Denis Gindre

Adrien Vesperini

Jean-Michel Nunzi

Hervé Leblond

Kokou D Dorkenoo

Refractive-index saturation-mediated multiple line emission in polymer thin-film distributed-feedback lasers

We report experimental and theoretical investigations of multiple laser-line emission in a distributedfeedback dye laser pumped by two coherent optical beams. We have used a Lloyd interferometer configuration with second-and third-order Bragg reflections to study the interaction between the two incident pumps in an organic thin film. We demonstrated theoretically that the number of laser emission lines can be interpreted with reference to the saturation effect in the refractive index.

Since the first papers by Kogelnik and Shank 1,2 about distributed-feedback lasers, many groups have successfully attempted, both experimentally 3,4 and theoretically, [START_REF] Eichler | Laser-Induced Dynamic Gratings[END_REF]6 the study of various configurations to obtain a laser emission. Optimized architectures have been obtained with high conversion efficiency and low divergence. 7 The Bragg reflector needed for laser emission is obtained in different configurations, some of them using a permanent grating [8][9][10][11] or a dynamic grating 12,13 obtained by means of an interferometer. If the pumps are perfectly superposed, many interactions can occur in the media. For example, Khan and Hall produced five and nine lines by using two and three pairs of pump beams, respectively. 14,15 Their experiments were performed in dye solutions, and their model, 16 based on the geometry of the system, gives the wavelengths of the laser lines as a function of the incidence angles. In this Letter we present the study of the emission spectrum of a distributed-feedback laser of a Rhodamine 6G-(Rh6G) doped polymer excited by two coherent pumps beams in a spin-coated thin film. By using the flexibility of the Lloyd mirror device, we have obtained multiple lines corresponding to second-and third-order Bragg reflections written in the thin film. In addition, we show that the content of the emission spectra can be explained by the saturation effect that occurs in the refractive index during the pumping step. A similar effect has been evidenced in a fiber Bragg grating. 17 The thin films were composed of a polymer (polymethylmethacrylate) and an organic dye (Rh6G) and were spin coated on glass substrates. The concentration of the polymer was 150 g / l and that of Rh6G solution in 1,2,3-trichloroethan solvent was 4 ϫ 10 -3 M / l. The 2 m thick film acts as a single-mode waveguide and avoids the spectral effects due to multimode propagation. 18 Figure 1 presents the experimental setup. The sample was excited by a pulsed Nd: YAG laser ( = 532 nm, ⌬t = 35 ps, 5 Hz). The temporary superposition of the two pumps was monitored by a delay line. The pair of pump beams forms a complex pattern in the thin film, and this intensity modulation creates a phase grating that induces more than two laser lines simultaneously. For each beam, the output wavelength i ͑i =1,2͒ is linked to the grating period by the Bragg condition as i = n eff p / ͑m sin i ͒, where n eff is the effective refrac- tive index, p is the pump wavelength, and m is the Bragg reflection order. This formula cannot describe the interaction between the coherent pumps, and further theoretical insight is needed to explain the spectral signature of this interaction. We consider two incident pumps, with wavelengths p at angles 1 and 2 , and the reflected beams (see Fig. 2). The electric field amplitude is

E = ͚ j=1 4 E j e ik j ជ •r ជ , ͑1͒
where k j ជ =2 / p ͑sin j , cos j ,0͒ for j = 1, 2, 3, 4, with 3 =-1 and 4 =-2 . The interference pattern is given by the intensity

I = ͉E͉ 2 = ͚ j,l=1 4 E j E l * e i͑k j ជ -k l ជ ͒• r ជ . ͑2͒
Because the emission is confined inside the polymeric film, only r ជ = ͑x ,0,0͒ and the x components k j of the wave vectors are involved. The nonlinear refractive index ⌬n, including its imaginary part that corresponds to nonlinear gain or absorption, is, in a first approximation, proportional to I. However, a saturation occurs above some value I s of the intensity. It can be accounted for by the expression ⌬n = n 2 I / ͑1+I / I s ͒.

Because of saturation, the Fourier spectrum of ⌬n involves all harmonics of the frequencies present in the spectrum of I. Since k 3,4 =-k 1,2 , the latter can be written as the sums ͑k j + k l ͒, with j, l = 1, 2, 3, 4. The nonlinear index can expanded as

⌬n = ͚ j,l=1
4 a j,l e i͑k j +k l ͒x + ͚ j,l,m,p=1 4 b j,l,m,p e i͑k j +k l +k m +k p ͒x + ͚ j,l,m,p,q,s=1 4 c j,l,m,p,q,s e i͑k j +k l +k m +k p +k q +k s ͒x + ¯, ͑3͒

where a j,l , b j,l,m,p , and c j,l,m,p,q,s are the amplitudes of the fundamental frequencies, second harmonics, and third harmonics, respectively. We now proceed to coupled-wave theory, as introduced by Kogelnik and Shank. 2 They consider a sinusoidal nonlinear index ⌬n = n 1 cos͑k 0 x͒, where n 1 is the amplitude of the modulation and k 0 characterizes its spatial frequency. This factor couples a rightpropagating wave S͑x͒exp͑ikx͒ to a left-propagating one R͑x͒exp͑-ikx͒. The parametric interaction between the two waves involves a spatial dependency exp͑2ikx͒ and must be matched with exp͑±ik 0 x͒ according to the expression of ⌬n. Therefore we must have k = k 0 / 2. In principle, a periodical perturbation of the medium generates an infinity of diffraction orders, but only two orders are phase matched around the Bragg frequency. When the second order of the Bragg reflection is considered, the coupling is due to the second harmonic terms in the Fourier expansion [Eq. ( 3)] of the nonlinear index ⌬n. The intensities of the higher harmonics are assumed to be negligible. Under these conditions, k 0 = k j + k l + k m + k p ; hence the wave vectors of the emitted light have the form k 0 = ͑k j + k l + k m + k p ͒ / 2, where j , l , m , p =1 ,2 ,3 ,4 .I no rder for k to belong to a spectral domain where gain is present, all k j must have the same sign. We get the corresponding wavelengths as = n eff 2 k = 2n eff p ͑sin j + sin l + sin m + sin p ͒ , ͑4͒

where j , l , m , and p , are either 1 or 2 . When the third order is considered, the coupling is due to the third harmonic, and the possible wave vectors are k 0 = ͑k j + k l + k m + k p + k q + k s ͒ / 2, while the corresponding wavelengths are = 2n eff p ͑sin j + sin l + sin m + sin p + sin q + sin s ͒ .

͑5͒

The Bragg reflection order is not explicitly involved in the formula. As mentioned above, we can observe laser lines only in the spectrum interval, which corresponds to the amplified spontaneous emission (ASE) of the dye. In our case, the Rh6G shows a maximum ASE near 560 nm, and the effective refractive index of the film is close to 1.5. Thus the laser lines in this interval are obtained with incident angles close to 45°, corresponding to the second Bragg order ͑m =2͒ or close to 20°for the third order ͑m =3͒. The obtained wavelengths originate in the coherent superposition for all Bragg reflection orders. We performed our experiments with m = 2 and m =3 to compare the experimental results with the theory. When the incident angle is 45°, the overlap between the reflected and the direct beam in the sample is perfect. Thus the size of the grating is optimal, and the efficiency of the spectral selectivity is maximal. The angle between the pumps is chosen smaller than 2.8°to fall in the region of highest gain of the ASE spectrum of Rh6G. Thus the pump beam is divided into two parts; a set of two mirrors allows the tuning of the incidence angle of the second pump, whereas the incidence angle of the first pump is fixed. We use s polarization for the pumps and an intensity just above the threshold to avoid multimode propagation. 19 The energy of each pump beam was approximately 100 J / pulse, and the conversion efficiency is typically 1%. The linewidth of the laser peaks is 0.5 nm, which is the limit of resolution of our spectrometer.

For m =2 ͑m =3͒, we used 1 = 44.92°, 2 = 44.47°͑ 1 = 28.16°, 2 = 27.90°͒, n eff = 1.485, and p = 532 nm. The spectrum of the output is shown as the dashed curve in Fig. 3(a) [Fig. 3(b)]. It appears that the theoretical values calculated with Eq. ( 4 Laser lines 1 and 7 for m = 3 (1 and 5 for m = 2) correspond to the emission of the gratings separately, while the lines 2-6 for m = 3 (2-4 for m = 2) were induced by the coherent interaction between the pump beams. Notice that Eqs. ( 4) and (5) give the positions of the laser lines, but not their amplitudes.

The basis of the theoretical approach is expansion (3), which requires a numerical confirmation. Indeed, in the case of a periodical pattern, the formation of harmonics due to the saturation of the nonlinear index follows from the properties of the Fourier series.

Here the situation is more complicated because the pattern is quasi-periodic with two fundamental frequencies corresponding to k 1 and k 2 . It is therefore interesting to investigate numerically which harmonics will or will not be produced. The electric field E, the intensity I, and the nonlinear index ⌬n are first computed; the spectrum of ⌬n is then obtained by using the fast Fourier transform algorithm. The coupled wavelength is related to the wave vector k as in the first equality in Eq. ( 4). We can then plot the spectrum against as shown in Fig. 3(a) (solid curve) for m = 2 and Fig. 3(b) (solid curve) for m = 3. Notice that the plots are expected to reproduce the locations of the emitted wavelengths but by no means the corresponding amplitudes. In the numerical simulation, other wavelengths appear corresponding to higherorder harmonics. They are not located between the peaks corresponding to the emission of the independent gratings 1-5 (or 1-7) but outside this range. Their amplitude is smaller and depends strongly on the value of the saturation intensity I s . It can also be reduced by taking into account some partial incoherences and amplitude dissymmetry.

In conclusion, we have shown that the multiple laser peaks observed in the distributed-feedback laser pumped by two beams are neither due to a nonlinear interaction between the fundamental lines nor to some higher-order wave coupling; they are due to the presence of higher-order spatial harmonics in the structure of the grating. Taking the higher-order harmonics into consideration allows one to predict the observed laser frequencies with great precision. The formation of the harmonics was interpreted as the consequence of the saturation of the nonlinear refractive index, as shown by the results of our numerical computations.

D. Gindre's e-mail address is denis.gindre@univangers.fr.

Fig. 1 .

 1 Fig. 1. (Color online) Experimental setup. The emitted beam is collected by an optical fiber (length 1 m) and dispersed by a 1200 lines/ mm grating spectrometer.

Fig. 2 .

 2 Fig. 2. (Color online) Configuration of incident wave vectors on the sample.

Fig. 3 .

 3 Fig. 3. (Color online) Comparison between theory and experimental results. (a) Experimental spectrum (dashed curve) obtained with 1 = 44.92°and 2 = 44.47°and plot of a part of the spectrum of the nonlinear index ⌬n (solid curve) against the corresponding coupled wavelengths with m =2. The peaks correspond to the wavelengths emitted by the laser. The circles assign an arbitrary number to the laser lines. (b) Experimental spectrum (dashed curve) and calculated spectrum (solid curve) with m = 3 obtained with 1 = 28.16°and 2 = 27.90°. The lines are superposed onto a slope pattern because of residual ASE. The imperfect spatial matching of direct and reflected beams in the sample is probably at the origin of this effect.

  ), 559.41, 560.52, 561.63, 562.74, and 563.87 [Eq. (5), 558.00, 558.79, 559.58, 560.38, 561.17, 561.98, and 562.78] almost perfectly match the experimental results: 559.40, 560.45, 561.54, 562.66, and 563.88 (557.96, 558.82, 559.56, 560.37, 561.08, 561.90, and 562.72).


Fouad Aouinti

Zeynep Sonat Baltaci

Mathieu Aubry

Aleksandre Guilbaud

Stavros Lazaris

Aouinti Fouad

Zeynep Baltaci

Sonat

Aubry Mathieu

Guilbaud Alexandre

Lazaris Stavros

Computer Vision and Historical Scientific Illustrations

Keywords: historical document analysis, IIIF, illustration detection, deep learning

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

Introduction and context

Illustrations and their evolution in the scientific corpora of the Middle Ages and modern Western cultures have only been partially studied. More generally, the role of image in the construction and dissemination of scientific knowledge raises complex questions that remain historically delicate to grasp and for which adapted analysis tools are lacking. To fill this gap, we aim to develop automated methods for extracting illustrations and then analyzing their similarities, which will lead to the constitution of an iconographic series that can be interpreted by historians. One major difficulty lies in the heterogeneity characterizing these corpora and the impossibility of carrying out large-scale annotations for every specific task. We thus develop approaches to handle the heterogeneity of historical data with limited annotation effort.

Various deep learning methods have been developed for the visual analysis of historical documents. In this paper, we focus on our dataset creation and illustrations extraction. We introduce the digital platform developed to allow historians to leverage deep learning methods, improve them by collaborative annotation, and historically analyze the results.

Dataset and illustration extraction

We started by building a platform to enable historians to upload 54 manuscripts and 704 printed volumes of interest in varying formats (IIIF manifests, JPEG, PDF), and the associated metadata. This resulted in 405,594 document images. To ease the work of annotating illustrations, we developed a two-stage workflow (Fig. 1). First, we exploited YOLOv5s (Jocher et al., 2020) trained on SynDoc [START_REF] Monnier | docExtractor: An off-the-shelf historical document element extraction[END_REF] and predicted illustrations bounding-boxed on a subset of 4451 representative images. Our platform then enabled historians to easily verify and if necessary use SAS1 to correct Tab. 1: Performance of YOLOv5s in precision (P), recall (R), F-score and average precision (AP).

and complete these detections. We split this curated 'VHS' data into train (60%), validation (20%), and test (20%) sets and use it to fine-tune YOLOv5s. We will release this data upon publication. We tested different pre-training and fine-tuning strategies and obtained a very significant boost over the network trained on synthetic data (see Table 1). We used the network pre-trained on SynDoc and fine-tuned on VHS to obtain 235,198 illustrations (including false positives) from our complete corpora.

Discussion

Based on this huge high-quality illustrations dataset, we work on developing similarity search algorithms and interfaces. We already obtained high-quality results using (Xi et al., 2021). Our next step is to equip the platform with tools that will enable historians to easily analyze and annotate the obtained image clusters, to assess their historical interest, and to refine the developed similarity search tools, in a spirit similar to the one we have presented for illustration detection.

1stFig. 1 :

 1 Fig. 1: Overview of our pipeline to refine illustrations detections in scientific historical documents.

https://github.com/glenrobson/SimpleAnnotationServer

Acknowledgements

This work was supported by the ANR (ANR project VHS ANR-21-CE38-0008). MA and SB were supported by ERC project DISCOVER funded by the European Union's Horizon Europe Research and Innovation programme under grant agreement No. 101076028. Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible for them.


Francesco Ruscelli
email: francesco.ruscelli@iit.it

Arturo Laurenzi
email: arturo.laurenzi@iit.it

Enrico Mingo Hoffman

Nikos G Tsagarakis
email: nikos.tsagarakis@iit.it

A Fail-Safe Semi-Centralized Impedance Controller: Validation on a Parallel Kinematics Ankle

This paper proposes the implementation of an impedance controller on the ankle level of COMAN+, a robot with parallel kinematics ankles actuated by a dual four-bar mechanism. The main contribution of the work is a realization of said control scheme that grants a less abrupt and safer robot response in case of system failures, that would cause the local joint torque controllers to lose their torque reference inputs. In particular, we propose a semi-centralized impedance control implementation which eliminates the instability of the pure joint torque control schemes used in the classical fully centralized methods when torque reference interruptions occur. Finally, we present experimental results, proving the effectiveness of our method and demonstrating how it ensures a safer behaviour compared to a fully centralized impedance control implementation when the communication to the ankle joints is interrupted. This paper is a follow-up work of [1], which presented and analyzed the parallel kinematics ankles.

I. INTRODUCTION

Robust locomotion of humanoid robots in complex environments is still an unsolved topic. Many solutions have been proposed, both on a higher behavioural level, e.g. footstep planning [START_REF] Herdt | Online Walking Motion Generation with Automatic Footstep Placement[END_REF], [START_REF] Chestnutt | Footstep Planning for the Honda ASIMO Humanoid[END_REF] or balance recovery [START_REF] Stephens | Humanoid push recovery[END_REF], and on a lower control level, such as active (admittance controller [START_REF] Zhou | A passivity based compliance stabilizer for humanoid robots[END_REF], [START_REF] Li | A Passivity Based Admittance Control for Stabilizing the Compliant Humanoid COMAN[END_REF], joint impedance control [START_REF] Hoffman | Upper limb compliant strategy exploiting external physical constraints for humanoid fall avoidance[END_REF], Virtual Model Control [START_REF] Pratt | Virtual Model Control: An Intuitive Approach for Bipedal Locomotion[END_REF]), and passive compliance (Series Elastic Actuators [START_REF] Pratt | Series elastic actuators[END_REF], [START_REF] Tsagarakis | COMpliant huMANoid COMAN: Optimal joint stiffness tuning for modal frequency control[END_REF], [START_REF] Tsagarakis | WALK-MAN: A High-Performance Humanoid Platform for Realistic Environments[END_REF]). Furthermore, in implementing any control strategy particular emphasis should be put on the robot reliability: reducing, in the presence of communication interruptions between the centralized high layer controllers and the local actuation regulators, the risk of abrupt and unwanted reactions that could damage the robot and the environment or create unsafe conditions in human-robot interaction tasks. Similarly to the intrinsic elasticity of human muscles, which allows small adaptations to external obstacles and partial energy absorption from impacts [START_REF] Pratt | Virtual Model Control: An Intuitive Approach for Bipedal Locomotion[END_REF], an effective solution to improve the interaction skills of robotics systems experiencing frequent dynamic interactions with the environment is the employment of impedance control schemes [START_REF] Hyeon | Hybrid control for biped robots using impedance control and computed-torque control[END_REF], [START_REF] Lim | Balance and impedance control for biped humanoid robot locomotion[END_REF]. Interacting effectively with the environment is very relevant to humanoids and in general to legged locomotion where the ability to conform and adapt to ground disturbances and impacts is a fundamental prerequisite for achieving stable locomotion. However, one of the limitations of the classical impedance control is its sensitivity to interruptions of torque reference signals for the local joint torque controllers. Due to the inherent instability of a pure joint torque regulator, if the communication between the centralized impedance control and the joint torque controllers is interrupted, unwanted torques are commanded to the motors. This could compromise the overall robot stability, creating dangerous situations for humans interacting with the robot and increasing the risk of damages both to the physical system and the environment.

In this work, we propose a method to realize a desired impedance using a partially centralized scheme that reduces abrupt motions when the torque reference of the actuators are interrupted. Specifically, our work is divided into two parts, both oriented towards the goal of devising an impedance controller capable of partially rejecting external perturbation while being more stable against interruptions in the actuation reference signals:

• the implementation of an impedance controller, on the humanoid COMAN+ (shown in Figure 1); • the design of the Semi-Centralized Impedance Controller (SCIC), a specific formulation of the impedance controller based on a semi-centralized scheme.

The proposed control was implemented on the ankle joint, which is based on a parallel kinematics mechanism [START_REF] Zhou | On the Comprehensive Kinematics Analysis of a Humanoid Parallel Ankle Mechanism[END_REF] and actuated by two motors through a dual four-bar transmission mechanism. Since the foot is the joint directly interacting with the terrain and hence subjected to impact forces and disturbances from the ground, compliant ankles could enhance the robot behaviour on uneven grounds, adapting its step to the terrain and rejecting minor disturbances on the joint.

The main motivation of this work is the realization and validation of the aforementioned controller so as to ensure the system stability and safety in the presence of failures, i.e. communication interruptions between the on-board centralized computer and actuation DSPs controllers. Such events could lead to unwanted torques commanded abruptly to the robot, leading to hardware damages or sudden motions, especially dangerous in human-robot interaction tasks. The Semi-Centralized Impedance Controller (SCIC) takes advantage of the local joint impedance controller implemented at each motor's DSP, sending the desired torque as two different components: the first depends only on the reference and the sensed position of the motor, while the second term is a bounded torque sent in feed-forward. This allows the robot, if the communication is interrupted, to keep the last motor position reference sent, while the feed-forward term becomes a bounded disturbance which is easily rejected.

Following the Introduction, in Section II we describe the transmission mechanism between joints and motors of the ankle. In Section III we present the impedance controller extended for the humanoid COMAN+. In Section IV we propose the semi-centralized impedance controller. Experimental results and relative discussion are presented in section V. Finally, conclusion and future works are discussed in section VI.

II. JOINT-TO-MOTOR TRANSMISSION

In this section, we describe the parallel ankle mechanism and the kinematics of the joint-motor transmission. A parallel ankle offers numerous advantages over its serial counterpart, though at the expense of more complex kinematics. By using a four-bar mechanism the motors of the ankle can be placed on the shin, decreasing the inertia of the tibia, improving the dynamic behaviour and reducing power consumption. Additionally, a parallel configuration allows the distribution of the torque provided by the actuators to the two degrees-of-freedom of each ankle: motors moving in common mode produce a pitch motion, while differential mode movement results in roll motion (see Figure 1). The transmission between the parallel ankle joint and actuators requires a transformation when passing from the motors variables to joint variables. In particular, joint and motor angles are related through a forward kinematic mapping, and differentiation yields the velocity mapping:

q = J(θ) θ. (1)
where J(θ) ∈ R 2×2 is the Jacobian for the parallel ankle mechanism of COMAN+, q ∈ R 2 are the joint-side velocities and θ ∈ R 2 the motor-side velocities. A complete computation of the Jacobian for this mechanism can be found in [START_REF] Zhou | On the Comprehensive Kinematics Analysis of a Humanoid Parallel Ankle Mechanism[END_REF]. From the principle of virtual work, we obtain a linear mapping between the motor and joint torques, which is required for the design of the ankle impedance controller:

τ m = J(θ) T τ j , (2)
where τ m ∈ R 2 is the torque sensed by the actuator and τ j ∈ R 2 is the torque at the joint (from now on we are going to drop the dependency on θ).

III. IMPEDANCE CONTROLLER

In this section, we describe the implementation of an impedance controller on the parallel kinematics ankles of the humanoid robot COMAN+. The impedance controller at the motor-side has the following formulation:

τ m,d = K m ∆θ + D m ∆ θ + τ m,ff , (3)
where K m ∈ R 2×2 and D m ∈ R 2×2 are respectively the motor-side stiffness and the damping matrices, ∆θ = θ d -θ ∈ R 2 is the error between the desired and measured motor positions, ∆ θ = θdθ ∈ R 2 is the error between the desired and measured motor velocities and τ m,ff ∈ R 2 are feed-forward reference torques at the motor-side. The same law can be defined at the ankle joint level:

τ j,d = K j ∆q + D j ∆ q + τ j,ff , (4)
where the quantities are defined w.r.t. the joint-side instead.

Our goal is to define a motor-side stiffness and damping matrix corresponding to a desired joint-side stiffness matrix:

{K m , D m } = T {K j , D j }, (5)
Since the ankle has two DoFs, pitch and roll, we design diagonal stiffness and damping matrices to decouple the effect of the impedance controller on the two axes:

K j = k p 0 0 k r , (6)
By using ( 2) and (3) we can write:

K m = ∂τ m ∂θ = ∂(J T K j ∆q) ∂θ = J T K j J + ¨¨¨¨B 0 ∂J T ∂θ K j ∆q. (7)
It can be seen how this formulation is position dependent, since the mapping T changes with the robot configuration. If we compute the stiffness matrix around the equilibrium point (∆θ = 0) the second term ∂J T ∂θ K j ∆q can be discarded, though neglecting the conservativeness of the mapping: an in-depth analysis including this term is found in [START_REF] Albu-Schäffer | Cartesian impedance control techniques for torque controlled light-weight robots[END_REF]. A similar approach can be used to obtain the damping D j :

D m = ∂τ m ∂ θ = ∂(J T D j ∆ q) ∂ θ = J T D j ∂ q ∂ θ = J T D j J, (8)
K j = k p 0 0 k r , D j = d p 0 0 d r , K m = k a k c k c k b , D m = d a d c d c d b . (9)
In general, from a diagonal joint stiffness matrix K j , a nondiagonal motor stiffness matrix K m is derived from the mapping. Finally, the term τ m,ff can be obtained from τ j,ff using (2).

It should be noted that ∂(J T Kj ∆q) ∂θ = J T K j J is valid only locally [START_REF] Albu-Schäffer | Soft robotics: what cartesian stiffness can obtain with passively compliant, uncoupled joints?[END_REF], [START_REF] Albu-Schäffer | Cartesian impedance control techniques for torque controlled light-weight robots[END_REF]. Therefore, large displacements from the equilibrium point lead to errors in the stiffness matrix. However, since the parallel ankle mechanism has a narrow range of motion, both on the pitch and the roll axes (i.e. the term ∆q always remains relatively small), the mapping never causes significant errors.

Remark: In this paper, a traditional notation was used to describe the motor-joint mapping q = f (θ) for the sake of clarity. However, in the robot we used the inverse mapping, since the kinematics were provided in this form by [START_REF] Zhou | On the Comprehensive Kinematics Analysis of a Humanoid Parallel Ankle Mechanism[END_REF]. Following the calculations, we require the inverse of the Jacobian to obtain the motor torque: τ m = J(θ) -T τ j . Thanks to the mechanical design of the ankle mechanism the matrix is square and has always full rank: therefore we don't run up against singularity issues.

IV. SEMI-CENTRALIZED IMPEDANCE CONTROL (SCIC)

In this section, we propose an implementation of an impedance controller to make the ankles less susceptible to hardware breakdown or disconnections, reducing potential movements that could compromise the hardware or create unsafe conditions in general. To achieve the desired behaviour in the ankle joint, the controller can be implemented in two different ways for a torque controlled robot.

Pure Torque Control. A conventional way is sending the computed torque directly to the torque controller implemented in each DSP. The on-board processor computes joint torques and transforms them into motor torques before sending them to the DSP of the actuator:

τ m,d = J T (K j ∆q + D j ∆ q + τ j,ff ) . ( 10
)
This method, although viable, is not robust against hardware failures that cause interruptions of the torque reference signals sent the joint local torque controllers. If the connection between the computer and the DSP is severed, the robot behaviour can become unstable, generating high speed motions that could damage the robot and the environment: having solely a torque reference and no position feedback, this impedance control implementation would directly send the last computed torque reference to the DSP causing the ankle to jolt against the mechanical limits with the risk of damaging the actuator transmission system or other mechanical components. Besides, such behaviour can severely compromise safety in any robotics application where humans are in close vicinity of the robot. Since most of the humanrobot collaboration applications rely on the use of impedance control to generate compliant interactions, bounding the robot reactions in the case of system failures is critically necessary.

Semi-Centralized Impedance Control. A way to grant a safer response is sending to the DSP a position reference and the computed stiffness and damping diagonals, and let the DSP calculate the torque using the feedback of the robot. This allows to bound the system motion reaction, resulting in a graceful and safer degradation if a system failure causing discontinuities of the torque reference signals occurs.

To achieve this scheme, we can arrange the matrices K m and D m , derived from the joint-motor transformation, in the following way:

K m = K m + K m , D m = D m + D m , (11)
where (•) m contains the diagonal elements and (•) m contains the off-diagonal elements of the matrix. Recalling equation (3), we can write:

τ m,d = K m ∆θ + D m ∆ θ + τ m,ff , (12)
where:

τ m,ff = K m ∆θ + D m ∆ θ + J T τ j,ff . (13)
This formulation enables the decentralization of the input signals. The first two terms in [START_REF] Hyeon | Hybrid control for biped robots using impedance control and computed-torque control[END_REF] are computed in each DSP, each row separately. Conversely, the τ m,ff component in ( 13) is computed on the on-board computer in a centralized way: since the i th DSP has information solely about the i th actuator to which is connected, it can only determine the torques dependent on position θ i and velocity θi , whereas both terms (θ i and θi , ∀i = 1, 2) are needed to obtain each motor total torque. Therefore, the on-board computer uses the off-diagonal terms K m and D m to compute the part of the torque that depends on the other joint, which is then sent directly to the DSP as a feed-forward component. Term τ j,ff is a torque directly commanded to the joint by the user. Thus, in case of a failure:

• the last position sent to the robot is tracked by the DSP which continuously receive the position feedback locally; • the last feed-forward torques received, provided by the remaining components of the stiffness and damping matrices, will act as a disturbance in the resultant position feedback loop implemented using the diagonal stiffness and damping elements. However, this disturbance is always bounded and not large enough to induce instabilities. Our controller, in case of system faults that could interrupt the references signals to the joint local controllers, is then reduced to a simple closed loop with a bounded disturbance at the input which can be rejected. This favors the hardware protection and results in slower and safer robot reactions, especially important in human-robot interaction. Remark: A similar semi-centralized method was implemented in [START_REF] Albu-Schäffer | Cartesian impedance control techniques for torque controlled light-weight robots[END_REF] under different circumstances, i.e. a Cartesian impedance controller used to implement a compliant behaviour for robotic manipulators. However, the rationale behind our implementation differs from the one in [START_REF] Albu-Schäffer | Cartesian impedance control techniques for torque controlled light-weight robots[END_REF], as we are exploiting it to demonstrate its benefits in ensuring a bounded robot reaction in case of system failures.

V. EXPERIMENTS AND DISCUSSION

The robot platform used in this work is the lower body of the humanoid robot COMAN+, built at IIT, composed of 13 series elastic actuated (SEA) joints, providing a first layer of passive compliance and the ability to sense external forces. The software architecture of COMAN+ is powered by XBotCore [START_REF] Muratore | XBotCore: A Real-Time Cross-Robot Software Platform[END_REF], a software framework which permits to run hard real-time control loops at 1 kHz. The ankles motors were controlled in impedance mode: the control loop in each DSP consists of an impedance controller which takes as inputs the desired torque together with desired values of stiffness and damping. The impedance controller sends the computed reference torque to the lower-level torque controller of the actuator.

To validate the SCIC we compared it to the conventional pure torque control (PTC). In particular, we designed two classes of experiments to estimate performance, detailed in V-A and V-B: the first one is a set of experiments to measure the accuracy of the SCIC in rendering a desired impedance, and the second is a qualitative comparison to prove its efficacy in reducing the system motion reaction when a failure occurs. Finally, an experiment was designed to show qualitatively the advantages of an impedance controller implemented on the ankle joint. While lifted by the crane, the robot was set in a fixed standing configuration. Subsequently, it was lowered on an inclined wooden board. The experiment was repeated using three different strategies on the ankle joint:

• impedance control mode;

• position control mode;

• no control. Videos of the experimental comparisons can be found at https://youtu.be/-QULP4xUD9A.

A. Accuracy validation

For the sake of clarity, we conducted two independent set of experiment, one for the stiffness alone and one for the damping, so that we could decouple the effects of damping and stiffness in (4). The feed-forward commanded torque was set to zero. Separate experiments were carried out for both ankle DoFs. The robot was lifted with a crane in order to raise it above the ground, starting with a nominal commanded stiffness (or damping). Then, an external force was applied along one DoF of the left ankle, while the commanded stiffness (or damping) was changed by the user. In particular, the ankle was kept in place while changing the stiffness reference, since its value depends on the position error. On the contrary, the damping depends on the velocity error; hence the ankle was moved while changing the damping reference to estimate its value. The sensed motor torque and position were translated into the sensed joint torque and position by using the Jacobian associated with the joint-motor mapping. A simple way to obtain the measured stiffness (or damping) is rearranging (4), but the results are way too noisy. Hence, we used a Recursive Least Square Filter (RLSF) to estimate the desired values. This adaptive algorithm allows to recursively evaluate stiffness and damping over time by minimizing a weighted least square cost function. We conducted two experiments for each degree of freedom, one for the stiffness and one for the damping.

B. Efficacy validation

To validate our controller compared to the traditional one in the case of communication interruptions, we used the same experimental setting of the previous experiments, keeping the robot above the ground. While the SCIC was controlling the ankle with a non-zero stiffness and damping, a force was applied along one axis. In this configuration, the EtherCAT cable connecting the on-board computer to the actuators was unplugged, simulating a hardware failure. The same experiment was conducted with a pure torque input controller.

C. Results

The first class of experiments validated the accuracy of our method compared to the traditional torque input approach. Figures 2 and3 shows the reference stiffness signal and the measured stiffness using the semi-centralized impedance controller (SCIC) and a pure torque input (PTC) for the roll and pitch DoF, respectively. Figures 4 and5 shows the same graphs for the measured damping. All the measured values are filtered with the Recursive Least Square (RLS) filter. Note that the differences between the two reference values in each graph are due to the fact that the data for the SCIC and for the PTC belongs to separate experiments and the commands sent to the robot by the user differ in time. The measures from the SCIC are close to the reference signal and comparable to the ones from the PTC, thus validating our method.

The second experiment demonstrates the effectiveness of the SCIC's design: as the disconnection between on-board computer and DSPs occurs, our approach allows the ankle's DSP to keep the last position sent by the SCIC, while using the traditional control the last torque sent before the disconnection is kept as reference input, resulting in a rapid motion of the ankle joint towards its limit. Figure 7 shows a sequence of frames comparing the behaviour of the robot controlled with the SCIC and the PTC. Due to the nature of this particular experiment, i.e. unplugging the etherCAT cable to simulate a hardware failure, we couldn't collect any data.

Finally, the qualitative experiment proved the importance of the impedance control in granting the robot ankles with the ability to conform to the terrain. The ankles adapt to the slope without yielding under the weight of the robot. This desired behaviour falls between two unwanted extremes, shown in the remaining trials of Figure 8: the stiff position control doesn't allow the ankle to comply with the terrain, which rigidly keeps its position when the robot is lowered on the inclined board, while a non-controlled ankle collapses under the weight of the robot. Note that in these experiments we don't implement any balancing strategy on the robot, which is standing entirely thanks to its initial configuration.

A video including the full sequence of the controllers comparison (SCIC-PTC) and the qualitative experiments is available at https://youtu.be/-QULP4xUD9A.

D. Offset discussion

In this section, we discuss the offset between the reference signal and the measured values in the graphs of Figures 2,3, 5 and 4. The offset error, which is particularly evident in Figure 4, arises from an offset in the torque tracking of the motor lower level control, which is shown in Figure 6. The corresponding stiffness and damping offset can be computed by solving the equation (4) for the stiffness or the damping. Note that the equation always has one addend only, since in the experiments either pure stiffness or damping was commanded to the robot and the feed-forward component was set to zero. The following expression is obtained for the stiffness offset: where δK meas j and δτ meas j are the offset error in the measured stiffness and torque, respectively. Since we measure the torque on the motor level, we use equation ( 2) to map the sensed torque to the joint level and compute δτ meas j . A similar equation can be written to estimate the offset in the damping:

δK meas j = diag (∆q) -1 J -T • δτ meas j ( 14
)
δD meas j = diag (∆ q) -1 J -T • δτ meas j ( 15
)
Both experiments to estimate stiffness and damping have an offset, but in the stiffness experiments the error is significantly smaller. This is due to the different nature of the two experiments: since the quality of the torque tracking is worse during torque reference variations, the experiment to estimate the damping, in which the ankle is moved, is affected by a higher torque tracking error compared to the experiment to estimate the stiffness, in which the ankle stays still. In the first graph of Figure 6 the two motors actuating the ankle joint are moving in common mode, and their offsets are concordant. In the second graph the two motors are moving in differential mode, and their offsets are discordant. This justifies both the offset in Figure 4, where the offsets from the two motor readings sum up, as well as the absence of an error in Figure 5, where the two errors cancel each other out. Finally, Figure 4 shows the damping measures corrected by compensating the offset: the resulting values follow closely the reference signal.

VI. CONCLUSION

The semi-centralized impedance controller was successfully implemented on the COMAN+ robot, extending it for the parallel kinematics mechanism of the ankles. We experimentally proved the effectiveness of our method, which performs closely to a conventional impedance control implementation based on pure local torque controllers: no (a) The SCIC allows a safe response of the robot when the external force exerted on the robot is removed, the robot tracks the last position sent by the controller, returning to the nominal configuration.

(b) The PTC makes the ankle jolt, reaching violently the joint limit, when the external force exerted on the robot is removed. This sudden reaction could push the robot off balance or cause damages to the ankle joint. evident difference is present between the experimental data of the SCIC and the PTC, as shown in Figures 2345. The offset error in the measurements, due to the torque tracking error in the lower control, was analyzed and effectively compensated. Furthermore, we showed how the implementation of the impedance controller improves the behaviour of the COMAN+ lower body, allowing the ankles to comply with the environment.

The main contribution of this paper is the design, implementation and testing of an impedance controller capable of withstanding potential system failures without generating rapid motion reactions, which could compromise the hardware integrity or the safety of humans interacting with the robot. This is a fundamental prerequisite for a robotic platform used in real-time operations such as human-robot interactions or rescue missions, where the effects of failure should be reduced as much as possible to protect both the hardware and the human operators. Our implementation allows the robot to fail gracefully, as opposed to a traditional impedance control implementation relying on pure local torque controllers, in which interruptions or breakdowns of the main on-board computer could lead to abrupt movements of the joints and dangerous behaviours. The SCIC was experimentally validated in the presence of disconnections: by simulating a hardware failure, we proved how the ankle doesn't drift violently, but returns to its nominal position. This paper follows the development of the new robot CO-MAN+: future work will build upon this work to further develop and apply the proposed scheme to the full humanoid robot.

Fig. 1 .

 1 Fig. 1. The legs of the humanoid robot COMAN+ on the left and close back view of its tibia (CAD) without covers on the right. When the ankle's motors situated on the shin move in common mode (both rotate in the same direction, as illustrated by the orange arrows), a pitch motion is produced. Conversely, the motors moving in differential mode result in a roll motion of the ankle (green arrows).

Fig. 2 .

 2 Fig. 2. Reference stiffness for the roll axis DoF commanded manually by the user (steps of 100, 50, 20 Nm/rad) compared to the measured ankle stiffness with the semi-centralized impedance control (SCIC) and the pure torque control (PTC).

Fig. 3 .

 3 Fig. 3. Reference stiffness for the pitch axis DoF commanded manually by the user (steps of 100, 50, 20 Nm/rad) compared to the measured ankle stiffness with the semi-centralized impedance control (SCIC) and the pure torque control (PTC).

Fig. 4 .

 4 Fig. 4. Reference damping for the pitch axis DoF commanded manually by the user (one step from 0 to 50 Nm s/rad) compared to the measured ankle damping with the semi-centralized impedance control (SCIC) and the pure torque control (PTC). The offset between measured and reference values is due to a torque tracking error. The compensated measures are shown, closely tracking the reference values of damping.

Fig. 5 .

 5 Fig. 5. Reference damping for the roll axis DoF commanded manually by the user (one step from 0 to 50 Nm s/rad) compared to the measured ankle damping with the semi-centralized impedance control (SCIC) and the pure torque control (PTC).

Fig. 6 .

 6 Fig. 6. Torque offset caused by a torque tracking error of the motor lower level control, filtered with a zero-phase digital filter. Above: the ankle motors are moving in common mode, resulting in a rotation of the ankle along the pitch axis. The torque offsets have the same sign. Below: the ankle motors are moving in differential mode, resulting in a rotation of the ankle along the roll axis. The torque offsets have opposite sign.

Fig. 7 .

 7 Fig. 7. Sequence of frames of the robot COMAN+'s behaviour when disconnected from the on-board computer. The red circles highlight the etherCAT, which is unplugged in the second frame simulating a hardware failure.

Fig. 8 .

 8 Fig. 8. Frames of the COMAN+ robot being lowered on a wooden board inclined of 16 • . In a) the robot is position-controlled. The ankles don't adapt to the slope, making the robot fall backward. In b) the robot is controlled with an impedance controller (stiffness: 30 Nm/rad, damping: 5 Nm s/rad), which allows the ankles to comply with the inclined board. Finally, in c) the robot is not controlled. The ankles collapse under the weight of the robot, making it fall forward.

ACKNOWLEDGEMENTS

The research leading to these results has received funding from the European Union Horizon 2020 robotics program [ICT-23-2014], grant agreements n.644727 (CogIMon).


A Ghzayel
email: alaa.ghzayel@univ-poitiers.fr

A Beaudoin
email: anthony.beaudoin@univ-poitiers.fr

Scour downstream of a sluice gate using SedFoam and FLOW-3D

Keywords: Sediment transport, fluid mechanics, numerical simulation. 1

 will be used as the reference to study this phenomenon. This will show the advantages and defaults of both models as well as reproduce numerically this phenomenon.

Introduction

Sediment bed erosion, near hydraulic structures, that is caused by the erosion due to water flow that interact with part of this structures is a natural phenomena called "scour". Scour has been an arising problem that can cause severe failure in fluvial and marine hydraulic structures. Calculating and forecasting the maximum depth that the scour can reach and the size and shape of this scour is interesting to predict points of failures and to protect against such failures in these hydraulic structures. Sediment scour due to a water wall jet downstream of a sluice gate in a channel with a high degree of confinement is not well documented, several experiments were made but to our knowledge, no numerical simulations were tested to describe the hydro-morphological phenomena of sediment erosion and deposition downstream of a sluice gate in a very confined channel. In 2019, Martino et al. tested for several flow rates, the sediment scour formation caused by a water jet downstream of a sluice gate in a very confined channel [START_REF] Martino | Experimental study on the scour due to a water jet subjected to lateral confinement[END_REF], where a sediment erosion and a dune formation were observed. Furthermore, the dune changed its shape from a triangle to a trapezoidal shape and for certain flow rates a phenomenon of cycles of erosion and filling of the scour was also observed, which is not observed in the literature for non-confined channels. The model was reproduced numerically using two different models, SedFOAM [START_REF] Chauchat | Sedfoam-2.0: a 3-D two phase flow numerical model for sediment transport[END_REF] and FLOW-3D [START_REF]FLOW-3D Version 12.0 User's Manual[END_REF], where the experimental results were compared to the numerical results and an explanation of the observed phenomena was deduced from the numerical results as well as the effect of the water jet on the erosion process and dune shape formation.

Numerical models

SedFOAM, is a two-phase flow RANS model for sediment transport based on the open-source toolbox, OpenFOAM, in which dynamical equations are applied on both the fluid and particles phases. This approach takes into consideration most of the physical processes. Two important closures are implemented, the flow turbulence and granular stress closures. For the turbulence closure there are two models, the k -ε model [START_REF] Bakhtyar | Two-phase hydrodynamic and sediment transport modeling of wave-generated sheet flow[END_REF] and the most recent k -ω model [START_REF] Amoudry | Extension of turbulence closure to twophase sediment transport modelling: Application to oscillatory sheet flows[END_REF] similar to that of Amoudary (2014), the k -ω model was used. While for the granular stress models, the dense granular flow rheology initially proposed by GDRmidi (2004) [START_REF] Gdrmidi | On dense granular flows[END_REF] is used. On the other hand, FLOW-3D is a commercial CFD software that use a RANS-VOF method, which uses a three-dimensional sediment transport and scour model that is not considered as a coupled multi-phase flow model. In this model, the Flow is calculated using the Navier-Stokes equations with turbulence closure but without taken the solid phase into consideration. The RNG k -ε model [START_REF] Yakhot | The renormalization group, the e-expansion and derivation of turbulence models[END_REF] is used and calibrated by the software to assure numerical stability. As for the sediment transport, a conventional equation of sediment transport like the Nielsen equation [START_REF] Nielsen | Coastal bottom boundary layers and sediment transport[END_REF] as well as the Meyer-Peter and Muller equation [START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF]

Results

Numerical simulations were performed using several flow rates Q, for a sluice gate aperture of 1 cm. The numerical results were then compared with the experimental data from the study of [START_REF] Martino | Experimental study on the scour due to a water jet subjected to lateral confinement[END_REF]. A comparison of the form of the dune and change of shape of the dune from triangle to trapezoidal shape was compared at first between the experimental data and the numerical results on both models, SedFoam and FLOW-3D (Figure 2). The three-dimensional numerical simulations were able to reproduce the form of the dune as well as the shape of both the scour and dune since the dune submitted a change in form, from a triangular shape to a trapezoidal shape on both numerical models. These numerical results showed the capability of these models to perform scour simulations with a high degree of physical accuracy.

Figure 2: Dune formation, Experimental (A), SedFoam (B) and FLOW-3D (C).

This shape shifting was only observed in a very confined channel. The numerical results were only able to reproduce this phenomenon using three-dimensional numerical simulations which indicate the importance of walls on giving accurate numerical results as well as their impact on the dune size and shape. More numerical results are to be analyzed during this research where the scour maximum depth of both models will be compared with the experimental data. The evolution of the maximum scour depth through time will also be studied and analyzed based on the experimental data to check for the accuracy of these models in reproducing every aspect of the confined scour phenomena. Moreover, the numerical results of the water velocity will be validated through the experimental data to assure a good agreement of the velocity profiles between the experimental and the numerical results. This will also allow us to check if the numerical models are capable of reproducing the velocity decay phenomena through a confined channel. This phenomenon occurs for wall jets where the maximum velocity of the jet decreases the more it gets far from the sluice gate aperture in the stream wise direction.

Conclusion

Two different numerical models, SedFoam and FLOW-3D Hydro, were used to study a particular phenomenon of scouring downstream of a sluice gate in a confined channel. The main objective of this

(A) (B) (C)
study is to reproduce the scour phenomena observed in the experimental research done by [START_REF] Martino | Experimental study on the scour due to a water jet subjected to lateral confinement[END_REF] which will show the capability of both models as well their advantages and defaults in such complex conditions. Moreover, the numerical results will be validated using the experimental data, especially the dune shape shifting as well as the maximum scour depth. This validation can serve as a numerical data base that can help complete the experimental data, where measurements were not made as well as a tool to explain the scouring phenomena observed in the experimental procedure which will help later on in predicting maximum scour depth and changes in the scour and dune characteristics.

  is implemented and used to calculate the transport of the sediments. The numerical results were validated using the experimental data by Martino et al. (2019).3 Experimental setup A hydro-sedimentary canal setup used by Martino et al. (2019) is reproduced numerically. It consists of a 140 cm long confined channel of a 2 cm width and lateral glass walls confinements. A pump is installed that can control the flow rate Q. Several flow rates were used to test the scour formation as a function of the flow rate. A 60 cm long sediment bed composed of 1.1 mm diameter granular glass beads and a density of 2.65 g/cm 3 was laid directly downstream of the sluice gate with a 10 cm depth. A 13 cm water depth is used. The used aperture of the sluice gate is 1 cm.

Figure 1 :

 1 Figure 1: Experimental setup sketch by Martino et al. (2019).


Matteo Parigi Polverini

Enrico Mingo Hoffman

Arturo Laurenzi
email: arturo.laurenzi@iit.it

Nikos G Tsagarakis
email: nikolaos.tsagarakis@iit.it

Sparse Optimization of Contact Forces for Balancing Control of Multi-Legged Humanoids

Keywords: Optimization and Optimal Control, Humanoid Robots, Force Control

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

T HE control of multi-legged humanoid robots is a re- cent research topic that combines balancing, flexible mobility and precise manipulation. Only few examples of these advanced research platforms are currently available, e.g. the Momaro robot [START_REF] Schwarz | Nimbro rescue: Solving disaster-response tasks with the mobile manipulation robot momaro[END_REF], RoboMantis by Motiv Robotics, and CENTAURO [START_REF] Baccelliere | Development of a human size and strength compliant bi-manual platform for realistic heavy manipulation tasks[END_REF] a hybrid wheeled-legged quadruped equipped with a bimanual upper-body developed at Istituto Italiano di Tecnologia (IIT). A distinctive feature of this type of platforms is their redundancy with respect to the number of end-effectors required to perform balancing and/or manipulation tasks. We hereafter refer to this feature as end-effector redundancy, in contrast to the common kinematic redundancy. Consider for example the balancing problem for a bipedal humanoid. As a consequence of end-effector redundancy for this task, a stable balance can be accomplished not only with two feet on the ground, but also with a single foot on the ground, or even combining foot and hand contacts. This applies even Fig. 1. The left sketch shows a possible use of multiple end-effectors and retargeting of the foot for a pushing task. The right sketch shows an example of singularity avoidance for the leg: one leg is lifted, instead of being stretched, during a balancing task. more to a multi-legged humanoid as CENTAURO, which can establish up to six independent contacts with the environment. Indeed, end-effectors redundancy opens up several possibilities of application, see Fig. 1, spanning:

end-effector re-targeting, e.g. think of a bipedal humanoid standing on one foot while the other foot is re-targeted to help the arms pushing an heavy object; -use of multiple end-effectors, e.g. push an object with one hand and additionally with another hand or a foot, if actuation constraints are reached; -singularity avoidance, see the right sketch in Fig. 1.

Although these possibilities seem attractive, in order to be able to exploit end-effector redundancy, it is crucial to find the smallest set of end-effectors required to perform a given primary task. For the balancing control of a torque-controlled humanoid, as the CENTAURO robot, this is equivalent to finding a sparse solution of a contact force distribution problem, i.e. a solution in which most of the contact forces are zero, thus minimizing the number of employed end-effectors. As for the authors' knowledge, this is the first work that exploits the concept of sparsity for realtime control purposes, whereas approaches employing sparse optimization with highly redundant robots can be found in the context of inverse kinematics [START_REF] Ho | Computing inverse kinematics with linear programming[END_REF], [START_REF] Gonc ¸alves | Parsimonious kinematic control of highly redundant robots[END_REF], and planning [START_REF] Ponton | A convex model of humanoid momentum dynamics for multi-contact motion generation[END_REF]. Two different methodologies are presented and discussed in this paper to achieve a sparse optimization of contact forces for the balancing control of a multi-legged humanoid robot. The proposed optimization methods are formulated and discussed as a general regularization approach to a Least Squares (LS) estimation problem and consequently extended to the balancing control. An extensive analysis performed on a simplified case study complements the discussion on the impact of the proposed approaches on the degree of sparsity of the solution.

Finally, the effectiveness of the proposed approaches has been validated in simulation using the CENTAURO robot.

The paper is organized as follows. Existing control approaches to humanoid balancing are reported in Sec. II, while background knowledge on sparse regularization methods is provided in Sec. III. The proposed sparse optimization methods are presented in Sec. IV and compared on a simplified case study. Gazebo simulation results performed on the CENTAURO robot are presented in Sec. V. Finally, concluding remarks and future work directions can be found in Sec. VI.

II. RELATED WORKS ON BALANCING CONTROL

State-of-the-art control approaches to the humanoid balancing problem are generally classified in two categories. A first category of methods performs a pre-optimization of contact forces with respect to the robot centroidal dynamics, and subsequently maps the computed contact forces to the actuated joint torques, see [START_REF] Ott | Posture and balance control for biped robots based on contact force optimization[END_REF]- [START_REF] Stephens | Dynamic balance force control for compliant humanoid robots[END_REF]. In contrast, the second category entirely exploits the full-body inverse dynamics, see e.g. [START_REF] Sentis | Synthesis of whole-body behaviors through hierarchical control of behavioral primitives[END_REF], [START_REF] Righetti | Optimal distribution of contact forces with inverse-dynamics control[END_REF] and the methods based on a hierarchical Quadratic Programming (QP) formulation of the inverse dynamics [START_REF] Saab | Dynamic whole-body motion generation under rigid contacts and other unilateral constraints[END_REF], [START_REF] Herzog | Momentum control with hierarchical inverse dynamics on a torque-controlled humanoid[END_REF].

As an alternative approach, see [START_REF] Koolen | Design of a momentum-based control framework and application to the humanoid robot atlas[END_REF], [START_REF] Kim | Computationally-robust and efficient prioritized whole-body controller with contact constraints[END_REF] and the authors' work in [START_REF] Laurenzi | Balancing control through post-optimization of contact forces[END_REF], the balancing problem for torque-controlled legged robots can be addressed through a post-optimization of contact forces. This way the the tricky specification of a desired momentum of rotation inherent in the first category of methods is circumvented, additionally reducing the number of decision variables and thus the computational load compared to the second category.

Regardless of the specific methodology, the aforementioned approaches share the common drawback (from the perspective of the present work) to employ the total number of endeffectors in contact with the environment to perform the balancing task. This is done, see e.g. [START_REF] Ott | Posture and balance control for biped robots based on contact force optimization[END_REF], by computing the Least Squares (LS) solution of the force distribution problem through pseudo-inversion. As it will be shown in the following Section, this solution minimizes the Euclidean norm of the contact forces, although resulting in non-zero contact forces associated to each end-effector, i.e. it is inherently not sparse. In contrast, our main objective is to propose a methodology to automatically find the sparsest solution to the force distribution problem.

III. BACKGROUND ON SPARSE REGULARIZATION

This Section provides background knowledge on sparse regression approaches to Least Squares (LS) estimation.

A. Problem statement

Consider the Least Squares (LS) estimation problem:

min x∈R n Ax -b 2 2
(1) where x ∈ R n is the vector of optimization variables, A ∈ R m×n is the design matrix, b ∈ R m is the target vector. In order to give preference to solutions with smaller norm, the l 2 -regularization of (1) is given by: min

x∈R n Ax -b 2 2 + λ x 2 2 (2)
where • 2 is the Euclidean norm and λ ∈ R + is a weighting parameter. In this respect, the Moore-Penrose inverse (also known as left pseudo-inverse, or simply pseudo-inverse) of the design matrix is probably the most widely adopted method ensuring the minimization of the Euclidean norm of the optimization variables, i.e.: min

x∈R n x 2 2 s.t. Ax = b. (3)
Nevertheless note that l 2 -regularization and pseudo-inversion do not encourage sparsity, and the resulting solution typically has non-zero values associated with all the variables.

B. LASSO regression

The Least Absolute Shrinkage and Selection Operator (LASSO) is a well known linear regression technique obtained by adding a l 1 -regularization term to a LS estimation [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. LASSO regression has the tendency of yielding a sparse solution, in which some optimization variables are equal to zero, and it can be formulated as the following unconstrained optimization problem:

min x∈R n Ax -b 2 2 + λ x 1 (4)
where • 1 is the l 1 -norm. Despite the continuous but nondifferentiable nature of (4), the LASSO can be converted to a tractable constrained Quadratic Programming (QP) problem by introducing additional variables t ∈ R n , e.g. as in [START_REF] Stellato | Osqp: An operator splitting solver for quadratic programs[END_REF]:

min x,t∈R n Ax -b 2 2 + λ n i=1 t i s.t. -t ≤ x ≤ t (5)
It can be finally proved that any optimal solution of (5) satisfies:

t * i = |x * i | ∀i = 1, . . . , n
hence showing the equivalence of (4) and [START_REF] Ponton | A convex model of humanoid momentum dynamics for multi-contact motion generation[END_REF].

C. MILP regression

Mixed Integer Linear Programming (MILP) is a Linear Programming (LP) class of problems in which a subset of variables can assume only integer values. A generic MILP problem can be expressed as:

min x,δ c T x + d T δ s.t. Ax + Dδ ≤ b, (6)
where x ∈ R n1 are continuous optimization variables, δ ∈ Z n2

+ are variables having integer non-negative components, c ∈ R n1 , d ∈ R n2 are the linear cost functions, A ∈ R m×n1 , D ∈ R m×n2 and b ∈ R m are the problem constraints.
In order to exploit MILP to perform sparse regression, it is necessary to introduce a set of n binary variables δ associated to each variable of x enforcing the following relation:

δ i = 0, x i = 0 1, x i = 0 , ∀i = 1, . . . , n (7)
This way δ represents an indicator of which one of the variables x is used, i.e. it is not zero. By applying the socalled big-M method, (7) can be formulated as the following big-M or indicator constraint:

-M δ i ≤ x i ≤ M δ i , ∀i = 1, . . . , n (8)
where M ∈ R + is a big enough value to guarantee that the constraint is not active when δ i = 1. Recalling that finding a sparse solution to the LS regression problem in (1) means to maximize the number of elements of x which are equal to zero, the resulting MILP optimization problem can be formulated as follows:

min x,δ n i=1 δ i s.t. Ax = b, -M δ ≤ x ≤ M δ, δ ∈ {0, 1} n (9)
Compared to LASSO, MILP regression ensures the maximization of the sparsity of the solution. Although hybrid optimization models are more difficult to solve than continuous models, there exist methods e.g. Branch-and-Bound algorithms that, by combining a partial enumeration strategy with relaxation techniques, can obtain the global optimum with much fewer sub-problems solved than with the complete enumeration method.

D. Comparison on a simplified example

To better understand the differences between these approaches, we here consider the following simplified LS estimation problem:

min x∈R 2 a T x -b 2 2 (10) with a = 2 T , b = 1
This analysis has been performed in MATLAB. The l 2regularization has been performed according to (3) using pinv function, while LASSO regression has been performed using linprog to solve the following equivalent LP constrained problem:

min x,t∈R 2 x 1 s.t. a T x = b (11)
Finally, the MILP problem in (9) has been solved using intlinprog. The solutions of the different optimization problems are geometrically illustrated in Fig. 2(a).

As it can be noticed, the solution of l 2 -regularization, i.e.

x * = 0.2 0.4 T , minimizes the Euclidean norm of all the variables, at the expense of resulting in non-zero values. On the contrary, both LASSO regression and MILP produce a sparse solution where one variable is equal to zero. Nevertheless, note that while LASSO solution is unique, two optimal solutions exist for the MILP problem. As a final comment, notice that, in case box constraints are further introduced so that the LASSO solution of ( 11) is ruled out, e.g.:

-0.25 ≤ x 2 ≤ 0.25 LASSO regression is unable to provide a sparse solution, whereas MILP does, see Fig. 2(b). As a matter of fact, this example highlights how, unlike the LASSO, MILP regression always returns the sparsest solution compatible with the problem constraints.

IV. BALANCING CONTROL OF A LEGGED ROBOT

The focus of this work is the application of ideas from the sparse regression field to the balancing control of endeffector-redundant legged robots. In this section, we introduce the main concepts about balancing from two recent works by the authors [START_REF] Laurenzi | Balancing control through post-optimization of contact forces[END_REF], [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF] that are later used in conjunction with sparse optimization techniques.

A. Floating-base dynamics

Fixed-base and floating-base robots differ in terms of two major aspects, namely under-actuation and, as a consequence, the prominent role played by contact forces for their control. Indeed, legged robots must always interact with the environment in order to be controlled in their full (6+n)-dimensional coordinate space. The dynamics of a floating-base robot is explained by the following equation:

B(q)q + h(q, q) = Sτ + J T C F C + J T F x , (12)
which highlights two major differences w.r.t. the fixed-base case. First, an actuation matrix S ∈ R (n+6)×n models the system under-actuation. Second, task dynamics becomes coupled with contact forces, which are taken into account by considering the Jacobian of all support links J C ∈ R k×(n+6) and the corresponding overall contact wrench F C ∈ R k , with k equal to the contact constraint dimension (e.g. k = 12 for a humanoid in double support).

In [START_REF] Laurenzi | Balancing control through post-optimization of contact forces[END_REF], the authors have shown that it is possible to deal with such a coupling by separating the inverse-dynamics problem into two independent stages. The first stage consists in a fixed-base hierarchical force control problem, as introduced by the authors in [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF]; subsequently in the second stage, the system under-actuation is tackled by a post-optimization of the contact forces. While some details about this procedure are given in the following Section, the reader can refer to [START_REF] Laurenzi | Balancing control through post-optimization of contact forces[END_REF] for a detailed description.

B. Post-optimization of contact forces

First and foremost, an equivalent fully-actuated torque vector τ is defined as follows:

τ = Sτ + J T C F C , (13)
such that the floating-base dynamics formally resembles the fixed-base one:

B(q)q + h(q, q) = τ + J T F x . (14)
Consequently, our algorithm in [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF] can be applied without any modification, yielding some optimized value τ * for the fullyactuated torque vector (13) that permits to achieve the desired hierarchical motion and interaction. Then, the contact force information is recovered by solving the system of equations ( 13) for τ and F C . Furthermore, the structure of the actuation matrix S can by exploited in order to improve efficiency by reducing the number of unknowns. Indeed, considering the topmost six rows, the following formula is obtained:

τu = J T C,u F C , (15)
where the subscript "u" indicates the sub-matrix corresponding to the unactuated virtual joints. Equation ( 15) is torqueindependent, and contains only the contact forces as variables.

Once that these have been determined, joint torques can be recovered by looking at the bottom n rows of (13) and solving for τ . The transpose of the unactuated part of the contact Jacobian acts as a grasp matrix G ∈ R 6×k :

G = J T C,u (16)
hence, a force distribution problem can be formulated according to several possible criteria. Assuming n C point contacts, i.e. no torque is exchanged, the following inequalityconstrained LS optimization has been explored in [START_REF] Laurenzi | Balancing control through post-optimization of contact forces[END_REF]:

min FC GF C -τu 2 2 s.t. b l ≤ D F C ≤ b u u l ≤ F C ≤ u u . (17)
whereas in the present work we aim to find a sparse solution.

Pyramidal linearized friction cones have been considered as inequality constraints in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], i.e.:

F C,n ≥ 0, |F C,t | ≤ √ 2 2 µ F C,n (18)
where µ ∈ R is the friction coefficient. Note that the contact forces, and the corresponding Jacobians, are expressed with respect to a local reference frame such that the z-axis represents the normal component of the contact force.

C. Block-symmetrical structure of the grasp matrix

We now discuss the structure of the grasp matrix G in ( 16), again assuming n c point contacts. To this aim, recall that the un-actuated part of the contact Jacobian J C,u encodes the relationship between the linear velocity of the contact points and the virtual joint velocities qu ∈ R 6 . Considering the i-th contact point velocity v C,i , the following equation expresses the relation between the linear and angular velocity of the robot floating-base, v fb ∈ R 3 and ω fb ∈ R 3 , respectively, and v C,i :

v C,i = v fb + ω fb × r C,i , (19)
where the vector r C,i ∈ R 3 points from the floatingbase origin to the i-th contact point. Let us further denote T (q u ) ∈ R 6×6 the floating-base Jacobian matrix realizing the linear mapping:

v fb ω fb = T (q u ) qu (20)
By now combining ( 19), ( 20) and ( 16), the resulting expression for the grasp matrix is given by:

G = T (q u ) T I 3 • • • I 3 S(r C,1 ) • • • S(r C,nC ) , ( 21
)
where S is a skew-symmetric matrix such that S(a) b = a × b.

At this point, it is important to highlight the blocksymmetric nature of G. In fact, if the contact points r C,i are located symmetrically w.r.t. one or more axis of the floatingbase local reference frame (a particularly recurrent situation on a multi-legged humanoid, like the quadruped CENTAURO), this eventually reflects on a block-symmetrical pattern of the last three rows of the grasp matrix, in addition to the inherently block-symmetric first three rows. As it will be shown in the following Section, this peculiar feature of G, i.e. the design matrix for the LS force distribution problem in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], directly affects the degree of sparsity of the solution.

D. Impact of symmetry of the design matrix on sparsity

In order to understand the impact of a symmetrical structure of the design matrix on the solution produced by sparse regression approaches, we hereafter resume the simplified example previously introduced in Sec. III-D. By now choosing:

a = 1 1 T , b = 1
the simplest symmetrical structure of the design vector a is considered; note that both optimization variables now contribute equally to the target value. The l 2 -regularization has been performed according to (3) using pinv function, while LASSO regression has been performed using linprog to solve the LP constrained problem in [START_REF] Righetti | Optimal distribution of contact forces with inverse-dynamics control[END_REF]. Finally, the MILP problem in (9) has been solved using intlinprog. The solutions of the different optimization problems are geometrically illustrated in Fig. 3. Once again, the solution of the l 2 -regularization, i.e.

x * = 0.5 0.5 T , tends to minimize the Euclidean norm of all the variables at the expense of sparsity, while two equally Fig. 3. Graphical interpretation for a = 1 1 T . Due to the parallelism between the line representing the equality constraint a T x = b and the contours of the l 1 -norm cost function in [START_REF] Righetti | Optimal distribution of contact forces with inverse-dynamics control[END_REF], an infinite number of optimal solution exist for the LASSO approach, of which only two are sparse.

optimal sparse solutions exist for the MILP regression. Particularly interesting in this case, is the behaviour of the LASSO, which now has an infinite number of optimal solutions, only two of which are sparse. As it can be noticed from Fig. 3, this is due to the fact that the line representing the equality constraint a T x = b is now parallel to the contours of the l 1 -norm cost function.

Following from this example, it can be inferred that, due to the block-symmetrical nature of G, multiple sparse solution to the force distribution problem can arise using LASSO and MILP regression. As a matter of fact, considering for instance the balancing problem for a biped, this is perfectly coherent with the fact that a stable balancing can be equivalently achieved by standing on the left or on the right foot. Although acceptable, this could result in high-frequency switches between supporting legs, thus affecting the applicability of the proposed approach for control purposes.

In order to improve this behaviour aiming at ensuring the uniqueness of a sparse solution, we propose to introduce a weighted penalty formulation. For the considered simplified example, being γ ∈ R 2 a weighting vector, the resulting weighted LASSO optimization problem can be expressed as:

min x,t∈R 2 γ T t s.t. a T x = b, -t ≤ x ≤ t ( 22
)
while the weighted MILP regression is given by: min

x,δ γ T δ s.t. a T x = b, -M δ ≤ x ≤ M δ (23)
where δ ∈ {0, 1} 2 . The solutions of the different optimization methods are geometrically illustrated in Fig. 4 for a different choice of the weighting vector. Thanks to the proposed weighted penalty formulation both LASSO and MILP now return a unique sparse solution, consistent with the chosen weighting vector. In particular, thanks to the weighted formulation of the LASSO, the squareshaped l 1 -norm contour lines are stretched in one direction, hence braking the parallelism with the line representing the equality constraint and thus effectively restoring the sparsity property of the LASSO.

E. Sparse post-optimization of contact forces

In the following, the proposed sparse optimization methods in ( 22) and ( 23), accounting for the block-symmetry of G, are extended to the balancing control.

In order to guarantee feasibility of the resulting optimization problem, we introduce two prioritized post-optimization stages: the one with higher priority is represented by the Quadratic Programming (QP) problem in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], while the second stage, i.e. the proposed spare regression methodology, acts in the null-space of the previous one. To this end, the equality constraint:

GF C = GF * C (24)
needs to be considered in the second stage, being F * C the optimal solution produced by the first stage [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

It is worth noticing that, thanks to the friction cone constraints in [START_REF] Stellato | Osqp: An operator splitting solver for quadratic programs[END_REF], it is inherently ensured that both the x and y components of the contact force for each contact point are simultaneously equal to zero, whenever the z-component is null. Following from this, it is sufficient to minimize the number of contact forces along the z Cartesian direction, aggregated in the vector F C,z ∈ R nC , in order to produce a solution which is sparse in all directions. As a matter of fact, this has the nice computational advantage of requiring a number of optimization variables equal to 4 n C instead of 6 n C . By now introducing the auxiliary variables t ∈ R nC , the binary variables δ ∈ {0, 1} nC and the weighting vector γ ∈ R nC , the resulting weighted LASSO formulation of the force distribution problem is given by: min

FC,t γ T t s.t. GF C = GF * C -t ≤ F C,z ≤ t b l ≤ D F C ≤ b u u l ≤ F C ≤ u u . (25)
while the weighted MILP formulation is given by:

min FC,δ γ T δ s.t. GF C = GF * C -M δ ≤ F C,z ≤ M δ, b l ≤ D F C ≤ b u u l ≤ F C ≤ u u (26)
respectively, where δ ∈ {0, 1} nC and M ∈ R + is a big enough value to guarantee that the box constraint is not active when δ i = 1.

As an additional technique to reduce the possibility of highfrequency switches in the solution of the force distribution problem, let us further introduce the following first order discrete-time dynamics:

F k+1 C = F k C + T s Ḟ k C ( 27
)
where superscript k refers to the value of the related vector at discrete time instant t k , while T s is the controller sampling time. By plugging ( 27) in ( 25) and ( 26), it is possible to obtain an equivalent formulation for the weighted LASSO and the weighted MILP force distribution problem, respectively, being now Ḟ k C the optimization vector.

V. SIMULATION RESULTS

In order to validate the proposed sparse optimization approaches to the balancing control problem, we set up a simulation scenario in Gazebo using the legged platform CENTAURO [START_REF] Baccelliere | Development of a human size and strength compliant bi-manual platform for realistic heavy manipulation tasks[END_REF], a 39 degrees-of-freedom (DoF) hybrid wheeled-legged quadruped equipped with a humanoid upper-body. The robot is fully torque-controlled by feeding back the measured linkside joint torques. The controller is written in C++ within the OpenSoT control library [START_REF] Hoffman | Robot control for dummies: Insights and examples using opensot[END_REF], and executed under the XBotCore framework [START_REF] Muratore | Xbotcore: A real-time cross-robot software platform[END_REF], within a control loop running at 1 kHz.

A. Gazebo simulation setup

The first optimization stage, that according to [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF] computes the fully actuated torque vector in [START_REF] Herzog | Momentum control with hierarchical inverse dynamics on a torque-controlled humanoid[END_REF], can be written using the Math of Task (MoT) formalism [START_REF] Hoffman | Robot control for dummies: Insights and examples using opensot[END_REF] as follows:

         i World TFoot i / World TWaist/ World TLHand + World TRHand / T Posture          << CJoint Torque Limits , (28)
where the symbol A T B denotes a Cartesian impedance task of the frame B relative to the frame A; such a task is implemented by commanding simple virtual wrenches as in the following expression:

F d = K d (x d -x) + D d (x d -ẋ) , (29)
where x, x d , ẋ and ẋd represent the actual and desired Cartesian poses and twists, respectively, whereas K d and D d denote the desired Cartesian impedance. These are combined by means of the operators "+" and "/", which are used to set aggregation and null-space relations, respectively. Finally, the symbol "<<" denotes insertion of constraints into the problem. Note that the Cartesian impedance tasks on the feet are considered at the first level in order to constrain the contact locations, i.e. to guarantee that whatever motion is produced, the feet will stay in contact with the ground. Conversely, this behaviour cannot be ensured by using soft priorities among tasks. The corresponding hierarchical QP problem has been solved using the qpOASES [START_REF] Ferreau | qpOASES: a parametric active-set algorithm for quadratic programming[END_REF] solver. According to Sec. IV-E, the fully actuated torque vector ( 13) is then mapped to an under-actuated torque vector via the proposed sparse postoptimization of contact forces. OSQP solver [START_REF] Stellato | Osqp: An operator splitting solver for quadratic programs[END_REF] has been used to solve the weighted LASSO problem in (25), while the weighted MILP regression in (26) has been solved using GLPK (GNU Linear Programming Kit) package. Finally, in all the simulations the state of the floating-base link is directly taken from the simulator.

B. Balancing under external disturbances

The simulation scenario consists on a balancing task for the CENTAURO robot under the disturbance of external forces applied on the robot waist. A constant force of 200 N is first applied downward along the z-direction for 2 s, while a second constant force of 90 N is subsequently applied sideways along the y-direction, again for 2 s, see Fig. 5. I. In the initial phase, all the approaches produce the same sparse solution, which exploits three supporting legs and is consistent with the selected weighting vector, i.e. the contact force on the rear left leg is null. When the downward external force is applied (left shaded area) the weighted LASSO, see Fig. 7(a), employs three legs until the bound on the maximum contact force is reached. At this point, a contact force is exerted also on the rear left leg, subsequently restoring the initial sparse configuration once the external disturbance stops. Note that, according to [START_REF] Stellato | Osqp: An operator splitting solver for quadratic programs[END_REF], the x and y components of the contact force are driven to zero at once with the z-component. In contrast, during the same phase, high-frequency switches between supporting legs occur with the weighted MILP, see Fig. 7(b), despite the weighted penalty formulation. Nevertheless, it is worth noticing that, thanks to the proposed first order dynamics formulation of the weighted MILP, see Fig. 7(d), it is possible to effectively remove this chattering, while ensuring sparsity. Another interesting behaviour arises during the application of the sideways external force (right shaded area). As a matter of fact, both the weighted LASSO, see Fig. 7(a), and its first order dynamics formulation, see Fig. 7(c), cannot find a sparse solution. Instead, with the first order dynamics formulation of the weighted MILP, see Fig. 7(d), a sparse solution is returned which exploits three contacts, consequently involving a switching between the rear supporting legs.

Aiming at further showing a possible application of endeffector redundancy and sparse optimization methods to singularity avoidance, in a second simulation, involving the application of a 90 N sideways force along the y-direction, the binary variables vector δ has been used to relax the feet Cartesian impedance tasks as follows:

F d Ci = δ i K d x d Ci -x Ci + D d x d Ci -ẋCi (30)
The first order dynamics formulation of the weighted MILP has been exploited for this purpose. As shown in Fig. 6, the rear right leg is automatically lifted (Fig. 6(b)), instead of being stretched (Fig. 6(a)), effectively preventing the leg from reaching a singular configuration. These results are also illustrated in the accompanying video.

VI. CONCLUSION AND FUTURE WORKS

In order to be able to exploit the inherent end-effector redundancy of a multi-legged humanoid platform, two sparse optimization approaches to the balancing control problem have been presented and discussed in this paper. The proposed optimization methods have been first formulated and compared as a general regularization approach to a Least Squares (LS) estimation problem and consequently extended to the balancing control. The discussion on the impact of the proposed approaches on the degree of sparsity of the solution has been carried out employing a simplified case study. Finally, the effectiveness of the proposed approaches has been finally validated on a simulated model of the CENTAURO robot, together with an example of a possible application of endeffector redundancy to a singularity avoidance situation.

The specification of the weighting vector γ deserve special attention and remains an important future work. The weighting vector has in fact the prominent role to bias the sparsity pattern towards a preferable one, e.g. think of a bipedal robot standing on the left foot instead of on the right one. Note also that the proposed sparse methodologies constitute a real-time reactive control layer. As such, they are meant to be used in conjunction with a high-level (not real-time) planner, which, based on tasks and environment knowledge, could be responsible for the specification of the weighting vector, as well as as for possible changes in the contact points' locations, e.g. in order to move a foot from the floor to a wall to help the arms pushing an heavy object. Finally, validation on the real robotic platform remains a forthcoming work.

Fig. 2 .

 2 Fig. 2. The black solid line represents the constraint a T x = b, while red and black dashed contour lines refer to the l 2 -norm cost function (3) and to the l 1 -norm cost function (11), respectively. The possible solutions of the corresponding optimization problem (circle markers) are located at the intersection of the equality constraint line with the innermost contour line.

Fig. 4 .

 4 Fig. 4. Graphical interpretation for a different choice of the weighting vector: γ = 1 1.5 T and γ = 1.5 1 T for the left and right figure, respectively. The l 1 -norm contour lines are stretched along one axis according to γ, breaking the parallelism with the line representing the equality constraint and thus effectively restoring the sparsity of LASSO.

Fig. 5 .

 5 Fig. 5. Screen shots from CENTAURO simulations in Gazebo: applied external forces on the robot waist. The weighting vector has been set equal to: γ = 1 1 1 1.5 , i.e. penalizing the contact force applied on the rear left leg. The desired Cartesian impedance for the waist position task is given by: K d = 500 500 500 N m and D d = 200 200 200 N s m , while a friction coefficient µ = 0.3 has been selected for the linearized friction cones. The following upper and lower bounds have been considered for the contact forces: u u = 500 500 500 N and u l = -u u . Finally, M has been set equal to 1000. The time histories of the contact forces, obtained from the proposed sparse optimization methods, are shown in Fig. 7, while the corresponding computation times are reported in TableI. In the initial phase, all the approaches produce the same sparse solution, which exploits three supporting legs and is consistent with the selected weighting vector, i.e. the contact force on the rear left leg is null. When the downward external force is applied (left shaded area) the weighted LASSO, see Fig.7(a), employs three legs until the bound on the maximum contact force is reached. At this point, a contact force is exerted also

Fig. 6 .

 6 Fig. 6. Application of end-effector redundancy to singularity avoidance: the rear right leg is lifted during a balancing task under the disturbance of a sideways external forces on the robot waist.

Fig. 7 .

 7 Fig. 7. Comparison among the proposed sparse optimization methods: time histories from CENTAURO balancing simulations under the disturbance of external forces applied on the robot waist. A constant force of 200 N is first applied downward along the z-direction (left shaded area), while a second constant force of 90 N is subsequently applied sideways along the y-direction (right shaded area).

TABLE I COMPUTATION

 I TIME OF THE PROPOSED SPARSE OPTIMIZATION METHODS (MEAN, STANDARD DEVIATION AND MAXIMUM VALUE).

	Method Computation time [ms]	mean	std	max
	Weighted LASSO	3.79	0.83 14.84
	Weighted LASSO first order dynamics	4.05	1.12 30.81
	Weighted MILP	3.41	0.94 15.00
	Weighted MILP first order dynamics	2.82	0.76 13.03

This paper was recommended for publication by Editor D. Song upon evaluation of the Associate Editor and Reviewers' comments. This work was supported by the European Union Horizon 2020 research and innovation program under grant agreements No. 644839 (CENTAURO) and No. 644727 (CogIMon).