L Alonso-Llanes

G Sánchez-Colina

A J Batista-Leyva

C Clément

E Altshuler

R Toussaint
email: renaud.toussaint@unistra.fr

Sink versus tilt penetration into shaken dry granular matter: The role of the foundation

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

I. INTRODUCTION

The Kocalei earthquake occurring on August 17, 1999 affected in various ways many constructions in the city of Adapazari, Turkey. Following observers, some buildings sank vertically into the soil, others tilted, and some even suffered lateral translation over the ground [START_REF] Bray | Damage patterns and foundation performance in adapazari[END_REF][START_REF] Sancio | Correlation between ground failure and soil conditions in adapazari, turkey[END_REF][START_REF] Sancio | Performance of buildings over liquefiable ground in adapazari, turkey[END_REF]. This case illustrates well the diversity of damage that earthquake fluidization of soils may cause to man-made structures [START_REF] Ambraseys | Engineering seismology: Part i[END_REF].

Liquefaction in the ground may be triggered dynamically by waves emitted during earthquakes, generating cyclic shear stresses that lead to the gradual build-up of pore water pressure [START_REF] Ben-Zeev | The combined effect of buoyancy and excess pore pressure in facilitating soil liquefaction[END_REF][START_REF] Ben-Zeev | Compaction front and pore fluid pressurization in horizontally shaken drained granular layers[END_REF]. The shaking produced by seismic events is a trigger for extensive liquefaction, as was observed recently in Belgium [START_REF] Vanneste | Late quaternary earthquake-related softsediment deformation along the belgian portion of the feldbiss fault, lower rhine graben system[END_REF]. Ground fluidization [START_REF] Council | Liquefaction of Soils During Earthquakes[END_REF][START_REF] Berrill | Energy dissipation and seismic liquefaction of sands: revised model[END_REF] has been investigated in different kinds of media like sand [START_REF] Berrill | Energy dissipation and seismic liquefaction of sands: revised model[END_REF], dry granular soils [START_REF] Clement | Fluidization of a bidimensional powder[END_REF] and sediments [START_REF] Wang | Earthquakes and Water[END_REF]. Of immediate interest for engineering and for the geosciences is to understand how man-made structures such as buildings, and massive rocks laying on granular soils respond to fluidization associated to seismic waves.

Granular matter itself displays a variety of puzzling phenomena [START_REF] Jaeger | Granular solids, liquids, and gases[END_REF][START_REF] Shinbrot | Reverse buoyancy in shaken granular beds[END_REF][START_REF] Altshuler | Sandpile formation by revolving rivers[END_REF][START_REF] Aranson | Patterns and collective behavior in granular media: Theoretical concepts[END_REF][START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Altshuler | Revolving rivers in sandpiles: From continuous to intermittent flows[END_REF][START_REF] Boudet | The granular jump[END_REF][START_REF] Vinningland | Granular rayleigh-taylor instability: Experiments and simulations[END_REF][START_REF] Johnsen | Pattern formation during air injection into granular materials confined in a circular hele-shaw cell[END_REF][21][START_REF] Niebling | Sedimentation instabilities: impact of the fluid compressibility and viscosity[END_REF][START_REF] Niebling | Mixing of a granular layer falling through a fluid[END_REF][START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF][START_REF] Turkaya | Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium[END_REF], but during the last decade or so, our understanding of the dynamics of objects penetrating into granular media has advanced quickly [START_REF] Uehara | Low-speed impact craters in loose granular media[END_REF][START_REF] Boudet | Dynamics of impact cratering in shallow sand layers[END_REF][START_REF] Katsuragi | Unified force law for granular impact cratering[END_REF][START_REF] Goldman | Scaling and dynamics of sphere and disk impact into granular media[END_REF][START_REF] Pacheco-Vázquez | Infinite penetration of a projectile into a granular medium[END_REF][START_REF] Kondic | Microstructure evolution during impact on granular matter[END_REF][START_REF] Torres | Impact dynamics in "hard" and "soft" granular matter[END_REF][START_REF] Clark | Granular impact model as an energy-depth relation[END_REF][START_REF] Ruiz-Suárez | Penetration of projectiles into granular targets[END_REF][START_REF] Brzinski | Depth-dependent resistance of granular media to vertical penetration[END_REF][START_REF] Altshuler | Settling into dry granular media in different gravities[END_REF][START_REF] Joubaud | Forces encountered by a sphere during impact into sand[END_REF][START_REF] Harich | Intruder mobility in a vibrated granular packing[END_REF][START_REF] Kolb | Reorganization of a dense granular assembly: The unjamming response function[END_REF][START_REF] Li | A terradynamics of legged locomotion on granular media[END_REF][START_REF] Aguilar | Robophysical study of jumping dynamics on granular media[END_REF][START_REF] Kang | Archimedes' law explains penetration of solids into granular media[END_REF][START_REF] Díaz-Melián | Rolling away from the wall into granular matter[END_REF]. While laterally shaken granular beds have received a certain degree of attention [START_REF] Medved | Modes of response in horizontally vibrated granular matter[END_REF][START_REF] Kruelle | Physics of granular matter: pattern formation and applications[END_REF], the performance of objects initially laying on the surface of a granular bed submitted to lateral shaking has been rarely studied [START_REF] Tennakoon | Onset of flow in a horizontally vibrated granular bed: Convection by horizontal shearing[END_REF][START_REF] Liu | Seismic response of shallow foundation on liquefiable sand[END_REF][START_REF] Sánchez-Colina | Note:"lock-in accelerometry" to follow sink dynamics in shaken granular matter[END_REF][START_REF] Clément | Shake and sink: liquefaction without pressurization[END_REF][START_REF] Clément | Sinking during earthquakes: Critical acceleration criteria control drained soil liquefaction[END_REF].

In this paper we perform systematic experiments associated to the latter scenario, which may help understanding the performance of human constructions and rocks laying on granular beds during earthquakes. In particular, using a cylinder as a simplified model, we study its settling dynamics on a granular bed submitted to lateral vibrations. Somewhat unexpectedly, we have found that, within a large range of lateral shaking powers, cylinders with flat bottoms sink vertically, while those with a "foundation"consisting in a shallow ring attached to their bottom, tilt besides sinking. The latter scenario seems to dominate independently from the nature of the foundation when strong enough lateral vibrations are applied. Quasi-2D simulations mimicking the experiments were also performed. The settling dynamics of the simulated intruders, with or without foundation, reproduces the corresponding experimental results. Our simulations also reveal how the difference in forcechain distributions between flat and non-flat bottom cylinders produces different torques justifying the two types of penetration. In addition, we present a simple phenomenological model that reproduces well the sinking dynamics and helps understanding how the tilting process influences the sinking one.

II. EXPERIMENTAL

The penetration experiments were performed on a granular bed contained in a test cell of approximately 25 × 25 × 25 cm 3 filled with Ugelstad spheres of non expanded polystyrene with a bulk density 1050 kg/m 3 , and diameter 140 µm (monodisperse within a 1 percent), type Dynoseeds ©, produced by Microbeads, Norway. The box was horizontally shaken at different amplitudes of motion (A), and a frequency (f ) of 5.0 Hz (a value commonly found in seismic waves), using a TIRA ©TV51120-M shaker, see figure 1. By controlling the voltage of the shaker input signal we varied the amplitude of the oscillations up to a maximum value corresponding to a peak ground acceleration of A(2πf ) 2 ≈ 12.2 m/s 2 [51]. This acceleration range covers most potentially damaging earthquakes, from weak to strong [START_REF]Shakemap scientific background[END_REF], though there has been reports of larger peak ground accelerations [START_REF] Goto | Extreme accelerations during earthquakes caused by elastic flapping effect[END_REF]. The time the shaker needs to reach the steady state depends on the dimensionless acceleration, being longer for larger accelerations. The time intervals can range from one (0.2 s) to three periods.

Two types of intruders were used in the experiments: (a) a hollow 3D printed cylinder of 44 mm diameter, 44 mm height h c (external dimensions), and 5 mm thick walls, and (b) the same cylinder with a ring of 5 mm height and 3 mm thickness glued to its bottom (illustrated in the upper right corner of Fig. 1). Intruders (a) and (b) will be called "No-ring"and "Ring", respectively, from now on. Their masses were adjusted with ballast in such a way that their densities matched the average effective density of the granular medium, which was measured as 430 kg/m 3 . As far as the ballast used has a density near the effective density of the granular material, it was almost evenly distributed inside the cylinder. Note that, using a flat bottom cylinder and a ring-like bottom cylinder, we are modifying the "foundation"of our intruder.

A digital camera Hero 2 made by GoPro © was fixed to the electromagnetic shaker, in such a way that it could take a video of the sinking process from an oscillating reference frame locked to the test cell, as proposed in [START_REF] Sánchez-Colina | Note:"lock-in accelerometry" to follow sink dynamics in shaken granular matter[END_REF]. This method allowed a much better quality of the cylinder's images, and made easier their digital processing. Videos were taken at a maximum rate of 120 frames per second, with a resolution of 1920 × 1080 pixels.

The images were processed as follows. We first converted the videos to image sequences binarized through an appropriate threshold. Using the tool regionprops from MatlabR2014a ©, we identified and assigned coordinates to several bright marks we had glued to certain points of the cylindrical intruder. The coordinates of the marks were used to calculate the position of the intruder's geometrical center and inclination relative to the vertical in each picture. In some experiments where the sinking was particularly big, it was difficult to obtain the tilt angle, since part of the marks sank below the level of the sand surface, and they were impossible to follow. In such cases the upper border of the cylinder was identified using the Matlab's tools find and bwtraceboundary, and then fitted to a polynomial using the function polyfit. The fit was used to find the inclination. In the case of experiments ending in a very inclined position, the reference to calculate the inclination was the cylinder's corner above the sand surface, that was identified as the intersection of the two polynomial fits of the upper and one lateral borders of the cylinder.

As the cylinder oscillates due to the vibration of the box, it is difficult to determine the final position, particularly when there is a big tilting. Then, in order to determine the sinking depth and tilting, we observe in the videos the onset of a cyclic movement of a reference point in the cylinder. Then, the final position could be measured in the frames filmed after the shaker was stopped.

The experimental protocol can be described as follows: (I) preparing the granular medium by stirring it evenly with a long rod, (II) gently depositing the cylinder in the upright position on the free surface of the granular bed, (III) turning ON the camera, (IV) switching ON the shaker after setting the desired frequency and amplitude (V) turning OFF the shaker and the camera after the penetration process had visibly finished.

In Fig. 2 we define the main parameters describing the sinking process of a No-ring cylinder (upper row), and the tilting and sinking of a Ring cylinder (bottom row), during the experimental lapse, defined as ∆t. As the figure indicates, in the following we will call h the penetration of the geometrical center at a time t and ∆h the final penetration at time t = ∆t. Note that both magnitudes are defined as the vertical displacement of the geometrical center of the cylinder (without taking the ring into account). In the same way we will call θ the inclination of the intruder at time t and ∆θ the final inclination at t = ∆t.

We also explored the phenomenology through numerical simulations. They were based on a discrete element method code (DEM) for the computation of granular systems [START_REF] Johnsen | Pattern formation during air injection into granular materials confined in a circular hele-shaw cell[END_REF][START_REF] Niebling | Sedimentation instabilities: impact of the fluid compressibility and viscosity[END_REF][START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF][START_REF] Clément | Shake and sink: liquefaction without pressurization[END_REF][START_REF] Clément | Sinking during earthquakes: Critical acceleration criteria control drained soil liquefaction[END_REF][START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Parez | Unsteady granular flows down an inclined plane[END_REF]. We modeled a quasi-2D granular medium, made of finite-sized hard spheres with radii between 1.0 and 1.5 mm, to avoid the crystallization effect. The medium contains 4000 particles and is prepared by placing the latter randomly in a space 30 cm wide and 25 cm high and then allowing them to settle under the action of gravity g = 9.81 m/s 2 .

Once the medium reaches equilibrium, it occupies a virtual space 30 cm wide and about 8 cm high, laterally delimited by flat walls that define the Hele-Shaw cell. The simulation box used was created narrow in order to have a single plane of particles in the direction perpendicular to the images shown in Figures 3 and4. The components of the velocities and forces along this direction are set to zero at each time step. To mimic the experimental conditions, we simulate particles of density 1050 kg/m 3 .

The two intruders are made of cohesive particles. One is a square of 40 mm side, made of N = 1681 particles with diameter 1 mm, placed in a quasi-2D square arrangement, which simulates the No-ring intruder of the experiments. The second one is also a square of 40 mm side to which two small feet are attached. Each foot has a size of 4.5 × 2.7 mm 2 , so the simulated Ring intruder contains a total of N = 1705 particles. The density of the spheres ρ p which form the intruders is adjusted so that the bulk density of the rigid body matches the effective density of the granular medium. The latter is calculated once the medium has settled down and is stable, and was always found to be around ρ m = 566 kg/m 3 . Then, the density of the particles forming the intruder is obtained as ρ p = ρ m V i /N V p , where V i is the volume of the intruder and V p the volume of a sphere.

Once our granular medium is created, we place the intruder 1 mm above the medium, with its bottom parallel to the horizontal direction. We release it, under the action of the force of gravity, and wait until the whole system becomes motionless (i.e. its total kinetic energy reaches a value under 10 -7 J). Then, we apply horizontal oscillations of different amplitudes and a frequency of 5 Hz (∆t ≈ 8s) to the walls of the cell and compute the time evolution of the position and tilting angle of the intruder. The amplitudes were chosen in such a way that the dimensionless acceleration Γ = A(2πf ) 2 /g had the values of 0.16, 0.25, 0.5, 1.0, 1.25 and 1.5.

The contact between spheres was modeled as a linear spring-dashpot [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Thompson | Lammps -a flexible simulation tool for particlebased materials modeling at the atomic, meso, and continuum scales[END_REF],

F ij = (k n δn ij - m ef f γ n v n ) -(k t ∆S t + m ef f γ t v t )
where k and γ are the elastic and viscoelastic damping constants, δn ij is the overlap distance along the line connecting the centers of the two spheres, and v their relative velocity. ∆S t is the tangential displacement vector between two spheres, which is truncated to satisfy a frictional yield criterion, and m ef f = m i m j /(m i + m j ) is the effective mass of two spheres of mass m i and m j . We considered normal (n) and tangential (t) forces components between the particles and, in order to model hard spheres that interact on contact (i.e. spheres whose deformation during collisions is less than a small fraction of their radii), we used the following parameter values:

k n = 1.2 × 10 7
N/m, k t = 2/7k n , γ n = 12 s -1 and γ t = 0.1γ n . The interaction force between the walls and the particles touching them is the same as the corresponding for two particles but considering the wall of infinite radius and mass (flat wall). The microscopic friction coefficient between spheres, and between spheres and boundaries, was taken as µ = 0.3. The time step dt was chosen to guarantee that there are at least 50 steps during one characteristic time of a collision dt = t c /50, where

t c = π/ (k n /m ef f ) -γ 2 n .
Figures 3 and4 show the initial and final positions of both types of intruders in two typical runs. Fig. 3 indicates that the No-ring cylinders are slightly inclined, while in Fig. 4 the large inclination of the Ring one becomes obvious.

We also performed an additional set of simulations aimed at elucidating the influence of the intruder's aspect ratio in the sink -tilt behaviour. Two intruders, one No-ring with aspect ratio 1.125, and a Ring with aspect ratio 1 (see Fig. 5), were submitted to the same range of dimensionless accelerations. Note that the new intruders (in dark gray in Fig. 5)

have the same dimensions of the former ones (in light gray) but the geometry of the bottom is interchanged.

The new intruders are also rigid bodies made of 1 mm diameter cohesive particles and their bulk densities also correspond to that of the granular medium. In light gray the intruders previously described and in dark gray the new ones.

III. RESULTS AND DISCUSSION

A. Sink vs. tilt penetration in experiments An important characteristic of the sinking process in this type of cylinder is that the intruder penetrates the granular medium with almost no tilting, and a final inclination smaller than 5°.

Figure 6(b) is similar to the previous one, but measurements were performed with Ring cylinders. The general features of both graphics are similar, but there is a difference, that will be better observed in the following figures: the dimensionless acceleration at which the cylinder sinks completely in the medium is bigger for the Ring cylinders than for the No-ring ones. The long-time creep process is not completely shown. The tilting angle of No-ring cylinders is not displayed, due to the fact that it oscillates around angles not larger than 5°relative to the vertical direction.

shaking, while cylinders with rings tend to tilt. Figure 7 quantifies the differences between the initial and final stages of the process, for almost all the range of accelerations our experimental setup was able to reach. granular bed. Vertical penetrations started to increase significantly above this acceleration, reaching a plateau around Γ ≈ 0.7. At the plateau, the cylinder has sunk completely, but stays "floating" into the fluidized granular medium, as expected for an object isodense relative to it, so there is no further sinking.

In Fig. 7(b) the sinking process of the Ring cylinders is summarized. Though the low acceleration part is similar to Fig. 7(a), now the plateau is not observed for the range of accelerations recorded. Notice that from the depth reached at Γ ≈ 1.2, approximately the height of the cylinder h c , it would not sink any further, and that from this value of Γ onwards a plateau would appear. The sinking process can be understood taking the experimental results in Ref. [START_REF] Tennakoon | Onset of flow in a horizontally vibrated granular bed: Convection by horizontal shearing[END_REF] into account. When the system is submitted to lateral shaking, a fluidized layer appears in the upper part of the granular cell. This layer reaches a depth h f that depends on the dimensionless acceleration Γ. Below this layer exists a "solid" layer. For accelerations in the range spanned in our experiments, h f varies almost linearly with Γ (see Fig. 3(a) in Ref. [START_REF] Tennakoon | Onset of flow in a horizontally vibrated granular bed: Convection by horizontal shearing[END_REF]), so we can write

h f (Γ) = α(Γ -Γ * ); Γ > Γ * (1)
where Γ * is the onset of fluidization and α is the slope of the linear dependence. If Γ ≤ Γ * the depth of the fluidized layer is zero.

Then, at low values of Γ the granular medium is not fluidized, and the cylinder almost does not sink (merely 5 mm at Γ = 0.27; see Fig. 7(a)). For accelerations above the fluidization threshold, the cylinder sinks until it gets in contact with the solid layer. The larger is the acceleration, the deeper is that layer, so the bigger is ∆h. But as soon as the solid layer appears at a depth larger than the cylinder's height, it does not sink further: instead, it "floats" due to isodensity with the sand, so a plateau is reached.

According to reference [START_REF] Toussaint | Soil liquefaction: how the granular medium evolves, macroscale and microscale study[END_REF], Γ * can be taken as proportional to the friction coefficient µ between the cylinder and the granular medium. In these experiments we can approximate µ ≈ 0.3, which is the value we use in the simulations. The authors also conclude that the final depth of intrusion depends on isostasy, and on the severity of shaking. It can be entirely determined by isostasy, when the shaking completely unjam the medium and suppresses the average friction around the intruder [START_REF] Clément | Sinking during earthquakes: Critical acceleration criteria control drained soil liquefaction[END_REF].

To better understand the differences in the dynamics of both types of intruders, we performed numerical simulations and their results are described below.

B. Sink vs tilt in quasi-2D numerical simulations Regarding the vertical sinking in Fig. 8(a), we do not observe major changes between Ring and No-ring intruders; both types of intruders sink less than in the experiments. This may be related with the lower dimensionality of the simulations relative to the real experiment.

Quasi-2D granular media allow less choices of readjustment than in 3D: they are easily jammed, which makes it more difficult for an object to sink. Moreover, the size ratio of the intruder over the particles is 8 times smaller in the simulations than in the experiments (experiments: 44 mm/0.140 mm ≈ 300; simulations: 40 mm/1 mm= 40), which means that if one particle is stuck under the intruder during the simulations, it will slow down the intruder significantly more than if the particle were 8 times smaller. values of Γ. The conclusions obtained from Fig. 8 are valid for all the dimensionless accelerations used in the simulations: there are no significant differences in the final sinking depth between both types of intruders.

In Fig. 10, on the contrary, the difference in tilting angles between the two types of intruders is clearly seen. For all the values of Γ the simulated Ring intruder tilts more than the No-ring one. For Γ=1.50, the No-ring tilts up to an angle that is closer to the Ring's one, corresponding to preliminary observations found in experiments with frequencies above 5 Hz for dimensionless accelerations Γ >1.25.

According to our simulations, the difference in tilting between intruders lies in that one type of intruder, the No-ring one, is somehow more capable of rectifying its rotation during sinking, while the other, the Ring intruder, is not. This rectification can be understood as the process of returning to, or recovering, the initial rotation angle once one oscillation of the cell has concluded and, as can be seen in the temporal evolution of θ (Fig. 11(a)), the difference in the rotation angles between the No-ring and Ring intruders is produced by a non-rectifying cumulative process taken by the latter.

To understand why the tilting dynamics is affected by the presence of the legs, which makes the No-ring intruder able to further rectify its rotation -at least for the values of gamma between 0.25 and 1.25-, we calculate from the simulations the torque about the center of mass and the angular velocity. Fig. 11 shows the time evolution of θ, the angular velocity, the torque about the center of mass and the oscillation of the cell for the two types of intruders in a simulation with Γ = 1.0. Focusing on the θ curve, the difference in tilting can be noticed during the first second of the simulation (as in Fig. 8 (b) for the averaged values) as well as the aforementioned Ring intruder non-rectification process. Unexpectedly, the values of the torques shown in Fig. 11 (c) are very similar for the two intruders, contradicting the intuitive idea that the sole presence of the legs would produce higher torques about the Ring intruder's center of mass making it rotate more. However, the slight differences may affect, to a greater or lesser extent, the rotational movement of the latter with respect to the No-ring intruder (see Fig.

(b)).

In general, during the first half of a cell oscillation, both intruders rotate in the same direction -counterclockwise due to the fact that the granular bed moves from left to right-.

But, during the second half, some forces appear in the bottom and/or the inside of the Ring intruder's right leg that do not allow it to rotate in the same way as the No-ring does.

In some cases, these forces completely prevent it from rotating clockwise, as in both first oscillations shown in Fig. 11. This process is prone to occur in each of the oscillations during the simulation and its repetition causes the differences in inclination observed for the two intruders after 8 simulated seconds (see Fig. 10). to the depth h, as expressed in Eq. (A.13). Though this model does not include the degree of freedom associated with tilting, the consideration of its influence in the lineal and surface dimensions of the intruder is enough to explain why a tilted intruder reaches a final depth smaller than that reached by a non tilted intruder for a given Γ. Now, we examine the influence of the intruder's aspect ratio on the penetration dynamics.

The results of the simulations performed with the intruders of different aspect ratios are summarized in Fig. increasing aspect ratio. Figures 14(c) and (d) compare intruders with the same foundation and different aspect ratios. Both figures lead to the same conclusion: the larger the aspect ratio, the larger the tilt angle. Summarizing the results of Fig. 14, the intruder with higher aspect ratio and ring placed on the bottom has the largest tilt angle for all dimensionless accelerations, and the intruder with flat bottom and lower aspect ratio has the smallest tilt angle. Interestingly, the other two intruders show approximately equal final tilt angle values for equal values of Γ.

We have studied the behavior of cylindrical objects as they sink into a dry granular bed fluidized by horizontal oscillations, as a model system to understand the effects of earthquake-related fluidization of soils on human constructions and other objects like rocks.

We have found that, within a relatively large range of lateral shaking amplitudes at a frequency of 5 Hz, cylinders with flat bottoms sink vertically, while those with a "foundation"consisting in a shallow Ring attached to their bottom, tilt laterally besides their vertical sinking.

We have been able to mimic the above described behaviors by quasi-2D numerical simulations. With their help we found that these differences are not necessarily due to the sole existence of the legs that generate bigger torques about the center of mass. Instead, they can be associated with the jamming of the particles in the region between the legs, which can increase the friction at the bottom of the Ring intruder generating force chains capable of preventing the total recovery of its initial angle of rotation. Numerical experiments also helped to clarify the influence of the intruder aspect ratio on the tilt dynamics: of two intruders with the same foundation, the one with higher aspect ratio will have a larger tilt angle. So, the aspect ratio and the foundation type cooperate to establish the penetration dynamics of the intruder.

We have also reproduced the vertical sink dynamics of cylinders with a flat base using a Newtonian equation of motion for an object penetrating a fluidized layer of granular matter, where the granular effective density increases with depth, eventually reaching a solid phase.

The same model allows to understand the sinking even in the present of tilting (Appendix 1).

Finally, it is worth noting that preliminary experimental data and quasi-2D numerical simulations suggest that, when strong enough lateral shaking is applied, the tilting scenario tends to dominate regardless the nature of the intruder's foundation.

where h, as previously, is the depth reached by the cylinder below the surface of the granular medium. In Eq. (A.2) we have made explicit that the density of the material varies with depth. Let us assume that it varies as a power law between zero and the density of the solid layer, ρ sl , that is reached at a depth h f :

ρ(z ) = ρ sl z h f p (A.3)
where p ∈ [0, 1]. The selection of the value of p is discussed below (see also Appendix 2).

By combining (A.3) and (A.2) and integrating, we find the hydrostatic buoyancy force acting on the cylinder with a length h under the (average) level of the granular bed, as:

(-P )ndS = - ρ sl Sg (p + 1)h p f h p+1 ĥ (A.4)
where S is the characteristic area of the intruder cross section, and ĥ is a unit vector pointing downwards. It is easy to see that the buoyancy force depends on the volume submerged into the granular medium.

Neglecting the inertial forces, which according to our simulations is typically two orders of magnitude smaller than the contact forces, the shear stress component goes as

σ s • ndS = -Dγv ĥ (A.5)
where γ has the dimensions of a viscosity, D is the characteristic size of the cross section of the intruder and v is its sinking speed [START_REF] De Bruyn | Penetration of spheres into loose granular media[END_REF][START_REF] Hou | Projectile impact and penetration in loose granular bed[END_REF]. By substituting Eq. (A.4) and Eq. (A.5) into Eq. (A.1), and only recovering the modular values, we get:

m d 2 h dt 2 + Dγ dh dt + ρ cs Sg (p + 1)h p f h p+1 = mg (A.6)
Before solving Eq. (A.6) we will assume that the sink velocity is constant, which follows quite well the behavior during the fast sink regime, as seen in Fig. 15 (i.e., we neglect the inertial term). So,

dh dt + ρ sl Sg Dγ(p + 1)h p f h p+1 = mg Dγ (A.7)
which can be written as if p = 0, and

h(t) = b a tanh( √ abt) (A.10) if p = 1.
It is easy to see that both expressions correspond to an exponential growth that saturates. Then, we solve this equation numerically. This procedure is repeated for values of p between 0 and 1, with a step of 0.1.

The inset in Fig. 15 shows some of the numerical solutions for the values of p in the legend. The main conclusion is that the density profile has small influence on the first (and most important) part of the sinking process. Of course, the final depth is influenced by the value of p, but due to experimental uncertainties, it is almost impossible to choose any particular value.

Let us now study the influence of the values of p in the quality of the fit of the solution of Eq. (A.6) to the experimental data. For doing this we notice that the values of a and b in Eq. (A.9) can be easily obtained from the experiments. Considering Eq. (A.8) in the first moments of motion, as h is small, h (t) b, so b can be evaluated as the initial slope. As at large times h(t) ∼ h eq (h eq = ∆h if h(0) = 0) then a = b/h p+1 eq . Then solving Eq. (A.6) for a given value of b, p, a(p) and naming the result h mod , the best value of p arises from the minimization:

p opt = arg min N i=1 [h mod (t i , p) -h exp (t i )] 2 (A.11)
where h exp (t) are the experimental values of h.

The result for Γ ≤ 1.0 is indifferent to p: the fit is equally good no matter which is the value of p ∈ [0, 1]. For Γ = 1.24 there are differences in the quality of the fits for various values of p, but Eq. (A.11) gives a minimum for p = 0, so, we will use p = 0 in the following (Appendix 2 supports the selection of p from the simulations). Then, Eq. (A.6) becomes:

m d 2 h dt 2 + Dγ dh dt + ρ sl Sg h f h = mg (A.12)
that can be taken as the simplest equation of motion describing the vertical sink dynamics of our cylinders. It is worth noticing that Eq. (A.12) reproduces quite closely the results reported in Fig. 15, and can be used to qualitatively describe the vertical sinking of Ring-cylinders while tilting, as we will see below. It is possible to demonstrate that Eq. (A.12), developed for a granular bed fluidized by shaking, is closely related with that proposed in [START_REF] Pacheco-Vázquez | Infinite penetration of a projectile into a granular medium[END_REF] to describe the penetration of an intruder into ultra-light granular material that eventually behaves like a fluid medium even in the absence of shaking.

In order to understand the influence of the tilting dynamics in the sinking process, it is useful to note that, when applying Eq. (A.12) to a tilted cylinder, the values of both D and S change. The reason is that when we calculate the surface integrals, the result will be proportional to the cylinder's immersed volume. As the cylinder tilts, both the immersed surface and linear dimensions increase more than in the case of sinking without tilting, so the drag force is bigger in the former case. Considering, for instance, the situation represented in the lower row of Fig. 2, when the cylinder sinks a distance ∆h, the surface and linear dimensions increase as the inverse of cos θ (of course, other intruder geometries may follow different laws).

To test it, let us assume a simplified model: the increase factor of S and D is proportional to the characteristic size of the cross section of the cylinder projected on the horizontal plane, i.e., it is proportional to the inverse of cos θ. Then, instead of D and S, we will solve Eq.

(A.12) using D/ cos θ(t) and S/ cos θ(t), where θ(t) is a function that grows from zero to the maximum angle θ max reached by the cylinder, mimicking Fig. 6(c), i. e. with the same time constant.

The consequences can be seen in Fig. 16. While in the beginning the sinking process in all situations occurs with the same dynamics, as the cylinder approaches the final angle, the behavior changes, being the final depth larger for the situations corresponding to low tilting.

The upper curve, calculated for θ = 0 coincides with the upper curve in the inset of Fig. 15 (calculated for p = 0). Subsequent curves are calculated for values of θ max varying in steps of π/15, the lowermost curve corresponds to θ max = π/3. As the inclination of the cylinder increases, both the buoyancy and the viscous drag do. The effect of these factors on the sinking process of Ring cylinders was already noted in Fig. 6(b): an immediate consequence is the decrease of the sinking depth (for a given Γ), compared with that of the No-ring ones, which can be easily observed in the experiments. From the inset it is possible to deduce that, for the larger angles, a small decrease in the depth is observed.

In spite of the simplifications assumed, it is worth noting that one of the basic differences Finally, there is another element that was neither considered by us: as the container shakes horizontally, it produces a horizontal drag that changes periodically its direction.

According to [START_REF] Zhang | Numerical investiga[END_REF], it creates an additional lift force, and also a dependence of the drag force with depth, which, of course, must influence the detailed penetration dynamics of the Ring cylinders. The results of Li et al [START_REF] Li | A terradynamics of legged locomotion on granular media[END_REF] also support these ideas.

Density profile of the granular medium

At each time step the simulation box is subdivided into a fixed number of rectangles for which the density is calculated as ρ = m/V , where m is the sum of the masses of all the particles of the granular medium within each rectangle and V is the volume of the rectangle calculated as its area multiplied by the average diameter of the particles of the medium inside it.

The density profile for each dimensionless acceleration Γ is obtained by calculating the average of the density profile at each time step. The latter is obtained by averaging the profiles of a set of rectangle columns near the intruder. These columns are chosen avoiding the presence of the intruder in them in order not to affect the density profile since only the particles that compose the granular medium are taken into account in its calculation. Γ. As can be seen in the region of the graph enclosed by dashed lines, the density of the medium saturates rapidly with increasing depth for all dimensionless accelerations. This fact supports from the simulations the use of p = 0 in expression A.3 of Appendix 1.

FIG. 1 .

 1 FIG. 1. Experimental setup. At the upper right, we have illustrated the intruder consisting in a cylinder with ring.

FIG. 2 .

 2 FIG.2. Sketch of sinking and tilting processes. The top row illustrates the sinking process of a No-ring cylinder in three moments during the experimental interval from t = 0 to a final time t = ∆t. The bottom row shows the same temporal sequence for a Ring cylinder, which tilts in addition to sink.

FIG. 3 .

 3 FIG. 3. (color online) Snapshots of the initial (a) and final (b) positions of a No-ring intruder in a typical quasi-2D simulation using a shaking frequency of 5 Hz.

  The new Ring intruder is composed by N = 1498 particles forming a 40 mm wide and 35.5 mm high rectangle to which two 4.5 × 2.7 mm 2 feet were added. Note that the size of these feet and those of the Ring intruder are the same. The new No-ring is a rectangle of 40 mm wide and 44.5 mm high composed of N = 1845 particles.

8 4FIG. 5 .

 85 FIG.5. Comparison of the geometry of the intruders used to test the influence of the aspect ratio.

Figure 6 (

 6 Figure 6(a) shows the time variation of the sinking depth for selected values of the dimensionless acceleration Γ = A(2πf ) 2 /g (where g = 9.81 m/s 2 is the gravitational acceleration and A(2πf ) 2 is the horizontal peak acceleration of the sand box) for No-ring cylinders. It is easy to see that the penetration of the No-ring cylinders follows a common pattern for all the accelerations. A first process of fast sinking is followed by a slow creep. Only the penetration depth increases with Γ. In this figure we do not show the total creep process, due to its long duration. As the height of the cylinder is h c = 44 mm, it is possible to check from Fig. 6(a) that, for a dimensionless accelerations of 1.24, the cylinder sinks completely.

Figure 6 (FIG. 6 .

 66 Figure 6(c) presents the time evolution of the tilting angle for a Ring cylinder, a process that occurs simultaneously with the sinking. The sinking and tilting dynamics of Ring cylinders is more irregular than that of the No-ring ones. This is illustrated in Fig.6(b) and (c), even after being submitted to a moving average process, to get a smoother graph.No-ring cylinders tend to sink vertically as the granular soil is fluidized by horizontal

Figure 7 (FIG. 7 .

 77 Figure 7(a) shows sink data for No-ring cylinders. As can be seen, for dimensionless accelerations up to Γ = 0.27, there was no significant penetration of the intruder into the

Figure 7 (

 7 Figure 7(c) shows the tilt data for Ring cylinders. No significant tilting is observed for Γ smaller than approximately 0.25. With the increase of the dimensionless accelerations, the cylinder significantly tilts, increasing abruptly the tilting angle with Γ, until it slows down at Γ ≈ 0.7. We do not show the tilting angle of No-ring cylinders, because it is always smaller than 5°, with a random distribution of values around the vertical direction.Figures7(b) and (c) are closely related, because they are two descriptions of the same process: the motion of Ring cylinders into the granular medium, that includes both sinking and tilting. The fact that at the accelerations shown in this figure the plateau in the sinking depth is not reached while for the tilting angle at higher values of Γ the inclination almost saturates, could be explained by the increase of the friction of the intruder with the granular medium when the tilting angle increases. At Γ ≈ 0.7 the intruder has reached a large inclination, but is not completely submerged into the medium. An increase in the acceleration does not increase significantly the angle, because the resulting torque has diminished due to the influence of both sinking and tilting, but the increase in fluidization helps further sinking, until most of the cylinder is submerged.

Figure 8

 8 Figure 8 shows the time dependence of the penetration depth (a) and tilting angle (b) for both types of intruders at the dimensionless acceleration Γ = 1.0. In both figures the thick curves represent the average value of six repetitions varying the initial conditions and the surrounding zone represents ±1 standard deviation.

Figure 8 (

 8 Figure 8(b) indicates that the presence of a foundation at the bottom of the intruder causes a large tilting. Indeed, for the shaking with no ring, the intruder tilting angle is around 10°, but during the shaking with ring, the intruder tilts fast, reaching an angle around 50°. This resembles what happens in the experiments (see Fig. 4): the intruder almost ends up lying on one of its sides. Of course, the tilting is also limited by the diminished dimensionality in the quasi 2D simulations.

Figure 9 FIG. 8 .FIG. 9 .FIG. 10 .

 98910 Figure9compares the penetration depth reached for both types of cylinders at different

FIG. 11 .

 11 FIG. 11. (color online) Simulations. Temporal evolution of the calculated magnitudes for the intruders during the first second of a simulation using Γ = 1.0. From top to bottom: the evolution of the rotation θ (a), the angular velocity Ω (b), the torque about the center of mass τ (c) and finally, the position of the cell (d) showing the oscillation described by it. The vertical line indicates the instant t = 1.84 T. τ max is maximum torque on a non-tilting but horizontally accelerated No-ring intruder.

FIG. 12

 12 FIG. 12. (color online) Simulations. Pressure field in the granular medium represented as a jet color map for No-ring (a) and Ring (c) intruders. The color of the intruders represents their angular velocity using a seismic color map. In (b) and (d) the outer segments of the intruders are represented as a whole with a color corresponding to the resulting torque about the center of mass that is generated on all grains belonging to them. Notice that for the Ring intruder each leg is divided into 3 segments. The grains in these segments are magnified for better viewing. The resulting torque about the center of mass is represented as the interior color of the intruder. For both angular velocity and torque, blue color represents clockwise.

Fig. 12 illustrates

 12 Fig.12illustrates in more detail what is described above. In (a) and (c) it shows the pressure field in the granular medium at t = 1.84 T (time indicated in Fig.11by the vertical line), where the force chains are represented using a sequential color map. In them, the color of the intruders represents their angular velocity, which in turn is displayed as a diverging color map where blue indicates clockwise rotations. Figures12 (b) and (d)show the contribution to the torque about the center of mass of each of the intruder segments as a result of the forces acting on them. In this case, the color of the intruders is associated with the resulting torque about the center of mass. Note that all the grains of each outer edge in

FIG. 13 .

 13 FIG. 13. (color online) Simulations. Grain tracking. Initial (left) and final (right) positions of the No-ring (top) and Ring (bottom) intruders as well as the grains of the granular medium next to them. The grains right below the intruders are represented in blue while the yellow and cyan colors represent the grains initially at the left and right side of the intruders, respectively. Note that for the Ring intruder, the blue grains initially located between the legs remain in the same region until the end of the simulation

FIG. 14 .

 14 FIG. 14. Simulations. Final tilt angles at different Γ for two intruders with (a) aspect ratio 1 and different foundations; (b) aspect ratio 1.1125 and different foundations; (c) a flat bottom and different aspect ratios and (d) a ring at the bottom and different aspect ratios.

FIG. 15 .

 15 FIG. 15. (color online) Time dependence of sinking depth for the No-ring cylinder from experiment, compared with that determined from Eqs. (A.9, A.10). The inset shows the solutions of Eq. (A.6) for different values of p (see text).

Figure 15 shows

 15 Figure 15 shows the experimental results (continuous line) obtained for a dimensionless acceleration Γ 0.3. It is possible to see in more detail the initial fast sinking process, followed by the slow creep. Fig. 15 also shows the fitting of equations (A.9, A.10) to experimental data. Both solutions reproduces well the main features of the sinking process.It is almost impossible to determine experimentally the exact density profile. But we do not need to know it in order to validate our model, if we use the following rationale.

FIG. 16 .

 16 FIG. 16. (color online) Time dependence of sinking depth as calculated solving numerically Eq. (A.12) considering the variation of S and D provoked by tilting (see text). Upper curve is for θ max = 0 while the lower one is for θ max = π/3. Between them, θ max varies in steps of π/15. The inset shows the last three seconds.

Figure 17 shows 25 FIG. 17 .

 172517 Figure17shows the density profiles as a function of depth for the different values of

ACKNOWLEDGMENTS

L. A. thanks the Photovoltaic Research Lab of the University of Havana for allowing the use of its facilities. E. A. drew inspiration from the late M. Álvarez-Ponte. We acknowledge support from Project 29942WL (Fonds de Solidarité Prioritaire France-Cuba), from the EU 262644. We thank Mustapha Meghraoui, Einat Aharonov, Knut Jørgen Måløy, Eirik G.

Flekkøy, for fruitful discussions. This research was made in the frame of the University of Havana's institutional project "Medios granulares: creando herramientas para prevenir catástrofes".

Appendix 1. Sink dynamics: a phenomenological Newtonian model

The model to be formulated should account for two related processes, the sinking in the vertical direction and the oscillations perpendicular to it. But as was shown above, the Noring intruders have only small oscillations that end fast, being the overall sinking process almost vertical. Regarding the Ring intruders, though they strongly oscillate, the tilting process ends first, so we will consider only the equation controlling the vertical sinking, figuring out how the tilting angle affects the sinking dynamics.

In order to formulate a model to describe analytically the sinking process, let us consider the forces acting on the cylinder. As soon as the shaking starts, if the dimensionless acceleration is above threshold, the upper part of the granular bed is fluidized, and the intruder sinks.

Let us assume that the cylinder just sinks vertically, and let us name the vertical downward axis as z. The force balance on the intruder can be written as

where P is the pressure, σ s the shear stress tensor, n is the vector normal to the intruder's surface, and the integrals run over the boundary of the intruder that is inside the granular material. Assuming a hydrostatic pressure profile, we can write:


Pierre-Jean Meyer
email: pierre-jean.meyer@univ-eiffel.fr

Reachability analysis of neural networks using mixed monotonicity

Keywords:

This paper presents a new reachability analysis approach to compute interval over-approximations of the output set of feedforward neural networks with input uncertainty. We adapt to neural networks an existing mixedmonotonicity method for the reachability analysis of dynamical systems and apply it to each partial network within the main network. This ensures that the intersection of the obtained results is the tightest interval over-approximation of the output of each layer that can be obtained using mixed-monotonicity on any partial network decomposition. Unlike other tools in the literature focusing on small classes of piecewise-affine or monotone activation functions, the main strength of our approach is its generality: it can handle neural networks with any Lipschitz-continuous activation function. In addition, the simplicity of our framework allows users to very easily add unimplemented activation functions, by simply providing the function, its derivative and the global argmin and argmax of the derivative. Our algorithm is compared to five other interval-based tools (Interval Bound Propagation, RELUVAL, NEURIFY, VERINET, CROWN) on both existing benchmarks and two sets of small and large randomly generated networks for four activation functions (ReLU, TanH, ELU, SiLU).

I. INTRODUCTION

A RTIFICIAL intelligence and particularly neural networks are quickly spreading in various fields including safetycritical applications such as autonomous driving [START_REF] Xiang | Verification for machine learning, autonomy, and neural networks survey[END_REF]. With it comes a growing need for replacing verification methods based on statistical testing [START_REF] Kim | A programmatic and semantic approach to explaining and debugging neural network based object detectors[END_REF] by more formal methods able to provide safety guarantees of the satisfaction of desired properties on the output of the neural network [START_REF] Liu | Algorithms for verifying deep neural networks[END_REF]. Such methods are also useful for the verification of a closed-loop system with a neural network controller [START_REF] Akintunde | Reachability analysis for neural agent-environment systems[END_REF], [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF].

When focusing on isolated neural networks, the rapidly growing field of formal verification has been categorized into three main types of objectives [START_REF] Liu | Algorithms for verifying deep neural networks[END_REF]: counter-example result to find an input-output pair that violates the desired property; adversarial result to determine the maximum allowed disturbance that can be applied to a nominal input while preserving the properties of the nominal output; reachability result to evaluate (exactly or through over-approximations) the set of all output values reachable by the network when given a bounded input set. While the first two types of results often rely on solving optimization problems [START_REF] Liu | Algorithms for verifying deep neural networks[END_REF], [START_REF] Salman | A convex relaxation barrier to tight robustness verification of neural networks[END_REF], reachability results (which are considered in this paper) naturally lean toward developing new reachability analysis approaches to be applied layer by layer in the network. The current reachability methods in the literature consider various set representations to overapproximate the network's output set with static intervals [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF], symbolic interval equations and linear relaxations of activation functions (RELUVAL [START_REF] Wang | Formal security analysis of neural networks using symbolic intervals[END_REF], NEURIFY [START_REF]Efficient formal safety analysis of neural networks[END_REF], VERINET [START_REF] Henriksen | Efficient neural network verification via adaptive refinement and adversarial search[END_REF], CROWN [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF]) or polyhedra and zonotopes in the tool libraries NNV [START_REF] Tran | NNV: The neural network verification tool for deep neural networks and learning-enabled cyber-physical systems[END_REF] and ERAN [START_REF]ERAN: ETH robustness analyzer for neural networks[END_REF].

The main weakness of existing neural network verifiers is that most of them are limited to a very small class of activation functions. The large majority of verifiers only handles piecewise-affine activation to focus on the most popular ReLU function [START_REF] Wang | Formal security analysis of neural networks using symbolic intervals[END_REF], [START_REF]Efficient formal safety analysis of neural networks[END_REF], [START_REF] Katz | The marabou framework for verification and analysis of deep neural networks[END_REF], [START_REF] Botoeva | Efficient verification of relu-based neural networks via dependency analysis[END_REF]. A handful of tools cover other popular activations, such as S-shaped functions sigmoid and hyperbolic tangent [START_REF] Henriksen | Efficient neural network verification via adaptive refinement and adversarial search[END_REF], [START_REF] Tran | NNV: The neural network verification tool for deep neural networks and learning-enabled cyber-physical systems[END_REF], [START_REF]ERAN: ETH robustness analyzer for neural networks[END_REF]. Several other tools claiming to handle general activation functions are either implicitly restricted to monotone functions in their theory or implementation [START_REF] Dvijotham | A dual approach to scalable verification of deep networks[END_REF], [START_REF] Raghunathan | Certified defenses against adversarial examples[END_REF], or their claimed generality rather refers to a mere compatibility with general activation functions but they require the user to provide their own implementation of how to handle these new functions [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF].

Although ReLU or sigmoid networks often perform well, the focus on these activation functions in the literature is also strongly tied to the fact that most current optimization and verification tools cannot handle more general activations. Proposing a new approach to deal with these general activation functions is thus an important step towards opening the field of formal verification of neural networks to new non-monotone activations that have already been shown to outperform classical ReLU or sigmoid networks (see e.g. Swish/SiLU [START_REF] Ramachandran | Searching for activation functions[END_REF], Mish [START_REF] Misra | Mish: A self regularized non-monotonic activation function[END_REF], PFLU and FPFLU [START_REF] Zhu | PFLU and FPFLU: Two novel non-monotonic activation functions in convolutional neural networks[END_REF]).

Contributions: This paper presents a novel method for the reachability analysis of feedforward neural networks with general activation functions. The proposed approach adapts the existing mixed-monotonicity reachability method for dynamical systems [START_REF] Meyer | Interval Reachability Analysis: Bounding Trajectories of Uncertain Systems with Boxes for Control and Verification[END_REF] to be applicable to neural networks so that we can obtain an interval over-approximation of the set of all the network's outputs that can be reached from a given input set. Since our reachability method is also applicable to any partial network within the main neural network and always returns a sound over-approximation of the last layer's output set, intersecting the interval over-approximations from several partial networks ending at the same layer can only tighten the approximation while preserving the soundness of the results. To take full advantage of this, we propose an algorithm that applies our new mixed-monotonicity reachability method to all L(L + 1)/2 partial networks contained within the considered L-layer network. This ensures that the resulting interval over-approximation is the tightest output bound of the network obtainable using mixed-monotonicity reachability analysis on any partial network decomposition, at the cost of a computational complexity of O(L 2 ).

Mixed-monotonicity reachability analysis is applicable to any system whose Jacobian matrix is bounded on any bounded input set [START_REF] Meyer | Interval Reachability Analysis: Bounding Trajectories of Uncertain Systems with Boxes for Control and Verification[END_REF]. In the case of neural networks combining linear transformations and nonlinear activation functions, the above requirement implies that our approach is applicable to any network with Lipschitz-continuous activation function. Since extracting the bounds of the Jacobian matrix of a system (or of the derivative of an activation function in our case) is not always straightforward, for the sake of self-containment of the paper we provide a method to automatically obtain these bounds for a (still very general) sub-class of activation functions whose derivative can be defined as a 3-piece piecewise-monotone function. Apart from the binary step and Gaussian function, all activation functions we could find in the literature (including all non-monotone functions reviewed in [START_REF] Zhu | PFLU and FPFLU: Two novel non-monotonic activation functions in convolutional neural networks[END_REF]) belong to this class.

Although most of the tools cited above provide a two-part verification framework for neural networks (one reachability part to compute output bounds of the network, and one part to iteratively refine the network's domain until a definitive answer can be given to the verification problem), this paper focuses on providing a novel method for the first reachability step only. The proposed approach can thus be seen either as a preliminary step towards the development of a larger verification framework for neural networks when combined with an iterative splitting of the input domain, or as a standalone tool for the reachability analysis of neural networks as is used for example in the analysis of closed-loop systems with a sample-and-hold neural network controller [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF].

In summary, the main contributions of this paper are:

• a novel approach to soundly bound the output of neural networks using mixed-monotonicity reachability; • a method that is compatible with any Lipschitzcontinuous activation function; • a framework that allows to easily add new activation functions (the user only needs to provide the activation function, its derivative, and the global arg min and arg max of the derivative).

The paper is organized as follows. Section II defines the considered problem and provides useful preliminaries for the main algorithm, such as the definition of mixedmonotonicity reachability for a neural network, and how to automatically obtain local bounds on the activation function derivative. Section III presents the main algorithm applying mixed-monotonicity reachability to all partial networks within the given network. Finally, Section IV compares our novel reachability approach to the bounding methods of 5 other interval-based optimization-free tools [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF], [START_REF] Wang | Formal security analysis of neural networks using symbolic intervals[END_REF], [START_REF]Efficient formal safety analysis of neural networks[END_REF], [START_REF] Henriksen | Efficient neural network verification via adaptive refinement and adversarial search[END_REF], [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF] on both existing benchmarks [START_REF] Bak | The second international verification of neural networks competition[END_REF] and randomly generated networks, to highlight the complementarity with these tools and the cases when our approach outperforms them.

II. PRELIMINARIES Given x, x ∈ R n with x ≤ x, the interval [x, x] ⊆ R n is the set {x ∈ R n | ∀i ∈ {1, . . . , n}, x i ≤ x i ≤ x i }.

A. Problem definition

Consider the L-layer feedforward neural network

x (i) = Φ(W (i) x (i-1) + b (i) ), ∀i ∈ {1, . . . , L} (1)
where x (i) ∈ R ni is the output vector of layer i, x (0) and x (L) being the input and output of the neural network, respectively. We assume that the network is pre-trained and all weight matrices W (i) ∈ R ni×ni-1 and bias vectors b (i) ∈ R ni are known. The function Φ is defined as the componentwise application of a scalar and Lipschitz-continuous activation function. For simplicity of presentation, the activation function Φ is assumed to be identical for all layers, but note that the proposed approach in this paper is compatible with having different activation functions for each layer (or even for each node of the same layer). In the literature, the activation function of the last layer is often omitted since a monotone activation (the most commonly used) would only change the output values but not their relative comparison (which is considered for classification problems). Since our approach introduces the ability to handle non-monotone activation, our network definition (1) needs to allow for the activation of all layers. Some verification problems on neural networks aim to measure the robustness of the network with respect to input variations. Since the output set of the network cannot always be computed exactly, we rely on over-approximating it by a simpler set representation, such as an interval. We can then solve the verification problem using this interval: if the desired property is satisfied on the over-approximation, then it is also satisfied on the real output set of the network. In this paper, we focus on the problem of computing an interval overapproximation of the output set of the network when its input is taken in a known bounded set, as formalized below.

Problem 1: Given the L-layer neural network (1) and the interval input set [x 0 , x 0 ] ⊆ R n0 , find an interval [x L , x L ] ⊆ R n L over-approximating the output set of (1):

{x (L) | x (0) ∈ [x 0 , x 0 ]} ⊆ [x L , x L ].
Naturally, our secondary objective is to ensure that the computed interval over-approximation is as tight as possible in order to minimize the number of false negative results in the subsequent verification process.

B. Mixed-monotonicity reachability

In this paper, we solve Problem 1 by iteratively computing over-approximation intervals of the output of each layer, until the last layer of the network is reached. These overapproximations are obtained by adapting to neural networks existing methods for reachability analysis of dynamical systems. More specifically, we rely on the mixed-monotonicity approach in [START_REF] Meyer | Interval Reachability Analysis: Bounding Trajectories of Uncertain Systems with Boxes for Control and Verification[END_REF] to over-approximate the reachable set of any discrete-time system x + = f (x) with Lipschitz-continuous vector field f : R n → R n . Since the neural network ( 1) is instead defined as a function y = f (x) with input x = x (0) and output y = x (L) of different dimensions, we present below the generalization of the mixed-monotonicity method from [START_REF] Meyer | Interval Reachability Analysis: Bounding Trajectories of Uncertain Systems with Boxes for Control and Verification[END_REF] to any Lipschitz-continuous function.

Consider the function y = f (x) with input x ∈ [x, x] ⊆ R nx , output y ∈ R ny and Lipschitz-continuous function f : R nx → R ny . The Lipschitz-continuity assumption ensures that the derivative f (also called Jacobian matrix in the paper) is bounded.

Proposition 2: Given an interval [J, J] ⊆ R ny×nx bounding the derivative f (x) for all x ∈ [x, x], denote the center of the interval as J * . For each output dimension i ∈ {1, . . . , n y }, define input vectors ξ i , ξ i ∈ R nx and row vector α i ∈ R 1×nx such that for all j ∈ {1, . . . , n x },

(ξ i j , ξ i j , α i j ) = (x j , x j , min(0, J ij )) if J * ij ≥ 0, (x j , x j , max(0, J ij )) if J * ij ≤ 0.
Then for all x ∈ [x, x] and i ∈ {1, . . . , n y }, we have:

f i (x) ∈ f i (ξ i ) -α i (ξ i -ξ i ), f i (ξ i ) + α i (ξ i -ξ i ) .
Intuitively, the output bounds are obtained by computing for each output dimension the images for two diagonally opposite vertices of the input interval, then expanding these bounds with an error term when the bounds on the derivative f spans both negative and positive values. Proposition 2 can thus provide an interval over-approximation of the output set of any function as long as bounds on the derivative f are known. Obtaining such bounds for a neural network is made possible by computing local bounds on the derivative of its activation functions, as detailed in Section II-C.

C. Local bounds of activation functions

Proposition 2 can be applied to any neural network whose activation functions are Lipschitz continuous since such functions have a bounded derivative. However, to avoid requiring the users to manually provide these bounds for each different activation function they want to use, we instead provide a framework to automatically and easily compute local bounds for a very general class of functions describing most popular activation functions and their derivatives.

Let R ∞ = R∪{-∞, +∞} and consider the scalar function ϕ : R ∞ → R ∞ , where ϕ(x) ∈ {-∞, +∞} only when x ∈ {-∞, +∞}. In this paper, we focus on 3-piece piecewisemonotone functions for which there exist z, z ∈ R ∞ such that ϕ is:

• non-increasing on (-∞, z] until reaching its global minimum min x∈R∞ ϕ(x) = ϕ(z); • non-decreasing on [z, z] until reaching its global maximum max x∈R∞ ϕ(x) = ϕ(z); • and non-increasing on [z, +∞). When z = -∞ (resp. z = +∞), the first (resp. last) monotone segment is squeezed into a singleton at infinity and can thus be ignored.

Although this description may appear restrictive, we should note that the vast majority of activation functions as well as their derivatives belong to this class of functions. In particular, this is the case for (but not restricted to) popular piecewise-affine activation functions (identity, ReLU, leaky ReLU), monotone functions (hyperbolic tangent, arctangent, sigmoid, SoftPlus, ELU) as well as many non-monotone activation functions (SiLU, GELU, HardSwish, Mish, REU, PFLU, FPFLU, which are all reviewed or introduced in [START_REF] Zhu | PFLU and FPFLU: Two novel non-monotonic activation functions in convolutional neural networks[END_REF]).

Example 1: Two examples of such functions are provided in Figure 1. The first one (in black) is the Sigmoid Linear Unit (SiLU, also called Swish [START_REF] Ramachandran | Searching for activation functions[END_REF]) activation function Φ(x) = x/(1+e -x ), whose global arg min is z = -1.2785. Its global arg max is z = +∞, which implies that the third monotone component is pushed toward +∞ since it is not needed for this function.

The second example (in dashed red) is the derivative of the SiLU activation function Φ (x) = (1 + e -x + xe -x )/(1 + e -x ) 2 , with global arg min and arg max defined as z = -2.3994 and z = 2.3994, respectively.

The only exceptions of activation functions whose derivative do not belong to this class that the author could find in the literature are: the binary step, which is discontinuous so its derivative is undefined in 0; and the Gaussian activation function Φ(x) = e -x 2 which belongs to this class, but not its derivative (although -Φ does belong to this class, so a very similar approach can still be applied).

Proposition 3: Given a scalar function ϕ as defined above and a bounded input domain [x, x] ∈ R, the local bounds of ϕ on [x, x] are given by:

min x∈[x,x] ϕ(x) = ϕ(z) if z ∈ [x, x], min(ϕ(x), ϕ(x)) otherwise, max x∈[x,x] ϕ(x) = ϕ(z) if z ∈ [x, x], max(ϕ(x), ϕ(x)) otherwise.
In this paper, Proposition 3 is used to obtain local bounds of the activation function derivatives in order to compute bounds on the network's Jacobian matrix for Proposition 2. However, Proposition 3 can also be useful on activation functions themselves to compare (in Section IV) our method with the naive interval propagation through the layers.

III. REACHABILITY ALGORITHM

Solving Problem 1 by applying the mixed-monotonicity approach from Proposition 2 on the neural network (1) can be done in many different ways. One possibility is to iteratively compute bounds on the Jacobian matrix of the network and then apply Proposition 2 once for the whole network. However, this may result in wide bounds for the Jacobian of the whole network, and thus in a fairly conservative over-approximation of the network output due to the term α i in Proposition 2. The dual approach is to iteratively apply Proposition 2 to each layer of the network, which would result in tighter Jacobian bounds (since one layer's Jacobian requires less interval matrix products than the Jacobian of the whole network), and thus tighter over-approximation of the layer's output. However, the loss of the dependency with respect to the network's input may induce another source of accumulated conservativeness if we have many layers. Any intermediate approach can also be considered, where we split the network into several consecutive partial networks and apply Proposition 2 iteratively to each of them.

Although all the above approaches result in overapproximations of the network output, we cannot determine in advance which method would yield the tightest bounds, since this is highly dependent on the considered network and input interval. We thus devised an algorithm that encapsulates all possible choices by running Proposition 2 on all partial networks within (1), before taking the intersection of the obtained bounds to tighten the bounds on each layer's output. The main steps of this approach are presented in Algorithm 1 and summarized below.

Given the L-layer network (1) with activation function Φ and input interval [x 0 , x 0 ] as in Problem 1, our goal is to apply Proposition 2 to each partial network of (1), denoted as NN(k, l), containing only layers k to l (with k ≤ l) and with input x (k-1) and output x (l) . We initialize the Jacobian bounds of each partial network to the identity matrix. Then, we iteratively explore the network (going forward), where for each layer l we first use interval arithmetics [START_REF] Jaulin | Applied interval analysis: with examples in parameter and state estimation, robust control and robotics[END_REF] to compute the pre-activation bounds based on the knowledge of the output bounds of the previous layer, and then apply Proposition 3 to obtain local bounds on the activation function derivative (ϕ = Φ ). Next, for each partial network NN(k, l) (with k ≤ l) ending at the current layer l, we compute its Jacobian bounds based on the Jacobian of layer l ([Φ , Φ ] * W (l) ) and the Jacobian of the partial network NN(k, l-1). We then apply Proposition 2 to the partial network NN(k, l) with the Jacobian bounds we just computed and input bounds [x (k-1) , x (k-1) ]. Finally, once we computed the l over-approximations of x (l) corresponding to each partial network ending at layer l, we take the intersection of all of them to obtain the final bounds for x (l) .

Theorem 4: The interval [x (L) , x (L) ] returned by Algorithm 1 is a solution to Problem 1. This is easily proved by the fact that we compute several interval over-approximations of the output of each layer using Proposition 2, which ensures that their intersection is still an over-approximation of the layer's output. Then, using these over-approximations as the input bounds of the next layers Input: L-layer network (1), input interval [x 0 , x 0 ], activation function Φ (with Φ defined as in Section II-C) ∀k, l ∈ {1, . . . , L}, J(k, l) ← I, J(k, l) ← I for l ∈ {1, . . . , L} do

[Φ , Φ ] ← Prop3(Φ , W (l) * [x (l-1) , x (l-1) ] + b (l) ) for k ∈ {1, . . . , l} do [J(k, l), J(k, l)] ← [Φ , Φ ] * W (l) * [J(k, l -1), J(k, l -1)] [x(k, l), x(k, l)] ← Prop2(NN(k, l), [x (k-1) , x (k-1) ], [J(k, l), J(k, l)]) [x (l) , x (l ) ] ← [x(1, l), x(1, l)] ∩ • • • ∩ [x(l, l), x(l, l)]
Output: Over-approximation [x (L) , x (L) ] of the network output Algorithm 1: Mixed-monotonicity reachability analysis of a feedforward neural network.

guarantees the soundness of the approach.

Note also that in Algorithm 1, Proposition 2 is applied to every partial network that exists within the main network [START_REF] Xiang | Verification for machine learning, autonomy, and neural networks survey[END_REF]. Although this implies a computational complexity of O(L 2 ), it guarantees that the resulting interval [x (L) , x (L) ] from Algorithm 1 is the least conservative solution to Problem 1 that could be obtained from applying the mixed-monotonicity reachability analysis of Proposition 2 to any decomposition of (1) into consecutive partial networks.

IV. NUMERICAL COMPARISONS

To evaluate the performances of our approach, we run Algorithm 1 on both established benchmarks and randomly generated neural networks, and compare the results with stateof-the-art tools for reachability analysis and verification of neural networks. 1 Since our method is an optimization-free reachability analysis approach using interval bounds, we focus these numerical comparisons to the most relevant toolboxes of this category, and we leave to future work the additional comparisons with less related verification tools relying on optimization-based approaches [START_REF] Liu | Algorithms for verifying deep neural networks[END_REF] or reachability analysis with more complex set representations (polyhedra, zonotopes) as in the ERAN [START_REF]ERAN: ETH robustness analyzer for neural networks[END_REF] and NNV libraries [START_REF] Tran | NNV: The neural network verification tool for deep neural networks and learning-enabled cyber-physical systems[END_REF].

We provide comparisons to the following 5 tools and methods for reachability analysis of neural networks. We first use as a baseline the naive interval bound propagation (IBP), as presented e.g. in [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF], since this is the simplest approach using interval arithmetic at each layer, but it is also very conservative due to losing dependency to the network's inputs at each step. Three of the other tools (RELUVAL [START_REF] Wang | Formal security analysis of neural networks using symbolic intervals[END_REF], NEURIFY [START_REF]Efficient formal safety analysis of neural networks[END_REF], VERINET [START_REF] Henriksen | Efficient neural network verification via adaptive refinement and adversarial search[END_REF]) rely on methods propagating symbolic intervals (i.e. bounding functions linear in x (0) ), sometimes combined with linear relaxations of the nonlinear activation functions.

The fifth tool (CROWN [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF]) instead relies on a backward propagation of these linear relaxations.

Since these tools are written in various programming languages (Matlab, C, Python), for the convenience of comparison we have re-implemented the bounding method of each tool in Matlab. The implemented approaches are those described in each of the original papers cited above, and may thus slightly differ from the latest updates of the corresponding public toolboxes. Note also that the comparisons in this section only focus on reachability analysis as formulated in Problem 1 (over-approximating the network's reachable set for a given input uncertainty), and we thus disregard the iterative refinement algorithms that some of the above tools use to verify or falsify properties on the network's behavior.

VNN benchmarks: We first compare our mixedmonotonicity approach with the above 5 reachability tools on the MNISTFC benchmark as presented in the 2021 VNN competition [START_REF] Bak | The second international verification of neural networks competition[END_REF]. This benchmark is composed of 3 feedfoward ReLU neural networks trained on the MNIST dataset for the recognition of handwritten digits. All 3 networks have 784 inputs (28 × 28 pixels of a gray-scale picture), an output layer with 10 nodes, and an additional 2, 4 or 6 hidden layers of 256 nodes each. We run all 6 reachability methods on these 3 networks for 250 different input pictures, each with an = 0.03 uncertainty around their nominal values. In terms of tightness of the output bounds, our approach is 10 to 1000 tighter than IBP and has a comparable order of magnitude with the other four tools. In the top half of Table I, the first three columns summarize the percentage of the 250 evaluated input pictures for which our mixed-monotonicity approach results in tighter (or at least as tight) output bounds than the other five methods. We can see in particular that our method outperforms all others in at least 70% cases on the benchmark with two hidden layers. In terms of computation times (bottom of Table I), our approach is significantly slower than IBP and VERINET and a bit slower than CROWN, but it is also up to 3 times faster than RELUVAL and NEURIFY.

We ran the same comparison on the only non-ReLU feedforward benchmark of the 2021 VNN competition [START_REF] Bak | The second international verification of neural networks competition[END_REF]: the ERAN-sigmoid benchmark is also trained on the MNIST dataset and has the same structure as the above MNISTFC networks, but with 6 hidden layers of 200 nodes each, sigmoid activation functions and input uncertainty = 0.012. On the same 250 input pictures, we observe that: IBP is always more conservative than our approach (as expected); RELUVAL and NEURIFY cannot be run on non-ReLU networks; VERINET and CROWN both fail for every 250 inputs because their results are so conservative on this benchmark that Matlab evaluates e -x = ∞ when the pre-activation value x goes too far in the negative, and the (theoretically horizontal) tangent of the sigmoid e -x /(1 + e -x ) 2 cannot be computed. VERINET and CROWN's failure due to excessive conservativeness on this sigmoid benchmark is also observed when using an input uncertainty = 1.2 * 10 -5 which is 1000 smaller than the one suggested in the VNN benchmark [START_REF] Bak | The second international verification of neural networks competition[END_REF].

Random networks: To get more representative comparisons on a wider variety of network dimensions (since the above Method mnistfc-2 mnistfc-4 mnistfc-6 ERAN-sig IBP [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF] 100% 100% 100% 100% RELUVAL [START_REF] Wang | Formal security analysis of neural networks using symbolic intervals[END_REF] 80% 74.8% 65.6% -NEURIFY [START_REF]Efficient formal safety analysis of neural networks[END_REF] 72.4% 72% 58.8% -VERINET [START_REF] Henriksen | Efficient neural network verification via adaptive refinement and adversarial search[END_REF] 70.4% 54% 15.6% 100% CROWN [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF] 70.4% 54% 15.6% 100% IBP [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF] 0.012 0.019 0.025 0.016 RELUVAL [START_REF] Wang | Formal security analysis of neural networks using symbolic intervals[END_REF] 42 62 85 -NEURIFY [START_REF]Efficient formal safety analysis of neural networks[END_REF] 43 62 84 -VERINET [START_REF] Henriksen | Efficient neural network verification via adaptive refinement and adversarial search[END_REF] 0.012 0.019 0.026 0.017 CROWN [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF] 0 benchmarks are only 4 specific networks), and more importantly to highlight the main strength of our approach handling general activation functions (since we could not find wellestablished benchmarks using any activation other than ReLU and Sigmoid), we run these comparisons over two large sets of randomly generated neural networks.

We first consider a set of 10000 small feedforward networks as defined in (1) and whose parameters are randomly chosen as follows: depth L ∈ {1, . . . , 5}; input and output dimensions n 0 , n L ∈ {1, . . . , 10}; width of hidden layers (each layer may have a different width) n 1 , . . . , n L-1 ∈ {1, . . . , 30}. The second set contains 1000 deeper and wider networks with: depth L ∈ {5, . . . , 10}; input dimension n 0 ∈ {500, . . . , 1000}; output dimension n L ∈ {10, . . . , 50}; width of hidden layers n 1 , . . . , n L-1 ∈ {100, . . . , 200}. The input bounds as in Problem 1 are defined as the hypercube [x (0) , x (0) ] = {x * } + [-0.1, 0.1] n0 around a randomly chosen center input x * ∈ [-1, 1] n0 . The comparison between our mixed-monotonicity method in Algorithm 1 and the five approaches from the literature cited above (when applicable) is done for each of these 4 activation functions:

• Rectified Linear Unit (ReLU) is the piecewise-affine function Φ(x) = max(0, x), handled by all 6 tools; • Hyperbolic tangent (TanH) is the S-shaped function Φ(x) = (e x -e -x )/(e x + e -x ), handled by all tools apart from RELUVAL and NEURIFY;

• Exponential Linear Unit (ELU) is the monotone function Φ(x) = e x -1 if x ≤ 0 and Φ(x) = x if x ≥ 0,
natively handled only by IBP and our method. ELU is compatible with the framework of VERINET and CROWN (although not implemented) at the condition of providing a method to compute linear equations bounding Φ on any bounded input set. We added this to our own Matlab implementation of these tools; • Sigmoid Linear Unit (SiLU) is the non-monotone function Φ(x) = x/(1 + e -x ), only natively handled by our method. Although the implementation of IBP in [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF] is restricted to monotone activation functions, our own Matlab implementation of IBP is extended to handle nonmonotone activation by using Proposition 3.

All these activation functions are Lipschitz continuous and their derivatives satisfy the desired shape described in Sec-tion II-C. Algorithm 1 can thus be applied to all of them. For the set of 10000 small networks, Table II gives the average computation times of each methods for each activation function. For these time averaging to be meaningful despite the varying sizes of the networks, we actually divide the computation times by the total number of neurons in the considered network, and only then we take the average. Table III gives the proportion of the 10000 small networks for which our mixed-monotonicity approach results in at least as tight output bounds as the other methods. Table IV gives the same result but only when our approach is strictly tighter. For the 1000 large networks, the same results are displayed in Tables V-VII. Note that the time tables are in microseconds for small networks (Table II), and milliseconds for large ones (Table V). In the next three paragraphs, we analyse these results and highlight the main take-away conclusions from these numerical comparisons in terms of generality of the methods, complexity, and tightness of the output bounds.

Among all 6 compared methods, our mixed-monotonicity approach is the most general one and the only one that could natively handle all considered activation functions. Indeed, RELUVAL and NEURIFY are limited to piecewise-affine functions (such as ReLU) and VERINET and CROWN cannot handle non-monotone functions (such as SiLU). In addition, although we added an implementation for IBP to handle SiLU, and for VERINET and CROWN to handle ELU, the original tools do not natively handle these activation functions. (The corresponding results are given in parentheses in Tables II-VII.) This last comment brings up another important aspect of the generality of our mixed-monotonicity approach: the ease of implementation to add new activation functions. Indeed, NEURIFY, VERINET and CROWN rely on linear relaxations of the nonlinear activation functions, which may require long and complex implementations to be provided by the user for each new function (e.g. finding the optimal relaxations of sigmoid-shaped functions takes several hundreds of lines of code in the implementation of VERINET). In contrast, for a user to add a new activation type to be used within our mixed-monotonicity approach, all they need to provide is the definition of the activation function and its derivative, and the global arg min and arg max of the derivative as defined in Section II-C. Everything else is automatically handled internally by Algorithm 1.

This generality and ease of use however comes at the cost of a higher computational complexity of O(L 2 ) since Algorithm 1 calls the mixed-monotonicity reachability method from Proportion 2 on all L(L + 1)/2 partial networks within the main L-layer network. On the set of smaller networks, this implies that our approach is the slowest, with a similar order of magnitude as CROWN (which also has a complexity of O(L 2 )), while the other four methods are 10 to 20 times faster (Table II). When the size of the network increases, IBP and VERINET are always the fastest methods and keep a consistent average computation time (per number of neurons in the network), while all other methods see their computation time per neuron increase, the worst being RELUVAL and NEURIFY which become slower than our mixed-monotonicity approach (Table V). Finally, in terms of tightness (measured as the 2-norm width of the output bounds), we first observe that, as expected, our approach is always at least as tight as the naive IBP method (Tables III andVI), and even strictly tighter in 70-80% cases in the set of small networks (Table IV). Apart from the very conservative IBP, it should be noted that none of the other 5 methods always outperforms the others. For example on the set of 10000 small ReLU networks, our approach is strictly tighter than the other four in 12-36% cases, equal in 31% cases, and strictly looser in the remaining 32-57% cases (Tables III andIV). The respective strengths of each method thus make them complementary. On average, our approach returns at least as tight output bounds as the state-of-the-art tools in less than half cases on smaller networks, but this percentage significantly increases on larger networks except for ReLU (Tables III andVI). When the size of the network increases, we also observe a larger proportion of cases where the compared methods return approximately equal output bounds, particularly on activation functions with approximate saturations, such as TanH (Tables VI andVII). However, one interesting case if with the ELU activation (for which the stateof-the-art tools where not designed and optimized), where our approach not only outperforms both VERINET and CROWN on almost all 1000 tested large networks, but it also obtains strictly tighter output bounds than VERINET in 79% cases, and CROWN is seen to fail to return a numerical result in all cases due to its excessive conservativeness and its formulation in [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF] being incompatible with horizontal linear relaxations (similarly to the ERAN-sigmoid benchmark).

V. CONCLUSIONS

This paper presents a new method for the sound reachability analysis of feedforward neural networks using mixedmonotonicity. The main strength of our approach is its very broad applicability to any network with Lipschitz-continuous activation functions, which is satisfied by the vast majority of activation functions currently in use in the literature. Another

Method

ReLU TanH ELU SiLU IBP [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF] 0.016 0.018 0.016 (0.018) RELUVAL [START_REF] Wang | Formal security analysis of neural networks using symbolic intervals[END_REF] 44 ---NEURIFY [START_REF]Efficient formal safety analysis of neural networks[END_REF] 44 ---VERINET [START_REF] Henriksen | Efficient neural network verification via adaptive refinement and adversarial search[END_REF] 0.034 0.05 (0.037) -CROWN [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF] 4 significant strength of our framework is the greater simplicity to implement new activation functions compared to existing tools relying on linear relaxations. We compared our approach with the naive Interval Bound Propagation and 4 state-ofthe-art tools using interval reachability analysis on neural network (RELUVAL, NEURIFY, VERINET, CROWN) on both benchmarks from the VNN competition [START_REF] Bak | The second international verification of neural networks competition[END_REF] and sets of randomly generated neural networks of a wide variety of depth and width. From these comparisons, we observed a good complementarity between our approach and the state-of-the-art tools (when applicable on the considered activation functions), where each outperforms the others on a significant percentage of the tested inputs and networks. Although this tool can already be used on its own for the reachability analysis of neural networks, our next objective is to address its main limitation of a O(L 2 ) complexity to make it more compatible with iterative refinement approaches similar

Method

ReLU TanH ELU SiLU IBP [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF] 100% 100% 100% (100%) RELUVAL [START_REF] Wang | Formal security analysis of neural networks using symbolic intervals[END_REF] 100% ---NEURIFY [START_REF]Efficient formal safety analysis of neural networks[END_REF] 2% ---VERINET [START_REF] Henriksen | Efficient neural network verification via adaptive refinement and adversarial search[END_REF] 0.1% 100% (98%) -CROWN [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF] 0.1% 100% (100%) -

Method

ReLU TanH ELU SiLU IBP [START_REF] Xiang | Reachable set estimation for neural network control systems: A simulation-guided approach[END_REF] 0% 0% 0% (0%) RELUVAL [START_REF] Wang | Formal security analysis of neural networks using symbolic intervals[END_REF] 0% ---NEURIFY [START_REF]Efficient formal safety analysis of neural networks[END_REF] 0% ---VERINET [START_REF] Henriksen | Efficient neural network verification via adaptive refinement and adversarial search[END_REF] 0% 0% (79%) -CROWN [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF] 0% 2% (100%) - to existing neural network verification algorithms.

Fig. 1 .

 1 Fig. 1. SiLU activation function (black line) and its derivative (red dashed).

  PROPORTION OF THE 250 MNIST INPUT PICTURES FOR WHICH OUR APPROACH GIVES OUTPUT BOUNDS TIGHTER THAN (OR EQUAL TO) OTHER METHODS. BOTTOM: AVERAGE COMPUTATION TIME IN s.

		.78	4.1	10	5.7
	Mixed-Mono	15	36	69	39
			TABLE I
	TOP:

  TABLE IV PROPORTION OF THE 10000 SMALL RANDOM NETWORKS FOR WHICH OUR MIXED-MONOTONICITY APPROACH RESULTS IN STRICTLY TIGHTER OUTPUT BOUNDS THAN OTHER METHODS.

	Method	ReLU TanH	ELU	SiLU
	IBP [5]	73%	71%	79%	(79%)
	RELUVAL [7]	36%	-	-	-
	NEURIFY [8]	16%	-	-	-
	VERINET [9]	12%	11%	(19%)	-
	CROWN [10]	12%	10%	(17%)	-

TABLE V AVERAGE

 V COMPUTATION TIME (PER NEURON IN THE NETWORK) IN ms FOR THE SET OF 1000 LARGE NETWORKS.

		.8	4.8	(4.8)	-
	Mixed-Monotonicity	33	29	34	33

TABLE VI PROPORTION

 VI OF THE 1000 LARGE RANDOM NETWORKS FOR WHICH OUR MIXED-MONOTONICITY APPROACH RESULTS IN TIGHTER (OR EQUAL) OUTPUT BOUNDS THAN OTHER METHODS.

TABLE VII PROPORTION

 VII OF THE 1000 LARGE RANDOM NETWORKS FOR WHICH OUR MIXED-MONOTONICITY APPROACH RESULTS IN STRICTLY TIGHTER OUTPUT BOUNDS THAN OTHER METHODS.

All computations are done on a laptop with 1.80GHz processor and 16GB of RAM running Matlab

2021b. The code used to generate all numerical comparisons described in this section is available at: https://gitlab.com/pj_meyer/MMRANN