V Thoréton

P Rozier

E Capoen

O Mentré

B Boukamp

S Daviero-Minaud

Evidence of the Current Collector Effect: Study of the SOFC Cathode Material Ca 3 Co 4 O 9+d Ã

Keywords: Ca 3 Co 4 O 9+d, Cathode, Current Collection, Electrochemical Impedance Spectroscopy, Screen Printing, Solid Oxide Fuel Cell -

Introduction

In a solid oxide fuel cell (SOFC), the fundamental part of the cathode model is the oxygen reduction pathway. The mechanisms involved at the cathode at the triple phase boundary, where the oxygen reduction reaction occurs, are quite complex and are often modeled using multi-step mechanisms. They imply interactions, which can take place on the surface, but also in the whole thickness of the cathode layer. The steps to consider are oxygen diffusion, oxygen adsorp-tion/dissociation, charge transfer, and diffusion of ionic species [1]. They involve specific requirements for the considered materials in terms of microstructure, porosity, electrocatalytic properties, electronic/oxide ion conduction, and electrochemical properties. The cathode composition and morphology plays the main role. But one must keep in mind that the current collector can also interfere with the oxygen reduction reaction, especially if platinum is used. Indeed, its excellent catalytic properties toward this reaction justified its use as cathode in the first SOFCs developments [1][2][3].

-

[ * ] Corresponding author, aurelie.rolle@ensc-lille.fr

Abstract

In the study of the performance of solid oxide fuel cell (SOFC) electrodes, the possible influence of the applied current collector is often not mentioned or recognized. In this article, as part of an optimization study of the potentially attractive Ca 3 Co 4 O 9+d cathode material (Ca349), special attention is paid to the choice of current collector. The influence of both gold and platinum paste or grid current collectors on pure and composite (Ca349 + 30 wt.% Gd-doped ceria) is studied, using electrochemical impedance spectroscopy (EIS). Although, platinum is catalytically active in the oxygen reduction reaction and then is often considered as current collector for SOFC cathodes, in combination with Ca349 cathodes, additional low frequency dispersion is observed, leading to a larger polarization resistance than found in the case of gold current collectors. A subsequent experiment revealed that Pt reacts with Ca349, forming undesirable phases: CaPt 2 O 4 , Ca 4 PtO 6 , Ca 3 Co 2 O 6 , and Co 3 O 4 . The impedance spectra were analyzed with ZView 3.3a and with EqCwin v1.2. One series equivalent circuit was deduced using ZView, whereas, two possible equivalent circuits (series and nested), leading to the same quality of fits, were evidenced in EqCwin. The circuits are closely related to interactions of the current collector and layer thickness effects of the cathodes.

Within the framework of the study of cathodes for SOFCs, electrochemical impedance spectroscopy (EIS) is a method of choice to describe the electrochemical properties. However, drastic effects can be evidenced in the cathode electrochemical response according to the measurement conditions and more precisely the collector choice. In many studies using EIS to characterize the cathode properties, either the choice of current collector is not justified, or the type of current collector is not always specified. Generally platinum paste is used as current collector in the characterization of SOFCs' cathodes. We have recently pointed out as part of an optimization study of Ca 3 Co 4 O 9+d (Ca349, a potential cathode material [4,5]), that particular attention should be paid on the choice of the current collector.

Ca 3 Co 4 O 9+d , a thermoelectric material, is expected to have two-dimensional ionic conduction due to oxygen vacancies associated with the mixed valence of cobalt in 2D-rock salt subunits. It is also a good electronic conductor (r = 100 S cm -1 at 700 °C). It shows no reactivity with the gadolinium doped ceria standard electrolyte (CGO) and possesses a thermal expansion coefficient of 9-10 × 10 -6 K -1 , which makes it compatible with CGO [2].

In the article, all the abbreviations considered, the cells characteristics and the equivalent circuit symbols are provided in the Supplementary Data Section (S1, S2, S3). In Figure 1, the recently published results on Ca349 are summarized [4,5]. The same description of samples as in ref. [5] will be used in this paper. The sample designations are described in more details in the Supplementary Data Section (S2). From these results, one trend has been observed: the polarization resistance is lower for Ca349/CGO composites than for the single-phase electrodes. Moreover, the lowest polarization resistance is obtained for the thicker layers. The lowest polarization resistance is currently reported in ref. [4] for the 70 wt.% Ca349-30 wt.% CGO composite cathode: approximately 1 X cm -2 at 700 °C measured on thick electrolyte support.

From the impedance diagrams presented in Figure 1, it is obvious that the use of platinum as current collector raises the question of a possible current collector effect.

With the aim to highlight the importance of the choice and on the effect of the current collector, this article focuses on the influence of several current collectors (gold grid, gold paste, platinum grid, and platinum paste) on the electrochemical response of the Ca349 cathodes. The impedance data has been analyzed with ZView 3.3a [START_REF] Johnson | Z View version 3.3a[END_REF] and with EqCwin v1.2 (Windows version of Equivalent Circuit [START_REF] Boukamp | [END_REF][8][START_REF] Boukamp | Equivalent Circuit (version 1.2)[END_REF]), in which more detailed equivalent circuit models have been extracted.

Experimental

The synthesis of the Ca349 powder, the Ca349-based inks and the preparation of the pellets for the microstructural and electrochemical characterizations were performed as presented in ref. [5]. Cathode layers were screen-printed using a DEK 248 screen-printer. Depending on the type of sample/ current collector assembly thermal treatments at 700, 800, or 900 °C were applied. More information on each cell is provided in the Table S2 given in the Supplementary Data Section. For current collector deposition, Au paste (Ferro couleurs, France SA), mixed with DH 43 (dmc 2 Limoges, France) was deposited by hand painting on the cell, followed by a thermal treatment at 650 °C in order to decompose the organic part. The Pt paste (dmc 2 Limoges) was deposited similarly and treated at 800 °C. EDS analysis has revealed the purity of the gold and platinum paste.

Microstructural characterizations were performed by scanning electron microscopy using a Hitachi S3400N SEM. X-ray diffraction was performed on powders and on deposited layers using a D8 Advance A25 Bruker with a Position Sensitive Detector Lynxeye or using a D8 Bruker AXS diffractometer (CuKa = 1.54 Å). In order to determine the oxygen mobility in Ca349, permeation measurements were carried out, using the same protocol as described in ref. [START_REF] Capoen | [END_REF], on dense pellets (97%) obtained by spark plasma sintering (SPS), performed in the PNF2 Platform (T max = 850 °C and P max = 50 MPa, densification starting at 600 °C until 850 °C) (Toulouse, France).

Pellets have been checked by XRD. Impedance measurements were performed on symmetrical cells, cathode/electrolyte/cathode, using a SI 1260 Solartron FRA, during a heating-cooling sequence. An optimized signal amplitude of 50 mV (pure Ca349) and 70 mV (composite Ca349) were imposed, with respect to the linearity of the electrical response. Four types of current collectors were applied: gold grid, platinum grid, gold paste, and platinum paste. The data collection procedure was identical to the one reported in ref. [5]. At each temperature, the thermal and time stability were checked. The data were collected when the change in the polarization resistance was less than 5%. After a temperature change, measurements were performed as function of time. It was found that stability is reached within 2 h after each temperature change. The impedance data were analyzed using two different software packages: ZView 3.3a [START_REF] Johnson | Z View version 3.3a[END_REF] and Equivalent Circuit 1.2. (Windows version EqCwin) [START_REF] Boukamp | [END_REF][8][START_REF] Boukamp | Equivalent Circuit (version 1.2)[END_REF]. For use of the latter programme, the stability and validity of the considered data was checked using the Kramers-Kronig (KK) test for Windows, version 1.01 [11,12].

Results and Discussion

Ionic Conduction Properties

Ca349 is known for its good electrical conductivity property, which is reported to be higher than 100 S cm -1 over the 20-800 °C temperature range [4,13]. However, the ionic conduction of this phase has never been checked before. For that purpose, permeation measurements were performed on 0.5 mm thick and 15 mm diameter SPS-dense pellets (relative density >97%) of the pure Ca349 material. Both faces of the pellets were polished with up to 2,400 grid and the experiments were carried out between 580 and 650 °C. The oxygen flux is deduced from the DpO 2 variation at the permeate side. The oxygen permeation measurements revealed a significant oxygen flux even at low temperatures (0.0014 lmol cm -2 s -1 at 600 °C), which unquestionably proves the existence of oxygen-ion migration (Table 1). For information, the oxygen fluxes reported for La 0.6 Sr 0.4 Co 0.8 Fe 0.2 O 3-d (LSCF) in the literature [14] are about 0.00562 lmol cm -2 s -1 at 850 °C, but one can notice, that the membrane thickness, pO 2 gradients, and temperature range are different, making a direct comparison difficult. Nevertheless the flux measured for Ca349 is about the same order, but measured at significantly lower temperatures. However, after the permeation measurements, a surfacial passivation of Ca349 into CaO and CoO, was observed due to the long-exposure time of the pellet under Ar atmosphere at 700 °C. This experiment confirms Ca349 as a mixed oxide ion and electronic conductor and shows its potentiality as cathode material in SOFC applications.

Current Collector Effect

Gold Paste Versus Gold Grid

The same composite sample CC800 20 characterized in ref. [5] with gold grid only (CC800 20 ) was also studied here with gold paste (CC800Au 20 ) as current collector. The comparison between these two kinds of current collectors revealed that the electrochemical properties are the same, independent of the type of gold current collector. For example, the comparison of the impedance spectra recorded at 700 °C is shown in Figure 2(a). A random shift of the series resistance is observed on the impedance spectra, reported Figure 2, which might be due to the global measuring set-up.

For this 20 lm thick cathode layer, polarization resistances of 1.39 and 1.41 X cm 2 were obtained for the gold grid and the gold paste, respectively. The cross-section of the cell after EIS measurements is presented in Figure 3(a). The gold layer is about 5 lm thick. After the EIS measurements no additional phase could be observed in the X-ray diffraction pattern [Figure 4(a)]. The reactivity test between Ca349 and Au at 700 °C for 72 h also revealed no interaction between these two phases, as can be noticed in the X-ray diffraction pattern of Figure 4(b). The impedance plots consist of a single depressed semicircle which does not change noticeably with time [Figure S4 in Supplementary Data Section]. From these results, it can be inferred that the impedance response represents the properties of Ca349 only.

Platinum Paste Versus Gold Grid

A comparison of the electrochemical properties of identical composite CC800 5 cathodes with either platinum paste or gold grid as current collector is presented in Figure 2(b). These results have been reproduced from ref. [5]. Activation energies determined for the case of gold grid and platinum paste current collectors (1.15 and 1.52 eV, respectively), reveal that, when platinum paste is used as current collector, the cathode behavior is more complex and the ASR values are less favorable [5]. Special attention is paid to the impedance properties at 700 °C. At this temperature, polarization resistances of about 5.6 and 12.9 X cm 2 are obtained, respectively, for gold grid and platinum paste used as current collector [5] with the same sample. Even if the platinum layer is quite thin [about less than 1 lm, Figure 3(b)], a clear growing contribution is observable in the impedance plot at low frequency, as compared to the gold grid collector case. This Temperature / °C 600-700 850-1000 J(O 2 ) / lmol cm -2 s -1 0.0014-0.0020 0.00562-0.0316 Thickness / mm 0.50 0.55-0.98 p′O 2 / atm 1.00 0.21 3 × 10 - 6 1 × 10 -3 growing contribution is observed in the impedance plot of all the samples (pure as well as composite ones) measured with a Pt-paste current collector, independently of the platinum thickness [5]. X-ray diffraction was performed on the CC800 5 pellets after EIS characterizations, in order to get a better understanding of the difference in behavior as observed between the gold grid and the platinum paste current collectors. For the Pt-paste, the X-ray analysis revealed the formation of undesirable phases [Figure 5 5(b andc)], confirming a chemical reaction between Ca349 and Pt. A study, whose results have never been published before [START_REF] Preux | A la Recherche de Nouveaux Matériaux D'électrolyte et de Cathode Pour SOFC : Weberite et Cobaltite[END_REF], has revealed that the Ca 3 Co 2 O 6 phase, which consists of a one-dimensional structure, has quite poor conducting properties and cathode behavior, as compared to Ca349, with a polarization resistance of 96.2 X cm -1 at 700 °C (Figure 6). The growing contribution observed in the impedance plot could thus be assigned to undesirable phases resulting from the reaction on the surface between Ca349 and the platinum paste.

Platinum Grid

On sample CC700 20 , of which the layer microstructure is presented in Figure 3(d), a platinum grid collector was tested. In comparison with the gold grid case, a small deviation at the low frequency end of the depressed semicircle can be discerned [Figure 2(c)]. The low frequency semi-circle starts to appear. This trend is intermediate between the electrochemical behavior as observed using gold grid and platinum paste. Again this may be related to the Ca349/Pt reaction. X-ray diffraction performed after the impedance spectroscopy measurements confirms a beginning of reaction between Ca349 and the platinum grid leading to Co 3 O 4 , CaPt 2 O 4 , and Ca 3 Co 2 O 6 formation, as indicated by additional peaks to those of the Ca349 and CGO phases (Figure 7).

Fit Procedures

Evidence of the Ca349 Equivalent Circuit. Comparison of Z View Versus EqCWin

We have first focussed on the pure compounds PC800c 10 , measured with a gold grid. In this article, two fitting softwares are considered and compared to analyze the same data set: Z View 3.3a [START_REF] Johnson | Z View version 3.3a[END_REF] and EqCwin v1.2 [START_REF] Boukamp | [END_REF][8][START_REF] Boukamp | Equivalent Circuit (version 1.2)[END_REF]. Using EqCwin v1.2, the approach consists of first checking the data quality using the KK test for Windows, version 1.01 [11,12]. This test, which should be performed directly after an impedance measurement, provides a check for the validity and stability of the data. The quality of the data can be observed in the residuals plot of the KK-test. The relative residuals, D real and D imag , are defined as [11]:

D real x i Z real x i À KK real x i Zx i j j ; D imag x i Z imag x i À KK imag x i Zx i j j
where KK(x i ) is the Kramers-Kronig test function.

Valid data should present a random distribution of the relative residuals around the log(x i ) axis, representing basically the noise distribution of the data. A systematic deviation from the frequency axis indicates non-KK behavior, which can be caused by several factors. The most important reasons are a nonequilibrium state of the cell, for example, the system has not completely adjusted to the last temperature change or a change in the ambient composition. But it can also be caused by a strong non-linear transfer function, which is inherent to electrode responses when a too large perturbation voltage is used [11,12]. In general, for a good CNLS-fit the residuals distribution should be similar in shape as observed for the KK-test [12]. Since the KK-test also indicates the frequency limits that will allow an optimal fit procedure, a significant difference between the residuals plot for the KK-test and the final CNLS-analysis can indicate that the data contains either an additional contribution or that the used equivalent circuit is not appropriate.

Whereas, ZView 3.3a uses a (semi) automatic search for the equivalent circuit with the "best fit" to the data set, EqCwin provides a partial fit and subtraction procedure with which a possible equivalent circuit is built and starting values for the fit parameters are obtained. These are subsequently used for a final overall CNLS-fit of the complete data set. The clear advantage of this procedure is that minor contributions, that are not obvious from visual inspection of the data, are readily detected [START_REF] Boukamp | [END_REF]12]. Final comparison of the residuals distributions between the KK-test and the CNLS-fit can confirm the viability of the derived equivalent circuit to model the data set. In Figure 8 (a andb) the KK-tests indicate for the data of sample PC800c 10 a perfect KK behavior.

Two possible equivalent circuits can be derived using EqCwin: a series circuit and a nested circuit. A very good match between data set and model function can be observed, as demonstrated in Figure 8(c andd), where the final CLNSfit of the data measured at 700 °C is presented. The corresponding circuit element values and error estimates are presented in Table 2. The nested arrangement assumes a double layer capacitance between the electrode and the electrolyte, which results in a parallel capacitance to the electrode response. Both circuits have been checked for various temperatures, various compositions (composite 70% Ca 3 Co 4 O 9 -30% CGO, pure Ca 3 Co 4 O 9 ), various layer thicknesses, various current collectors (gold grid, gold paste, and platinum paste), various oxygen partial pressure, and leads to good fits of similar quality. The extracted circuits correspond to the minimal number of contributions necessary for a good fit quality. Adding an extra circuit does not automati-cally improve the fit. Inappropriate additions to the EqCwin will generally result in a shift of the element parameters to either zero or infinity during the fit procedure. Particular focus should be paid on the errors obtained in ZView and EqCwin which are not calculated in the same way, then, the obtained values cannot be compared. The errors in ZView are based on the number of experimental points considered, whereas, the error estimates in EqCwin are based on the assumption that the residual errors can be assumed to represent the parent distribution [8,12]. This last point requires that the residuals are randomly distributed around the log(x) axis. A systematic deviation from 0 around the log(x) axis invalidates the calculation of the error estimates.

Here, as example, we present Figure 8 and Table 2 the results concerning only the impedances data measured at 700 °C under air atmosphere. The specific resistances associated to a frequency range are of the same order of magnitude, whatever the type of circuit, series or nested. That means, with the chosen notation: (R 1 ) series = (R 1 ) nested ; (R 2 ) series = (R 4 ) nested ; (R 3 ) series = (R 2 ) nested ; (R 4 ) series = (R 3 ) nested ; and (R 5 ) series = (R 5 ) nested . Both circuits involve the same number of parameters. The series circuit imply successive steps in the cathode reaction, whereas, the nested one, consider one step occurring at the same time as the other ones.

Considering the ZView 3.3a program, an equivalent circuit consisting of a series association of a finite Warburg (W) and (R//CPE) elements is extracted, corroborating the results of previous studies [4]. A comparison between the fitted model and the experimental data is provided in the Supplementary Section S5. The corresponding circuit elements are given Table 2.

The two fit procedures, ZView and EqCwin, lead in the case of PC800c 10 to quite different model circuits. But for a comparison purpose, the series circuits extracted from EqCwin has been applied in ZView (Supplementary data S6). However, a small improvement of the fit is observed because of the higher number of fitted elements. But the additional part of the circuit is associated to a too big error to be considered. By constructing progressively the circuit, using various data in function of the temperature, the first proposed ZView circuit with the corresponding errors is more justified. Moreover, if the nested circuit is applied in ZView, the fit is not stable (Supplementary data S6).

As the series and the nested circuits extracted from EqCwin lead to the same fit quality and as the application of the nested circuit in ZView is unstable, for an easier comparison between the two softwares, the series circuit only will be considered for the rest of the study. Knowing that CCO is a mixed ionic-electronic conductor, we can expect several intermediate mechanisms, involving the surface reaction, the electrode-electrolyte interface diffusion, and the diffusion in the cathode itself. They cannot be clearly identified using the elementary ZView circuit with two (even three) contributions, which provide a too rough description. In that case, the subtraction technique developed in EqCwin leads to a more complex circuit consisting of four contributions corresponding to the best fit quality with the minimal number of contri-

R 1 L 1 W s (R 2 //CPE 1 ) [5] R 1 L 1 W s (R 2 //CPE 1 ) (R 3 //CPE 2 ) L 1 R 1 (R 2 Q 1 )(R 3 Q 2 )(R 4 Q 3 )(R 5 Q 4 ) L 1 R 1 (Q 1 [R 2 (R 3 Q 2 )(R 4 Q 3 )(R 5 Q 4 )]) t = 696 °C t = 696 °C t = 696 °C t = 696 °C Circuit elements
Circuit elements Circuit elements Circuit elements R 1 = 3,557 X (0.1%) R 1 = 3,557 X (0.1%) L 1 = 6.5 × 10 -7 Henri (0.1%) L 1 = 6.5 × 10 -7 Henri (0.1%) L 1 = 6.5 × 10 -7 Henri (0.1%) L 1 = 6.5 × 10 -7 Henri (0.1%) butions, which may provide additional information on the cathode behavior.

R 1 = 3.6 X (0.4%) R 1 = 3.7 X (0.2%) W s -R = 3.64 (2%) W s -R = 2.10 (3%) R 2 = 1.7 X (35.2%) Q 1 -Y 0 = 2.2 × 10 -4 mho (104.8%) W s -T = 0.0046 (5%) W s -T = 0.0018 (5%) Q 1 -Y 0 = 2.3 × 10 -2 mho (30.3%) Q 1 -n = 7.4 × 10 -1 (10.3%) W s -P = 0.336 (0.5%) W s -P = 0.309 (1%) Q 1 -n = 9.3 × 10 -1 (4.1%) R 2 = 3.1 × 10 -1 X (29.6%) R 3 = 6.2 X (79.2%) R 3 = 2.6 X (40.0%) R 2 = 3.33 X (2%) R 2 = 2.9 X (8%) Q 2 -Y 0 = 2.
To have a better understanding of the current collector and of the layer thickness, the same analysis approach will be applied on selected composite cathodes. We have focused on a cell with a composite Ca349 cathode, measured successively with a gold grid and with platinum paste (cathodes: CC800 5 and CC800Pt 5 ) and second cell with a thicker cathode layer (20 lm), measured successively with a gold grid and gold paste (cathodes: CC800 20 and CC800Au 20 ). A detailed com-parison is performed on the CNLS-fits of the impedance data collected at 700 °C.

Using ZView, the fit procedure considers the whole data set and leads to the following series circuit: R 1 L 1 W s (R 2 // CPE 1 ) [4] and R 1 L 1 W s (R 2 //CPE 1 ) (R 3 //CPE 2 ). In the case of Pt paste as current collector the contribution (R 3 //CPE 2 ) obtained in the low frequency range is widely increased, compared to Au grid collector. No particular information of the cathode behavior can be extracted. The ZView and CNLS-fit results are provided in Table 3 and Figure 9. One Table 3 CC800 5 , CC800Pt 5 , CC800 20 , and CC800Au 20 : fitted circuit elements using ZView 3.3a software.

Fit results 700 °C-ZView 3.3a CC800 5 [5] Fig. 9 Experimental impedance plots and the corresponding fits using ZView 3.3a for CC800 5 , CC800Pt 5 , CC800 20 , and CC800Au 20 .

must notice, in case of cells measured with the gold paste or gold grid as current collector, that the last part of the circuit (R 3 //CPE 2 ) is not necessarily justified, the associated errors being high for a slightly improved R p value. The equivalent circuit evidenced in ZView is in good agreement with the circuits reported in refs. [4,5] for the pure and composite compounds.

For the subtraction procedure, the quality of the data is first checked using the KK test. The data are of quite good quality, except for CC800Pt 5 , where a larger noise contribution is observed in the low frequency region. The KK-test residuals and CNLS-fit results are presented in Figures S7 andS8 of the Supplementary Data Section. The circuits extracted pre-viously on the pure compound, measured with a gold grid as current collectors, were applied to the composite compounds studied in this part. For a better comparison with ZView fits, the series equivalent circuit will be considered, as it is equal to the other one with respect to the fit quality. The series equivalent circuit

L 1 R 1 (R 2 Q 1 )(R 3 Q 2 )(R 4 Q 3 )(R 5 Q 4 ) describes
properly CC800 5 , CC800 20 , and CC800Au 20 . One can notice that the same circuit can be used for CC800Pt 5 . The fit results are presented in Table 4 and the corresponding CNLS-fits compared to the experimental data in Figure 10. The first point to notice is that the pure and the composite compounds show the same behavior, as evidenced by the same equivalent circuits. This means that the same physical phenomena are involved independently of the composition of the cathode: pure Ca349 or composite 70 wt.% Ca349-30 wt.% CGO.

= 702 °C L 1 R 1 (R 2 Q 1 )(R 3 Q 2 )(R 4 Q 3 )(R 5 Q 4 ) L 1 R 1 (R 2 Q 1 )(R 3 Q 2 )(R 4 Q 3 )(R 5 Q 4 ) L 1 R 1 (R 2 Q 1 )(R 3 Q 2 )(R 4 Q 3 )(R 5 Q 4 ) L 1 R 1 (R 2 Q 1 )(R 3 Q 2 )(R 4 Q 3 )(R 5 Q 4 ) L 1 = 6.
In the following sections, the focus is on the subtraction procedures using EqCwin, considering the series equivalent circuit. Some large errors can be noticed. But they are just a consequence of the cross-correlation between parameters. From a general point of view, the errors in the ZView fit appear to be lower, but this is due to the systematic deviation in the residual plots. Knowing that, a correlation between the circuit contributions and physical phenomena has been suggested.

Behavior of CC800 5 Versus CC800Pt 5

As presented above, a comparison of the same CC800 cathodes, either with gold grid or with platinum paste current collectors, reveals that two different equivalent circuits: a series and a parallel circuit, could be considered.

As frequency decreases, the resistances appear in the order: R 1 , R 3 , R 4 , R 2 , and R 5 . The resistance parallel to a capacitor chosen to describe the low frequency contribution is related to surface phenomena associated either with the platinum paste or the degradation phases. The high and medium frequency semicircles correspond to the electrochemical behavior associated with the combination of Ca349 and the degradation phases, which constitutes a new cathode. This changed cathode composition leads to quite different electrochemical properties. One can notice that R 1 , R 2 , and R 4 are of the same order of magnitude, whereas, using platinum paste as current collector leads to an increase of R 3 and R 5 .

Behavior of CC800 20 Versus CC800Au 20

CC800 20 and CC800Au 20 cathodes responses are both described by a series equivalent circuit. The characteristic resistances appear, with decreasing frequency, in the order R 1 , R 3 , R 4 , R 2 , and R 5 . Special attention paid to the fitted values of the circuit elements (Table 4) reveals that the same order of magnitude is observed for both cathodes. This confirms that the same electrochemical reaction applies when using gold paste or gold grid current collectors.

Behavior of CC800 5 Versus CC800 20

These two samples differ in their layer thickness. CC800 5 and CC800 20 have, respectively, a thickness of 5 and 20 lm. Assuming that each contribution in the equivalent circuit can be associated with different steps in the overall electrochemical cathode reaction, one could assign R 3 and R 4 to diffusion phenomena in the cathode layer according to their frequency location on the plot in the 13-10 000 Hz range, approximately, corresponding to the "C arc" mentioned in the ref. [START_REF] Jorgensen | [END_REF]. Similarly R 5 can be assigned to surface phenomena and R 2 to the electrode-electrolyte interface impedance. For both samples, the R 5 values and associated time constants are of the same order of magnitude. For R 2 , differences are evidenced for R values, but the associated time constants are similar. For R 3 and R 4 , differences are noticed for R values and time constant. It seems to indicate that the surface and interface phenomena are not affected by the thickness effect. However, the diffusion phenomenon is affected by the thickness effect. The same physical phenomenon is involved and the thickness affects highly the resistances values. The thick layer shows the lowest polarization resistance. This indicates the importance of the mixed conduction properties of the material and the increased triple phase boundary length which is essential for the oxygen reduction reaction.

Conclusion

Ca 3 Co 4 O 9-d (Ca349) is a mixed ionic, electronic conductor with interesting electrochemical properties, which makes it a good candidate cathode for application in a SOFC. Continued investigations are underway to improve the properties and to understand the electrochemical behavior of this material. Complex mechanisms are involved at the cathode which are much more complicated in case of the composite cathode/ electrolyte. EIS is a key technique to describe the cathode behavior. However, special attention should be paid to the choice of the current collector, which depends on the composition under study. This study, focused on Ca349, revealed that in this case the best current collector is a gold grid. The commonly used platinum paste or grid is here inappropriate as current collector for the SOFC cathode. Reactions between Ca349 and platinum lead to the formation of several undesirable phases: CaPt 2 O 4 , Ca 4 PtO 6 , Ca 3 Co 2 O 6 , and Co 3 O 4 . With platinum paste as current collector, the electrode impedance clearly shows two contributions, whereas, with a gold grid current collector only one contribution is observed. Two fit procedures, ZView 3.3a and EqCwin v1.2 were applied for the impedance data analysis. The subtraction procedure provided in EqCwin yields more detailed information on the equivalent circuits than ZView. The CNLS-analysis has made it possible to study the effect of current collector and layer thickness, as part of the study of electrode materials. To confirm the correlation between the circuit elements and the physical phenomena, further electrochemical studies considering various electrode compositions and microstructural properties will be performed under various conditions (temperature, atmosphere, and current). A perspective of the work could be also to describe the actual electrode structure with finite element modeling, as reported by Kornely et al. [17].

R 1 //CPE 1 : R1 CPE1 L 1 R 1 (R 2 Q 1 )(R 3 Q 2 )(R 4 Q 3 )(R 5 Q 4 ): L1 R1 R2 QPE1 R3 QPE2 R4 QPE3 R5 QPE4 L 1 R 1 (Q 1 [R 2 (R 3 Q 2 )(R 4 Q 3 )(R 5 Q 4 ]):

Supplementary Data Section

Ca349 and LSCF refer to Ca

Fig. 1

 1 Fig.1Recently published data on Ca349[4,5].

  (a)]. To confirm this result, a reactivity test was then performed on pellets prepared by mixing Ca349 and Pt powder and firing at 700 °C for 72 h. Additional peaks corresponding to phases, such as CaPt 2 O 4 , Ca 4 PtO 6 , Ca 3 Co 2 O 6 , and Co 3 O 4 are highlighted [Figure

Fig. 2

 2 Fig. 2 Effect of the current collector on the electrochemical property: (a) gold paste-gold grid, (b) platinum paste-gold grid, and (c) platinum grid. Data are reported at 700 °C.

Fig. 3

 3 Fig. 3 SEM cross-section of (a) CC800Au 20 , (b) CC800Pt 5 , (c) PC900cPt 10 , and (d) CC700 20 after EIS characterization.

Fig. 4 (

 4 Fig. 4 (a) Comparison of the X-ray diffractograms of Ca349 powder, CC800 20 , and CC800Au 20 before gold painting (bef.) and after EIS measurements (aft.). (b) Reactivity test of Ca349 and gold at 700 °C for 72 h.

Fig. 5 (Fig. 6

 56 Fig. 5 (a) Comparison of the X-ray diffractograms of Ca349 powder, CC800 5 , and CC800Pt 5 before Pt painting (bef.) and after EIS measurements (aft.). (b and c) Reactivity test of Ca349 and platinum at 700 °C for 72 h.

Fig. 7 X

 7 Fig. 7 X-ray diffractograms of CC700 20 after EIS measurements.

Fig. 8 (

 8 Fig. 8 (a) KK tests for PC800c 10 . The red plot corresponds to the real part differences and the blue plot to the imaginary part differences. (b) Ideal fits possible to get using the considered data set for PC800c 10 . Fit is represented in red and impedance data in black. (c and d) Comparison of the experimental and fitted impedance plots of PC800c 10 , using the subtraction procedure on Equivalent Circuit 1.2 software, considering the series and the nested equivalent circuit, respectively. The corresponding residual plots are reported in insert. Fit are represented in red cross and impedance data in blue square.

  3 × 10 -2 mho (45.4%) Q 2 -Y 0 = 6.9 × 10 -3 mho (18.9%) CPE 1 -T = 0.0124 (2%) CPE 1 -T = 0.0071 (7%) Q 2 -n = 4.1 × 10 -1 (11.1%) Q 2 -n = 5.4 × 10 -1 (4.4%) CPE 1 -P = 0.797 (1%) CPE 1 -P = 0.71 (2%) R 4 = 2.6 X (76.9%) R 4 = 2.2 X (62.6%) Q 3 -Y 0 = 5.7 × 10 -3 mho (17.6%) Q 3 -Y 0 = 8.3 × 10 -3 mho (9.1%) R 3 = 1.8 X (11%) Q 3 -n = 6.0 × 10 -1 (12.3%) Q 3 -n = 7.4 × 10 -1 (14.5%) CPE 2 -T = 0.021 (7%) R 5 = 2.0 X (105.0%) R 5 = 1.7 X (25.8%) CPE 2 -P = 0.91 (2%) Q 4 -Y 0 = 9.4 × 10 -3 mho (18.3%) Q 4 -Y 0 = 2.1 × 10 -2 mho (19.0%) Q 4 -n = 7.8 × 10 -1 (21.2%) Q 4 -n = 9.3 × 10 -1 (3.6%) v 2 = 6.3 × 10 -5 v 2 = 6.3 × 10 -6 v 2 = 3.6 × 10 -7 v 2 = 5.5 × 10 -7

7 Fig. 10

 710 Fig. 10 Experimental impedance plots and the corresponding fits using EqCwin 1.2 for (a) CC800 5 , (b) CC800Pt 5 , (c) CC800 20 , and (d) CC800Au 20 . The corresponding residual plots are in insert.

S4

  Evolution of gold grid cell CC800Au 20 , before and after 7 days measurements. S5 Comparison of the experimental and fitted impedance plots of PC800c 10 using ZView 3.3a sofware. Fit are represented in red line and impedance data in black square. S6 Test in ZView of the series and the nested equivalent circuits extracted using EqCwin. S7 Kramers-Kronig tests for (a) CC800 5 , (b) CC800Pt 5 , (c) CC800 20 , and (d) CC800Au 20 . The red plot corresponds to the real part differences and the blue plot to the imaginary part differences. S8 Ideal fits possible to get using the considered data set of (a) CC800 5 , (b) CC800Pt 5 , (c) CC800 20 , and (d) CC800Au 20 .

Table 1

 1 Oxygen permeation measurements: comparison of Ca349 and LSCF.

Membrane

Ca 3 Co 4 O 9+d (Ca349) La 0.6 Sr 0.4 Co 0.8 Fe 0.2 O 3 (LSCF)

[13]

Table 2 PC800c 10

 210

	ZView 3.3a

: comparison of the fitted circuit elements using ZView 3.3a and Equivalent Circuit 1.2 softwares.

Table 4 .

 4 CC800 5 , CC800Pt 5 , CC800 20 , and CC800Au 20 : fitted circuit elements using Equivalent Circuit 1.2 software.

	CC800 5	CC800Pt 5	CC800 20	CC800Au 20
	t = 714 °C	t = 701 °C	t = 701 °C	t

  3 Co 4 O 9+8 and La 0.6 Sr 0.4 Co 0.8 Fe 0.2 O 3 , respectively. EIS refer to electrochemical impedance spectroscopy. S1 Legend of the abbreviations considered.

	Composition	Treatment t / °C Support	Current collector Thickness / lm	Name	Reference
	Pure Ca 3 Co 4 O 9+d	700 °C	Dense CGO	Au grid	~40	PC700 40	[4]
	Composite 70 wt.%	700 °C	Dense CGO	Au grid	~15	CC700 15	[4]
	Ca 3 Co 4 O 9+d -30 wt.% CGO
	Pure Ca 3 Co 4 O 9+d + 20 wt.%	800 °C	Dense CGO	Au grid	5-10	PC800c 10	[5]
	cellulose
	Pure Ca 3 Co 4 O 9+d + 20 wt.%	900 °C	Dense CGO	Pt paste	~10	PC900cPt 10	[5]
	cellulose
	Composite 70 wt.%	800 °C	Dense CGO	Au grid	~5	CC800 5	[5]
	Ca 3 Co 4 O 9+d -30 wt.% CGO
	Composite 70 wt.%	800 °C	Dense CGO	Pt paste	~5	CC800Pt 5	[5]
	Ca 3 Co 4 O 9+d -30 wt.% CGO
	Composite 70 wt.%	800 °C	Dense CGO	Au grid	~20	CC800 20	[5]
	Ca 3 Co 4 O 9+d -30 wt.% CGO
	Composite 70 wt.%	800 °C	Dense CGO	Au paste	~20	CC800Au 20	This work
	Ca 3 Co 4 O 9+d -30 wt.% CGO
	Composite 70 wt.%	700 °C	Dense CGO	Pt grid	~20	CC700 20	This work
	Ca 3 Co 4 O 9+d -30 wt.% CGO
	S2 Cell characteristics.

Acknowledgements

The PHC Alliance exchange programme is acknowledged for funding the collaboration with Prof. Irvine's group, at the University of St Andrews. We thank the group for performing screen-printing depositions. K. Nagasawa and H. Nakatsugawa are thanked for initiating the work on the Ca349 material. Nicolas Preux is acknowledged for the results on Ca 3 Co 2 O 6 , which are part of his PhD.


Baptiste Lafabregue

Jonathan Weber

Pierre Gançarski

Germain Forestier
email: germain.forestier@uha.frgancarski@unistra.fr

Deep constrained clustering applied to satellite image time series

Keywords: time series, satellite images, remote sensing, clustering, constraints, deep embedded clustering

The advent of satellite imagery is generating an unprecedented amount of remote sensing images. Current satellites now achieve frequent revisits and high mission availability and provide series of images of the Earth captured at different dates that can be seen as time series. Analyzing satellite image time series allows to perform continuous wide range Earth observation with applications in agricultural mapping, environmental disaster monitoring, etc. However, the lack of large quantity of labeled data generally prevents from easily applying supervised methods. On the contrary, unsupervised methods do not require expert knowledge but sometimes provide poor results. In this context, constrained clustering, which is a class of semi-supervised learning algorithms, is an alternative and offers a good trade-off of supervision. In this paper, we explore the use of constraints with deep clustering approaches to process satellite image time series. Our experimental study relies on deep embedded clustering and the deep constrained framework using pairwise constraints (must-link and cannot-link). Experiments on a real dataset composed of 11 satellite images show promising results and open many perspectives for applying deep constrained clustering to satellite image time series.

Introduction

Deep learning is widely used in a large panel of domains, and sees a high interest in the remote sensing community [START_REF] Zhu | Deep learning in remote sensing: A comprehensive review and list of resources[END_REF]. It achieves great results but is highly dependant on the amount of available data, and more specifically labeled data. Remote sensing produces a large amount of data which often lacks of annotations. This problem is even more pregnant for times series of satellite images. Satellite images can be freely acquired every 5 days, however, due to the complexity of the data and the lack of a well defined typology, it is difficult to use supervised approaches. Therefore multiple clustering methods are applied in the domain [START_REF] Khiali | Detection of spatio-temporal evolutions on multi-annual satellite image time series: A clustering based approach[END_REF][START_REF] Rey | Investigating lake-area dynamics across a permafrost-thaw spectrum using airborne electromagnetic surveys and remote sensing time-series data in Yukon Flats, Alaska[END_REF]. But even if fully labelled data are not available, there are still a lot of background knowledge available from experts.

Constraints are a good way to exploit this knowledge, in this paper we will focus only on pairwise constraints, must-link (ML) and cannot-link (CL). They are widely used and well studied [START_REF]Constrained clustering: Advances in algorithms, theory, and applications[END_REF], which makes it easier to have a comparison baseline, as they are implemented by a lot of different methods. These constraints indicate that two instances should be assigned to the same cluster (must-link) or that they should be assigned to different clusters (cannot-link).

In this paper, we want to study if we can benefit from the advances in deep learning through a constraint-based framework that seems more appropriate to the remote sensing domain. We first introduce related works on clustering for time series and constrained clustering, respectively in section 2.1 and 2.2, then we present the deep constrained clustering framework and its adaptation to time series in section 3, then we compare the results to classic constrained clustering on remote sensing data in section 4. Finally, we discuss our results and the future works in the section 5.

Related Work

Clustering for time series

Different approaches for clustering of time series have been proposed through time, mostly based on some representation methods like Discrete Wavelet Transform [START_REF] Chan | Efficient time series matching by wavelets[END_REF] or similarity measures like Dynamic Time Warping [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF]. Those methods are then usually used to be incorporated into some standard clustering algorithms from the k-means, k-medoid, spectral or hierarchical clustering families, as shown in the review by Aghabozorgi et al. [START_REF] Aghabozorgi | Time-series clusteringA decade review[END_REF]. In this field, some improvements are still made, with, e.g., the k-shape algorithm [START_REF] Paparrizos | k-shape: Efficient and accurate clustering of time series[END_REF], which is based on an iterative refinement procedure that uses a normalized version of the crosscorrelation measure. One of the difficulties with time series is the heterogeneity of the related topics and the types of data, starting from the number of features or the sequences length, to the type of correlation between elements, based on shape or structure, with different amplitudes.

This issue has been tackled in the supervised learning by the rise of representation learning through deep neural networks. Recently some deep clustering approaches have been proposed. They are essentially based on some end-to-end architectures that simultaneously learn an embedding for the data and a clustering assignment, through an autoencoder and a clustering layer plugged on the encoder output [START_REF] Xie | Unsupervised deep embedding for clustering analysis[END_REF][START_REF] Guo | Improved Deep Embedded Clustering with Local Structure Preservation[END_REF]. A derived architecture have been developed for time series, that uses a 1D-CNN followed by a Bi-LSTM as autoencoder, to preserve the time dimension in the encoded features, they are then clustered by a similarity metric as a clustering layer [START_REF] Sai Madiraju | Deep Temporal Clustering: Fully Unsupervised Learning of Time-Domain Features[END_REF].

Constraints in clustering

A lot of works have been done on constraints integration for clustering. Most of them rely on some extension of standard clustering algorithm like k-means [START_REF] Wagstaff | Constrained kmeans clustering with background knowledge[END_REF] or spectral approaches [START_REF] Li | Constrained clustering via spectral regularization[END_REF], but there are also some dedicated methods like constrained programming clustering [START_REF] Duong | Constrained clustering by constraint programming[END_REF]. A comparative study has been done on the subject on time series [START_REF] Lampert | Constrained distance based clustering for time-series: a comparative and experimental study[END_REF].

In the deep learning domain, most of the semi-supervised methods refer to self-learning approaches or other way to include knowledge in supervised task. But recently some works have been proposed to include pairwise constraints in deep clustering networks [START_REF] Zhang | Deep Constrained Clustering-Algorithms and Advances[END_REF][START_REF] Ren | Semi-supervised deep embedded clustering[END_REF]. Both of these papers use constraints as an input for an extended loss function that tends to maximize the similarity of the encoded values between elements of a must link constraints and respectively minimize it for a cannot link constraints. Our work is based on the paper of Zhang et al. [START_REF] Zhang | Deep Constrained Clustering-Algorithms and Advances[END_REF]. Most of our work consists in adapting this method to time series and study its results on satellite image time series. This method handles different types of constraints, but, as indicated above, we focus on ML and CL in this paper.

Deep Constrained Clustering Framework and its adaptation to time series

The Deep Constrained Clustering framework (DCC) presented by Zhang et al. [START_REF] Zhang | Deep Constrained Clustering-Algorithms and Advances[END_REF] is based on a deep clustering method, the Deep Embedded Clustering (DEC) [START_REF] Xie | Unsupervised deep embedding for clustering analysis[END_REF] and its improved version (IDEC) [START_REF] Guo | Improved Deep Embedded Clustering with Local Structure Preservation[END_REF]. We first describe the IDEC framework, then how it is extended with constraints through the DCC framework and finally how we adapted it to time series.

Improved Deep Embedded Clustering

Deep Embedded Clustering (DEC) [START_REF] Xie | Unsupervised deep embedding for clustering analysis[END_REF], in the initial step, trains an autoencoder (x i = g(f (x i ))) and then removes the decoder. Then, it fine-tunes the remaining encoder (z i = f (x i )) by optimizing the Kullback Leiber divergence between two distributions Q and P . Q is a soft cluster assignment, where for each instance i we compute a vector q i of length k, k being the number of expected clusters, where q ij is the degree of belief that the instance i belongs to cluster j. P is the target distribution that is defined from Q as a "hard" assignment vector that assign i to only one cluster. We have the following loss, L c , as clustering loss:

L c = KL(P |Q) = i j p ij log p ij q ij (1)
where q ij is the similarity between the embedded point z i and cluster centroid µ j measured by Student's t-distribution [START_REF] Maaten | Visualizing data using t-SNE[END_REF]:

q ij = (1 + ||z i -µ j || 2 ) -1 j (1 + ||z i -µ j || 2 ) -1 (2)
and p ij is the target distribution defined as:

p ij = q 2 ij / i q ij j (q 2 ij / i q ij ) (3)
The set of centroids µ is initialized using a k-means on z. The improvement from IDEC [START_REF] Guo | Improved Deep Embedded Clustering with Local Structure Preservation[END_REF] is to keep the decoder and the reconstruction loss L r even after the initialization. The intuition behind it is that the clustering loss, by distorting the embedding space, may alter the representativeness of embedded features and thus the clustering performance, clusters being no more meaningful. Therefore, the clustering loss improves separability of clusters, while the reconstruction loss keeps clusters matched to the features learned by the autoencoder in the first step. The reconstruction loss is computed as the mean squared error between the input time series and the output of the autoencoder. So we have a combined loss that is defined as follow:

L = L r + γ * L c (4)
where γ > 0 is a coefficient that controls the degree of distorting embedded space.

Constraints integration

The extension of DEC to incorporate constraints is based on the Deep Constrained framework (DCC) [START_REF] Zhang | Deep Constrained Clustering-Algorithms and Advances[END_REF]. They propose four types of constraints, but we only took in consideration pairwise constraints, because they are supported by various types of constrained clustering algorithm. The loss function used for must-link constraints set ML is:

l M L = L r -γ M L * (a,b)∈M L log j q aj * q bj (5)
In the same way, the loss function for cannot-link constraints set CL is:

l CL = - (a,b)∈CL log(1 - j q aj * q bj ) (6)
In an intuitive way, ML loss prefers instances with same soft assignments and the CL loss prefers instances with opposite assignments. The ML loss is mitigated by a coefficient γ M L > 0 and the addition of the reconstruction loss L r in a similar way that the clustering loss L c to prevent the method to assign all elements to only one cluster.

Application to satellite image time series

The main purpose of these experiments is to see if this new type of constrained clustering can be applied on satellite image time series and how it competes with the state of the art. We tested the original DCC framework, which is composed of fully connected layers. The network architecture remains unchanged, the time series is only flattened before being fed to the encoder. We also proposed a modified version of DCC with 1D-convolutional layers, as it proved to work well on time series supervised classification [START_REF] Fawaz | Deep learning for time series classification: a review[END_REF]. In this new version, we keep the original input dimension and the network is composed only of 1D-convolutional layers followed each time by a batch normalization layer and the embedding layer is preceded by a global average pooling layer. The embedding layer remains a fully connected layer. We did not include the method in [START_REF] Sai Madiraju | Deep Temporal Clustering: Fully Unsupervised Learning of Time-Domain Features[END_REF] yet, as it relies on the choice of a similarity metric for the target distribution.

4 Experiments and results

Dataset and experimental setup

For these experiments, we apply these methods on crop classification which is an important field of research in remote sensing that has seen numerous study [START_REF] Sicre | Determination of the crop row orientations from Formosat-2 multi-temporal and panchromatic images[END_REF][START_REF] Garnot | Time-Space tradeoff in deep learning models for crop classification on satellite multi-spectral image time series[END_REF]. The dataset is composed of 12 class of different kind of crops (wheat, corn silage, irrigated corn, ect, see Fig. 1c), located near Toulouse (Southwest France). The original set of images3 is composed of 11 multispectral (green, red, and near-infrared) 1000 × 1000 pixel images non-uniformly sampled from 15/02/07-20/10/07 and captured by the Formosat-2 satellite. One of the images is presented in Fig. 1a. The dataset is composed of pixels randomly selected within the annotated areas (see Fig. 1b) that were then split into a train and test sets composed of 1974 and 9869 pixel time series respectively. The algorithms are evaluated using the test set. The train set is only used to set hyperparameters for the spectral clustering method (Spec). Constraints are generated from the test set by randomly sampling pairs of pixels and creating an ML or CL constraint depending upon their labels. The reference data is based on farmer's declaration to the EEA's Common Agricultural Policy. To test how methods benefit from constraints, we define three levels of constraints size: 5%, 15% and 50% of the cardinality of the dataset N = 9869 (a very small fraction of the number of possible constraints, 1 2 N [N -1]). For the evaluation we used the Adjusted Rand Index (ARI) and the constraints satisfaction rate (Sat.) averaged over ten runs. For each level of constraints, ten random sets of constraints were generated to be used at each run. This ensures that each method benefits from the same constraints.

Methods compared and parametrization

To have a comparison baseline, we add to the two deep constrained clustering methods, four standards constrained clustering methods. We use a constrained k-means algorithm (COP-KMeans) [START_REF] Wagstaff | Constrained kmeans clustering with background knowledge[END_REF], a spectral method (Spec) [START_REF] Li | Constrained clustering via spectral regularization[END_REF], a declarative method (CPClustering) [START_REF] Duong | Constrained clustering by constraint programming[END_REF] and a collaborative method with 3 k-means agents (SAMARAH) [START_REF] Forestier | Collaborative clustering with background knowledge[END_REF]. To highlight the variability caused by the choice of metric, we use the Euclidean and the DTW metrics [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] Note that CPClustering does not require any initialization parameters. The others need at least the number of clusters, otherwise default parameters were used. The only exception is the Spec method which needs some hyperparameters that are defined by grid search on a the train set. For deep clustering, we follow the settings in DEC [START_REF] Xie | Unsupervised deep embedding for clustering analysis[END_REF] and IDEC [START_REF] Guo | Improved Deep Embedded Clustering with Local Structure Preservation[END_REF] but as the results were not stable we had to do some minor changes (see section 4.3 for more details). For the layers dimensions we set the embedding layer to a dimension of 2 instead of 10, because it seems to give more stability. For the DCC, the encoder network is set to dimensions d -500 -500 -2000 -2, where d = l * f and l is the length of the input time series and f the number of features per time step, and for DCC-conv the dimensions are l * t -128 -256 -128 -2, the corresponding 1D-filters have a dimension of 8 -5 -3, following Wang et al. [START_REF] Wang | Time series classification from scratch with deep neural networks: A strong baseline[END_REF]. For each of them the decoder is a mirror of the encoder. For the optimizer both of them are trained using a SGD with momentum of 0.9 and a decay value of 1e -6, to compensate the variability mentioned before. γ and γ M L are both set to 0.1 as described in the original papers, these values should be set as small as the learning rate used is high, otherwise the effect of the distortion is too important if both are high, or negligible if both are low. In our case, the learning rate is set high (0.1) so the clustering loss should have less weight (0.1) than the reconstruction loss, more explanations can be found in [START_REF] Guo | Improved Deep Embedded Clustering with Local Structure Preservation[END_REF]. 5

Results

The results with and without constraints are presented in Table 1. Spec gives the best overall results, but that must be mitigated, as mentioned before, by the fact that it needs a training for its hyperparameters (for unconstrained its average results is at 0.367 against the retained configuration at 0.737). We can also point out that the result depends on the chosen similarity metric. This is also the case for the other metric-based methods at the exception of COP-KMeans which performs badly in both cases. The deep clustering only proves to be more efficient when the number of constraints are very high and only with the convolutionnal architecture. But the most intriguing point is the partially opposite behavior of the two architectures. DCC gives relatively good results in unconstrained configuration, but the constraints have a strong negative effect. This effect has already been studied in constrained clustering, and can be observed for all other methods at the exception of SAMARAH. It was observed in Lampert et al. [START_REF] Lampert | Constrained distance based clustering for time-series: a comparative and experimental study[END_REF] that if the algorithm already captures well the structure of the data it will not benefit, or even have a negative effect, from constraints addition. This seems to be the case for DCC. This goes in an opposite conclusion than Zhang et al. [START_REF] Zhang | Deep Constrained Clustering-Algorithms and Advances[END_REF], which concludes to no negative effect. In our case this can be explained by the noisy ground truth used (presence of trees or road across the field, problem of frontier pixels).

DCC-Conv, in the other hand, obtains good results with constraints, but only if the number of constraints is high enough. Indeed the way constraints are used, in the backpropagation, does not enforce the algorithm to respect constraints, and needs a lot of information to be sure that it propagates to the weights. On this point DCC-Conv really shows its specificity, it was observed in Lampert et al. [START_REF] Lampert | Constrained distance based clustering for time-series: a comparative and experimental study[END_REF] that methods does not necessarily benefits from a higher number of constraints, but mostly from informative and coherent constraints. In this case, the ARI seems linked to the number of constraints which may let think that the network starts to learn the dataset itself and not the structure, the algorithm being trained and learn on the same set.

Two other points have to be highlighted, that are not displayed in Tab 1. First the standard deviation strongly increases with constraints and stabilizes again when the number of constraints is higher, for DCC and DCC-Conv (i.e. for DCC-Conc, standard deviation is respectively for unconstrained, 5%, 15%, 50% at 0.005, 0.069, 0.010 and 0.015). This seems to show that the quality of constraints is important, the impact of noise might therefore be smoothed by the number of constraints. The second one is that the network is not stable during the training, this is the case for both architectures but in a higher amplitude for DCC. This tends to stabilize, due to the decay added in the optimizer, but the variability remains (i.e. for DCC, ARI goes from 0.55 to 0.74 to go down to 0.45). An illustration of this instability can be seen in figure 2. In a first step, the ARI increase smoothly, but then it starts to oscillate and struggle to converge again. This is mostly the case for constrained runs and but also, in a smaller amplitude for unconstrained runs. This seems to come mainly from the fact that the target distribution (the hard assignment) is updated regularly, so the objective function does not aim the same target.

Conclusion

The deep clustering demonstrates that it can achieve good results on remote sensing time series without and with constraints. Moreover, it does not require to choose a representation method as it is learned by the autoencoder, which makes it easier to handle for a domain expert. However the role played by the hyperparameters (i.e. dimensions of layer, especially the embedding layer, optimizer) needs some further investigations, as they seem to influence the output quality. There are also two main problems that we plan to study. First, DCC is not as robust as expected from the result of the previous studies. Second, the instability in the training and in the impact of constraints lowers the average result. We plan to study which factors may induce these problems, and further investigate the effect of noisy constraints, the size of the dataset and how constraints can be differently integrated. Finally, we also plan to study how the learned network behave on transfer learning.

4

 4

  .

Fig. 1 .

 1 Fig. 1. An image from the time series: 12 classes, and 11 time points (t4 displayed here).

Fig. 2 .

 2 Fig. 2. Evolution of ARI measure through the training process compared between DCC and DCC-Conv.

Table 1 .

 1 Unconstrained and Constrained ARI and constraint satisfaction. The best performances for each measure, constraint fraction, and distance measure are highlighted in bold. Unconstrained satisfaction was measured using the 50% constraint sets.

	Method Distance Unconstrained	5%	15%	50%
		ARI	Sat.	ARI Sat. ARI Sat. ARI Sat.
	COP-KMeans [14]	DTW 0.426 0.812 0.416 1.00 0.407 1.00 0.436 1.00 Eucl. 0.420 0.807 0.406 1.00 0.443 1.00 0.369 1.00
	Spec [15]	DTW 0.531 0.840 0.683 0.867 0.725 0.888 0.786 0.911 Eucl. 0.737 0.885 0.671 0.854 0.702 0.875 0.781 0.916
	CPClustering [16]	DTW 0.437 0.803 0.469 1.00 0.510 1.00 0.589 1.00 Eucl. 0.681 0.413 0.650 1.00 0.542 1.00 0.510 1.00
	SAMARAH [23]	DTW 0.406 0.802 0.597 0.870 0.637 0.867 0.681 0.878 Eucl. 0.463 0.817 0.691 0.884 0.714 0.890 0.702 0.885
	DCC [13]	0.703 0.885 0.550 0.852 0.448 0.816 0.615 0.862
	DCC-Conv	0.508 0.833 0.497 0.844 0.491 0.819 0.820 0.936

This is the author's version of an article published at ECML/PKDD Workshop MACLEAN 2019.

Provided by the Centre d Études Spatiales de la Biosphère (CESBIO) Unité Mixte de Recherche CNES-CNRS-IRD-UPS, Toulouse, France.

Implementations for compared method available at from https://icube-forge. unistra.fr/lampert/TSCC

The source code used for this paper: https://github.com/blafabregue/ DeepConstrainedClustering

Acknowledgments. This work was supported by the ANR TIMES project (grant ANR-17-CE23-0015) of the French Agence Nationale de la Recherche.


A Durand Petiteville

S Hutchinson

V Cadenat

M Courdesses

2D visual servoing for a long range navigation in a cluttered environment

In this paper, we consider the problem of realizing a vision-based long range navigation in a cluttered environment. To this aim, we couple a topological environment representation and a supervision algorithm. The latter organizes several elementary tasks which allow to perform image based visual servoing, obstacle avoidance and to estimate on line several required parameters. Well known and new methods based on the tasks function approach are introduced.The dynamical sequencing formalism is used to guarantee the global control law smoothness. Finally, simulation results validate our approach.

I. INTRODUCTION

We address the problem of the navigation of a nonholonomic mobile robot based on visual informations. In the survey [START_REF] Bonin-Font | Visual navigation for mobile robots : a survey[END_REF], there exists two main approaches: map-based navigation and mapless navigation. The first method is subdivided into: metric map-using navigation systems (MMUNS), metric map-building navigation systems (MMBNS) and topological map-based systems (TMBS). MMUNS require a complete map of the environment before the navigation starts. MMBNS build the map by themselves then use it during the navigation. Finally, TMBS build and/or use topological maps which consist of nodes linked by lines where each node represents the most characteristic places of the environment [START_REF] Gaspar | Vision-based navigation and environmental representations with an omnidirectional camera[END_REF]. Mapless navigation mostly include reactive techniques that use visual clues derived from the segmentation of an image [Chaumette and Hutchinson, 2006] or optical flow [START_REF] Santos-Victor | Visual-based obstacle detection: a purposive approach using the normal flow[END_REF].

We aim to realize a long range navigation in a cluttered environment using the minimal a priori and a posteriori knowledges. The MMUNS and MMBNS need and/or provide a too rich information to be considered in our case. Several appropriate approaches based on reactive techniques coupled with a topological map have been proposed. In [START_REF] Blanc | Indoor navigation of a wheeled mobile robot along visual routes[END_REF], an indoor teleoperated prenavigation is realized in order to record the image stream. Then, a visual route, comparable to a topological map, is built thanks to the key frames extracted from the stream. The key frames correspond to robot close configurations to ensure the local visual servoing stability. Finally, to reach the goal the robot has to navigate with the successive key frames. Thus, the path followed during the pre-navigation step is replayed thanks to the images. A Durand Petiteville, Cadenat and Courdesses are with CNRS; LAAS : 7 avenue du colonel Roche, F-31077 Toulouse, France; Université de Toulouse ; UPS, INSA, INP, ISAE ; LAAS ; F-31077 Toulouse, France.

[adurandp, cadenat, courdess]@laas.fr S. Hutchinson is a professor of Electrical and Computer Engineering at the univerrsity of Illinois, USA. seth@illinois.edu similar approach is proposed in [START_REF] Kralnìk | A simple visual navigation system with convergence property[END_REF] in an outdoor environment. Nevertheless, the obstacle avoidance is not taken into account in these methods. In [START_REF] Cherubini | A redundancy-based approach to obstacle avoidance applied to mobile robot navigation[END_REF], a controller allowing the robot to perform an obstacle avoidance is proposed. The previous methods involve two main drawbacks. First, numerous key frames compose the image data base to perform a long range navigation. Second, the navigation method requires to develop new necessary controllers, such as the obstacle avoidance in [START_REF] Cherubini | A redundancy-based approach to obstacle avoidance applied to mobile robot navigation[END_REF], to perform a safe navigation.

The image based visual servoing presented in [Chaumette and Hutchinson, 2006] is a well known and studied technique. It consists in controlling the robot to make the image features converge to the desired ones. These latter, which correspond to a robot configuration with respect to a landmark, can be obtained either during a pre-navigation work or estimated on-line [Durand Petiteville et al., 2010]. In addition to the controller realizing the navigation, several methods have been developed to avoid obstacles and to handle occlusions [START_REF] Folio | Folio and Cadenat[END_REF]. Nevertheless, image based visual servoing is usually used to perform short range navigation. Indeed, it requires to see the considered target to compute the control law during all the navigation and especially at the beginning. This condition can be unsatisfied because of the presence of occulting obstacles. If the target to reach cannot be seen at the inital robot position, a solution consists in dividing the global navigation mission into elementary navigation tasks with respect to a set of landmarks. Each intermediate target can be seen from the neighborhood of the previous one. So, in order to extend the navigation range, we propose to add a topological map whose each node corresponds to a landmark. Moreover, a supervision algorithm is developed to organize the tasks. Finally, each task is performed using the task function approach [START_REF] Samson | Robot control : The task function approach[END_REF] and the transition between them are managed using the dynamical sequencing formalism proposed in [START_REF] Souères | Dynamical sequence of multi-sensor based tasks for mobile robots navigation[END_REF].

The paper is organized as follows. Section 2 is dedicated to the topological map and supervision algorithm description. Section 3 details the several controllers computation. Finally, simulation results validate the proposed approach.

II. CONTROL STRATEGY

In this paper, we propose to realize a long range navigation based on a visual reactive controller and a topological map. The considered robot is a cart-like vehicle equipped with a camera mounted on a pan-platform and a laser. The system is detailed in III-A.

A. The topological map

We propose to model the environment by a set T G = (T 1 , T 2 , ..., T n ) of n landmarks. For an indoor navigation different kinds of landmarks can be used : corridors, doors or wall edges [START_REF] Vassalo | Visual servoing and appearance for navigation[END_REF]. Now, we define a graph G (Fig 1 ) where each node corresponds to a landmark T i (for i ∈ [1, ..., n]). A node T i is linked to a node T i+1 if T i+1 can be seen from the neighborhood of T i . If the target to reach cannot be seen from the initial robot configuration, the graph G is used to provide a set T P of successive targets allowing to perform the whole navigation. T P only contains m ≤ n elements of T G . To compute the set T P , we propose to use the Dijkstra's algorithm. For example, in figure 1, the environment is depicted by the whole graph. The computed path allowing to reach T 7 from T 1 is represented in grey. In this case

T P = [T P1 , T P2 , T P3 , T P4 , T P5 ] = [T 1 , T 2 , T 4 , T 6 , T 7 ].
Thanks to the topological map, the set of landmarks to reach in order to perform a long range navigation is known.

B. The global strategy

To perform a long range navigation in a cluttered environment guaranteeing the robot safety, several tools, such as visual servoing, obstacle avoidance and visual features estimation, are available. These tools are coupled to define a set of useful tasks, such as the visual servoing using estimated visual features to handle occlusion during the navigation. The tasks are placed in an supervision algorithm (Fig 2) in which the events involving a transition between two tasks are described. Thus the task to perform is defined using the robot perception of the environment and the knowledge of the previous task.

The different tasks have to be defined in order to allow the robot to reach successively the targets included in T P while preserving the system and navigation integrity. The main task [B] is the visual servoing which is performed thanks to a reactive controller allowing the current image data, denoted by s j , to converge towards the reference ones, noted s * j (for j ∈ [1, ..., m]). The latter corresponds to the sensor data values at the desired system configuration with respect to T P j . During the navigation, if an obstacle is detected in the robot neighbourhood, the algorithm switches to the Once the value of s * 1 has converged, the visual servoing can be started.

For the following visual servoing, the next target T P j+1 to reach needs to be found by the camera during the navigation wrt T P j . To this aim, during the tasks [F] and [G], the mobile base is controlled to perform a visual servoing or an obstacle avoidance, based on the estimated visual data as in [D] and [E], whereas the camera is looking for the next target by scanning the environment from one of its mechanical bound to the another one. If the target T P j+1 is not found at the end of the scan, the visual servoing wrt the current target T P j or the obstacle avoidance restarts. Otherwise the camera is controlled to keep the target T P j+1 in its field on view. Then, the obtained visual data are used to estimate s * j+1 . If the target is not found at the end of the visual servoing task, the map is updated and the planning task is restarted.

Once the convergence of s * j+1 is detected, T P j+1 becomes the task to realized and task [H] is launched. This task is defined in section III and allows to guarantee the continuity of the control law when switching between two visual servoing with respect to two different targets.

By applying this algorithm for the m targets included in T P , the robot is able to reach the final target which cannot be seen by the camera at the initial position. Moreover, the robot safety and the navigation success are guaranteed thanks to the use of several tasks allowing the system to overcome the problems of occlusion and collision. However, the continuity problems between each task have to be addressed during the controllers definition.

III. CONTROL ASPECTS

In this section, a set of controllers allowing to perform each task of the supervision algorithm are presented. We choose to compute them thanks to the task function approach in order to use the dynamical sequencing formalism. Thus, the continuity between two tasks is guaranteed. Moreover the different estimation processes are introduced.

A. System modelling

The considered system is a cart-like robot equipped with a camera mounted on a pan-platform. Moreover a laser has been added. Fig 3(a) introduces the robot model, which requires to define the successive frames : F O attached to the world, F M linked to the robot, F P attached to the platform, and F C linked to the camera. Let θ be the direction of the robot wrt. x O , ϑ the direction of the pan-platform wrt. x M , P the pan-platform centre of rotation, C x and C y the coordinates of C in F P , and D x the distance between the robot reference point M and P. Defining vector q = (l, θ , ϑ ) T where l is the robot curvilinear abscissa, the control input is given by q = (υ, ω, ϖ) T , where υ and ω are the cart linear and angular velocities, and ϖ is the pan-platform angular velocity wrt. A landmark, which can be characterized by a number k of points P i (for i ∈ [1, ..., k]) in the world frame, is represented by a 2k-dimensional vector s made of the coordinates (X i ,Y i ) in the image plane. For such a robot and visual features, the relation between ṡ and the control inputs q is given by : ṡ = L (s,z) J r q

(1)

where J r is the robot jacobian [START_REF] Folio | Folio and Cadenat[END_REF]]

J r =   -sin(ϑ (t)) D x cos(ϑ (t)) +C x C x cos(ϑ (t)) D x sin(ϑ (t)) -C y -C y 0 -1 -1   (2)
and L (s,z) = [L T (P 1 ) , ..., L T (P k ) ] T the interaction matrix. L (P i ) is classically given by [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]:

L (P i ) = L x (s i , z i ) L y (s i , z i ) = 0 X i z i X i Y i f -f z i Y i z i f + Y 2 i f (3)
B. Controller sequencing using the task function approach

In this section we propose to define the various tasks necessary for the long range navigation, using the task function approach [START_REF] Samson | Robot control : The task function approach[END_REF]. With this method, each task is defined by a vector e, called task function, to vanish. To ensure the continuity between e i the current task and e i+1 the next one, we suggest to use the dynamical sequencing formalism [START_REF] Souères | Dynamical sequence of multi-sensor based tasks for mobile robots navigation[END_REF]. To ensure the continuity, the tasks have to be admissible at the switching time t s . A task is admissible if its jacobian can be invertible. For a transition, the control law expresses as follows:

q(t) = qi (t) ∀t ≤ t s qi+1 (t) ∀t ≥ t s (4)
where qi (t) and qi+1 (t) are respectively the control inputs allowing to make e i and e i+1 vanish. To guarantee the control law smoothness between these controllers, [START_REF] Souères | Dynamical sequence of multi-sensor based tasks for mobile robots navigation[END_REF] propose to use the dynamical sequencing formalism. The control law is defined as:

qi+1 (t) = J -1 i+1 ė * i+1 (t) (5)
where ė * i+1 (t) is computed to ensure the global control input continuity at the first order. We have chosen the dynamic presented in [START_REF] Mansard | Task sequencing for sensor-based control[END_REF], which allows to set independently the error decreasing speed and the transition time duration. So we obtain :

ė * i+1 (t) = -λ e i+1 (t) + ρ(t) (6)
where

ρ(t) = [J i+1 J -1 i ėi (t s ) + λ e i+1 (t s )] exp -τ(t-t s ) (7)
with λ > 0 and τ > 0. Choosing ėi (t s ) = J i qi (t s ) ensures the control law smoothness at t = t s , as the dynamics of the task function e i+1 at t = t s depends on the control law qi (t).

C. The controllers

In this section, the necessary controllers to perform the different previously introduced tasks are presented. The visual servoing and obstacle avoidance controllers come from former works, whereas the reorientation one is developed for the first time.

1) Visual servoing: The vision-based navigation task consists in positioning the camera with respect to a given static landmark using visual data. The value of the visual features s depends on the relative camera position with respect to the landmark. So the goal of the image based visual servoing is to make the current visual signals s converge to their reference values s * . s * then corresponds to the value of s obtained at the desired camera pose with respect to the target.

To perform the desired vision-based task, we apply the visual servoing technique given in [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] to mobile robots as in [START_REF]Applying visual servoing techniques to control a mobile handeye system[END_REF]. The proposed approach consists in expressing the visual servoing task by the following task function to be regulated to zero:

e vs = C(s -s * ) (8)
where matrix C, called combination matrix, allows to take into account more visual features than available degrees of freedom. For C = L + (s * ,z * ) , the task jacobian is equal to J V S = CL (s * ,z * ) J r = J r . The task is admissible as det(J r ) = D x = 0.

By imposing an exponential decrease on e, a controller making e vanish is obtained in [START_REF] Folio | Folio and Cadenat[END_REF] :

qvs = -(CL (s,z) J r ) -1 λ vs C(s -s * ) (9)
where λ vs is a positive scalar or a positive definite matrix.

2) Obstacles avoidance: To perform a vision based navigation task in a cluttered environment, we have to deal with the problems of collision. In this section we first present the mobile base controller allowing to avoid obstacles, then the pan-platform one to perform the target tracking. The chosen strategy consists in making the vehicle avoid the obstacle by following the security envelope ξ 0 as shown in Fig. 4. With this method, only the mobile base (i.e υ and ω) is taken into account. To design the desired controller, the following task function has been defined [START_REF] Folio | A controller to avoid both occlusions and obstacles during a vision-based navigation task in a cluttered environment[END_REF]: where Y 0 is the abscissa of the visual pattern gravity center. We can now obtain a task function e oa = (e mb e pp ) T for the obstacle avoidance. Using classical visual servoing techniques, a controller making e oa vanish is obtained :

e mb = (l -v r t δ + kα) T (
qoa = J -1 oa (-λ oa e oa -A oa ) (12)
where λ oa is a positive scalar, A oa = ( -v r 0 0 ) T and:

J oa =   1 0 0 sin(α) -kχ cos(α) k 0 L Y 0 J r   (13) We define χ = 1 R 1+ σ
R δ , with R the curvature radius of the obstacle and σ = {-1, 0, 1} depending on the sense of the robot motion around the obstacle. Moreover,

L Y 0 = (-f z 0 , Y 0 z 0 , f + Y 2 0 f
), and z 0 is the depth of the visual pattern gravity center. It should be noticed that the obstacle avoidance task is admissible as the task function jacobian J oa is invertible:

det(J oa ) = -κ C x z 0 + C y Y 0 z 0 + 1 +Y 2 0 = 0 (14)
3) The reorientation phase: The classical image based visual servoing presented in III-C.1 allows to compute the camera kinematic screw. Some of the approaches proposed in [Chaumette and Hutchinson, 2006] allow to decouple the camera velocities. Thus, if the camera trajectory appears to be correct, it is not necessary the case of the one of the mobile base. For example, as shown in figure 5, where q init = [ 0 -2π 3 5π 6 ] T , the system performs the task whereas this global behaviour is incongruous. To overcome this problem we propose to define a controller allowing to turn toward the mobile base toward the next target without performing a reverse movement. Moreover, the camera has to be controlled to keep the target in its field of view. To this aim, we define the following task function :

e tr = (l tr -v tr t θ tr -ω tr t Y 0 ) T (15)
where v tr and ω tr are the reference linear and angular velocities, l tr is the curvilinear abscissa of point M and θ tr the robot heading. v tr is chosen > 0 to preserve the robot of a reverse movement while ω tr is the same sign as ϑ (t s ), the pan-platform orientation in the robot frame at the switching time. Indeed, before the transition phase, the camera is controlled to preserve the target visibility. So, the pan-platform orientation at t s informs us about the suitable mobile base heading. Controlling the mobile base to regulate to zero the camera orientation means directing the former towards the new target. Making e tr vanish implies that the robot linear and angular velocities are the right ones, while the camera is looking to the landmark. We propose to impose an exponential decrease: ėtr = -λ tr e tr (16)

The expression of ė is now required:

ėtr = (v -v tr ω -ω tr Ẏ0 ) T = J tr qtr + A tr (17)
where

J tr =   1 0 0 0 1 0 L Y 0 J r   A tr =   -v tr -ω tr 0   (18)
where L Y 0 is given in III-C.2. It shoud be noted that:

det(J tr ) = -f (D x cos(ϑ ) -C x ) +Y 0 (D x sin(ϑ ) -C y ) z 0 -f - Y 2 0 f = 0 (19)
Now, thanks to ( 15), ( 16) and ( 17), we obtain the controller to realize a transition:

qtr = J -1 tr (-λ tr e tr -A tr ) (20)
The transition phase is considered to be finished when ϑ = 0. At this moment the mobile base is oriented toward the target.

D. Parameters estimation for visual servoing

Several parameters are needed to be able to compute the control inputs (9) : the current visual features s, the desired visual features s * and the depth z i of each point p i which characterizes the landmark. In the following sections, the estimation process for each parameter is presented.

1) The current visual features: The controller (9) can only be used if the visual data are available. If they are not the task cannot be realized anymore. To remedy this critical situation, [START_REF] Folio | Folio and Cadenat[END_REF] has recently proposed to compute the visual data evolution for any t ∈ [t k-1 ,t k ], with t k = t k-1 + T s and T s the sampling period. To this aim, the values of X i (k),Y i (k) and z i (k) are computed using X i (k -1),Y i (k -1), z i (k -1) and q(k -1). This method requires the initial values of X i (t occ ),Y i (t occ ) and z i (t occ ), where t occ correspond to the occlusion time. The two first ones can be easily obtained thanks to the last image before occlusion, while the depth has to be estimated.

This method is used to realize the visual servoing [D] or the obstacle avoidance [E] tasks using the estimated visual features.

2) The visual feature depth: To provide the necessary depths for III-D.1, the estimator proposed in [Durand Petiteville et al., 2010] is used. The estimation process is based on a predictor/corrector pair using a number m of previous images. In the proposed algorithm, a high number of data is used to improve the depth estimation. Thus, the estimation process is robust to the errors due to the image process and odometry.

3) The desired visual features: In [Durand [START_REF] Petiteville | [END_REF], a method is proposed to automatically compute the desired visual features. Thanks to the estimated depth of the visual data and the images, the landmark dimensions are computed. Then a desired camera position with respect to the landmark is defined. Thus, using the camera pinhole model, the desired visual features can be obtained.

IV. SIMULATIONS

We have simulated a long range navigation in a cluttered environment using the software Matlab TM . The scene is composed of four occulting obstacles, a non-occulting one and a set of nine landmarks (see figure 6 ). The landmarks are made of four circles, whose centers C k are extracted by an image processing. So, for a target T j ,

s j = [X 1 ,Y 1 , X 2 ,Y 2 , X 3 ,Y 3 , X 4 ,Y 4 ]
T where X k and Y k are the centers C k coordinates. The topological map described in figure 1 is used to model the navigation environment. In this simulation the robot has to reach the target T 7 from its initial configuration q init = [1, 1, 1.81, 0] T . As T 7 cannot be seen from q init because of an occulting obstacle, the set T P = [T 1 , T 2 , T 4 , T 6 , T 7 ] is computed using the Dijkstra's algorithm, as shown in figure 1. In order to be closer of the real conditions, a one pixel noise has been added to the visual data whereas a 3% error has been introduced in the odometry process. Finally the sampling time is T s = 0.1s

As shown in figure 6, the long range navigation (about 20 meters) has been successfully performed. The robot safety is guaranteed and the different targets are sequenced until the goal is reached. It should be noticed that each task shown in the supervision algorithm has been realized at least one time. Indeed, in addition to the classical [B] and [C] tasks, the reorientation task [H] is executed when the considered target changes, unless an obstacle avoidance is performed. Thus, the robot behaviour presented in III-C.3 is not found during the navigation. From a global point of view, the long range navigation is successfully and automatically performed using the unique topological map knowledge. Indeed at the end of the navigation, the camera has reached the desired final position q f inal = [3.9, 7, π] T . However, the trajectory made by the robot is not the shortest one. As the path is realized using only local and reactive techniques, the global trajectory is not improved. Moreover, it should be noticed that the system does not switch to the next landmark when it can be seen but only when it has been detected by the camera. The search task cannot be improved because the topological map does not contain any informations about the landmarks positions. Thus, it is impossible to privilege any direction of research for the next target.

V. CONCLUSION AND FUTURE WORKS

In this paper we have proposed a method allowing to perform a vision based long range navigation. The latter relies on a supervision algorithm a set of elementary tasks coupled to a topological map. This map is made of a small number of landmarks, which is the only knowledge available about the environment. The environment is then represented using a reduced database contrary to the other approaches of the literature. The supervision algorithm allows to select the right elementary task at the right instant. Each of these tasks is performed thanks to a local controller, based either on real or estimated sensory data. In this way, it becomes possible to perform long range navigation mission using well-known reactive techniques.

We will validate the proposed approach using artificial landmarks on our robot Pekee II, developed by Wany Robotics. Finally, it will be also necessary to consider the topological map construction. The use of advanced pattern recognition and image processing techniques will allow to take into account a wider variety of landmarks. In this way, it will become possible to perform a vision based long range navigation in a natural indoor environment.

Fig. 1 .

 1 Fig. 1. The topological map

Fig. 2 .

 2 Fig. 2. Supervision algorithm

  Fig. 3. The robotic system The pinhole model is used to represent the camera (Fig 3(b)). In this case, a point P, with coordinates (X,Y ) in the image plane, is the projection of a point p, with coordinates (x, y, z) in F O . Moreover we define z as the depth in F C of the projected point p and f as the camera focal.A landmark, which can be characterized by a number k of points P i (for i ∈ [1, ..., k]) in the world frame, is represented by a 2k-dimensional vector s made of the coordinates (X i ,Y i ) in the image plane. For such a robot and visual features, the relation between ṡ and the control inputs q is given by :

Fig

  Fig. 4. Collision detection

  10) where l is the curvilinear abscissa of point M (Fig 4), k a positive gain to be fixed and v r the desired linear velocity. The first component of this task function allows to regulate the linear velocity to v r . The second component can be seen as a sliding variable whose regulation to zero makes both δ and α vanish (Fig 4). Thanks to the task (10), only the mobile base is controlled. To keep the target in the camera field of view during the avoidance phase, we have to define a specific task function for the pan-platform : e pp = Y 0 (11)

Fig. 5 .

 5 Fig. 5. Example of incongruous behaviour during a visual servoing

Fig. 6 .

 6 Fig.6. The robot trajectory

Fig. 7 .

 7 Fig. 7. Estimated and real reference visual data with respect to T 7