Jérôme Sueur
email: sueur@mnhn.fr

Bernie Krause

Almo Farina

Climate change is breaking Earth's beat

Keywords: biodiversity, climate, sound, ecoacoustics

teaching and research institutions in France or abroad, or from public or private research centers.

Text

From pole to pole, numerous animal species transmit sounds for communication, foraging, and navigation, in particular. Sound is a major element of animal behavior that is used, for example, in territory defense, mate attraction, parent-offspring identification, orientation, prey localization, predator escape. Natural sounds may also come, in part, from plants, initiated both internally and as a consequence of external abiotic forces such as rain, wind, and moving water. Here, we emphasize that climate change is transforming two components of the global soundscapes (Figure 1): biophony which includes biotic sounds; and geophony that groups abiotic, but natural, sounds coming, for example, from running water, wind, and earth movements.

Weather properties play a role in the propagation of animal sound in the surrounding environment as sound speed depends on a range of parameters including ambient temperature, relative humidity, wind, and rain intensity. Variations in these parameters can drastically change the range of communication between individuals [START_REF] Larom | The influence of surface atmospheric conditions on the range and area reached by animal vocalizations[END_REF] unusual weather conditions might prevent sound propagation, for example, the temperature might be too high or the air too dry. In addition, there are abundant examples that demonstrate how temperature directly governs the acoustic behavior of terrestrial, aquatic, and marine animals through chemical and physiological regulation, in particular neuron activity and muscle contraction. It is well known that the ambient temperature controls the vocalization repetition rate, pitch, volume, and the auditory processes of ectothermic species, such as amphibians, fish, and arthropods [START_REF] Narins | Climate change and frog calls: long-term correlations along a tropical altitudinal gradient[END_REF]. The dependence of animal vocalization on climate is actually so strong that it is apparent in past and present popular culture: animals are thought to be able to predict weather as the thermometer crickets, the forecasting frogs, and the summer cicadas.

Individual sensitivity to meteorological conditions induces cascading effects at larger ecological scales. There is growing evidence that climate change is affecting the phenology of plant and animal populations [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF], a property that has a major role in both ecology and evolution. Acoustic phenology can be viewed as the seasonal acoustic activity of animals, a well-known and popular phenomenon for bird and whale songs for instance. Climate change can induce shifts in acoustic phenology either directly through a change in the phenology of signal-producing organisms or indirectly through a change in their resources (e.g., food plants). Similarly, the population spatial distribution can be modified by a change in climatic conditions with movements toward higher latitudes or altitudes [START_REF] Chen | Rapid range shifts of species associated with high levels of climate warming[END_REF], remodeling locally the population acoustics and soundscapes. Changes in population phenology and biogeography may then have relevant consequences in terms of nature acoustics.

The signal-producing capacity and distribution of a significant number of species may be then altered by climate change. Novel weather characteristics and dynamics might reduce, or increase, vocal emissions of soniferous species with consequences for other energy demanding activities. Some populations face the risk of extinction leading to silencing of the biophony. Alternatively, other populations may benefit from climate change, extend their distribution, and dominate the acoustic space. The loss and gain of soniferous species could thus strongly rearrange the natural soundscapes of large areas. This reorganization could lead to a loss of acoustic diversity among sites and across time and an acoustic homogenization with ubiquitous species dominating the soundscapes. This would be a reduction of beta diversity, a crucial component of biodiversity.

Acoustic changes might be even more pronounced for the species living in both cold and tropical areas that have a low thermal tolerance and a limited acclimation capacity [START_REF] Burivalova | The sound of a tropical forest[END_REF].

The acoustic niche hypothesis, derived from the ecological niche concept, stipulates that at a given location each species would occupy non-overlapping specific acoustic space to avoid interference. This partitioning observation is often noted in pristine environments but might not hold in perturbed environments, and may not prove to be viable with climate change due to the rearrangement of community composition, shifts in acoustic phenology, and to modifications of signal parameters as time pattern or frequency range. These changes would then lead to asynchrony and niche partition obfuscation.

Extreme climate events are already occurring and are predicted to increase in frequency and severity over time. These will not only damage the habitats of soniferous species [START_REF] Gordon | Habitat degradation negatively affects auditory settlement behavior of coral reef fishes[END_REF] but also deeply affect natural soundscapes, in particular, geophony. Cyclones, heat waves, droughts, floods, and wildfires generate soundscapes that will be more frequent than before.

For humans, changes in biophonic patterns can be a factor of anxiety and stress as natural soundscapes have been proven to be a source of tranquility and serenity, an acoustic refuge where the human mind can recover, and an important factor for green tourism and local economies [START_REF] Francis | Acoustic environments matter: Synergistic benefits to humans and ecological communities[END_REF].

Humans should be prepared for a change in acoustic environments, from remote wilderness to urban areas, and efforts in education should alert citizens to these future changes. The seasonal acoustic signatures will shift with noticeable changes in species density, and diversity. In particular, the populations of birds may show different song dynamics [START_REF] Oliver | Eavesdropping on the Arctic: automated bioacoustics reveal dynamics in songbird breeding phenology[END_REF], the populations of frogs may call earlier or less [START_REF] Gibbs | Climate warming and calling phenology of frogs near Ithaca, New York, 1900-1999[END_REF], and the populations of nocturnal insects may decrease inducing a reduction of stridulations [START_REF] Jeliazkov | Large-scale semi-automated acoustic monitoring allows to detect temporal decline of bush-crickets[END_REF]. If animal sounds are the material of several research programs in ecology and biological conservation, they have been also a source of inspiration for artists since the origin of music, in all styles. Sound travels across oceans and passes over country boundaries; as the acoustic medium of the natural world, it knows no boundaries. Such cultural importance and universality make natural soundscapes a fundamental means of increasing popular awareness about climate change. Anthropogenic noise and land use have been identified so far as major threats to animal acoustics [START_REF] Buxton | Noise pollution is pervasive in U.S. protected areas[END_REF]. Nonetheless, climate change also acts as a major risk for changing Earth's natural soundscapes through direct, indirect, and combined effects. If singing animals are proto-musicians interacting in a proto-orchestra, we can then predict that changing the thermal and moisture conditions of the concert hall will detune the ensemble and affect the performance capacities of the players. Furthermore, it will likely introduce a perceived dissonance in the music. The same is happening for Earth's orchestra: new atmospheric conditions are detuning natural sounds and only major mitigation actions will help preserve Earth's beat.

Box 1. The challenges of monitoring soundscapes over the long-term

For those researchers attuned to habitat-defining ecoacoustic indicators, there are many scattered reports of shifts in the density and diversity of organism voices due to climate change [START_REF] Sugai | Terrestrial passive acoustic monitoring: review and perspectives[END_REF].

However, these fluctuations are difficult to quantify, in part, because this subject is so new resulting in a paucity of quality data sets and analytic tools. But also, the older institutional sound libraries have primarily focused on individual species' acoustic samples rather than the collective expressions inherent in a more holistic biophonic approach; this last method is one that has evolved independently over the last half century. These biophonic collections must be mined for their baseline potential, especially where they have been recorded under strict protocols. In addition, ecoacousticians now have access to recording technologies that store layers of long-term metadata not previously possible with earlier systems. Because of their inherent value, we have reached a critical preservation turning point for such data acquisition and research.

  While solutions may exist to control anthropogenic sound by reducing the sound pressure level of sound sources or by protecting the receiver, there are unfortunately no straightforward solutions to reduce the effects of climate change on natural sound. The only upstream solutions are those advocated for climate mitigation. Nonetheless, current research in ecoacoustics can help in monitoring the changes affecting natural sounds. Noninvasive, long-term, and large-scale techniques are able to record and analyze natural soundscapes for ecological research, introduced in Box 1.

Figure 1 .

 1 Figure 1. Climate change and soundscape change. Changes in the soundscape of a natural environment over time. Top is a spectrogram revealing sound intensity in color according to time (year) and frequency. Bottom is an oscillogram showing sound amplitude according to time (year). This time-lapse is a series of 10-s samples extracted from recordings conducted from 2005 to 2019 with the exception of the year 2012. The recordings were obtained at the same location (Sugarloaf State Park, CA, USA, 38°26′23.00″N/122°29′53.44″W), at the same date (15 April), at the same time (06.35 h), with the same calibrated equipment. In red is overlaid the trend of the average temperature of April recorded at the same hour. This trend and 95% confidence intervals (broken lines) were obtained with a polynomial regression fitted over 30 years (1989-2019). Located about 50 m from a year-round stream, the site is an ecotone between riparian and oak savannah which did not suffer main land use changes. The lower half of the spectrogram shows the signature of the nearby stream with different flow regimes under different precipitation patterns. Weather-related changes reduced bird diversity from 2011. Recordings by Bernie Krause; temperature data from weatherspark.com.

Acknowledgements

We thank Sylvain Haupert for his advice, and the Parc Naturel Région du Haut Jura and the CNRS Guyane LEEISA for their support in conducting long-term acoustic research related to climate change. We are also grateful to two anonymous referees and to Kevin F.A. Darras for their valuable comments on the manuscript. This paper is dedicated to Stuart Gage (Emeritus, Michigan State University), a pioneer whose research, dedication, and insights were seminal to the emerging science of ecoacoustics.


Yannick Descantes
email: yannick.descantes@ifsttar.fr

Fabien Tricoire

Patrick Richard

Classical contact detection algorithms for 3D DEM simulations: drawbacks and solutions

The present paper sheds new light on ongoing drawbacks of three classical algorithms dedicated to initial contact detection between overlapping convex polyhedra, (i) Cundall's common plane, (ii) Nezami's fast common plane and (iii) Gilbert-Johnson-Keerthi's algorithm (GJK). Solutions to these drawbacks are suggested and implemented into revised versions of those three algorithms, which are further benchmarked for accuracy and speed using nine overlapping contact situations. The benchmarking results show that the revised version of GJK, called GJK -T D, and Nezami (revised) return values of the contact normal components and overlap depth which are identical to machine precision, whereas Cundall (revised) results differ beyond the ninth decimal place. Furthermore, for a given contact situation, GJK-TD returns those values within a few tens of microseconds on average, whereas Nezami (revised) and Cundall (revised) are respectively 6 and 65 times more computationally intense. It is believed that the robustness and efficiency of GJK-TD will boost its use into DEM simulations, all the more that this versatile algorithm may easily be customized to detect contact between convex polyhedra and spheroid particles.

Introduction

Discrete Element Methods (DEM) are widely used techniques to perform numerical simulations of granular assemblies aiming at investigating their flow behaviour [START_REF] Silbert | Granular flow down an inclined plane[END_REF][START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF] or packing properties [START_REF] Silbert | Geometry in Frictionless and Friction Sphere Packings[END_REF][START_REF] Farr | Close packing density of polydisperse hard spheres[END_REF] to mention just a few. Among existing literature, most studies simulate granular assemblies of spherical particles [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Makse | Packing of compressible granular materials[END_REF][START_REF] Torquato | Is Random Close Packing of Spheres Well Defined ?[END_REF][START_REF] Silbert | Granular flow down an inclined plane[END_REF][START_REF] Silbert | Geometry in Frictionless and Friction Sphere Packings[END_REF], whereas nonspherical particles have only recently began to be simulated using ellipsoids [START_REF] Donev | Underconstrained Jammed Packings of Nonspherical Hard Particles : Ellipses and Ellipsoids[END_REF], spherocylinders [START_REF] Williams | Random packings of spheres and spherocylinders simulated by mechanical contraction[END_REF], tetrahedra [START_REF] Zhao | Shape Effects on the Random-Packing Density of Tetrahedral Particles[END_REF], or other polyhedra [START_REF] Torquato | Dense Packings of Polyhedra: Platonic and Archimedean Solids[END_REF][START_REF] Azéma | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF][START_REF] Camenen | Effect of confinement on dense packings of rigid frictionless spheres and polyhedra[END_REF][START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF], despite their large presence in natural (e.g. ice debris, rocks) or manufactured state (e.g. sugar grains, crushed sand particles).

A well-known limitation of DEM applied to polyhedra is the lengthy calculation time needed to simulate large granular assemblies, say over a few tens of thousands particles, and still such assemblies only represent negligible quantities of granular materials (e.g. a 5 ml teaspoon already contains about 50000 powder sugar grains of size 0.5 mm). A significant share of this calculation time is dedicated to contacts detection. Indeed, though simple between two spheres -contact occurs when their centers are separated by a distance smaller than the sum of their radiicontacts detection between non-spherical particles is generally much more complex, hence more time-consuming. From the authors' experience, the calculation CPU time required for contacts detection between convex polyhedra generally lies in the range 10% to 20% of the total calculation CPU time. It may be lower for granular assemblies undergoing quasi-static strains, since new contact situations can generally be deduced from known contact situations at the previous time step. It may also reach higher values when simulating granular flows using highly faceted polyhedra.

The present paper focuses on initial contact detection between two overlapping convex polyhedra A and B having respectively v A and v B vertices, e A and e B edges, and f A and f B faces. Otherwise stated, we assume that the contact situation between these polyhedra has not previously been determined (e.g. during an earlier timestep). Initial contact detection is critical, since inappropriate detection may yield unacceptable overlap between the contacting particles that will be hard to eliminate during the simulation. Nevertheless it should be stressed that a small overlap between two contacting particles is often required by DEM methods in order to calculate repulsive contact forces (such a force is generally proportional to the overlap depth, e.g. see [START_REF] Wachs | Grains3d, a flexible dem approach for particles of arbitrary convex shape part i: Numerical model and validations[END_REF]). To avoid the lengthy testing of all interaction possibilitiesvertex/vertex, vertex/edge, vertex/ f ace, edge/edge, edge/ f ace, f ace/ f ace -, a number of algorithms has been designed and published, namely in the frameworks of rock engineering [START_REF] Cundall | Formulation of a three-dimensional distinct element modelpart i. a scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF][START_REF] Nezami | A fast contact detection algorithm for 3-d discrete element method[END_REF][START_REF] Boon | A new algorithm for contact detection between convex polygonal and polyhedral particles in the discrete element method[END_REF], robotics [START_REF] Lin | A fast algorithm for incremental distance calculation[END_REF] and computer graphics [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF][START_REF] Mirtich | V-clip: Fast and robust polyhedral collision detection[END_REF]. Unfortunately, these algorithms have drawbacks which more or less strongly affect their use into DEM simulations.

The objectives of the present paper are to address the ongoing drawbacks of three of the most widely used algorithms dedicated to initial contact detection between overlapping convex polyhedra, and then to identify the most accurate and efficient revised algorithm among them.

The paper is organized as follows: section 2 reviews initial contact detection algorithms and highlights their drawbacks while suggesting solutions when possible. Then section 3 focuses on one of the most widely used algorithms and introduces a revised version which addresses its drawbacks. Eventually, section 4 benchmarks this revised algorithm against two other revised algorithms over a set of overlapping contact situations, in order to rank them in terms of accuracy and calculation speed.

Review and discussion of contact detection algorithms

Pionneering work by Cundall [START_REF] Cundall | Formulation of a three-dimensional distinct element modelpart i. a scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF] has introduced the idea of a common plane (CP) separating two contacting polyhedra, which deflects to bisect the space between them. Basically, his idea is that the polyhedra are in contact with each other provided that each of them is in contact with the CP. Starting from the perpendicular bisector plane (PB) of the centroids of the particles (Fig 1), Cundall's algorithm consists in identifying the closest vertex of each particle to this plane, then translating the plane to the mid-point of the segment linking the identified closest vertices, and finally perturbating step by step the orientation of the unit normal of this plane until the gap between the closest vertex and the plane is maximum. Note that the gap may take negative values to account for situations in which particles overlap the CP -hence particles overlap -, so that maximizing the gap means minimizing the overlap and closest vertex means most deeply burried into the plane. Overall, Cundall's algorithm detects contacts between two particles and, in this case, yields the coordinates of contact point(s) and CP normal vector as well as gap width, all this with a reported complexity of order O(n iter .(v A + v B )), where n iter denotes the number of perturbation steps. In fact, a large perturbation step improves the algorithm complexity by decreasing the number of perturbation steps required to maximize the gap, but this is done at the expenses of the CP orientation accuracy and may hence yield significant unwanted interpenetration. More embarrassing, depending on the shape and relative location of the contacting particles, initializing Cundall's algorithm with the PB of the particle centroids may end up with an erroneous CP as shown on Fig. 2, insofar as the common plane fails to minimize the overlap between the particles. This problem is caused by the limited range of CP normal orientations explored when applying Cundall's perturbations, thus leading to a local maximum of the gap rather than the global one. To overcome this problem, a simple solution consists in repeating Cundall's algorithm three times, first upon initializing the algorithm with the PB of the centroids and then with two orthogonal planes which are perpendicular to the PB and intersect at the centroids mid-point. Among the calculated CP candidates, the one which minimizes the overlap should eventually be retained. Naturally, this solution triples the algorithm complexity from existing one. The revised Cundall algorithm is summarized as algorithm 1.

Nezami and coworkers [START_REF] Nezami | A fast contact detection algorithm for 3-d discrete element method[END_REF] have sped up Cundall's algorithm by observing that the number of common plane candidates is finite. Their algorithm begins similarly by (step 1) finding the perpendicular bisector plane of the centroids of the particles and (step 2) identifying the closest vertex of each particle to this plane. In case several vertices minimize the distance of a particle to the plane, the authors state that any of them may be chosen. Then (step 3), the CP is sought among the PB of the closest vertices and the planes parallel to one of the edges or faces of each particle which pass through its closest vertex, the latter planes being further translated to intersect at the midpoint of the closest vertices (see Fig. 3). The set of closest vertices of the particles to each CP candidate is identified, and, among all these, Figure 1: 3D-contact detection using the common plane strategy [START_REF] Cundall | Formulation of a three-dimensional distinct element modelpart i. a scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF]. The perpendicular bisector plane (PB) of the particle centroids C 1 and C 2 allows to identify the closest vertices V 1 and V 2 . The PB is further translated to the mid-point of the [V 1 ; V 2 ] segment and rotated until maximisation of the plane-particle gap (gap = ----→ V 1 V 2 . n) to achieve the common plane (CP). Translate the CP-candidate to mid-point of these vertices 7:

Initialize p 1 and q 1 so that ( n, p 1 , q 1 ) is an orthonormed frame 8:

for ( p, q) in {( p 1 , q 1 );(rot π/4 ( p 1 ),rot π/4 ( q 1 ))} do

n 1temp ← ( n + k p)/(1 + k 2 ) 12: n 2temp ← ( n -k p)/(1 + k 2 ) 13: n 3temp ← ( n + k q)/(1 + k 2 ) 14: n 4temp ← ( n -k q)/(1 + k 2 )
15:

for i = 1 to 4 do 16:

Find closest vertex of each particle to CP-candidate along n itemp 17:

Calculate interparticle gap i along n itemp end while 29: end for Figure 3: Two polyhedra and some of their common plane candidates (CP) according to [START_REF] Nezami | A fast contact detection algorithm for 3-d discrete element method[END_REF]. CP candidates are the perpendicular bissector plane (PB) of the closest vertices V 1 and V 2 of the polyhedra identified during the previous iteration as well as planes which are parallel to at least one edge (highlighted in bold red/dark color) of a polyhedron that contains either V 1 or V 2 . Note that each CP candidate passes through point M, the midpoint of V 1 and V 2 . Here, the actual common plane (CP) coincides with the PB of vertices V 1 and V 2 .

the CP candidate achieving the maximum distance between its set of closest vertices is retained. Finally (step 4), the retained set of closest vertices is compared to the previous one. If the two sets of closest vertices are equal, then the CP has been obtained, otherwise the algorithm iterates from step 3 using the last set of closest vertices and their CP candidate. Note that unlike Cundall's algorithm, Nezami's algorithm requires that no particle is in contact with the CP candidate prior to applying step 3. In case such a contact exists, the authors prescribe to translate both particles perpendicular to this plane and away from it to achieve a positive and small gap. Obviously, the opposite translation is applied to the particles once the CP has been identified. According to Nezami and coworkers, their algorithm yields the coordinates of contact points and common plane normal vector as well as overlap distance with a complexity of order smaller than O(v A + v B ). However, as shown in Fig. 4, the translation prescribed by Nezami in case of an existing contact between the CP candidate and a particle prior to applying his algorithm may yield an erroneous CP. In fact, such an error occurs when the most deeply penetrated vertex of a particle into the other particle changes with the translation direction. In Fig. 4, V 2 is the most deeply penetrated vertex into particle 1 along direction n 1 , hence it becomes the closest vertex to the CP candidate once the translation along n 1 is performed and CP 1 is found to be the common plane. Nevertheless, V 2 is the most deeply penetrated vertex into particle 1 along direction n 2 and this latter direction minimizes the overlap, hence CP 2 is the correct common plane and it should be observed that both CP candidates differ significantly. Steps 2, 3 and 4 of Nezami's algorithm only require minor modifications to appropriately account for overlapping particles without translating them. In case of a negative gap and as suggested by Cundall [START_REF] Cundall | Formulation of a three-dimensional distinct element modelpart i. a scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF], the closest vertices simply need to be redefined as the most deeply burried into the plane and the CP as the plane which minimizes the overlap, hence maximizes the algebraic distance between the sets of most deeply burried vertices (or equivalently minimize their absolute distance). The revised Nezami algorithm is summarized as algorithm 2.

Other algorithms are focused on finding the closest features (point, edge or face) between convex polyhedra A and B. Among these, the GJK algorithm [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF] was designed upon observing that finding the closest points between noncontacting polyhedra A and B is equivalent to finding the point of another convex set, consisting of their Minkowski difference A -B (A -B = {x Ax B , x A ∈ A, x A ∈ B}), which minimizes the distance to the coordinates origin O. Interestingly, the coordinates origin always lies outside the outer boundary of the Minkowsky difference A -B of noncontacting polyhedra, so that the point minimizing the distance of A -B to the coordinates origin may be redefined as the point V of the outer boundary of A -B located closest to the coordinates origin and || OV|| is the interparticle distance. When the polyhedra are in contact or overlap, the coordinates origin is located respectively on or inside the outer boundary of the Minkowski difference A -B, so that the point V of the outer boundary of A -B located closest to the coordinates origin defines both the overlap depth -|| OV|| (negative distance) and contact normal OV of the overlapping polyhedra (pointing from polyhedron A to polyhedron B).

To find the point V whilst avoiding the lengthy calculation of the entire Minkowski difference, the algorithm builds a sequence of simplices whose distance to the coordinates origin decreases monotonically to achieve the intended Figure 4: Common plane according to ref. [START_REF] Nezami | A fast contact detection algorithm for 3-d discrete element method[END_REF] (CP 1 , normal n 1 ) and correct common plane (CP 2 , normal n 2 ) between two interpenetrated polyhedra. n 2 is the direction of minimum relative translation to remove any interpenetration between the particles. V 1 and V 2 denote the closest vertices of particles 1 and 2 respectively to CP 1 , while V 1 and V 2 denote the closest vertices of these particles to the correct common plane (CP 2 ).

Algorithm 2 Revised Nezami algorithm

Require: Normal n 12 is known, with n 12 and mid-point of particle centroids defining their PB, CP-candidate initialized to PB 1: Find vertices V 1 and V 2 of particles located closest to CP-candidate along n 12 2: Calculate corresponding gap 3: gap max ← gap 4: Initialize n 12save so as to enter the loop 5: while n 12 n 12save do 6:

n 12save ← n 12 7: Determine mid-point MP of vertices V 1 and V 2 8: n 12 ← ---→ V 1 V 2 / ---→ V 1 V 2 9:
Find all faces F 1 and edges

E 1 of particle 1 comprising vertex V 1 10:
Find all faces F 2 and edges E 2 of particle 2 comprising vertex V 2 11:

for all couples of edges E 1 and E 2 do

12:

if E 1 non-parallel to E 2 then for all planes PB, F 1 , F 2 , P 12 , P 1 , P 2 do 23:

Find vertices V 1 and V 2 of particles located closest to plane along its normal n plane

24:

Calculate interparticle gap along n plane end for 31: end while minimum. Here, a simplex is a convex feature, e.g. a vertex, segment, triangle or tetrahedron, at most comprising r + 1 affinely independent vertices of the outer boundary of the Minkowski difference in a r-dimensional space. Given a simplex τ k at the end of the k-th iteration, the subsequent simplex τ k+1 is obtained by addition of a new vertex w k+1 so that τ k+1 is located closer to the coordinates origin than τ k , thus building iteration after iteration the closest simplex to the coordinates origin. Note that if τ k already contains r + 1 affinely independent vertices, then the vertex not belonging to the segment or the triangular face located closest to the coordinates origin respectively for r = 2 and r = 3 is deleted. Practically, w k+1 is calculated as the difference between the outermost vertex of A in the oriented direction v k , which minimizes the distance from previous simplex τ k to the origin and points towards the origin, and the outermost vertex of B in the opposite direction (Fig. 5). The GJK algorithm stops either when the origin lies inside the new simplex τ k+1 or in the absence of a new vertex w k+1 generating a closer simplex to the coordinates origin. Last, the coordinates of the point V k+1 of the final simplex τ k+1 located closest to the coordinates origin are calculated using the Johnson algorithm [START_REF] Johnson | The optimization of robot motion in the presence of obstacles[END_REF]. This algorithm recurses on the subsets of the final simplex, starting with those having the smallest number of vertices (e.g. vertices, then segments, then triangular faces, then tetrahedron) until one contains the point V k+1 . Upon calling S 1 to S p the p vertices of this subset, V k+1 is uniquely defined as the barycentre of vertices S 1 to S p which verifies the following equations:

p i=1 λ i = 1 and ∀ i ∈ {1, ...p} λ i > 0, (1)
p i=1 λ i OS i .( OS j -OS 1 ) = 0 for j ∈ {2, ..., p}, (2)
with OV k+1 = p i=1 λ i OS i .

Observe that equations ( 1) and ( 2) provide a small system of the form Mλ= L with p equations (p ≤ 4) and p unknowns λ i , which may be inverted to determine the λ i values and calculate both the distance || OV k+1 || of V k+1 to the origin and the orientation of OV k+1 . Furthermore, note that if the particles are not in contact, V k+1 is also the closest point of the Minkowsky difference A -B to the origin and || OV k+1 || the distance between the two polyhedra. Conversely, when the polyhedra are in contact or overlap, || OV k+1 || denotes the overlap depth and OV k+1 defines the contact normal.

According to Gilbert and coworkers [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF], their algorithm has a complexity of order O(v A + v B ).

Despite its widespread use due to its low computational costs and high versatility, the GJK algorithm was reported to suffer numerical unstabilities resulting from rounding errors [START_REF] Van Den Bergen | A fast and robust gjk implementation for collision detection of convex objects[END_REF][START_REF] Montanari | Improving the gjk algorithm for faster and more reliable distance queries between convex objects[END_REF]. These rounding errors would either cause the simplex generation algorithm to loop forever [START_REF] Van Den Bergen | A fast and robust gjk implementation for collision detection of convex objects[END_REF] or the Johnson algorithm to be numerically unstable as the combination of dot products and differences may yield ill-conditionned system of equations ( 1) and (2). As an example, ill-conditioning would occur when the simplex generation algorithm converges towards a degenerate simplex, in other words a simplex comprising vertices which are nearly affinely dependent such as an elongated triangular face having two vertices lying very close to each other. This situation may be encountered when a large flat face of a polyhedron interacts with another polyhedron smaller by several orders of magnitude [START_REF] Van Den Bergen | A fast and robust gjk implementation for collision detection of convex objects[END_REF], a situation that may be somewhat frequent in DEM simulations of convex polyhedra assemblies. Solutions to circumvent GJK issues have been suggested and tested with reasonable success by several authors [START_REF] Van Den Bergen | A fast and robust gjk implementation for collision detection of convex objects[END_REF][START_REF] Ericson | Real-Time Collision Detection[END_REF][START_REF] Tereshchenko | Algorithm for finding the domain intersection of a set of polytopes[END_REF][START_REF] Montanari | Improving the gjk algorithm for faster and more reliable distance queries between convex objects[END_REF], starting with Gilbert et al themselves [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF] who designed a backup procedure to cache the above-mentioned dot products which would otherwise be calculated several times by the Johnson algorithm with the risk of rounding errors being introduced. This backup procedure was further improved by Van den Bergen [START_REF] Van Den Bergen | A fast and robust gjk implementation for collision detection of convex objects[END_REF] who suggested a caching algorithm which minimizes the caching overhead. This author also sped up the GJK algorithm by designing an early termination procedure to detect non-contacting particles, focused on detecting the existence of a separating axis between them. More recently, Montanari et al [START_REF] Montanari | Improving the gjk algorithm for faster and more reliable distance queries between convex objects[END_REF] have suggested the signed volume method to circumvent the numerical instabilities of the Johnson algorithm. The authors' idea is to disembed the determination of the simplex closest to the coordinates origin from the calculation of the closest point coordinates. Once the GJK algorithm has converged to the simplex located closest to the coordinates origin, their method consists in (step 1) identifying the smallest set of vertices of this simplex which contains the closest point to the coordinates origin and (step 2) determining the barycentric coordinates of this closest point by solving a small system of equations comprising equation ( 1) and a revised form of equation ( 2) which writes:

p i=1 λ i OS i = P 0 , (3)
where P 0 contains the coordinates of the projection of O on the affine-hull of the vertices S i followed by 1. Again, equations ( 1) and (3) provide a small system of the form Mλ= P 0 , though this time the i-th column of M no more comprises scalar products likely to cause numerical instabilities when their values tend to zero but simply the coordinates of vector OS i followed by 1. When the coordinates origin is inside the simplex, meaning that the particles overlap, the authors state that λ i values may be determined by the following equation:

λ i = C 4,i /det(M), (4)
where det denotes the determinant of M and C 4,i its i-th cofactor. Observe that equation ( 4) is erroneous since, in accordance with Cramer's rules, C 4,i cannot be the minor of M in absolute value but should comprise the coordinates of vector P 0 as its i-th column to read:

λ i = det[ M 1 ... M i-1 P 0 M i+1 ... M p ]/det(M), (5)
where M j denotes the j-th column of matrix M.

The Lin-Canny algorithm [START_REF] Lin | A fast algorithm for incremental distance calculation[END_REF] uses another strategy to find the closest features, which may be summarized as follows: define criteria (step 1) to partition the space surrounding each convex polyhedron into regions, each region being closer from a given feature than from any other feature making up the surface of the polyhedron. Then (step 2), select two features, one on each polyhedron, and determine their closest points. Next (step 3), check whether the closest point on one feature belongs to the region of the other feature and conversely. If either check fails, then (step 4) by identifying which criteria failed, select a new set of features on the polyhedra which more likely yields successful checks and iterate from step 3. According to the authors, this algorithm terminates in order O(v A + v B ) iterations and yields the distance and closest features between two polyhedra. Unfortunately, this algorithm gets stuck in a loop when the two polyhedra are in contact and forcing its termination relies on arbitrary iteration threshold value yet does not allow determining an interpenetration depth if any.

More recently, algorithms based on convex optimization have been designed [START_REF] Boon | A new algorithm for contact detection between convex polygonal and polyhedral particles in the discrete element method[END_REF][START_REF] Boon | A new contact detection algorithm for three-dimensional non-spherical particles[END_REF]. Basically, such an algorithm consists in (step 1) writing linear inequalities which define the interior of each polyhedron, then (step 2) establishing whether there exists at least one point verifying simultaneously all the previous inequalities. Next (step 3), in case such a point exists, meaning that both polyhedra intersect, the contact point -where the interparticle contact forces apply -is determined as the analytical center of the intersection area using an adapted version of the logbarrier method [START_REF] Boyd | Convex Optimization[END_REF]. This method expresses the burrial depth -with reference to each plane defining the intersection area -of any point located into this area and returns the point maximizing the sum of burrial depth log-function. Last (step 4), the contact normal is calculated as the weighed average of the gradients of two polynomial inner potentiels, one for each polyhedron, each potential being defined inside one of the particles as a positive and decreasing function of the distance to the particle surface.

Since GJK appears to be the most versatile contact detection algorithms -unlike others, it may easily handle polyhedron/spheroid contact situations -, the next section sheds a new light on its numerical instabilities and introduces a more robust GJK algorithm which fixes them.

New light on GJK drawbacks and revised algorithm

New light on GJK drawbacks

The contact situation depicted in Fig. 5 is of particular interest to shed a new light on the numerical instabilities of GJK algorithm. In this case, the GJK algorithm converges to a simplex whose face located closest to the coordinates origin belongs to the convex hull of the Minkowsky difference of particles A and B (see Fig. 5e). However, if Particle B is slightly translated along the -X and -Z axes with reference to Particle A as shown on Fig. 6a, then the GJK algorithm will converge to the simplex depicted in Fig. 6b. Observe that this tetraedral-shaped simplex has two interesting particularities: first, unlike its vertices, none of its triangular faces is located on the convex hull of the Minkowsky difference; second, the coordinates origin is located strictly inside this tetrahedron. As a consequence, applying the Johnson algorithm to determine the closest point of this simplex to the coordinates origin would lead to point V A -V B (Johnson) (see zoomed Fig. 6c). Unfortunately, this point is burried into the Minkowsky difference of particles A and B and not located on its convex hull, hence the vector linking the coordinates origin to this point cannot be representative of the overlap direction or distance between the two particles. The true closest point V A -V B (correct) of the Minkowsky difference of particles A and B to the coordinates origin, also depicted in Fig. 6b andc, differs significantly from the previous one. In fact, it is not unusual that the GJK algorithm converges to a simplex such as the one depicted in Fig. 6b, with no face located on the convex hull of the Minkowsky difference of the two particles. Furthermore, the response returned by the GJK algorithm may be even more unsatisfactory. Indeed, upon slightly translating particle B of Fig. 5a along -X as depicted in Fig. 8a, the GJK algorithm will converge to the new tetrahedral-shaped simplex depicted in Fig. 8b. Further to having none of its triangular faces located on the convex hull of the Minkowsky difference of the particles, the particularity of this simplex is that the coordinates origin lays on one of its faces. As a consequence, the closest point of this simplex to the coordinates origin as determined by the Johnson algorithm will be the origin itself, hence the distance between the particles will erroneously be found equal to zero and no overlap will be detected.

A more robust algorithm: GJK-TD

Observe that, whether polyhedra A and B overlap or not, the point V located closest to the coordinates origin always lies on the outer boundary of their Minkowsky difference A -B (see section 2). Furthermore, for overlapping polyhedra, GJK always converges to a tetrahedral-shaped simplex and the point of this simplex located closest to the coordinates origin always lies on one of its triangular faces (e.g. see points V A -V B (Johnson) in Fig. 5e, 6c and8b). A straightforward consequence is that, in order to circumvent the drawbacks of the GJK algorithm evidenced in sec-tion 3.1, additional steps are needed to ensure that, at the end of the GJK algorithm, the simplex triangular face located closest to the coordinates origin lies on the outer boundary of the Minkowsky difference A -B of the polyhedra. These steps are depicted in the next paragraph using the notations of section 3.1.

Consider the contact situation depicted in Fig. 8a between Particle A and Particle B (shi f ted), for which applying the GJK algorithm yields the tetrahedral-shaped simplex of Fig. 8b. The basic idea of the steps taken to relocate the closest triangular face to the coordinates origin on the outer boundary of the Minkowsky difference A -B of the polyhedra may be summarized as follows : 1) determine the distance of each triangular face of the tetrahedral-shaped simplex to the coordinates origin, then 2) expand the located closest to the coordinates origin along its normal towards the outer boundary of the Minkowsky difference A -B of the polyhedra, and 3) repeat this expansion with the subsequent closest triangular face among the resulting set of triangular faces until the condition that the closest triangular face is located on the outer boundary of the Minkowsky difference A -B of the polyhedra is fulfilled. Both the expansion principle and face location condition may simply be explained using the contact situation depicted in Fig. 8.

To describe the expansion principle, consider the tetrahedral-shaped simplex (A 1 , B 1 , C 1 , D 1 ) of Fig. 9a, which is the same simplex as in Fig. 8b but observed from another perspective. Once the triangular face (A 1 , B 1 , C 1 ) minimizing the distance d 1 to the coordinates origin has been identified, its normal n 1 pointing outside the initial simplex (or, equivalently, in the direction opposite to the coordinates origin) is calculated. Next, vertex E 1 is determined as the Minkowsky difference between the outermost vertex of particle A in the oriented direction n 1 and the outermost vertex of particle B in the opposite direction (see Fig. 9b). Last, triangular faces (A 1 , B 1 , E 1 ), (B 1 , C 1 , E 1 ) and (A 1 , C 1 , E 1 ) are substituted for triangular face (A 1 , B 1 , C 1 ) in the set of faces in which the closest face to the coordinates origin is sought. Fig. 9c to e depict the application of the expansion principle to the subsequent face located closest to the coordinates origin and so on. It should be emphasized that, until the closest triangular face to the coordinates origin is found on the Minkowsky outer boundary of the polyhedra (see triangular face (B 1 , E 4 , E 3 ) in Fig. 9e), each triangular face is deleted right after being visited (e.g. face (A 1 , B 1 , C 1 )). Furthermore, records are kept to ensure that, during the subsequent iterations of the expansion process, deleted faces are never rebuilt, thus ensuring that the algorithm will not get trapped in a non-ending loop.

The face location condition used to ensure that the current triangular face is located on the outer boundary of the Minkowsky difference A -B of the polyhedra simply consists in 1) determining the outermost vertex of particle A in the oriented direction n normal to the current face (see Fig. 9e) and the outermost vertex of particle B in the opposite direction, and 2) checking that the resulting point E 4 , which is a vertex of the outer boundary of the Minkowsky difference A -B of the polyhedra, belongs to the current triangular face.

Once the triangular face located closest to the coordinates origin (B 1 , E 4 , E 3 ) has been determined on the Minkowsky outer boundary of the polyhedra (see Fig. 9e), the closest point V of this face to the coordinates origin O is calculated as the barycentre of vertices (B 1 , E 4 , E 3 ) for which vector OV is normal to the face:

3 i=1 λ i = 1 and ∀ i ∈ {1, ..., 3} λ i ≥ 0, (6) 3 i=1
λ i OS i = 0V, (7)
with

S i in {B 1 , E 4 , E 3 } and 0V = ( OS 1 . n) n.
Equations ( 6) and ( 7) provide a small system of the form Mλ= OV, whose unknowns λ i are easily determined by the following equations: Similarly, the last step (e) identifies (B 1 , C 1 , E 3 ) as the closest triangular face to the coordinates origin, with n 4 its outward normal allowing to identify outermost vertex E 4 = E 2 . This new vertex yields triangular face (B 1 , E 4 , E 3 ) which is both located closest to the coordinates origin and on the convex hull of the Minkowsky difference (represented in light-grey), hence the expansion is stopped.

λ 1 = det[ OV OS 2 OS 3 ]/det[ OS 1 OS 2 OS 3 ], (8)
λ 2 = det[ OS 1 OV OS 3 ]/det[ OS 1 OS 2 OS 3 ], (9)
λ 3 = det[ OS 1 OS 2 OV]/det[ OS 1 OS 2 OS 3 ]. (10)
Last, these barycentric coordinates λ i are used to determine the closest points V c A and V c B of polyhedra A and B as the barycentres of the particle vertices forming the Minkowski vertices {B 1 , E 4 , E 3 } (see Fig. 8a). As already mentioned, -|| OV|| (negative distance) defines the overlap depth of the overlapping particles and OV their contact normal (pointing from particle A to particle B).

The new algorithm that has just been presented is summarized as Algorithm 3. It was called GJK -T D to emphasize the face expansion steps (T D stand for T etrahedron Distortion) performed at the end of a basic application of the GJK algorithm. Several comments shall be made about this new algorithm when applied to two overlapping convex polyhedra made of triangular faces. First, like the Signed Volumes method [START_REF] Montanari | Improving the gjk algorithm for faster and more reliable distance queries between convex objects[END_REF], the matrix M of the small system defined by equations ( 6) and ( 7) involves no scalar product likely to cause ill-conditioning of the system. Hence, λ i values defined by equations ( 8), ( 9), [START_REF] Zhao | Shape Effects on the Random-Packing Density of Tetrahedral Particles[END_REF] may always be calculated. Second, the expansion algorithm always converges to a triangular face located on the outer boundary of the Minkowsky difference between the polyhedra, and the convergence occurs in less than 3(v A + v B ) -8 steps. Indeed, the algorithm Algorithm 3 GJK-TD algorithm Require: GJK applied to convex polyhedra A and B has converged to a tetrahedral-shaped simplex for each new triangle made of E q and two vertices (S r1 , S r2 ) among (S i , S j , S k ) do 21:

if Triangular face (S r1 , S r2 , E q ) not in ( visited U ) then

22:

Calculate normal n q+4 and distance d q+4

23:

← n q+4 , d q+4 and face connectivity 24:

q ← q + 1 25:

end if end if 30: end loop consists in finding a sequence of triangular faces whose distance to the coordinates origin is strictly increasing. Furthermore, the vertices of these triangular faces belong to the Minkowsky difference of the polyhedra. Hence, the distances to the coordinates origin of the sequence of triangular faces are bounded by that of the closest triangular face of the outer boundary of the Minkowsky difference. As a consequence, like for any strictly increasing and bounded mathematical sequence in a convex set, the convergence is always achieved. Next, each convex polyhedron consisting of triangular faces, the outer boundary of their Minkowsky difference is also made of triangular faces. Besides, the triangular face to which the expansion algorithm converges is always located on the outer boundary of the Minkowsky difference of the polyhedra as a consequence of the face location condition. Finally, the complexity of the algorithm may be assessed as follows:

• According to the Descartes-Euler theorem, the total number of faces f , edges e and vertices v of any convex polyhedron relate by the formula fe + v = 2;

• Now, the outer boundary of the Minkowsky difference of the polyhedra is made exclusively of triangular faces, each comprised of three edges shared by exactly two faces, hence the total number of edges and faces relate according to e = 3 f 2 ; • The two preceeding relations yield f = 2v -4;

• Besides, the expansion algorithm may visit internal triangular faces, that is faces embedded into the Minkowsky difference of the polyhedra and not located on its outer boundary (the latter being external faces), whose total number is f -4 (this may easily be shown by recurrence);

• Last, the outer boundary of the Minkowsky difference between two polyhedra made of v A and v B vertices respectively contains exactly v A + v B vertices. As a consequence, the expansion algorithm will visit at most all 2(v A + v B ) -4 external and v A + v B -4 internal triangular faces, and these faces will be visited at most once since any visited face is deleted and never restored until convergence.

Efficiency of the GJK-TD algorithm

In order to examine the efficiency of the GJK-TD algorithm, one shall address both its calculation accuracy and speed. This is achieved upon considering nine contact situations with various overlap, which are used as a benchmark to compare GJK-TD with the revised versions of Cundall's and Nezami's algorithms as depicted in section 2. Given that these three algorithms implement significantly different strategies to determine a contact situation, it seems reasonable to consider that the level of agreement between the results returned by these methods reflects their accuracy. Furthermore, upon selecting more or less overlapped contact situations, which are representative of all common contact types, namely vertex-face, edge-face, edge-edge, or face-face contacts, comparing the calculation speed between these three algorithms should yield representative orders of magnitude. The nine contact situations used are depicted in Fig. 2,4, 5f, 7, 8a, 10a, 10b, 10c, 10d respectively.

All calculations are performed on a double precision Ubuntu 16.04 -LT S machine equipped with an intel R Core TM i7 -6700HQ CPU@2.60GHz processor having 8 GB RAM. The calculation routines, GJK-TD as well as the revised versions of Cundall's and Nezami's algorithms, are written in Fortran90 and compiled using GNU Fortran 5.4.0.

Computing accuracy

For each of the nine contact situations, both the normal components and overlap depth are calculated to the highest possible precision level. All calculation routines use double precision values, hence results are returned with 15 significant digits (machine precision). For each contact situation, Table 1 summarizes the results returned by the GJK-TD algorithm and Fig. 11 displays the accuracy with which the revised versions of Cundall's and Nezami's algorithms agree with the results returned by GJK-TD. Note that two separate charts are displayed, Fig. 11a depicting accuracy in terms of mean absolute difference between the normal components calculated by any of the revised algorithm and GJK-TD, and Fig. 11b depicting accuracy in terms of overlap depths absolute difference with reference to GJK-TD. Observe that, whatever the contact situation, GJK-TD and the revised version of Nezami's algorithm return results which agree with eachother to machine precision. This is consistent with the finite and limited number of possible common plane candidates when two convex polyhedra overlap. By contrast, the revised version of Cundall's algorithm returns results which agree with those of the other two algorithms at best to 10 -10 . Note that a better agreement could be achieved upon setting the Cundall's minimum orientational perturbation step to a value below the requested 10 -7 π/180 radians, but this would strongly affect the calculation speed.

Next, the level of accuracy of these results is compared with that of the results reported by Montanari et al. [START_REF] Montanari | Improving the gjk algorithm for faster and more reliable distance queries between convex objects[END_REF] on gear teeth with various number of vertices (see their figure 15). Montanari's results appear to be at least one order of magnitude less accurate than the results reported in the present paper.

Computing speed

For each of the nine contact situations depicted in Fig. 2,4, 5f, 7, 8a, 10a, 10b, 10c, 10d, Table 2 summarizes the calculation speed of all three contact detection algorithms both in terms of absolute and relative to GJK-TD values. These calculation speed figures are obtained using the S YS T EM CLOCK Fortran function, which allows a level of precision up to 1 nanosecond, and the reported values are averaged over one million repetitions of each algorithm. According to this table, GJK-TD is able to return the contact normal and overlap depth within 6 to 95 microseconds, depending on the contact situation, whereas the revised versions of Cundall's and Nezami's algorithms are respectively 6 and 65 times more computationally intense on average, and up to 12 and 135 times slower respectively. Furthermore, it should be pointed out that, due to the extra loop encompassing lines 4 to 28 of Algorithm 1, the revised Cundall algorithm has a complexity which is three times that of the classical Cundall algorithm. In other words, the computational cost of the former is three times that of the latter. Yet, as mentioned in section 2, in some configurations, the classical Cundall algorithm fails to return correct results for the contact plane.

A first attempt was made to confront these figures with those reported by Montanari et al. [START_REF] Montanari | Improving the gjk algorithm for faster and more reliable distance queries between convex objects[END_REF] on overlapping polygonal spheres with various number of vertices tested over one million cycles (see their figure 17c). Their figure 17c suggests that the speed of the tested algorithms does not significantly vary with the number of vertices, hence confronting their results in terms of algorithm velocity to those of the present paper obtained with different polyhedralshaped particles makes sense. However, their figure raises questions about the CPU time unit (ns): indeed, repteating a simple assignment to 1 of a unique integer variable for one million cycles already takes 6.153152 10 -3 s with the machine used in the present paper, hence about 6 ns per assignement, and about 2.5 ns per assignment on an Intel R Xeon R CPU E5 -2630 v4 2.20GHz processor having 28 GB RAM. As a consequence, it seems doubtfull that the run of a full algorithm such as JB, BK or S V as reported by Montanari et al. only takes a few nanoseconds. Then we confronted our contact detection speed figures to those reported by Wachs et al. [START_REF] Wachs | Grains3d, a flexible dem approach for particles of arbitrary convex shape part i: Numerical model and validations[END_REF] on overlapping cubes or tetrahedra. These authors have implemented a tailored version of the GJK algorithm in which overlapping particles are first reduced through a homothety, so that they are no more in contact, prior to determining the set of closest points. The contact detection speed reported by these authors on an Intel4 R Pentium CPU 3GHz processor is 74 µs and 106 µs for overlapping cubes and tetrahedra respectively. These figures are consistent with those of the present paper.

Conclusion

Algorithms dedicated to initial contact detection between overlapping convex polyhedra have been reviewed, focusing on three of the most prominent ones [START_REF] Cundall | Formulation of a three-dimensional distinct element modelpart i. a scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF][START_REF] Nezami | A fast contact detection algorithm for 3-d discrete element method[END_REF][START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF]. Drawbacks affecting the use of these three algorithms into DEM simulations have been evidenced and solutions to these drawbacks have been suggested. In particular, new light was shed on GJK drawbacks and a new algorithm supplementing the original one was introduced. Finally, revised versions of these three algorithms implementing the suggested solutions have been benchmarked over nine contact situations. These contact situations involve two overlapping polyhedra which are representative of all common contact types, namely vertex-face, edge-face, edge-edge and face-face contacts. The benchmarking results show that GJK-TD and Nezami (revised) return identical results to machine precision, whereas Cundall (revised) returns identical results to the tenth decimal place. Furthermore, these results show the GJK-TD returns both the normal components and overlap depth within a few tens of microseconds, whereas Nezami (revised) and Cundall (revised) are respectively 6 and 65 times more computationally intense on average. It is believed that the robustness and efficiency of GJK-TD will boost its use into DEM simulations, all the more that this versatile algorithm may easily be customized to detect

Figure 2 :

 2 Figure 2: 3D-common plane candidates CP 1 (Cundall) and CP 2 (correct) between two parallelepiped interpenetrated by an edge. Interpenetration distance d 1 to CP 1 is significantly larger than interpenetration distance d 2 to the correct common plane (see inset). Note that corresponding normals n 1 and n 2 are also significantly different.

4 :

 4 while closest vertices to CP-candidate change do 5:Find closest vertex of each particle to CP-candidate along n 6:

Figure 5 :

 5 Figure 5: Principle of the GJK algorithm [20] demonstrated using (a) Overlapping parallelepiped and prism generalized from 2D Fig. 7a of ref. [23] (resp. particle A with centroid C 1 and eight cross-marked vertices A 1 to I 1 , and particle B with centroid C 2 and six cross-marked vertices A 2 to G 2 ), with A 1 the outermost vertex of particle A in the oriented direction v 1 and D 2 the outermost vertex of particle B alongv 1 . (b) The difference of preceeding vertices A 1 -D 2 belongs to the convex hull of the Minkowsky difference of particles A and B and defines initial simplex τ 1 as well as new oriented direction v 2 pointing from A 1 -D 2 to coordinates origin O, along which outermost points of particles A and B define new point D 1 -A 2 of the convex hull (fully represented here in light-grey together with crosses representing its vertices). (c) v 3 is the oriented direction minimizing the distance between new simplex τ 2 containing two points A 1 -D 2 and D 1 -A 2 and coordinates origin, along which outermost points of particles A and B define new point E 1 -B 2 of the convex hull. (d) Following the same principle, v 4 is the oriented direction minimizing the distance between new simplex τ 3 containing three points A 1 -D 2 , D 1 -A 2 and E 1 -B 2 , and coordinates origin, and B 1 -G 2 is the new point identified on the convex hull of the Minkowski difference. (e) Following the same principle, the gjk algorithm has converged to the tetrahedra (D 1 -A 2 , E 1 -B 2 ,H 1 -B 2 ,D 1 -F 2 ) which encloses the coordinates origin and V k+1 = V A -V B the closest point of the Minkowsky difference A -B to the origin. (f) Overlapping polyhedra and corresponding common plane (CP) and normal n calculated using the GJK algorithm, with V A and V B the closest points identified between particles A and B ( V B V A and n are parallel).

Fig. 7

 7 reflects this difference in terms of common plane orientation and overlap distance between particles A and B.

Figure 6 :

 6 Figure 6: Particles shift (a) and new simplex resulting from the GJK algorithm (c) is a magnified view of (b). The shift applied is a simple translation of particle B along the -X and -Z axes. (b) and (c) display the closest point V A -V B (Johnson) of the tetrahedral-shaped simplex to the coordinates origin as calculated by the Johnson algorithm, the true closest point of the Minkowsky difference of particles A and B to the coordinates origin V A -V B (correct), as well as the corresponding overlap distances d 1 and d 2 respectively.

Figure 7 :Figure 8 :

 78 Figure 7: Common plane orientations (CP), closest points V β α and overlap distances d β between particles A and B, superscript J refers to Johnson's algorithm and c to true values.

Figure 9 :

 9 Figure 9: Tetrahedron expansion principle. Starting from tetrahedral-shaped simplex (A 1 , B 1 , C 1 , D 1 ) (a), the closest triangular face (A 1 , B 1 , C 1 ) to the coordinates origin is identified and its outward normal n 1 (oriented opposite to the coordinates origin) is calculated. Then, the outermost vertex E 1 along n 1 is identified (b), triangular face (A 1 , B 1 , C 1 ) is deleted and replaced by triangular faces (A 1 , B 1 , E 1 ), (B 1 , C 1 , E 1 ), (A 1 , C 1 , E 1 ) so that the former tetrahedral-shaped simplex becomes the new polyhedron (A 1 , D 1 , B 1 , C 1 , E 1 ). Similarly, the closest triangular face (D 1 , B 1 , C 1 ) to the coordinates origin of the new polyhedron is identified (c) and its outward normal n 2 is calculated, in order to determine vertex E 2 and substitute triangular faces (B 1 , C 1 , E 2 ), (C 1 , D 1 , E 2 ) and (D 1 , B 1 , E 2 ) for triangular face (B 1 , C 1 , D 1 ) to achieve a new polyhedron. Following the same principle, the closest triangular face (E 1 , B 1 , C 1 ) to the coordinates origin of the new polyhedron is identified (d) and its outward normal n 3 is calculated, in order to determine vertex E 3 and substitute triangular faces (B 1 , C 1 , E 3 ), (C 1 , E 1 , E 3 ) and (E 1 , B 1 , E 3 ) for triangular face (B 1 , C 1 , E 1 ).Similarly, the last step (e) identifies (B 1 , C 1 , E 3 ) as the closest triangular face to the coordinates origin, with n 4 its outward normal allowing to identify outermost vertex E 4 = E 2 . This new vertex yields triangular face (B 1 , E 4 , E 3 ) which is both located closest to the coordinates origin and on the convex hull of the Minkowsky difference (represented in light-grey), hence the expansion is stopped.

2 : 3 : 9 :

 239 for face p in simplex do Calculate d p , n p 4: ← n p , d p , face connectivity 5: end for 6: Search for triangular face (S i , S j , S k ) with smallest d p 10: if Triangular face location condition = TRUE then 11: Calculate OV = ( OS i . n p ) n p 12:Calculate λ i , λ j and λ k

Figure 10 :

 10 Figure 10: Various 3D contact situations used to benchmark GJK-TD with the revised versions of Cundall's and Nezami's algorithms, (a) vertexedge contact, (b) edge-edge contact, (c) edge-face contact and (d) face-face contact. n is the normal vector of the common plane (CP), pointing from Particle A to Particle B.

Figure 11 :

 11 Figure 11: Comparison of calculation accuracy against GJK-TD for algorithms Cundall (revised) and Nezami (revised), for nine contact situations, (a) normal mean components, (b) overlap, (c) edge-face contact and (d) face-face contact.

  Algorithm 1 Revised Cundall algorithm Require: Orthonormed frame ( n 12 , p 12 , q 12 ) is known, with n 12 and mid-point of particle centroids defining their PB, CP-candidate initialized to PB, gap max and orientation perturbation step range initialized 1: gap max ← -1000. 2: [k min ; k max ] ← [tan(10 -7 );tan(5 • )] 3: for n in { n 12 , p 12 , q 12 } do

Table 1 :

 1 Normal components and overlap calculated to machine precision(10 -15 ) by GJK-TD for nine contact situations.

	Test-case	n x	n y	n z	overlap
	Fig. 2	0.000000000000000	1.000000000000000	0.000000000000000	-0.250000000000000
	Fig. 4	0.992277876713668	-0.124034734589208	0.000000000000000	-4.026012719933860
	Fig. 5f	0.894427190999916	-0.447213595499958	0.000000000000000	-0.223606797749979
	Fig. 7	0.894427190999916	-0.447213595499958	0.000000000000000	-1.118033988749900
	Fig. 8a	0.894427190999916	-0.447213595499958	0.000000000000000	-1.118033988749900
	Fig. 10a	0.000000000000000	0.000000000000000	-1.000000000000000 -0.100000024000000
	Fig. 10b	-0.168309951548722	0.950616325200044	-0.260768791217108 -0.572082006256242
	Fig. 10c	0.832050294337844	-0.554700196225229	0.000000000000000	-0.554700196225229
	Fig. 10d	0.970142500145332	0.000000000000000	-0.242535625036333 -0.485071250072666

Table 2 :

 2 Comparison of calculation speed between algorithms Cundall (revised), Nezami (revised) and GJK-TD, for nine contact situations.Test case GJK-TD Nezami (revised) Nezami (revised)/ Cundall (revised) Cundall (revised)/

		(x10 -6 s)	(x10 -6 s)	GJK-TD	(x10 -6 s)	GJK-TD
	Fig. 2	10.896	83.920	7.7	540.656	49.6
	Fig. 4	94.688	121.456	1.3	1225.261	12.9
	Fig. 5f	6.425	77.768	12.1	870.992	135.6
	Fig. 7	19.095	91.640	4.8	938.476	49.1
	Fig. 8a	18.975	77.156	4.1	671.476	35.4
	Fig. 10a	10.697	67.356	6.3	1288.066	120.4
	Fig. 10b	27.189	85.672	3.2	2582.582	95.0
	Fig. 10c	22.239	161.700	7.3	1056.096	47.5
	Fig. 10d	19.208	85.068	4.4	719.848	37.5
	mean	25.490	94.637	5.7	1099.273	64.8

contact between convex polyhedra and spheroid particles.

References References


Mariana Ballesteros
email: mballesteros@ctrl.cinvestav.mx

Andrey Polyakov
email: andrey.polyakov@inria.fr

Denis Efimov
email: denis.efimov@inria.fr

Isaac Chairez

Alexander Poznyak

Differential Neural Network Identification for Homogeneous Dynamical Systems ⋆

Keywords: Differential Neural Network, Nonlinear Systems, Homogeneous systems, Identification

In this paper, a non parametric identifier for homogeneous nonlinear systems affine in the input is proposed. The identification algorithm is based on the neural networks using sigmoidal activation functions. The learning algorithm is derived by means of Lyapunov function method and homogeneity theory. A numerical example demonstrates the performance of the proposed identifier.

INTRODUCTION

Analysis and design of dynamic control systems need a valid mathematical model of the plant. However, most control systems have model uncertainties. Therefore, the tools looking for a valid mathematical description of the dynamic system are widely demanded and applied.

Nonlinear system identification is a field of control theory [START_REF] Ljung | Some aspects on nonlinear system identification[END_REF], which develops algorithms of mathematical modeling of control systems based on input and output signals measured on-line or/and during some experiments. Identification problem has been tackled using different approaches due to a large class of system models and its inherent complexity. The identification process in general needs to handle the available input and output data in order to postulate a model and validate it somehow. Some of the most popular identification techniques are: functional series methods, frequency domain approaches, fuzzy models and neural networks (NNs) (see e.g. [START_REF] Billings | Identification of nonlinear systems-a survey[END_REF], [START_REF] Haber | Nonlinear system identification-Input-output modeling approach[END_REF], [START_REF] Nelles | Nonlinear system identification: from classical approaches to neural networks and fuzzy models[END_REF] and [START_REF] Haykin | Neural networks: a comprehensive foundation[END_REF] for more details).

NNs are a special kind of approximation algorithms having the theoretical capability to approximate a large class of nonlinear mappings, (see [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF]). Differential neural networks (DNN) are utilized for approximation of dynamical systems [START_REF] Poznyak | Differential neural networks for robust nonlinear control: identification, state estimation and trajectory tracking[END_REF], [START_REF] Lewis | Neural network control of robot manipulators and non-linear systems[END_REF], [START_REF] Sontag | Some topics in neural networks and control[END_REF], since they can be trained on-line (in a real time). DNNs can process many inputs and outputs, so they are applicable to multi-variable systems.

⋆ This work was partially supported by the Government of Russian Federation (Grant 08-08) and the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031). Mariana Ballesteros is sponsored by Mexico scholarship from the National Council for Science and Technology, (Scholar reference: 550803) The DNN identification admits the selection of different activation functions, which represent a certain basis for the model of the system in an admissible space, for example, sigmoidal functions, polynomials or radial basis functions [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF], [START_REF] Diaconis | On nonlinear functions of linear combinations[END_REF], [START_REF] Hubbert | Radial basis function interpolation on the sphere[END_REF]. The adjustment of the time-varying parameters (weights) in the DNN according its structure and the set of activation functions should be adjusted, for example, by a stability analysis based on Lyapunov procedure (see e.g. [START_REF] Chairez | Adaptive neural network nonparametric identifier with normalized learning laws[END_REF], [START_REF] Jagannathan | Identification of nonlinear dynamical systems using multilayered neural networks[END_REF]). In this work the DNN identification algorithms are developed for a specific class of systems: homogeneous nonlinear systems affine in the input.

Homogeneity is a symmetry-like property under which an object remains consistent with respect to a certain scaling or dilation. Homogeneous systems can be utilized for local approximations [START_REF] Hermes | Nilpotent approximations of control systems and distributions[END_REF], [START_REF] Andrieu | Homogeneous approximation, recursive observer design, and output feedback[END_REF] or setvalued extensions [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF], [START_REF] Levant | Homogeneity approach to high-order sliding mode design[END_REF] of nonlinear control systems. In particular, some models of process control [START_REF] Zimenko | A note on delay robustness for homogeneous systems with negative degree[END_REF], nonholonomic mechanical systems [START_REF] Pomet | Time-varying exponential stabilization of nonholonomic systems in power form[END_REF] and systems with frictions [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF] are homogeneous or at least locally homogeneous.

To the best of authors knowledge, identification problem of homogeneous systems is not well studied in the literature. One of the main features of homogeneous system is that an analysis of its behavior in a whole state space can be reduced to a similar analysis on a unit sphere [START_REF] Hermes | Homogeneous feedback controls for homogeneous systems[END_REF], [START_REF] Bernuau | On homogeneity and its application in sliding mode control[END_REF], [START_REF] Polyakov | On homogeneous distributed parameter systems[END_REF]. This feature implies a specific structure of the DNN identifier in this work. The activation functions are selected to approximate the systems on the unit sphere, next due to homogeneity the system model can be expanded to the whole space. This paper is organized as follows: In Section 2, the class of uncertain nonlinear systems considered in the identification scheme are described. In Section 3, the approximation property of NNs for homogeneous systems is studied. Identification algorithms as well as their convergence are discussed in Section 4. Section 5 presents the numerical results in order to demonstrate the performance of the identifiers. Section 6 concludes this manuscript with some remarks.

Notation: R + = {x ∈ R : x ≥ 0}, where R denotes the set of real numbers. For any ϑ ∈ R + and ∀x ∈ R we set ⌈x⌋ ϑ = sign(x)|x| ϑ . The Euclidean norm is denoted by .

. For M ∈ R m×n , M = [m 1 , m 2 , . . . , m n ], we set vec(M) = [m ⊤ 1 , m ⊤ 2 , . . . , m ⊤ n ] ⊤ .
The Kronecker product is denoted by ⊗.

PROBLEM STATEMENT

Let us consider a nonlinear affine control system:

ẋ = f (x, u) := f 0 (x) + m ∑ i=1 f i (x)u i , (1)
• x ∈ R n is the state vector of the system;

• u = (u 1 , .., u m ) ⊤ ∈ R m is the control input, m ≤ n; • f i : R n → R n , i = 0, 1, .
.., m are unknown nonlinear vector fields;

The system is studied under the following basic assumptions: Assumption 1. The vector fields f i , i = 0, 1, ..., m are continuous on the unit sphere:

S = {x ∈ R n : x = 1} . (2)
Assumption 2. The vector fields f i are homogeneous in the standard sense with known degrees

ν i ∈ R + , i.e. f i (λx) = λ ν i f i (x), ∀x ∈ R n , ∀λ > 0, ( 3
) where i = 0, 1, ..., m.
The standard homogeneity means that the function f is symmetric with respect to dilation x → λx of its first argument. The generalized concepts of homogeneity have been developed for other types of dilation D : R → R n×n , x → D(λ)x, and for both finite and infinite dimensional systems, (see e.g. [START_REF] Zubov | On systems of ordinary differential equations with generalized homogenous right-hand sides[END_REF], [START_REF] Khomenuk | [END_REF], [START_REF] Kawski | Geometric homogeneity and stabilization[END_REF], [START_REF] Polyakov | On homogeneous distributed parameter systems[END_REF]). The matrix D(λ) denotes a dilation in R n . In this manuscript we consider the standard homogeneity property since, in the finite dimensional case, any generalized homogeneous system is topologically equivalent to a standard homogeneous one [START_REF] Grüne | Homogeneous state feedback stabilization of homogenous systems[END_REF], [START_REF] Polyakov | Sliding mode control design using canonical homogeneous norm[END_REF]. Assumption 3. It is assumed that the whole state vector of (1) is bounded and sufficiently excited, i.e. 0 < x(t) < +∞, ∀t ≥ 0 and it can be measured, i.e. h(x) = x.

The main goal of the identification problem for the control system (1) under assumptions 2 and 3 is to find the model representation for the vector field f and the parameters for that model, such that, the error between the states of the real system and the proposed structure is bounded and small enough. In other words, the first step is to represent (1) with a valid model for its approximation. In this work we use a DNN structure for this purpose. Then, the second step is to use a DNN identifier structure and design the adaptive laws for the adjustment weights, such that, the error e := xx (4) between the system states and the identifier state x tends to zero or, at least, it is bounded, i.e. lim sup t→∞ e ≤ β < +∞.

Homogeneity allows local properties of vector fields to be extended globally. For example (see e.g. [START_REF] Bhat | Geometric homogeneity with applications to finite-time stability[END_REF] and Polyakov (2018)), a vector field f i : R n → R n satisfying Assumption 2, is Lipschitz continuous on R n \ {0} if and only if it satisfies Lipschitz condition on the unit sphere (2). Similarly, the system (1) satisfying the Assumptions 1 and 2 has the right-hand side continuous (on the first argument) in R n \{0}. For ν i = 0 the function f i may have discontinuity at the origin. Therefore, we need to develop an identifier which can deal with discontinuous models under considerations. Classical approximation theorems does not work in this case, since they usually assume (at least) continuity of the model on a compact.

NEURAL NETWORKS APPROXIMATION PROPERTY FOR HOMOGENEOUS SYSTEMS

The universal approximation property of NN has been used in many works, which says that any continuous function can be approximated arbitrarily closely on a compact set using superposition of nonlinear functions such as polynomials, radial basis functions and sigmoidal functions [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF], [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF]. Theorem 4. [START_REF] Haykin | Neural networks: a comprehensive foundation[END_REF]). Let σ(•) be a bounded monotoneincreasing continuous function and let C(I c p ) be the space of continuous functions defined in the unit hypercube of p dimensions I c p := [0, 1] p . Then, given any function f ∈ C(I c p ) and ε > 0, there exists an integer N and sets of real constants α i , b i , ω i j , where i = 1, 2, ..., N and j = 1, 2, ..., p, such that the following representation:

G(x) := N ∑ i=1 α i σ p ∑ j=1 ω i j x j + b i , (5)
is an approximate realization of the function f (•), i.e.

|G(x) -f (x)| < ε, ∀x ∈ I c p . (6)
Notice that Theorem 4 is formulated for a two layers static NN structure, the proof of this theorem is based on the Stone-Weisstrass and the Kolmogorov approximation theorems, this property is used in DNN structures for the representation of dynamic systems. Remark 5. In [START_REF] Lewis | Neural network control of robot manipulators and non-linear systems[END_REF] and [START_REF] Igelnik | Stochastic choice of basis functions in adaptive function approximation and the functional-link net[END_REF], it is stated that the parameters ω i j , and b i in the structure presented in (5) can be selected randomly with a uniform distribution. Then, in the aforementioned works, it is shown that the universal approximation property holds by finding only the parameters α i .

The next result uses the universal approximation property of NN structures for the case of the dynamic nonlinear system (1) satisfying assumptions 2 and 3. Corollary 6. Let the system (1) satisfy the assumptions 1, 2 and 3. Then, for any ε i ∈ R + and for any matrix A ∈ R n×n there exist W i ∈ R n×N i , i = 0, 1, .., m such that:

f (x, u) -F(x, u) ≤ ε 0 x ν 0 + m ∑ i=1 ε i x ν i |u i | ∀x ∈ R n , ∀u ∈ R m , where F(x, u):= x ν 0 Ax x + W 0 σ 0 x x + m ∑ i=1 x ν i W i σ i x x u i (7)
and the elements of the vector fields σ i (•) are sigmoidal activation functions:

(σ i (x)) j = α i j 1 + c i j e -b ⊤ i j x -1 , ( 8
)
where α i j ∈ R + , c i j ∈ R + and b i j ∈ R n are some properly selected parameters with i = 0, 1, ..., m and j = 1, ..., N i .

Proof. Obviously that under Assumption 2 the system (1) can be rewritten as:

ẋ = x ν 0 f 0 x x + m ∑ i=1 x ν i f i x x u i , (9)
i.e. the right-hand side of the system (1) is uniquely identified by its values on the unit sphere (2).

Let us consider the functions fi :

R n → R n , i = 0, 1, ..., m defined as follows fi (x) = f i x x
.

Obviously, fi is continuous on R n due to continuity of f i on the unit sphere.

Applying Theorem 4 to each component of the vector fi , i = 1, 2, ...m and to f0

x x -A x x on the unit sphere S, we obtain that

W i σ i x x approximates fi (x) = f i x x
with an error ε i , and this error can be made arbitrary small by means of a proper selection of parameters W i , α i j , b i j and c i j . Finally, taking into account homogeneity of nonlinear functions f i , i = 0, ..., m we complete the proof.

In the view of Remark 5 below we assume that appropriate parameters α i j , b i j and c i j , are somehow selected, and we just need to find matrices W i in order to identify the model.

IDENTIFICATION OF AFFINE HOMOGENEOUS CONTROL SYSTEMS

The case of known control gains

Let us consider initially the case when the nonlinear maps f i : R n → R n , i = 1, ..., m associated with the inputs are known and we need to identify only the vector field f 0 .

In the most straightforward case, system (1) is represented with an exact NN matching (ε 0 = 0). The following modeling considers that for a selected Hurwitz matrix A ∈ R n×n , there exists a weight matrix W * 0 ∈ R n×N 0 such that f 0 can be exactly presented by the NN structure described in Corollary 6. Notice that

W * 0 σ 0 x x := Σ 0 x x w * 0 , (10)
where

Σ 0 (z) = I n ⊗ σ ⊤ 0 (z) ∈ R n×nN 0 , z ∈ R n , I n ∈ R n×n is the identity matrix and w * 0 = vec (W * 0 ) ⊤ ∈ R nN 0 .
The identification problem can be seen as finding w * 0 by using some learning law such that x can be reproduced by x, where x ∈ R n represents the state vector of the DNN identifier (whose equations are given below in ( 11)).

The following theorem introduces the first result for the convergence of the identification error. Theorem 7. Let assumptions 1, 2 and 3 be satisfied and the homogeneous vector field f 0 admits the exact representation (7) and the DNN identifier is defined as follows:

d dt x = x ν 0 A x x + Σ 0 x x w 0 + ΩKΩ ⊤ e + m ∑ i=1 f i (x) u i , (11
) where A ∈ R n×n is a Hurwitz matrix, e := xx is the identification error, w 0 ∈ R nN 0 is the vector of the weights to be adjusted as follows

d dt w 0 = -x ν 0 KΩ ⊤ e (12)
and Ω ∈ R n×nN 0 is an auxiliary variable satisfying

d dt Ω = x ν 0 A Ω x -Σ 0 x x . ( 13
)
If K ∈ R nN 0 ×nN 0 is a positive definite matrix and t+ℓ t Ω ⊤ (τ)Ω(τ)dτ ≥ υI nN 0 , ∀t ∈ R + , (14)
for some ℓ > 0 and υ > 0. then, the vector of weights w 0 converge to w * 0 asymptotically and the identification process is asymptotically consistent, i.e. lim t→∞ e(t) = 0.

Proof. To prove the stability of the DNN as well as the stability of the identification error e, let us define the following auxiliary equation: δ = e + Ω ω, (15) where ω := w * 0w 0 . Hence, the dynamics of the auxiliary equation is:

δ = ė + dΩ dt ω + Ω d ω dt , ( 16
)
the dynamics of the identification error is:

ė = x ν 0 Ae x + Σ 0 x x ω -ΩKΩ ⊤ e , (17)
the time derivative d ω dt = -dw 0 dt , by substituting ( 12), ( 13) and ( 17) in ( 16), we have:

δ = x ν 0 A e x + Σ 0 x x ω -ΩKΩ ⊤ e + x ν 0 A Ω x -Σ 0 x x ω + x ν 0 ΩKΩ ⊤ e = x ν 0 A e x + A Ω x ω = x ν 0 A x δ. ( 18
)
Since the matrix A is Hurwitz and x = 0, then the condition (14) implies δ → 0 and e → -Ω ω as t → +∞. Since the matrix A is Hurwitz and the state x(t) is bounded due to Assumption 3, it is straightforward to observe that the variable Ω(t) is also bounded.

Using the selected learning law ( 12) and the fact on the stability of the dynamic auxiliary equation ( 16), we obtain:

d ω(t) dt = x ν 0 KΩ ⊤ (Ω ω -δ) , (19)
hence, since the condition of persistence excitation ( 14) is satisfied, the latter system is input-to-state stable with respect to the input δ (see e.g. [START_REF] Efimov | Design of impulsive adaptive observers for improvement of persistency of excitation[END_REF]). Therefore, for asymptotically converging signal δ, the asymptotic convergence of ω (and e, respectively since both δ and ω are converging) to zero can be established.

Notice that the DNN identifier ( 11) is not regular [START_REF] Poznyak | Differential neural networks for robust nonlinear control: identification, state estimation and trajectory tracking[END_REF], [START_REF] Chairez | Adaptive neural network nonparametric identifier with normalized learning laws[END_REF] because it has a direct injection of the identification error e.

The case of unknown control gains

The previous design considers the nonlinear part associated to the input as known. In this section we design a DNN identifier for the considered system assuming that functions f i are unknown but representation (7) admits the exact match (ε 0 = 0 and ε i = 0).

Similarly to (10) we introduce the vectors w * i and the matrixvalued function Σ i such that

W * i σ i x x := Σ i x x
w * i . Theorem 8. Let the system (1) admits an exact representation of the form (7), assumptions 1-3 be satisfied and the control input u be selected as follows:

u i (t) = ũi (t) x(t) ν i -ν 0 , i = 1, . . . , m, ( 20
)
where ũi are continuous uniformly bounded functions and the DNN identifier is defined as follows

d x dt = x ν 0 -1 A x + m ∑ i=0 x ν 0 Σ i x x w i ũi + Ω i K i Ω ⊤ i e , ( 21
)
where u 0 ≡ 1, A ∈ R n×n is a Hurwitz matrix and e := xx is the identification error, w i ∈ R nN i are the vectors of the weights to be adjusted as follows:

d dt w i = -x ν 0 K i Ω ⊤ i e (22)
and Ω i ∈ R n×nN i are auxiliary variables satisfying:

d dt Ω i = x ν 0 -1 AΩ i -ũi x ν 0 Σ i x x , (23)
If K i ∈ R nN i ×nN i , i = 0, 1, .
.., m are positive definite matrices and:

t+ℓ W t H ⊤ (τ)H(τ)dτ ≥ υ W I ∑ m i=0 N i ( 24
)
for all t ∈ R + , some ℓ W > 0, υ W > 0, where H ∈ R n×∑ m i=0 N i is a block matrix H = {Ω i }.
Then, the vectors of weights w i converge to w * i asymptotically and the identification process is asymptotically consistent, i.e. lim t→∞ e(t) = 0.

Proof. To prove the stability of the DNN as well as the stability of the identification error e, let us define the following auxiliary equation:

δ = e + m ∑ i=0 Ω i ωi , (25)
where ωi := w * iw i . Hence, the dynamics of the auxiliary variable is:

δ = ė + m ∑ i=0 dΩ i dt ωi + Ω i d ωi dt = x ν 0 -1 Ae + m ∑ i=0 x ν 0 Σ i x x ω ũi -Ω i K i Ω ⊤ i e + m ∑ i=0 x ν 0 -1 AΩ i -ũi x ν 0 Σ i x x ωi + x ν i Ω i KΩ ⊤ i e = x ν 0 -1 (Ae + A ∑ m i=0 Ω i ωi ) = x ν 0 -1 Aδ,
Since the matrix A is Hurwitz and x = 0, then the condition (24) implies δ → 0 and e → -∑ m i=0 Ω i ωi as t → +∞. The variables Ω i (t) are bounded due to the structure of ( 23), Hurwitz property of the matrix A, boundedness and separation from zero of u i (t) and x(t) .

Using the selected learning law ( 22) and the fact on the stability of the dynamic auxiliary equation ( 25), we obtain:

d ωi (t) dt = x ν i K i Ω ⊤ i m ∑ j=0 Ω j ω j -δ , (26)
hence, since the condition of persistence excitation ( 24) is satisfied, the latter system is input-to-state stable with respect to the input δ (see e.g. [START_REF] Efimov | Design of impulsive adaptive observers for improvement of persistency of excitation[END_REF]). Therefore, for asymptotically converging signal δ, the asymptotic convergence of ωi (and e, respectively) to zero can be established.

NUMERICAL RESULTS

Let us consider the three tank system depicted in figure 1. This nonlinear system can be modeled as in [START_REF] Join | Control of an uncertain three-tank system via on-line parameter indentification and fault detection[END_REF][START_REF] Seydou | Actuator fault diagnosis for flat systems: a constraint satisfaction approach[END_REF][START_REF] Zimenko | A note on delay robustness for homogeneous systems with negative degree[END_REF] by:

ẋ1 = 1 S tank -a 13 ⌈x 1 -x 3 ⌋ 0.5 + u 1 , ẋ2 = 1 S tank a 32 ⌈x 3 -x 2 ⌋ 0.5 -a 20 ⌈x 2 ⌋ 0.5 + u 2 , ẋ3 = 1 S tank a 13 ⌈x 1 -x 3 ⌋ 0.5 -a 32 ⌈x 3 -x 2 ⌋ 0.5 , (27)
where x 1 , x 2 and x 3 represent the liquid level of each tank respectively, S tank is the diameter of the three tanks, the input flows u 1 and u 2 are the control signals and the constant parameters a 13 , a 32 and a 20 are coefficients related with the outflow rate according to Torricelli's rule.

Checking assumptions 2 and 3 for the system (27), it is straightforward to check that the system is of homogeneity degree ν = -0.5. In table 1, the values used to simulate the three thank system are presented, this values are not used on the identifier.

Parameter Description Value

S tank

Tank diameter 1m

a 13 Outflow rate coefficient 3 m 3 2 /s 2 a 32 Outflow rate coefficient 2 m 3 2 /s 2 a 20 Outflow rate coefficient 1 m 3 2 /s 2 x 0 Initial conditions [3, 1, 2] ⊤ λ Real positive value 0.4
Table 1. Parameters for the simulation of the three tank system

The numerical simulations for the identifiers were made in Simulink Matlab ® by using the Runge Kutta integration method with a step of 1 ms.

For the first identifier, let us consider the vector of activation functions σ 0 (•) in ( 8) with N 0 = 3 and the constant parameters α 0,1 = 2.7, α 0,2 = 1.8, α 0,3 = 3.6, c 0,1 = 8, c 0,2 = 16, c 0,3 = 8, the constant vectors:

the matrix A and the functions vectors f i are :

A := 0 1 0 0 0 1 -2 -3 -1 ; f 1 (x) = 1 0 0 ; f 2 (x) = 0 1 0 , f 0 (x) = 1 S tank -a 13 ⌈x 1 -x 3 ⌋ 0.5 a 32 ⌈x 3 -x 2 ⌋ 0.5 -a 20 ⌈x 2 ⌋ 0.5 a 13 ⌈x 1 -x 3 ⌋ 0.5 -a 32 ⌈x 3 -x 2 ⌋ 0.5 , x = (x 1 x 2 , x 3 ) ⊤ .
The initial conditions for the adjustment law are w 0 (0) = [1 2 1 1 3 2 1 2 3] ⊤ and the gain matrix is selected as K := I 3 ⊗ K, where the matrix K is defined as: K := 100 -20 30 -20 30 10 30 10 100

.

For the second identifier with an unknown nonlinear term associated to the input was used the same integration method and matrix A. The activation functions were selected as in ( 8) with all N i = 3. The constant parameters α 0, j , c 0, j and the vectors b 0, j are the same as the used in the first identifier. The constant parameters for σ i (•) were selected as α i, j = α 0, j and c i, j = c 0, j and the vectors b i, j = b 0, j .

The initial conditions for the adjustment law were selected as w 0 (0) is the same as the used for the first identifier and w 1 (0) = w 0 (0) and w 2 (0) = w 0 (0) and the gain matrices K i := I 3 ⊗ K, Figures 2, 3, 4 and 5 show the obtained results. The first identifier (-blue) uses a known gain of the input. Its state estimate converges quickly to the state of the uncertain system (-black). The convergence rate of the second identifier (• • • red) is much slower but it does not assume that functions f 1 and f 2 are known. In this paper, DNN identifiers for affine standard homogeneous (possibly discontinuous) uncertain control system are developed. Their convergence is proven theoretically and their precision confirmed by numerical simulation for homogeneous model of three tank system. The main assumption required for the design of a homogeneous DNN identifier is knowledge of the homogeneity degree. The problem of identification of homogeneity degree is planned to be studied in a recent future.

Fig. 1 .

 1 Fig. 1. Three tank system

Fig. 5 .

 5 Fig. 2. Identification result for the first state (x 1 )