G Amiard-Hudebine

G Tison

E Freysz
email: eric.freysz@u-bordeaux.fr

Study of a compact nanosecond laser system for the ignition of aeronautic combustion engines

We have studied and developed a compact nanosecond laser system dedicated to the ignition of aeronautic combustion engines. This system is based on a nanosecond microchip laser delivering 6 µJ nanosecond pulses, which are amplified in two successive stages. The first stage is based on an Ytterbium doped fiber amplifier (YDFA) working in a quasi-continuous-wave (QCW) regime. Pumped at a 1 kHz repetition rate, it delivers TEM 00 and linearly polarized nanosecond pulses centered at 1064 nm with energies up to 350 µJ. These results are in very good agreement with the model we specially designed for a pulsed QCW pump regime. The second amplification stage is based on a compact Nd:YAG double-pass amplifier pumped by a 400 W peak power QCW diode centered at =808 nm and coupled to a 800 µm core multimode fiber. At a 10 Hz repetition rate, this system amplifies the pulse delivered by the YDFA up to 11 mJ while preserving its beam profile, polarization ratio, and pulse duration. Finally, we demonstrate that this compact nanosecond system can ignite an experimental combustion chamber.

Introduction

Nanosecond laser-induced spark ignition of combustion engines is usually based on nonresonant gas breakdown [START_REF] Morsy | Review and recent developments of laser ignition for internal combustion engines applications[END_REF][START_REF] Dale | Application of high energy ignition systems to engines Progress[END_REF][START_REF] Weinrotter | Laser ignition of engines[END_REF][START_REF] O'briant | Review: Laser ignition for aerospace propulsion[END_REF][START_REF] Phuoc | Laser-induced spark ignition fundamental and applications[END_REF]. The breakdown process is initiated by few free electrons in the interaction region. For nanosecond laser pulses breakdown in air, multi-photon ionization of air molecules is initially responsible for the generation of these seed electrons. The optical breakdown process is subsequent to avalanche ionization which results in a cascade growth of electron density in the interaction volume. Laser beam energy is absorbed via an inverse bremsstrahlung process and hot and dense plasma is created. The plasma generation and expansion is started and continued within the pulse duration time. The breakdown threshold power density, which depends on many parameters and thermodynamic properties, has been evaluated between and in air [START_REF] Tambay | Laser induced breakdown study of laboratory air at 0[END_REF]. Nowadays, this power density is relatively modest and can be delivered by commercial laser systems and a nanosecond pulse in near infrared with peak power ranging from 0.1 to 1 MW can easily reach this threshold.

However, inducing laser plasma within a combustion chamber does not ensure its ignition. The plasma has to bring a minimum energy to sustain the initial combustion within the ignition kernel [START_REF] Lorenz | Pulse train ignition with passively Q-switched laser spark plugs[END_REF][START_REF] Srivastava | Effect of focal size on the laser ignition of compressed natural gas fueled engine[END_REF][START_REF] Srivastava | Effect of laser pulse energy on the laser ignition of compressed natural gas fueled engine[END_REF]. The minimum ignition energy depends upon the nature of the fuel (gas or liquid droplets). Different laser systems have been developed to ignite the gas mixture of an automobile engine [START_REF] Kofler | An innovative solid-state laser for engine ignition[END_REF][START_REF] Sakai | MW peak power single-mode high-brightness passively Q-switched Nd 3+: YAG microchip laser[END_REF][START_REF] Tsunekane | High peak power, passively Q-switched microlaser for ignition of engines[END_REF][START_REF] Pavel | Composite, all-ceramics, high-peak power Nd:YAG/Cr 4+ :YAG monolithic micro-laser with multiple-beam output for engine ignition[END_REF][START_REF] Lorenz | Characterization of energy transfer for passively Q-switched laser ignition[END_REF][START_REF] Manfletti | Laser ignition of a cryogenic thruster using a miniaturised Nd:YAG laser[END_REF][START_REF] Ma | A novel miniaturized passively Q-switched pulse-burst laser for engine ignition[END_REF][START_REF] Ma | Multiple-beam, pulse-burst, passively Q-switched ceramic Nd:YAG laser under micro-lens array pumping[END_REF][START_REF] Dascalu | High-temperature operation of a diode-pumped passively Qswitched Nd:YAG/Cr 4+ :YAG laser[END_REF][START_REF] Pavel | Ignition of an automobile engine by high-peak power Nd:YAG/Cr 4+ :YAG laser-spark devices[END_REF]. However, the environment where the laser source has to be placed is often very hostile (vibrations, high temperature variations, oil projections). For aeronautic applications where the ignition of kerosene liquid droplets is required, we will show that a laser device delivering short (~1 ns) and powerful pulses (~10 mJ), which is a very good spatial beam shape (M 2 ~1.2), can ignite an experimental combustion chamber. To satisfy all the criteria mentioned above, we have developed a two-stage amplifier, which is injected by nanosecond pulses yielded by a microchip laser. The first fiber amplification stage acts as a pre-amplifier and it is a means of transportation for the pre-amplified pulse. The second amplification stage brings the laser pulse energy delivered by the first amplifier above the spark ignition threshold. This latter compact device, built around a Nd:YAG rod, is placed very close to the motor chamber. It can work under relatively large temperature variations and is mechanically robust.

In a compact combustion engine, the transportation of the laser pulses cannot be in free space; therefore, the first amplification stage that we developed is a fully integrated master oscillator power amplifier, which is based on an ytterbium doped fiber amplifier (YDFA)

preserving the polarization of the seeded pulses. It has been shown that YDFA can cover broadband amplification regimes ranging from continuous wave (CW) to femtosecond pulses [START_REF] Jauregui | High-power fiber lasers[END_REF]. The former seeks for the highest average power with narrow spectral linewidth [START_REF] Theeg | Core-pumped singlefrequency fiber amplifier with an output power of 158 W[END_REF], while the latter aims at the highest peak power. Amplification of nanosecond pulses lies in between, targeting the highest energy per pulse. In the past decade and in the nanosecond regime, many research teams have developed YDFA delivering nanosecond pulses with pulse energy over 1 mJ [START_REF] Lago | High-energy temporally shaped nanosecond-pulse master-oscillator power amplifier based on ytterbium-doped single-mode microstructured flexible fiber[END_REF][START_REF] Brooks | 1 mJ energy, 1 MW peak-power, 10 W average power, spectrally narrow, diffraction-limited pulses from a photonic-crystal fiber amplifier[END_REF][START_REF] Teodoro | High-power ns-pulse fiber laser sources for remote sensors[END_REF] and even 25 mJ using Yb-doped large mode area fibers [START_REF] Stutzki | 26 mJ, 130 W Q-switched fiber-laser system with near-diffraction-limited beam quality[END_REF][START_REF]Nonlinear fiber optics[END_REF]. All those results were achieved by cascading multiple YDFA pumped by laser diodes centered at 976 nm and working at a repetition rate around 10 kHz. However, most of these YDFA could not fit the standard engine constraints mainly because, without dedicated temperature regulation, the large temperature changes in the environment of the engine (T~50 K) lead to a large shift in the central wavelength of the pumping laser diodes (~15 nm for T 50 K) and, therefore, to a drastic reduction of the amplified pulse energy. Moreover, YDFA, which delivers powerful nanosecond pulses, also suffers other limitations associated with amplified stimulated emission (ASE) and nonlinear optical effects, such as four-wave mixing (FWM), and stimulated Brillouin and Raman scattering (SBS and SRS) [START_REF] Fève | Four-wave mixing in nanosecond pulsed fiber amplifiers[END_REF]. To limit these side effects, one generally cascades multiple stage amplifiers separated by dedicated filters and increases the amplification repetition rate up to 10 kHz to further filter and delay ASE and nonlinear effects. To further limit SBS, one can also work with a pulse duration below 2 ns [START_REF] Zhang | 1.5-kW Yb-Raman combined nonlinear fiber amplifier at 1120 nm[END_REF]. The YDFA we developed is based on a 1.1 meter single mode 25/250 µm (radius of the core/clad) Yb-doped double clad (DC) large mode area (LMA) polarization maintaining (PM) fiber. It is seeded by nanosecond pulses centered at 1064nm whose duration is =1.00±0.05 ns delivered by a Nd:YAG microlaser. Its repetition rate can be varied by the user from a single pulse up to 10 kHz. To overcome the limitation imposed by the temperature variations sensed by the pumping laser diodes within the engine, this YDFA is pumped at 915 nm. Compared to the absorption cross section of Yb 3+ ions at 976 nm, the absorption cross section of the Yb 3+ ions at this latter wavelength is reduced; however, its broader spectral bandwidth and smoother structure makes it rather less sensitive to temperature variations. This first amplification stage brings the microchip pulse energy from 6 µJ up to 350 µJ. While this system yields nanosecond pulses with no measureable SBS, SRS, or FWM appearing in the spectra of the laser pulse and with a negligible ASE (less than a few percentages of the total spectral density), the delivered output pulse energy is not able to ignite a combustion engine. Therefore, to further boost the pulse energy, we designed a compact rod type Nd:YAG amplifier. The latter was built around a 0.6% doped Nd:YAG rod crystal. The diameter of this 10 mm long crystal is 4 mm. It is pumped by a quasi-CW (QCW) 800 µm fiberized diode laser centered at 808 nm. By carefully designing this amplifier and choosing the pump pulse duration, we were able to bring the energy of the pre-amplified pulse from 350 µJ up to 11 mJ at a 10 Hz repetition rate. The output pulses preserve the polarization, spatial, and spectral profile of the pulses delivered by microchip seed laser. We demonstrate that such laser pulses can ignite a spark in air and ignite an experimental combustion chamber.

As already mentioned, for aerospatial applications, energy saving is a crucial issue.

Therefore, it is important to minimize, as much as possible, the time of use of all the optical amplifiers. In this context, we have specifically studied the transient amplification regime reached during the turn on of the pre-amplifier and the amplifier. Such studies, which are rather important in our context, are often neglected in most of the simulations so far developed [START_REF] Koester | Amplification in a fiber laser[END_REF][START_REF] Giles | Modeling Erbium-doped fiber amplifiers[END_REF][START_REF] Wang | Optimization of pulse amplification in Ytterbium-doped double-clad fiber amplifiers[END_REF].

Therefore, we have developed a simulation that takes into account the transient behavior of the pre-amplifier and the amplifier. The results we obtained are in good agreement with the experiments we performed.

The paper is arranged as follows. In Sec. II, IV, V and VI, we describe the first amplification stage we built and the design we used to avoid the limitations introduced by ASE and nonlinear effects. The experimental data obtained in steady state and transient regimes will then be compared to the model we developed. In Sec. VII, we describe the design of the Nd:YAG amplifier stage that we developed. The evolution of the output pulse energy versus the input pulse energy and the pump pulse energy will then be presented and discussed. In Sec. VIII, we demonstrate the ability of this amplified microchip laser to create plasma in air and to ignite an aeronautic combustion chamber. We conclude by presenting some of the prospects for this device.

II.

Pre-amplifier design and performances.

The YDFA we developed is sketched in Figure 1. It uses a 110 cm long, single mode, 25/250 DC-LMA-PM active fiber doped at 5%. It was pumped with two 25 W 105/125 fibered diodes operating at 915 nm (Oclaro BMU25). The use of two 25 W fiberized pump diodes ensured a large enough amplification of the seeded nanosecond pulses. As already mentioned, the pumping quantum efficiency of Yb 3+ ions is less important at 915 nm compared to 976 nm.

However, the broader and smoother absorption of Yb 3+ ions at 915 nm makes it less sensitive to temperature variations, thereby avoiding the use of an active control of the temperature of the laser diode. This YDFA is seeded by nanosecond pulses yielded by a microlaser (Crylas DPSS 1064-Q2). It delivers linearly polarized pulses centered at the wavelength=1064 nm. At 10 kHz, which is the maximal repetition rate of this microlaser, the pulses had a duration of 1.00±0.05 ns and an energy of 6 µJ. The beam delivered by this microlaser had a divergence smaller than 3 mrad. To prevent any harmful back-reflection, an optical isolator was placed in between the amplifying fiber and the microlaser. To ensure the proper injection of the polarization of the seed pulses along the neutral axis of the amplifying fiber, a half-wave plate was placed in front of the 18 mm aspheric lens, which coupled the seeding pulses within the amplifying fiber. The pulses were coupled into a passive 25/250 PM fiber (DC-LMA-PM). The output of this passive fiber and the two fibered pump diodes were spliced with the input of the active fiber by a 2+1→1 combiner. Due to the small cross section absorption of Yb 3+ ions at 915 nm and the small length of our fiber power, we found that only ~20% of the pump power was absorbed by the amplifier. The output of the active fiber was spliced to an end-cap, which lowered the peak power at the glass-air interface. This end-cap and the input of the passive fiber were cleaved at a large angle to avoid any free lasing. The microlaser and the laser diodes were synchronized by a homemade electronic system that made it possible to remotely adjust the parameters, such as the repetition rate of the microlaser, the duration and peak power of the pump, and the time delay between the pumping and microchip seeding pulses. At the exit of this YDFA, the amplified nanosecond pulse was mixed with the residual pump power and ASE. Most of the remaining pump power was rejected by a dichroic mirror (DM). A diaphragm was used to block the few amounts of pulse energy and ASE propagating in the cladding modes. A polarizing beam-splitter cube transmitted the amplified pulses with the desired polarization. Since the ASE was unpolarized, ~50% of its intensity was reflected by this polarizing cube. The latter was measured and used to correct the effective amplified pulse intensity transmitted by the cube. We pumped the amplifier in the CW regime and recorded the evolution of the amplified output pulse energy versus the pump power. We already mentioned that only ~20% of the pump power was absorbed by the amplifying fiber. The repetition rate of seeding pulses generated by the microlaser was adjusted from 1 kHz to 4 kHz. Whatever the pump power, the far field spatial profile of the output pulse was found to be almost TEM 00 (Fig. 2) and its degree of polarization was always better than 100:1.

III. Continuous pumping regime

The evolution of the output pulse energy versus the pump power absorbed by the amplifying fiber and recorded for different repetition-rates of the microlaser is displayed in We noticed that whatever the repetition rate, and for intensities of the pump power larger than 5W, the amplifier acted as a laser and oscillated at =1030 nm. This phenomenon was due to the fluorescence within the fiber, which was centered at =1030 nm. Above a pump threshold, the fluorescence was strong enough to compensate for the losses at the end of the fiber, and the fiber was able to oscillate at this wavelength. By keeping the pump power below 5 W, the maximal output pulse energy was achieved at 1 kHz repetition rate and was ~185 µJ. As displayed in Figure 3, as the repetition increased the pulse energy decreased. This behavior was also recorded when the pulse repetition rate decreased below 1 kHz. As already mentioned, above 5 W of pump power we noticed the apparition of a free lasing effect. This indicated that during its propagation in the amplifying fiber, the seed pulses did not sufficiently lower the residual gain of amplifier. Indeed, the energy of the seeding pulses was way below the saturation energy of the amplifier medium (E sat =75 µJ). Therefore, the latter was insufficient to extract a large amount of the stored energy in the fiber amplifier. For 4.5 W of pump power, when we increase the repetition rate of the microlaser from 1 kHz to 4 kHz, the output power increased from 170 mW to 480 mW indicating that the pump to seed conversion efficiency increased from ~19% up to ~53%. Simultaneously, we noticed that the amplitude of the spontaneous emission centered at =1030 nm decreased. This clearly stressed that, as expected, at =1064 nm the energy stored in the amplifier could be extracted by increasing the microlaser repetition rate.

However, the energy of the amplified pulses decreased from 170 µJ to 120 µJ. The following sections will show how to further limit the free lasing regime of this amplifier experimentally and numerically, and therefore increase the energy of the amplified pulses.

IV. Quasi-continuous wave pumping

To avoid the apparition of free lasing in our amplifier, we decided to pump it in a QCW regime. In this regime, the amplifier is pumped during a time that is lower or comparable to the fluorescence life time  rad of its long living excited state ( 2 F 5/2  2 F 7/2 ), which, for Yb 3+ ions in a germanosilica glass, is about  rad ~1 ms. New pumping of the amplifying fiber was repeated after an off waiting time of T rad . Hereafter, we report on the influence of the pump duration and the pump power on the output energy of the seeded pulse. Figure 4 displays the evolution of the pulse energy at the output of the amplifier versus the pump pulse duration when the repetition rate of the microlaser is fixed to 1 kHz. We stopped to increase the pump power each time we reached the lasing threshold of the amplifier. As expected, at a given peak power of the pump, as the pulse duration increased, the energy stored in the fiber and the output pulse energy increased. Energies up to 230 µJ were recorded when the amplifier was pumped with 9 W of pump power during 500 µs. In this regime, which corresponded to a duty cycle of 50%, we recorded a slight modification of the temporal envelope of the seed pulse, and the ASE was rather limited (Fig. 5a-b). The pulse duration at half width increased from ~1.00±0.05 ns to ~1.33±0.05 ns. It is important to stress that under this QCW regime, while the output pulse energy increased, the overall efficiency of this amplifier reduced. It was, at most, about 24%. Nevertheless, the beam shape of the amplified pulse remained very good as did its polarization ratio (~100:1). The measurement of the beam diameter for different positions indicated M 2 ~1.1. These measurements, performed with the seeding pulses polarization along the two neutral axes of the fiber, gave very similar results. Figure 5: a) Temporal envelope of the input and amplified pulse: input (amplified) pulse in red (black). b) Spectrum of the output beam. The amplified pulse is centered at =1064 nm. At =1030 nm the apparition of ASE was noticeable. These measurements were performed with 9 W of pump power, and the energy of the amplified pulse was 230 µJ.

Due to the relatively weak absorption (~20%) of the pump power by the amplifying fiber, we did not restrict the length of our amplifying fiber to 110 cm. Using a fiber with the length of L=1.5 m, we amplified the signal to energies up to 350 µJ (Fig. 6). However, occasionally, we noticed that when the amplified pulse energies were ~350 µJ, the generation of nonlinear effects, such as four-wave mixing, was taking place. The latter phenomenon was easily noticeable since it resulted in the generation of visible spectral components. The observation of four-wave mixing regularly foreshadowed the occurrence of damage within the fiber. Therefore, it was safer to keep the length of the amplifier fiber to ~110 cm. Under these conditions, we did not notice the occurrence of any damage within the fiber.

V. Simulations

To simulate the behavior of our fiber amplifier in the QCW regime, we developed a simple model. We considered that, at room temperature, Yb 3+ ions in silica fibers behaved as a quasi-three-level system. However, in silica fiber, the high phonon energy led to the fast decay of the higher excited states to the metastable level. We therefore assumed that the amplifier media could be described by a simplified two-level rate equation [START_REF] Siegman | Laser pumping and population inversion[END_REF]. To account for the behavior of our amplifier, we decomposed its evolution in different sequences. The first was the pumping sequence when the microchip laser was switched off. We therefore derived the population inversion evolution induced by the sole pumping of the amplifying medium. During this sequence, we easily derived the evolution of the pump power at each point within the fiber as well as the amplified stimulated emission. The fiber was split in a 1D grid with meshes of 1 cm.

Smaller meshes did not improve the results. For the mesh at position n, the absorbed pump energy is:

(1) (2)
where is the remaining pump energy, which was not absorbed; ( ) is the absorption (emission) cross section at the pump wavelength (=915 nm); and and are the population density within the fundamental and excited state level. Accordingly, we derived the inversion of population for each mesh following an Eulerian integration:

, ( 3
)
where dt is our time increment and is the energy of the pump photons.

Knowing the excited population, we derived the number of photons emitted by the spontaneous emission SE(n) at each position of the mesh:

For Yb 3+ , this spontaneous emission was mainly centered at =1030 nm. Spontaneous emission occurring over 4 steradians, only the photons emitted in the numerical aperture of our fiber were guided in the core fiber and could contribute efficiently to ASE. The other photons not guided in the core of the fiber or leaking through the cladding were considered lost. To compute the ASE, we assumed that for each mesh, ASE resulted from photons emitted within the numerical aperture of the fiber at mesh n and were added with photons emitted and amplified in the neighboring meshes:

,

where ASE(n) is the number of emitted photons by ASE in the mesh n and is the spontaneous emission at =1030 nm emitted in the numerical aperture of the fiber. According to the experimental results presented in Fig. 5b, we considered that the value of was averaged over the ASE spectral range from =1020 nm to =1040 nm. To compute the overall ASE, we considered photons emitted in the forward direction from mesh 1 to the end of the fiber.

According to Eq. 5, this ASE was overestimated since, doing so, we did not consider the impact on the excited state population of ASE generated in the backward direction. To account for this overestimation, we computed the ASE in the backward direction. We repeated the computation of ASE in the forward and backward directions back and forth. With this method, and after about four iterations, we promptly converged towards a solution, where the difference between ASE between two iterations was less than 1%. The set of equations (1-5) made it possible to compute the population in the excited state N 2 (n) at each point of the fiber and allowed us to evaluate the evolution of the energy of the seed pulse E out (n) within each mesh of the fiber. According to Estable et al. [START_REF] Estable | amplification régénérative et multipassage d'impulsions lumineuses dans des milieux solides (YAG dopé néodyme, alexandrite, saphir dopé titane)[END_REF], E out (n) is:

(6)
where S is the surface of the fiber core; is the emission cross section of ytterbium at 1064 nm; dz is the length of a mesh (1 cm); ( ) is the energy of the seed at the input (output) of the n th mesh with ; and is the saturated energy density which is J sat =h/ . With this latter equation, we were able to compute the energy of the seed pulse at each position of the fiber. Hereafter, we report on the comparisons between our model and our experimental results. As evidenced in Figure 7, theoretical and experimental data are in disagreement. Figure 7 revealed that our simulations were always underestimated by at least a factor of 2 in the experimental results! This meant that our simulations neglected a major phenomenon. As already mentioned, the experimental data were recorded while the repetition rate of the microlaser was fixed at 1 kHz. Therefore, the experimental output pulse energy displayed in Figure 7 was computed according to the average power measured at the exit of the fiber. In practice, this meant that each millisecond, the fiber was pumped during 0.5 ms before the seed pulse was injected. Our simulations were performed considering, first, that this phenomenon took place once; and, second, neglecting the population that remained in the excited state when another pump pulse excited the amplifying fiber. In other words, our simulation neglected the transient phenomena and, therefore, the transient evolution of the population of the Yb 3+ excited state towards its steady state value.

VI. Comparison between simulation and experimental results

To demonstrate that, indeed, such transient states have to be considered, we recorded the evolution of the intensity of a pulse yielded during ~2 ms by the microlaser working at 4 kHz repetition rate while pumping the amplifying fiber with a pump pulse of ~2 ms. As shown in Figure 8, the pulses yielded by the microlaser were increasingly amplified as time went on. In fact, to account for this phenomenon, we also had to consider that both the pump and the input pulse were way below the saturation density energy. Hence, after each amplification, many Yb 3+ ions remained in their excited state. The latter had to be added to the population brought by the pump beam during the following pumping sequence. To account for this transient behavior in our simulation, after the amplification of the microlaser pulse, we introduced the natural decay of the Yb 3+ that remained in the excited state just after the amplification of the injected pulse. This decay stopped immediately upon the application of another QCW pump pulse. The excited population that was computed just prior to this was considered as the new starting point. Accordingly, we could readily compute the transient evolution that should be recorded. Figure 9 shows the plotted evolution of the energy of a 6 µJ nanosecond pulse centered at =1064 nm as it propagates in the amplifying fiber when it is seeded at a 1 kHz repetition rate. During each cycle and prior to seeding, the fiber is QCW pumped during 500 µs by a laser diode centered at =915 nm with a peak power of 9 W. First, we noticed that the steady state regime was reached for the eighth pulse; and second, that, in very good agreement with the discrepancies between the simulation and experimental results displayed in Figure 7, the first seeded pulse was about two times less amplified than the third one (Fig. 9a). In our simulation, we considered that the steady state regime was reached when the energy difference between two consecutive pulses was less than 5%. Figure 9b also showed that a stacking gain effect was computed, considering that the Yb 3+ ions remaining in the excited state emitted at 1064 nm. Indeed, we noted that along the first 40 cm of the fiber, the available gain stocked in the fiber was higher. This observation indicated that the energy of the pulse injected in the fiber was too small compared to the saturation energy; therefore, it was too weak to efficiently empty the population of the Yb 3+ excited state level. Each millisecond, the remaining energy that was stacked up within the fiber was added to the energy brought by the next pumping cycle. Figure 9b also made it possible to compare the energy stocked all along the fiber that was available for the first pulses and in the steady state regime. It clearly indicates that, in the steady state regime, at the entrance (resp. the exit) of the fiber, more energy is stored in the fiber. This latter phenomenon accounted for the increase in energy of the amplified pulse. Figure 9: a) Transient evolution of the energy of a 6 µJ, nanosecond pulse delivered at 1 kHz along a 110 cm long Yb 3+ amplifying fiber. Each millisecond, the fiber is pumped during 500 µs by a 9 W peak power pump pulse centered at 915 nm. The nanosecond pulse is injected at z=0 in the fiber immediately after each pumping pulse. b) Associated evolution of the energy stored the fiber along the 110 cm fiber.

With this transient regime added to our simulations, the experimental and numerical results were in very good agreement (Fig. 10). To further optimize the amplification within our fiber, we used our simulations to study the influence of the pump duration on the evolution of the output energy pulses in the transient regime. In Figure 11, we plotted the evolution of the energy of the first eight pulses exiting the fiber. The fiber was pumped at a repetition rate of 1 kHz by QCW delivered by a laser diode centered at 915 nm where the peak power was 9 W. After each QCW pumping sequence, a 6 µJ nanosecond pulse centered at 1064 nm was injected in the fiber. As expected, the amplified pulse energy increased with the pump pulse duration. The steady state regime was reached after about five pulses. This result was important for the characterization of time ignition of a combustion chamber.

VII. Nd:YAG crystal amplification stage

Once the pre-amplified seed pulse was guided close to the combustion chamber, it needed to be further amplified in order to reach the energy ignition threshold and the breakdown threshold power density. As previously mentioned, the peak power required to induce a plasma ranged between 0.1 and 1 MW. To deliver nanosecond pulses with both adequate peak power and energy, we developed another amplification stage. It had to be inserted in the close vicinity of the combustion chamber; therefore, it would have to deal with relatively large temperature changes and high vibration constraints. To design a compact system, we used a double-pass amplification stage. The chosen amplifying medium was a 10 mm long Nd:YAG rod crystal doped at 0.6% with a diameter of 4 mm.

This crystal that had a lower quantum defect (~0.76) compared to Yb 3+ fiber (~0.89), behaved like a four-level system. Therefore, it was noticeably less sensitive to temperature variations. The radiative lifetime of the exited state was independent of temperature but depended on the concentration of Nd 3+ ions. It was ~230 µs for 1% concentration and decreased to 150 µs for 2% [START_REF] Dong | Temperature-dependent stimulated emission cross section and concentration quenching in highly doped Nd 3+ :YAG crystals[END_REF]. Therefore, in order to store a maximum of energy in the Nd:YAG rod, we chose to fix the concentration to 0.6% and to pump it in a QCW regime. Contrary to its excited state lifetime, the emission cross section of Nd:YAG at =1064 nm-hereafter labeled as  em (1064nm)-decreased as the temperature increased [START_REF] Rapaport | Temperature Dependence of the 1.06 µm stimulated emission cross section of neodymium in YAG and in GSGG[END_REF][START_REF] Sato | Temperature dependencies of stimulated emission cross section for Nd-doped solid-state laser materials[END_REF]. A simple empiric law made it possible to infer the evolution of the emission cross section versus the temperature (Equ. 7):

 em (1064 nm)(T)= (-3.9 10 -22 T + 3.445 10 -19 ) cm -2 (7)

Hence, if the temperature increases from 300 K to 450 K,  em (1064 nm) decreases by ~26%. It was noted that, simultaneously, the absorption cross section of this crystal at =808 nm increased [START_REF] Sato | Temperature dependencies of stimulated emission cross section for Nd-doped solid-state laser materials[END_REF]. Therefore, the latter phenomena could partly compensate for the decrease of  em (1064 nm) with temperature. Moreover, the amplification spectral width of Nd:YAG at =1064 nm was relatively small. This was useful in our system since all the ASE centered at =1030 nm produced in the fiber amplifier could not be further amplified in Nd:YAG.

Furthermore, it should be noted that the peak of emission cross-section band shifted with temperature. The latter phenomenon, which resulted in the decrease of amplification gain, was also partly compensated by the temperature increase of the amplification bandwidth [START_REF] Pavel | Enhancing performances of a passively Q-switched Nd:YAG/Cr 4+ :YAG microlaser with a volume Bragg grating output coupler[END_REF].

Regarding the vibration constraints, it is important to underline that the small dimension of our amplifying crystal increased the vibration modes to higher frequencies and limited their amplitudes to ensure a weaker sensitivity to the latter parameter. Hence this crystal presented many advantages for our particular application. The design of our double-pass amplifier is depicted in Figure 12. At the exit of the fiber amplifier, the beam collimated to a diameter of ~800 µm; it was sent through an optical isolator (OI) and entered the amplification stage where its size was kept constant. To ensure an optimal transmission of the pre-amplified seed pulse, a /2 plate was inserted in between the collimator and the OI. Then, the linearly polarized seed pulse passing through the thin film polarizer (TFP) was amplified. After a first pass in the Nd:YAG crystal, the polarization of the seed pulse was rotated by /2 by means of the /4 plate and was back reflected by a high reflection DM. After its second pass in the Nd:YAG crystal, the beam was reflected by the TFP and was focused in the combustion chamber by an aspheric short focal lens. The DM reflected more than 99% of the signal pulse and transmitted more than 80% of the beam delivered by the pump, a DILAS QCW 800 µm fiberized diode laser centered at 808 nm. The duty cycle of this diode was 2%. To prevent the occurrence of any thermal effect induced by optical pumping of the Nd:YAG rod, we decided to fix the pumping repetition rate at 10 Hz. To optimize the amplification process within the crystal, we had to optimize both the focusing of the pump pulse within the crystal and its spatial overlap with the seed beam. Therefore, good imaging with a very compact optical system of the output of the multimode 800 µm the pump fiber in the Nd:YAG rod was rather difficult but important. The imaging system was computed with ZEMAX. It was made of two antireflection coated aspherical lenses with 8 mm and 13.5 mm focal length. As depicted in Figure 13, the pump beam was gently focused to a minimum size of 920 µm within the crystal.

This pump beam size changed by less than 10% over 90% of the length of the crystal. It was slightly larger than the seed beam size (~800 µm). Figure 15 displays the beam profile and the experimental evolution of the output pulse energy of our double-pass amplifier versus the peak power of the pump pulse when the energy of pre-amplified nanosecond seed pulse was 350 µJ. Here again, we used a home-made synchronization electronic device that insured the Yb 3+ doped fiber and the Nd:YAG rod were pumped during 500 µs and 250 µs, respectively, prior to being injected by the pulse delivered by the microchip laser. This device also made it possible to separately control the repetition rate of the Yb 3+ fiber and the Nd:YAG rod amplifier. In the data displayed in Figure 15 the Yb 3+ fiber and Nd:YAG amplifier were pumped at 1 kHz and 10 Hz, respectively. We noticed that the output pulse energy increased almost linearly versus the peak pump power, which reached up to 11 mJ for 400 W of peak pump power (Fig. 15). This indicated an optical (pump) to optical (seed) conversion efficiency of 11%. The numerical and experimental data presented in Figure 16 displays the evolution of the output pulse energy delivered by the double pass amplifier versus the energy of the pulses delivered by the fiber amplifier when the peak power of the pump beam is fixed at 400 W. We noted that the output pulse energy increased almost linearly with the injected pulses energy. However, as the input pulse energy increased, the factor of amplification in energy decreased from ~50 to ~25. Since the intensity of the pulse was far below the saturation intensity in Nd:YAG (J sat =0.6 J cm -2 ), we attributed this behavior to phenomena occurring in the fiber amplifier when the pulse energy was over 250 µJ, which was related partly to polarization instabilities and partly to the increase of ASE at =1030 nm, which both increased when the energy of the seeding pulses increased. Indeed, as we increased the energy of the injected pulse, the contribution of ASE became less negligible and could not be perfectly corrected. Moreover, the polarization instabilities affected the input energy of the seed pulse, which was transmitted by the thin film polarizer and amplified in the Nd:YAG rod. We also did not record any significant change of the amplified pulses duration that remained about =1.33±0.05 ns. The polarization ratio of the beam was better than 100:1.

This whole set of data indicated that the Nd:YAG amplifier that we designed complied with our experimental requirements: a very compact system that delivered a nanosecond pulse with energy higher than 10 mJ and a beam profile with an M 2 <1.5. It confirmed that the combination of Yb 3+ amplifying fiber and compact double pass Nd:YAG bulk amplifiers was well adapted for our application.

VIII. Laser-induced plasma generation and ignition of a combustion chamber

An interesting review of laser ignition phenomena can be found in the literature [START_REF] Weinrotter | Laser ignition of engines[END_REF][START_REF] O'briant | Review: Laser ignition for aerospace propulsion[END_REF].

Ignition of aeronautic combustion engines is usually realized using conventional spark plugs.

The spark plug is inserted in a metallic jacket and fixed on the combustion chamber. The anode and cathode of the spark plug are located in the near surface of the combustion chamber.

Therefore, when the spark plug is powered, the electrical sparks are produced near the surface of the combustion chamber, which is drilled with hundreds of holes through which flows compressed air. This air flow ensures the efficient cooling of the surface of the combustion chamber as well as the anode and cathode of the spark plug. However, when the combustion chamber is at full power, the air flow near its surface prevents a good interaction between the electric spark and the fuel. Hence, at the full regime, this phenomenon strongly impacts the probability of re-ignition of the engine. This problem can be overcome using laser ignition by focusing the laser pulse a few centimeters away from the surface of the combustion chamber.

Moreover, the combustion chamber can be ignited at positions where optimal ignition conditions may be found. To do so, we removed the conventional spark plug and replaced it with a hollow cylindrical metallic tube on which a sapphire lens was fixed. The diameter of the lens was limited to 8 mm by the diameter of the metallic jacket of the spark plug. Antireflection-coated sapphire lenses with different focal lengths ranging from 20 mm to 7.5 mm were tested. The beam delivered by the Nd:YAG amplifier was directly injected into the metallic tubing containing the focusing lens. For pulse energies larger than 1.5 mJ, we could easily record the creation of plasma in air with all the lenses we tested. The typical plasma generated by a 2.5 mJ nanosecond pulse focused by a 7.5 mm focal lens is displayed in Figure 17. Finally, we tested this amplifying system on a combustion chamber prototype developed at the Office National d'Etudes et de Recherches Aérospatiales (ONERA). The experiment was carried out at an atmospheric pressure using the experimental combustion chamber previously described [START_REF] Linassier | Application of numerical simulations to predict aircraft combustor ignition[END_REF]. The fuel/air mixture of 2.1 was injected during 5 seconds in the chamber. When focusing a 10 mJ pulse produced by our system in the combustion chamber with a 10 mm focal lens, we achieved successful ignition. Under this experimental condition, we found that the probability of ignition was 20% for 10 mJ and it rapidly reached 27% as the pulse energy was increased to 11 mJ. The details of these experiments and the different parameters we explored to improve the ignition of this experimental combustion chamber will be presented in a forthcoming paper.

IX. Conclusion

We studied and developed a compact nanosecond laser system dedicated to the ignition of combustion engines. This system was based on a nanosecond microchip laser that was injected into two successive amplifiers. The first amplification stage was a YDFA working in a QCW mode. At a 1 KHz repetition rate, it delivered TEM 00 and linearly polarized nanosecond pulses with a duration of ~1 ns, centered at 1064 nm with energies up to 350 µJ. These results

were in very good agreement with the model we specially designed for the pulsed-pump regime.

This model was also able to compute the transient state that occurred for the early pulses, which was crucial to enable the design of a lasing ignition device. The second amplification stage was based on a compact Nd:YAG double pass amplifier pumped by a 400 W QCW diode coupled to a 800 µm multimode fiber. At a 10 Hz repetition rate, this system was able to amplify the pulse delivered by the YDFA up to 11 mJ while preserving its beam profile, polarization ratio, and pulse duration.

Finally, we demonstrated that this compact nanosecond system could ignite an experimental combustion chamber. This work offers many prospects that we are now investigating. One is directly related to the increase in the energy of the pulse delivered by the Nd:YAG amplifier as well as the increase in its repetition rate. Another interesting prospect is related to the transient regimes described in this paper. One important problem in the ignition of aeronautic engines is related to the creation of a core combustion that satisfies the Ballal and

Lefebvre criterion [START_REF] Ballal | Flame propagation in heterogeneous mixtures of fuel droplets, fuel vapor and air[END_REF], which indicates that a combustion core can only ignite the chamber when its decrease in temperature, induced during its expansion, is compensated by a temperature increase induced by combustion. Under this criterion, a too-small combustion core cannot survive. To circumvent this problem, the use the picket-fence strategy proposed for ignition of plasma fusion can be envisioned [START_REF] Rothenberg | Ultrafast picket fence pulse trains to enhance frequency conversion of shaped inertial confinement fusion laser pulses[END_REF]. A short and intense pulse creates hot plasma, which is subsequently stabilized and heated by a sequence of less powerful pulses. Prior to the injection of a pulse sequence produced by the microlaser, the fiber amplifier and the Nd:YAG rod need to be pumped. Upon injection of the pulse sequence in the amplifiers, the first pulse within the pulse sequence produced by the microlaser will be strongly amplified while the subsequent pulses will be less amplified, but will take advantage of the remaining population inversion within these amplifiers. While less intense, these laser pulses should be absorbed by the plasma generated by the first laser pulse heat; therefore, they will contribute to its stabilization. This strategy has already been explored for automotive combustion engines by different authors [START_REF] Lorenz | Pulse train ignition with passively Q-switched laser spark plugs[END_REF][START_REF] Dumitrache | Laser generated plasma using a dual pulse approach with application to laser ignition[END_REF].

Figure 1 :

 1 Figure 1: Experimental setup of the YDFA amplifier (DM: dichroic mirror).

Figure 2 :

 2 Figure 2: Spatial profile of the beam amplified in the ytterbium fiber.

Figure 3 Figure 3 :

 33 Figure 3.

Figure 4 :

 4 Figure 4: Evolution of the energy of the amplified pulse versus the peak pump power for different pump pulse durations.

Figure 6 :

 6 Figure 6: Evolution of the energy of the amplified pulse at the exit of a 1.5 meter length fiber versus the peak pump power for a pump duration of 500 µs.

Figure 7 :

 7 Figure 7: Evolution of the energy of the amplified pulse versus the peak pump power and for the different durations of the pump pulse.

16 0Figure 8 :

 168 Figure 8: Temporal spectrogram of 8 amplified pulses injected in the fiber at 4 kHz repetition rate. The amplifying fiber is pumped at the time t=0 during 2 ms with 9 W of peak power.

Figure 10 :

 10 Figure 10: Comparison between simulation and experimental data for different pump durations versus the peak pump power. The energy of the seed pulse is 6 µJ.

Figure 11 :

 11 Figure 11: Evolution of the energy of the amplified versus the duration of the pump pulse for a sequence of eight pulses. The pulses delivered at a repetition rate of 1 kHz are injected immediately after each cycle of pumping.

Figure 12 :

 12 Figure 12: Design of the Nd:YAG amplifier. The seed beam in red (green) is p (s) polarized. /2: half-wave plate; /4: quarter-wave plate; DM: dichroic mirror; OI: optical isolator; TFP: thin film polarizer.

Figure 13 :Figure 15 :

 1315 Figure 13: Simulation of the focusing performed by the two lenses within the Nd:YAG crystal of the pump beam delivered by the 800 µm optical fiber.

Figure 14b and Figure 15 Figure 16 :

 1516 Figure 14b and Figure 15 agree well, indicating that with our simulations we were able to

Figure 17 :

 17 Figure 17: Image of plasma created in the air by a 2.5 mJ pulse focused by a lens of 7.51 mm. The pulse propagates from the right to the left. The plasma is induced near the lens.

Acknowledgments

This project was partly financed by the French National Research Agency (ANR) under the project ECLAIR, the Conseil Régional d'Aquitaine under the project ParaSources. We thank the center Alphanov for the technical help they provided to us at the beginning of this project.

This system was design to be very compact: the distance between the output of the pumping fiber and the DM was ~7 mm. The DM, the quarter-wave plate, and one of the end surfaces of the crystal were almost in contact. To compute the pulse energy after single and double passes in the Nd:YAG crystal, we used Eq. ( 6) where the parameters (J sat , ) of a Nd:YAG crystal doped at 0.6% were considered. Typical results of our simulation are depicted in Figure 14. In this simulation, we considered that 350 µJ pulses delivered by the fiber amplifier were injected in a 10 mm long Nd:YAG crystal with a diameter of 4 mm. Prior to amplification, the latter was considered to be homogenously pumped at 808 nm during 250 µs by a pump contra-propagating with respect to the injected nanosecond pulse. In these latter simulations, the pump repetition rate was fixed at 10 Hz. Due to the fast relaxation time of the excited state (~250 µs), we considered that the steady state regime of this amplifier was immediately reached. To compute the evolution of the population inversion, the pump beam profile computed with software ZEMAX was considered.

The evolution of the inverted population indicated that an inversion of ~45% with 400 W pump power could be reached (Fig. 14a). Figure 14b presents the evolution of the pulse energy during


Joseph Lehec

On the Yao-Yao partition theorem

The Yao-Yao partition theorem states that for any probability measure µ on R n having a density which is continuous and bounded away from 0, it is possible to partition R n into 2 n regions of equal measure for µ in such a way that every affine hyperplane of R n avoids at least one of the regions. We give a constructive proof of this result and extend it to slightly more general measures.

Introduction

In [START_REF] Yao | A general approach to d-dimensional geometric queries[END_REF], Yao and Yao show that for any probability measure µ on R n having a density which is continuous and bounded away from 0, it is possible to partition R n into 2 n regions of equal measure for µ in such a way that every affine hyperplane of R n avoids at least one of the regions. This theorem was designed for computational geometry purposes but it turned out to be useful in other areas of mathematics. For instance, the authors of [START_REF] Alon | Crossing patterns of semi-algebraic sets[END_REF] use it to prove Ramsey type theorems and in [START_REF] Lehec | Partitions and functional Santalo inequalities[END_REF], we show a connection with the Blaschke-Santaló inequality. More precisely, using a partition à la Yao and Yao, we reduce a general functional form of the Blaschke-Santaló inequality (due to Fradelizi and Meyer [START_REF] Fradelizi | Some functional forms of Blaschke-Santaló inequality[END_REF]) to an easy inequality between functions dened on R + . The proof of Yao and Yao is by induction on the dimension and uses the Borsuk-Ulam theorem. The purpose of this paper is to show that the Yao-Yao theorem can be obtained in a much more concrete way, by applying the real intermediate values theorem again and again. Of course this proof is longer, but we believe that it gives a better understanding of the structure of the Yao-Yao partition. Also we are able to get rid of the technical assumptions (µ having a continuous density bounded away from 0) which are annoying for applications. The article [START_REF] Yao | A general approach to d-dimensional geometric queries[END_REF] being very sketchy, we intend to give (almost) every detail.

The paper deals with finite dimensional real affine spaces; if E is such a space, E denotes the associated vector space. We say that P is a partition of E if ∪P = E and if the interiors of two distinct elements of P do not intersect. For instance, with this definition, the set

P = {(-∞, a], [a, +∞)} is a partition of R. Definition 1. If E = {x}
is an affine space of dimension 0, we say that P is a Yao-Yao partition of E if P = {x} and we define the center of P to be x. Let E be an affine space of dimension n ≥ 1. We say that P is a Yao-Yao partition of E if there exists an affine hyperplane F of E, a vector v ∈ E\ F and two Yao-Yao partitions P 1 and P -1 of F having the same center x such that

P = A + R -v | A ∈ P -1 ∪ A + R + v | A ∈ P 1 ,
and we say that x is the center of P.

If E has dimension n, then a Yao-Yao partition of E has 2 n elements and we shall see in the following section that every hyperplane of E avoids at least one of the elements of a Yao-Yao partition. Let us state our main theorem. Let M(E) be the set of non-negative Borel measures µ on E which are finite and satisfy µ(F ) = 0 for any affine hyperplane F .

Definition 2. Let µ ∈ M(E), a Yao-Yao equipartition P for µ is a Yao-Yao partition of E satisfying ∀A ∈ P, µ(A) = 2 -n µ(E). (1)
We say that x ∈ E is a Yao-Yao center of µ if x is the center of a Yao-Yao equipartition for µ.

Theorem 3. Let µ ∈ M(E), there exists a Yao-Yao equipartition for µ. Moreover, if µ has a center of symmetry x ∈ E, then x is a Yao-Yao center for µ.

Main properties

If A is a subset of E we denote by pos(A) the positive hull of A, that is to say the smallest convex cone containing A. A Yao-Yao partition P of an n-dimensional space E has 2 n elements and for each A in P there exists a basis v 1 , . . . , v n of E such that

A = x + pos(v 1 , . . . , v n ), ( 2
)
where x is the center of P. Indeed, assume that P is defined by F, v, P 1 and P -1 (see Definition 1). Let A ∈ P 1 and assume inductively that there is a basis v 1 , . . . , v n-1 of F such that A = x + pos(v 1 , . . . , v n-1 ). Then

A + R + v = x + pos(v, v 1 , . . . , v n-1 ).
Proposition 4. Let P be a Yao-Yao partition of E and x be its center. Any affine half-space containing x contains an element of P.

Proof. When E has dimension 0, the result is obvious. Let E have dimension n ≥ 1 and assume that the proposition holds for any affine space of dimension n -1. Let be an affine form on E such that (x) ≥ 0, and let H = {y ∈ E | (y) ≥ 0}. We use the notations of Definition 1. By the induction assumption, there exists A + ∈ P 1 and A -∈ P -1 such that ∀y ∈ A + , (y) ≥ 0 and ∀y ∈ A -, (y) ≥ 0.

Let be the linear form on E associated to , if (v) ≥ 0 then (x + tv) ≥ 0 for all x ∈ A + and t ∈ R + , thus

A + + R + v ⊂ H. Similarly if (v) ≤ 0 then A -+ R -v ⊂ H.
In both cases H contains an element of P.

Proposition 5. Let µ ∈ M(E) and let K be a convex subset of E satisfying µ(E\K) < 2 -n µ(E). Then any Yao-Yao center x of µ is contained in K.

Proof. Assume on the contrary that there is a center x of µ outside K, and let P be an equipartition with center x. By Hahn-Banach there is a halfspace containing x and disjoint from K. By Proposition 4, this half-space contains an element A of P. So on the one hand µ(A) = 2 -n µ(E) and on the other hand A ⊂ E\K, thus we get a contradiction.

Center with respect to a basis

Let E be an affine space of dimension n and L = ( 1 , . . . , n ) be a family of affine forms on E such that the map

x ∈ E → 1 (x), . . . , n (x) ∈ R n
is one to one. We say that L is a system of coordinates. If x ∈ E and v ∈ E we write x i and v i for i (x) and i (v), respectively. We also let e 1 , . . . , e n be the basis of E satisfying e i j = δ ij . Let us introduce a restricted notion of Yao-Yao partition. Definition 6. A Yao-Yao partition P of E given by F, v, x, P 1 and P -1 is adapted to L if F = {y ∈ E | y 1 = x 1 } and if P 1 and P -1 are adapted to ( 2|F , . . . , n|F ) (which is a system of coordinates on F ).

In the sequel, the vector v is called the axis of P, and it is said to be normalized when v 1 = 1. Since P has 2 n elements, there is a one to one map between P and the discrete cube {-1, 1} n . Let us construct such a map.

Definition 7. Let P be a Yao-Yao partition of E adapted to L, defined by x, P 1 , P -1 and v (with v 1 = 1). Let P(∅) = E and let (ε 1 , . . . , ε k ) be a sequence of ±1 of size k ∈ {1, . . . , n}. Recall that P ε 1 is a partition of F , hence the notation P ε 1 (•) is relative to F , for instance P ε 1 (∅) = F . We assume inductively that we have defined P ε 1 (ε 2 , . . . , ε k ) and we let

P(ε 1 , . . . , ε k ) = P ε 1 (ε 2 , . . . , ε k ) + R + (ε 1 v).
An easy induction shows that

P = P(ε) | ε ∈ {-1, 1} n . (3)
We now give some basic properties of the sets P(•). It is easy to prove by induction that if P is a Yao-Yao partition of E adapted to the basis L, then for all ±1 sequence ε of size k, there exists a sequence of vectors v 1 , . . . , v k satisfying v i j = 0 and v i i = 1 for all j < i ≤ k (we call sub-diagonal such a sequence hereafter) such that

P(ε) = x + pos(ε 1 v 1 , . . . , ε k v k ) + span(e k+1 , . . . , e n ), (4)
where x is the center of P. Besides, the first vector v 1 is equal to the axis of P (in particular, it does not depend on ε). Let x, y ∈ E and v 1 , . . . , v n and w 1 , . . . , w n be two sub-diagonal sequences of E. Observe that if

x + pos(v 1 , . . . , v n ) = y + pos(w 1 , . . . , w n ) then x = y and v i = w i for all i = 1, . . . , n. Let P and Q be Yao-Yao partitions of E satisfying P(ε) = Q(ε) for some ±1 sequence ε of length k ≤ n. Let P be the projection of E with range span(e 1 , . . . , e k ) and kernel span(e k+1 , . . . , e n ). Using (4), the equality P P(ε) = P Q(ε) and the observation above, we get (x 1 , . . . , x k ) = (y 1 , . . . , y k ) and

(v 1 , . . . , v k ) = (w 1 , . . . , w k ), (5)
where x and y are the centers of P and Q, respectively; and v and w are their normalized (v 1 = w 1 = 1) axes. If S is a set of affine functions on E (possibly empty), we let B E (S) be the smallest σ-algebra of subsets of E making all elements of S measurable. When a set A belongs to B E (S) for some S, we say that A depends only on S. Another simple consequence of ( 4) is that for any ε of length k the set P(ε) depends only on ( 1 , . . . , k ). Let α ∈ R and F = { 1 = α}. For any u ∈ E\ F we let π F,u : F + R + u → F be the projection which, for any x ∈ F and t ∈ R + , maps x + tu to x. Let A ⊂ F , we let

A u = (π F,u ) -1 (A) = {x ∈ E | x 1 ≥ α and x -(x 1 -α) u u 1 ∈ A}. (6)
Let P be a partition of E. Let k < n and ε 1 , . . . , ε k be a ±1 sequence. There exists an affine form ∈ span(1, 1 , . . . , k ) such that

P(ε 1 , . . . , ε k , 1) = P(ε 1 , . . . , ε k ) ∩ { k+1 ≥ } (7a) P(ε 1 , . . . , ε k , -1) = P(ε 1 , . . . , ε k ) ∩ { k+1 ≤ }. (7b)
Let us prove (7) by induction on k. When k = 0 we have P(1)

= F + R + v = { 1 ≥ x 1 } and P(-1) = { 1 ≤ x 1 }. If k ≥ 1, let v be the normalized (v 1 = 1)
axis of P and assume inductively that

P ε 1 (ε 2 , . . . , ε k 1) = P ε 1 (ε 2 , . . . , ε k ) ∩ { k+1 ≥ }
for some ∈ span(1, 2 , . . . , k ). Let us assume that ε 1 = 1, then

P(1, ε 2 , . . . , ε k , 1) = (π F,v ) -1 P 1 (ε 2 , . . . , ε k ) ∩ { k+1 ≥ } = P(1, ε 2 , . . . , ε k ) ∩ A v ,
where A = F ∩ { k+1 ≥ }. Using (6), we obtain

A v = {x 1 ≥ α} ∩ { k+1 ≥ } where = + ( 1 -α)(v k+1 -(v)
). Since additionnally {x 1 ≥ α} contains P(1, ε 2 , . . . , ε k ), we get (7a) when ε 1 = 1. The proof for ε 1 = -1 and the proof of (7b) are similar. From (7), we get in particular P(ε 1 , . . . , ε k ) = P(ε 1 , . . . , ε k , 1) ∪ P(ε 1 , . . . , ε k , -1).

Applying this equality again and again, we obtain, for any ε of size k

P(ε) = ∪ P(ε, τ ) | τ ∈ {-1, 1} n-k
where (ε, τ ) is the sequence obtained by concatenation of ε and τ . Therefore, if P is an equipartition for µ then

µ P(ε 1 , . . . , ε k ) = 2 -k µ(E). (8)
4 Uniqueness

In this section we prove that, under reasonable assumptions, the Yao-Yao center of a measure, with respect to a given basis, is unique. In the sequel, the space E is equipped with a system of coordinates 1 , . . . , n and all Yao-Yao partitions are adapted to this system.

Lemma 8. Let α ∈ R and F = {z ∈ E | z 1 = α}. (i) Let A ⊂ F depend only on 2 , . . . , k , let v, w ∈ E satisfying v 1 = w 1 = 1 and (v 2 , . . . , v k ) = (w 2 , . . . , w k ), then A v = A w .
(ii) Let be an affine form on E, non-constant on F and let A = F ∩ { ≥ 0}. Let v, w satisfy v 1 = w 1 = 1 and (v) > (w), then A v A w .

(iii) Again let A = F ∩ { ≥ 0}, and let (v p ) be a sequence satisfying v p 1 = 1 for all p and (v p ) ↑ +∞. Then for any µ ∈ M(E) we have

lim p→+∞ µ(A v p ) = 0.
Proof. Point (i) follows easily from (6). For (ii), observe that when A = F ∩ { ≥ 0}, equation (6) becomes

A v = {x ∈ E | x 1 ≥ α and (x) ≥ (x 1 -α) (v)}, (9)
and similarly for A w . The inclusion A v ⊂ A w follows immediately. Besides, since is non-constant on F , we can find x satisfying x 1 = 1 + α and (w) < (x) < (v), so inclusion is strict. By (ii) the sequence (A v p ) is decreasing, and clearly by (9) the intersection ∩A v p is included in F , hence (iii).

The key step is the following Lemma 9. Let P and Q satisfy

∀ε 1 , . . . , ε k , P(ε 1 , . . . , ε k ) = Q(ε 1 , . . . , ε k ) (10) ∃ε 1 , . . . , ε k+1 , P(ε 1 , . . . , ε k+1 ) = Q(ε 1 , . . . , ε k+1 ) (11)
for some k < n. If x k+1 ≥ y k+1 , then there exists δ 1 , . . . , δ k such that

P(δ 1 , . . . , δ k , 1) is strictly included in Q(δ 1 , . . . , δ k , 1).
Proof. The proof is by induction on k. If k = 0 then (11) becomes

P(1) = { 1 ≥ x 1 } = { 1 ≥ y 1 } = Q(1)
.

So x 1 = y 1 , thus x 1 > y 1 , hence P(1) Q(1).
Assume that k ≥ 1. Let v and w be the normalized (v 1 = w 1 = 1) axes of P and Q, respectively. Recall that (10) implies ( 5): (x 1 , . . . , x k ) = (y 1 , . . . , y k ) and (v 1 , . . . , v k ) = (w 1 , . . . , w k ). Also, intersecting (10) with

{ 1 = x 1 }, we get ∀ε 1 , . . . ε k , P ε 1 (ε 2 , . . . , ε k ) = Q ε 1 (ε 2 , . . . , ε k ). ( 12
)
There are three cases. If v k+1 = w k+1 , since for all ε 1 , . . . , ε k+1 , the set A = P ε 1 (ε 2 , . . . , ε k+1 ) depends only on 2 , . . . , k+1 , Lemma 8 (i) implies that

P(ε 1 , . . . , ε k+1 ) = P ε 1 (ε 2 , . . . , ε k+1 ) + R + (ε 1 v) = P ε 1 (ε 2 , . . . , ε k+1 ) + R + (ε 1 w). ( 13
) Also Q(ε 1 , . . . , ε k+1 ) = Q ε 1 (ε 2 , . . . , ε k+1 ) + R + (ε 1 w).
Therefore there must exist (δ 1 , ε 2 , . . . , ε k+1 ) such that

P δ 1 (ε 2 , . . . , ε k+1 ) = Q δ 1 (ε 2 , . . . , ε k+1 ), otherwise (11)
would fail. Recalling (12), we remark that we can apply the induction assumption to P δ 1 and Q δ 1 : there exists δ 2 , . . . , δ k such that P δ 1 (δ 2 , . . . , δ k 1) is stricly included in Q δ 1 (δ 2 , . . . , δ k 1). Then (13) shows that P(δ 1 , . . . , δ k 1) is strictly included in Q(δ 1 , . . . , δ k 1), which concludes the proof in this case. If v k+1 > w k+1 , then by (7) and Lemma 8 (ii) we obtain

P 1 (ε 2 , . . . , ε k 1) + R + v P 1 (ε 2 , . . . , ε k 1) + R + w, (14)
for all ε 2 , . . . , ε k . Therefore, it is enough to prove that there exists δ 2 , . . . , δ k such that

P 1 (δ 2 , . . . , δ k 1) ⊂ Q 1 (δ 2 , . . . , δ k 1).
This holds by the induction assumption applied to P 1 and Q 1 : either we have equality for all δ 2 , . . . , δ k , or there exists δ 2 , . . . , δ k such that the inclusion is strict.

The case v k+1 < w k+1 is similar, we just need to deal with P -1 and Q -1 instead of P 1 and Q 1 .

We now let M * (E) be the set of measures µ ∈ M(E) which satisfy µ(U ) > 0 for every open set U . Here is the main result of this section.

Proposition 10. Let µ ∈ M * (E), then µ has at most one center with respect to the system L. Moreover the k first coordinates of the center of µ depend only on the restriction of µ to B E ( 1 , . . . , k ).

Proof. Let µ, ν ∈ M * (E) satisfy ν(A) = µ(A) for all A ∈ B E ( 1 , . . . , l ). We can assume that µ(E) = 1(= ν(E)). Let P and Q be equipartitions for µ and ν, respectively, and let us call x and y the respective centers of P and Q. We want to prove that (x 1 , . . . , x l ) = (y 1 , . . . , y l ). If P = Q, there exists k satisfying (10) and (11), we have in particular (x 1 , . . . , x k ) = (y 1 , . . . , y k ) (

and we may assume that x k+1 ≥ y k+1 . By the previous lemma, there exists (ε 1 , . . . , ε k+1 ) such that P(ε 1 , . . . , ε k+1 ) is strictly included in Q(ε 1 , . . . , ε k+1 ). Since these sets are polytopes, their difference Q(ε 1 , . . . , ε k+1 )\P(ε 1 , . . . , ε k+1 ) has non-empty interior. Recall that the set Q(ε 1 , . . . , ε k+1 ) depends only on 1 , . . . , k+1 . So if k + 1 ≤ l, applying (8), and the fact that µ(U ) > 0 for any open set U , we get

2 -k-1 = µ P(ε 1 , . . . , ε k+1 ) < µ Q(ε 1 , . . . , ε k+1 ) = ν Q(ε 1 , . . . , ε k+1 ) = 2 -k-1 ,
which is absurd. Therefore l ≤ k, and the result follows from (15). The uniqueness of the center is obtained by letting ν = µ and l = n.

Continuity

We now deal with the continuity of the center, for this we need a topology on M(E). Let (µ p ) p∈N be a sequence of elements of M(E) and µ ∈ M(E).

The sequence (µ p ) converges narrowly to µ when

φ dµ p → φ dµ (16)
for any function φ continuous and bounded on E. A subset F of M(E) is tight if for every ε > 0 there exists a compact subset K of E such that µ(E\K) < ε for all µ ∈ F. A converging sequence of measures is obviously tight.

Recall that when T : Ω 1 → Ω 2 is a measurable map, and µ is a measure on Ω 1 , the image measure T # µ of µ by T is defined by

T # µ(A) = µ T -1 (A) ,
for every measurable subset A of Ω 2 .

Let F be an affine hyperplane of E. For any µ ∈ M(E) and v ∈ E\ F we let µ

F,v = (π F,v ) # µ, hence ∀A ⊂ F, µ F,v (A) = µ(A + R + v). (17)
Let (ν p ) be a sequence of elements of M(E) converging narrowly towards ν ∈ M(E). Notice that if A ⊂ E is a polytope, then ν(∂A) = 0 and thus ν p (A) → ν(A). Let (v p ) be a sequence of elements of E\ F converging to v ∈ E\ F . Then

ν p F,v p → ν F,v , narrowly.
Indeed, by tightness of the set {ν, ν p | p ∈ N}, we can assume that all these measures are supported by a compact set K. We can also assume that F + R + v p = F + R + v for all p. Then, setting µ p = ν p F,v p and µ = ν F,v , it is easy to see that the supports of µ and µ p for all p are contained in a compact subset L of F . Let φ be continuous on F , then φ is uniformly continuous on L, from which we get φ • π F,v p → φ • π F,v , uniformly on K. The result follows easily. In the same spirit, if (u p ) is a sequence of elements of E converging to 0 and if T p is the translation of vector u p , then (T p ) # ν p → ν, narrowly. The following lemma is an easy consequence of Proposition 5.

Lemma 11. Let F ⊂ M(E) be tight and such that {µ(E) | µ ∈ F} is bounded. Then there exists a compact subset of E containing any center of any element of F.

We now reformulate the definition of a Yao-Yao center. Let µ ∈ M(E). Then an element x ∈ E is a center for µ according to ( 1 , . . . , n ) if and only if letting F = { 1 = x 1 } there exists v ∈ E\ F such that -µ(F + R + v) = 1 2 µ(E).

and similarly for x (-1) . Therefore T satisfies (a).

Let be an affine form on E, and (v p ) be a sequence of elements of F such that (v p ) → +∞. Then it follows from Lemma 8 (iii) that for any m > 0 and any big enough p we get µ u+v p (F ∩ { ≥ -m}) < 2 -n µ(F ), which by Proposition 5 implies that (x (1) (v p )) ≤ -m. Hence

x (1) (v p ) → -∞.

Similarly x (-1) (v p ) → +∞. Thus T satisfies (b). Therefore T is onto.

There exists v ∈ F such that T v = 0, then x (1) (v) = x (-1) (v) which concludes the proof.

Lemma 15. If µ ∈ M * (E) has a center of symmetry z then z is the unique Yao-Yao center of µ, whatever the basis L.

Proof. Let x be the Yao-Yao center of µ and s : E → E be the symmetry of center z. Then clearly s(x) is a center for s # µ, with respect to L. Since s # µ = µ and by uniqueness of the center, we obtain s(x) = x. Therefore x = z.

We are now in a position to prove Theorem 3. Let E be an affine space of dimension n and ( 1 , . . . , n ) be a system of coordinates. Let γ be an arbitrary element of M * (E). Let µ ∈ M(E) and for p ≥ 1 let µ p = µ + 1 p γ. Then obviously µ p ∈ M * (E) and µ p → µ, narrowly. Let x p be the center of µ p with respect to ( 1 , . . . , n ). By Lemma 11 and Corollary 13, the sequence (x p ) is bounded and the limit of any of its converging subsequence is a center for µ, so µ has a center. If µ is symmetric with respect to x, then we let γ be an element of M * (E) symmetric with respect to x. Then so is µ p , by the preceding lemma we get x p = x for any p, therefore x is a center for µ.

-the point x is a center for both µ F,v and µ F,-v , according to the system of coordinates ( 2|F , . . . , n|F ).

In the sequel, such a vector v, is called an axis for (µ, F ), and we say that

Here is a first continuity property of the center.

Lemma 12. Let (µ p ) be a sequence of elements of M(E) converging narrowly towards µ ∈ M(E). If all the measures µ p share a common center x with respect to ( 1 , . . . , n ); then x is also a center of µ.

Proof. The proof is by induction on the dimension n of E. When n = 1, the result is obvious: if x is a median for all the measures µ p then it is a median for µ. We assume that n ≥ 2 and that the result holds for any affine space of dimension n -1. We can also assume that µ and the measures µ p are probability measures. Let F = { 1 = x 1 }. For all p there exists a normalized axis v p for (µ p , F ). We claim that the sequence (v p ) is bounded. Indeed otherwise there exists affine on E such that (v p ) → +∞. Let H = { ≥ 0} ∩ F . Let > 0, by Lemma 8 (iii), there exists w such that w 1 = 1 and µ(H w ) < . For p big enough (v p ) > (w), applying Lemma 8 (ii) we get H v p ⊂ H w . Also for p large, we have µ p (H w ) ≤ µ(H w ) + . Thus

Taking small enough, it follows from Proposition 5 that for any such p the center x of µ p F,v p does not belong to H. Hence (x) ≤ 0. Let m > 0, the same holds if we replace by + m, so we get (x) ≤ -m for all m > 0, which is absurd.

Up to an extraction we can assume that (v p ) has a limit, say v (which satisfies v 1 = 1). Then µ p F,v p → µ F,v narrowly. By the induction assumption, x is a center for µ F,v with respect to the basis ( 2|F , . . . , n|F ), and the same holds for µ F,-v . Therefore x is a center for µ.

Corollary 13. Let (µ p ) be a sequence of elements of M(E) converging narrowly towards µ ∈ M(E). If every measure µ p has a center x p with respect to ( 1 , . . . , n ) and if the sequence (x p ) has a limit x, then x is a center of µ with respect to ( 1 , . . . , n ).

Proof. Let v p = x -x p and T p be the translation of vector v p . Since v p → 0 we have (T p ) # µ p → µ, and clearly x = T p (x p ) is a center for (T p ) # µ p . Then the result follows from the previous lemma.

Proof of Theorem 3

Let us start with an easy fact. Let V be a vector space of finite dimension n, with a given basis (e 1 , . . . , e n ). Let T : V → V be continuous and satisfy the following properties:

Then T is onto. Indeed, let u in V and u 1 , . . . , u n be its coordinates. By (b) and the continuity of T , we can find v 1 ∈ R such that the first coordinate of T (v 1 e 1 ) is u 1 . Then we find v 2 ∈ R such that the second coordinate of T (v 1 e 1 + v 2 e 2 ) is u 2 , and by (a) the first coordinate of T (v 1 e 1 + v 2 e 2 ) is still u 1 . And so on. Proposition 14. Let E be an affine space and ( 1 , . . . , n ) be a system of coordinates. Any element of M * (E) admits a unique center with respect to ( 1 , . . . , n ).

Proof. We have already proved the uniqueness of the center, we shall prove its existence by induction on the dimension. If E has dimension 1 then we just have to show that any µ ∈ M * (E) has a median, which is clear. If n ≥ 2, we assume that the proposition holds for any affine space of dimension n-1. Let α ∈ R satisfy µ{ 1 ≥ α} = 1 2 µ(E) and let F = { 1 = α}. Let u ∈ E satisfy u 1 = 1. By the induction assumption, for all v ∈ F and ε ∈ {-1, 1}, the measure µ F,ε(u+v) admits a unique center with respect to ( 2 , . . . , n ) which we call x (ε) (v). If we can prove that there exists v such that x (1) (v) = x (-1) (v), then we are done. We define

If a sequence (v p ) goes to v, then µ F,u+v p converges narrowly to µ F,u+v . Then by Lemma 11, the sequence x (1) (v p ) p is bounded, and by Corollary 13, any of its converging subsequences goes to the unique center of µ F,u+v . Therefore x (1) (v p ) → x (1) (v), and similarly for x (-1) , hence the continuity of T . If (v 2 , . . . , v k ) = (w 2 , . . . , w k ) then, by Lemma 8 (ii), we have A u+v = A u+w for all A ∈ B F ( 2 , . . . , k ), hence µ F,u+v (A) = µ F,u+w (A). By Proposition 10, this implies that for i = 2, . . . , k i x (1) (v) = i x (1) (w) ,


J-P Merlet
email: jean-pierre.merlet@sophia.inria.fr

A Local Motion Planner for Closed-loop Robots

Global motion planners have been proposed for closed-loop robot based on the same paradigm than has been proposed for serial chains. First a sparse representation of the configuration space of the robot is constructed as a set of nodes. This is somewhat more complicated than for serial chain as the closure equations of the mechanism should be satisfied. Then a motion planning query consists simply in connecting the start and goal points through an appropriate set of nodes (usually minimizing the length of the trajectory). But such motion planner should be complemented by a local motion planner that addresses the following issues:

1) ensure that two successive nodes belong to the same robot kinematic branch (otherwise connecting these nodes will require to disassemble the robot) 2) verify that all poses between nodes satisfy the robot constraints (if possible taking into account the uncertainties in the robot modeling) 3) eventually try to shorten the trajectory length We present such a local motion planner that addresses all three issues and illustrates its use on a Gough parallel robot.

I. INTRODUCTION

A. Global motion planner

Motion planning is a classical problem in robotics and has been largely addressed for serial chains [START_REF] Canny | The complexity of robot motion planning[END_REF], [START_REF] Latombe | Robot Motion planning[END_REF], [START_REF] Lavalle | Planning algorithm[END_REF]. Among the most successful method we may mention the roadmap approach: a representative, but limited, set of reachable poses, called the sampling tree, is pre-computed and a planning query consists in connecting the reachable poses so that, for example, the length of the trajectory is minimal. An interest of this approach is that the construction of the sampling tree may be done once off-line, while finding a trajectory relies on determining a shortest path within a graph, a task for which there are efficient and fast algorithms. For serial chains the sampling is performed usually in the joint space (a point in this space leads to a unique configuration of the robot). Unfortunately adapting this approach to closed-loop chains is not easy as the joint variables should satisfy the closure equations and cannot be arbitrary chosen.

B. Closed-loop robots and kinematic constraints

The planning problem for serial chains is usually related to avoiding obstacles and self-collision, while for closedchains kinematic constraints become preponderant.

A specific kinematic constraint for closed-loop structure is that their operational configuration space may have different components that are not connected. This may be illustrated simply on a four-bar mechanism (figure 1). Fig. 1. A 4-bar mechanism with an operational space having two not connected components. Going from a pose in a component to another in the second component requires to disassemble the mechanism.

A consequence is that even if a sampling method is able to determine poses that satisfy the closure equations, they may still belong to different kinematic branches and therefore cannot be connected by a trajectory.

Other kinematic constraints are joint limits, collision avoidance between the robot's bodies and the absence of singularity on the trajectory. All these constraints are functions of the pose of the robot and of parameters that describe the geometry of the robot and must be satisfied although there are uncertainties on these parameters.

The complexity of the kinematic constraints and the uncertainties leads to complicated configuration space. In some cases this configuration space can be constructed either analytically [START_REF] Trinkle | Complete path planning for closed kinematic chains with spherical joints[END_REF] or numerically [START_REF] Merlet | Determination of 6D workspaces of Gough-type parallel manipulator and comparison between different geometries[END_REF] but manipulating these representations for motion planning is a complicated task. One of the most successful motion planning method for closed-chains is called the Probabilistic Roadmap [START_REF] Cortés | Probabilistic motion planning for parallel mechanisms[END_REF], [START_REF] Xie | A kinematics-based probabilistic roadmap method for high dof closed chain systems[END_REF], [START_REF] Yakey | Randomized path planning for linkages with closed kinematic chains[END_REF]. Here the sampling poses are chosen at random although general geometric algorithms are used to guide the sampling so that the poses satisfy the closure equations.

C. Local planner

In all cases these global motion planners must be complemented by a local motion planner that addresses the following issues:

1) ensure that two successive nodes in the trajectory belong to the same kinematic branch of the robot. 2) find a path that connect two successive nodes of the trajectory, such that any pose on the path are guaranteed to satisfy the kinematic constraints 3) eventually improve the trajectory in terms of some optimality criteria (e.g. its length) and provide a solution that is "close" to the global optimum of the criteria (we will precise this closeness issue later on)

II. BASIC INGREDIENTS

The local planner we have developed relies on the following assumptions:

• being given a pose of the end-effector we are able to state if the kinematic constraints are satisfied • all the kinematic constraint may be written as C(X, Θ) ≤ 0, where X is a n-dimensional vector constituted of the parameters that describes the endeffector pose and Θ the joint variables (active and/or passive) • the workspace of the robot is bounded On the other hand we don't impose any particular form for C or any choice on the parameters X.

The output of the planner will be a list of way points in the operational space, the start point S being the first element of the list while the goal point G will be the last element (S, G will be provided by the global planner). A trajectory between two way points W j , W j+1 will be defined by a set F of n continuous, analytical functions F k such that any pose X t on the trajectory may be written as X t = F(W j , W j+1 , t), where t is the time and lie in the range [0,1] and such that F(W j , W j+1 , 0) = W j , F(W j , W j+1 , 1) = W j+1 . At time t the k-th component of X t will be obtained as F k (W j k , W j+1 k , t). Any choice of F k is possible as soon as it is possible to formulate analytically the optimality criteria for the trajectory (if any is used). For the sake of simplicity we will use here linear functions i.e. X t = W j + t(W j+1 -W j ).

The number of way points may not be defined in advance. If an optimality criteria is used, the final number of way points of the trajectory will be such that adding a new way point will not improve significantly the value of the optimality criteria. Note however that this local planner is only able to provide a limited number of way points (typically 3). But requiring more way points means that the global planner has performed poorly.

Any pose belonging to the trajectory between two way points will be guaranteed to satisfy the kinematic constraints. This is evidently a key issue that is addressed in the next section.

A. Checking a trajectory between two way points

We have already proposed a method to check if a given trajectory, defined by analytical function of the time, satisfies the kinematic constraints [START_REF] Merlet | A parser for the interval evaluation of analytical functions and its applications to engineering problems[END_REF] and we will just outline its principle. We rely on interval analysis to determine if a kinematic constraint C i is satisfied or violated over a given time range. Being given an analytical function f (x 1 , . . . , x m ) and ranges [x i , x i ] for the unknowns x i , interval analysis allows to determine simply an interval [a, b], called an interval evaluation of f , such that for all x i in [x i , x i ] we have a ≤ f (x 1 , . . . , x m ) ≤ b i.e. a, b are lower and upper bound of the minimum and maximum of f over the ranges. Usually a, b overestimate the minimum and maximum but the overestimation decreases with the width x i -x i of the range. A major interest of interval analysis is its robustness with respect to round-off errors: even if such errors occurred, the range [a, b] is still guaranteed to include the real minimum and maximum of f .

Being given a time interval [t 0 , t 1 ] this method will be used to determine ranges for the pose parameters along the trajectories defined by the set F. In turn these ranges allows one to determine a range [a i , b i ] for the value of each kinematic constraint C i . If a i > 0, then the constraint is violated at any time in [t 0 , t 1 ], while if b i < 0, then the constraint is satisfied all over the time range. It may occur that due to the overestimation of interval analysis we get a i < 0, b i > 0, so that we cannot state if the constraint is violated or satisfied. In that case we will bisect the time range in two ranges [t 0 , (t 0 + t 1 )/2], [(t 0 + t 1 )/2, t 1 ] that will be stored in a list. All ranges in this list will be processed until we may determine for all boxes in the list that either all constraints are satisfied or that for some time range a constraint is violated. It may also occur that for a given time interval analysis does not allow to determine if one (or more) kinematics constraint is satisfied or violated because round-off errors do not allow to determine the sign of the constraints.

It must be noted that although we have assumed an analytical form for C, this is not strictly necessary. Indeed our approach requires only to be able to calculate an interval evaluation of an index stating if the constraint is satisfied or violated. For example for closed-chain robot a singularity occurs if the Jacobian matrix is singular. To detect such occurrence we use the sign of the determinant without having its analytical form: only analytical forms of the components of the Jacobian are required [START_REF] Merlet | On the regularity of the inverse jacobian of parallel robot[END_REF].

A first extension of this algorithm is to consider uncertainties in the geometric modeling of the robot. In mechanical engineering they correspond to the bounded manufacturing tolerances. We use a worst case scenario by assuming that the real value of the geometrical parameters may be any value within known ranges. If there width are relatively small their ranges are used as it for the interval evaluation of the kinematic constraints. In the worst case it may happen that for a given fixed time that the interval evaluation of one (or more) kinematic constraint has a lower negative bound and positive upper bound, thereby not allowing to determine if the constraint is satisfied: in that case the parameters uncertainties are added as new variables and submitted to the bisection process. Note that the same process may be used if we assume control errors, i.e. the robot will not follow exactly the planned trajectory, the differences between the trajectory and the robot pose being still bounded.

At this point we are thus able to design an algorithm A P (W j , W j+1 ) that returns 1 if the trajectory between these two way points satisfies the kinematic constraints, -1 if it violate them or 0 if the trajectory is unsafe (i.e. at some time we cannot determine if a kinematic constraints is satisfied or violated).

We may extend this algorithm to deal with set of way points defined by a box in the operational workspace. This new algorithm, A P (W j , W j+1 ), takes as input two boxes of the operational workspace that include the potential way points W j , W j+1 and returns 1 if the trajectory between any poses in the boxes W j , W j+1 satisfies the constraint, -1 if all trajectories violate them. Furthermore the algorithm returns 0 if it is not able to complete its task after a limited number of iterations.

B. The optimality criteria

We will assume that we can calculate an interval evaluation of the eventual optimality criteria H (without loss of generality we will assume that H should be positive and minimal). Our algorithm will determine a trajectory with a criteria that is at most , that is defined by the end-user, away from its minimal value.

Now assume that we have determined a trajectory for which the value of H is l: using the interval evaluation of the criteria we are able to design an algorithm O(W 1,...,m ) that takes as input a set of boxes for the m way points and returns 1 if the lower bound of the interval evaluation of H is lower than l -(i.e. there may be trajectories having a value of H that is lower than lwith way points in the boxes ), -1 if the lower bound of the interval evaluation of H is greater than l-(i.e. even if there are valid trajectories with way points in the boxes, they will have a higher value of H than the current trajectory). Otherwise the algorithm will return 0.

III. THE LOCAL PLANNER

The first objective of the local planner is to determine if there is valid trajectory with one way point. This eventual way point must be be located in the workspace, hence we have bounds for its coordinates. A way point box is a box of the operational space that may contain the way point W 1 . The robot workspace will be decomposed into way point boxes and for each box we will use algorithm A P on the trajectory components SW 1 , W 1 G to test if the trajectory with the center of the box as W 1 is valid. If not we will use algorithm A P to determine if the current box does not contain any valid trajectory, in which case we will discard it, otherwise the box will be bisected. We will not elaborate on the bisection strategy (i.e. which variable will be bisected) but it plays an important role in the efficiency of the algorithm.

If an optimality criteria has been defined and a valid trajectory has been found with H = l, then we use O to determine if the box may include trajectories with H < l -. If this is not the case the box will be discarded.

A. Incremental addition of way points

It may occur that there is no valid trajectory with only one way point or that the purpose of the planner is to determine a trajectory whose value of H is optimal up to the accuracy . In that case we will incrementally add way points until either a trajectory (or a better value of H) is found or we are able to determine that the current trajectory is optimal. The dimension of the way point boxes will be k n where k is the number of way points.

B. The algorithm

We present here the algorithm for the case where an optimality criteria has been defined. The number of way points will be denoted N and the algorithm starts with N = 1. A flag T , initially set to 0, will be set to 1 if a valid trajectory has been found. The value of H for this trajectory will be denoted H T . The algorithm will process a list L of way point boxes B = {W 1 , . . . W N }. The i -th element of the list will be B i and r will denote the total number of way point boxes in the list. The mid point of a way point box is obtained by taking the mid point of all the ranges of the variables in the box: it correspond to a set of fixed poses and will be written as M id(B i ). The value of the optimality criteria obtained for a trajectory with N way point will be denoted l N o (l 0 o will be equal to +∞). When starting the algorithm we have T = 0 and there is one way point box B 1 in L, whose ranges correspond to the ranges of the robot workspace E. The Cartesian product of N box E will be denoted E N .

r = i = 1, L = {E N }, T = 0 while true do if T = 1 then l N o = l N -1 o while i ≤ r do if T = 1 and O(B i ) = -1 then i = i + 1, next if A P (B i ) = -1 then i = i + 1, next if A P (M id(B i )) = 1 then l w = H(M id(B i )) if T = 1 and l w < l N o then l N o = l w if T = 0 then T = 1, l N o = l w bisect B i , r = r + 2, i = i + 1, next end-do if T = 1 and l N o < l N -1 o -then N = N + 1, next if T = 1 and l N o = l N -1 o then break if T = 0 then N = N + 1, next end-do
Efficiency of the algorithm will be improved by adding filtering methods, i.e. heuristics that reduce the size of a way point box without using bisection. There are standard interval analysis filtering methods [START_REF] Hansen | Global optimization using interval analysis[END_REF] but also adapted methods as will be shown in the examples.

Note that exiting the algorithm as soon as adding only one way point does not allow to improve significantly l o (i.e. leads to a decrease larger than ) does not allow to guarantee that a close approximation of the optimal trajectory has been found. In practice however this will imply that the global planner has performed poorly.

A first drawback of this algorithm is that its complexity quickly increase with the number of way points. This will be confirmed by the example: determining a trajectory with one or two way points is fast but the computation time drastically increase for 3 way points. But we may reduce the computation time for 2 (or more) way points with an incremental use of the local planner, an approach that we will call the intermediary step method. Assume that a trajectory with one way point has been found and let S, W 1 , G be the current trajectory. We use the local planner with a threshold equal to, for example, /2 to determine if there is a trajectory with one way point W 1 between S, W 1 that reduces the value of H compared to its value in the current trajectory. If this is the case we get the trajectory with 2 way points S, W 1 , W 1 , G that has already a better H than the one way point trajectory. The process may be repeated with W 1 , G and we will select the two way points trajectory that leads to the best improvement of H. We get therefore a better initial trajectory for the local planer when looking for a trajectory with 2 way points.

A second drawback of this algorithm is its lack a memory. Indeed when checking a trajectory with two or more way points we will examine possible location for W 1 , W 2 that have already been determined as non valid when we have looked for a trajectory with one way point.

IV. EXAMPLES

Four our tests we have used a 6 d.o.f. parallel robot that is presented in the next section.

A. Parallel robot

As example of complex 6-dof closed-loop robot we will consider a Gough platform (figure 2). The fixed frame (O, x, y, z) will be called the base frame while a mobile frame (C, x r , y r , z r ) attached to the platform will be called the mobile frame. The pose of the platform will be parametrized by the location of C in the base frame and 3 angles will be used to define the orientation of the mobile frame with respect to the base frame. Motion of the platform are obtained by changing the length of the 6 legs, that are attached on the base at A i and on the platform at B i . The coordinates of the vector OA i in the reference frame and the coordinates of CB ir in the mobile frame are known. The legs are restricted to have a minimal and maximal length, thereby restricting the workspace of this robot. Being given the location of the center C of the platform in the reference and the rotation matrix R between the mobile and reference frame, the length ρ of a leg is determined as

ρ 2 = ||AO + OC + RCB r || 2 (1)
and we should have

ρ 2 min ≤ ρ 2 ≤ ρ 2 max
The shape of the workspace due to this limitation is quite complicated [START_REF] Merlet | Determination of 6D workspaces of Gough-type parallel manipulator and comparison between different geometries[END_REF] and planning a trajectory to fully lie in this workspace is a complex task. Other kinematic constraints may be considered as well: limitation of the passive joints motion at A, B [START_REF] Su | Trajectory planning for constrained parallel manipulators[END_REF], singularity avoidance [START_REF] Bhattacharya | Comparison of an exact and an approximate method of singularity avoidance in platform type parallel manipulators[END_REF], [START_REF] Dasgupta | Singularity-free path planning for the Stewart platform manipulator[END_REF], [START_REF] Dash | Workspace generation and planning singularity-free path for parallel manipulators[END_REF], [START_REF] Jui | Path tracking of parallel manipulators in the presence of force singularity[END_REF], [START_REF] Sen | Variational approach for singularity-path planning of parallel manipulators[END_REF] or leg interference [START_REF] Chablat | Moveability and collision analysis for fully-parallel manipulators[END_REF], [START_REF] Merlet | Legs interference checking of parallel robots over a given workspace or trajectory[END_REF]. We will consider here only the limitation on the leg lengths and limits on the motion of the passive joints.

The optimality criteria that we will use is the length of the trajectory of the center C of the platform. This allows to use a dedicated filtering method that is available as soon as a trajectory of length l with one way point has been determined. Indeed any potential way point that will lead to a smaller trajectory length (i.e. smaller than l -) will lie within an ellipsoid which is defined as the set of points M such that the sum of the distances from M to S, G is equal to l -. A bounding box E of this ellipsoid may be found and its intersection I with the way point box B is calculated. If I is empty, then the box is discarded otherwise I is substituted to B.

B. Implementation and tests

To test the motion planning algorithm we have used the C++ interval arithmetic package BIAS/Profil and components of our interval analysis library ALIAS.

We assume first that the geometry of the robot is perfectly known (i.e. the location of the A i , B i are exact as indicated in table I). The minimal and maximal leg lengths The various poses are defined by the 3 coordinates of C followed by the three rotation angles. The start point S is (0,0,52.1,0,0,0) and the goal G is (11,5,52.1,0,0,0) (having pose with the same orientation and in an horizontal plane allows for an easy representation of the workspace). The distance between S, G is 12.083. First we consider only the limitation on the leg lengths as kinematic constraints. It may been that on figure 3 that the workspace has voids in its interior. We first try to determine a trajectory that keep constant the orientation of the platform. Trajectories with 1, 2 and 3 way points are presented in figure 3: their lengths are respectively 19.5373, 17.1118, 16.7887 and these trajectories are established respectively in 1 second, 17 seconds and 11 mn 44 seconds on a DELL D400 laptop. As expected the computation time grows quickly with the number of way points. But the computation time does not significantly increase if we add as kinematic constraint that the determinant of the inverse Jacobian matrix has a value larger than a given threshold for any pose on the trajectory in order to ensure that the trajectory lies on a given kinematic branch. If we now relax the constraint that the trajectory should lie in the plane z = 52.1 by allowing a z value in the range [50,55] while keeping the same orientation we found out in 0.26 second that the trajectory with one way point at (5.562,2.5,52.5351,0,0,0) has a length of 12.1144 (figure 4). As the minimal value of H is 12.083, the algorithm stops as adding new way point cannot make the length decrease by more than 0.3. Setting the threshold to 0.01 allows to determine a trajectory with length 12.0917 in a computation time of 1.2 seconds. We may also keep the trajectory in the same plane while allowing to change the platform orientation. For the same trajectory and allowing the orientation angles to lie in the range [-5,5] degrees, we find a one way point trajectory of length 12.3967 in 15 seconds and the planner immediately determine that adding a new way point will not shorten the trajectory length by more than 0.3.

S G W 1 1 W 2 1 W 3 1 W 2 2
We may also have to deal with uncertainties in the robot modeling. We assume a tolerance error of ± 0.01 on each coordinates of the A i , B i meaning for example that the x coordinate of A 1 may have any value in the range [-9.01,-8.99]. A trajectory with one way point and length 21.2389 is found in 8 seconds, with two way points and length 17.8501 in 5 minutes and 29 seconds. For three way points a trajectory of length 17.4191 is found in 4h27mn49s. In the later case by using the intermediary step method we may reduce this time to 3h41mn and the intermediary path allows to get an initial three way point trajectory of length 17.6235 in 16mn10s. These trajectories are presented in figure 5. If we allow z to lie in the range [50,55] a minimal Fig. 5. Trajectories with one to three way points when the coordinates of A i , B i have a tolerance of ± 0.01 length trajectory is found in 2.65s.

Assume now that we have limits on the passive joint motion. To model this limits for the joints at A i we assume that the angle between the i -th leg and a fixed direction defined by the unit vector n i should be lower than a threshold µ i . This constraint may be written as

| A i B i .n i ρ i | ≤ cos(µ i ) (2)
We consider the trajectory between S (0,0,52.2) and G (-8,5,52.2) and assume that the angle between the vertical direction and the legs should not exceed 17 degrees. Without considering these constraints we find a trajectory of length 16.3498 with one way point and of length 13.9694 with two way points (figure 6). But if the passive constraints are taken into account the algorithm found out in 3 seconds that there is no trajectory with only one way point. The previously calculated trajectory with two way points is also not satisfactory as part of it violate the joint limits. But a two way point trajectory that takes into account the joint limits may be found with a length of 14.1181 with

V. CONCLUSION

Safe motion planning of closed-loop robot requires to verify complex kinematic constraints apart of the classical obstacle avoidance problem. Global motion planner have the advantage to provide a draft trajectory but a local planner must then be used to provide a certified trajectory (i.e. such that the kinematic constraints are satisfied all along the trajectory). The local planner proposed in this paper is based on interval analysis and allows both to certify the trajectory and to manage possible uncertainties in the robot modeling.

Fig. 2 .

 2 Fig. 2. A Gough-Stewart platform

Fig. 3 .

 3 Fig. 3. Trajectories with 1, 2 (squared way points) and 3 (circled way points) way points.

Fig. 4 .

 4 Fig. 4. A 3D trajectory. Note that the scales are different on the x, y, z axis for a better understanding of the trajectory.

Fig. 6 .

 6 Fig.[START_REF] Dash | Workspace generation and planning singularity-free path for parallel manipulators[END_REF]. Trajectories with one (1) and two (2) way points when the constraints on the passive joint limits are not considered.

Acknowledgment: this work has been partly supported by the European ARES project N • 015653.