J Erôme Berjonneau

Laetitia Colombel

Jacques Poirier

Michel Pichavant

Franc -Oise Defoort

Jean-Marie Seiler

Determination of the Liquidus Temperatures of Ashes from the Biomass Gazification for Fuel Production by Thermodynamical and Experimental Approaches

During the gasification of biomass, inorganic species are produced and constitute an important obstacle in this process. The aim of this study is to understand the high temperature behavior of the inorganic species particularly the ashes close to the miscanthus and the straw ashes. A thermodynamic approach is used to obtain this characteristic. The database used for the calculation is first validated by comparison with systems available in the literature and by using an experimental approach. This study reveals the lack of optimization and the lack of data concerning alkali species and particularly K 2 O-rich compositions. An experimental approach has been used to determine the liquidus temperature of ashes by a classical quenching method. Experiments using sealed capsules largely prevent the volatilization of K 2 O. The liquidus temperature of three compositions with alkali and one without alkali were measured. A liquid miscibility gap is observed experimentally and by calculation for the straw ash.

Introduction

Biomass is one of the most important renewable energy sources. The use of biomass energy has the potential to greatly reduce greenhouse gas emissions. It is also a way to increase fuel diversification and reduce dependence on foreign oil.

The context of this research is the gasification of biomass for fuel production in an entrained flow reactor. 1 During the gasification of biomass, inorganic species are produced and constitute an important obstacle in this process. Inorganics are the mineral elements and compounds present in biomass apart from the main gas and organic species. They play an important role on fluidized bed agglomeration, on the corrosion of the reactor wall, on the gas-phase pollution, or on the corrosion of the structure.

The overall objective of this work is understanding the behavior of inorganics under entrained-flow gasification conditions. Once gasification has taken place, the remaining solid, mainly ashes, will soften, liquefy and deposit on the walls, hereby forming the slag. A layer of solidified and molten ash will be condensed on a refractory wall. This wall is cooled externally. This ash layer will play the role of a thermal insulating lining to minimize the heat losses and to protect the refractory layer against external attacks. 2 The slag has to be liquid to be able to run down the wall.

This paper presents the results of the characterization of slag behavior, particularly the liquidus temperature. The study will be carried out for two typical compositions of ashes (Table 1).

The first one, which is rich in SiO 2 (C1), is close to the miscanthus ash, and the second one, with the same amount of SiO 2 and CaO (C2), is close to the straw ash contents. Apart from CaO and SiO 2 ,K 2 O is an important component of the ash. The selection of these two types of ashes is justified by the great potential of miscanthus and straw for biomass valorization.

In the first part of this paper, the thermodynamic approach is described. The thermodynamic tool used is described and is validated by comparison with existing data from the literature. The liquidus temperatures calculated for the two typical ashes are then presented.

Second the experimental approach is described. The results for the ternary system (SiO 2 -CaO-K 2 O) and for the two typical compositions (C1 and C2) are presented. Finally, the thermodynamic and experimental approaches are compared.

The same kind of research combining experimental and thermodynamic modeling studies has already been successfully done on coal combustion slag. 3 2. Thermodynamic Approach 2.1. Description of the Thermodynamic Calculations. The calculation is based on the minimization of the Gibbs energy of the system. It is thus possible to deduce the nature of the solid, liquid and gaseous phases as well as their chemical compositions at the thermodynamic equilibrium. This knowledge makes it possible to deduce the liquidus and solidus temperatures and the amount of the solid and liquid phases. *To whom correspondence should be addressed. Telephone: þ33.2.38.25.55.14. Fax: þ 33.2.38.63.81.03. E-mail: Jacques.Poirier@ univ-orleans.fr.

(1) Jinsong, Z; Qing, C; Hui, Z; Xiaowei, C; Qinfeng, M; Zhongyang, L; Kefa, C. Biomass-oxygen gasification in a high-temperature entrained-flow gasifier. Biotechnol. Adv. 2005, 25 (5), 606-611.

(2) Poirier, J; Colombel, L; Prigent, P. The corrosion mechanisms of SiC refractory lining in waste incineration plants and in reactors of biomass gasification. Proc. Int. Colloq. Refract., 52nd 2009,49-52.

(3) Jak, E; Degterov, S; Zhao, B; Pelton, A. D.; Hayes, P. C. Coupled experimental and thermodynamic modeling studies for metallurgical smelting and coal combustion slag systems. Metall. Mater. Trans. B 2000, 31 (4), 621-630.

Calculations have been carried out with Factsage (version 5.4.1).

Factsage is a fully integrated database and software package developed jointly between Thermfact/CRCT (Montreal) and GTT-Technologies (Aachen). It consists of a series of modules that access and manipulate thermodynamic databases and perform various calculations. Factsage is composed of two types of thermo chemical databases: compounds (i.e., pure substances, over 4,400 compounds) and real solutions.

• Compound databases are mainly for stoichiometric solid, liquid and gaseous species. Allotropes (for example cristobalite, tridymite) are included. Nonstoichiometric compounds such as FeO x (x =0.947)arealso stored. • Solution data are for solid and liquid alloys, ceramics, salts, mattes, slags, etc. The data are stored in the form of Gibbs energy functions for the phase constituents and parameters for the excess Gibbs energy of mixing between the phase constituents. Factsage supports several different solution models including simple polynomial models such as the Redlich-Kister or the Legendre polynomial combined with different higher order extrapolations (Muggianu, Kohler, Toop), the quasichemical model, Pitzer parameters, and sublattice models. 4 The associated databases for this study are: ELEM (elements thermodynamic database), FACT 53 (gas species, solid and liquid compounds thermodynamic database), and FT oxid (compounds and solutions for oxides database).

The expected equilibrium phases were predicted by using the "Equilibrium" module of Factsage considering as possible phases: the slag, several solid-solution phases and all liquid and solid compounds between the inputs. Oxide compositions are entered as initial compositions in weight percent in the submenu "reactants" of "Equilib".

The FToxid database contains data for pure oxides and oxide solutions of 20 elements. All the binary and ternary subsystems have not been evaluated and optimized, nor are all composition ranges covered. Subsystems which have not been evaluated and optimized have been assumed ideal or have been approximated.

Thus, this database contains the following major oxides: Al 2 O 3 ,CaO,Fe 2 O 3 ,MgO,SiO 2 , and all the binary and many ternary and quaternary subsystems more or less optimized. Interaction between these five oxides should be correctly calculated.

When Na 2 O and K 2 O are present, the Factsage documentation specifies that optimization is less accurate. Major improvements are in progress for these oxides. 5,6 To conclude, the use of the thermodynamic FToxid database for the calculation of the liquidus temperature for the two compositions will be marred by some uncertainties due to the alkalis.

To ensure the reliability and the pertinence of the calculations, it seems necessary to validate the thermodynamic calculations for simpler systems with fewer oxides than the multi component ashes. Consequently, the main subsystems (binary and ternary) have been tested.

2.2. Validation of the Thermodynamic Calculations. The thermodynamic database FToxid is constructed from the optimization of several binary and ternary systems. A phase diagram is said to be optimized when all the data from the literature have been analyzed, evaluated, selected, and used to determine the main parameters to carry out the modeling of the phase diagram. The modeling is not always carried out for the complete compositional range of a system. To validate the FToxid database, calculations have been done in the most pertinent binary, ternary and quaternary systems.

To assess the validation of the FToxid database, the systems have been checked in two ways:

• For the optimized diagrams, the optimization has been checked in the studied domain of ashes C1 and C2. If the range of composition is not optimized, bibliographic research is done to check the existence of other measurements and to reproduce them by calculation from the database. • For the nonoptimized diagram, bibliographic research is performed to find experimental results that are compared (5) Pelton, A. D. Thermodynamic calculations of chemical solubilities of gases in oxide melts and glasses. Glastech. Ber. 1999, 72 (7), 214-226.

(6) Pelton, A. D.; Wu, P. Thermodynamic modeling in glass-forming melts. J. Non-Cryst. Solids. 1999, 253, 178-191. with the calculation. This comparison gives information about the validation of the database.

For the four main oxides SiO 2 , CaO, Al 2 O 3 and MgO, all the binary systems have been checked with the same proportions as the two ashes. [7][8][9][10][11][12] This has not been done for minor oxides such as Fe 2 O 3 and Na 2 O.

For K 2 O, which is also a major oxide, K 2 O-SiO 2 is the only optimized binary system. 13 It has been successfully checked for the two ashes compositions C1 and C2. The two other main binary systems K 2 O-CaO and K 2 O-MgO are not optimized in the FToxid database. Experimental data have not been found in the literature. It makes the validation of the database in this domain impossible.

On the other hand, K 2 O-Al 2 O 3 is optimized 7 but not in the composition range relative to ashes C1 and C2 (Figure 1). Besides, experimental data has not been encountered in the literature to validate the database in this domain.

Concerning the ternary systems resulting from the major oxides SiO 2 , CaO, K 2 O, MgO, only SiO 2 -CaO-MgO is optimized 14 and successfully checked for the two ashes C1 and C2. This result is in agreement with the literature. 15,16 No data has been published for CaO-K 2 O-MgO system. Experimental data has been encountered in the literature for SiO 2 -K 2 O-MgO 17 and SiO 2 -CaO-K 2 O. 18 Validation of the database has been poorly obtained for these two ternary systems. Particularly the liquidus temperature is higher 200-400 °C for C1 and C2 as illustrated, for example, for C1 in SiO 2 -CaO-K 2 O (figure 2).

2.3. Calculation Results for the Typical Ashes. The liquidus temperatures calculated with Factsage are presented in Table 2 for the two ash compositions with and without Na 2 O and K 2 O.

2.3.1. Composition C1. Without Na 2 O and K 2 O, the liquidus temperature is 1560 °C (Figure 3). The primary crystallization phase is SiO 2 . With Na 2 O and K 2 O, the liquidus temperature (1324 °C) is lower by around 230 °C (16) Kondratiev, A; Jak, E. Prediction coal ash slag flow characteristics (viscosity model for the Al 2 O 3 -CaO-FeO-SiO 2 system). Fuel 2001, 80, 1989-2000. (17) Roedder, E. W. in agreement with the existence of a eutectic in the binary system with silica (Figure 4). The calculation predicts the formation of two immiscible liquid phases called slag 1, which is rich in SiO 2 , and K 2 O and slag 2, which is rich in SiO 2 and CaO. The crystallization phases are composed of melilite (Ca 2 MgSi 2 O 7 ), clinopyroxene (CaMgSi 2 O 6 ), and wollastonite (CaSiO 3 ).

2.3.2. Composition C2. Without Na 2 O and K 2 O, the liquidus temperature is 1420 °C (Figure 5).

The primary crystallization phase is wollastonite. With Na 2 O and K 2 O, the liquidus temperature (1256 °C) is lower by about 170 °C (Figure 6). The calculation predicts the presence of two immiscible liquid phases called slag 1, which is rich in SiO 2 , and CaO and slag 2, which is rich in SiO 2 and K 2 O. The primary crystallization phase is a spinel phase close to Fe 3 O 4 in composition.

For both compositions, the modeling of Na 2 O and K 2 O with CaO, which is a major oxide, has not been optimized. Thus the liquidus temperature is quite uncertain.

Besides it is important to remember that ternary systems with K 2 O (SiO 2 -CaO-K 2 O and SiO 2 -MgO-K 2 O) are poorly validated by Factsage. The calculations are probably less reliable because of the lack of optimization for this oxide and both liquidus temperatures from calculation are probably too high.

To test the validity of these results it is necessary to compare them with experimental measurements.

Experimental Approach

3.1. Experimental Procedure. Experiments have been carried out in a vertical quench furnace at ISTO (Institut des Sciences de la Terre; CNRS laboratory, Orl eans, France) (Figure 7). The objectives of the experiments are to determine an experimental liquidus temperature and to identify the phases crystallizing below the liquidus. Experiments were performed using the quenching method. 19 The starting sample is made of ground glass mixed with polyvinyl alcohol to ensure its cohesion during the installation in the furnace and which is later vaporised at low temperature. Two different procedures were used: one in an open system and the second one in sealed capsules. In the first one, an open system, the slurry obtained is then spread over a platinum hook (O e = 0.2 mm) hung on a wider platinum wire (O e = 0.5 mm) inserted into an alumina rod. A thermocouple (Pt-Pt Rh 10 wt.%) inserted into another alumina rod is located close to the sample to measure the temperature. For the second one, sealed capsules, the samples have been enclosed in sealed capsules to minimize the volatilization of the alkali oxides (K 2 O and Na 2 O) during the experiment. 20 Around 40 mg of powder is prepared. Capsules are made from a platinum tube of 2 mm in diameter and 10 mm in length. These capsules are cleaned, dried, and then arc welded (Figure 8).

For both procedures, the temperature distribution is uniform.

Daily fluctuation is less than to 2 °C and temperature is known to better than ( 5 °C. The whole assembly is inserted into a vertical furnace at a steady temperature for 3 h, under air.

Quenching is obtained by applying a high voltage to the platinum wires. This fuses the thinner wires (O e =0.2 mm). As a result, the sample falls into the cold part of the furnace, allowing the melt phases to quench into a glass or a metastable quench crystallized phase. Once the sample has been quenched, the capsule is opened (if sealed capsule procedure). The glass is extracted and then polished for SEM (scanning electron microscopy) observations. WDS (wavelength dispersive spectroscopy) electron microprobe analysis is used for the determination of the composition of the experimental products.

The possible formation of crystals indicates if the sample was completely liquid at the temperature of the experiment. The liquidus temperature is obtained by bracketing the disappearance of crystals by means of experiments at different temperatures.

The experimental approach has been first validated with the ternary system K 2 O-CaO-SiO 2 . Second this approach has been applied to both ash compositions (C1 and C2).

3.2. Validation of the Method with the Ternary System K 2 O-CaO-SiO 2 . To validate this method, one composition called KCS was studied in the ternary system K 2 O-CaO-SiO 2 ,in oxide proportion similar to the C1 composition (Table 3).

This ternary diagram has not been optimized in the FToxid database. However experimental data exist in the literature. 18 For this composition, the liquidus temperature, calculated with Factsage, is 1418 °C (Figure 4). It is higher than the one measured in the publication, 18 which is around 1175 °C.

The ternary composition has been prepared with K 2 CO 3 , CaCO 3 , and SiO 2 oxides. Before the experiments, the powder was previously decarbonated by step heating at 700 and 900 °C. At the end of the experiment, a hollow ball with a diameter of about one millimeter is obtained, then polished and observed with SEM.

The cross section of a glass ball, quenched at 1175 °C, reveals a ring of glass in which crystals, with a size of about 4 μm, occur (Figure 9).

For samples quenched between 1177 and 1225 °C, crystals are present for this specific composition. The crystals grow in size as the temperature increases, while their number decreases (Figure 10). At 1230 °C, no crystals are observed. The liquidus temperature is bracketed between 1225 and 1230 °C. The experimental temperature is close to the literature value 18 and differs from the calculated value (Figure 2). This result confirms the inability of the FToxid database to reproduce phase equilibria involving K 2 O.

The chemical analysis reveals that the crystals are wollastonite (CaSiO 3 ), in agreement with the literature and the thermodynamic results.

3.3. Composition C1. This composition has been measured by the two procedures (open system and sealed capsule) 3.3.1. without the Sealed Capsule (Open System). This composition has been prepared by Schott (Table 4). The experimental composition (analyzed by electron microprobe) is in agreement with the miscanthus ash C1.

An iterative approach has allowed estimation of the liquidus temperature within a range of 5 °C.

Below the liquidus temperature, crystals in the shape of needles are observed (Figure 11).

The number of crystals decreases as temperature rises while their size increases. The chemical analysis shows that crystals are wollastonite CaSiO 3 . The glass quenched at 1240 °C does not contain any crystals. The liquidus temperature is therefore bracketed between 1235 and 1240 °C.

The evolution of crystal and glass phases is followed with the electron microprobe as a function of the temperature. Data were obtained using a microscope image acquisition and analysis software package (Table 5).

The amount of crystal decreases as the temperatures rises, consistent with a closer approach to the liquidus temperature. This analysis shows the volatilisation of K 2 O with increasing temperature (Table 5 and Figure 12). Up to 1240 °C, the volatilization of K 2 O is weak. At 1240 °C, 1% wt. of the glass is volatilized which represents a loss of 10% wt. of the K 2 O charge. Even if the experiment could be improved to avoid volatilization, this measure gives a good estimate of the liquidus temperature.

3.3.2. with Sealed Capsule. As in the former experiment without the sealed capsule, an iterative approach, by means of several experiments at different temperatures, has made it possible to estimate the liquidus temperature within a range of 5 °C. The liquidus temperature is around 1230 °C (Figure 13).

The evolution of crystal and glass phases is followed with the electron microprobe as a function of the temperature (Table 6) and K 2 O is reported on the same graph as the experiments without capsules (Figure 12). The loss of K 2 O is less important with the use of sealed capsule (less than 0.5 wt. % of the total K 2 O).

In both cases, the loss of K 2 O is quite weak. The liquidus temperature is confirmed and is around 1230 °C. This value is one hundred degrees lower than the thermodynamic calculation for alkali species which is 1324 °C (table 2). The composition of the liquid phase at the liquidus temperature is presented in table 7.

The primary crystallization phase predicted by thermodynamics, the melilite (Ca 2 MgSi 2 O 7 ) followed with clinopyroxene (CaMgSi 2 O 6 ) and then wollastonite (CaSiO 3 ), is only in agreement with experimental results for the formation of the wollas-tonite. In addition, the presence of two immiscible liquid phases predicted by the calculation is not observed.

3.4. Composition C2. The composition C2 has been prepared by Schott (Table 8). This composition has been measured by the two procedures (open system and sealed capsule) 3.4.1. without Sealed Capsules (Open System). The composition C2 is richer poorer in SiO 2 and richer in K 2 O and CaO. In this range of temperature and without a capsule, the experiments have shown significant vaporisation of K 2 O at a temperature bellow 1280 °C. For the real composition, that is, without alkali species, as total K 2 O vaporization has occurred, the liquidus temperature has been measured around 1375 °C (Figure 14). This value is in reasonable agreement with the liquidus temperature predicted from the calculation without alkali species: 1420 °C (Table 2) and the primary crystallization phase is wollastonite CaSiO 3 (Figure 3). The evolution of crystal and glass phases is presented in Table 9.

3.4.2. with Sealed Capsules. At 1202 °C and at 1230 °C, two immiscible liquid phases are present together with a crystal phase From 1230 °C, the two immiscible liquid phases remain, but the crystal phase is no longer observed (Figure 15). A mass balance, performed by microprobe analyses, shows that the quench crystallized phases observed on the figure 15 have been formed upon quenching of the melt which makes the determination of the liquidus temperature less straightforward.

The compositions of these two immiscible liquid phases and the crystal phase are presented in Table 10 andFigure Figure 17 shows the evolution of the ratio K 2 O/K 2 O þ CaO of the three phases between 1200 and 1400 °C. The melting temperature ranges between 1230 and 1295 °C.

This result is confirmed by DSC (differential scanning calorimetry; Figure 18). These results are in agreement with the thermodynamic calculations which predict that this temperature is 1256 °C with the disappearance of the solid phase.

Discussion

Measurements and calculations of liquidus temperatures of two typical ashes from biomass gasification for fuel production has been compared in this publication and are summarized in Table 11.

The liquidus temperature has been calculated with the thermodynamic software Factsage for two typical ash compositions. To assess the reliability of this calculation a validation step has been used. It reveals the reliability of the database concerning the studied system without alkalis. This excellent agreement between F*A*C*T equilibrium calculations and measurements has already observed by Jack. 15,16 However, with alkalis, there is a lack of thermodynamic data and the database is very uncertain, which clearly is indicated in the F*A*C*T instructions. Indeed, with K 2 O, only one binary system K 2 O-SiO 2 13 is optimized. For the ternary system resulting from the major oxides with K 2 O, no optimization is performed in the database, but measurements are published in the literature for SiO 2 -K 2 O-MgO 17 and SiO 2 -K 2 O-MgO. 18 For the others, the comparison with the literature exhibits strong differences which reveal the limitations of the FToxid database. Thus, the liquidus temperatures obtained from the calculation are not quite certain. To make up for the lack of thermodynamic data, some updating is required. However, this procedure depends on the experimental measures available in the literature which are unfortunately very scarce for these systems and particularly for K 2 O.

An experimental approach has been used in this paper to get a liquidus temperature measurement and to compare it with the calculation. The experiments have been carried out with and without sealed capsules.

Without capsules, the first system studied is the ternary K 2 O-CaO-SiO 2 (KCS), which made it possible to validate the method used thanks to the comparison with existing literature data. Very little K 2 Oi sv o l a t i l i z e d .H o w e v e ra 200 °C discrepancy with calculation (Table 11) confirms the lack of optimization of the FToxid database as previously discussed. This experimental method (without capsule) has been applied to the multi component composition C1 containing seven oxides. This result has given a first approximation of the liquidus temperature for C1, which has been confirmed by measurement with the sealed capsule method. Because of the large differences with the calculation (∼100 °C), the lack of optimization (for FToxid database) is again confirmed. Besides, the primary crystallization phase predicted is not in agreement with the one experimentally obtained and no miscibility gap in the liquid phase is observed experimentally contrary to calculations.

This procedure without capsule could not be applied for the multicomponent C2 composition, which is richer in K 2 Othan C1. Because of K 2 O volatilization with increasing temperature, all the K 2 O is volatilized at the liquidus temperature of C2 (1375 °C). Nevertheless, these experimental results are in agreement with thermodynamic predictions for composition without alkali species. The liquidus temperatures are reasonably close (Table 11). There is a difference of less than 50 °C, and the primary crystallization phase is wollastonite as calculated. This result validates the FToxid database without alkali species as already observed by Jack. 15,16 For the composition C2, the use of capsule gave reliable results. The liquidus temperature ranges between 1230 and 1295 °Ca n dn oK 2 O loss has occurred. The samples exhibit also the formation of two immiscible liquid phases. This is in agreement with thermodynamic predictions (1256 °C) taking into account alkali species. The good accordance between thermodynamic calculations and measurements for this C2 multicomponent composition is very surprising because this composition is richer in K 2 Oa n d because the FToxid database has been found to be poorly validated when alkalis are present. Thus, it should not be concluded that the FToxid is well validated for this C2 ash. Indeed, one must conclude that an opposite effect might be compensated in the database.

Conclusion

This paper has shown a good agreement between thermodynamic calculations and measurements of liquidus temperatures and primary phases of crystallization for ashes without alkali as already observed in the literature. 15,16 However the discrepancy between the measured and calculated liquidus temperatures of systems containing alkalis (K 2 O) is between several hundred and 0 °C. Thus, if thermodynamic calculations need to be used to predict the liquidus temperature of typical ashes from biomass gasification for fuel production, it is necessary to improve the FToxid database where alkalis are present as it is clearly indicate in the F*A*C*T instructions. However, the literature provide too few measurements that can be taken into account by the FToxid database developers. Such measurements with alkalis are not easy to perform because of their high volatility. Nevertheless an experimental set up has been presented in this paper which is able to produce reliable results with binary or ternary systems and multicomponent systems (7 oxides). Only the ternary result presented in this paper could be taken into account to improve the database. The two multicomponent compositions C1 and C2 measured in this study can only be used to validate a thermodynamic database as too many oxides are present and some compensation error may occur as is probably the case with the composition C2.

Acknowledgment. The authors thank ANR (Agence National de la Recherche, France) and the PNRB (Programme National de Recherche sur les Bio energies) for their helpful financial contribution in this research.

Figure 1 .

 1 Figure 1. Binary diagram K 2 O-Al 2 O 3 [4].

Figure 2 .

 2 Figure 2. Ternary diagram SiO 2 -CaO-K 2 O [13].

  (14) Jung, I. H.; Decterov, S. A.; Pelton, A. D. Critical thermodynamic evaluation and optimization of the CaO-MgO-SiO 2 system. J. Eur. Ceram. Soc. 2005, 25 (4), 313-333. (15) Jak, E. Prediction of coal fusion temperature with the F*A*C*T thermodynamic computer package. Fuel 2002, 81, 1655-1668.

Figure 4 .

 4 Figure 4. Weight percent of the phases formed as a function of the temperature for C1 with Na 2 O and K 2 O.

Figure 3 .

 3 Figure 3. Weight percent of the phases formed as a function of the temperature for C1 without Na 2 O and K 2 O.

Figure 5 .

 5 Figure 5. Weight percent of the phases formed as a function of the temperature for C2 without Na 2 O and K 2 O.

Figure 6 .

 6 Figure 6. Weight percent of the phases formed as a function of the temperature for C2 with Na 2 O and K 2 O.

(

  20) Matousek, J. Chemical reactions taking place during vaporization from silicate melts. Ceram.-Silik. 1998, 42 (2), 74-80.

Figure 7 .

 7 Figure 7. Quenching method.

Figure 8 .

 8 Figure 8. Sealed capsule method.

  16 (K 2 O-CaO-SiO 2 ternary diagram).

Figure 10 .

 10 Figure 10. Glass balls quenched at 1225 and 1230 °C for the ternary composition.

Figure 9 .

 9 Figure 9. Glass balls quenched at 1175 °C for the ternary composition.

Figure 11 .

 11 Figure 11. Glass balls quenched at 1235 and 1240 °C for composition C1 without sealed capsule.

Figure 12 .

 12 Figure 12. Evolution of K 2 O with and without capsule as a function of the temperature in comparison with the initial composition C1.

Figure 13 .

 13 Figure 13. Glass balls quenched at 1228 and 1234 °C for composition C1 with sealed capsule.

Figure 14 .

 14 Figure 14. Glass balls quenched at 1372 and 1376 °C for composition C2 without sealed capsule.

Figure 15 .

 15 Figure 15. Glass balls quenched between 1202 and 1385 °C for composition C2 with sealed capsule.

Figure 16 .

 16 Figure 16. Composition (% wt.) of the phases 1, 2, and 3 (K 2 O-CaO-SiO 2 ternary diagram).

Figure 17 .

 17 Figure 17. Compositions of three phases in a K 2 O-CaO-SiO 2 ternary diagram.

Figure 18 .

 18 Figure 18. Differential scanning calorimetry (DSC) of C2 (under argon, T max = 1600 °C).

Table 1 .

 1 Typical Composition of Ashes

		C1 (wt. %)	C2 (wt. %)
	SiO 2	60	34
	Al 2 O 3	51
	Fe 2 O 3	35
	CaO	13	34
	MgO	4	3
	K 2 O1 4	2 2
	Na 2 O1	1
	total	100	100

Table 2 .

 2 Liquidus Temperatures for C1 and C2 with and without Na 2 O and K 2 O

		C1 rich in silica (miscanthus)	C2 rich in K 2 O (straw)
		with Na 2 OK 2 O
	T liquidus	1324 °C	1256	°C
		without Na 2 OK 2 O
	T liquidus	1560 °C	1420	°C

Table 3 .

 3 Studied Composition in the K 2 O-CaO-SiO 2 System

		ternary composition (% wt.)	composition C1 (% wt.)
	SiO 2	69	60
	CaO	15	13
	K 2 O1 6	1 4
	Al 2 O 3	05
	Fe 2 O 3	03
	MgO	0	4
	Na 2 O0	1

Table 4 .

 4 Theoretical and Experimental Compositions of C1

	composition C1 siliceous slag close	composition C1 given by
		to miscanthus ash (% wt)	Schott and analyzed (% wt)
	SiO 2	60	61.3
	Al 2 O 3	5	4.9
	Fe 2 O 3	3	2.6
	CaO	13	12.8
	MgO	4	3.9
	K 2 O	14	13.5
	Na 2 O	1	1.1

Table 5 .

 5 Evolution

	% wt. glass	89.1	96.7	99.8	100	100
	% wt. crystal	10.9	3.3	0.2	0	0
	% wt. K 2 O	13	11.75	12.2	11.8	9.6
	(glass þ crystal)

of the Amount of Glass and Crystal As a Function of the Temperature Composition C1 without the Sealed Capsule (Open System) 1170 °C 1216 °C 1235 °C 1240 °C 1278 °C

Table 6 .

 6 Evolution of the Amount of Glass and Crystal As a Function of the Temperature Composition C1 with the Sealed Capsule

		1198 °C	1224 °C	1228 °C	1233 °C	1234 °C	1240 °C	1254	°C
	% wt. glass	90.9	95.9	97.0	100	100	100	100
	% wt. crystal	9.1	4.1	3	0000
	% wt. K 2 O (glass þ crystal)	13.41	13.31	13.38	13.29	13.36	13.29	13.31

Table 7 .

 7 Composition of Glass at the Liquidus Temperature (Experimental Values)

		K 2 O	SiO 2	FeO	CaO	Na 2 OA l 2 O 3	MgO
	% wt.	13.29	61.63	2.47	12.71	1.02	4.92	3.96

Table 8 .

 8 Theoretical and Experimental Compositions of C2

		composition C2 close to	composition C2 given by Schott and
		straw ash (% weight)	analyzed (% weight)
	SiO 2	34	34.8
	Al 2 O 3	1	1.7
	Fe 2 O 3	5	4.3
	CaO	34	33.2
	MgO	3	5.1
	K 2 O	22	20.3
	Na 2 O	1	0.6

Table 9 .

 9 Evolution of the Amount of Glass and Crystal As a Function of the Temperature Composition C2 without the Sealed Capsule (Open System) 1296 °C 1316 °C 1348 °C 1357 °C 1362 °C 1372 °C 1376 °C

	% wt. glass	66.7	72.5	77.5	82.5	84.1	87.8	100
	% wt. crystal 33.3	27.5	22.5	17.5	15.9	12.8	0

Table 10 .

 10 Composition of the Two Immiscible Liquid Phases and the Solid Phase C2 with Sealed Capsules

	temperature		K 2 O (% wt.)	SiO 2 (% wt.)	FeO (% wt.)	CaO (% wt.)	Na 2 O (% wt.)	Al 2 O 3 (% wt.)	MgO (% wt.)
	1230 °C	liquid phase 1	4.5	34.2	0.4	57.9	0.3	0.04	0.9
		liquid phase 2	40.9	30.9	7.7	1.4	0.7	4.4	12.9
		solid phase 3	24.7	37.4	7.4	19	1	1.3	3.7
	1295 °C	liquid phase 1	3.9	34.4	0.35	58.4	0.2	0.03	1.05
		liquid phase 3	25.4	36.9	6.93	19.3	1.1	1.5	4.5

Table 11 .

 11 Liquidus Temperatures for KCS, C1, and C2 T liquidus (°C) < 1240, without capsule 1228 < T liquidus (°C) < 1234, with capsule C2 1420 °C, without alkalis 1372 < T liquidus (°C) < 1376, without capsule, that is, without alkalis no miscibility gap in the liquid phase no miscibility gap in the liquid phase 1256 °C, with alkalis 1230 < T (°C) < 1295, with capsule, that is, with alkalis presence of miscibility gap disappearance of the solid phase, presence of miscibility gap

		thermodynamic calculations	experimental measurements
	KCS	1418 °C	1225 < T liquidus (°C) < 1230
	C1	1342 °C, with alkalis	1235 <


Antonin Dallard

Mehdi Benallegue

Fumio Kanehiro

Fellow, IEEE Abderrahmane Kheddar

Synchronized Human-Humanoid Motion Imitation

We present a tele-operation control framework that (i) enhances the upper motion synchrony between a user and a robot using the minimum-jerk model coupled with a recursive least-square filter, and (ii) synchronizes the walking pace by predicting user's stepping frequency using motion capture data and a deep learning model. By integrating (i) and (ii) in a task-space whole-body controller, we achieve full-body synchronization. We assess our humanoid-to-human wholebody synchronized motion model on the HRP-4 humanoid robot in forward, lateral and backward walks with concurrent upper limbs motions experiments.

Whole-Body

Control

I. INTRODUCTION

Synchronized human motion such as those witnessed in artistic swimming or dancing choreography performances come after an intensive technical training. Those motions are not only precisely predefined but they are strictly timed with auditive or visual rhythms to ensure synchrony. Much looser forms of synchrony occur in other social contexts [START_REF] Wiltermuth | Synchrony and cooperation[END_REF].

On the other hand, making robots produce synchronized motions is 'simpler'. This has been demonstrated with a large number of small size humanoids 1 or quadrupeds2 or both 3 . Moreover, it can be done even without predefined trajectories with good network communication protocol and dedicated controllers, see e.g., [START_REF] Chung | Cooperative robot control and concurrent synchronization of lagrangian systems[END_REF], [START_REF] Bechon | Synchronization and quorum sensing in a swarm of humanoid robots[END_REF].

Yet when it comes to synchronizing motions between a human and a (humanoid) robot, the problem becomes more challenging, e.g., [START_REF] Miyake | Interpersonal synchronization of body motion and the walk-mate walking support robot[END_REF]. When the trajectory is predefined, we need to account for various discrepancies such as dissimilar anthropomorphic features of all kind (i.e., in terms of size, range of motion, degrees of freedom, torque limitations, etc.). Finally the human has to be able to perform motions perfectly to ensure synchrony during the execution, see the example of the dance performance reported in [START_REF] Nakazawa | Imitating human dance motions through motion structure analysis[END_REF], [START_REF] Nakaoka | Learning from observation paradigm: Leg task models for enabling a biped humanoid robot to imitate human dances[END_REF].

However, to our best knowledge none existing work has tackled the problem of whole body online motion synchronization, which is certainly a type of imitation but where the start and end of critical motion timing prevails and motion lag of both actors is to be brought close to zero. Yet, we report in Sec. II few exceptions having links to this problem.

We suggest that human-humanoid synchronized motion imitation requires the following ingredients:

• knowledge on the human motion, which can be provided by any low-lag human tracking technology 4,5 ; • anticipating on the start of a new motion, i.e., from a zero velocity state limbs; because of the causality of the motion, it is impossible to bring the start-motion lag between the human and the robot to zero, unless an early intention of the motion can be observed, e.g., using brain computer interfaces [START_REF] Gergondet | Using brain-computer interface to steer a humanoid robot[END_REF]; • anticipating on the end-time of an ongoing motion; here a good prediction can theoretically put the lag in the motion synchronization to nil. To realize human-humanoid synchronized motion, we devise two separate imitation-mapping strategies. The first deals with reproducing arm and attitudes synchronized imitation, Sec. III. The second focuses on steps synchronized motion, that is to say, humanoid footstep and walking are to be synchronized to those of the human, Sec. IV. Experiments with the HRP-4 humanoid robot assess the effectiveness of our approach together with its current limitations, see Sec. V.

II. BACKGROUND AND PROBLEM STATEMENT A. Whole Body Imitation

Achieving motion imitation between a human and a robot relies first on obtaining a similar posture between both entities. Kinematics retargeting, in the tasks space or in the joint space (or hybrid) [START_REF] Darvish | Whole-body geometric retargeting for humanoid robots[END_REF] allows a humanoid robot to mimic a human posture. For floating base robots, balance must be considered and applied in the control scheme, either as a balance constraint [START_REF] Penco | A multimode teleoperation framework for humanoid locomanipulation: An application for the iCub robot[END_REF] or as an imitation of the human state. In [START_REF] Benallegue | A computational model for synchronous motion imitation by humans: The mirror controller applied on stepping motions[END_REF], walking imitation is modeled by estimating the dynamics of a demonstrator from measurements to match the foot location; the driving idea is that the demonstrator and imitator walks use similar controller which input (for the robot pace) is determined from observing the same controller output (from the human). Recently, [START_REF] Ramos | Dynamic locomotion synchronization of bipedal robot and human operator via bilateral feedback teleoperation[END_REF] also used this assumption by analyzing a simplified dynamics of a human in order to make a robot track it.

B. Motion Prediction

The motion of the human has been extensively studied [START_REF] Flash | Models of human movement: Trajectory planning and inverse kinematics studies[END_REF] and allowed to derive from the motion of the humans limbs some invariant properties. One of those properties gives models for hands trajectory generation such as the minimum jerk model. This latter state that the hand motion always tend to minimize the mean square jerk (derivative of acceleration).

In addition to those properties, the motion prediction has also been tackled as a whole body motion generation using the past movement. This kind of prediction can rely on Deep Learning methods using time-sensitive model [START_REF] Lyu | 3d human motion prediction: A survey[END_REF] which are commonly used in robotics to plan a path [START_REF] Kutsuzawa | Sequence-to-sequence model for trajectory planning of nonprehensile manipulation including contact model[END_REF] or adapt it for safety interactions [START_REF] Liu | Deep learningbased human motion prediction considering context awareness for human-robot collaboration in manufacturing[END_REF].

Finally, if the motions to be predicted tend to be more specific and can be pre-recorded offline as sets of trajectories,

Walking Gait Generation

Step phase estimation Minimun Jerk QP

Steps Timings and speed estimation

RLS Filter

Minimum-Jerk Polynome

Past motion parameter estimation Fig. 1: Proposed control scheme for human-humanoid synchronized motion; the superscript * refer to references forwarded to the robot. Motion Data refer to the recorded data detailed in Fig. 4; p c is the center of mass reference for the robot; p swing represents the swing foot trajectory. it is possible to convert each set into a probability distribution and therefore create a so called "motion primitive" [START_REF] Paraschos | Probabilistic movement primitives[END_REF]. This solution allows to anticipate motions that belong to the recorded set over a long period of time [START_REF] Penco | Prescient teleoperation of humanoid robots[END_REF].

C. Contribution

Our control scheme is summarized in Fig. 1; it aims at imitating and anticipating an operator's intention at two levels:

• The upper body motion by anticipating the hand motion of the operator in the cartesian space to achieve synchrony, this motion prediction algorithm requires only the past and present motion data of the user limb; • The locomotion by evaluating and anticipating the operator walking pace and forward it to the robot considering balance constraints. These anticipatory strategies are integrated in a task-space whole-body controller, with experiments using the HRP-4 humanoid robot.

III. ARM MOTION SYNCHRONIZATION

A. Trajectory estimation model

We estimate the hand motion using a polynomial interpolation that match a minimum jerk model and correct it through an RLS (Recursive Least Square) filter.

We can derive from the minimal jerk model an analytic solution [START_REF] Viviani | Minimum-jerk, two-thirds power law, and isochrony: converging approaches to movement planning[END_REF]. It states that to go from an initial point P 0 = (p 0x , p 0y ) at t 0 to a point P 2 = (p 2x , p 2y ) at t 2 passing by an intermediate point P 1 = (p 1x , p 1y ) at t 1 , the trajectory can be expressed as a fifth-order polynomial functions of time:

x(t) = 5 ∑ k=0 c kx t k + p 1x (t -t 1 ) 4 + + p 2x (t -t 1 ) 5 + (t -t 1 ) + = t -t 1 if t ≥ t 1 0 if t < t 1 (1)
And similarly for y(t). [START_REF] Fligge | Minimum jerk for human catching movements in 3D[END_REF] extended, by analyzing human catching movements, the minimum Jerk model to 3D. Therefore we propose to apply the model formulated in [START_REF] Viviani | Minimum-jerk, two-thirds power law, and isochrony: converging approaches to movement planning[END_REF] to 3D. Using this assumption, the hand motion can be interpolated between P 0 and P 2 through an interval [t 0 ;t 2 ] with a waypoint P 1 at t 1 through the parameters C = (c 0 , . . . , c 5 , p 1 , p 2 ) T for each axis. Therefore, we aim to identify the parameters (C x ,C y ,C z ) which define the polynomial regressing the trajectory in the interval [t 0 ;t 2 ].

We evaluate those parameters by using the past hand trajectory considered to be in the interval [t 0 ;t 1 ]. For each axis, C can be determined on the basis of the position, velocity, and acceleration at P 0 and P 2 , and the position and velocity at P 1 . Therefore, identifying C is equivalent to identifying ψ = (x 0 , ẋ0 , ẍ0 , x 1 , ẋ1 , x 2 , ẋ2 , ẍ2 ) T . C can then be computed using a 8 × 8 matrix A such that ψ = AC.

We can compute ψ by formulating an optimization problem that minimizes, within [t 0 ;t 1 ], the norm between:

• The polynomial and the recorded position;

• The second derivative of the polynomial and the recorded acceleration. The recorded trajectories are obtained from a motion capture system composed of IMUs which records angular velocity and linear acceleration that are exploited to reconstruct a user's posture. Motion data are provided at fixed frequency f m . We record b pos as the past position and b acc as the past acceleration of the operator's hand of a duration t 1t 0 . Setting t 0 = 0 and t 2 > t 1 , we formulate the problem as:

ψ = argmin ψ b pos -H pos ψ 2 2 + b acc -H acc ψ 2 2 ψ min ≤ ψ ≤ ψ max (2)
with : H pos = T 0 A -1 and H acc = T 1 A -1 and

T 0 =      1 0 0 0 0 0 0 0 1 1/ f m 1/ f 2 m 1/ f 3 m 1/ f 4 m 1/ f 5 m 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 1 t 1 t 2 1 t 3 1 t 4 1 t 5 1 0 0      T 1 =      0 0 2 0 0 0 0 0 0 0 2 6/ f m 12/ f 2 m 20/ f 3 m 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 2 6t 1 12t 2 1 20t 3 1 0 0     
ψ min and ψ max are set as safety margin. Symbol ψ represent an initial guess of the predicted trajectory, using a RLS filter, we estimate a parameter θ ∈ R 64 to guess a more accurate value ψ following:

ψ = Ψθ + ω (3)
where: Ψ is a 8 × 64 diagonal block matrix of the vector ψT and ω is a zero-mean noise with covariance R. The parameter θ can be estimated online with the past data; it is updated using the real value ψ with the following two equations:

P -= P + R; K = P -Ψ T (ΨP -Ψ T + R). (4)
We then update P, R and θ as follows:

P ← (I -KΨ)P -; θ ← θ + K(ψ -ψ); R ← (λ -1 -1)P;
with λ ∈ [0; 1] often refereed as a "forgetting factor" and provides an approximated exponential decaying weighting on the past data [START_REF] Simon | Optimal State Estimation: Kalman, H-Infinity, and Nonlinear Approaches[END_REF]. P defines the estimation covariance and K the estimator gain matrix.

We train this filter on past known data and use ψ to obtain a motion prediction of the user hand. Then, we can compute from ψ any set of position acceleration and velocity in the interval [t 1 ;t 2 ]. However, we rather use (x 2 , ẋ2 , ẍ2 ) as they represent the furthest prediction available for a value of ∆t.

Finally, we can forward the estimated (x 2 , ẋ2 , ẍ2 ) to an inverse kinematics control scheme using an acceleration based tracking law:

pd = -K p (p -p r ) -K d ( ṗ -ṗr ) + pr , (5)
where p r represent the reference end-effector pose and pd the robot desired end-effector acceleration. The gain K d is set to the critical value K d = 2 K p , therefore, K p reflects the tracking behavior w.r.t the input target.

B. Determining the tracking latency to compensate

The main external parameters of the presented model are the timings t 1 and t 2 (having set t 0 = 0); t 1 is the duration of the input trajectory sequence and t 2t 1 is how far ahead in time we estimate (x 2 , ẋ2 , ẍ2 ).

Therefore, the optimal value t 2t 1 = ∆t opt should compensate two sorts of delays:

1) The first delay is inherent to communication lag between the instant the user performs a motion and the instant the measured position is transmitted to the robot. This delay can be estimated using communication average delays; 2) The second delay comes from the robot control itself.

The robot motion has an unfixed latency to reach the target depending on the control parameters. We can measure the delay between a user motion and the robot using the cross correlation between two discrete signals f (t) and g(t) using:

f * g(n) = ∞ ∑ m=-∞ f (m)g(m + n)
A value of ∆ opt can then be estimated beforehand using pre-recorded values of the robot hand velocity magnitude v r and of the user v h , and applying:

∆t opt = arg max t∈R v h * v r (t).
However, it is important to keep in mind that the prediction error increases with higher values of t 2t 1 , so ∆t cannot always be set at will.

To illustrate this, the proposed model is evaluated offline using the recorded motion of a user's hand. Figure 3 a prediction sample in one axis with two values of ∆t = t 2t 1 , and Fig. 2 shows the average error with the range 2σ depending of ∆t. The limit of the model are presented from the choice of ∆t. We note in Fig. 2 that the predicted trajectory has an increasing average error during uniform motion in the interval [1.0, 3.0] as ∆t increases. Moreover, notice that stroke motions induce a high prediction error due to two reasons:

• A stroke motion is not 'grasped' by the minimum jerk model; • As the estimation is performed from past data, a stroke motion cannot be anticipated.

The limit of the prediction of the stroke motions can be seen in Fig. 3a where the estimated pose during the stroke motion is seen 0.1 seconds ahead instead of 0.2 seconds. The lack of precision can also be evaluated by monitoring the filter estimation covariance matrix P; ∆t is chosen to account for the control delay and the limit of the model itself.

IV. WALKING IMITATION

Whereas arm motion synchrony can be achieved by imitating the kinematics of the human limbs, walking synchrony requires to meet balance constraints that are paramount for both humans and humanoids. Besides, one must estimate relevant human walking parameters to be conveyed to the humanoid robot in order to achieve a synchronized walking. There are of two kinds:

• The kinematics that represent the support foot and the future step location; • The timings which indicate the take off and landing time of the swing foot. We believe this is the most critical feature to reach a high level of embodiment [START_REF] Aymerich-Franch | The second me: Seeing the real body during humanoid robot embodiment produces an illusion of bi-location[END_REF]. It is important that human's and humanoid's legs take-off and land synchronously [START_REF] Benallegue | A computational model for synchronous motion imitation by humans: The mirror controller applied on stepping motions[END_REF].

In the case of teleoperated humanoid locomotion, it is very likely that the operator will have a limited space to control the humanoid robot; it is important to consider also the case where the operator uses on-site stepping.

A. Walking frequency estimation

Synchronizing the user's steps with the humanoid means providing to the robot the timing of foot lift-off and landing ahead of time. We also need to tackle both cases where the user is stepping in place and actually walking. For that, we propose to identify those stepping parameters using a Deep Learning model in order to extract from the user's motion data the stepping phase s(t) defined as:

s(t) = t -t st T s -t st ( 6
)
where t st is the double support starting time and T s is the single support ending time, we define (t st , T s , s ds ) as the step phase parameters.

As it can be unclear to define t st (especially when the user's motion starts), we set t st such that s reaches a value s ds when the user starts its single support phase, hence:

t st = s ds T s -t ds s ds -1 (7)
In fact, s ds can be set in order to match the human double support period during one step which is s ds = 1 φ 3 ≈ 0.24, where φ is the golden ratio [START_REF] Iosa | The golden ratio of gait harmony: repetitive proportions of repetitive gait phases[END_REF]. 1) Dataset: We use as a dataset the motion data of a user walking recorded with the Neuron motion capture mockup. We record the angular velocity of the lower body parts and the velocity of the upper body parts. Figure 5 maps the sensor location and the recorded data where the displayed frame is attached to one of the user limb (in our case the hip center). We do not apply an additional filter to the built-in software. The proprietary software can provide the contact state of each foot with the floor using the acceleration data that we use to compute the user step phase parameters.

We split the data in sequences containing the user's motion data as inputs and sequence of the future user step phase as expected output. The sequences on which s(t) transitions from a constant 0 value to non-zero indicate the transition from a standing phase to a walking motion in a close future. However this sequence cannot be predicted and therefore must be discarded. We also include motion data where the user is static or is doing random motions with the upper body, e.g., rotating, bending forward... where s(t) is always equal to 0. This is to reduce, as much as possible, false positive estimations of steps. Figure 4 shows an example of a motion data and the output step phase.

2) Deep learning model using Recurrent Neural Networks: When dealing with temporal data, our model is expected to see each time values of a sequence as a consequence of the previous ones. Recurrent Neural Networks (RNN) models proved high performance in various problems such as language translation [START_REF] Cho | Learning phrase representations using RNN encoderdecoder for statistical machine translation[END_REF]. A RNN cell process each elements of a sequence one by one and use the previous computation of the sequence input to get the next sequence value. Therefore, if we set as an input a temporal sequence (x 1 , • • • , x n ) T , then a RNN layer will compute a sequence (h 1 , . . . , h n ) T and an output sequence (y 1 , . . . , y n ) T as follows:

h i = σ (W h x i +Uh i-1 ) (8) y i = W y h i + b y (9)
(h 1 , ..., h n ) T is called the hidden states, W h ,W y , b y ,U are the RNN weights and σ is a smooth, bounded function such as a hyperbolic tangent or a logistic sigmoid function.

The concept of RNN has been extended to more complex models to handle the long term dependencies in the input sequence. Long Short Term Memory (LSTM) or Gated Recurrent Unit (GRU) are two RNN models that belongs to the group of Gated Recurrent Neural Network and are widely used in the deep learning community. Their differences lies in the method to control the long term past data in the input sequence. [START_REF] Chung | Empirical evaluation of gated recurrent neural networks on sequence modeling[END_REF] showed some empirical equivalence between both models, however, as GRU requires less parameters to be trained. Therefore, we use GRU for our training model.

A GRU [START_REF] Cho | On the properties of neural machine translation: Encoder-decoder approaches[END_REF] updates the hidden state value h i using the previous hidden state h i-1 and a candidate hidden state hi such as:

h i = (1 -z i )h i-1 + z i hi (10)
with the "update gate" z i and the candidate hidden state hi :

z i = σ (W z x i +U z h i-1 ) hi = tanh(W h x i +U(r i h i-1 )) (11)
with the element-wise multiplication and r i the reset gate:

r i = σ (W r x i +U r h i-1 ) (12)
Our model tends to extract from a motion data sequence the future step phase sequence of the user. It can therefore be set as a Seq2Seq (Sequence to Sequence architecture). [START_REF] Cho | Learning phrase representations using RNN encoderdecoder for statistical machine translation[END_REF], [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF] presented Seq2Seq which is divided into two parts: the first, Encoder, is a RNN layer that converts the input sequence into an internal representation of the motion, a context, (which will be the last element of the hidden state). The last part is the Decoder that converts the computed context into the expected sequence using a second RNN layer. Then, we apply Eq. 9 on the output hidden state to match the sequence output shape.

This decoder is set as the next input to the previous decoder output. In the model proposed in [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF], the context is set as the first decoder hidden value h 1 . This model is illustrated in Fig. 6. 3) Learning: We train our model using the Seq2Seq model [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF], the loss function is the mean absolute error between the prediction and the model output which is optimized by the Adam (Adaptive Moment Estimation) algorithm. The dataset characteristics are shown in Table I. The samples in Fig. 7 represent the average estimation of the stepping phase on validation data. We obtain this average by computing the predicted step phase sequence at the model frequency. We show that the model is able to anticipate the user's single support phase but it represents the prediction with the biggest variance, whereas the single support phase appears to be properly identified. It appears necessary, especially in double support, to use multiple estimations of the model to provide a more accurate value of the stepping parameters.

Input

4) Prediction: Our model outputs a prediction sequence of s = (s 0 , • • • , s n s -1 ) T supposed to be linear if we do not consider sequences corresponding to an end and start of new steps. Therefore, for each prediction we can extract a predicted slope and offset x = (a, b) T using a linear regression such that:

x = argmin x s -T x 2 T = 0 1/ f md • • • (n s -1)/ f md 1 1 • • • 1 T (13)
And then compute the step phase parameters

t st = -b/a (14)
T s = 1 -b a -t st (15)
Finally, using multiple computations of the step phase parameters we can forward an average estimation to the robot and use the standard deviation of this latter as a safety metric for false positive step detection or unreliable estimations.

B. Steps location

When the user is performing a walk-in-place mode, we can estimate an equivalent forward walking speed that will match the user step frequency. The relation between reference velocity and the steps frequency is set from empirical observation of the human walking which lead to Dean's relation [START_REF] Dean | An analysis of the energy expenditure in level and grade walking[END_REF] linking the human's step frequency f s , speed magnitude |v| and height h as follows: [START_REF] Paraschos | Probabilistic movement primitives[END_REF] with h in cm and |v| in cm.s -1 . We can therefore use the estimated step frequency to generate steps location that will match the computed speed under the humanoid robot kinematics constraints.

|v| = f s 0.157 × h 172 2

V. EXPERIMENTS AND RESULTS

We assessed whole-body synchronized walking using the HRP-4 humanoid from Kawada Robotics. The human operator is equipped with the Neuron Perception 2 IMU wholebody tracking system. The precision of Neuron has been evaluated w.r.t the VICON motion capture system in [START_REF] Choo | Validation of the perception neuron system for full-body motion capture[END_REF]. Performances shows a RMSE below 4 degrees for mosts joint angles that are in the same range of a the wearable multi-modal system in [START_REF] Fang | Robotic teleoperation systems using a wearable multimodal fusion device[END_REF].

We processed the experiment as follows: (i) we use the Neuron tracking system to collect data from which human walking features prediction is learned; (ii) we integrate upper-body and walking synchronizations using the taskspace mc rtc quadratic programming framework [START_REF] Bouyarmane | On weight-prioritized multi-task control of humanoid robots[END_REF], [START_REF] Bouyarmane | Quadratic programming for multirobot and task-space force control[END_REF], that is used also to control the robot; (iii) we assess the walking and terminal point synchronization on simulation using Chorenoid; finally (iv) mc rtc switches into realtime whole body walking with motions of hands and head. We show a human-humanoid side-by-side synchrony and a human-humanoid face-to-face mirror synchrony motions. The computation of the commanded angles and the logged data are done at a frequency of 200 Hz.

A. Experimental constraints

Unlike simulations where torque and joint speed limits can be set at will, the real HRP-4 robot has limited capabilities; e.g., velocity limits and range of motions, lower degrees of freedom, less graceful motions for the available degrees of freedom, and lower range of motion for the existing DoFs as compared to their human counterpart. This is the reason (also considering safety) why the following restrictions have to be considered in real humanoid experiments:

• the learning process applies to a specific human, because of the anthropomorphic and gait differences, the learning process does not generalize to any human operator. Expanding the learning to a large range of human operator is possible, but requires a huge amount of data which is not the core part of our work; • in order for the synchrony to occur, the operator is asked to perform slow motions to match the speed capacity of the robot. This also allows to achieve experiments without tethering the robot; • the robot is powered through an outside connected cable and controlled from a wired external computer, therefore, the range of motion of the robot is limited and the reachable space are limited; • an additional game-joystick is used by another operator for safety purposes. It allows to trigger on and off the terminal points for synchronized tracking on demand. Fig. 8: Robots and operator's hands linear velocity magnitude and the computed delay using cross-correlation.

B. User hand control

We apply the model developed in Sec. III such as the robot's hands track the operator's predicted hands' Cartesian position and velocity and the current hands orientation and rotational speed, that are updated at a rate of 60 Hz.

However, task space retargeting leads to dissimilarity between the human and the robot posture because of the human/robot shape and degrees of freedom mismatch. Moreover, if the robots arm has kinematics redundancy, the arm posture, tracking only a 6d reference could differ from the user's. Hence, to achieve a convincing synchrony and motion between the human and the robot arm motion, we must ensure that the robot's links between the shoulder and the hand are also following an orientation similar to the user.

Therefore, each of the robot's two upper limbs are controlled using two tasks. The first is a position task on the robot hand. The second is a task that tracks the operator's upper arm orientation, this has the benefit to provide a similar posture between the operator at the cost of the precision on the hand tracking.

The accuracy of the predicted hand position is shown in Figs. 2 and3, therefore to assess the motion prediction performance in terms of synchrony, we compare the operator and the robot linear velocity magnitude. More specifically, as the robot's hand motion is also constrained by the arm orientation task, the velocity between the user and the robot will not be identical, we focus on the instant when user's and robot's hand speed reach a velocity close to zero. We show in Fig. 8 the velocity magnitude of the user and the robot during the experiment with and without using the predicted hand position and velocity. Note that anticipating the motion of the operator and forwarding it to the controller enables the robot motion to be synchronized to that of the operator w.r.t the instant when the user's hand comes to a stop. Finally, Fig. 8 shows that the measured delay is reduced from 58 ms to 3 ms (a typical robot control loop is 5 ms) which is under our data logging frequency. However, to obtain this behavior, the model has been set for a value ∆t = t 2t 1 = 100 ms which is almost twice the measured anticipation. Nevertheless, this is not a problem as ∆t is a tunable parameter that can be empirically set to minimize the delay.

C. Walking parameters online estimation

The robot is using a walking pattern generator that provides a center of mass reference trajectory (position, velocity, acceleration), see [START_REF] Scianca | MPC for humanoid gait generation: Stability and feasibility[END_REF], and that requires timings parameters and a reference walking velocity.

Therefore we must provide: (i) the single and double support duration; (ii) the walking direction; (iii) the user standing foot. The two first inputs come from the steps timing estimator described in Sec. IV using the same frequencies described in Table I. During the experiments we provide a double support duration time only once; whereas the single support duration is updated regularly during the swing foot phase. On the other hand to estimate the walking direction, one solution could be to use the measured velocity from the IMU on the user's hip; however, if this data can be used for forward walking it cannot be reliable when it comes to side walking as we can't distinguish sway from locomotion. Therefore we rather rely on the velocity measured at the user's swing foot, the drawback of this estimation is that it is available only during the single support phase, therefore, we notice a kinematic discrepancy between the robot and the human swing phase as the robot will swing its leg once the walking direction will be estimated (see, Fig. 9). To estimate the user's support foot, we rely on the sign of the user's hip velocity in the y-axis (cf. Fig. 5) during the double support phase. We assume that the user does not use the same support twice in a row, therefore during a walking sequence, we estimate the first support foot only once when the walking sequence starts. Figure 10 gathers the contact state of the operator (also obtained from the mocap software) and the robot with the current foot height of both entities. We show that the robot is able to follow the user's pace and swing foot, we still notice some errors in the phase estimation which can be due to a bad estimation from the model itself and from the implemented safeties that can prevent the robot to change its step-timing when it is close from landing w.r.t the previous input step timing.

VI. CONCLUSION

We proposed a novel control framework that enables a humanoid robot to synchronize its steps and arms motions with those of a human in real-time. Using a whole-body motion capture system base on IMUs, we used an invariant of the human hand motion, the minimum jerk model, coupled with a RLS filter to anticipate the hand's cartesian motion. Using recorded data of human walking, we trained a deep learning model capable of predicting the human duration of double and single support phase. Those predicted data are integrated in a task-space whole-body controller with constraints handling to synchronize the robot with the operator. Our approach is successfully assessed in real experiments with the HRP-4 humanoid robot synchronizing motions with an operator.

The training of the model was performed offline and therefore dependent of the past training data, it would be interesting to switch to a self-learning model approach which might also require an improvement in the detection of the contact state. Another major improvement could be to use a brain computer interface to anticipate on the motion instead of relying solely on the body motions [START_REF] Gergondet | Using brain-computer interface to steer a humanoid robot[END_REF] (EEG detects intentional motion about 300 ms ahead of time).

We used the minimum jerk model to anticipate Cartesian translation of the hands of an operator. However a similar method cannot be applied on the operator's hand orientation since (i) to our knowledge there is no equivalent model that study the human hands orientation; (ii) our prediction model relies on the acceleration data of the user's hand, our device fail in estimate properly the hand angular acceleration. Exploring the behavior of the hand's and head orientation in future work would open perspectives in embodiment [START_REF] Aymerich-Franch | The second me: Seeing the real body during humanoid robot embodiment produces an illusion of bi-location[END_REF].

Fig. 2 :

 2 Fig.2: Average error between the estimated x 2 and the real state depending of ∆t for the motion in Fig.3for t ∈ [1; 4] using t 1t 0 = ∆t, f m = 120 Hz.

Fig. 3 :

 3 Fig. 3: Evaluation of the prediction model in one axis depending on ∆t = t 2t 1 . Using t 1t 0 = ∆t, f m = 120 Hz.

Fig. 4 :

 4 Fig. 4: Example of step phase profile correlated with the step height.

Fig. 5 :

 5 Fig. 5: Experimental set-up: operator walking with the robot (left) and location of sensor and measured data type (right).

Fig. 6 :

 6 Fig.6: Seq2Seq model in[START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF].

Fig. 7 :

 7 Fig.7: Phase estimation with validation data compared to the ground truth value, the grey area represent the confidence interval of 95%.

  Using measured hand position and velocity, measured delay: 58 ms.

  Using predicted hand position and velocity. t 2t 1 = 0.1s, t 1t 0 = 0.1s, f m = 60 Hz, measured delay: 3 ms.

Fig. 9 :

 9 Fig. 9: Delay in the walking reference velocity. The user has already placed his feet backward meanwhile the robot is vertically lifting the feet.

Fig. 10 :

 10 Fig. 10: Support state of the user and the robot with the vertical position of the swing foot; the estimations errors can be seen at the instants when the support state change.

  shows

		0.0175	Average Prediction error
		0.0150
		0.0125
	Error (m)	0.0075 0.0100
		0.0050
		0.0025
		0.0000
			0.06	0.08	0.10	0.12 Δt	0.14	0.16	0.18	0.20

TABLE I :

 I Model parameters

				Data stream frequency f d	60 Hz
				Model frequency f md		30 Hz
				Input sequence length		30
				Output sequence length n s	10
				Number of hidden layers	128
				Motion data component	8
				TABLE II: Training results
				Trained samples		7811
					Batch size		125
				Number of epochs		50
				Training data loss		0.0023
			Training data mean absolute error	0.0206
				Validation data loss		0.0052
			Validation data mean absolute error 0.0385
					sds
		1.0			Validation step phase Estimation average Left Foot	0.16
					Right Foot		0.14
		0.8						0.12
	Phase	0.4 0.6						0.06 0.08 0.10	Height (m)
								0.04
		0.2						0.02
		0.0						0.00
		0	1	2	3	4	5	6
					Time (s)

https://www.youtube.com/watch?v=InYsET-v7wM

https://www.youtube.com/watch?v=7atZfX85nd4

https://www.youtube.com/watch?v=G9p9jdmJQOQ

https://neuronmocap.com/

https://www.xsens.com/motion-capture


Clément Cosserat

Ben Gabrielson

Emilie Chouzenoux

Jean-Christophe Pesquet

Tülay Adali

A PROXIMAL APPROACH TO IVA-G WITH CONVERGENCE GUARANTEES

Keywords: blind source separation, IVA-G, proximal algorithms, joint ISI, non-convex optimisation

Independent vector analysis (IVA) generalizes independent component analysis (ICA) to multiple datasets, and when used with a multivariate Gaussian model (IVA-G), provides a powerful tool for joint analysis of multiple datasets in an array of applications. While IVA-G enjoys uniqueness guarantees, the current solution to the problem exhibits significant variability across runs necessitating the use of a scheme for selecting the most consistent one, which is costly. In this paper, we present a penalized maximum-likelihood framework for the problem, which enables us to derive a non-convex cost function that depends on the precision matrices of the source component vectors, the main mechanism by which IVA-G leverages correlation across the datasets. By adding a quadratic regularization, a block-coordinate proximal algorithm is shown to offer a suitable solution to this minimization problem. The proposed method also provides convergence guarantees that are lacking in other state-ofthe-art approaches to the problem. This also allows us to obtain overall slightly better performance, and in particular, we show that our method yields better estimation in average than the current IVA-G algorithm for various source numbers, datasets, and degrees of correlation across the data.

INTRODUCTION

Joint blind source separation (JBSS) enables recovery of information across multiple datasets through latent variables, and has proven powerful for extracting useful features from related datasets in a large number of application areas. In particular, by enabling multiple datasets to fully interact, it leverages the joint information across them and has been attractive for data fusion, [START_REF] Lahat | Multimodal data fusion: A methodological overview methods, challenges and prospectives[END_REF][START_REF] Adalı | Diversity in independent component and vector analyses: Identifiability, algorithms, and applications in medical imaging[END_REF]. Independent vector analysis (IVA) [START_REF] Adalı | Diversity in independent component and vector analyses: Identifiability, algorithms, and applications in medical imaging[END_REF][START_REF] Kim | Independent vector analysis: Definition and algorithms[END_REF] generalizes independent component analysis (ICA) to multiple datasets. It captures the joint information across the datasets through selection of an appropriate multivariate probability density model, for the source component vectors (SCVs), the key element of the IVA formulation. When a multivariate Gaussian density model is used, it is termed IVA-G, and the method enables effective fusion of multiple datasets using only second-order statistics [START_REF] Adalı | Diversity in independent component and vector analyses: Identifiability, algorithms, and applications in medical imaging[END_REF][START_REF] Anderson | Joint blind source separation with multivariate Gaussian model: Algorithms and performance analysis[END_REF], and also generalizes multiset canonical correlation analysis (MCCA), which has been widely used in an array of applications [START_REF] Kettenring | Canonical analysis of several sets of variables[END_REF]. In addition, IVA provides identifiability under rather general conditions [START_REF] Anderson | Joint blind source separation with multivariate Gaussian model: Algorithms and performance analysis[END_REF][START_REF] Via | Joint blind source separation from second-order statistics: Necessary and sufficient identifiability conditions[END_REF]. However, one of its main difficulties is that the cost function is non convex, which makes it necessary to use of a criterion to evaluate the reproducibility (stability) of the solutions and pick one as the best run among many obtained with random initializations [START_REF] Adali | Reproducibility in matrix and tensor decompositions: Focus on model match, interpretability, and uniqueness[END_REF]. The variability in the solution space becomes a bigger issue when operating in the neighborhood of the local stability conditions, increasing the variability in the solution space. Hence, a solution with convergence guarantees is highly desirable.

In this paper, we propose and analyze a new estimator for performing IVA-G by using a penalized maximum likelihood approach. A proximal alternating minimization algorithm is used to minimize the resulting cost function, which alternates over updates on the covariance matrices of the SCVs and on the mixing coefficients. Convergence guarantees are provided based on recent tools from non smooth analysis. As we demonstrate through simulations in a number of relevant scenarios, the new algorithm exhibits a greater level of accuracy than IVA-G in the estimation of the ground truth sources.

The paper is organized as follows. In Section 2, we formulate the problem and describe the data model and the parameters involved. We derive a cost function starting from the maximization of a likelihood criterion. Section 3 introduces our proximal alternating algorithm, details its updates, and discusses its convergence properties. Finally, we describe our experimental protocol, and discuss the results we obtained in Section 4, before concluding in Section 5.

JBSS PROBLEM FORMULATION

We first formulate the problem of JBSS, and consider the case when the latent sources are modeled as random vectors. In this model, a signal is a sequence of T ∈ N realisations of its source. There are K observed datasets (k ∈ {1, . . . , K}). Within each dataset, there exist N observed signals (n ∈ {1, . . . , N }). We denote the N observed signals for a sample t ∈ {1, . . . , T } within the k-th dataset by the vector

x [k] (t) = [x [k] 1 (t), . . . , x [k] N (t)] ⊤ ∈ R N (
here, (.) ⊤ denotes the transpose). These signals are obtained by the linear mixing of N latent source signals s [k]

(t) = [s [k] 1 (t), . . . , s [k] N (t)] ⊤ ∈ R N ,
which we aim to jointly estimate.

The linear mixing of sources within each dataset is represented by unknown invertible matrix

A [k] = a [k] i,j 1≤i,j≤N ∈ R N ×N .
The generative model is given as

x [k] (t) = A [k] s [k] (t), k ∈ {1, . . . , K}. (1)
Since the K datasets are observed over T samples of the data, the datasets are themselves represented by matrices

X [k] = [x [k] 1 , . . . , x [k] N ] ⊤ ∈ R N ×T
where, for every k ∈ {1, . . . , K} and n ∈ {1, . . . , N }, x

[k] n = [x [k] n (1), . . . , x [k] n (T )] ⊤ ∈ R T .
Using a similar notation, the corresponding source signals are given by S

[k] = [s [k] 1 , . . . , s [k] N ] ⊤ ∈ R N ×T .
The goal of JBSS is to estimate the source signals within each dataset, which is done by estimating demixing matrices

W [k] ∈ R N ×N that demix each dataset X [k] into the source sig- nals estimates Y [k] = W [k] X [k] = W [k] A [k] S [k] , with Y [k] = [y [k] 1 , . . . , y [k] N ] ⊤ ∈ R N ×T .
Source signals are perfectly estimated when the mixing-demixing matrices W [k] A [k] can be expressed as permuted diagonal matrices, indicating perfect estimation subject to permutation and scale ambiguities. The nth row of W [k] is given by the vector w

[k] n ∈ R N ,
which is used to obtain the nth estimated source signal in the kth dataset, given by y

[k]⊤ n = w [k]⊤ n X [k] , with y [k] n ∈ R T .
JBSS exploits dependence across the datasets by modeling dependence through the sources. Across the K datasets, sources of the same index n are modeled as dependent to other sources of index n, thus forming N sets of K dependent sources. In IVA terminology, each of these sets is referred to as a "source component vector" (SCV). The nth SCV is given by sn = [s [START_REF] Lahat | Multimodal data fusion: A methodological overview methods, challenges and prospectives[END_REF] n , . . . , s

[K] n ] ⊤1
, which is a K-variate random vector. The observed signals are thus realisations of the random vector xn = [x

[1] n , . . . , x [K] n ] ⊤ where ∀k ∈ {1, . . . , K} : x [k] n = N m=1 a [k] n,ms [k] m . Sources sn are es- timated by yn = [y [1] n , . . . , y [K]
n ] ⊤ which are similarly modeled by the mixing of the xn by the demixing matrices. Over T samples of data (i.e., T realizations), the nth SCV sn is represented through the matrix Sn = [s [START_REF] Lahat | Multimodal data fusion: A methodological overview methods, challenges and prospectives[END_REF] n , . . . , s

[K] n ] ⊤ ∈ R K×T , and is estimated from Yn = [y [1] n , . . . , y [K] n ] ⊤ ∈ R K×T .
Each SCV is modeled as independent from all other SCVs, thus any two sources across the datasets are modeled as dependent only if they correspond to the same source index n (nth SCV).

For each of the SCVs sn, n ∈ {1, . . . , N }, we denote covariance Σn = E sn s ⊤ n ∈ R K×K as the covariance of sn, and Cn = Σ -1 n as the corresponding precision matrix of sn. We similarly define Σn ∈ R K×K as the covariance of yn, and Ĉn = Σ-1 n as the precision matrix of yn.

IVA with multivariate Gaussian model

IVA is a method of JBSS that estimates W [k] , for k ∈ {1, . . . , K}, such that the corresponding SCV estimates are maximally independent. IVA does this by minimizing the mutual information across the SCVs, which is dependent on the probability distribution function (PDF) assumed on each of the sources. Typical implementations model each SCV by a single multivariate distribution, allowing mutual information of the SCVs I{yn} to be written with respect to the marginal entropy of each individual source H{y n } -H{yn}.

IVA assuming a multivariate Gaussian PDF for each SCV (IVA-G) [START_REF] Anderson | Joint Blind Source Separation With Multivariate Gaussian Model: Algorithms and Performance Analysis[END_REF] can be considered the simplest of the IVA algorithms, in that it exploits only correlation between sources as the measure of dependence. By denoting W as the collection of matrices W [k] for k ∈ {1, . . . , K}, the IVA-G cost function is given by:

JIVA-G (W) = N K log(2πe) 2 + 1 2 N n=1 log det Σn - K k=1 log det W [k] . ( 2
)
1 throughout this paper, bold characters represent random variables, and upper cases letters represent matrices To derive our own variant of IVA-G for our proximal approach, we revisit the general cost for maximum likelihood IVA where now we additionally estimate the corresponding parameters that parametrize the PDFs of the SCVs. For the Gaussian case, these parameters are the covariances Σn, or alternatively represented by the corresponding precision matrices Ĉn. It is straightforward to show that by applying a change of variable to the general IVA cost, assuming a multivariate Gaussian PDF for each SCV, then applying independence across the SCVs and across the T samples, the IVA-G cost function can alternatively be expressed by the following up to additive constant:

JIVA-G (W, C) = - K k=1 log det W [k] + N n=1 1 2 tr ĈnWn ΛW ⊤ n -log det Ĉn . (3)
where we use the following notation:

• C is the collection of all Ĉn, for n ∈ {1, . . . , N }.

• Wn ∈ R K×KN is a matrix corresponding to the nth SCV, defined as

Wn = W [1] n , . . . , W [K] n
, a block matrix obtained by horizontal concatenation of the

W [k] n ∈ R K×N where W [k] n has w [k]
n as its k-th line and zeros elsewhere.

• Λ ∈ R KN ×KN is the empirical covariance matrix of vector vec(x) ≜ [(x [1] ) ⊤ , . . . , (x

[K] ) ⊤ ] ⊤ ∈ R N K , where x [k] ≜ [(x [k] 1 ) ⊤ , . . . , (x [k]
N ) ⊤ ] ⊤ is the random vector whose realisations produce the columns of the dataset X [k] and vec(x) is the vertical concatenation of these random vectors.

This yields an alternative but equivalent cost to [START_REF] Adalı | Diversity in independent component and vector analyses: Identifiability, algorithms, and applications in medical imaging[END_REF] where we also measure cost as a function of Ĉn. Here we assume that log|det( Ĉn) | = ∞ if Ĉn is not positive definite.

Both the cost functions in (2) and ( 3) are non-convex, which makes minimization challenging. In [START_REF] Anderson | Joint Blind Source Separation With Multivariate Gaussian Model: Algorithms and Performance Analysis[END_REF], the authors minimize (2) by proposing a block-Newton descent method for iteratively updating W [k] . Despite this method's practical performance, it does not benefit from proven convergence guarantees, and may suffer from instabilities dependent on the data. We thus introduce a regularized version of the problem by instead optimizing over (3), and build a structured optimization method to tackle the problem with the explicit goal of proven convergence.

Taking into account ambiguities

As noted in Section 2, there exists scaling ambiguity that is unresolved by demixing matrices W [k] , that is, for any α

[k] n ∈ R \ {0}, simultaneously multiplying s [k] n and w [k] n by α [k]
n does not change the cost. This scaling ambiguity ultimately results in an unbounded set of minima of the cost function [START_REF] Kim | Independent vector analysis: Definition and algorithms[END_REF]. This may introduce some misbehaviour of minimization algorithms as the norm of the estimated W [k] may tend to zero or infinity. A way to reduce this to only a sign ambiguity thus consists of constraining the diagonal components of Ĉn to be equal to a given value (e.g. 1). To this end, we add a quadratic regularization term in our cost function, thus leading to a regularized version of (3):

JIVA-G (W, C) = - K k=1 log det W [k] + N n=1 1 2 tr ĈnWn ΛW ⊤ n -log det Ĉn + K k=1 α Ĉn(k, k) -1 2 (4)
where we denote Ĉn(k, k) as the kth diagonal element of Ĉn, and α > 0 is a regularization hyperparameter.

OPTIMIZATION ALGORITHM

We now present our optimization scheme to minimize (4). Our method relies on the PALM algorithm (Proximal Alternating Linearized Minimization) [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. We first present the general form of PALM, and then specify it applied to (4).

General form

The PALM algorithm aims at finding a critical point of a function of the form:

(x1, . . . , xp) → Ψ (x1, . . . , xp) = H(x1, . . . , xp) + p k=1 f k (x k )
where H is a smooth function defined on R d 1 ×...×dp and each f k is a function defined on R d k . An advantage of PALM is that none of the involved functions need to be convex. Given an initialization (x 0 1 , . . . , x 0 p ) ∈ R d 1 ×...×dp , the PALM algorithm iterates, for iteration index i, over the blocks of variables x k for k ∈ {1, . . . , p}. Specifically, it generates x i = (x i 1 , . . . , x i p ) such that, for k ∈ {1, . . . , p},

x i+1 k ∈ prox f k c i k x i k - 1 c i k ∇ k H(x i+1 1 , . . . , x i+1 k-1 , x i k , . . . , x i p )
where

c i k = γ k L i k , L i k is a Lipschitz modulus of the gradient func- tion x ∈ R d k → ∇ k H(x i+1 1 , . . . , x i+1 k-1 , x, x i k+1 , . . . , x i p ),
and γ k > 1. Under some mild technical assumptions [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], if Ψ is regular enough, the sequence generated by PALM has two interesting properties: the sequence (x i ) i≥1 converges toward a critical point of Ψ and the sequence of their cost decreases. Note that variable metrics variants of the PALM algorithms have been proposed in [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF] and shown to be effective to solve various inverse problems [START_REF] Repetti | Euclid in a taxicab: Sparse blind deconvolution with smoothed l1/l2 regularization[END_REF][START_REF] Abboud | An alternating proximal approach for blind video deconvolution[END_REF].

Application to IVA-G

PALM can be applied to the minimization of (4), by defining the functions within (4) operating on W and C:

H(W, C) = N n=1 1 2 tr ĈnWn ΛW ⊤ n + K k=1 α Ĉn(k, k) -1 2 , (5)
f (W) = - K k=1 log det W [k] , (6)
g(C) = - 1 2 N n=1 log det Ĉn . (7)
Let us discuss the practical implementation of PALM in our context. First, we can express the gradient of ( 5) with respect to W, or ∇ W H(W, C), by considering the W [k] blocks separately as:

∂H(W, C) ∂W [k] = ∂H(W, C) ∂w [k] n,m 1≤n,m≤N =    [ Ĉ1W1 Λ] k,(k-1)N +1 . . . [ Ĉ1W1 Λ] k,(k-1)N +N . . . . . . [ ĈN WN Λ] k,(k-1)N +1 . . . [ ĈN WN Λ] k,(k-1)N +N    . (8)
Similarly, ∇ C H(W, C) is calculated for each Ĉn block as:

∂H(W, C) ∂ Ĉn = 1 2 Wn ΛW ⊤ n + α (diag( Ĉn) -IK ). (9)
From these expressions, we can derive the Lipschitz moduli

L i W = ∥ Λ∥S max 1≤n≤N ∥ Ĉ i n ∥S and L i C = α
where Ĉ i n is the estimate of Ĉn at iteration i and ∥ • ∥S denotes the spectral norm, that is the largest spectral value of the symmetric definite positive matrix it applies to.

Let us discuss the calculation of the proximal operators of ( 6) and [START_REF] Adali | Reproducibility in matrix and tensor decompositions: Focus on model match, interpretability, and uniqueness[END_REF], respectively. We can separate the blocks into sub-blocks and calculate, for every k and n, the proximal operator of f k (W [k] ) = -log | det W [k] | and gn( Ĉn) = -1 2 log det Ĉn. Both these functions are spectral, i.e., they only depend on the singular values of the matrices. Function gn corresponds to a convex function of a symmetric matrix for which the proximity operator is known [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. Function f k is non-convex, but its proximity operator can be derived very similarly to the previous one by performing a singular value decomposition of W [k] . We thus have all the analytical formulas needed to implement the iterations of PALM algorithm.

Regarding the convergence analysis, we can easily establish the existence of a minimizer. In addition, the above formulas combined with Hadamard and Jensen inequalities allow to prove that (4) satisfies assumptions for the convergence of the PALM algorithm. We can thus guarantee that the generated sequence, assuming it is bounded, will converge to a critical point of the cost function.

EXPERIMENTS

We tested the performance of our algorithm on simulated data over K equal to 2 or 5 datasets, and a number of sources N equal to 3, 5, or 10. For each pair (K, N ), we randomly generated 100 simulations of data defined by N covariance matrices and K mixing matrices. The models were used to simulate SCVs sn, and sample the source datasets S [k] , which were then mixed by A [k] ∈ R N ×N to generate the observed datasets X [k] . Each datasets was composed of T = 10 4 samples. The coefficients of the A [k] were generated according to random realizations of the zero mean unit variance normal distribution. The SCV covariance matrices Σn were generated according to random realizations from the Wishart distribution Wp(V, n), where n = K +10 and V ∈ R K×K is a matrix whose diagonal coefficients are 1 and off-diagonal coefficients are all equal to ρ ∈ [0, 1]. Those forms for V in Wishart distribution allow us to control the typical correlation within the SCVs across the datasets, which is high if ρ is close to 1 and low if ρ is close to 0.

For each simulation of the data, we compared the performance of PALM-IVA with α = 1 to IVA-G with a vector-gradient update (IVA-G-V), and IVA-G with a block-Newton update (IVA-G-B), both of which are introduced in [START_REF] Anderson | Joint Blind Source Separation With Multivariate Gaussian Model: Algorithms and Performance Analysis[END_REF] and have matlab implementation available at https://mlsp.umbc.edu/resources.html. We use joint inter-symbol-interference (joint-ISI) used in [START_REF] Anderson | Joint Blind Source Separation With Multivariate Gaussian Model: Algorithms and Performance Analysis[END_REF] as the metric to compare separation performance of JBSS algorithms' demixing matrices W [k] when the true mixing matrices A [k] are known. Joint-ISI is normalized in [0 1], and jointly measures distance of each dataset's "mixing-demixing" matrix W [k] A [k] to a permuted diagonal matrix, with 0 joint-ISI representing perfect separation performance. The three algorithms were used with a stopping criteria of 10 -4 to ensure that the results were not biased, and we also added a limit of 5000 iterations.

We repeated twice the experiment described above: once with ρ = 0.2 to simulate low correlation within an SCV, representing a lower level of dependence across the datasets, and once with ρ = 0.8 to simulate high correlation within an SCV (and across the datasets). Indeed, it is known that IVA-G algorithms usually perform better if the datasets are strongly correlated, that is why we also tested the influence of the correlation across dataset on the performance of PALM-IVA. These experiments have been implemented in Python 3.10.8.

The results of the 100 simulations are displayed in Tables 1 and2. We indicated for each algorithm and for each pair (K, N ) the average computation time t in seconds, the mean ISI score µISI and the standard deviation of the scores σISI. For the mean and the standard deviation, we noted the best performance in bold characters.

  E.C. acknowledges support from the ERC Starting Grant MAJORIS ERC-2019-STG-850925. Part of the work of J.-C. P. was supported by the ANR Chair in AI, BRIDGEABLE. This work was also supported in part by NSF-NCS 1631838, and NIH grants R01 MH118695, R01 MH123610, R01 AG073949.

  [k]n } and the joint entropy of the entire SCV H{yn}, via I{yn} = K k=1 H{y[k]

Table 1: Results for low correlation across datasets (ρ = 0.2). [k] . Therefore, our method enables to find a proper solution to the IVA-G problem when there are only 2 datasets which are barely correlated. The standard deviation is also smaller for PALM-IVA, which confirms that this method is more consistent than IVA-G-B. PALM-IVA also achieves overall better performance than IVA-G-V, but these algorithms perform similarly when K = 5. The computation time however is the downside of PALM-IVA which, in our Python implementation, is about 10 times slower than IVA-G-B and 3 times slower than IVA-G-V. Despite this, computation time using PALM-IVA may be significantly improved by initializing PALM-IVA with an efficient JBSS algorithm, such as MCCA or the faster variants of IVA-G.

IVA-G-B

Table 2: Results for high correlation across datasets (ρ = 0.8). In the case of high correlation where JBSS is expected to be most useful, PALM-IVA still performs better than IVA-G-B and IVA-G-V in most scenarios. Namely, the mean ISI score is 2 to 3 times smaller with PALM-IVA than with IVA-G-B if K = 2.

CONCLUSION

We proposed a new algorithm based on alternating block proximal descent, which provides convergence guarantees to a critical point of the proposed IVA-G cost function, and performs competitively with respect to previous state-of-the art approaches to IVA-G.

In our future work, we would like to improve the computational efficiency of our approach, and try to adapt the PALM algorithm to versions of IVA with other assumed multivariate distributions, such as the multivariate Laplacian distribution assumed within IVA-L. Another lead would be to consider more general constrained formulations as alternative solutions to the regularization on the diagonal coefficients of the precision matrices we used. Finally, we would also like to consider applications of PALM-IVA to real data, such as fMRI datasets typically modeled by IVA.


Victor Vallejo-Otero

Arnaud Valour

Stéphanie Reynaud

Émilie Gamet

Christophe Donnet

Nicolas Crespo-Monteiro

Yves Jourlin

E Gamet

à la diffusion de documents scientifiques de niveau

Keywords: properties Keywords (max. 4-5): Titanium nitride, Micro-nanostructure, Nitridation, Optical lithography Disciplinary fields involved: Material Sciences, Optics

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

Innovative process of micro-nanostructure of titanium nitride (TiN) by association of a photostructurable titanium oxide (TiO2) and a rapid thermal annealing (RTA)

The titanium nitride (TiN) is usually used for its thermal, mechanical, electrical and optical properties. Its micro-nanostructure allows to use it in a wide range of applications especially in the area of optic [1] and photocatalysis [2]. These properties combined make this material a promising alternative of gold and silver for optical applications at high temperature or in extreme environment.

The fabrication of microstructured or nanostructured surfaces in TiO2 and in TiN on large substrates with plane or complex geometry is a challenge. Due to their high thermal, chemical and mechanical stability, TiN layer are difficult to micro-nanostructure by etching or engravement after their deposition. An alternative route is to micro-nanostructure titanium oxide film (TiO2) before converting it into TiN by a nitridation process [3,4].

In this presentation, a new method of elaboration of micro-nanostructured TiN layer from a TiO2 photostructurable sol-gel by rapid thermal annealing (RTA) will be presented. This sol-gel is photostructurable by optical lithography and nano-imprint lithography (NIL), allowing an easy and quick production of micro-nanostructure TiO2 layers (figure 1). These layers can then be converted into micronanostructure TiN layers thanks to an RTA of few minutes under ammonia flux. The micro-nanostructured films have been characterized by X-ray diffraction, Scanning Electron Microscopy (SEM), Spectroscopy RAMAN... The optical and electrical properties have been studied and a plasmonic application will be presented.