Arthur Marmin
email: arthur.marmin@centralesupelec.fr

Marc Castella
email: marc.castella@telecom-sudparis.eu

Jean-Christophe Pesquet
email: jean-christophe.pesquet@centralesupelec.fr

Laurent Duval
email: laurent.duval@ifpen.fr

Signal reconstruction from sub-sampled and nonlinearly distorted observations

come

Signal Reconstruction from Sub-sampled and Nonlinearly Distorted Observations Arthur Marmin * , Marc Castella † , Jean-Christophe Pesquet * , Laurent Duval ‡

I. INTRODUCTION

Accurate data acquisition of sparse signals from real-world measurements remains an open challenge. Our present motivation can be illustrated by an actual situation, encountered in analytical chemistry, especially in chromatography. Let us consider a mixture of chemical compounds in different concentrations. An optimal measurement would produce a set of isolated peaks whose locations and amplitudes characterize each compound and its relative concentration. A discrete version is illustrated with the three-sample parsimonious black-dot signal modeled in Figure 1. Usually, such a sparse information is subject to undesirable fluctuations and degradations, requiring to employ restoration methods from observed responses. Here, the linear part of the system assumes the convolution by a limited-support kernel. In Figure 1, the kernel is a five-sample discrete approximation of a binomial filter. Such models are commonly used in physico-chemical data processing [START_REF] Vendeuvre | Characterisation of middle-distillates by comprehensive two-dimensional gas chromatography (GC × GC): A powerful alternative for performing various standard analysis of middledistillates[END_REF]- [START_REF] Ning | Chromatogram baseline estimation and denoising using sparsity (BEADS)[END_REF]. The resulting peak signal is further affected by a nonlinear distortion, representing a saturation, that flattens higher amplitudes. This is depicted on the blue-cross signal. Those smoothed data are further degraded by an additive Gaussian noise as shown in solid red lines.

This paper extends results in [START_REF] Castella | Optimization of a Geman-McClure like criterion for sparse signal deconvolution[END_REF]- [START_REF]A global optimization approach for rational sparsity promoting criteria[END_REF] and provides a novel recovery method for sparse signals from sub-sampled observations obtained through a nonlinear and noisy model. This stands in contrast with the vast majority of available methods, which only apply in a linear context [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]- [START_REF] Soubies | A continuous exact 0 penalty (CEL0) for least squares regularized problem[END_REF]. Among the few methods that have tried to weaken the linearity assumption, one can mention [START_REF] Schetzen | Nonlinear system modelling and analysis from the Volterra and Wiener perspective[END_REF]- [START_REF] Deville | An overview of blind source separation methods for linear-quadratic and post-nonlinear mixtures[END_REF]. Our method builds on the minimization of a sum of rational functions, for which recent methodologies can be applied. Here, we specifically focus on the possibility to acquire the data at a reduced rate. This feature is important in applicative areas, since it permits faster acquisitions for high-throughput experiments and analysis. In addition, the structure induced by the sub-sampled acquisition reduces the overall complexity, which was a major difficulty in the previous work in [START_REF]A global optimization approach for rational sparsity promoting criteria[END_REF]. We provide an analysis of the benefit of the sub-sampling on the complexity of the problem.

Our paper is organized as follows: Section II introduces our model for the observed signal. Section III explains how our method recovers the target signal. It also provides a study of the computational difficulties encountered. Simulation results can be found in Section IV and Section V concludes our work.

We introduce the following notation: for any nonnegative integers n and p, S n + is the cone of n × n real positive semidefinite matrices, n p is the binomial coefficient "among n choose p" and • (resp. • ) is the greatest (resp. smallest) integer lower (resp. bigger) than its argument.

II. PROBLEM STATEMENT

A. Observation model

We consider the reconstruction from sub-sampled measurements of an unknown sparse discrete-time signal x of length T , comprising few peaks. The measurement process deteriorates x in the following way: peak enlargement, nonlinear distortion due to sensor saturation, and noise. Finally, the measured signal is sub-sampled during the acquisition. This is modelled by a convolution with a finite impulse response filter given by vector h followed by a nonlinear function Φ, the addition of a noise w with samples drawn from an i.i.d. zero-mean Gaussian distribution, and a decimation operator D. Defining the observation vector y of size U after sub-sampling, the corresponding modelling equation reads

y = D Φ(h * x) + w .
We will be mainly interested in regular decimations D α where all the elements indexed with a multiple of α are deleted, namely D α (v t ) 1≤t≤T = (v ∆(u,α) ) 1≤u≤U where ∆ is defined as ∆(u, α) = u + u-1 α-1 . We denote by D ∞ the identity operator that preserves the entire signal.

B. Signal recovery criterion

To estimate the original signal x, we minimize a criterion J composed of two terms. The first one is a fit measure between the recorded measurements y and the output of the noiseless model for a given estimate x of the original signal x. The second term is a sparsity-promoting penalization approximating the 0 pseudo-norm weighted by a positive parameter λ, i.e.

J (x) = y -D α (Φ(h * x)) 2 2 data fidelity + λ T t=1 Ψ(x t ) penalization .
In the following, we choose Φ and Ψ as rational functions and Φ acts component-wise. To simplify our notation, we set the components of x with nonpositive index to be identically zero. Then J can be written as a sum of rational functions

J (x) = U u=1 y u -Φ L l=1 h l x ∆(u,α)-l+1 2 + λ T t=1 Ψ(x t ) = U u=1 p u (x ∆(u,α)-L+1 , . . . , x ∆(u,α) ) q u (x ∆(u,α)-L+1 , . . . , x ∆(u,α) ) + T t=1 r(x t ) s(x t ) .
The polynomials (p u ) 1≤u≤U and (q u ) 1≤u≤U are functions of L variables whereas the polynomials r and s depend on one variable. We define also d q = degree(q u ) 2 and d s = degree(s) 2 .

III. METHOD

A. Lasserre's relaxation

We here extend the methodology developed in [START_REF]A global optimization approach for rational sparsity promoting criteria[END_REF] to the case of decimated observations. Using a recent approach introduced by Lasserre [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF] and related to the sum of squares relaxation, the solution to the rational optimization problem

J * = min x∈R T J (x) , (1)
can be found by solving a hierarchy of convex semi-definite programming (SDP) problems in their standard dual form. Defining a vector z whose components are related to the monomials up to degree 2k in the variables of Problem ( 1), the relaxation of order k takes the following form

P * k = min z∈R m b z (2) subject to C - m i=1 z i A i ∈ S n + (3) f -G z = 0 , (4)
where C and

(A i ) 1≤i≤m are symmetric matrices, b ∈ R m , f ∈ R n l , and G is a matrix of R m×n l .
The size of the SDP, expressed by the integers n, m and n l , is related to Problem (1) and to the order k of the relaxation. This is detailed in Section III-B, where it will appear that the size of the SDP grows tremendously as k increases. In addition, it has been proved [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF], [START_REF] Bugarin | Minimizing the sum of many rational functions[END_REF] that the SDP relaxations P * k converge to the optimum J * of the rational optimization problem when the relaxation order goes to infinity, that is lim k→∞ P * k = J * . Fortunately, low relaxation orders suffice in practice to provide correct solutions. The order of relaxation yet needs to be at least equal to half of the maximal degree of the polynomials involved in the criterion to allow all the polynomials in the original problem to be represented. Finally, the convergence results require the original variables to belong to a compact subset of R T . We thereby restrict the domain of our problem to a closed bounded box chosen such that a solution to Problem (1) belongs to it.

B. Structure and complexity of the relaxation

This section studies the structure and the size of the SDP relaxation. For simplicity, we do not take into account here that some variables in the expression of J are set identically to zero. The SDP size presented here will hence be slightly overestimated.

In the SDP relaxation, each rational fraction in J is associated with a measure and the components of the variable z are moments (truncated to the order 2k) of these measures. To ensure that each sub-sequence corresponding to a given rational term is a valid moment sequence for a set of say p variables, 1 + 2p positive semi-definite constraints must be enforced (one for the whole set of variables, and 2p for the two bound constraints coming from the lower and upper bounds imposed to each variable). Since our criterion is a sum of U rational functions over L variables and T rational functions over 1 variable, the total number of such semi-definite positive constraints is equal to

n s = U (1 + 2L) + 3T .
These constraints can be expressed more concisely under the form (3) provided that matrices (A i ) 1≤i≤m and C have a block diagonal structure involving n s blocks.

In addition, since any given variable x t may appear in multiple rational functions, many moments must be identical in z. These moment equality constraints correspond to the linear equality constraints in [START_REF] Castella | Optimization of a Geman-McClure like criterion for sparse signal deconvolution[END_REF]. Let us denote by θ u the overlap parameter, defined as the number of variables in common between pu+1 qu+1 and pu qu when u ∈ {1, . . . , U -1}. θ u depends on u, but also on L and α.

According to the general methodology [START_REF] Jibetean | Global optimization of rational functions: a semidefinite programming approach[END_REF], we have to consider the equality of the moments of all monomials of these θ u overlapping variables up to degree 2(k -d q ).

As a consequence, this gives θu+2(k-dq)

2(k-dq)
equality constraints for every u ∈ {1, . . . , U -1}. Similar arguments hold for the T rational functions r s appearing in J , but since it depends on one variable only, the above binomial coefficient reduces to 2(k -d s ) equality constraints on moments. Hence,

n l = U -1 u=1 θ u + 2(k -d q ) 2(k -d q ) + 2(k -d s )T .
Concerning the size m of the vector z, it is given by the sum of the number of monomials up to degree 2k for all rational terms in J . A total of U of these terms depend on L variables, and T terms depend on one variable. Noting that

1+2k 2k = 1 + 2k, it follows that m = U L + 2k 2k + T (1 + 2k) .
Furthermore, n is the sum of the block sizes in (A i ) 1≤i≤m and C. The blocks can be viewed as matrices whose rows and columns are indexed by monomials. Each of their entries is a sum of moments whose maximal degree is the product of the row index, the column index and, if any, a polynomial defining a bound constraint. For a maximal degree 2k, and for the first U terms in J the block sizes are thus L+k k for the whole set of variables and L+k-1 k-1 for the upper and lower bound constraints. The last T terms in J can be considered similarly. Finally, we obtain

n = U L + k k + 2L L + k -1 k -1 + T (1 + 3k).
The order of relaxation k appears as a binomial coefficient in both SDP size variables n and m, so that increasing k quickly blows up the size of the SDP. The previous remark can be extended to L. Actually, for a given order of relaxation k, the size of the SDP asymptotically becomes of the order:

m = O(U L 2k ), n = O(U L k ) . (5)
In our context, decimation yet reduces efficiently the size of the SDP relaxation since it decreases U and (θ u ) 1≤u≤U -1 .

Both play a prominent role in the values of n, m, n s , and n l which determine the complexity of the optimization problem. Table I compares the size of SDP relaxations without decimation (D ∞ ) and with the two-fold D 2 decimation. As discussed above, the size variables n and m increase quickly with the order of relaxation k and the length of the filter L. In the decimated case, the variables n and m are decreased about one half while n s and n l are respectively reduced about one third and two thirds.

C. Exact relaxation of 0

In this section, we show how to theoretically address an exact 0 penalty with our relaxation. We explain why this relaxation is intractable in practice by relying on our previous complexity estimation. This justifies our choice of a rational approximation of 0 as a regularization.

An equivalent form of the minimization of J , where Ψ is replaced by the exact 0 penalization is as follows minimize

(x,ξ)∈(R T ) 2 y -D α (Φ(h * (x ξ))) 2 + λ T i=1 ξ i subject to (∀i ∈ {1, . . . , T }) ξ i = ξ 2
i , where the operator denotes the element-wise Hadamard product and λ is a nonnegative regularization parameter. Unfortunately, in the above formulation, the introduction of the ξ i necessary to formulate the 0 penalization in a polynomial form results in twice more variables than in the original problem. As a consequence, the variables n and m defining the size of the SDP relaxation become much higher. A rough and underestimated idea of the SDP size is obtained by replacing L with 2L. Our method may then become intractable.

The substantial increase in the computation time has been confirmed in our experiments. Even for a signal of length 50 with a filter of length 2, we needed about three hours of computations to reconstruct the solution compared to a few seconds using our rational approximation of the 0 penalty.

IV. SIMULATIONS AND RESULTS

A. Simulation setup

We have run our simulations for nonnegative signals x with 50, 100 and 200 samples and for filters of length 3, 4 and 5 whose coefficients are binomial coefficients on the corresponding line of Pascal's triangle. These specific filter coefficients as well as the nonegativity of input signals are of practical importance for chemistry problems, especially in chromatography. Indeed the concentrations of compounds are always positive and the coefficients of the filters correspond to peak enlargement.

We have set the approximation Ψ of 0 to the following Geman-McClure-like potential [START_REF] Castella | Optimization of a Geman-McClure like criterion for sparse signal deconvolution[END_REF]

(∀t ∈ R) Ψ(t) = |t| δ + |t| ,
where the parameter δ controls the quality of the 0 approximation. Smaller δ leads to better approximation, but also to a harder optimization problem. We have set δ to 0.01. Note that this approximation of 0 is well-chosen because it is rational with low degree, which allows a low order of relaxation k.

The saturation function Φ acts component-wise and is defined as (∀t ∈ R) Φ(t) = t χ+|t| where the parameter χ has been set to 0.3 and controls the level of the saturation.

For each test, we have run 100 simulations and we show here the average results. We have set the order of relaxation to 3 and have used GloptiPoly [START_REF] Henrion | GloptiPoly 3: moments, optimization and semidefinite programming[END_REF] coupled with the SDP solver SDPT3-4 [START_REF] Tütüncü | Solving semidefinitequadratic-linear programs using SDPT3[END_REF] to solve the optimization problem. We have run the simulations on a standard computer with an Intel I7 CPU running at 3.60 GHz and 32 GB of RAM.

We have compared our method to an improved LASSO approach (named iLASSO). It consists of applying first the least absolute shrinkage and selection operator (LASSO) [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] using a linearization of Φ. More precisely

x lasso = argmin x∈[0,+∞[ T y -D α (χ -1 (h * x)) 2 + λ x 1 .
The solution x lasso is then used to initialize a proximal gradient algorithm based on iterative hard thresholding (IHT) [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF] extended to the nonlinear model.

The value of the relaxation parameter λ has been optimized empirically and set to 0.10 for our method and 0.05 for LASSO and IHT. Decimation D ∞ represents the case when the signal is not sub-sampled and the corresponding results are provided for comparison with results discussed in [START_REF]A global optimization approach for rational sparsity promoting criteria[END_REF].

B. Results

In this section, we study the effect of the different decimation patterns first on the reconstructed signal quality and then on the computational time.

1) Reconstruction quality: Figure 2 illustrates, from top to bottom, a sub-sampled observed signal y, the corresponding original signal x, the signal reconstructed with our method and the signal reconstructed using iLASSO. Within the framework of our nonlinear model, iLASSO performs poorly whereas our estimation of x is close to the original signal.

To assess peak detection, we introduce a varying threshold on the estimated amplitude of x. peaks but is also more prone to false detection. Conversely, our method detects a high number of true peaks for the different thresholds and make few mistakes. In accordance with our observations in Figure 2, our method has a very low false positive rate which maintains its ROC curve close to the ordinate axis. Moreover, the amplitudes of the peaks are sensibly closer to the real ones with our method than with iLASSO. Amplitudes are highly important in chromatography since they generally correspond to the concentration of the different compounds that we are analyzing.

Comparing the ROC of our method for the different decimation patterns in Figure 3, we notice that the higher the decimation is, the fewer peaks detected in the true position. We observe that our model keeps a very low false positive rate which is related to the fact that very few peaks are detected at a shifted position. Decimation mainly causes the system to miss some peaks and weakly affect their positions.

2) Computational load: Table II shows the computational time corresponding to the different test cases for the three decimation operators D ∞ , D 4 and D 2 .

The computational time explodes when the length of the filter increases. This is due to the presence of L in the binomial coefficients of the expressions of n and especially of m in [START_REF]Recovery of nonlinearly degraded sparse signals through rational optimization[END_REF]. Indeed, state-of-the-art SDP solvers exhibit high computational time when the dimension m is large.

However, the higher the decimation, the lower the computational cost. The effect is especially beneficial for high values of L and it allows us to deal with filter lengths that are intractable without decimation. This is of high practical interest, because filter lengths from 5 to 7 can be sufficient for dedicated chemistry needs.

V. CONCLUSION

We have considered the challenging problem of reconstructing sparse signals degraded by a decimated non-linear model. We have defined a least-squares fitting criterion penalized with an 0 approximation that can be written as a sum of rational functions. Thereby we have applied a recent methodology for rational minimization and relaxed the original problem into a hierarchical sequence of convex SDP. One of the main advantages of the proposed approach is that the SDP size is significantly reduced when the decimation factor increases, which may be of high interest in applications such as chromatography. The reduction of the computation time also allows us to consider filters with a larger length as encountered in signal processing applications. Finally, we have shown that our method compares favorably to LASSO followed by IHT.

Fig. 1 .

 1 Fig. 1. Example of chromatography signals

Fig. 2 .

 2 Fig. 2. Comparison between iLASSO and our method using D 2 decimation


Anne Atlan
email: anne.atlan@univ-rennes1.fr

Catherine Boyen

Laurence D Hurst

How

HOW TO LIVE WITHOUT SEX: THE ANSWER OF CYTOPLASMIC GENOMES

Keywords: organelle, sex evolution, mutation, Muller's ratchet

whether they are published or not. The documents may come

Introduction

The prevalence of sex is an enigma for evolutionary biologists. Because sexual reproduction is costly compared to asexual reproduction [START_REF] Smith | The evolution of sex[END_REF], it is a challenge to understand why most extant species are sexual. The debate around this problem has generated over 20 different hypothesis [START_REF] Kondrashov | Classification of hypotheses on the advantage of amphimixis[END_REF]. Indeed, there are so many good reasons why organisms ought to have sex [START_REF] Hurst | Recent advances in understanding of the evolution and maintenance of sex[END_REF], it is considered that nothing could survive for a long time without sex. This led Maynard [START_REF] Smith | Contemplating life without sex[END_REF] to describe as an "evolutionary scandal" the persistence of asexuality in Bdelloid rotifers, a whole taxon in which no male has ever been observed.

In the long term, asexual species are expected to go extinct because of the accumulation of deleterious mutations. Muller's ratchet [START_REF] Muller | The relation of recombination to mutational advance[END_REF]) is one formulation of this mutational meltdown. It is envisaged that in an asexual population the number of slightly deleterious mutations can only increase. Selection may remove the individuals with the lowest fitness but if the deleterious effect of a single mutation is relatively weak then selection will be unable to remove it. Over time the net accumulation of slightly deleterious mutations will have a large effect leading to the extinction of the species. By allowing recombination between two different genomes, sexuality permits the creation of individuals with fewer mutations than either of the parents and hence mutational meltdown is not inevitable. In a genome with clonal inheritance, however, weakly deleterious mutations progressively accumulate in an irreversible manner, just as the turn of a ratchet is irreversible.

For these reasons, non recombining genomes are often thought to be an evolutionary dead-end. In sexual species, nuclear genomes recombine regularly during meiosis. Similarly, prokaryotes, though not having sex in a conventional sense, do undergo related parasexual processes in which DNA from other cells is incorporated. Even viruses have the capacity to recombine with other viruses as well as with their host. In fact, very few species have been shown to be obligately clonal [START_REF] Judson | Ancient asexual scandals[END_REF]. The rarity of such species ensures that they can be viewed as evolutionary curiosities. One type of asexual genome is however so widespread that its case cannot simply be so easily dismissed. These genomes are the uniparentally inherited cytoplasmic genomes of eukaryotes. As most cytoplasmic genomes have uniparental inheritance they do not have the opportunity to recombine with the genome of a different cell. As pointed out by [START_REF] Gabriel | Müller's ratchet and mutational meltdowns[END_REF], if Muller's ratchet is really a problem for asexuals, the DNA from organelles should suffer from its effect. This raises two questions: is mutational meltdown a problem for these genomes and if so how do they manage to avoid the accumulation of deleterious mutations?

With the aim of answering these questions we review the structure and properties of the genomes of different types of mitochondria and chloroplasts. The parameters that we consider are the genome structure, copy number, and rate of substitution. We then argue that these parameters appear not be independent, and that the link between them may slow down the ratchet. Nevertheless, in some cases deleterious mutations do accumulate. Other mechanisms such as paternal contribution, the mode of transmission and RNA editing may also preserve the functionality of theses genomes.

Molecular and Structural Properties of Cytoplasmic Genomes

There exist a number of similarities between genomes of all organelles. Nearly all are circular DNA molecules (or at least have circular maps) and exist in multiple copies. Their endosymbiotic origin is now widely accepted (Gray, 1993); and phylogenetic analysis demonstrates that mitochondria are related to alpha-proteobacteria, while chloroplasts are related to cyanobacteria [START_REF] Loiseaux-De Goër | Plastid lineages[END_REF]. Contrary to what was previously thought, it is now thought that all mitochondria derive from the same endosymbiotic event [START_REF] Leblanc | Origin and evolution of mitochondria: what have we learnt from red algae?[END_REF][START_REF] Gray | Mitochondrial evolution[END_REF]. Genetic transfer to the nucleus has strongly reduced the size of the mitochondrial genome, with the smallest mitochondrial genome being observed within parasistic organisms. Because this implies that constraints on parasites' cytoplasmic genomes are relaxed, we will not incorporate these organelles in this review.

An increasing set of data on organelles from various independent eukaryotic lineages show that genome size, structure, organisation, and substitution rate may be very different, even among closely related taxa [START_REF] Leblanc | Origin and evolution of mitochondria: what have we learnt from red algae?[END_REF]. On the other hand, some similarities are striking e.g. the curious structural resemblance between the mitochondrial genome of oomycetes and the chloroplast genome of angiosperms. The amount of data available on the different lineages is, however, variable. Among animals, much more is known about the genomes of mammals than about those of ectotherms and invertebrates. Among plants, angiosperm cytoplasmic genomes have been extensively studied, but almost nothing is known about those of gymnosperms. Fungi and oomycetes are also well studied, but our knowledge of other lineages is very recent and limited to few species. Table 1 summarises the properties of the best known cytoplasmic genomes.

Table 1: Some properties of the best studied cytoplasmic genomes. Data for the number of copies come from concordant estimations given by specialists of each lineages. All other data come from references listed in the text. For the rate of point substitution, we arbitrary choose to give the value of 1 for the nucleus (which taking current averages roughly corresponds to 2 to 3 substitutions per site and per billion years, although no doubt there is variation between eukaryotes). a Up to 10 5 copies of the genome may be found in the oocytes. b The number of copies given here corresponds to copies of the "master circle", that is the molecule, often virtual, that contains at least one copy of each gene. c The estimation of the number of copies is highly variable depending on the tissue. It may vary from few dozen in the roots to hundred thousands in the green leafs. The number in the female gametes is unknown. d The number of copies is variable depending on the species. cd the substitution rate is smaller in the inverted repeat than in the single copy region. e yeasts are the only organisms present in this table that exhibit biparental inheritance of mitochondria.

Genome

Genome structure

Animal mitochondria: Mammalian mitochondria have small genomes with minimal variability in genome size. They vary from 16 Kb for the smallest to 18 Kb for the largest. Furthermore, this genome is remarkably compact [START_REF] Grivell | Small, beautiful and essential[END_REF][START_REF] Boore | Animal mitochondrial genome[END_REF]. There are no duplicated sequences and non-translated 5' and 3' sequences are virtually absent. For instance, in humans, no more than 87 of the 16,569 base pairs can be regarded as lying in intercistronic regions. Five of the reading frames even lack a stop codon and terminate with either UA or U. However because mRNA must be polyadenylated prior to translation these sequences become UAA and hence a stop codon is produced without having to be coded for. This is clearly a base saving device. The same explanation can be mounted to explain the fact that some genes have an overlap between them. Often the start of one gene is the termination base of the previous gene. This compactness of the genome strongly suggests that selection on genome size has been an important force. Mammalian mitochondrial DNA does not normally exhibit recombination and gene arrangement is highly conserved [START_REF] Boore | Animal mitochondrial genome[END_REF].

Whilst the above appears to be also true for birds, i.e. the other homeothermic animals, it appeared that ectothermic vertebrates and invertebrates have a much larger mt size and size variation [START_REF] Rand | Endotherms, ectotherms, and mitochondrial genome-size variation[END_REF]. Ectothermic vertebrates may have genomes over 10kB larger than that of the largest homeotherm, and some invertebrates have a genome over twice the size of homeotherms (e.g. 35 Kb for Armadillidium vulgare or 42 Kb Placopecten magellanicus). Large insertions or duplications, or large scale-variation of tandemly repeated sequences are not uncommon. Recombination through these sequences may generate heteroplasmy, as was first observed within crickets [START_REF] Rand | Molecular population genetics of mtDNA size variation in crickets[END_REF] or scallops [START_REF] Gjetvaj | Repeated sequences and large-scale size variation of mitochondrial DNA: a common feature among scallops (Bivalvia: Pectinidae)[END_REF], and are associated with the presence of short recombining repeats in direct orientation. However, changes in the order of the genes are rare among animal species as such are unlikely to arise independently in separate evolutionary lineages [START_REF] Boore | Animal mitochondrial genome[END_REF].

Fungal mitochondria have considerable variation in genome size, ranging from 18 to 174 kb [START_REF] Clark-Walker | Evolution of mitochondrial genomes in fungi[END_REF]. The increase in size partly results from the presence of introns and noncoding sequences. This suggest that selection for genome size is not as strong as in animal mitochondria. In fungal mitochondria, recombination is often facilitated by the presence of short directly repeated sequences, and it is noticeable that larger genomes are more prone to rearrangements than small ones [START_REF] Clark-Walker | Evolution of mitochondrial genomes in fungi[END_REF]. In Saccharomyces cerevisiae recombination has been extensively studied not least because it is involved in the generation of "petite" mutants. Recombination may be so frequent than up to 10% of petite mutants may be generated per generation [START_REF] Dujon | Mitochondrial instabilities and rearrangements in yeasts and fungi[END_REF]. As a consequence, small differences in gene order may sometimes be observed within the same species.

Oomycete mitochondria and angiosperm chloroplasts have striking similarities. The size of the former varies from 30 to 80 Kb [START_REF] Mcnabb | Uniformity of mitochondrial DNA complexity in oomycetes and the evolution of the inverted repeat[END_REF] while the size of the later vary from 120 to 217 Kb [START_REF] Palmer | Conservation of chloroplast genome structure among vascular plants[END_REF]. In both cases, most of the genome size variation comes from variation in the length of an inverted repeat, that may make up nearly three quarters of the genome. Apart from this variation, the order of the genes is generally conserved. When present, the inverted repeat leads to the coexistence of two isomers generated by recombination [START_REF] Palmer | Chloroplast DNA exists in two orientations[END_REF][START_REF] Boyd | Evolutionary stability of DNA organization in Achlya[END_REF]. In both lineages, the inverted repeat always contain ribosomal genes, if nothing else.

Plant mitochondria have the largest genome and with it considerable genome size variation. It varies between 200 and 2500 Kb, of which a variable proportion is coding [START_REF] Schuster | The plant mitochondrial genome: physical structure, information content, RNA editing, and gene migration to the nucleus[END_REF]. A comparison of different angiosperm species reveals that most, but not all genome size variation is due to variation in the number of large duplications. Size variation probably often appears in evolution, since large variations were found among closely related species, and the whole range of variation was even found within the same genus, namely the Brassicae. This suggests an absence of selection on genome size, or at least that if such a selection exists, it is of little importance.

Angiosperm mitochondria have the highest rate of rearrangements, due to the presence of small repeated sequences with a very high rate of recombination. Recombination between these sequences is so frequent that a totally different order of the genes is often found between closely related species, within different individuals of the same species, or even within the same individual. This generation of diversity through recombination is well documented [START_REF] Schuster | The plant mitochondrial genome: physical structure, information content, RNA editing, and gene migration to the nucleus[END_REF].

Rate of point substitution

Animal mitochondria have an extremely high rate of substitution. In mammals, it is about ten times more than in the nucleus [START_REF] Wolfe | Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast and nuclear DNAs[END_REF]. The substitution rate of vertebrate ectotherms, is highly variable and dependant on the species. The substitution rate is much slower, five to ten times less than in homeotherms, in some marine turtles [START_REF] Avise | Mitochondrial DNA evolution at a Turtle's pace: evidence for low genetic variability and reduced microevolution rate in the Testudines[END_REF], and all studied amphibians and fishes, from sharks to salmonidae [START_REF] Martin | Body size, metabolic rate, generation time, and the molecular clock[END_REF]. By contrast, the substitution rate is much higher (up to ten times higher than in homeotherms) in snakes, alligators [START_REF] Kumazawa | The complete nucleotide sequence of a snake (Dinodon semicarinatus) mitochondrial genome with two identical control regions[END_REF] and the turtle Pelomedusa subrufa (D. Casane, pers. com.). The diversity of rates appears to covary with the thermal habits of the above species. While species of the first group do not control their temperature, species of the second group control their temperature by behaviour (warming themselves under the sun). The relationship between the rate of mitochondrial substitution and life history traits has already been studied [START_REF] Martin | Body size, metabolic rate, generation time, and the molecular clock[END_REF][START_REF] Rand | Thermal habit, metabolic rate and the evolution of mitochondrial DNA[END_REF], but discarded by the discovery of ectotherms with high mutation rate (e.g. [START_REF] Kumazawa | The complete nucleotide sequence of a snake (Dinodon semicarinatus) mitochondrial genome with two identical control regions[END_REF]. Further analysis including these new data and taking into account the behaviour of the animals may be of great interest.

Plant mitochondria have the lowest known substitution rate among living organisms: about ten times less than the nuclear rate, that is 100 times less than mammalian mitochondria! The fact that this low substitution rate is found even in the case of synonymous mutations [START_REF] Wolfe | Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast and nuclear DNAs[END_REF] indicates that the low substitution rate is not the result of selection on genomes with a low mutation number but rather that, above all, the mutation rate is low. This very low substitution rate, especially when interpreted as the consequences of a low point mutation rate, suggests that strong selection exists to lower the mutation load in this genome.

Yeast and oomycetes mitochondria and angiosperm chloroplasts do not differ very much in their rate of substitution, which is slightly higher than that of the nucleus. In oomycete mitochondria [START_REF] Boyd | Evolutionary stability of DNA organization in Achlya[END_REF] and in angiosperm chloroplasts [START_REF] Wolfe | Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast and nuclear DNAs[END_REF] the substitution rate is significantly lower in the inverted repeat than in the single copy region. As a consequence, genes located in the inverted repeats, namely those coding for ribosomal RNA, are the least variable.

Variation in substitution rates may result either from variation in the rate at which deleterious mutations accumulate, or from variation in the rate at which they are generated. Comparison between synonymous and non-synonymous mutations allows one to differentiate these two mechanisms. The existence of variation in the rate at which mutations (synonymous or not) accumulate is suggested in mammals by [START_REF] Lynch | Mutation accumulation in transfer RNAs: molecular evidence for Muller's ratchet in mitochondrial genomes[END_REF], and in aphid endosymbionts by [START_REF] Lambert | Deleterious mutations destabilize ribosomal RNA in endosymbiotic bacteria[END_REF]. These variations may come either from difference in the fidelity of replication, or from differences in the efficiency of repair.

Cytoplasmic and nuclear genomes are replicated by different specific polymerases. In vertebrate, plant and fungi, a polymerase of the gamma type is devoted to mitochondrial DNA synthesis. This polymerase is capable of proof reading activity and exhibits a high in vitro fidelity [START_REF] Lecrenier | Mitochondrial DNA polymerases from yeast to man: a new family of polymerases[END_REF]. Therefore the difference in in vivo mutation rate is thought to result from differences in the post-replication repair machinery, however, a DNA polymerase of the beta type that involves a DNA repair activity has been identified. The additional presence of a beta polymerase is suspected in trypanosomid Crithidia fasciculata and in yeast. The variation of the replication faithfulness might be more related to the absence or deficiency of a repair system, such as the mutS-like gene [START_REF] Pont-Kingdon | A coral mitochondria mutS gene[END_REF] than to the efficiency of different DNA polymerases. The involvement of a post-replication repair mechanism is also consistent with the observed link between recombination and the mutation rate.

The relation between genome parameters and Muller's ratchet

It is clear from table 1 and from the above descriptions, than the different types of genome structure and evolution don't follow the phylogenetic relationship between the lineages. However, the different parameters are not independent: the smaller the genome, the smaller the genome size variation and the frequency of rearrangements, the higher the rate of substitution and the number of copies per cell. Animal mitochondria are at one extreme with a small genome with high substitution rate, and plant mitochondria at the other extreme, with a large genome and a low substitution rate. Plant chloroplasts, along with oomycete and fungal mitochondria, lie in between.

Let us first consider the parameters involved in Muller's formulation of the problem of the accumulation of deleterious mutations [START_REF] Muller | The relation of recombination to mutational advance[END_REF]. The rate at which the deleterious mutations accumulate is dependent upon the number of genomes in the population, or population size (N), the number of mutable loci per genome (l) and the rate of mutation at each locus (). Only slightly deleterious mutations are taken into account, because the ratchet will not advance if the deleterious effect of each mutation (s) is sufficiently large for the corresponding genotype to be eliminated immediately by selection. The ensemble of these parameters determine No, the number of genomes which have the lowest number of mutations, which at equilibrium will be given by : No = N exp -(l/s) The smaller No the greater the probability that all mutation free genomes are expected to be lost by drift. Hence, as No increases so Muller's ratchet is expected to slow. That is, the ratchet can be slowed by either increasing N, the population size, or decreasing l, the number of loci, or decreasing , the mutation rate, or increasing s, the deleterious effect of a single mutation. This prediction holds even without considering the formal expression of Muller's equation. Furtheremore, other formulations of mutational meltdown, such as those involvinf background selection, also link to the same prediction than for asexual organisms to persist, either the population size must be high, or the mutation rate low, or the genome must have few loci. For instance, asexual protists probably avoid mutational meltdown because their nuclear genome is relatively small and their population size is very high [START_REF] Bell | Sex and death in protozoa[END_REF].

As a consequence of the necessity to avoid Muller's ratchet, a negative relation may be inferred between population size and genome size, and a positive relation may be inferred between genome size and mutation rate. We can relate these parameters to the properties of cytoplasmic genomes by making the assumption that the population size is related to the number of copies per cell, that the mutation rate is related to the synonymous substitution rate, and that the number of loci per genome is related to the genome size. Then a negative relation is expected between the number of copies per cell and genome size, and a positive relation may be inferred between genome size and substitution rate. This prediction, deduced from the functioning of Muller's ratchet, appears to be consistent with the qualitative pattern between the above parameters as described in table 1.

For animal mitochondria, and especially for mammals, the large number of copies per cell and the small genome size is sufficient to slow the ratchet, as first suggested by Bell (1888). For plant mitochondria, by contrast, the number of copies of the genome is at least one order of magnitude smaller than in animals and there is no evidence of selection favouring the smallest genome size, which is one to two orders of magnitude greater than in animals. With such a genome and such a number of copies, it is to be expected that Muller's ratchet should click very rapidly. It would seem, in this case, that the low rate of mutation of this genome (100 times smaller than in homeotherms' mitochondria) is central to the prevention of the accumulation of deleterious mutations. The fact that chloroplasts and fungal mitochondria lie in between for both genome size and mutation rate may indicate that they follow the rule, although they have a different strategy. It is also possible that the existence of an inverted repeat, by allowing recombination and reducing the rate of point substitution of the genes within the repeat, is part of this strategy. Locating genes in this inverted repeat may have been selected because it protects these genes against mutational meltdown. The fact than the genes present in the inverted repeats (rRNA genes) are essential genes known to undergo strong selective pressure, are in favour of this hypothesis. It is also possible to interpret as the result of such a selection the fact that vertebrate ectotherms with a very high substitution rate (such as snakes) have a duplicated control region with concerted evolution [START_REF] Rand | Thermal habit, metabolic rate and the evolution of mitochondrial DNA[END_REF].

We therefore propose that the rules underpinning the molecular evolution of cytoplasmic genomes (i.e. genome size, substitution rate and number of copies per cell), is in agreement with the necessity to prevent the accumulation of deleterious mutations within a genome with uniparental inheritance. To test this hypothesis consider one organism with biparental inheritance with clear evidence of frequent recombination between parental genomes. This is the case for yeast, where we may not consider the mitochondrial genome as asexual. Yeasts do not follow the above rule in sense that they have a genome size and substitution rate comparable to that of oomycetes with uniparental inheritance, but they have an especially low number of copies per cell: only a few dozen. This supports the hypothesis that selective pressures acting on the above parameters is dependent of the necessity to avoid mutational meltdown.

This does not mean that the relation between these parameters was directly selected for the avoidance of Muller's ratchet. Selection for metabolic efficiency may also produce comparable relationship. Consider an organism where, for one reason or another (maybe a high energetic demand), there is a selection for rapid replication and transcription of the cytoplasmic genome. This may be achieved by both decreasing genome size (for faster transcription), increasing copy number (for multiple simultaneous replication or transcription events) and increasing the rapidity of the DNA polymerase. As a consequence of rapid replication, the efficiency of repair is decreased and the mutation rate goes up. This alone may explain a non-independence between genome size, mutation rate and copy number. That animal mitochondria properties seem to be related with homeothermy or thermal habits of ectotherms would support such a metabolic explanation. Further, if allocation of DNA to a given organelle is a fixed parameter then the small size and copy number would be expected to covary.

Both mechanistic and mutational decay models predict that only a limited number of parameter combinations could be successful and that alteration of one parameter will select for a co-alteration of others. If Muller's ratchet is a long term process only those genomes with an appropriate parameterization will survive in the long term. Alternatively, mutational decay, constraints on DNA replication and adaptation to metabolic demands may be short term forces. We remain agnostic as to whether Muller's ratchet is a strong selective force promoting adaptation or a long term filter.

Other ways to avoid mutational meltdown

Increasing the efficiency of selection between gametes One way to increase the efficiency of selection against deleterious mutations is to allow their grouping in the same individual cell or organelle. In an asexual lineage, this can be achieved by the mode of transmission during cell division or by the germline bottleneck that occurs in mammals. This was first examined by [START_REF] Takahata | Evolutionary dynamics of extranuclear genes[END_REF], and further modelled by [START_REF] Bergstrom | Germline bottlenecks and the evolutionary maintenance of mitochondrial genomes[END_REF], who show that "while a bottleneck indeed increases the rate of genetic degradation within a particular (cell) lineage, it also serves to strengthen selection among lineages, and hence has a net effect of resisting genetic decay". A different model simulating the functioning of angiosperm mitochondrial genome [START_REF] How | Dynamics of plant mitochondrial genome: Model of a three-level selection process[END_REF] showed that selection at the mitochondrial level may prevent the fixation of deleterious mutations (in this case, deletions). Most of work on this topic is however theoretical, since experimental data on the mode of transmission of the different cytoplasmic genome are rare [START_REF] Boore | Transmission of mitochondrial DNA-playing favorites?[END_REF][START_REF] Yaffe | The machinery of mitochondrial inheritance behavior[END_REF]. Another way to increase the efficiency of selection is to decrease the number of surviving cells. A comparative analysis made among animals by [START_REF] Krakauer | Mitochondria and germ-line death[END_REF] led them to conclude that the high level of germ-cell death in female germ lines may be considered as a "developmental solution" to avoid Muller's ratchet in mitochondria.

Paternal transmission

The more obvious way to escape the ratchet is, of course, to have some level of paternal inheritance and recombination. Indeed, numerous models show that a low rate of recombination may be adequate to prevent mutational meltdown. Uniparental inheritance is not equally displayed through eukaryotic organisms. Some degree of biparental inheritance may be relatively frequent among sexual isogamous protists, where two haploid cells of equal size fuse to give a diploid. In angiosperms, maternal transmission is the rule, but there are some cases of rare paternal transmission, especially for chloroplasts [START_REF] Reboud | Organelle inheritance in plants[END_REF]. In higher animals transmission is usually strictly maternal although rare paternal transmission may be found [START_REF] Wallis | Do animal mitochondrial genome recombine?[END_REF]. In gymnosperms cytoplasmic genes are usually uniparentally inherited, but paternal transmission of chloroplasts seems to be more frequent that in angiosperms [START_REF] Chesnoy | La reproduction sexuée des gymnospermes[END_REF].

These paternal transmission events could play an important role in stopping the ratchet. It is however unclear just what level is necessary to have any effect. Futhermore, the parameter that has an importance here is the rate of paternal transmission that results in recombination between maternal and paternal mitochondria. In most lineages, this rate is almost impossible to estimate. However, we can predict that the higher the rate of paternal transmission/recombination, the lower the evolutionary constraints induced by Muller's ratchet.

Move of cytoplasmic genes to the nucleus

The ancestral prokaryotic mitochondrial genome has undergone massive gene transfer from the primary symbiont to the nucleus. If the species is sexual, this results in reducing the number of genes that are not submitted to sexual recombination.

Post-replicative mechanisms

In mammals, the ratchet is not totally stopped. As demonstrated by [START_REF] Lynch | Mutation accumulation in transfer RNAs: molecular evidence for Muller's ratchet in mitochondrial genomes[END_REF], mitochondrial tRNA of mammals accumulated deleterious mutations that reduce their stability compared to their nuclear counterpart. This author, however, suggests that the functionality of the mitochondrial tRNA may not be necessarily affected, because compensatory mutations could prevent a net loss of fitness. Arguments in favour of this hypothesis came from the study of [START_REF] Lambert | Deleterious mutations destabilize ribosomal RNA in endosymbiotic bacteria[END_REF]. By comparing aphid endosymbionts and their closest free leaving bacteria, Moran found that the 16S RNA accumulate more mutations in five independently evolving symbiont lineages than in their free leaving relatives. Moran also found that some endosymbionts contain a very high concentration of a chaperone protein that may help stabilise mutationally decayed proteins.

RNA editing has been found in almost all lineages, although it is particulary common among angiosperms [START_REF] Schuster | The plant mitochondrial genome: physical structure, information content, RNA editing, and gene migration to the nucleus[END_REF] and trypanosomes [START_REF] Arts | Mechanism and evolution of RNA editing in kinetoplastida[END_REF]. By modifying the messenger RNA to restore the ancestral sequence, it may also be viewed as a way to maintain the functionality of a genome, once a deleterious mutation has been fixed at the DNA level. The fact that editing results in reconstituting the original structure of a protein was first described in land plants [START_REF] Gualberto | RNA editing in wheat mitochondria results in the conservation of protein sequences[END_REF]. In Kinetoplastida [START_REF] Arts | Mechanism and evolution of RNA editing in kinetoplastida[END_REF], it is even suggested that extensively edited RNA may be converted in DNA through retrotransposition. In that case, RNA editing not only allows the proteins to stay functional despite the ratchet, but would even allow the ratchet to go in reverse!

Conclusion

Most current studies on cytoplasmic genome evolution are based on sequence data. This kind of data is powerful for estimating mutation and substitution rates, as well as estimating the number of loci under selection. However, cytological data on the number of copies per cell, their partitioning among organelles, and the way they segregate during cell division are necessary to understand the link between the above parameters. Furthermore, life history traits such as the rate of paternal transmission and recombination, metabolic demand, thermal habit, parasitic vs. free-living, and others may be of importance. It would also be useful to study cytoplasmic DNA of a phylogenetically broad representation.

Despite these limitations, and the imperfect correspondence between measurable data and the theoretical parametrization, it appears that molecular evolution of cytoplasmic genomes is, at the very least, enlighted by current theories of the evolution of sex. Conversely a comparative study of the structure and evolution of cytoplasmic genomes allows a better understanding of the evolutionary constraints leading to the prevalence of sex.

Acknowledgements

We wish to thank Professors John Maynard Smith and Bengt Bengtsson for discussion. We also thank Professor Graham Bell for suggesting the problem.


Sylvain Devie

Pierre-Philippe Robet

Yannick Aoustin

Maxime Gautier

Estimation of the External Force for a Co-manipulation Task Using the Drive Chain Robot

Keywords: Control, Identification, Robot, Co-manipulation, Sensor-less

The aim of this paper is to show that the observation of the external effort and the sensor-less control of a system is limited by the mechanical system. First, the model of a one-joint robot with a prismatic joint is presented. Based on this model, two different procedures were performed in order to identify the mechanical parameters of the system and observe the external effort applied on it. Experiments have proven that the accuracy of the force observer, based on the DC motor current, is limited by the mechanics of the robot. The sensor-less control will be limited by the accuracy in estimation of the mechanical parameters and by the maximum static friction force, that is the minimum force which can be observed in this case. The consequence of this limitation is that industrial robots without specific design are not well adapted to perform sensor-less precision tasks. Finally, an efficient control law is presented for high effort applications.

I. INTRODUCTION

F ORCE control is an important part of the mechatronic control [START_REF] Ohnishi | Motion control for advanced mechatronics[END_REF]. This type of control system is currently made by using force sensors in order to measure the external effort. However, in particular cases, this solution cannot be applied and a sensor-less solution by using motor currents has to be chosen. For example, in [START_REF] Li | Current-sensor-based feed cutting force intelligent estimation and tool wear condition monitoring[END_REF], the electrical current measurement is used to estimate the feed cutting force for a manufacturing operation. In [START_REF] Qi | A foot force sensing approach for a legged walking robot using the motor current[END_REF], the current of the DC motors of a legged walking robot was used to estimate the ground reaction force based on the dynamic equation in the joint space for each joint. The inner friction were identified by using a neural network to improve the identification. Recently, Kambara et al. have proposed, in [START_REF] Kambara | A design method of observers for bilateral control using dc brushed motor[END_REF], an innovative disturbance force observer based on the model of their minor current loop. Then this observer is used to carry out acceleration-based control for tele-operation tasks without external force sensor and have proven its efficiency by comparing it to more classical observation-based controls.

In the case of force control, and more specifically in case of co-manipulation tasks, the external perturbations are mostly due to the mechanical part of the system. These perturbations have a direct effect on the observation of the effort and have to be identified or mechanically compensated. For example, the Haption master arm [START_REF] Gosselin | Design of a High Fidelity Haptic Device for Telesurgery[END_REF] has been designed in order to be mechanically transparent. This means that all the mechanical losses are compensated and the effort given by the actuator is exactly equal to the effort applied on the environment by Sylvain Devie, Pierre-Philippe Robet, Yannick Aoustin and Maxime Gautier are with the LS2N (Digital Sciences Laboratory of Nantes) (Sylvain.Devie, Pierre-philippe.Robet, Yannick.Aoustin, Maxime.Gautier)@univ-nantes.fr the end-effector. With this design, sensor-less force control can be performed for surgical tasks. In the case of industrial robots, the power from the actuator is shared between the useful work, dynamic effects like inertia or the Coriolis effect, and perturbation effects like friction. In this spirit, in [START_REF] Canudas De Wit | A new model for control of systems with friction[END_REF] the authors proposed a friction observer in order to control industrial robots by using an innovative estimation of the friction perturbation.

In the case of the co-manipulation control, the perturbation are be divided into two categories: the effort from the operator and the friction perturbation. The effort from the operator is one of the inputs of the control law. In the classical application, we desire the velocity of the robot to be proportional to this effort. The friction perturbation consists of mechanical losses in the transmission which consumes energy between the actuator and the end-effector. The identification of these perturbation presents some difficulties. For example, the breakaway, between the state of rest and motion is difficult to define, and several empirical models have been proposed ( [START_REF] Armstrong-Helouvry | Stick-slip arising from Stribeck friction[END_REF], [START_REF] Hess | Friction at a Lubricated Line Contact Operating at Oscillating Sliding Velocities[END_REF], [START_REF] Bo | The friction-speed relation and its influence on the critical velocity of stick-slip motion[END_REF]). Another difficulty is the non-repeatability of the friction identification, which makes the friction compensation difficult ( [START_REF] Gilbart | Adaptive compensation for an optical tracking telescope[END_REF], [START_REF] Walrath | Adaptive bearing friction compensation based on recent knowledge of dynamic friction[END_REF]). This difficulty can be minimised with online correction of the friction estimation. For example, in [START_REF] Roveda | Iterative Learning Procedure with Reinforcement for High-Accuracy Force Tracking in Robotized Tasks[END_REF] the friction estimation is corrected on-line in order to carry out impedance control.

More globally, the mechanical parameters of an industrial robot are estimated with a general identification method, based on the inverse dynamic model (IDM) and the least squares (LS) method. This method has been successfully applied to identify the inertial and friction parameters of many prototypes and industrial robots ( [START_REF] Khosla | Parameter identification of robot dynamics[END_REF], [START_REF] Gautier | Numerical calculation of the base inertial parameters of robots[END_REF]). The identification consists of moving a robot without a load (or external force) or with a constant given payload [START_REF] Khalil | Identification of the payload inertial parameters of industrial manipulators[END_REF], by following a specific trajectory [START_REF] Gautier | Exciting trajectories for the identification of base inertial parameters of robots[END_REF]. More complicated models can also be solved in order to consider the flexibility of the system [START_REF] Gautier | Dynamic Identification of flexible joint manipulators with an efficient closed loop output error method based on motor torque output data[END_REF].

The aim of this paper is to illustrate the limit for the force observation for an industrial robot in the case of a comanipulation task. This task is similar to the tele-operation task presented in [START_REF] Kambara | A design method of observers for bilateral control using dc brushed motor[END_REF]. The difference is that in the case of comanipulation a single robot is considered and the force applied at its end-effector by the operator is the reference effort that the robot has to follow, while two robots are considered in the case of tele-operation: the force applied on one robot is the reference for the force applied by the other one. For the experimental part, a one-link robot with a prismatic joint is presented and is assumed to be rigid. This robot testbed has already been used to test several identification processes over the last few years ( [START_REF] Hamon | Dynamic Modeling and Identification of Joint Drive with Load-Dependent Friction Model[END_REF], [START_REF] Gautier | Dynamic Identification of flexible joint manipulators with an efficient closed loop output error method based on motor torque output data[END_REF]) and is a good tool to perform this kind of study, as recently shown in [START_REF] Yashiro | Fast Stiffness Estimation using Acceleration-based Impedance Control and its Application to Bilateral Control[END_REF]. There, an acceleration-based impedance control is used in a tele-operation scenario for fast environmental stiffness estimations with a time delay. Or in [START_REF] Seki | Improvement of Force Transmission Performance Considering Nonlinear Friction in Bilateral Control Systems[END_REF] where the bilateral control is improved by considering a non-linear model for the perturbation observer.

From the work of Onishi [START_REF] Ohnishi | Motion control for advanced mechatronics[END_REF] a force observer is first implemented for our robot. An identification of its parameters was a priori made. The aim of this force observer will be to check if it is possible to perform a precise task with a low force observer. After that, an impedance based control law based on Hogan's ( [START_REF] Hogan | Impedance Control: An Approach to Manipulation[END_REF], [START_REF]Stable execution of contact tasks using impedance control[END_REF]) work is implemented. This control law allows us to successfully perform co-manipulation tasks with large external forces. In [START_REF] Devie | Impedance control using a cascaded loop force control[END_REF], an equivalence was shown between this control law and the one already used in the harware used in this paper ( [START_REF] Devie | Accurate force control and co-manipulation control using hybrid external command[END_REF], [START_REF] Devie | Optimized force and co-manipulation control using stiffness of force sensor with unknow environnement[END_REF]).

This paper is outlined as follows. Section 2 focuses on the model of the robot and presents its Newton-Euler equation in order to define a linear form of the Inverse Dynamic Model. Section 3 uses this model to perform the identification of the mechanical parameters of the robot and the observation of the effort of the operator while performing a co-manipulation task. Section 4 proposes an efficient sensor-less control law for high forces and Section 5 offers our conclusion and perspectives.

II. MODELING OF A RIGID INDUSTRIAL ROBOT

A. Experimental set-up

The EMPS is a high-precision linear Electro-Mechanical Positioning System, presented in Fig. 1. It is a standard configuration of a drive system for the prismatic joint of robots. It is actuated by a Maxon DC motor which drives a carriage by a ball screw. The carriage moves a force sensor and a gripper in translation. The motor rotor and the ball screw are connected by a flexible coupling. Two incremental encoders are presented on the robot. The first one measures the motor position q m (rad) while the second one measures the position of the ball screw q 1 (rad).

For the applications considered in this paper, the robot is assumed to be rigid in the frequency range of the system harmonics. This means that the encoder measurements are equal and q m = q 1 . In the following, we will consider only the displacement of the carriage q = q m /r (m), where r is the pitch of the ball-screw (rad/m).

An inner current loop is physically integrated at the input of the DC motor and allows it to control the motor without using the current sensor. This helps the hardware of the system becomes simpler. This current regulation cannot be easily modified by the operator and has a frequency ω I 20 times greater than the velocity loop. In the following, the current loop will be assimilated into a simple gain.

In the following, all the variables on the load side will be given in SI units. The force applied by the motor on the carriage is linearly proportional to the motor torque with a factor of the reduction ratio of the ball screw r. This torque is linearly proportional to the current sent to the motor I m (A), with a factor of the torque constant k t (N/A). A current amplifier is present on the input of the current loop v I (V ), with I m = G i v I .

In the following, we will consider:

τ m = G i r k t v I = G τ v I (1)
The parameter G τ (N/V ) will be considered to be known and independent of the experimental conditions.

B. Inverse dynamic model of the robot

The mechanical system, presented in Fig. 2 has one prismatic degree of freedom q (m). It is affected by the actions of the motor τ m (N) and of the environment τ e (N). In the considered case, the body of the robot is not affected by gravitational force.

The Inverse Dynamic Model (IDM) calculates the motor force as a function of the joint position and its derivatives. Newton-Euler equations give the following IDM [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]:

τ m + τ e = M q + F v q + F c sign( q) (2)
Here:

• q (m/s) and q (m/s 2 ) are respectively the carriage velocity and acceleration, • τ m (N) is the drive force from the motor, • M (kg) is the total load side equivalent mass, including the inertia of all rotating elements in the drive chain and the mass of all the translating elements, • τ e (N) is the external force, applied by the environment,

• F v (N/m/s) is the viscous damping coefficient, • F c (N) is the Coulomb friction force,
• sign(u) denotes the sign function. The friction considered here, as shown in Fig. 3, is a classical model that takes into account only the kinetic (Coulomb) friction and viscous damping:

F f c ( q) = F v q + F c sign( q) (3)
However, this model does not take into account the static friction: in the static case, a break-away force F s (N) is needed to overcome the static friction and move the robot, with F s > F c . This break-away is very difficult to model efficiently, the most efficient empirical model to approximate it is the Stribeck model, taking into account the non-linearity due to adherence at the null velocity:

F f s = F v q + F c sign( q) + (F s -F c )e -| q/ qs | δs sign( q) ( 4
)
where qs (m/s) is the Stribeck velocity constant and δ s is a coefficient between 1 and 2. The shape of this model is also presented in Fig. 3. The classical model restricts us to experiments in areas far from the non-linear area, away from null velocities.

Another model problem would be the unsymmetrical behaviour of friction [START_REF] Hamon | Dynamic Modeling and Identification of Joint Drive with Load-Dependent Friction Model[END_REF]. In order to deal with this, a two quadrant model will be considered for the Coulomb and viscous friction, depending on the sign of the speed. This implies the following:

• F + c , F + v and q+ will respectively be used for F c , F v and q when q > 0 • F - c , F - v and q-will respectively be used for F c , F v and q when q < 0 Encoder on the motor side Ball screw Force sensor

Carriage Flexible coupling

Encoder on the load side Motor 2) can be re-written as function of the model parameters, gathered in vector X.

τ m + τ e = WX (5)
Here :

• X = M F + v F - v F + c F - c is a 5 × 1 vector • W = q q+ q-sign( q+ ) sign( q-) is a 1 × 5 ma-
trix. This equation will be used in the following in order to identify the parameters in X.

C. Control design

In a previous study [START_REF] Devie | Accurate force control and co-manipulation control using hybrid external command[END_REF], a simple control law was studied for the force control using an integral proportional (IP) cor-rection for the velocity loop and a proportional (P) correction for the force loop. Then this control law was adapted in order to perform the co-manipulation task.

In a co-manipulation case, the environmental impedance depends on the impedance of the environment and of the operator's hand. This environmental impedance is supposed to be unknown. In this case, a classical control system is illustrated in Fig. 4. Here, an inner velocity loop is used to control the performance of the system while an outer force loop is used to control its transparency. Increasing the gain k e1 , we increase the transparency. However, a too high gain can lead to instability. The aim is to get a simple comanipulation controller with τ e = 0 when the robot is not moving ( q1 = q2 = 0). It leads us to choose the force reference τ re f = 0 which is the offset of the force the operator has to apply on the system.

Thanks to this reference, the system allows a linear relationship between the external force τ e and the velocity reference v q. If the external perturbation has a low frequency, the relation τ e = k c q1 can be used in order to calculate the correction gains of the closed loops.

In order to calculate the correction and the performance of the system, we consider that the frequency range (< 20 rad/s) of the exogenous disturbance τ e is small compared to the bandwidth of the velocity loop (100 rad/s), in order to provide a linear relation between the velocity and the external force: τ e = k c q1 with k c = 1/k e1 . However, in the case of low external forces, the external force can be neglected for the calculation of the velocity loop correction. According to Fig. 4 and considering this last assumption, for co-manipulation applications the open loop transfer function of the velocity loop is given by the equation:

T vo (s) = q1 (s) v q(s) -q1 (s) = k v1 G τ t v s 1 M 1 s + F v1 + k v1 G τ (6)
Imposing the phase margin φ v at a frequency ω v leads to T vo ( jω v ) = 1e j(-π+φ v ) and gives the values of k v1 and t v [?].

k v1 = M 1 ω v tan(φ v ) -F v1 G τ , t v = k v1 G τ M 1 ω 2 v cos(φ v ) (7)
The closed loop transfer function for the velocity loop is: with

T vc (s) = q1 (s) v q(s) = 1 1 + (t v + F v1 K v )s + M 1 K v s 2 (8)
K v = k v1 G τ /t v .
Let us write this transfer function in the canonical form:

T vc (s) = 1 1 + 2z v ω 0v s + s 2 ω 2 0v ( 9
)
with

ω 0v = K v /M 1 and z v = (t v + F v1 K v ) ω 0v 2 .
The outer open force loop has the following equation:

T eo (s) = τ e τ re f + τ e = k e1 k c 1 + 2z v ω 0v s + s 2 ω 2 0v ( 10
)
According to the definition of k c , which is equivalent to a viscous friction coefficient, the equation can be written in the following form:

T eo (s) = 1 1 + 2z v ω 0v s + s 2 ω 2 0v (11)
Here, the coefficient k c is chosen to provide good transparency.

III. REACTION FORCE OBSERVER FOR SMALL EXTERNAL

FORCES

A rigorous observation of the reaction force needs a perfect estimation of the electro-mechanical parameters of the system. Hence, the observation of the reaction force needs two steps: the identification of the system parameters, and the estimation of the reaction force from the parameters. In the considered case, the electrical parameters are perfectly known and used to control the inner velocity loop, however the mechanical parameters have to be identified.

A. Identification of the parameters

Considering an acquisition free of reaction force, with a position control of the robot allowing it to follow a specific trajectory, equation ( 5) becomes:

τ m1 = W 1 X 1 (12)
Where τ m1 is a n × 1 matrix of motor forces, X 1 is a p × 1 matrix of mechanical parameters and W 1 is n × p matrix of velocity and acceleration signals calculated numerically from the measured positions. Here p = 5 is the number of independent parameters and n is the number of samples. In this section, the subscript 1 is used to describe data from the first set of experiments related to identification.

According to this model, for a robot following a dynamic trajectory [START_REF] Gautier | Exciting trajectories for the identification of base inertial parameters of robots[END_REF], the matrix X 1 can be estimated with the least squares method, for example by using the Matlab function:

X 1 = W 1 \ τ m1
Standard deviations are estimated using classical and simple results from statistics, considering the matrix W 1 to be a deterministic one, and ρ to be a zero mean additive independent noise, with standard deviation σ ρρ such that :

C ρρ = E(ρρ ) = σ 2 ρ I N ( 13
)
where E is the expectation operator and I n is the (n × n) identity matrix. An unbiased estimation of σ ρ is used and given by the expression :

σ 2 ρ = τ m1 -W 1 X 1 2 n -p (14)
The variance-covariance matrix of the estimation error and standard deviations can be calculated by :

C X 1 X 1 = E[(X 1 -X 1 )(X 1 -X 1 ) ] = σ ρ 2 (W 1 W 1 ) -1 (15)
The relative standard deviation is given by the expression:

%σ X 1 (i) = 100| σ X 1 (i) X 1 (i) | (16)
where

σ X 1 2 (i) = C X 1 X 1 (i, i) is the i th diagonal coefficient of C X 1 X 1
. Because of the perturbations due to the noise and the modeling errors, the actual force τ m1 differs from the estimated one τ m1 by an error e 1 , such that:

τ m1 = τ m1 + e 1 = W 1 X 1 + e 1 (17)
Calculating the least squares solution [START_REF] Gautier | Dynamic Identification of flexible joint manipulators with an efficient closed loop output error method based on motor torque output data[END_REF] from the data in W 1 and τ m1 can lead to a bias if W 1 is correlated to e 1 . Hence, it is essential to filter the data in τ m1 and W 1 before computing W 1 \ τ m1 .

Velocities are directly estimated by a backward derivative of the position, and the acceleration is estimated by a central derivative of the filtered velocity. All the variables were corrected with the medfilt1 and decimate functions of Matlab R Ṫhe measured and filtered values of the effort from the motor τ m1 , the velocity q and the acceleration q are presented in Fig. 5.

The velocity reference trajectory is a mixture of trapezoidal segments and constant velocity segments. Specific shapes of the trajectory highlight specific dynamic parameters. A linear variation of the velocity highlights the inertial parameters while a constant velocity highlights the viscous friction parameters. The reference trajectory is optimized in order to have a linear variation of velocity half the time, and a constant velocity in the other half of the time, in order to excite all the parameters in order to identify each of them with a good accuracy. Ten experiments were performed in order to identify the parameters of the robots. Fig. 6 shows the evolution of the motor effort τ m1 calculated thanks to the motor current measurement and its estimated values without taking the error into account W 1 X 1 .

The identified values for the experiment shown in Fig. 6 are given in Table I. It is important to remark that for the experiment presented in Fig. 6, the estimated Coulomb frictions is F + c1 = 17 N. Another static experiment allowed us to estimate the static friction as F + s = 17 N. This means that the classical model used for the observation is well adapted for a ball screw system. That also means that, in the case of sensor-less control, the operator has to apply a break-away force τ e = 17 N on a static robot in order to allow the controller to detect its action. = m1 W 1 X Fig. 6: Force from the motor τ m1 in blue and W 1 X 1 in red.

The root mean square (RMS) error between the motor force τ m1 and W 1 X 1 is RMS 1 = mean(e 1 2 ) = 3.6 N. The errors from the other experiments varies between 3.6 N to 5.3 N. This identification will be used for the observation in the next step. This means that the observer cannot be more accurate than 3.6 N in the best case.

In the following, X 1 represents the mean value of the mechanical parameters of X 1 and which were calculated with the ten experiments.

B. Observation of the reaction force

Once the model is known, the goal is to observe the reaction of the external forces. The co-manipulation control law presented in the previous section is considered. The aim of this operation is to control a robot with a transparent behaviour. This means that the operator should hold a gripper at the endeffector of the robot and move it without feeling the dynamics of the robot. This can be achieved thanks to high values of the gains of the outer force loop. However, we want to estimate the force applied by the operator on the robot. This interaction effort can be improved by reducing the transparency of the control law. The first step to perform a sensor-less task: if the estimation of the force τ e is not efficient enough, the control of the force τ m is impossible with any control law.

The control law presented in section II-C is used with the FN3280 sensor from Hoskin. The data collected thanks to this sensor is used to control the robot and will be compared to the estimated effort τ e1 from the observer.

Here, the aim is to observe the external force τ e1 . If the control system is transparent, the operator does not feel the force from the robot, and its action is compensated by the control system. In this case, the observation of the motor torque does not allow us to estimate this force. To avoid this issue, the gain k e1 is chosen to have a poor transparency of the system. The effort measured by the force sensor will be considered equal to the real effort from the operator, as the sensing error is insignificant compared to the identification error.

This effort is estimated with the help of ( 5), as follows:

τ e2 = τ m2 -W 2 X 2 (18)
In this section, the subscript 2 is used to describe data from the second set of experiments related to the observer. The vector X 1 from the mechanical parameters identified in the previous experiments will be considered instead of the vector X 2 . Fig. 7: Co-manipulation experiments, the operator apply a specific force on the force sensor of the robot presented into Fig. 1.

Ten experiments were performed in a co-manipulation scenario, they are illustrated by Fig. 7. During these experiments, a co-manipulation task was carried out, considering a non-null reference force in order to ensure an effort at any time. For each of these experiments, the force applied by the operator on the system was estimated by:

τ e2 = W 2 X 1 -τ m2 (19)
With τ m2 calculated with the help of equation [START_REF] Ohnishi | Motion control for advanced mechatronics[END_REF]. Again, because of the perturbations due to the noise, the modeling errors, and the identification errors, the actual force τ e2 applied by the operator differs from τ e2 by an error e 2 depending on e 1 , such that:

τ e2 = τ e2 + e 2 (20)
Two main cases of the calculation of τ m2 can be considered: the on-line calculation and the off-line calculation. The on-line calculation is done during the experiment, with a simple firstorder low pass filter on the measured velocity and current in order to obtain a real time estimation of the force, which can be used to perform a sensor-less co-manipulation operation. The off-line calculation is done after the experiment and allows better filtering of the data.

Considering the off-line calculation, the force applied by the operator is calculated by equation [START_REF] Yashiro | Fast Stiffness Estimation using Acceleration-based Impedance Control and its Application to Bilateral Control[END_REF]. The force measured by the sensor τ e2 and the estimated force τ e2 are presented on Fig. 8.

According to this figure, for a small force from the operator, the noise from the observer is bigger than the actual effort.

For the ten experiments, the RMS error of the observation RMS 2 = mean(e 2 2 ) varies from 3.9 N to 5.4 N. This result proves that the force observation with a simple 2 quadrant mechanical model is not well adapted for precision tasks.

However, the proposed model uses a cross validation of the observer -the observed effort is calculated thanks to the data identified with another experiment. This means that the error e 2 from the second experiment accumulates upon the error e 1 from the first experiment. This problem is cumulated with In our case, the used workspace was sufficiently small and the time between the two steps was sufficiently small to consider the estimation to be correct. However, in the case of industrial robots, the parameters identified in the first step could differ from the parameters which should be used in the second step.

Another approach to limit these errors is to carry out the two steps on the same set of data off-line.

C. Direct validation observation

Let us now consider only the co-manipulation experiment from the previous section. The efforts applied by the operator were specifically chosen in order to present particular behaviours like ramps and floors allowing the identification of the system. Here, the aim is not to propose a strategy for sensor-less control, but to use an identification and observation protocol based on force measurement to check the efficiency of the previous calculations.

Let us now consider [START_REF] Gosselin | Design of a High Fidelity Haptic Device for Telesurgery[END_REF]:

Y = τ m + τ e = WX (21)
Here, τ m1 is calculated thanks to the current measurement and equation ( 1), τ e1 is given by the force sensor and W is calculated off-line thanks to the motor encoder and time derivative computations.

According to the previous section, the dynamic parameters in X can be estimated by X with the help of a least squares method in order to approximate Y, with:

Y = W X + e i (22)
where e i is the vector of the identification error for the experiment, due to measurement noise and modelling errors.

The effort from the operator can be estimated similar to the previous section, by:

τ e = W X -τ m ( 23
)
The evolution of this effort is presented in Fig. 9, for the same experiment as in Fig. 8.

The identified values for the experiment shown in Fig. 6 are given in Table II The RMS error between the observed force and the actual effort from the operator is 2.8 N. For the ten experiments, this error varies from 2.7 N to 3.4 N, which is of the same order of magnitude as the value of the forces applied by the operator, but better than the observations of the previous section.

These results are mostly due to the fact that most of the power from the motor is used to compensate for the dynamics of the robot, and not to assist the operator. This means, in [START_REF] Devie | Impedance control using a cascaded loop force control[END_REF], τ m is approximately equal to W X. This assumption is illustrated in Fig. 10.

According to this figure, the motor forces and the forces due to the mechanical losses have a maximum value of 30 N, while the maximum value of the effort from the operator is equal to 5 N, which is the range of the sensing device. In this particular case, the effort applied by the operator has the same order of magnitude as the precision of the observation.

Several techniques can be proposed to reduce this problem. The most common technique is to design a specific robot with a mechanical compensation of the mechanical losses. In this case, the force needed to move the robot W X would be far smaller, and the force from the motor would be equal to the force from the operator. Another way to reduce the problem would be to perform large force operations instead of precision operations in order to have a force τ e higher than the precision of the observation. In this configuration, we can assume a good reconstruction of the force τ e . In any case, it appears that a sensor-less force control is not possible in the case of small forces for co-manipulation applications with industrial robots. However, the next section will present an efficient control law allowing for these kinds of applications for large external forces.

IV. SENSOR-LESS CO-MANIPULATION CONTROL FOR LARGE EXTERNAL FORCES

The aim of this section is to present an efficient control law in order to carry out a sensor-less co-manipulation task with an industrial robot, for large external forces based on on-line estimation of the operator force. Unfortunately, the available hardware does not allow the measurement of large forces on the system, so it is not possible to compare the evolution of the external force and the movement of the system. Simulation results are presented and external videos have proven the efficiency of the control law on a real device.

In this section, in order to simplify the expression, we consider only a one quadrant model: we consider the friction coefficients to be symmetric and not dependant on the sign of the velocity.

A. Control design

The control law used in this section is different from the one used in the section II-C. The previous control law is more efficient for systems controlled with a force sensor. In this case, the external force usually varies at a low frequency, which allows the use of an external force loop giving the reference of the inner velocity loop as shown in Fig. 4. However, in the case of a sensor-less control law, the estimated effort depends on the velocity for the estimation of the friction and the electric current for the estimation of the motor force. In our specific hardware, an inner electrical current loop is implemented and has a bandwidth 20 times higher than that of the bandwidth of the velocity loop.

In the case of large external forces, it is possible to apply an impedance control law on the system and control the motor torque to achieve the desired dynamic performance. Fig. 11 describes the implementation of the impedance control law and its model. The impedance control is used to control the motor torque according to the current velocity q (m/s), the motor electric current I m (A), and a reference position q 0 (m), as presented in Fig. 11a. The aim is to control the apparent impedance, described by Fig. 11b with M a , B, and K, around the position q 0 . τ m + τ e = M q + Fv q + Fc sign( q) [START_REF] Devie | Accurate force control and co-manipulation control using hybrid external command[END_REF] with M, Fv and Fc identified in the previous section. According to Fig. 11b, the desired mechanical equation is the following:

τ e = M a q + B( q -q0 ) + K(q -q 0 ) (25)
This equation gives:

q = 1 M a (τ e -B( q -q0 ) -K(q -q 0 )) (26)
According to [START_REF] Devie | Accurate force control and co-manipulation control using hybrid external command[END_REF] and ( 26), we get:

τ m = ( M M a -1)τ e + ( Fv - M M a B)( q -q0 ) + Fc sign( q) - M M a K(q -q 0 ) (27)
(27) is a classical impedance control law with a prediction term correcting the friction effects. In the following, in order to perform a co-manipulation task, we consider the following simplifications: K = 0 and q0 = 0. Thus, (27) becomes:

τ m = ( M M a -1)τ e + ( Fv - M M a B) q + Fc sign( q) ( 28
)
Let us now define the following external force observer: τe = M q + Fv q + Fc sign( q) -G τi I m (29)

where I m (A) is the current in the current loop of the motor.

For our model, we consider τm = G τi I m . With G τi = k t r, (28) and (29) give:

τ m = ( M M a -1) M q + M M a (F v -B) q + M M a F c sign( q)
-2( M M a -1)G τi I m (30)

B. Simulation result

Simulations were performed using the impedance control law (30). The parameters M a and B are calculated from the parameters of the control law of section II-C thanks to an equivalence criterion. The parameters M, Fv and Fc are the ones identified in section III-A, when q > 0.

The simulation is performed considering an identification error of 2% between the parameters used in the control law and the parameters used in the mechanic model. This error simulates the uncertainty of the identification.

The evolution of external force applied on the robot, its velocity and its position are plotted in Fig. 12. According to this figure, the control law allows a stable evolution of the velocity, proportional to the external force. Fig. 12: Simulation of the external force τ e (N), the velocity of the robot q(m/s) and its position q(m) function of the time (s), for a co-manipulation application without force sensor.

This result proves the efficiency of the force observer under these conditions. However, according to the hardware, the real conditions can be worse than the ones simulated. However, experiments were done with M a = 50 kg and B = 100 N/(m/s). The result can be seen on the videos.

It appears that the operator has to apply an effort higher that in the case , where a force sensor is used. It is coherent with the results of the previous section.

V. CONCLUSION

The aim of this paper was to show the limits of force observation in the case of co-manipulation tasks with an industrial robot. Firstly, the model of the robot was defined, and presented in a linear form. This model was used to perform identification of the robot parameters and observation of the external effort based on the current in the DC motor. It appears that the identification of the mechanical parameters has a large variance. The consequence of this dispersion is a large error on the force observation. This error coupled with the effects of static friction makes it unsuitable for precise sensor-less applications.

Then, an efficient control law was presented in order to control the robot for sensor-less applications with large efforts. Even if this application is less efficient that the one using a force sensor, due to noise and identification errors, it still gives acceptable results.

Fig. 1 :

 1 Fig. 1: Components of the EMPS robot.

Fig. 2 :

 2 Fig. 2: Components of the EMPS robot.

Fig. 3 :

 3 Fig. 3: Shape of the classical friction model in red, and the Stribeck friction model in blue, function of the velocity q

Fig. 4 :

 4 Fig. 4: Cascaded closed loop of speed and force in the case of co-manipulation.

Fig. 5 :

 5 Fig. 5: Measured values (blue) and filtered values (red) of the effort from the motor (top), the velocity (middle) and the acceleration (bottom)

Fig. 8 :

 8 Fig. 8: Force measured by the sensor τ e2 in red and estimated force τ e2 in blue

Fig. 9 :

 9 Fig. 9: Force measured by the sensor τ e in red and estimated force τ e in blue

Fig. 10 :

 10 Fig.10: Force applied by the motor τ m in blue, and estimated force W X in red and effort from the operator τ e in black.

Fig. 11 :

 11 Fig. 11: Definition of the impedance control law

TABLE I :

 I Identified values

	parameters	X 1	2 σ X 1	%σ X 1
	M 1 (kg) F + v1 (N/m/s) F -v1 (N/m/s) F + c1 (N) F -c1 (N)	115.7315 474.9727 563.8568 17.3811 11.5846	0.5765 17.5686 19.3088 0.3259 0.3363	0.4981 3.6989 3.4244 1.8750 2.9030

TABLE II :

 II Identified values

		parameters	X 1	2 σ X 1	%σ X 1
		M 1 (kg) F + v1 (N/m/s) F -v1 (N/m/s) F + c1 (N) F -c1 (N)	109 707 998 18.9 9.8	1.24 27.5 26.4 0.21 0.21	1.13 3.89 2.65 1.15 2.22
	0	2	4	6	8	10	12	14
				time (s)

ACKNOWLEDGMENT

The project is co-funded by the European Union. Europe is now investing in Pays de la Loire with the European Regional Development Fund.

Sylvain Devie received his Master's Degree in Advanced Robotics from the École Centrale de Nantes, and his Magister's Degree in Mechatronics from the École Normale Suprieure de Rennes in 2015. Since 2015, he has been a Ph.D. student in the Université de Nantes, working on the force control of robots in collaborative tasks.

Pierre-Philippe Robet received the Agrégation in electrical engineering from the Ecole Normale Suprieure, Cachan, France; the M.S. degree in satellite communication engineering from Surrey University, Guildford, United Kingdom; the Ph.D. degree from the Universit de Nantes et Ecole Centrale de Nantes, Nantes, France; in 1989, 1990, and 1996, respectively. Following receipt of the M.S. degree, he joined the International Maritime Satellite Organization to work in the eld of channel transmission modeling. Then, he was a teacher prior to begin his current research in automatic control and robotics at the University of Nantes and Ecole Centrale of Nantes. His current research topics include modeling, identication, control of electrical drives, and control of robots with force feedback.

Yannick Aoustin is Professor at University of Nantes, France. He earned his Ph.D. degree in Automatic in 1987 and his research degree for leading research and PhD students in 2006. He has been engaged in research for more than 25 years and his research interests include mechanical systems under actuated legged robots, bipedal, nonlinear observers, and biomechanics. He has publications on journals circulated internationally and premier conferences. He is also a topic editor of "International Journal of Advanced Robotic Systems" and member of the editorial board of the journal "Multibody system dynamics".

Maxime Gautier received the Doctorat d'Etat degree in robotics and control engineering from University of Nantes, Nantes, France, in 1990. Since 1991, he has been a Professor of automatic control with University of Nantes. He carries out his research with the Robotics Team, Institut de Recherche en Communications et Cybernétique de Nantes. His research interests include modelling, identification, and control of robots.