Rafael Ferro
email: rafael.ferro@ircam.fr

Nicolas Obin
email: nicolas.obin@ircam.fr

Axel Roebel
email: axel.roebel@ircam.fr

CycleGAN Voice Conversion of Spectral Envelopes using Adversarial Weights

Keywords: voice conversion, cycleGAN, GAN stability, adversarial weights

HAL is

I. INTRODUCTION

A. Related works

Voice identity conversion (VC) consists in modifying the voice of a source speaker so as to be perceived as the one of a target speaker. Over the past few years, VC has largely gained in popularity and in quality [START_REF] Toda | The voice conversion challenge 2016[END_REF], [START_REF] Lorenzo-Trueba | The voice conversion challenge 2018 : Promoting development of parallel and nonparallel methods[END_REF], in particular with the development of neural voice conversion algorithms [START_REF] Tanaka | AttS2s-VC : Sequence-to-Sequence Voice Conversion with Attention and Context Preservation Mechanisms[END_REF], [START_REF] Kameoka | ConvS2s-VC : Fully convolutional sequence-to-sequence voice conversion[END_REF]. VC consists in learning a conversion function between the acoustic space of a source and a target speaker. This conversion function is generally learned from a pre-aligned database (parallel VC) in which the source and the target speakers pronounce the same set of sentences, so that a direct correspondence between the frames of the source and target speakers can be established. Unfortunately, this constraint reduces the amount of available recordings of source and target speakers.

Modern VC is mainly based on neural architectures with the particular objective to extend the VC from parallel to nonparallel speech databases. The main advantage of non-parallel VC is that it provides the flexibility to learn the conversion from "on-the-fly" speech databases, which can more easily handle large amount of data and accommodate multiple speakers. These architectures includes Variational Autoencoders (VAEs) [START_REF] Kameoka | ACVAE-VC : Non-parallel many-to-many voice conversion with auxiliary classifier variational autoencoder[END_REF]- [START_REF] Huang | Refined WaveNet Vocoder for Variational Autoencoder Based Voice Conversion[END_REF], Generative Adversarial Networks (GANs) [START_REF] Kameoka | StarGAN-VC : Non-parallel many-to-many voice conversion with star generative adversarial networks[END_REF]- [START_REF] Fang | High-quality nonparallel voice conversion based on cycleconsistent adversarial network[END_REF], Phonetic PosteriorGrams (PPGs) [START_REF] Sun | Phonetic posteriorgrams for many-to-one voice conversion without parallel data training[END_REF] and sampleRNNs

The research in this paper has been funded by the ANR project TheVoice : ANR-17-CE23-0025 [START_REF] Zhou | Voice Conversion with Conditional SampleRNN[END_REF] among others. The use of GAN architectures [START_REF] Goodfellow | Generative Adversarial Networks[END_REF] for VC is inspired by advances conducted in the fields of image generation and manipulation. The cycleGAN is a particular configuration of GANs which has been specifically formulated to learn transformations between two different domains or between unaligned or unpaired datasets with application to image-to-image translation [START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF].

The cycleGAN-VC is the extension of the cycleGAN to the VC task, which has become a standard in non-parallel VC [START_REF] Kameoka | StarGAN-VC : Non-parallel many-to-many voice conversion with star generative adversarial networks[END_REF], [START_REF] Kaneko | Parallel-Data-Free Voice Conversion Using Cycle-Consistent Adversarial Networks[END_REF]. The main idea behind the application of the cycleGAN to VC is that the cycle-consistency encourages the preservation of the phonetic content through the cycle while learning to modify the speaker identity. Despite the advances accomplished in non-parallel VC, cycleGAN-VC still suffers from important limitations which conduct to conversion of mitigated quality. One main limitation is due to the wellknown stability issues of the GAN [START_REF] Salimans | Improved Techniques for Training GANs[END_REF], [START_REF] Radford | Unsupervised Representation Learning with Deep Convolutional Generative Adversarial Networks[END_REF]. This issue, combined to the limited amount of data available in the VC task, can lead to severe degradation of the voice conversion quality and the naturalness of the converted speech.

Tackling GAN stability issues has motivated the development of multiple ideas and heuristics [START_REF] Salimans | Improved Techniques for Training GANs[END_REF], such as the popular Wasserstein GAN [START_REF] Arjovsky | Wasserstein Generative Adversarial Networks[END_REF] or the LSGAN [START_REF] Mao | Least Squares Generative Adversarial Networks[END_REF]. More recently, so as to tackle stability issues, the weighted GAN has been introduced [START_REF] Pantazis | Training Generative Adversarial Networks with Weights[END_REF]. This novel approach, based on game theory, instead of equally weighted "fake" samples, more attention is given to samples that fool the discriminator. A particular form of weighted GAN has been recently applied to VC. [START_REF] Paul | Non-Parallel Voice Conversion Using Weighted Generative Adversarial Networks[END_REF].

B. Contribution of the paper

This paper proposes two contributions. First, we propose to use spectral envelopes as inputs instead of using cepstral coefficients as in [START_REF] Kaneko | Parallel-Data-Free Voice Conversion Using Cycle-Consistent Adversarial Networks[END_REF]. The use of the spectral envelope is assumed to simplify the task for the convolutional networks that will be required to only slightly move the spectral formants. Thus we assume to not need an extremely deep network and we will show that we can achieve similar performance with a significantly smaller network. Second, we propose to introduce a constraint as an additional training loss that enforces to preserve energy contour of the converted speech signal. This allows us to add a generator loss not depending on adversarial training, which is expected to contribute to the stability of the training. Third, we present a novel method to tackle the stability issue of GAN training, exploring a novel weighted GAN approach. We achieve this by adding a weight to the loss of the discriminator, giving more weight to "true" samples rather than to "fake" ones.

II. PRELIMINARY WORKS ON CYCLEGAN VC

A. Generative Adversarial Networks

A Generative Adversarial Network (GAN) [START_REF] Goodfellow | Generative Adversarial Networks[END_REF] is a neural network system composed by a generator G and a discriminator D, in which the discriminator is trained to discriminate real samples from generated samples, while the generator is trained to generate real-like samples, using the discriminator as the decision rule. The objective can be written as :

min G max D L(D, G) =E x⇠p(X) [log D(x)] +E z⇠p(Z) [log(1 D(G(z)))] (1)
where : x is a sample from a distribution p(X) to be modeled, z is a sample generated from a random distribution p(Z), and E x⇠p(X) represents the expected value of x given the distribution p(X).

B. Cycle Generative Adversarial Networks

In the cycleGAN architecture [START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF], [START_REF] Yi | DualGAN : Unsupervised Dual Learning for Image-to-Image Translation[END_REF], a generator G X!Y reads data from a dataset X and learns to map it into its respective position in a dataset Y, and vice versa for a generator G Y !X . If X and Y represent languages, this system should be analogous to two translators. To train these generators, the cycleGAN framework uses two adversarially trained discriminators to discriminate respectively any x 2 X in relation to G Y !X (y) for any y 2 Y and any y 2 Y in relation to G X!Y (x) for any x 2 X. Since G Y !X (G X!Y (x)) should be equal to x, and G X!Y (G Y !X (y)) should be equal to y, a loss named cycle-consistent loss is added to enforce this constraint. In the following equations we state the total objective of the cycleGAN, where E y⇠P Data (y) represents the expected value for the distribution Y and E Y ⇠P Data (X) represents the expected value for the distribution X.

The following equation describes the adversarial loss for the discriminator D Y (the equation for D X is analogous) :

L adv (GX!Y ,DY ) = E y⇠P Data (y) [log(DY (y))] +E x⇠P Data (x) [log(1 DY (GX!Y (x)))]. (2)
The following equation describes the cycle-consistency loss, using L1 norm :

Lcyc(GX!Y ,GY !X ) = E x⇠P Data (x) [||GY !X (GX!Y (x)) x||1] +E y⇠P Data (y) [||GX!Y (GY !X (y)) y||1]. (3)
The following equation describes the total objective of the cycleGAN, where c represents the weight for the cycleconsistency loss :

L full =L adv (GX!Y , DY ) + L adv (GY !X , DX ) + cLcyc(GX!Y , GY !X ) (4)
The cycleGAN-VC, introduced by [START_REF] Kaneko | Parallel-Data-Free Voice Conversion Using Cycle-Consistent Adversarial Networks[END_REF], is trained to convert a source speaker Mel Frequency Cepstral Coefficients (MFCCs) into a target speaker MFCCs, so as to perform VC. Their discriminators task is therefore to discriminate whether the conversions belong to their respective target speaker identity or not. In particular, so as to adapt the original cycleGAN framework, they used Gated CNNs as well as an identity-mapping loss, which is reported to encourage phonetic invariance. However, we did not find these two additional ideas beneficial.

III. GAN WITH ADVERSARIAL WEIGHTS

Generative Adversarial Networks are reported to be difficult to train. One problem are vanishing gradients when the discriminator achieves perfect discrimination, or when the generator is able to perfectly fool the discriminator, though it is producing nonsense. Another problem is the instability that is due to the fact that the discriminator is trained to systematically reject generated examples independent of their quality. In the case when the generator generates target samples covering only a small part of the target space the discriminator will improve its objective by means of pushing the generator out of the target space even if it has to wrongly classify some of the real samples as well. As a result the discriminator will push the generator away from the target space hindering the generator to converge.

To solve these issues, many ideas have been proposed, such as the DCGAN architecture [START_REF] Radford | Unsupervised Representation Learning with Deep Convolutional Generative Adversarial Networks[END_REF]. [START_REF] Gulrajani | Improved Training of Wasserstein GANs[END_REF] discussed mini-batch discrimination, historical averaging, one-sided label smoothing and virtual batch normalization. Also, new losses have been proposed, such as the Wasserstein GAN [START_REF] Arjovsky | Wasserstein Generative Adversarial Networks[END_REF] and the LSGAN [START_REF] Mao | Least Squares Generative Adversarial Networks[END_REF].

A. Weighted GAN System

Recently, so as to tackle CycleGAN optimization, Paul et al. developed the weStarGAN [START_REF] Paul | Non-Parallel Voice Conversion Using Weighted Generative Adversarial Networks[END_REF], implementing the weighted GAN idea [START_REF] Pantazis | Training Generative Adversarial Networks with Weights[END_REF] to the starGAN-VC architecture [START_REF] Kameoka | StarGAN-VC : Non-parallel many-to-many voice conversion with star generative adversarial networks[END_REF]. They do so by multiplying sample-wise a coefficient to the generator loss. Their idea is to give less weight to samples poorly produced by the generator, so that the generator has a stronger motivation to produce samples similar to the real data. For each sample j they compute its respective normalised weight w j,g , while introducing a hyper-parameter ⌘ gen : wj,g = e ⌘min(0,D j ) wj,g = wj,g P m j=1 wj,g

, j = 1, ..., m. (5)
This coefficient is then multiplied sample-wise by the generator loss. Note that samples that the discriminator does not correctly discriminate have a higher weight in the generator loss. A hypothesis assumed in [START_REF] Paul | Non-Parallel Voice Conversion Using Weighted Generative Adversarial Networks[END_REF] is that the discriminator is "faithful" : that quantitatively it returns on average values above 0.5 when the samples come from the real distribution and below 0.5 when fake samples are fed to the Discriminator. This structure is our starting point in our research on the stability and optimization of the GAN.

C. Generalized Adversarial Weights

Finally, we introduce a novel weighted GAN technique, also encouraging the discriminator "faithfulness". In the proposed approach, the idea is to encourage the training of the discriminator by weighting samples produced by the generator by a constant generator impact hyperparameter ⇢ gen , so that the discriminator is rather encouraged by the "real" distribution than by the generated one. We name this system the generator impact GAN (gimGAN), since we can control the impact of the generator on the discriminator loss. The idea is equivalent to replacing the "hard" label 0 scored by the discriminator on the samples produced by the generator D(G(x)), by a "soft" label calculated as:

D(G(x)) + ⇢ gen (0 D(G(x))) 2 [0, 1] (6)
where: ⇢ gen (0 D(G(x))) represents the current distance of the discriminator relatively to the target "hard" label 0, weighted by generator impact hyperparameter ⇢ gen . Once again, the motivation behind is to tackle the fact that in the vanilla GAN, the discriminator is trained to consider samples produced by the generator as "fake" ones, whatever its quality is, which might have the effect to discourage the generator to converge to the correct distribution.

was superior to 2s.

Our architecture was inspir cycleGAN-VC [START_REF] Sun | Phonetic posteriorgrams for many-to-one voice conversion without parallel data training[END_REF], [START_REF] Zhou | Voice Conversion with Conditional SampleRNN[END_REF]. Fo two-layer encoder with convo bottleneck with convolutions with transposed convolutions. decoder, we applied a kerne so as to avoid the checkerb consistently better results were effect was avoided. For the bo 3 and and stride size 1. The o 256, 512, 512 , 512, 256 and normalization, followed by a Since spectral envelopes expl tion, the task becomes much a rather small generator, com implementations. [START_REF] Sun | Phonetic posteriorgrams for many-to-one voice conversion without parallel data training[END_REF], [START_REF] Zhou | Voice Conversion with Conditional SampleRNN[END_REF] F 4 convolutional layers with a with respectfully 64, 128, 2 layers were followed by two f and 1 neurons respectively. We followed by a LeakyReLU a for the last one. Our inputs

B. Generalized Weighted GAN

Since in practice the discriminator might not be "faithful", in order to encourage its "faithfulness", we propose to train the discriminator rather by fairly produced samples than by poorly produced ones. This helps to reduces the above mentioned instability caused by training the discriminator to judge generated data, no matter how good it is, as equally wrong. To do so, we further develop Paul et al. idea, by also applying weights to the discriminator loss. We followed the exact same procedure : we introduce a hyper-parameter ⌘ dis and for each sample we compute its respective normalised weight w j,d , following the same procedure as in Eq. ( 5). This coefficient is then multiplied sample-wise by the discriminator loss. Thus samples perceived as poor by the discriminator have less weighting in the discriminator loss. We name this system as generalized weighted GAN (geweGAN).

C. Generator Impact GAN

Finally, we introduce a novel weighted GAN technique, also encouraging the discriminator "faithfulness". In the proposed approach, the idea is to encourage the training of the discriminator by weighting samples produced by the generator by a constant generator impact hyperparameter ⇢ gen , so that the discriminator is rather encouraged by the "real" distribution than by the generated one. We name this system the generator impact GAN (gimGAN), since we can control the impact of the generator on the discriminator loss. The idea is equivalent to replacing the "hard" label 0 scored by the discriminator on the samples produced by the generator D(G(x)), by a "soft" label calculated as :

D(G(x)) + ⇢gen(0 D(G(x))) 2 [0, 1] (6)
where : ⇢ gen (0 D(G(x))) represents the current distance of the discriminator relatively to the target "hard" label 0, weighted by generator impact hyperparameter ⇢ gen . Once again, the motivation behind is to tackle the fact that in the vanilla GAN, the discriminator is trained to consider samples produced by the generator as "fake" ones, whatever its quality is, which might have the effect to discourage the generator to converge to the correct distribution. The principle is illustrated in Figure 1.

D. Energy constraint

Additionally, we further apply an energy constraint to the cycleGAN architecture. This constraint enforces preservation of the energy contour of the original source speech signal during conversion, and avoids incoherence between source and converted envelope in the converted speech signal. Furthermore, by means of providing stable feedback to the generator, this constraint is expected to reduce the instability of the GAN training. To achieve this, we impose a reconstruction loss on the amplitude mean for each frame, on both generators. Since we work with spectral envelopes as inputs, we just add the following term to the total loss, where c represents the weight for the cycle-consistency loss :

Le = eEx⇠P Data (x) [|| T X t=0 GX!Y (x) x||1] + eEx⇠P Data (x) [|| T X t=0 GY !X (y) y||1] (7)
This constraint enforces preservation of the energy contour of the original source speech signal during conversion, and avoids incoherence between source and converted envelope in the converted speech signal. Further, by means of providing stable feedback to the generator, this constraint is expected to reduce the instability of the GAN training.

IV. EXPERIMENT

The proposed CycleGAN architectures have been trained and evaluated, using the VCC2018 database [START_REF] Lorenzo-Trueba | The voice conversion challenge 2018 : Promoting development of parallel and nonparallel methods[END_REF]. The VCC2018 training corpus contains 80 short sentences per speaker, sampled at 16 kHz and quantified on 16 bits. For the evaluation set, we used the first 5 sentences, whose length was superior to 2s.

Our architecture was inspired by the DCGAN and by the cycleGAN-VC [START_REF] Kaneko | Parallel-Data-Free Voice Conversion Using Cycle-Consistent Adversarial Networks[END_REF], [START_REF] Kaneko | CycleGAN-VC2 : Improved CycleGAN-based Non-parallel Voice Conversion[END_REF]. For both generators, we used a two-layer encoder with convolutions, followed by a two-layer bottleneck with convolutions and then a two-layer decoder with transposed convolutions. For both the encoder and the decoder, we applied a kernel size 2 and a stride size 2, so as to avoid the checkerboard effect [START_REF] Odena | Deconvolution and checkerboard artifacts[END_REF], noticing that consistently better results were obtained when this undesirable effect was avoided. For the bottleneck, we applied kernel size 3 and and stride size 1. The overall generator has respectfully 256, 512, 512 , 512, 256 and 1 filters. We applied instance normalization, followed by a ReLU at the end of each layer. Since spectral envelopes explicitly contain formant information, the task becomes much easier allowing us to implant a rather small generator, compared to traditional cycleGAN implementations. [START_REF] Kaneko | Parallel-Data-Free Voice Conversion Using Cycle-Consistent Adversarial Networks[END_REF], [START_REF] Kaneko | CycleGAN-VC2 : Improved CycleGAN-based Non-parallel Voice Conversion[END_REF] For both discriminators, we used 4 convolutional layers with a filter size 2, a kernel size 2, with respectfully 64, 128, 256 and 512 filters. These four layers were followed by two fully connected layers, with 512 and 1 neurons respectively. We applied instance normalization, followed by a LeakyReLU at the end of each layer, except for the last one. Our inputs were of size 32 and 128, for frequency bins and time frames respectively. We chose c and e to respectfully be 0.3 and 1. By lowering the traditional values for the cycle-consistency loss, we rather force GAN learning than cycle-consistency learning. This means that more weight is given to the identity learning, the task of the GAN, than to the phoneme reconstruction learning, the cycle-consistency task. In fact, otherwise results were found to be closer to a reconstruction, rather than a conversion. We applied a batch size 1 and we used least squares error for the discriminator loss, introduced in the LSGAN [START_REF] Mao | Least Squares Generative Adversarial Networks[END_REF] and optimized it with the Adam algorithm, as in [START_REF] Kaneko | Parallel-Data-Free Voice Conversion Using Cycle-Consistent Adversarial Networks[END_REF], [START_REF] Kaneko | CycleGAN-VC2 : Improved CycleGAN-based Non-parallel Voice Conversion[END_REF] for a total of 800k iterations, with a generator learning rate of 0.0002 and a discriminator learning rate of 0.0001. We implemented three GAN variations that all share the same network structure and input representation. The first variation, uses the Weighted GAN paradigm (weGAN) from [START_REF] Paul | Non-Parallel Voice Conversion Using Weighted Generative Adversarial Networks[END_REF] with ⌘ gen = 0.1. For the second variation, we implemented the generalized weighted GAN paradigm (geweGAN) for both the discriminator and the generator losses, as discussed above, with ⌘ gen = 0.9and ⌘ dis = 0.9. For the third variation, we implemented the generator impact weighted GAN paradigm (gimGAN), as discussed above, with ⇢ gen = 0.9. Finally, we also implemented cycleGAN-VC [START_REF] Kaneko | Parallel-Data-Free Voice Conversion Using Cycle-Consistent Adversarial Networks[END_REF] as a baseline.

Similarly to previous research on VC, the proposed VC is focused on spectral voice conversion only. The VC is based on a source/filter decomposition of the speech signal, in which the excitation of the source speaker is preserved during conversion and only the spectral envelope conversion is learned and modified. The analysis/synthesis engine relies on superVP, an extended phase vocoder developed by IRCAM 1 . The spectral envelope is estimated from the short-term Fourier transform (STFT) by using the True Envelope algorithm [START_REF] Röbel | On cepstral and all-pole based spectral envelope modeling with unknown model order[END_REF]. The Mel spectral envelope is then computed by integrating the 1. www.forumnet.ircam.fr/product/supervp-max-en/ estimated spectral envelope over 32 Mel filters in which the energy of each Mel filter is normalized to unity.

A. Experimental setups

The experiment consisted into the judgment by listeners of singing voice samples, based on the similarity to the target singer and the naturalness of the singer, as used for the voice conversion 2018 challenge [START_REF] Lorenzo-Trueba | The voice conversion challenge 2018 : Promoting development of parallel and nonparallel methods[END_REF]. Conversion were processed for all sentences contained in the test set. For the perceptual experiment, short excerpts were used and presented to the participants (around 5s.). We chose SF4 and TM4 as the source and TF3, TF4, TM3 and TM4 as the target. We evaluated the naturalness and speaker similarity of the converted samples, with a mean opinion score (MOS) test.

During the experiment, 15 short speech samples, original source and target speakers, and converted source-to-target speaker (each having duration of about 5s) were randomly selected from the test set, and presented to the participant in a random order. For each speech sample, the participant has the possibility to listen to an excerpt of the original target speaker. Then the participant is asked to rate the naturalness of the converted speech sample and its similarity to the target speaker. The experiment was conducted on-line, encouraging the use of headphones and quiet environment. 15 individuals participated in the experiment.

V. RESULTS AND DISCUSSION

The results of the perceptual evaluation are presented in table I. An overall result is that the original target speaker is consistently qualified to have high similarity and quality. In looking into the average results for all speaker pairs it is easy to see that the gimGAN and geweGAN variants perform significantly better than weGAN variant. This shows that the proposed GAN modifications used for training geweGAN and gimGAN have a positive impact on similarity and naturalness. Note that the energy constraint had an important contribution to these results. Without this constraint the similarity and naturalness rating are about 0.5 points lower (results not shown). Finally, the cycleGAN-VC baseline achieves approximately the same performance in similarity, while gimGAN and wegeGAN are perceived as more natural with gimGAN achieving overall the best results where differences in naturalness are more pronounced and differences in similarity remain marginal. Note however, that the geweGAN and gimGAN networks achieve this performance with about 3 times less parameters, which results in considerably shorter training and inference times, which is expected to be a result of directly working on the spectral envelope. Looking into the sub groups for MTM and FTF conversions, one can notice that the best network changes with gender. While cycleGAN-VC baseline is best for MTM conversion, wegeGAN has been evaluated best in the FTF case. These differences in the two gender groups indicate a relatively strong dependency of the conversion performance on the speaker.

VI. CONCLUSION

The present paper investigates Voice Conversion with Generative Adversarial Networks (GAN), particularly with the cycleGAN paradigm, addressing optimization and stability issues. First, we use spectral envelopes as inputs. Second, so as to optimize and to address the stability issues of the GAN training we propose a generalization of the weighted GAN, the geweGAN, and a similar approach, the gimGAN. Our conducted experiment shows first that both proposed methods are able to score a better performance than the previously proposed weighted GAN. Second, it shows that the proposed method performs similarly to the cycleGAN-VC baseline on similarity and considerably outperforms it on quality. Furthermore, these results were achieved using a significantly smaller network, which significantly reduces training time. Finally, an additional energy constraint to the loss was found to be essential for similarity and naturalness learning. For future work, we plan to combine the proposed generalized weighted GAN with the generator impact GAN, so as to further improve the stability of the GAN training procedure.

R ÉF ÉRENCES

FIGURE 1 .

 1 FIGURE 1. Comparison of GAN with and without generator impact. On left : vanilla GAN, on right : proposed GAN with weighted labels.

TABLE I MOS

 I AND 95% CONFIDENCE INTERVAL OBTAINED FOR THE DIFFERENT VC SYSTEMS.

		Male-to-Male (MTM)	Female-to-Female (FTF)	TOTAL
	Speech Signal Class	Similarity	Naturalness	Similarity	Naturalness Similarity	Naturalness
	orig : target conv : weGAN conv : geweGAN conv : gimGAN conv : cycleGAN-VC (baseline)	4.96 ± 0.09 2.52 ± 0.42 2.80 ± 0.76 3.21 ± 0.37 3.33 ± 0.46	4.91 ±0.12 2.15 ±0.23 2.20 ±0.30 2.19 ±0.42 2.43 ±0.42	5.00 ±0.00 2.29 ±0.40 3.32 ±0.35 3.09 ±0.46 2.75 ±0.33	5.00 ±0.00 2.13 ±0.34 3.00 ±0.40 3.41 ±0.36 2.06 ±0.28	4.98 ±0.04 2.42 ±0.30 3.10 ±0.47 3.15 ±0.30 3.05 ±0.31	4.96 ±0.05 2.16 ±0.20 2.56 ±0.28 2.81 ±0.33 2.27 ±0.27


est destinée au dépôt

Evaluation of power dissipation in a granular bed disturbed by an intruder a tool to design kneading operations

Julien Lehuen, Jean-Yves Delenne, Abdelkrim Sadoudi, Duri Agnès, Thierry

Ruiz

To cite this version:

Julien Lehuen, Jean-Yves Delenne, Abdelkrim Sadoudi, Duri Agnès, Thierry Ruiz. Evaluation of power dissipation in a granular bed disturbed by an intruder a tool to design kneading operations. 24th International Congress of Chemical and Process Engineering, 2021, Prague, Czech Republic. hal-02929270

Evaluation of power dissipation in a granular bed disturbed by an intruder a tool to design kneading operations J. Lehuen 1,2 , J-Y. Delenne 1 , A. Sadoudi 1 , A. Duri 1 , T. Ruiz 2 1 UMR IATE 1208 CIRAD/INRA/Montpellier SupAgro/Université Montpellier -2 Place Pierre Viala, 34060 Montpellier cedex 5, France, e-mail : julien.lehuen@supagro.fr, jeanyves.delenne@inra.fr, abdelkrim.sadoudi@inra.fr, agnes.duri@inra.fr 2 UMR QualiSud 95, CIRAD/Université Montpellier -15 avenue Charles Flahault, 34093, Montpellier cedex 5, France, e-mail : thierry.ruiz@umontpellier.fr

The study of stress transmission and typologies of motions in a powder bed under low shear condition constitutes a challenging issue to achieve a monitoring of processes which involve particle mobility like kneading or agglomeration process [1]. This work deals with the quantification of the ability of a probe to ensure the particle mobilities in a granular bed disposed in the tank of a reactor. An original experiment had been developed (i) to allow the visualization of the behavior of particles in the neighborhood of an intruder in ascendant vertical motion, and (ii) to measure the drag force applied to an intruder during its extraction from an ensiled granular medium. In order to identify the force propagation in a granular ensiled medium, experiments and Discrete Element Method (DEM) simulations [2] are compared to investigate the extent of force imposed by the blade at (i) different speeds, (ii) different geometric configurations and (iii) different intruder designs.

Methods

Experimental device. The granular medium consists of a population of glass beads of 1-2 mm diameter with a small span value. The experimental set-up consists of an open glass cell container (51x100x160 mm) made of transparent glass walls. The filling is provided by a funnel and the flow rate is modulated to reach a controlled static state of the granular packing. A steel probe is especially designed (shape: flat and rectangular, hemispherical, dihedral…and dimensions in regards the particle size and the size of the cell) to fit into the glass container. It is screwed on a rod that is linked to a load cell of a force sensor. The probe is then removed from the granular packing and the drag force is measured.

Numerical simulation. An "own homemade code" in c++ using the DEM is developed to calculus the force of a similar ensiled granular medium under gravity. We increment our code with the coarse graining method and we are now able to give different fields: compactness, stress tensor components, particle velocities. We can explore all the contacts between the particles by the calculus of (i) normal and tangential forces, (ii) moment. Besides DEM gives access to the mobility of the particles (position, velocity and acceleration). With the force and the velocity, the power dissipated by friction and collisions between each particle can be determined. Power dissipation was specially calculated at each time by the component of the dissipated forces (normal and tangential) and the velocity of each particle. The evolution of the repartition of the power dissipated inside the granular bed during the vertical ascension of the intruder, constitutes a criteria to quantify the efficiency of the process (single kneading with high mobilities, strong friction between particle). With this method, influences of different velocities, shapes of the intruder and size of the cell, on the power dissipation, could be determined.

Results and discussion

Experimental and numerical simulation result are in good agreement. Drag force and global powder kinetics are fitted. These results show the establishment of different regimes inside the granular bed during the vertical ascension of the intruder. These flowing regimes are correlated with characteristic lengths and allow to distinguish different zones in the granular packing. As a consequence, the power dissipation is not homogeneous inside the granular medium. Besides the repartition of this power depend also on the shape of the intruder, and on the size beneath wall and intruder. We can now purpose an analyze to design different shape of the intruder and used it for different geometrical conditions to correspond at the configuration linked with the specific industrial applications. This study at the particle scale can now be expand at the reactor scale for process of mixing in the food, pharmaceutical industries but also in civil engineering in the case of the cement truck for example.

Figure 1 .

 1 Figure 1. a) drag stress for different velocities of an intruder b) vertical velocity field obtained by coarse graining c) power dissipation inside granular medium d) vertical stress field obtained by coarse graining


Mehdi Benallegue

Rafael Cisneros

Abdelaziz Benallegue

Yacine Chitour

Mitsuharu Morisawa

Fumio Kanehiro

Lyapunov-Stable Orientation Estimator for Humanoid Robots

est destinée au dépôt

Lyapunov-stable orientation estimator for humanoid robots

Mehdi Benallegue, Rafael Cisneros, Abdelaziz Benallegue, Yacine Chitour, Mitsuharu Morisawa, Fumio Kanehiro Abstract-In this paper, we present an observation scheme, with proven Lyapunov stability, for estimating a humanoid's floating base orientation. The idea is to use velocity aided attitude estimation, which requires to know the velocity of the system. This velocity can be obtained by taking into account the kinematic data provided by contact information with the environment and using the IMU and joint encoders. We demonstrate how this operation can be used in the case of a fixed or a moving contact, allowing it to be employed for locomotion. We show how to use this velocity estimation within a selected two-stage state tilt estimator: (i) the first which has a global and quick convergence (ii) and the second which has smooth and robust dynamics. We provide new specific proofs of almost global Lyapunov asymptotic stability and local exponential convergence for this observer. Finally, we assess its performance by employing a comparative simulation and by using it within a closed-loop stabilization scheme for HRP-5P and HRP-2KAI robots performing whole-body kinematic tasks and locomotion.

I. INTRODUCTION

A humanoid robot is a system that is intended to move through the environment. However, this robot does not possess a dedicated actuator for this displacement. Instead, it has to use its joint actuators to generate contact forces with the environment, and the reaction forces create a centerof-mass acceleration. However, these contacts are limited by nature, mainly because of their unilaterality constraint. This constraint implies that the feet can only push on the environment and not pull. The consequence on a flat plane is that the center of pressure of the contact has to lie strictly inside of the convex hull of the support area for a motion to be physically feasible [START_REF] Wieber | Modeling and Control of Legged Robots[END_REF]. The estimation of the position of this center of pressure and the prediction of its dynamics require a precise estimation of the robot kinematics relative to the support area, especially the orientation of a root link referred to as the floating base. In an ideal case, this orientation could be obtained from the contact information and the joint encoders data. However, in practice, this is not guaranteed because of imprecise environment models and the presence of structural compliance in the robot or the contact surfaces. Therefore, this orientation is usually obtained using inertial measurement units (IMUs).

M. Benallegue, R. Cisneros, A. Benallegue, M. Morisawa, F. Kanehiro are with the CNRS-AIST JRL (Joint Robotics Laboratory), IRL, National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba, Japan. A. Benallegue is also Laboratoire d'Ingénierie des Systèmes de Versailles, France, and Y. Chitour is with Université Paris-Saclay, CentraleSupélec, CNRS, France. Email mehdi.benallegue@aist.go.jp, rafael.cisneros@aist.go.jp, benalleg@lisv.uvsq.fr, yacine.chitour@l2s.centralesupelec.fr, m.morisawa@aist.go.jp, f-kanehiro@aist.go.jp 1. Model of the anchor point definition in the case of non fixed anchor. In the top an example of single support phase and in the bottom a double support. The origin of the control frame C is located at the position of the anchor point. The origin of the world frame W can be located anywhere. The rotation R C between the world frame W and the control frame C is not available with encoders.

R c R c L L R L R L Fig.
These sensors traditionally include gyrometers (also called gyroscopes) and accelerometers measuring respectively, in the sensor frame, angular velocities, and linear accelerations including the gravitational one. Magnetometers cannot be used due to the presence of ferromagnetic materials and powerful electric actuation. Therefore, we can only partially reconstruct the orientation, since the yaw angle is not observable. Only the inclination with regard to the gravitational field, referred to as tilt, can be extracted. However, even with perfect measurements, it is not possible to distinguish the gravitational acceleration from the linear one using solely the sensors signals.

Several solutions to obtain this tilt have been developed. Most of them consider that the linear accelerations are negligible and use traditional estimators such as Kalman Filters [START_REF] Kajita | Biped walking stabilization based on linear inverted pendulum tracking[END_REF] or complementary filtering [START_REF] Allgeuer | Robust sensor fusion for robot attitude estimation[END_REF], [START_REF] Mao | Eliminating drift of the head gesture reference to enhance google glass-based control of an nao humanoid robot[END_REF]. However, a humanoid robot is a system subject to large linear accelerations, mostly due to impacts in high frequencies, structural oscillations in middle frequencies, and bipedal locomotion in low frequencies. Therefore, there is a need for an accurate tilt estimation that would not be sensitive to these accelerations.

One important progress has been made in this topic by taking explicitly into account the contact with the environment [START_REF] Bloesch | State Estimation for Legged Robots -Consistent Fusion of Leg Kinematics and {IMU}[END_REF]- [START_REF] Wawrzyński | Robust estimation of walking robots velocity and tilt using proprioceptive sensors data fusion[END_REF]. Indeed, this information provides an anchor between the robot's kinematics and the world frame and allows to algebraically link the rotations with the translations of the floating base. This model has also been extended to take into account multiple pivot points due to flexibilities [START_REF] Vigne | Estimation of multiple flexibilities of an articulated system using inertial measurements[END_REF], to be merged with other sources of measurements such as LIDAR [START_REF] Nobili | Heterogeneous sensor fusion for accurate state estimation of dynamic legged robots[END_REF], [START_REF] Maurice F Fallon | Drift-free humanoid state estimation fusing kinematic, inertial and lidar sensing[END_REF] or to take into account the dynamics in order to improve estimation accuracy [START_REF] Kuindersma | Optimizationbased locomotion planning, estimation, and control design for the atlas humanoid robot[END_REF].

However, to our best knowledge, none of these methods has been shown to give theoretically proven global (nor almost global) stability. One important contribution was implemented, in an open-source library, on the 2KHz loop of Cassie biped robot [START_REF] Hartley | Contact-aided invariant extended kalman filtering for robot state estimation[END_REF], [START_REF] Hartley | Contactaided invariant extended kalman filtering for legged robot state estimation[END_REF] and was based on invariant Kalman Filtering [START_REF] Barrau | The invariant extended kalman filter as a stable observer[END_REF] which allowed the estimation to have local stability. In our previous paper we have shown that if we model the robot as a simple inverted pendulum with a fixed contact point, we can use an almost global Lyapunov stable state estimator for the tilt [START_REF] Benallegue | Tilt estimator for 3D nonrigid pendulum based on a tri-axial accelerometer and gyrometer[END_REF]. The principle of the solution is to take advantage of the kinematic coupling due to the contact to build algebraically a measurement of the linear velocity expressed in the sensor frame. This additional measurement allows us to separate the gravitational acceleration and the linear acceleration within a simple state observer [START_REF] Hua | Stability analysis of velocity-aided attitude observers for accelerated vehicles[END_REF]. However, the limitation of this method to a fixed pivot point makes either the solution constrained on fixed contacts or it creates discontinuities if the contact is translated.

In this paper we extend our previous work by presenting three novel contributions:

1) The new pendulum model allows us to have a nonrigid kinematic chain with a moving pivot point. This model then can be used for humanoid robots during locomotion and does not create discontinuities in the estimated velocities. 2) We use a new state estimator with better convergence properties. This estimator is adapted from our prior work where a more general set of estimators has been presented [START_REF] Benallegue | Velocity-aided IMU-based Attitude Estimation[END_REF] 1 . They are based on a general two-stage observer scheme, one having efficient exponential global convergence driving the second one which is more robust to noise and disturbance errors. However, in [START_REF] Benallegue | Velocity-aided IMU-based Attitude Estimation[END_REF], only theoretical developments were made, assuming the presence of an additional velocity measurement and without studying their application to humanoid robots. In this paper we extract one of these estimators, adapt it to the case of humanoid robots, and for self-containment of the paper, we present a new specific proof of almost global asymptotic stability and local exponential convergence. 3) We run simulations of a biped in a variable time-step dynamical simulator to increase the accuracy of the simulation and test the sensitivity of the observer to noise and modeling error of the non-rigid pendulum model and to compare it to other methods. We also implement the estimator on a large scale humanoid robot and perform successful closed-loop stabilization while bending and walking. The next section will formally state the problem and the following sections will present the aforementioned contributions in the same order.

II. PROBLEM STATEMENT A. Frames and measurements

We denote W the world frame and L the local frame of the sensor. To simplify notations the symbol for the world frame W will be omitted. The attitude estimation has to rely on an IMU consisting in an three-axial accelerometer and a gyrometer. The accelerometer provides the sum of the gravitational field and the linear acceleration of the sensor and the gyrometer provides the angular velocity ω L of the IMU. Both signals are expressed in the sensor local frame L. The measurements are then

y g ∆ =ω L , (1)
y a ∆ =R T L pL + g 0 R T L e z , (2)
where y g and y a are the gyrometer and accelerometer measurements given by the IMU. R L , p L are the orientation and the position of the IMU in the world frame. And g 0 and e z are respectively the standard gravity constant, and a unit vector collinear with the gravitational field and directed upward. Finally ω L is the angular velocity of the sensor expressed in L such that

ṘL = R L S(ω L ), (3)
where S is the skew-symmetric operator allowing to perform cross products. The tilt can be entirely obtained through the estimation of the gravitational field in L. In other words, we need to build an estimation of R T L e z .

B. Linear velocity

The estimator that we use in this paper requires an additional measurement to the IMU, providing the vector of the linear velocity of the local frame L in the world W expressed in L [START_REF] Benallegue | Velocity-aided IMU-based Attitude Estimation[END_REF]. In other words, we need an estimation of v L = R T L ṗL . Note that with this measurement we can rewrite the model of the augmented IMU as

y v ∆ =v L (4)
y g ∆ =ω L , (5)
y a ∆ =S(ω L )v L + vL + g 0 R T L e z . ( 6
)
where y v is the measurement of v L . This measurement could be obtained with different kinds of sensors such as Doppler effect sensors or by the derivation of the position given by absolute kinematic estimators such as SLAM or LIDAR. However, most humanoid robots are not equipped with such sensors. Nevertheless, we show in the next section that this value can be obtained only with contact information.

III. RECONSTRUCTING THE LINEAR VELOCITY

A. The model of the contact as an anchor point

The main specificity of legged robots is that the locomotion relies on contacts with the environment that generate external forces. As we said before, the contacts are not perfectly rigid and deformations may occur. However, the contact points can be approximated to a fixed point in the environment as long as no contact slipping occurs, and this approximation is likely accurate to the order of the centimeter. Since the assumption is made, the presence of an anchor point in the environment creates a coupling between the orientation of the IMU and its position, which gives also a relation between the angular and the linear velocities. We show hereinafter how this coupling allows us to provide the signal of linear velocity.

Assume we have a robot equipped with an IMU, such that there is a perfectly-known joint-connected kinematic chain of rigid limbs linking the IMU to the feet of the robot. There are joint encoders that can provide the position and the velocity of each of these joints.

Let's take the case where contact occurs between a foot and the environment. We assume that the contact point is at a constant position in the environment, and we set this position at the origin of the world frame W without loss of generality. If the contact and the robot were completely rigid and the ground inclination was known, this situation would provide us directly with the full orientation of the IMU. This could be done thanks to the joint encoders and the kinematic model of the robot. Indeed, simple direct kinematics would provide the position and the orientation of the IMU as C p L ∈ R 3 and C R L ∈ SO(3) where C stands for the control frame having its origin at the anchor point as well and an orientation attached to any specific body of the robot (e.g. feet). Therefore, in the ideal case of a rigid and known environment, the control frame C and the world frame W are identical. Furthermore, since the values of C p L and C R L are obtained through the forward kinematics, their time-derivatives are also available as 2).

C ṗL ∈ R 3 and C ω L ∈ R 3 such that C ṘL = C R L S C ω L (see the left part of Figure
However such an ideal situation does not happen and there is a difference between the real position and orientation of the sensor in the world frame p L ∈ R 3 and R L ∈ SO(3) and their values in the control frame C p L and C R L . This difference is usually a combination of translations and rotations, but in the absence of slipping, we can approximate the contact point position as being fixed. This means that the transformation between the control frame C and the world frame W is a 3D rotation R C ∈ SO(3) around the origin, which means that 2).

p L = R C C p L and R L = R C C R L (see the right part of Figure
With this model, let's develop the expression of the desired signal of linear velocity y v which is entirely reconstructable from available signals.

y v = R T L ṗL (7) = R T L R C C ṗL + R T L R C S(ω C ) C p L (8) = C R T L C ṗL + S(y g -C ω L ) C R T L C p L (9) R c R C R L L R L L L L L L

B. The model with moving anchor point

The model of a fixed anchor can provide the required signals to build the attitude in the case of contact between one steady foot and the environment. However, a humanoid robot is required to perform bipedal walking, where the contact point transitions from one foot to the other, with a double support phase during the transition. One solution to this problem is to switch the contact position from one pivot point to the other during this double support phase. However, a discontinuous change in the pivot point generates a discontinuous y v signal which introduces artificial disturbances to the estimators. Therefore, a better solution is one where the transition is performed continuously. In that case we assume that we know the velocity v a of the anchor point p a in the world frame expressed in the control frame, i.e. such that R C v a = ṗa . This velocity can be simply the derivative of the center of pressure reconstructed with the force sensor or any virtual point approximating this transition.

We attach the position of the control frame origin p C at the position of the anchor point p a to keep the transformation as a pure 3D rotation around it. Note that varying this origin positions varies the values of variables expressed in the control frame, since they are translated by p a . This model is sketched in Figure 1, where we show an example for a single and a double support. With this new frame definition, the position of the sensor in the world frame can be written as

p L = R C C p L + p C = R C C p L + p a .
Thanks to this expression we can develop the expression of the linear velocity y v as

y v = R T L ṗL (10) = R T L R C C ṗL + R T L R C S(ω C ) C p L + R T L R C v a (11) = C R T L C ṗL + S(y g -C ω L ) C R T L C p L + C R T L v a (12
) which is constructed algebraically from available signals.

Note that this formula may produce inconsistencies between the trajectory of the anchor point p a and its time derivative v a = R T C ṗa . This is because of the estimation error in the orientation of the control frame and the lack of yaw measurements, therefore we do not recommend using it for navigation. Nevertheless, what we need is only a good enough approximation of the real sensor velocity v L . We believe that the disturbance caused by this small modeling error has much less impact than the disturbance caused by the assumption of no linear acceleration made in classic approaches. Moreover, we show in the next section that this measurement is used through low-pass processing within a nonlinear complementary filter, to produce a reliable and theoretically grounded tilt estimation.

IV. TILT ESTIMATION A. State and measurement definitions

Let's define the following state variables

x 1 ∆ = v L , (13)
x 2 ∆ = R T L e z , (14)
with x 1 ∈ R 3 and x 2 ∈ S 2 , where the set S 2 ⊂ R 3 is the unit sphere centered at the origin, and defined as

S 2 = x ∈ R 3 | x = 1 .
The variable x 1 is considered known using y v , even if noisy. On the contrary, x 2 is the tilt we aim to estimate and cannot be obtained algebraically from the measurements.

From equations ( 5) and ( 13) we get

ẋ1 = -S(y g )x 1 -g 0 R T L e z + y a . (15)
This, together with the time-differentiation of x 2 , provide us with the following state dynamic equations

ẋ1 = -S(y g )x 1 -g 0 x 2 + y a , ẋ2 = -S(y g )x 2 . ( 16
)
The system ( 16) is suitable for the observer synthesis presented hereinafter.

B. Two-stage observer designed in R 3 × R 3 × S 2

and error dynamics

We define the observer for the state defined in Equations ( 13) and [START_REF] Kuindersma | Optimizationbased locomotion planning, estimation, and control design for the atlas humanoid robot[END_REF]. It is designed in R 3 × R 3 × S 2 and is given by

     ẋ1 = -S (y g ) x1 -g 0 x 2 + y a + α 1 (y v -x1 ) ẋ 2 = -S (y g ) x 2 -α2 g0 (y v -x1 ) ẋ2 = -S (y g -γS (x 2 ) x 2 ) x2 (17)
where α 1 , α 2 and γ are positive scalar gains, x1 and x2 are estimations of x 1 and x 2 respectively and x 2 is an intermediate estimation of x 2 .

Assuming perfect measurements and using the estimation errors defined as

x1 ∆ = x 1 -x1 , x 2 ∆ = x 2 -x 2 and x2 = x 2 -x2 , we get the error dynamics as      ẋ1 = -S(y g )x 1 -α 1 x1 -g 0 x 2 ẋ 2 = -S(y g )x 2 + α2 g0 x1 ẋ2 = -S(y g )x 2 + γS 2 (x 2 )x 2 -γS 2 (x 2 )x 2 (18)
To run the analysis of errors, we set

z i = R L xi and z 2 = R L x 2 . Noticing R L x2 = e z -z 2 , one gets      ż1 = -α 1 z 1 -g 0 z 2 ż 2 = α2 g0 z 1 ż2 = γS 2 (e z -z 2 ) (z 2 -z 2 ) (19)
This new dynamics is autonomous and defines a time-invariant ordinary differential equation (ODE) which simplifies the stability analysis. Indeed, if one defines the error state ξ ∆ = (z 1 , z 2 , z 2 ) and the state space

Υ := R 3 × R 3 × S ez with S ez = z ∈ R 3 | (e z -z) ∈ S 2
, one can write [START_REF] Hua | Stability analysis of velocity-aided attitude observers for accelerated vehicles[END_REF] as ξ = F (ξ) where F gathers the right-hand side of ( 19) and defines a smooth vector field on Υ .

Note that the first two lines of ( 17) constitute a separate tilt estimator defined in R 3 × R 3 . A comparable tilt estimator can be found in [START_REF] Martin | A global exponential observer for velocity-aided attitude estimation[END_REF]. We show hereinafter the convergence and the performances of this intermediate estimation. We show, afterward, that even if this intermediate estimator has very good performances, the estimator of ( 17) is an extension allowing to ensure that the tilt estimation respects the normality constraint and reject more disturbances.

1) Global exponential stability of the intermediate estimator: The two first lines of the observer given by equation ( 17) constitute an independent observer of the tilt which is designed in R 3 × R 3 where the vector x 2 ∈ R 3 is the estimation of x 2 = R T L e z . In this case the error dynamics is given by ż1

= -α 1 z 1 -g 0 z 2 ż 2 = α2 g0 z 1 ( 20
)
The dynamics is linear and autonomous. If one define the state ξ := ( α2 /g0z 1 , z 2 ) and the state space Υ := R 3 × R 3 , one can write (20) as ξ = A 1 ξ where A 1 is constant and given by

A 1 = -α 1 I 3×3 -α 2 I 3×3 I 3×3 0 (21)
This matrix allows us to state the following theorem. Theorem 1: The system (20) is globally exponentially stable with respect to the origin (0, 0).

The proof comes from the stability of the matrix A 1 . The intermediate estimation x 2 is not constrained to be normalized, despite x 2 = R T L e z being a unit vector. This constitutes a strength of this estimator, it will converge through the unit sphere from any initial value. This constitutes also two weaknesses. The first one is that a simple normalization of this estimation would risk undefined output and unbounded time-derivatives when the norm is close to zero. The second one is the sensitivity noise including the one which leads to the violation of normality constraint. The solution we propose is to add the third line of ( 17) to take profit from the convergence of x 2 by tracking it while maintaining the normality constraint of the tilt estimation.

2) Almost global asymptotic stability of the full estimator: Let's now consider the full estimator of Equation [START_REF] Barrau | The invariant extended kalman filter as a stable observer[END_REF]. We state the following theorem, Theorem 2: The time-invariant ODE defined by ( 19) is almost globally stable with respect to the origin in the following sense: there exists an open dense set Υ 0 ⊂ Υ such that, for every initial condition ξ 0 ∈ Υ 0 , the corresponding trajectory converges asymptotically to (0, 0, 0). The proof is presented in the appendix.

In addition to this asymptotic stability, we can also show good local performances around the desired equilibrium. It is indeed really important for humanoid robots to maintain precise tilt observation after the convergence of the estimation since it is involved in high gain closed-loop stabilization.

3) Local exponential stability: The good local performances of this estimator can be stated with the following theorem Theorem 3: The time-invariant ODE defined by ( 19) provides local exponential convergence of the error ξ to (0, 0, 0). The proof is presented in the appendix.

These theoretical guarantees are provided in the presence of precise measurements. But since the sensors are noisy and the velocity measurement has been produced with an approximation, we need to assess the performances and compare them with state-of-the-art approaches. Therefore in the next section, we present the implementation of this observer on a dynamical simulation and on the full-size humanoid robots HRP-5P and HRP-2KAI.

V. IMPLEMENTATION

A. Simulations

The approach validity can be assessed through a precise simulation. It is performed under the Matlab Simscape Multibody environment and uses a model of a simple 12 degrees of freedom biped robot weighing 42.6Kg with two rectangular feet. Figure 3 shows simulation screen captures. The contact with the environment is simulated with point contact at the four corners of each foot. With the effect of the weight, each contact point gets strictly inside the ground and generates a viscoelastic reaction force to simulate a compliant environment with stiffness of 10 5 N/m and damping of 10 3 Ns /m. These contact dynamics change the orientation of the contact when the forces change. The robot is set to rest configuration, with a little oscillation during convergence, and it is pushed twice, at 4 s and 14 s with 100 N and 300 N forces respectively, each lasting 0.1 s. The first one leads the robot to tilt and oscillate until convergence and the second makes the robot tip and fall. Note that this oscillation does not have a fixed point in the environment and therefore breaks assumptions of a fixed anchor point.

The signals of an IMU consisting only of a gyrometer and an accelerometer were simulated with a uniform Gaussian noise with SD of 0.02 rad/s and 0.5 m/s 2 respectively. We use the position of the ZMP as a position of the anchor point. This allows us to generate the signal of the velocity x 1 using the method of Section III-B. Three state estimators were tested: (i) our estimator in R 3 × R 3 × S 2 of Sec. IV-B, (ii) the estimator in R 3 ×S 2 presented in our previous work [START_REF] Benallegue | Tilt estimator for 3D nonrigid pendulum based on a tri-axial accelerometer and gyrometer[END_REF], and the estimator of Hua 2016 from [START_REF] Hua | Stability analysis of velocity-aided attitude observers for accelerated vehicles[END_REF], which is also based on velocity-aided IMU and is comparable with our previous work, but with different dynamic properties. The parameters of the three estimators were

α 2 = k r 1 = 20, α 1 = k v 1 = k v 2
= 100 and γ = 3 and they were initialized to an error of 0.2 rad around the X axis. The simulations are very slow (about 25 minutes) because of the variable time-step duration allowing for high precision in the simulation.

The plots of Figure 4 show the result of this simulation. The beginning includes the convergence phase of the estimators. We note that our estimator in R 3 × R 3 × S 2 is the fastest to converge and the estimator of Hua 2016 is the slowest. The external forces are then applied at 4s where we can see (in the bottom) that the error has increased. This is only due to the modeling error arising from the assumption of an anchor point. However, by comparing the actual tilt with the estimation we can see that the tracking quality is good, especially the estimator in R 3 × R 3 × S 2 . Interestingly the estimator of Hua 2016 has good robustness to this modeling error compared to the estimator in R 3 × S 2 . This could be specifically due to the slowness of Hua's estimator providing it with stronger filtering of disturbances than our previous method. But the estimation of our method in R 3 × R 3 × S 2 is still better. However, the most interesting part is the fall induced at 14s where all estimators have rather good tracking thanks to their convergence properties. Nevertheless, If we focus on the error on the bottom, we can see that the previous estimators have an increased error while our estimator did not suffer a lot from this important disturbance. It is interesting to observe also, in the estimation on the Y axis, the noise sensitivity of the methods. We can see that our method in R 3 × R 3 × S 2 has the best filtering, followed by Hua 2016 and then our previous method in R 3 These results confirm the good behavior of the estimator in R 3 × R 3 × S 2 which is the fastest to converge without being too sensitive to noises and disturbances.

B. Open-source implementation on HRP-5P and HRP-2KAI

The estimator in R 3 × R 3 × S 2 , coupled with the linear velocity reconstruction exposed in the previous section has been implemented in the control framework of two robots, HRP-5P [START_REF] Kaneko | Humanoid robot HRP-5P: An electrically actuated humanoid robot with high-power and wide-range joints[END_REF] and HRP-2KAI [START_REF] Kaneko | Humanoid robot HRP-2Kai -improvement of HRP-2 towards disaster response tasks[END_REF]. Both are fullscale humanoid robots of 182 cm, 101 kg and 170cm, 65Kg respectively, and are equipped with 37 and 32 degrees of freedom, developed at the National Institute of Advanced Industrial Science and Technology (AIST). Both robots are controlled in a real-time framework based on Open-RTM middleware developed in AIST as well. Nevertheless, the implementation of the tilt estimator has rather been performed as a standalone state observation open-source C++ library [START_REF]State-estimator -open-source c++ library source code[END_REF], called simply tilt estimator.

The tilt estimation is used in the real-time balance control of the robots described in [START_REF] Morisawa | Biped locomotion control for uneven terrain with narrow support region[END_REF], which requires an estimation of the position and the velocity of the center of mass. To do that, let us assume that the positions of the feet are known. They can be for example provided by a footstep planner. Note that although a precise location is required for the navigation, it is not necessary for the balance control. From the feet position, the control frame C is built. The anchor p a , at the origin of C, can be described by a weighted ratio of the feet locations, the weights being the desired vertical forces at each foot. This definition allows us to keep a continuous trajectory of the origin while successfully transiting from one foot to another, and to know the velocity of this origin v a . The position C p L and orientation C R L of the IMU in C are known from the model and the encoders, and the linear and angular velocities C ṗL and C ω L can be estimated by using finite differences. This allows implementing the estimator as described in Section III-B. The desired yaw coming from the bipedal locomotion planner is combined with the estimated tilt using the TRIAD method [START_REF] Shuster | Three-axis attitude determination from vector observations[END_REF], where we construct a virtual magnetometer by using the desired yaw angle without disturbing the estimation of the tilt. This attitude estimation provides us with the orientation of the floating-base, for instance, the waist, and the position of the anchor point gives us the floating-base position. The angular and linear velocities of the base are obtained with the gyrometer and the anchor point. Having estimated the state of the floating-base it is straightforward to estimate the state of the center of mass. This is the necessary feedback that is required to implement the balance control.

The balance control is realized by a PID controller of the error between the desired and the estimated Divergent Component of Motion (DCM). The DCM is calculated as ξ = c + 1 ρc ċ, where c and ċ are the CoM position and its velocity, and ρ c is the natural frequency of the dynamics of the equivalent inverted pendulum [START_REF] Englsberger | Three-dimensional bipedal walking control using divergent component of motion[END_REF]. The DCM-based controller modifies the reference center of pressure, from which the reference forces and moments of the feet are computed. Then, a damping control, described in [START_REF] Kajita | Biped walking stabilization based on linear inverted pendulum tracking[END_REF], modifies the relative position of the feet and their orientation to realize the foot reference force and moment. Finally, the joint angles are calculated by using the prioritized inverse kinematics scheme proposed by [START_REF] Kanoun | Kinematic control of redundant manipulators: Generalizing the task-priority framework to inequality task[END_REF] and sent to the joint servos. Figure 5 illustrates the previous description.

Since there is no ground truth value for the real robot orientation, we test the performance of the estimator by assessing the stability in two experiments. The first one consists of leaning the chest forward by 25 degrees and coming back to the vertical, all in 2 seconds. This allows us to check the reactivity of the estimator to fast tilt changes, similarly to what happens in the falling simulation above. The second test is walk forward, this will assess the behavior concerning impacts and modifications of the contact point. All these experiments were successful and the balance of both robots could be maintained. Screen captures of these tests are available in Figure 6. For the case of HRP-5P, we show, in Figure [START_REF] Eljaik | Multimodal sensor fusion for foot state estimation in bipedal robots using the extended Kalman filter[END_REF], the measurements of the IMU together with the tilt estimation during each experiment. There we can see the tracking of the inclination of the chest and the quality of estimation during walking. We can evaluate the gait impacts by looking at the accelerometer measurements. The accelerometer measurements are on the top, the gyrometer measurements on the middle, and the tilt estimation on the bottom. The Z axis has been removed from the tilt estimation for legibility and can be reconstructed knowing that it is the positive solution to the vector normality constraint.

VI. CONCLUSION

We have presented a Lyapunov stable tilt estimator that can be applied to humanoid robots. We have shown how the presence of an anchor position, even if it is moving, can be exploited to algebraically produce approximations of the linear velocity of the sensor. This new signal allows us to use estimators dedicated to velocity-aided IMUs. We have shown the theoretical guarantees of almost global asymptotic stability and local exponential convergence of the presented two-stage estimator. Finally, we have studied the performance of the estimator on a dynamic simulator of a simple biped and have shown the stability of real robots when used within a closed-loop balance control framework. APPENDIX A. Proof of theorem 2

Proof: We define the positive-definite differentiable function V : Υ → R + radially unbounded over Υ

V ∆ = α 2 8α 1 g 0 z 1 2 + 1 8α 1 g 0 α 1 z 1 + g 0 z 2 2 + α 1 g 0 4α 2 z 2 2 + 1 2γ z 2 2 . ( 22
)
Using the error dynamics given by ( 18), the time derivative of V is then given by It means that V ≤ 0 and more specifically V < 0 if (z 1 , z 2 , z 2 ) is not an equilibrium point.

V = - α 2 1 4g 0 z T 1 z 1 - g 0 4 z T 2 z 2 + 1 γ z T 2 ż2 , (23)
The linearized system around the equilibrium (0, 0, 2e z ) is given by the following dynamics where I 3×3 is the 3 × 3 identity matrix. The characteristic polynomial of the matrix A 2 is given by

Ẋ = A 2 X, (25)
P (λ) = λ 2 + α 1 λ + α 2 3 λ (λ -γ) 2 . ( 27
)
This polynomial has the following double real positive root

λ = γ > 0, (28)
It can be verified that the linearized system at (0, 0, 2e z ) admits positive real eigenvalues. Hence the origin is almost globally asymptotically stable. This completes the proof of the theorem.

B. Proof of Theorem 3

Proof: The estimator has a unique stable equilibrium point defined by ξ = (0, 0, 0)

T . The linearization of the error dynamics of ( 19) around this point is defined in the tangent space to the state space Υ , which happens to be T (Υ ) =

R 3 × R 3 × R 2 .
We define the state ξ p ∆ = (z 1 , z 2 , z 2p ) where z 2p ∈ R 2 is the orthogonal projection of z 2 on the horizontal plane. The linearized dynamics is given by ξp = A 3 ξ p with A 3 being a constant matrix having the form

A 3 =   -α 1 I 3×3 -g 0 I 3×3 0 α2 g0 I 3×3 0 0 0 I 2×3 -γI 2×2   (29)
where I 2×3 = 1 0 0 0 1 0 . The matrix A 3 is Hurwitz giving the stability of the linearized dynamics. QED.

Fig. 2 .

 2 Fig. 2. Two different configurations of the robot model, with fixed contact point, represented on the top and on the bottom. The IMU of each robot is represented with the circle in the head of the robot. For both top and bottom the left figure represents the position and orientation of the robot in the control frame C and the right one represents the same robot in the world frame with the unknown rotation R C . Note that the rotation R C on the top and on the bottom is the same.

Fig. 3 .

 3 Fig. 3. Screenshots of the simulation of the simple biped. The first imageshowing it in rest position, the second one after the small push and the third one after the big push.

Fig. 4 .

 4 Fig. 4. Plot showing the evolution of estimations from three tilt observers for the simulated biped. The top three figures show each component of the tilt estimation x2 compared to the real value. The bottom shows the error angle for the three estimators.

Fig. 5 .

 5 Fig. 5. A simplified diagram of the bipedal locomotion control framework. Elements in red are components related to the estimation. Elements in yellow are components used to generate desired values. Elements in blue are components related to the control. Elements in green belong to the robot.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Screen capture of the robots (HRP-5P on top and HRP-2KAI on the bottom) performing the test motions. The two leftmost images show the leaning task and the remaining images show the walking motion.

z 2 z 2 - 2 + w z 2

 2222 + z T 2 S 2 (e z )z 2 -z T 2 S(e z )S(e z -z 2 )z 2(24)If we define the new vector w = S(e z )z 2 , we can write w T w + w T S(e z -z 2 )z 2We define the vector = [ z 1 z 2 w ] T , givingV ≤ -T H

with X = (z 1 )T(z 2 )T(z 2 -

 122 2e z ) T Tand A 2 being a constant matrix having the formA 2 =   -α 1 I 3×3 -g 0 I 3×3 0 α2 g0 I 3×3 0 0 0 -γS 2 (e z ) -γS 2 (e z )

  × S 2 .

	X
	0	2	4	6	8	10	12	14	16	18	20
						Time (s)
	Y
	0	2	4	6	8	10	12	14	16	18	20
						Time (s)
	Z
	0	2	4	6	8	10	12	14	16	18	20
						Time (s)
	0	2	4	6	8	10	12	14	16	18	20
						Time (s)

The paper is currently still under review, but the manuscript is available on HAL


F Thiesset

L Danaila

Focus on Fluids

Keywords: Turbulent Flows, Wakes/Jets, Turbulence theory

The theory of Kolmogorov, enunciated for very large Reynolds numbers, has progressively been shown to be inoperative for characterizing flows of practical relevance. Yet, in a recent study by Alves Portela et al. (2020), the turbulence statistics in the very near wake of a square prism at modest Reynolds numbers, reveal a significant portion of scales complying with a cascade of Kolmogorov-type. By resorting to a generalized version of the Kármán-Howarth-Kolmogorov equation, this intriguing observation is shown to be an illusion, hiding a measurable influence of coherent structures and statistical inhomogeneity. This striking conclusion highlights that a complete statistical theory of turbulence cannot dispense with the influence of large scales, possibly coherent, motions.

Introduction

Despite its ubiquity, turbulence continues to challenge mathematicians, physicists, and engineers because of the complex hierarchy of eddies of fluid motion, which interact in a puzzling manner. Kolmogorov (1941a) was the first to hypothesize that for asymptotically large Reynolds numbers, a clear separation is expected between the processes acting at large, intermediate and dissipative scales. The information conveyed by the large scales is likely to decline after the first non-linear interactions so that there should exist a scale beyond which turbulence becomes universal and therefore locally isotropic. Under the same assumption of sufficiently large Reynolds numbers, Kolmogorov (1941b) carried out a statistical analysis of the velocity field by considering the turbulent kinetic energy probed at fixed scale defined as the distance between two points. His analysis showed that, for locally isotropic turbulence, the energy transfer between scales in the intermediate range is constant and proportional to the mean energy dissipation rate. This result, known as the 4/5th law, is often considered the only exact (yet asymptotic) theory of turbulence.

However, the Reynolds numbers achievable in laboratory experiments are far smaller than would be required for Kolmogorov's 4/5th law to apply. In such situations, scale separation does not occur. The effect of finite Reynolds number, i.e. the combined influence of small-scale (viscous) and large-scale effects, has been accounted for in the generalized Kármán-Howarth-Kolmogorov equation, which is derived from first principles without the need to invoke any particular symmetry of the flow (e.g. [START_REF] Danaila | Turbulent energy scale budget equations in a fully developed channel flow[END_REF][START_REF] Hill | Equations relating structure functions of all orders[END_REF]. Extensive research on this theory in the last three decades (see, e.g., [START_REF] Qian | Inertial range and the finite reynolds number effect of turbulence[END_REF][START_REF] Danaila | An overview of the effect of large scale inhomogeneities on small-scale turbulence[END_REF][START_REF] Marati | Energy cascade and spatial fluxes in wall turbulence[END_REF][START_REF] Antonia | Approach to the 4/5 law in homogeneous isotropic turbulence[END_REF][START_REF] Valente | The energy cascade in grid-generated non-equilibrium decaying turbulence[END_REF] has demonstrated the predominant role of non-stationarity in decaying flows, statistical inhomogeneities, turbulent kinetic energy production, turbulent/pressure diffusion, and the general concept of forcing.

To complicate matters further, the large scales sometimes reveal an identifiable and regular behaviour related to the type of hydrodynamical instability that has led to their creation. These are generally referred to as coherent structures (see for example [START_REF] Reynolds | The mechanics of an organised wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments[END_REF]. They occur within most flows but they remain poorly understood because coherent motion possesses non-universal characteristics associated with initial and boundary conditions. [START_REF] Townsend | Early days of turbulence research in Cambridge[END_REF] synthesized this conundrum in the following terms "at this distance [50's] it may seem to be nearly the end of purely statistical approaches before concepts of eddy structure began to assume importance". Eighty years since Kolmogorov's theory, there remain open questions such as the way coherent motion interacts with randomly fluctuating turbulent eddies and how the former can be incorporated into a statistical theory of turbulence. Does the organized motion influence the cascade process? Since coherent structures generally emerge from velocity shear, does the cascade process become sensitive to inhomogeneities? The recent work by Alves Portela et al.

(2020) sheds new light on these long standing questions and proves that the effect of coherent motion needs to be explicitly accounted for.

Overview

The turbulence community owes a great deal to [START_REF] Reynolds | The mechanics of an organised wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments[END_REF] who first highlighted the dynamical influence of coherent structures or 'organized waves'. They came up with the astute idea of decomposing the fluctuating field into mean, coherent and random components and found that a portion of the energy conveyed by the organized structures could be transferred to random motion. However, they considered one-point statistics and thus could not gain insight into the scales affected by the coherent motion.

This issue has been addressed more recently by [START_REF] Thiesset | Scale-by-scale turbulent energy budget in the intermediate wake of two-dimensional generators[END_REF] who showed that coherent motion may act to force the cascade between scales. They found that in contrast with classical decaying turbulence, the transfer of energy in wakes behind twodimensional generators was systematically enhanced, and therefore closer to the asymptotic Kolmogorov scaling. This observation was further confirmed by invoking structure functions conditioned by the phase of the coherent motion [START_REF] Thiesset | Dynamical interactions between the coherent motion and small scales in a cylinder wake[END_REF], which is the two-point extension of the one-point phase-averaged statistics introduced by [START_REF] Reynolds | The mechanics of an organised wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments[END_REF]. Within this framework they showed that both the energy distribution and its transfer at a given scale are modulated (in amplitude) by the time-periodic dynamics of the coherent motion. Using two-point statistics, they derived the transport equations for the second-order structure functions of both coherent and random components of the velocity field to characterise the turbulent cascade process in the presence of coherent structures. They found that at a fixed location in the flow, energy transfer between scales is significantly enhanced when coherent motion is present, and diminished in its absence.

Alves Portela et al. (2020) apply a similar methodology to derive fully inhomogeneous and anisotropic Kármán-Howarth-Kolmogorov equations for the kinetic energy in the presence of quasi-periodic coherent structures, accounting for inhomogeneities of the mean flow. The novelty of their approach is to use this theoretical framework to explore the very near-field of the wake behind a prism where the scale/space features remain largely unknown, in contrast with intermediate and far-field flows. They use direct numerical simulations (DNS) to calculate the contributions of all processes including turbulence decay and production, pressure/turbulent-diffusion and the forcing due to coherent structures.

The total (coherent + random) energy transferred is restored, because all terms are computed from the first principles. Alves Portela et al. (2020) then discovered that there exists a reasonable portion of scales which surprisingly comply with the Kolmogorov constant energy transfer in the cascade. The authors qualify this result as "a clearly non-Kolmogorov yet Kolmogorov sounding [cascade]". Indeed, Kolmogorov's theory should not apply for the Reynolds numbers investigated by Alves [START_REF] Alves Portela | The role of coherent structures and inhomogeneity in near-field interscale turbulent energy transfers[END_REF] in the presence of the strong inhomogeneity and anisotropy of the flow. Could coherent motion and/or inhomogeneity be responsible for this apparent constant energy transfer? Using the previously introduced triple decomposition for the velocity structure functions, Alves [START_REF] Alves Portela | The role of coherent structures and inhomogeneity in near-field interscale turbulent energy transfers[END_REF] showed that the constant total energy transfer hides a contribution of the forcing associated with the coherent motion. By further decomposing the total transfer term into an inhomogeneous and homogeneous contribution, they find that inhomogeneities also contribute directly to the energy transfer in addition to manifesting themselves through the appearance of additional terms in the generalized Kármán-Howarth-Kolmogorov equations.

Impressively, they characterise how energy is transferred from the injection, through the mean velocity gradients, to the coherent motion and finally to the random fluctuations. While this general picture is believed to apply to most flows, this is the first time that this complex process is unravelled in a quantitative way, enabled uniquely by the combination of theory and DNS. Moreover, the Kármán-Howarth-Kolmogorov equation for turbulent motion in the presence of coherent structures appears to be the most appropriate framework to quantify the extent and range of scales over which each effect influences the total energy transfer, including the forcing associated with the coherent motion. The work of Alves Portela et al. (2020) provides a significant step forward in characterising the effects of inlet/boundary conditions in a complex, inhomogeneous and anisotropic flow, with the ultimate aim of predicting them at Reynolds numbers typically encountered in the laboratory. Finally, it is worth stressing that in some flows, coherent do not possess a periodic signature or multiple periodic motions can be at play (e.g. farms of wind turbines). In such situations, phase averages are precluded and classical two-point statistics (correlation and structure functions) may not remain appropriate. One is then likely to rely on some other statistical tool such as low-pass or band-pass filters, Proper-Orthogonalor Dynamic-Mode-Decomposition, or wavelets to name a few. By pursuing work in this direction, one should be able to explicitly identify the different processes (transfer, production, dissipation) in such a new space of scales, flow positions and time.

Future

  The work by Alves Portela et al. (2020) opens new avenues of investigation of the intricate interactions between coherence and turbulence. Will the organized motion continue to influence the cascade process when the inlet velocity is increased, and/or when travelling through the intermediate and far field? How are each individual velocity components and cross correlations affected? How do different directional components of energy flux perceive the coherent motion, beyond the circularly averaged transfer term considered by Alves Portela et al. (2020)? How can this information be used to build physics-informed predictive tools such as Large-Eddy-Simulation models? Answers to these open questions may be provided by generalized Kármán-Howarth-Kolmogorov equations in the presence of periodic forcing. In addition, the decomposition of the transfer term proposed by Alves Portela et al. (2020) into homogeneous and inhomogeneous contributions may help to ascertain the role of statistical inhomogeneity. The question of the uniqueness of such a decomposition which is addressed in Alves Portela et al. (2020) should be further examined.

Conflicts of interest

The authors report no conflict of interest.