Mohammed Kechadi

Jean Gamby
email: jean.gamby@upmc.fr

Lila Chaal

Hubert Girault

Boualem Saidani

Bernard Tribollet

Polymer microchip impedance spectroscopy through two parallel planar embedded microelectrodes: Understanding the impedance contribution of the surrounding polymer on the measurement accuracy

Keywords: Microchip, Contact free impedance, Constant-phase element, Polyethylene terephthalate, Photoablation

The present work describes a new methodology for contact free impedance of a solution in a polymer microchip taking into account the role played by the surrounding polymer on the impedance accuracy. Measurements were carried out using a photoablated polyethylene terephthalate (PET) microchannel above two embedded microband electrodes. The impedance diagrams exhibit a loop from high frequencies to medium frequencies (1 MHz to 100 Hz) and a capacitive behaviour at low frequencies (100 Hz to 1 Hz). The impedance diagrams were corrected by eliminating from the global microchip response the contribution of the impedance of the PET layer between the two microband electrodes. This operation enables a clear observation of the impedance in the microchannel solution, including the bulk solution contribution and the interfacial capacitance related to the surface roughness of the photoablated microchannel. Models for the impedance of solutions of varying conductivity showed that the capacitance of the polymer-solution interface can be modeled by a constant phase element (CPE) with an exponent of 0.5. The loop diameter was found to be proportional to the microchannel resistivity, allowing a cell constant around 4.93•10 5 (m -1 ) in contactless microelectrodes configuration.

INTRODUCTION

The development of miniaturization technology is increasing, due to the well-known advantages and the crucial need for micrometer scale measurements with improved detection sensitivity [START_REF] Rossier | Microchannel networks for electrophoretic separations[END_REF][START_REF] Yakovleva | Microfluidic Enzyme Immunoassay Using Silicon Microchip with Immobilized Antibodies and Chemiluminescence Detection[END_REF][START_REF] Manz | Miniaturized total chemical analysis systems: A novel concept for chemical sensing[END_REF]. Hence, micro-total analysis systems (µTAS) are usually designed for a well-defined analytical task [START_REF] Guber | Microfluidic labon-a-chip systems based on polymers fabrication and application[END_REF]. A wide range of detection techniques has been employed in microfluidic devices, such as optical detection, laser induced fluorescence, luminescence, electrochemistry, surface plasmon resonance, and mass-spectroscopy [START_REF] Cheek | Chemiluminescence Detection for Hybridization Assays on the Flow-Thru Chip, a Three-Dimensional Microchannel Biochip[END_REF][START_REF] Epstein | Fluorescence-based nucleic acid detection and microarrays[END_REF][START_REF] Dill | Immunoassays and sequence-specific DNA detection on a microchip using enzyme amplified electrochemical detection[END_REF]. However, electrochemical detection (ECD) offers specific advantages for microfluidic systems due to its small size, portability, low cost, high sensitivity, and high selectivity using a proper choice of detection potential and/or electrode material [START_REF] Dill | Electrochemical Detection on Microarrays[END_REF]. An additional advantage of ECD is the simplicity of the instrumentation, resulting in low electrical power requirements for in-field use [START_REF] Martin | In-Channel Electrochemical Detection for Microchip Capillary Electrophoresis Using an Electrically Isolated Potentiostat[END_REF][START_REF] Gamby | Supercapacitive Admittance Tomoscopy[END_REF][START_REF] Kjeang | Integrated electrochemical velocimetry for microfluidic devices[END_REF][START_REF] Ordeig | On-Chip Electric Field Driven Electrochemical Detection Using a Poly(dimethylsiloxane) Microchannel with Gold Microband Electrodes[END_REF]. Electrochemical impedance spectroscopy (EIS) is a sensitive analysis method for biological samples [START_REF] Gamby | Nanomosaic network for the detection of proteins without direct electrical contact[END_REF][START_REF] Dittami | A Multilayer MEMS Platform for Single-Cell Electric Impedance Spectroscopy and Electrochemical Analysis[END_REF][START_REF] Sun | On-chip epithelial barrier function assays using electrical impedance spectroscopy[END_REF][START_REF] Jagtiani | High throughput microfluidic electrical impedance flow cytometry for assay of micro particles[END_REF]. Indeed, the possibility to integrate impedance analysis into existing lab-on-a-chip designs offers great potential to obtain new insights and to achieve a profound understanding of biological processes, which are not accessible through conventional analysis. With this method, the bias of the sensing electrode potential can be controlled, even though some unwanted problems such as bubble generation, microelectrode passivation, and/or contamination can occur. To overcome these limitations, use of electrodes that do not contact electrolyte has been investigated, mainly in capillary electrophoresis with electrical detection. It is interesting to stress that non-contact conductivity measurements using external tubular electrodes all around the capillary was introduced in 1998 by Zemann et al [START_REF] Zemann | Contactless Conductivity Detection for Capillary Electrophoresis[END_REF] due to a problem encountered with electrical coupling between the high voltage applied in the capillary and the conductivity cell detection located at the end-column. However, in the literature, the signal excitation is high (5 -10 V) due to high impedance of silica capillary with pasted electrodes. In this case the linearity of the system is not maintained during measurement.

The growth of microfluidic technologies with widely varying substrate requirements has led to an increased interest in polymeric materials due to their good bio-chemical performance and low cost. Recently, Gamby et al. [START_REF] Gamby | Supercapacitive Admittance Tomoscopy[END_REF][START_REF] Gamby | Nanomosaic network for the detection of proteins without direct electrical contact[END_REF] described microchannels and contact-free microelectrodes using a microfabrication process based on photoablation of polyethylene terephthalate (PET). This process is an attractive alternative to other expensive, laborious, and time-consuming methods for microchannel fabrication such as photolithography approaches. However, when choosing a polymer-based substrate, the properties of the material are critical for both the fabrication process and the successful application of the device. Therefore, we propose in the presented work, a thorough study based on a sensitive non-contact technique in order to understand the physical phenomena taking place within the PET microsystem. Measurements were undertaken of dielectric response to a perturbation of voltage between two microelectrodes embedded in the dielectric polymer. Simulation and experimental results were used throughout the paper to illustrate the effectiveness of such an approach. This approach can be employed to detect any kind of molecule. Further development to identify the type of adsorption that governs the protein-PET system and to a better understanding of the factors that govern this adsorption can be also foreseen.

MATERIALS AND METHODS

The microchip device used in this work was described elsewhere [START_REF] Gamby | Supercapacitive Admittance Tomoscopy[END_REF][START_REF] Gamby | Nanomosaic network for the detection of proteins without direct electrical contact[END_REF]. Some modifications have been made and will be outlined hereafter. The set-up is displayed in Fig. 1. Impedance measurements were carried out through a photoablated polyethylene terephthalate (PET) microchannel having a trapezoidal cross-section shape with a depth of 45 µm, a top width of 100 µm and a length of 1.4 cm. The separation distance between both microchannels was 200 µm center to center and 120 µm edge to edge. The electrode fabrication was achieved using a carbon ink loaded with gold nanoparticles, thermally laminated at 135°C and a pressure of 2 bar by a polyethylene (10 µm)/polyethylene terephthalate (25 µm) (PE/PET) layer. The total thickness of the laminated layer is equal to 35 µm. The separation distance in the PET band between the two planar microelectrodes and the main microchannel was equal to 5 µm and the detection surface area per microelectrode was 80 µm x 100 µm.

Electrochemical impedance spectroscopy measurements were performed by applying an ac voltage with constant amplitude (0.1 V) through the microelectrodes with a frequency ranging from 1 MHz to 1 Hz. The current measured is thus related to the total impedance of the sensor. These measurements were performed using a frequency response analyzer (FRA 1255B, Solartron U.K.) together with a dielectric interface 1296 (Solartron, U.K) which extends the frequency range from 1 MHz to 10 Hz. Data analysis was performed using the commercial software SMaRT. The effect of NaCl solutions on the resistance in the microchannel was studied at various concentrations (0 -10 -2 M). To analyze these results, a physical model was used as is described in Fig. 2. The data obtained were treated using a home-made software SIMAD which allowed fitting both frequency-dependent analytical expressions and equivalent electrical circuits.

RESULTS AND DISCUSSION

The results of impedance measurements between the two band microelectrodes depend on the surface state of the flow channel and on adsorption taking place on the walls of the flow channel.

However, in the previous study [START_REF] Gamby | Supercapacitive Admittance Tomoscopy[END_REF], the impedance between the microelectrodes and the microchannel through the bulk of the PET substrate was not considered since it has been reported that the resistance is large enough that it could be disregarded. Thus, impedance measurements performed here are undertaken on the device when it is empty or filled with a sodium chloride solution in order to reach a better understanding of the PET impedance contribution in the global response.

EIS measurements in an empty microchannel

Polyethylene terephthalate is a dielectric material [START_REF] Schaumburg | On the Accuracy of dielectric measurements[END_REF]. Its capacitive behavior between the two microelectrodes can be evaluated using

2 2 0 , 2 d S C r PET    (1)
where C 2,PET is the geometric capacity in the PET thickness, ε 0 is the vacuum permittivity (8.854•10 -12 F m -1 ), ε r is the relative permittivity of the PET film (3.4), d 2 is the average distance separation between the two electrode edge to edge (120 µm) and S 2 is the vertical surface of the microelectrodes (2•10 -7 m 2 ). This derivation gives a value of 4.5•10 -14 (F). Equation 1

  PET PET Q j Z PET , 2 , 2 , 2 1 ) (     (2)
where Q 2,PET is the CPE constant (Ω -1 s α ) and α 2,PET is the CPE exponent.

As described by Orazem et al. [START_REF] Orazem | Enhanced Graphical Representation of Electrochemical Impedance Data[END_REF], the log-log representation of the imaginary part of impedance vs frequency is suitable for a CPE analysis. The α 2,PET value corresponding to the slope of the curve plotted in fig 3d is 0.99. The extracted value of the CPE element (Q 2,PET ) in the PET thickness of 120 µm is estimated to 8.5•10 -13 (Ω -1 s α ).

The slight difference in value between the CPE element Q 2,PET and C 2,PET is due to the capacitive effect of the horizontal surface of the microelectrodes as mentioned above. The average value of 0.85 pF value could be considered to be the overall background (or stray) capacitance.

In the following, the impedance Z 2,PET will be determined by using equation 2 with α 2,PET = 0.99 and Q 2,PET = 8.5•10 -13 Ω -1 s α .

From the knowledge for an empty microchannel Z 2 (ω), the impedance in the microchannel (Z 1 (ω) in Fig. 2) can be deduced from the global impedance measurement Z G (ω). Then in the following work, all results will be represented by the impedance Z 1 (ω).

EIS measurements in microchannel filled with NaCl solution:

The global impedance Z G (ω) between the two microelectrodes, in the case of the microchannel being filled with a solution, depends on the different parts which constitute the microchip. In this configuration, Z G (ω) can be expressed as the combination between Z 1 (ω) and Z 2 (ω) (Eq. 3):

) ( 1 ) ( 1 ) ( 1 2 1    Z Z Z G   (3)
where Z 1 (ω) corresponds to a combination of impedances in the microchannel with the microchannel interface/PET and the impedance of the 5-µm thick. The impedance

Z 1 (ω) is equal to ) ( 2 1 ) ( 2 ) ( int , 1 1     Z R C j R Z Z S S S PET     (4)
In the microchannel: R S with units of  denotes the solution resistance in the microchannel which depends on the solution conductivity and on the geometry of microchannel, and C S in F represents the solution capacitance, which depends on its relative permittivity,  r,solution , usually equal to water permittivity in the case of a dilute solutions ( r,water = 78). In a typical measuring cell (d = 1 cm), the distance between the two electrodes is large enough to neglect the capacitance of the solution (C S ). However, this is not the case in the present microchip configuration, where the distance between the microelectrodes is only 120 µm.

In the PET layer of 5 µm: the geometric capacitance calculated by equation 5 for the PET layer of 5 µm (C 1,PET ) is equal to 1.4•10 -14 F. This capacitance is depends on the PET permittivity, PET thickness (d 1 = 5•10 -6 m), and surface contact area (S 1 = 8•10 -9 m 2 ), (see in Fig. 1):

1 1 0 , 1 d S C r PET    (5)
Taking into account the results obtained in the case of the PET thickness of 120 µm, the impedance for the PET thickness of 5 µm is given by

  PET PET Q j Z PET , 1 , 1 , 1 1 ) (     (6)
where Q 1,PET in  1 s  is the modeled capacitive behavior of the PET thickness of 5 µm and α 1,PET is the CPE exponent.

At the PET/ microchannel interface:

The impedance in the PET/ microchannel interface has a CPE behavior given by

  int int int 1 ) ( Q j Z     (7)
where Q int is the CPE constant at the interface and α int is the CPE exponent.

Experimentally, Z 1 (ω) is obtained from equation 3 where Z 2 (ω) is the experimental values of the impedance for an empty microchannel and Z G (ω) is the experimental values of the impedance for a filled microchannel with sodium chloride solution.

The real and imaginary parts of Z 1 () are calculated from

             2 2 2 2 2 2 2 2 2 2 1 " " ' ' ' " ' " " " ' ' " " ' ' ' G G G G G G G G Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z          (8) and               2 2 2 2 2 2 2 2 2 2 1 " " ' ' " " ' ' " " ' ' ' " ' " " G G G G G G G G Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z          (9)
where The highest resistance corresponding to the largest diameter loop obtained was recorded in the case of water (water conductivity measured: 3 μS cm -1 ); whereas, the lowest value was obtained for the 10 -2 M NaCl solution (conductivity measured: 1194 μS cm -1 ). This effect is also clearly shown in Fig. 4c. For a frequency range from 10 kHz to 1 Hz, the real part of the measured impedance decreased linearly with the microchannel conductivity.

Z' G , Z " G , Z' 1 , Z " 1 , Z'
The evolution of the log-log plot representation of the imaginary part (Fig. 4c) shows that, in the high frequency range from 1 MHz to 50 kHz (filled with 10 -2 M NaCl solution) and from 1 MHz to 1 kHz (filled with water), all slopes are estimated to 0.99, which is close to a pure capacitive behavior. Moreover, the cut-off frequency varies and depends on solution conductivity.

It is important to stress that, at the intermediate frequency (5 kHz to 300 Hz for NaCl 10 -2 M and 1 kHz to 300 Hz for water), the slope of imaginary part was found equal to 0.5. This behavior is characteristic of diffusion impedance or of the impedance of a rough interface.

The diffusion impedance analysis could be complicated in case of anomalous behavior [START_REF] Kumar | Theory of Generalized Gerischer Admittance of Realistic Fractal Electrode[END_REF][START_REF] Sunde | Impedance analysis of nanostructured iridium oxide electrocatalysts[END_REF][START_REF] Batalioto | Dielectric Dispersion of Water in the Frequency Range from 10 mHz to 30 MHz[END_REF][START_REF] Lenzi | Anomalous-diffusion approach applied to the electrical response of water[END_REF].

But in the present case, the hypothesis of diffusion impedance could be disregarded due to the fact that the microelectrodes were not in contact with electrolyte and the current was not faradaic.

It is important to notice that diffusion impedance is a part of faradaic impedance. Only resistivity, permittivity and species adsorption may be considered for the present system, and, as a consequence, only the interfacial impedance of rough surface was used in the impedance analysis. As a proof of that, scanning electron microscopy (SEM) micrographs of the photoablated PET surface presented in Fig. 5 highlight a microstructured surface constituted by the well-known microhills and microvalleys region on the microchannel surface due to photoablation of crystalline and amorphous regions, respectively [START_REF] Lippert | Chemical and Spectroscopic Aspects of Polymer Ablation: Special Features and Novel Directions[END_REF][START_REF] Gamby | Electroacoustic miniaturized DNA-biosensor[END_REF][START_REF] Gamby | Electroacoustic Polymer Microchip as an Alternative to Quartz Crystal Microbalance for Biosensor Development[END_REF]. This observation supports the conclusion that the interfacial impedance including the photoablated PET microchannel should not be represented as a pure capacitance. To take into account this characteristic, the De Levie model [START_REF] Levie | Advances in electrochemistry and electrochemical engineering[END_REF] for rough surface was used and adapted by different authors [START_REF] Barcia | Application of the impedance model of de Levie for the characterization of porous electrodes[END_REF][START_REF] Candy | The characterization of porous electrodes by impedance measurements[END_REF][START_REF] Fievet | Impedance measurements for determination of pore texture of a carbon membrane[END_REF] in order to determine impedance in rough or porous surfaces. In our case, this model could be also adapted for a microstructured surface.

In order to study the influence of the CPE exponent (α parameter) and to validate the hypothesis that the De Levie model applies in the present case, the Z 1 (ω) impedance model was calculated for different values of the α CPE parameter. A comparison between the experimental data and Z 1 (ω) impedance model is presented in Figure 6 with the assumed value of α CPE as a parameter.

The value of α CPE has a strong influence on the impedance in a frequency range that extends from 1 MHz to 1 Hz for the imaginary and real parts. The best correspondence between model and experiment was found with α CPE = 0.5.

As demonstrated by Keiser et al [START_REF] Keiser | Abschätzung der porenstruktur poröser elektroden aus impedanzmessungen[END_REF], for any pore shape, the impedance of a rough surface at high frequencies approaches 45° in the Nyquist plot (corresponding to a slope of 0.5 in the loglog representation of the imaginary part of the impedance as a function of frequency). The geometry of the microstructured surface is close to the shape numbered 5 in reference [START_REF] Keiser | Abschätzung der porenstruktur poröser elektroden aus impedanzmessungen[END_REF] and can be seen to be a good approximation for a photoablated microchannel. If the impedance corresponds to the high frequency model in the entire frequency range, its representation is a CPE as in equation 7 with α CPE = 0.5. This indicates that the CPE element for a photoablated surface microchannel is proportional to the square root of the capacitance of a non photoablated one.

Then, we shall have to determine the value of Q int for NaCl at various concentrations.

In the Nyquist diagram, we observed at low frequencies a capacitive branch due to the increase of the imaginary part of the impedance for frequencies lower than 100 Hz. This behavior is influenced by the electrical conductivity of the microchannel solution as displayed in Fig. 4.

Referring to the proposed electrical circuit, this behavior defines the interface PET/ microchannel. In the knowledge that the PET surface has functional groups that can form electrostatic interactions with ionic species, it is obvious that high concentrations of NaCl increase the capacitive effect recorded at low frequencies. The use of this impedance model in Eq. 4 provides an excellent fit to the EIS spectra shown in Fig. 4.

A Set of fitted parameters is reported in table 1. The model provided a very good fit to the experimental data, as was confirmed by different graphical representations of the model and data

values. The values of the CPE element (Q 1,PET ) obtained for water and for a NaCl solution displayed in table 1 show a slight difference that can be related to the microchannel composition.

Moreover, the fitted value of Q 1,PET is much higher than the geometric capacitance C 1,PET .

Conversely to the slight difference observed between Q 2,PET and C 2,PET , in the case of CPE element Q 1,PET , the difference arose from the value of the surface contact area S 1 which was underestimated, particularly on the side of the microchannel. The estimated value of C 1,PET could be seen as a smaller value due to the fact that only the vertical streamlines from the S 1 to microchannel section displayed in fig. 1 were considered.

The R S values obtained with the fitting procedure show that there is a relationship between resistivity and solution resistance, given as

 S R K  (10)
where K in m 1 is the cell constant microchannel.

A plot of the solution resistivity, , against the length of loop diameter, R S , gives rise to straight lines with a slope that is comparable to a cell constant, κ . The cell constant was found to be equal to 4.93•10 5 m -1 (see Fig. 7). This latter observation supports the model used in this paper.

Conclusion

In this study, we have demonstrated that contact-free impedance measurements can be successfully implemented in an insulated polymer for flexible microchips. The impedance technique has the advantage that it can carry out in situ characterization to obtain several physical parameters. Moreover, a new approach involving a new data processing for the analysis of the contact-free impedance through two parallel planar microelectrodes is demonstrated. It is pointed out that the resistivity through the PET layer cannot be neglected and must therefore be taken into account in the global microchip response. Indeed, the photoablated PET surface should be considered to be a microstructured surface, characterized by a constant phase element with a typical exponent of 0.5. The impedance model proposed here is inspired by the classical equation established by De Levie which leads to accurate and efficient fitting of our experimental data.

This study allowed a superior understanding of the global response of the PET microchip impedance. These steps are crucial in order to choose and optimize the experimental conditions for biosensor development in any kind of flexible polymer microchips [START_REF] Miserere | Fabrication of thermoplastics chips through lamination based techniques[END_REF].

  does not take into account the horizontal part of the two microelectrodes shown in figure 1. Impedance measurements in the case of empty microchannel in the frequency range from 1 MHz to 1 Hz and with 0.1 V are presented in figure 3. The Nyquist diagram shows a capacitive behavior over the entire frequency range, but the representative line is not exactly vertical. It is interesting to stress that the conventional RC circuit is not convenient for modelling impedance response through the dielectric PET layer. The impedance ) ( 2  Z between the two microelectrodes for the empty microchannel has a capacitance behavior represented by constant phase element (CPE), defined as

  2 and Z " 2 denote the real and imaginary parts of the impedances Z G (ω), Z 1 (ω) and Z 2 (ω), respectively. Electrochemical impedance spectroscopy measurements with a voltage of 0.1 V and frequency ranging from 1 MHz to 1 Hz were performed for microchannels filled with water and various concentrations of NaCl (5•10 -4 M to 10 -2 M). The goal was to extract the values of the solution resistance in the microchannel included in the Z 1 (ω) expression by varying the solution conductivity in the microchannel.The corrected impedance diagrams plotted by eliminating the PET layer impedance (of 120 µm) contribution from the global microdevice response are shown in Fig.4. Nyquist plots exhibit a loop from very high frequencies to medium frequencies (1 MHz to 100 Hz) and a capacitive behavior at low frequencies (see Figs.4a and b). One can observe that the diameter of the loop is inversely related to the microchannel conductivity.

Table 1 :

 1 Fitted parameters values according to the equivalent circuit described in Fig.2and results of Fig.4. The values of the resulting χ 2 statistic were obtained with errors of fitted values, σ = 0.01, for frequency ranging from 1 MHz to 1 Hz.

FigureFig. 1 :Fig. 3 :

 13 Figure

Fig. 4 :

 4 Fig. 4: Experimental impedance diagrams and modeling diagrams Z 1 (ω) in the microchip measured for an empty microchannel, in the presence of water in the microchannel and in the presence of NaCl solution at various concentrations, between 1 MHz to 1 Hz at 0.1V. (a) Nyquist plots. (b) Nyquist plots magnification. (c) Real part vs. frequency. (d) Imaginary part vs frequency.

Fig. 5 : 3  CPE = 1  CPE = 0. 9  CPE = 0. 8  CPE = 0. 7  CPE = 0. 6  CPE = 0. 5  3  CPE = 1  CPE = 0. 9  CPE = 0. 8  CPE = 0. 7  CPE = 0. 6  CPE = 0. 5  3  CPE = 1  CPE = 0. 9  CPE = 0. 8  CPE = 0. 7  CPE = 0. 6  CPE = 0. 5  3  CPE = 1  CPE = 0. 9  CPE = 0. 8  CPE = 0. 7  CPE = 0. 6  CPE = 0. 5 Fig. 6 :Fig. 7 :

 5319876531987653198765319876567 Fig. 5: (a) SEM micrograph of the photoablated PET microchannel. (b) Microstructured aspest of the bottom microchannel.

Acknowledgments

We thank Stephan Borensztajn and Françoise Pillier at "Laboratoire Interfaces et Systèmes Electrochimiques -Université Pierre & Marie Curie, UPMC Paris 6" for technical assistance on SEM micrographs acquisition. This research was supported by the French research organism "Institut de Recherche pour le développement" (IRD).


Olivier Roussel

Andy Borum

Michel Taïx

Timothy Bretl

Manipulation planning with contacts for an extensible elastic rod by sampling on the submanifold of static equilibrium configurations

We consider the manipulation planning problem of an extensible elastic rod in collision-free or contact space. We assume the rod can be handled by grippers either at both or at only one of its extremities and during the manipulation, the grasped end may change. We show that the use of both quasi-static and dynamic models can be coupled efficiently with sampling-based methods. By sampling directly on the submanifold of static equilibrium and contact-free configurations, we can take advantage of the dynamic model to improve the exploration in the state space. We show the necessity of considering contacts for this type of problems with several simulation experiments on various scenarios.

I. INTRODUCTION

Motion planning plays an essential role in assembling and disassembling industrial cases. Approaches have been proposed to the extension of this problem to robots manipulating movable rigid objects [START_REF] Siméon | Manipulation planning with probabilistic roadmaps[END_REF]. Moreover, in the automotive or aeronautical industry, manipulating deformable parts is necessary. In this context, most of these consist in circular Deformable Linear Objects (DLOs) (cable, hose, pipe,...) and are usually referred as rods in the mechanical literature. In the classical motion planning problem formulation, the objective is to compute a path in the collision-free space, without contacts, but in realistic conditions for this context the contact is actually a necessary condition.

In this paper, we consider the manipulation planning problem of a rod in free or contact space (as in Figure 1). We assume the rod is handled by a gripper at one of its extremities and the grasped end may change during the manipulating phase. As rod can be geometrically interpreted as an infinite-dimensional continuous curve, planning in the discretization of this curve may lead to high finitedimensional configuration spaces. Furthermore, identifying the manifold of feasible configurations, i.e. which satisfy mechanical constraints, on this configuration space is a challenging task. This is even more difficult if also considering both dynamics of the rod and contacts. Some work has been done using the reasonable assumption of considering only the collision free space of quasi-static configurations. These approaches are based on the numerical minimization of the total elastic energy for given gripper placements [START_REF] Lamiraux | Planning paths for elastic objects under manipulation constraints[END_REF], [START_REF] Moll | Path planning for deformable linear objects[END_REF], [START_REF] Wakamatsu | Static modeling of linear object deformation based on differential geometry[END_REF].

An approach to the manipulation planning considering dynamics of the rod and contacts is to plan in the finite space 1 CNRS, LAAS, Univ. de Toulouse, UPS, 7 avenue du Colonel Roche, F-31400 Toulouse, France, oroussel@laas.fr, taix@laas.fr 2 Department of Aerospace Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA, borum2@illinois.edu, tbretl@illinois.edu of controls for a given number of grippers. In this case, the state transition function is assumed to be known and can be delegated, for example, to a simulator. For a physically realistic rod model, these might use complex finite elements methods resulting in computationally expensive local planning schemes. For instance, the motion planning problem for deformable objects has been addressed in [START_REF] Rodríguez | Planning motion in completely deformable environments[END_REF] which coupled a deformable dynamics model with a kinodynamic motion planning algorithm. However, the use of fully deformable environments prevents the robot to be stuck in local minima and bypasses the local control problem.

In [START_REF] Saha | Manipulation planning for deformable linear objects[END_REF], the manipulation planning problem of a DLO is adressed from a topological perspective by focusing on knot tying application with the help of many grippers.

More generally, the motion planning problem for deformable objects has already been investigated, especially in the case of simplified visually realistic deformation models. For example, [START_REF] Bayazit | Probabilistic roadmap motion planning for deformable objects[END_REF] extendeded the Probabilistic Roadmap Methods (PRMs) [START_REF] Amato | A randomized roadmap method for path and manipulation planning[END_REF] for deformable objects by reducing the deformation space to a one dimensional one. [START_REF] Gayle | Path planning for deformable robots in complex environments[END_REF] used the Constraint Based Motion Planning framework to simulate a deformable robot along a guide path computed for a pointlike robot. This work has then been extended to the specific class of DLOs in [START_REF] Kabul | Cable route planning in complex environments using constrained sampling[END_REF].

Recent results based on the local solution of a geometric optimal control problem enabled to define this configuration space [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF], where it is shown this space defines a finite dimensional manifold that can be parameterized by a single chart. In previous work [START_REF] Roussel | Efficient motion planning for quasistatic elastic rods using geometry neighborhood approximation[END_REF], we presented a motion planning algorithm for a quasi-static inextensible elastic rod in complex environments and we showed that quasi-static solution paths in the collision-free space can be computed efficiently by assuming the rod as handled by both extremities. However, slight changes in grippers wrenches involve a change of the whole rod shape which causes greater probability of collision. Consequently, in constrained environments, the roadmap suffers from a slow growth of rate.

On the other hand, the physical engine eXtended Dynamic Engine (XDE) offers a realistic multi-body dynamics simulation with contact for deformable bodies. Thus, XDE can compute dynamic motions but at the cost of a high computational time, conflicting with random motion planning needs. However, it seems interesting to allow sliding along the contact space to improve significantly the exploration.

In this context, we can define a manipulation planning approach based on sampling-based methods. By sampling on the submanifold of quasi-static contact-free configurations, it is then possible to use dynamic simulation to extend the exploration of the state space. Using a parameterization to sample directly on the submanifold, we can limit the number of calls to the simulator to extension steps and to take advantage of provided information such as contact forces. By allowing contacts between the deformable object and obstacles, we will show we increase the deformation space and we can efficiently handle constrained environments, i.e. having a very poor -goodness [START_REF] Kavraki | Randomized query processing in robot path planning (extended abstract)[END_REF], and narrow passages by sliding along the contact space.

As we use two different models for dynamic simulation and static equilibrium configurations, it is necessary to fit models parameters in order to ensure their convergence at quasi-static configurations.

The paper is organized as follows. Section II presents both models used for dynamic simulation or static equilibrium configurations. Section III describes how to take advantage of both models using sampling-based methods to solve the motion planning problem for an extensible elastic rod. Finally, we evaluate in Section IV the method on several planning scenarios and discuss about the approach.

II. ELASTIC RODS MODELING

A. Extensible elastic rods Static Equilibrium Model (SEM)

The work from Bretl and McCarthy [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF] offers a single global chart to describe the manifold of equilibrium configurations of an elastic rod. As coordinates in this chart are a subset of a low dimensional Euclidean space, it is especially well suited for sampling-based methods. However, this model only considers inextensible and non-shearable elastic rods. This section presents here an extension of [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF] that enables us to derive a global parameterization for extensible and shearable elastic rods.

Consider a thin, naturally straight elastic rod in static equilibrium (i.e. the rod is motionless) held fixed at each end by robotic grippers. We may assume without loss of generality that the rod has unit length. With t ∈ [0, 1] denoting arclength along the rod, the position and orientation of the rod at arc-length t are described by an element q(t) of the special Euclidean group SE(3). The shape of the entire rod is described by a continuous map q : [0, 1] → SE(3). The rod is allowed to bend and twist, and the rod experiences axial extension and shear deformations. Therefore, the map

q : [0, 1] → SE(3) must satisfy q = q 6 i=1 u i X i + X 4 (1)
for some u : [0, 1] → R 6 , where X i with i ∈ {1, . . . , 6} is a basis for the Lie algebra se(3).

The functions u i : [0, 1] → R with i ∈ {1, . . . , 6} are strains along the rod, where u 1 is the twisting strain, u 2 and u 3 are the bending strains, u 4 is the axial strain, and u 5 and u 6 are the shear strains. The functions q and u are referred to together as (q, u).

The ends of the rod q(0) and q(1) are held by robotic grippers. Denote the set of all possible q(0) by B 0 = SE(3) and the set of all possible q(1) by B 1 = SE(3), i.e. B 0 and B 1 are the set of all possible placements of the robotic grippers holding the ends of the rod. Static equilibrium configurations of the rod are those which locally minimize the elastic potential energy stored in the rod. Assuming the elastic potential energy is quadratic in the strain functions and that twist, bending, stretching, and shear deformations are not coupled, the elastic potential energy is

1 2 1 0 6 i=1 c i u 2 i dt
where c 1 > 0 is the torsional stiffness, c 2 > 0 and c 3 > 0 are the bending stiffnesses, c 4 > 0 is the axial stiffness, and c 5 > 0 and c 6 > 0 are the shear stiffnesses of the rod. We say that (q, u) is a static equilibrium configuration of the rod if it is a local optimum of minimize q,u 1 2

1 0 6 i=1 c i u 2 i dt subject to q = q 6 i=1 u i X i + X 4 q(0) = b 0 , q(1) = b 1 (2)
for some b 0 ∈ B 0 and b 1 ∈ B 1 .

A similar optimal control formulation was used by Bretl and McCarthy to model an unshearable and inextensible elastic rod [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF]. We now apply the procedure for deriving necessary and sufficient conditions for optimality outlined by [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF] to the optimal control problem (2). Define the constants 6}. Theorem 1 provides necessary conditions for (q, u) to be a static equilibrium configuration of the rod.

c ij = c -1 i -c -1 j for i, j ∈ {1, . . . ,
Theorem 1: Define A = R 6 . A trajectory (q, u) is a normal extremal of (2) for b 0 = q(0) and b 1 = q(1) if and only if there exists µ :

[0, 1] → R 6 satisfying μ1 = c 32 µ 2 µ 3 + c 65 µ 5 µ 6 μ2 = µ 6 + c 13 µ 1 µ 3 + c 46 µ 4 µ 6 μ3 = -µ 5 + c 21 µ 2 µ 1 + c 54 µ 4 µ 5 μ4 = c -1 3 µ 3 µ 5 -c -1 2 µ 2 µ 6 μ5 = c -1 1 µ 1 µ 6 -c -1 3 µ 3 µ 4 μ6 = c -1 2 µ 2 µ 4 -c -1 1 µ 1 µ 5 (3)
u i = c -1 i µ i for i ∈ {1, . . . , 6} (4)
and ( 1) with the initial condition µ(0) = a for some a ∈ A. Theorem 1 results from applying Theorem 3 of [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF] to the optimal control problem (2). Theorem 3 of [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF] relies on a geometric statement of the Pontryagin maximum principle, and the function µ : [0, 1] → R 6 is the costate trajectory provided by the maximum principle. The function µ can also be interpreted as the vector of internal forces and torques acting along the rod [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF].

Theorem 1 provides conditions under which (q, u) is a normal extremal of (2). Let C denote the set of all smooth maps (q, u) : [0, 1] → SE(3) × R 6 which satisfy the necessary conditions given in Theorem 1. Then any (q, u) ∈ C is completely defined by the choice of a ∈ A, as is the corresponding µ. Denote the resulting maps by Ψ(a) = (q, u) and Γ(a) = µ.

The trajectories (q, u) ∈ C are extrema of the optimal control problem (2). Theorem 2 provides a computational test to determine which of these extrema are local minima. This test relies on the theory of conjugate points in optimal control problems [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF]. This defines the matrix G by

G = diag(c -1 1 , c -1 2 , c -1 3 , c -1 4 , c -1 5 , c -1 6
) Also define the matrix functions F and H : R 6 → R 6×6 by

F(µ) =     0 c32µ3 c32µ2 0 c65µ6 c65µ5 c13µ3 0 c13µ1 c46µ6 0 1+c46µ4 c21µ2 c21µ1 0 c54µ5 -1+c54µ4 0 0 -c -1 2 µ6 c -1 3 µ5 0 c -1 3 µ3 -c -1 2 µ2 c -1 1 µ6 0 -c -1 3 µ4 -c -1 3 µ3 0 -1 1 µ1 -c -1 1 µ5 c -1 2 µ4 0 c -1 2 µ2 -c -1 1 µ1 0     H(µ) =       0 c -1 3 µ3 -c -1 2 µ2 0 0 0 -c -1 3 µ3 0 c -1 1 µ1 0 0 0 c -1 2 µ2 -c -1 1 µ1 0 0 0 0 0 c -1 6 µ6 -c -1 5 µ5 0 c -1 3 µ3 -c -1 2 µ2 -c -1 6 µ6 0 1+c -1 4 µ4 -c -1 3 µ3 0 c -1 1 µ1 c -1 5 µ5 -1-c -1 4 µ4 0 c -1 2 µ2 -c -1 1 µ1 0      
Theorem 2: Let (q, u) = Ψ(a) and µ = Γ(a) for some a ∈ A. Let M and J : [0, 1] → R 6×6 be the solutions of the linear, time-varying matrix differential equations

dM dt = F(µ(t))M dJ dt = GM + H(µ(t))J (5)
with initial conditions M(0) = I and J(0) = 0. Then, (q, u) is a local optimum of (2) for b 0 = q(0) and b 1 = q(1) if and only if det(J(t)) = 0 for all t ∈ (0, 1]. Theorem 2 results from applying Theorem 4 of [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF] to the optimal control problem (2). This Theorem provides a test to determine which extremals provided by Theorem 1 are local optima of (2), i.e which a ∈ A produce local optima Ψ(a) ∈ C. Let A stable ⊂ A be the subset of all a ∈ A which satisfy the conditions in Theorem 2, and let C stable = Ψ(A stable ) ⊂ C.

We now have a characterization of the set of all static equilibrium configurations of the rod. This allows us to use a sampling-based planning algorithm in which we directly sample static equilibrium configurations. This is done by sampling points in A = R 6 , using the test in Theorem 2 to check if these points are members of A stable , and then evaluating the map Ψ| Astable : A stable → C stable at these points.

B. The physical engine XDE

XDE offers a realistic multi-body dynamics simulation with various kinematics constraints (e.g. joints, kinematic loops) and with real-time performances. In addition to rigid bodies and kinematic chains, it can also handle DLOs, modeled as geometrically exact 3D beams [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF]. This model enables large displacements thanks to Reissner kinematics and uses geometrically exact finite elements.

Furthermore, XDE can handle non-smooth contacts accurately with friction with constraints based methods. Also, it provides smooth bodies interactions mechanisms as Proportional-Derivative coupling for a given body position. This will be used for controlling rod grippers in simulation in the local planning method, as detailed in section III-B.

C. Model fitting between the SEM and XDE

For the purposes of this paper, we consider rods that are made from a homogeneous isotropic linear elastic material. As is common in mechanical engineering, XDE defines rod material elasticity parameters through Young's modulus E and shear modulus G. SEM stiffness coefficients c i can then be deduced from elasticity parameters by c 1 = GJ, c 2 = c 3 = EI, c 4 = EA and c 5 = c 6 = GA where I represents the second moment of area, J the polar moment of inertia and A the cross-section area of the rod.

Using these parameters, we can bring simulated dynamic rods to quasi-static configurations described by the SEM by fitting rod grippers positions. As there exists a countable number of quasi-static configurations for grippers placements, this method does not give us the guarantee to fall into the desired local minimum of elastic energy. However, this will be sufficient to our needs and for keeping good exploration properties as explained in section IV-C. Attached video shows simulation results for models convergence.

III. MOTION PLANNING FOR EXTENSIBLE ELASTIC RODS

In this section, we show how the use of both models can improve the efficiency of solving our manipulation planning problem. By taking advantage of the parameterization provided by the SEM, we can sample quasi-static contact-free configurations that guide the exploration in the space of all admissible cable states reachable by the simulator.

A. Problem formulation

Consider a rod parameterized by t ∈ [0, 1], its configuration q(t) can be represented by the mapping q : [0, 1] → SE(3). The resulting configuration space would be a subset of the infinite-dimensional function space C ∞ ([0, 1], SE(3)). In our case, the rod is discretized by the simulator using finite elements into N -1 elements (thus N nodes), thus the resulting configuration space is a subset of SE(3) N , where each configuration q is given by q = q 1 q 2 ... q N T .

However, as we consider rod dynamics provided by the simulator, the phase space of the rod must be explored. Let X ⊂ (T SE(3)) N , i.e. the tangent bundle associated to each node configuration, be the state space of all rod states x satisfying elasticity constraints defined by x = q q T .

Let now X adm be the closed set of all admissible states that can be reached by the dynamic simulator, which includes the set of states with contacts, and X obs the closed set of all states in collision with an obstacle. Note that X adm = X \ int(X obs ). The simulator, which relies on a state transition function, must guarantee that from any initial state in X adm , all states resulting for any input will be also in X adm .

However, the main drawback of considering the space of states reachable by a dynamic simulator is that we do not have a direct way to sample states on X .

On the other hand, the result in section II-A has shown that all static equilibrium configurations of the rod can be parameterized using single chart mapping from a 6dimensional Euclidean space A stable for a fixed rod base q 1 . Let now consider a free-flying quasi-static rod, its configuration space can be expressed using the same chart by A stable = SE(3)×A stable . The 12-dimensional submanifold of X parameterized by this mapping is the open set Y that describes all free-flying quasi-static elastic rods held by its both extremities, i.e. without contacts along the rod. Alternatively, the set Y can be defined by

Y = {x ∈ int(X adm ) | q = 0}
Fig. 2. Illustration of the tree exploration in the state space X adm while sampling states on the submanifold of quasi-static contact-free states Y. We can see that the resulting state of an extension step can lie on this submanifold as x i , be in contact as x k or be any state in X adm as x j .

Note that initial x start and goal x goal states of the problem lie in Y. As shown in Figure 2, the key idea is to explore the space of all reachable states X adm while sampling states on the submanifold Y using coordinates in A stable .

B. Local planning using dynamics and contacts

In our context, the local planning problem consists in finding the trajectory from a given state x f rom ∈ X adm to a sampled state x to ∈ Y. There is various possible approaches and we describe in this subsection our controller that address that sub-problem.

We emphasize that taking advantage of contact information provided by the dynamic simulator by using contact sliding motions improves significantly the exploration of the search space for highly constrained cases (see section IV-C).

In this direction, we chose to exploit contact motions by not constraining the local planner trajectory on the submanifold Y. Our method, detailed in Algorithm 1, consists in a two steps approach as illustrated in Figure 3. First, the rod is manipulated by one of its gripper until the considered extremity reaches its goal position. Then, the other gripper is manipulated the same way while keeping the first one fixed.

Manipulation is done using a Proportional-Derivative controller on the grippers position. We also consider a time limit t max , chosen at random between given bounds, for which the controllers are being applied for a given local planning instance. Fig. 3. Illustration of the two steps approach for the local planner from a state x f rom to a quasi-static state xto. First, a randomly chosen gripper is manipulated to its respective goal (left). Once its goal is reached, the first gripper is fixed and the other one is manipulated the same way (right).

Algorithm 1 EXTEND(x f rom , x to ) Require: Time bounds [T low , T up ], simulation step time ∆t

1: Select at random first gripper g 1 2: Select at random a time limit t max ∈ [T low , T up ] 3: x new ← x f rom 4: while gripperDistance(x new , x to , g 1 ) > and t < t max do 5:

Apply position control on x new at g 1 for ∆t 6: end while 7: Let g 2 be other gripper 8: while gripperDistance(x new , x to , g 2 ) > and t < t max do 9:

Apply position control on x new at g 2 for ∆t 10: end while 11: return x new

C. Kinodynamic planning with sampling on a submanifold

Our approach to the manipulation planning problem for elastics rods is based on classical kynodynamic samplingbased approaches [START_REF] Lavalle | Randomized kinodynamic planning[END_REF] as detailed in Algorithm 2. Some work has been done in sampling-based planning while projecting configurations on a constraint manifold as sampling strategy [START_REF] Berenson | Manipulation planning on constraint manifolds[END_REF]. In our case, the complexity of the constraints based on the elasticity of the rod material make the formulation of a projector to the sampling submanifold a difficult task. However, the single global chart derived from the SEM that parameterize the submanifold Y provides us a direct way to sample on the manifold and in particular offers the guarantee that it can be covered entirely through random sampling.

The pseudo-metric considered in the NEAREST returns an approximation of ρ(a, b) is defined by

ρ(x i , x j ) = 1 0 p i (t) -p j (t) dt (6)
Algorithm 2 Kinodynamic RRT for an elastic rod with sampling on the quasi-static states submanifold Require: Environment model, x start , X goal 1: Initialize the tree T with x start 2: while ¬ solved and iter < N max do 3:

x rand ← random quasi-static state ∈ Y or goal state 4:

x near ← NEAREST(T , x rand ) 5:

x new ← EXTEND(x near , x rand ) 6:

Add current state x new to T 7:

Add edge (x near , x new ) 8: end while where p(t) ∈ R 3 is the translation part of the rod position at q(t) ∈ SE(3). Intuitively, this pseudo-metric represents the swept surface area in the workspace between the two rod states as if all nodes could move in straight line which gives us a lower bound to the real swept surface area.

IV. RESULTS AND DISCUSSION

In this section, we will present our motivations to use dynamics and contacts and we will show and discuss some experimental results obtained in simulation.

A. A critical case: the double funnel

Consider a rod of length l and radius r, and the case of an environment with a very narrow passage for the given length. In the extreme case, this narrow passage would tend to have the shape of the rod, i.e. here a tunnel of length l + and of radius 2r + as illustrated in Figure 4. Let X crt be the set of rod configurations that are collision free when the rod is centered in the tunnel. As tends to zero, the set X crt tends to be a unique state of with a zero volume. Consequently, the probability of sampling a state in X crt tends to zero when the tunnel size tends to fit perfectly the rod shape. This is an extreme narrow passage where sampling-based methods typically fail. The intuition behind taking advantage of contacts for this type of scenarios would be to:

• allow states to be in contact, i.e. x ∈ ∂X obs ,

• being able to slide along contact boundary ∂X obs ,

• being able to get off the contact. Considering this, we tested our previously described approach on a double funnel scenario, with a tunnel of length 2l and a radius 3r (see Figure 5). This type of case illustrate the advantage of our approach over classical collision-free sampling-based methods. On this scenario, a classical collision-free sampling-based method would typically fail (low probability to sample inside the tunnel). Our approach can solve the problem using dynamics/contacts with good performance results (see Tab. I).

B. Experimental results

A typical application to manipulation planning of elastic rods consists in assembly and disassembly studies of flexible parts such cables. In the case of disassembly study, the initial state is generally very constrained and the goal is a free region. Assembly studies can be solved by a symmetrical process. In addition to the double funnel case previously described, the following scenarios have been tested:

• engine (Figure 7) represents an industrial disassembly.

• grid (Figure 6) represents a disassembly operation for an elastic rod that is winded up into a grid. • powerplant (Figure 1) where both initial and goal states for the elastic rod show a poor -goodness. Solution paths on these scenarios are available in the attached video. Each scenario has been run 50 times with a limit in time of 20 minutes. All the benchmarks were run on a PC with 16GB of main memory and using one core of an Intel Core i7-2720QM processor running at 2.2GHz. Implementation has been done in C++ using the OMPL [START_REF] Sucan | The open motion planning library[END_REF].

Results of the manipulation planning on previously presented scenarios are shown in Table I.

Assembly and disassembly scenarios showing difficult narrow passages are solved efficiently using our approach. In addition, we provide a solution path with accurate rod dynamics thanks to the simulation based on finite element analysis. This solution includes the sequence of controls that should be applied to the rod to solve the problem, and this information could be used in realistic conditions. Note that the planner performance significantly decreases with the powerplant scenario, as it differs from typical assembly studies and consists in a complex and highly constrained case from initial to goal states.

C. Discussion about completeness

As our approach relies on sampling on a submanifold of the state space, our algorithm is probabilistically complete in the sampled submanifold Y. Indeed, our direct sampling method enables us to cover the entire submanifold Y of quasi-static contact-free configurations through the sampling the process. Furthermore, each infinitesimally small quasistatic motion that lies entirely in the submanifold Y can be achieved through dynamic simulation as both models fit.

Also, there exists a countable number of static equilibrium rod configurations for given grippers placement, the local planner may fail to reach the desired state by uniquely controlling grippers positions and fall into another local minimum of elastic energy. For the purpose of our local method, this does not invalidate the probabilistic completeness of the approach as there is a non null probability to reach the desired state. Moreover, these configurations are guaranteed to be close with respect to our metric, which enables us to keep a good Voronoi bias in Y. For this reason, we first limited our goal region X goal to the set of all static equilibrium configurations for given grippers positions.

V. CONCLUSION In this paper, we proposed an approach to the manipulation planning problem of an extensible elastic rod in collisionfree or contact space. We showed that by using a static equilibrium model, we can efficiently sample on the submanifold of its quasi-static contact-free states whereas the kinodynamic planning algorithm uses dynamic simulation and explore the full state space of the rod. Thanks to dynamic simulation, the planning algorithm can take advantage of contacts by allowing sliding motions in order to solve difficult narrow passages, where collision-free typical samplingbased method would fail.

Fig. 1 .

 1 Fig. 1. Solution path for an elastic rod to the "powerplant" scenario. Both sides (left) and top (right) views have been clipped for clarity.

Fig. 4 .

 4 Fig.4. A critical narrow passage for an elastic rod represented in the workspace W (left) and its equivalent illustrated in rod state space X . As tends to zero, the probability of sampling on Xcrt also tends to zero.

Fig. 5 .

 5 Fig. 5. The double funnel case where the elastic rod must go trough a very narrow passage (left) and a solution path (right) where the rod must go from the left side to the right side (view has been clipped for clarity).

Fig. 6 .

 6 Fig.6. Solution path on the "grid" scenario where the rod must go from a winded up state to a free state.

Fig. 7 .

 7 Fig. 7. Start (in green at clipped view a) and goal (in red at view b) states for the disassembly scenario "engine". The solution path makes use of contacts (views c and d).

TABLE I KINODYNAMIC

 I PLANNING FOR ELASTIC RODS WITH CONTACTS RESULTS

	Scenario	Number of vertices / faces	Success rate	Resolution time (sec) avg ± std. dev
	double funnel	198/388	100%	77.2 ± 36.0
	grid	704/1320	96%	34.4 ± 25.5
	engine	64885/131583	94%	119.6 ± 109.3
	powerplant	24870/20053	31%	402.5 ± 293.6

ACKNOWLEDGMENT

This work was supported by the French National Research Agency under the project Flecto (ANR-Digital Models). The engine industrial model is courtesy of Siemens-KineoCAM.