Ahmed Guellil
email: guellilamed@yahoo.fr

Tewfik Kernane

A NEW APPROACH OF POINT ESTIMATION FROM TRUNCATED DATA

Keywords: Dissection of mixture normals, EM algorithm, Minimum distance, Point estimation, Truncated data

We propose a new approach for estimating the parameters of a probability distribution from truncated data. We introduce …rst a new distance which measures the di¤erence between the variations of two discrete functions having the same support. We also introduce an auxiliary probability distribution depending on the same parameters of the considered distribution under study. We estimate the parameters of the later from truncated data. In a …rst method, we proceed by usual methods and in the second we use the minimum of the new distance. We illustrate the approach by some examples for discrete and continuous distributions. We apply the approach for dissection of mixture of two normal distributions.

Introduction

Point estimation is the most popular forms of statistical inference (see [START_REF] Lehmann | Theory of point estimation[END_REF]). We introduce in this paper a new statistical point estimation approach which found be useful in special practical situations such as truncated data and merged populations. The data are said to be truncated when measuring devices fail to report observations below and above certain readings. For example, truncated data frequently arise in the statistical analysis of astronomical observations ( see [START_REF] Efron | Nonparametric methods for doubly truncated data[END_REF]) and in medical data, and if the truncation is ignored this can cause considerable bias in the estimation (see [START_REF] Klein | Statistical challenges in comparing chemotherapy and bone marrow transplantation as a treatment for leukemia, Lifetime Data: Models in Reliability and Survival Analysis[END_REF]). There exist in the literature many approaches of estimation from "incomplete data" such as maximum likelihood based approach of the EM algorithm [START_REF] Hartley | Maximum likelihood estimation from incomplete data[END_REF], [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]), or nonparametric methods such as Kaplan-Meyer or Lynden-Bell estimators ( see [START_REF] Woodroofe | Estimating a distribution function with truncated data[END_REF], [START_REF] Efron | Nonparametric methods for doubly truncated data[END_REF]). The purpose of the present paper is to investigate other approaches based on variational aspects of probability distributions.

In estimation problems, we deal in general with three functions: a theoretical probability law f ( ; ) of a random variable X; depending on a parameter (real or vector valued), an empirical distribution f constructed from a sample of observations drawn from the random variable X; and an estimation b f ( ; b ) (where b is an estimation of ) obtained through the empirical law f . Usually, problems related to quality of the estimation b f are viewed through the theoretical properties of the estimator of : f is considered as a representative distribution of f; but in practice it is reduced to only few of its characteristics such as the mean and variance. The variational aspect of f is often neglected while its importance in our viewpoint. We can easily …nd, for instance, two distributions having the same support, mean and variance while their variations di¤er signi…cantly, or conversely having the same variations but their supports and characteristic parameters are di¤erent. But two distributions with same support and same variations in each subset of the support are necessarily the same. This suggest to look for a measure of proximity between variations of probability distributions having the same support.

One of the main tools used in parametric estimation is the method of minimum distance. We introduce here a new distance which proceed di¤erently from classical distances. It measures the di¤erence between variations of two discrete functions having same support. For estimation purposes, we build up a discrete distribution, with …nite support, associated to a sample of observations. The variations of this new distribution are assumed to be an empirical representation of those of the parent distribution restricted to the same support. In discrete case, we have only to use the empirical law of the observations. In continuous case, this could be done by considering the empirical distribution of the representative values (mid-classes or class means) after grouping the observations. The approach of estimation is then based on considering as a probability law the parent distribution restricted to the support that captures the data. In discrete case, we have simply to take the conditional distribution with respect to the observed values and we proceed analogously for the continuous case. We have to deal with two discrete probability distributions having the same …nite support, a theoretical distribution and its empirical representation with respect to the observations. The parameters of the former are those of the parent distribution and the aim is to estimate them from the …rst instead of the parent one as commonly used. We use …rst classical methods such as the method of moments or maximum likelihood principle. We proceed also by minimum distance using the new measure introduced here. We have then another way to parameter estimation giving new insights and perspectives as we will show in this paper. The setting that seems to us most suitable for illustrating our approach is the one of the estimation from truncated data. In usual practical problems, truncation can be on left or right or in either situations, and the "cut o¤" can be deterministic or random. In our approach, the truncation may be on any part of the range of the distribution so that the setting is more general. Also, classical approaches for truncated data are in general custom-made or subjective based methods. Instead, our approach is quite general and might be used in any situation where the underlying complete data come from a known family of distributions. We con…ne ourselves as a …rst presentation to deterministic truncations. Our intention is to present this new approach without exploring all possible areas of applications and extensions that it may generate.

In the subsequent section, we propose a variational distance between probability distributions. In Section 3, we de…ne a truncation of data and associated empirical and theoretical distributions and we use two di¤erent methods for estimation from truncation, a …rst method using traditional tools of estimations such as the method of moments and a second method based on the minimum of the new distance introduced in this paper. We illustrate the procedures by three examples: a binomial probability law, a normal distribution and a Gamma density function and we apply the two methods for dissection of two merged normal distributions. We present also an application to the estimation of the …rst trial value in the celebrate EM algorithm for incomplete data in the case of truncation. In Section 4, we present a basic feature of the new approach which prove the accuracy of the method. Finally, concluding remarks are made some pointing to possible extensions and applications.

Distance Based on Proportional Variations

As is usual is practice, given a sample of n independent and identically distributed observations, (x 1 ; x 2 ; :::; x n ) ; drawn from an unknown discrete random variable X falling in a discrete family of probability laws P = fp ; 2 R r g depending on a parameter (real or vector valued), i.e., p (x) = P(X = x) one can easily summarize the sample into k couples (y 1 ; f 1 ) ; (y 2 ; f 2 ) ; :::; (y k ; f k ) ; k n; where the y i are the di¤erent values taken by the sample and f is the empirical law

f j = n j n ; (1)
where n j represents the absolute frequency of the value y j ; j = 1; :::; k: Usually, it is hoped that f j p (y j );

in a certain probabilistic sense. However, we claim a more general assumption which is akin to the proportional allocation procedure for missing values, and take into account the variational aspects of probability distributions, so our aim is to have approximately

f i f j p (y i ) p (y j ) : (3)
Introduce the following distance of proportional variations between p and f d v (f ; p ) = X i;j2f1;:::;kg

f i f j p (y i ) p (y j ) : (4)
It turns out that this new distance, as we will show, measures the variations between probability distributions. One of its main powerful feature is that when using traditional distances of type P i2f1;:::;kg jf i p (y i )j we have to use the sample size n through the expression of f i = n i =n, but sometimes as for truncated data situations where measuring devices fail to report even the number of sample points in certain ranges, then the real size n is not known, but a truncated sample size n t is instead used. Using the ratios f i =f j will clear the e¤ect of the truncated sample size which can lead to considerable bias in the estimation. Furthermore, the new approach is related on the quality of the selected data and not mainly on large samples.

In continuous case also, any sample (x 1 ; x 2 ; :::; x n ) is summarized into k couples (y 1 ; f 1 ) ; (y 2 ; f 2 ) ; :::; (y k ; f k ) ; k n, which we hope contain all necessary informations about the variational aspects of the distribution of the continuous random variable X: We assume that this can be done uniquely, by grouping for example the sample in classes where the y i are the mid-classes (or class means) and the f i are the corresponding relative frequencies. This is justi…ed by the assumption that there exists an optimal number of classes which allows the empirical distribution function expressing the best the variations of the parent probability distribution. The proportional variational distance d v in this case between the density f (x; ) of X and its empirical law f is thus de…ned by

d v (f ; f ( ; )) = X i;j2f1;:::;kg f i f j f (y i ; ) f (y j ; ) : (5)
Let be given two functions f and g; d v associates for points x and y from the intersection of their domain of de…nitions, the quantity

d v (f; g) (x; y) = f (x) f (y) g(x) g(y) + f (y) f (x) g(y) g(x) : (6)
This measures the di¤erence between the variations of f and g on points x and y: Note that d v possesses the properties of symmetry and triangle inequality. But in the identity property

d v (f; g)(x; y) = 0 () f g; (7)
the equality between f and g must be understood in the sense that f and g have the same variations on the points x and y: From now on f will represent a theoretical probability law in both discrete or continuous cases and f shall represent the corresponding empirical law in both cases. Denote by = fx 2 R; f (x; ) > 0g the set of atoms of f or support. Let F be the algebra generated by the sets A = B \ ! where the ! are the Borel sets of R and B : For all A 2 F; we have

P (A) = Z A f (x; ) (dx); (8)
where is the Lebesgue measure on R.

In discrete case, we have

P (A) = X x2A f (x; ): (9)
For all i 1; we set i = ; F i = F and P i = P:

Let (n) = 1 2 ::: n ; F (n) = F 1
F 2 ::: F n and P (n) = P 1 P 2 ::: P n : The probability space (n) ; F (n) ; P (n) represents the space of samples of size n from the random variable X:

Let be given a classical distance d between two functions f and g which associates for points x and y from the intersection of their domain of de…nitions, the quantity

d (f; g) (x; y) = jf (x) g(x)j + jf (y) g(y)j (10)
Proposition 1 We have the following properties for the distance d v :

1. d(f; g)(x; y) = 0 =) d v (f; g)(x; y) = 0; the converse is not always true. 2. lim n!1 d v (f ; f ) = 0 in probability.
3. Let f and g be two functions de…ned on R and E R satisfying:

8 (x; y) 2 E E; d v (f; g)(x; y) = 0: (11) If Z R f d = Z R g d = 1; (12)
where is the Lebesgue measure on R, then

E = 0 =) f = g almost surely on R. (13)
Proof. 1. Follows directly from the de…nitions of d and d v :

2. Follows from the fact that whenever

lim n!1 d(f ; f ) = 0 then lim n!1 d v (f ; f ) = 0 in the same probabilistic notion of convergence. 3. Fix y 0 2 E; we have f (x)=f (y 0 ) = g(x)=g(y 0 ) for all x 2 E: This implies that Z E f (x)dx = 1 () Z E f (y 0 ) g(x) g(y 0 ) dx = f (y 0 ) g(y 0 ) Z E g(x)dx = 1: (14)
We deduce that f (y 0 ) = g(y 0 ) and the result follows.

Application to Truncated Data

Consider a sample of observations (x 1 ; x 2 ; :::; x n ) drawn from a theoretical probability law f ( ; ); depending on a parameter 2 R r : As usual, the data are summarized, in discrete or continuous cases (as shown in Section 2), into k couples (y 1 ; f 1 ) ; (y 2 ; f 2 ) ; :::; (y k ; f k ) ; k n; and let 4 = fu 1 ; u 2 ; :::; u m g a part from the set fy 1 ; y 2 ; :::; y k g ; m k; which we will call truncation.

De…ne the empirical distribution e f corresponding to a given truncation 4 by:

e f (x) = e f i if x = u i ; i = 1; :::; m; 0 otherwise, (15)
where the e f i satisfy the following set of proportional allocation equations

e f i e f j = f i f j ; (16)
for i; j = 1; :::; m and e f 1 + e f 2 + ::: + e f m = 1: De…ne the following auxiliary distribution from f ( ; ); which is based on proportional allocation as follows

h (x; ) = 8 < : f (x; ) f (u 1 ; ) + f (u 2 ; ) + ::: + f (u m ; ) if x = u i ; i = 1; :::; m; 0 otherwise (17)
Remark 2 If the truncation is random, that is, there exists a random variable T such that we observe, for example, the random variable X only if X > T or X < T; then the probability law used in (17) is replaced by the conditional law of X with respect to fX > T g or fX < T g respectively.

The distance d v and the auxiliary density h were found be useful for estimation problems in truncated data. Indeed, it is well known in classical estimation from truncated data (see [START_REF] Hartley | Maximum likelihood estimation from incomplete data[END_REF]) that missing values could be recovered by "proportional allocation" procedures, then the auxiliary density h; which is already based on proportional allocation, and the distance d v , which is based on proportional variations, will be intuitive and natural tools for estimation purposes from truncated data. The function h is a theoretical probability distribution depending on the same parameters of those of f . It has also the same support (atoms) as that of e f :

De…nition 3
We call e f and h( ; ) the empirical and theoretical distributions of a given truncation 4 = fu 1 ; u 2 ; :::; u m g from a sample of observations (x 1 ; x 2 ; :::; x n ) :

Methods of Estimation

As well known, two distribution having the same parameters (of central tendency for example) must fall in the same family to be similar. Classical distances do not give informations about variational aspects of probability distributions. So that, the distance d v will give an additional information about a sample of observations. This help making decisions about the accuracy of the estimation and having a measure of quality of the observations themselves. Indeed, the parameters are only a way for determining the theoretical distribution which …ts well the data. Whereas, probability distributions may also be characterized by both variations and supports. We notice that in general the classical and new approaches lead to approximately the same estimations. But in some cases they di¤er signi…cantly and this seems to be related to the situation where the estimation distribution obtained from classical approach does not …t well in variational viewpoint the parent distribution of the data. We will use mainly two methods of estimation. The …rst approach is similar to traditional ones such as the method of substitution or maximum likelihood principle, by considering e f as an empirical estimation of h( ; ): The second approach is a minimum distance estimation using the metric d v between e f and h: It was noticed, as will be shown in the examples below, that the two methods lead to approximately the same estimation results except few cases. Nevertheless, if the two estimations are signi…cantly di¤erent, it seems related to the quality of the selected data. However, We do not discuss the properties of the estimation by the minimum of distance d v ; since it is a complex problem which needs a separate thorough study. Our aim here is to validate new approaches that may be useful in some cases such as truncation in data or merged populations (see section 4). In order to test the performance of the proposed approach, we use synthetic datasets which are generated by simulation from three examples of probability law: binomial law, normal density and a Gamma distribution. The examples were selected from various simulation studies from di¤erent family of probability distributions and the two approaches have shown their e¤ectiveness and never deviate signi…cantly from the true parameter provided that the information contained in the truncation about the variational aspects of the parent distribution is su¢ cient. The reason for using synthetic datasets is that the true parameters for the synthetic datasets are known and the accuracy of the results obtained by using the two new methods can be compared.

Examples

Binomial distribution We generate a synthetic dataset of size 500 from a binomial law B(n; p) with n = 10 and p = 0:3: The data are summarized in the following table Our aim is to estimate the parameter p; with the knowledge of n = 10, from di¤erent truncation of the data.

For illustrating the two methods, consider the truncation 4 = f2; 3; 4; 5g with truncated sample size n t = 386: We have then a truncation proportion of Q = (n n t )=n = 22; 8 % in data. For the …rst method, the empirical distribution e f given the truncation = f2; 3; 4; 5g is given by e

f (2) = 108=386; e f (3) = 134=386, e f (4) = 97=386; e f (5) = 47=386 and e f (x) = 0 if x = 2 : Let f (x; p) = C x n p x (1 p) n x
the binomial probability law. The auxiliary distribution h( ; p) is given by:

h(x; p) = f (x; p) f (2; p) + f (3; p) + f (4; p) + f (5; p) if x = u i ; u i 2 f2; 3; 4; 5g 0 otherwise. ( 18
)
By the method of substitution, the estimation of p is obtained by solving the equation:

X u i 2f2;3;4;5g u i h(u i ; p) = X u i 2f2;3;4;5g u i e f (u i ) (19) that is 2 C 2 10 p 2 (1 p) 8 + 3 C 3 10 p 3 (1 p) 7 + 4 C 4 10 p 4 (1 p) 6 + 5 C 5 10 p 5 (1 p) 5 C 2 10 p 2 (1 p) 8 + C 3 10 p 3 (1 p) 7 + C 4 10 p 4 (1 p) 6 + C 5 10 p 5 (1 p) 5
= 1241 386 :

(20) Using a computer algebra package we obtain the result e p 1 = 0:3001.

For the second method, we have to search the value of the parameter p which minimizes the distance d v ; that is: Using computer algebra package, we obtain the result e p 2 = 0:2991: In the following table we present the estimations e p 1 from the …rst method using the auxiliary distribution, and e p 2 from the minimum distance approach using the distance d v ; of the parameter p; for known n; according to the truncation 4 = fu 1 ; :::; u m g considered. As previously said, the two estimations by the new approach, e p 1 and e p 2 ; are very accurate in all cases and close to each other, except few cases, since they are based on the same spirit. In practice, when the two estimations are close then we can take one of them or their mean. On the other hand, if they are signi…cantly di¤erent we can assert that the sample does not restore in a coherent way all the properties of the distribution from which it emanated. We then have to take the …rst estimation and use the second as a tool of decision on whether the estimation is credible or not.

min p d v ( e f ; h) = min p C 2
Normal distribution Consider a sample of size 400 drawn from a normal population with mean m = 0 and standard deviation = 1: Consider the data falling in 10 …xed class intervals as shown in the following table, with mid-classes u i and absolute frequencies n i Table 3 It should be noted that for the normal distribution the two estimations e m 1 and e m 2 by the new approach are somewhat close to each other for all the truncations considered, except the truncation n 12; however, for the latter, they are also close to the true value m = 0. For the …rst method ( e b 1 ), we remark that left truncations (n 1 10) are more accurate than right ones (n 11 21). Also, in general, it seems that the …rst method is better than the minimum of distance d v in this case of a continuous gamma density, but the distance d v may give a judgment about the accuracy and e¤ectiveness of the estimation from a given truncation. It means that if the two estimations are close then this will reassure us about the quality of the estimation.

Gamma probability density

Application to Dissection of Merged Normal Populations

We shall present in this section an application of the previous method used for truncated data in the situation where we have a mixture population of normal distributions.

Assume we have merged two samples of observations of sizes n 1 and n 2 generated from two normal distributions f 1 = N (m 1 ; 1 ) and f 2 = N (m 2 ; 2 ), with m 1 6 = m 2 . By assuming that 1 and 2 are known, our aim is to estimate the means m 1 and m 2 and also the merging proportion of each population.

In classical methods, we use the merged distribution f = f 1 + ( 1) f 2 and we estimate the parameters ; m 1 and m 2 using for example the EM algorithm. We will use a method based on truncations. The main idea being to split the range of the merged sample into three suitably chosen parts. A central part where the observations are highly merged, a left and right truncated parts where the observations become mainly from one of the distributions considered. If for example m 1 < m 2 , then to estimate m 1 we have to use the chosen right truncated part (left truncation 4).

The procedure is summarized as follows:

1. We compute the sample mean m g of the merged observations.

2. For determining the location of the two means m 1 and m 2 , we compute the empirical standard deviation S l of the observations less than m g ; and S r for those that are greater. Assume that S l < S r , in this case if 1 < 2 then we deduce that m 1 is situated on the left of m g . Otherwise, it will be assumed to be on its right. We follow the same idea for the case S l > S r : If 1 = 2 we pass directly to the third step.

3. Assume that m 1 is on the left. It is well known that for a normal distribution N (m; ) we have P (]m ; m + [) ' 0:68: We hope that on the left of sup l = m g 2 the number of observations generated from N (m 2 ; 2 ) is negligible, and on the right of min r = m g + 1 the number of observations generated from N (m 1 ; 1 ) is also negligible. Hence, to estimate m 1 , we consider only the part of observations situated on the left of m g 2 , and to estimate m 2 we consider the part situated on the right of m g + 1 :

The following example will provide some feel for the accuracy of the procedure.

Example

We consider the case where 1 = 2 : consider two samples of observations generated from N (m 1 ; 1 ) and N (m 2 ; 2 ); where m 1 = 1:3 and m 2 = 2:4, with known 1 = 2 = 1 and sizes n 1 = 300 and n 2 = 200: We combine them to obtain a merged sample of size n = 500: We have chosen the distributions in such a way that the histogram (Fig. 1) of the merged sample does not show directly the existence of a mixture of two distributions. When the histogram of the merged population is bimodal the situation is more easier, since when taking a suitably left (or right) part we get more accurate estimation from the situation that this part will have a negligible number of observations from the second distribution. It should be stressed that the histogram is one modal and does not show at …rst glance any mixture situation. Following the steps of the procedure we begin by calculating the mean of the resulting merged sample and we obtain m g = 1:8046. Since the standard deviations are assumed to be equal then we compute directly sup l = m g 2 = 0:8046: By grouping the observations on the left of sup l (which constitute the chosen right truncated part) in 7 classes we obtain the following table: The sample mean of the observations on the left of sup l is given by u l = 0:1483: Using the …rst method we have to solve on m the following formula

u 1 exp h (u 1 m) 2 2 2 i + y 2 exp h (u 2 m) 2 2 2 i + ::: + y k exp h (u k m) 2 2 2 i exp h (u 1 m) 2 2 2 i + exp h (u 2 m) 2 2 2 i + ::: + exp h (u k m) 2 2 2 i = u l : (22)
we obtain the estimation e m 1 = 1:2646: By deleting the …rst value u 1 which has a weak frequency n 1 = 1; that is using the truncation 4 = fu 2 ; u 3 ; u 4 ; u 5 ; u 6 g ; (we compute again u l = 0:1734) we obtain a better estimation e m 1 = 1:3011; which is very close to the true value m 1 = 1:3. = 2:3971: The mixture proportion can easily be estimated using the formula: e m 1 + ( 1) e m 2 = m g :

Estimation of the initial value in EM Algorithm

The initial starting value is of great importance in convergence behavior of algorithms such as EM Algorithm. Usually, as for the latter, the initial trial value is guessed. Surprisingly, we will show that our procedure gives an estimation of the starting value instead of having to guess. The approach will be illustrated by the following classical example which was the basis of the EM algorithm.

Example of Hartley (1958) revisited

Hartley (1958) used an algorithmic procedure to estimate the parameter of a Poisson distribution from data on the pollution of seeds of Phleum pratense by the presence of noxious weed seeds quoted from [START_REF] Snedecor | Statistical Methods[END_REF] and truncated them by missing the frequencies of the values 0 and 1 as shown in the following table (

Hartley [START_REF] Hartley | Maximum likelihood estimation from incomplete data[END_REF] has guessed the frequencies of the missing values 0 and 1 by taking n 0 = 4 and n 1 = 14 and after 4 steps of his algorithmic procedure, which has been the basis of the well known EM algorithm for incomplete data [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], has reached the estimation b = 3:026 (see table 1 p.177 [START_REF] Hartley | Maximum likelihood estimation from incomplete data[END_REF]). Using the …rst method we get the estimation e 1 = 3:1149: And by proportional allocation procedure we can see that the frequencies we get are n 0 = 4:29 and n 1 = 13:38 which are close to the guessed values.

A Basic Feature of the New Approach

The preceding results have shown the e¤ectiveness of the new approach and worked well in simulation experiments. Furthermore, the proposition below will give an insight of a major feature of the new approach by considering the one parameter exponential family. We will prove that for all truncation considered formed by more than two points, from a sample of observations , if the ratios of the relative frequencies of the u i are exact, then we may obtain the true value of the parameter. We may conjecture that when considering an arbitrary law of probability depending on r parameters, such that we have a truncation composed by r + 1 points having exact empirical ratios of the relative frequencies then we obtain the true values of the r parameters.

Proposition 4 Consider a probability distribution f from the one-parameter exponential family with density

f (x j ) = k(x) exp[ T (x) A( )]; (23)
where 2 R is the parameter, T a statistic, k(x) a function of x and A is a function of the parameter : Assume that we wish to estimate the parameter : If we consider a truncation having two points x and y with empirical frequencies f 1 and f 2 satisfying f 1 =f 2 = f (x j )=f (y j ), then, using the approach considered here, we obtain the true value of :

Proof. 1. If we consider the minimum of distance d v the result is immediate. 2. Consider now the second method to estimate m. Consider two values x and y from the exponential family with density given by (23); with e being the estimation by the new approach, and assume that their empirical frequencies f 1 and f 2 are such that

f 1 f 2 = f (x j e ) f (y j e ) : (24)
We obtain Remark 5 Note that the frequencies f 1 and f 2 need not be exact, that is f 1 may be di¤ erent from f (x j ) and also f 2 ; but we require only that their ratio is equal to the theoretical one f (x j )=f (y j ):

and also we get the true value e p 2 = 0:3.

Gamma distribution

Consider the Gamma probability distribution G(a; b) with a = 10 and b = 5: Assume that a is known and we wish to estimate b: Consider in table 5 the truncation 4 = fu 3 ; u 8 g with u 3 = 30:13 and u 8 = 60:02. We have the following value of the ratio f (u 3 j b) =f (u 8 j b) 0:7993 (the result is an approximate result since for probability density functions it is di¢ cult to get an exact rational value but we will show that the estimations are very close to the true value). Consider the absolute frequencies n 3 = 79:93 (or 80) and n 8 = 100 for the values u 3 = 30:13 and u 8 = 60:02 respectively. We have then n 3 =n 8 f (u 3 j b) =f (u 8 j b) : From the …rst method, we compute u = 46:7438 and solve on b the following equation (30:13 46:7438) 30:13 9 exp( 30:13=b)

parameters a = 10 and b = 5 :

 5 Consider a sample of size 400 drawn from a Gamma distribution G(a; b) with density given by f (x j a; b) Consider the data falling in 13 …xed class intervals as shown in the following

4 size e b 1 e b 2

 42 show the estimations e b 1 and e b 2 ; by the two methods, of b; with known a = 10; according to the truncation 4 considered Table 6.The estimations e b 1 and e b 2 by the new approach of the parameter b of the Gamma probability distribution G(a; b) with b = 5 and known a = 10. n

Fig 1 .

 1 Fig 1. Merged histogram of two normal distributions N (1:3; 1) and N (2:4; 1):

1 = 1 :

 11 Using the distance d v we obtain for all the truncation e m (dv) 244 and by deleting u 1 we get the value e m To estimate m 2 we consider the part situated on the right of inf d = m g + 1 = 2:8046: Grouping the observations on the right of inf d (which constitute the chosen right part) in 7 classes we obtain the following table:

u i 2

 2 The sample mean of the observations on the right part is given by u d = 3:5227: Using formula (22) with u d , we obtain the result e m 2 = 2:245: Deleting the extreme values u 1 and u 6 we obtain e m 2 = 2:4124: Using the distance d v for all the truncation we get e m (dv) 2

u = xf 1 + 1 A 1 A

 111 yf 2 = xk(x) exp e T (x) + yk(y) exp e T (y) k(x) exp e T (x) + k(y) exp e T (y) : Then, we solve on the following equation: 0 @ x xk(x) exp e T (x) + yk(y) exp e T (y) k(x) exp e T (x) + k(y) exp e T (y) k(x) exp ( T (x)) + 0 @ y xk(x) exp e T (x) + yk(y) exp e T (y) k(x) exp e T (x) + k(y) exp e T (y) k(y) exp ( T (y)) = 0; after straightforward algebra we obtain (x y) exp e T (y) + T (x) + (y x) exp e T (x) + T (y) = 0; yielding the true value e = : The proof is complete.

4. 1 we obtain the exact result e p 1 = 0: 3 :

 113 Examples4.1.1 Binomial distributionConsider again the binomial distribution B(n; p) with n = 10 and p = 0:3 and assume n is known and we wish to estimate p: Assume we have the following truncation with only two points 4 = f0; 1g : The exact ratio of their probability distribution is given by f (0; p)=f (1; p) = 7=30; which is a rational value that will simplify the example. Choose the absolute frequencies of the two values considered as being n 1 = 7 and n 2 = 30 for the values u 1 = 0 and u 2 = 1 respectively, in order for having f 1 =f 2 = f (x; p)=f (y; p) = 7=30: Using the …rst method we have to solve the following equation on p Using the second approach, that of the minimum of distance d v ; we have to solve min p d v ( e f ; h) = min p

Table 1 .

 1

	y i 0	1	2	3	4	5	6 7
	n i 15 71 108 134 97 47 23 5

Table 2 .

 2 The estimations e p 1 and e p 2 by the new approach of the parameter p of the binomial probability law B(n; p) with p = 0:3 and known n = 10.

		Truncated	Proportion of
	n 1 2 3 f0; 1; 2; 3; 4; 5; 6; 7g 4 f0; 1; 2; 3; 4g f0; 1; 2; 3; 4; 5g 4 f1; 2; 3; 4; 5g 5 f1; 2; 3; 4g 6 f2; 3; 4; 5g 7 f2; 3; 4g 8 f3; 4; 5g 9 f0; 3; 4; 5g 10 f1; 3; 4; 5g 11 f0; 2; 3; 4; 5g 12 f0; 1; 3; 4; 5g 13 f0; 1; 4; 5; 6; 7g 14 f2; 4; 5; 6; 7g 15 f0; 4; 5; 6; 7g 16 f0; 2; 3; 4g 17 f0; 5; 6; 7g 18 f0; 5g 19 f1; 4; 5; 6; 7g 20 f0; 1; 2; 5; 6; 7g	sample size truncation Q (%) 425 15 472 5:6 500 0 457 8:6 410 18 386 22:8 339 32:2 278 44:4 293 41:4 349 30:2 401 19:8 364 27:2 258 48:4 280 44 187 62:6 354 29:2 90 82 62 87:6 243 51:4 269 46:2	e p 1 0:2928 0:2952 e p 2 0:2926 0:2952 0:2978 0:3053 0:2922 0:2877 0:2922 0:2868 0:3001 0:2991 0:3052 0:3049 0:2906 0:2901 0:2932 0:2952 0:2872 0:2870 0:2983 0:2952 0:2892 0:2952 0:2945 0:3013 0:3071 0:3077 0:3018 0:3071 0:3009 0:2952 0:3011 0:3014 0:2937 0:2937 0:2935 0:3028 0:2985 0:2981

Table 4 .

 4 In the following table we show the estimations of m according to the truncation 4 considered, where e m 1 represents the estimation by the …rst method and e m 2 that of the minimum of the variational distance d v The estimations e m 1 and e m 2 by the new approach of the parameter m of the standard normal distribution with known standard deviation = 1.

		.
	i	1	2	3	4	5	6	7	8
	u i	2:5808	2:057	1:5331	1:0092	0:4853 0:0386 0:5625 1:0863
	n i	5	8	23	48	71	89	72	43
	9	10	11
	1:6102 2:1341 2:658
	25	10	6
	n 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 fu 1 ; u 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 g 394 size 4 161 fu 6 ; u 7 g 119 fu 4 ; u 5 g 79 fu 2 ; u 3 ; u 4 g 78 fu 8 ; u 9 ; u 10 g 232 fu 5 ; u 6 ; u 7 g 204 fu 6 ; u 7 ; u 8 g 142 fu 3 ; u 4 ; u 5 g 166 fu 3 ; u 5 ; u 7 g 140 fu 7 ; u 8 ; u 9 g 150 fu 2 ; u 3 ; u 4 ; u 5 g 231 fu 3 ; u 4 ; u 5 ; u 6 g 229 fu 6 ; u 7 ; u 8 ; u 9 g 300 fu 5 ; u 6 ; u 7 ; u 8 ; u 9 g 348 fu 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 g 371 fu 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 g	e m 1 0:1039 0:0003 6 10 6 e m 2 0:1040 0:1631 0:1769 0:2436 0:2181 0:0518 0:052 0:1284 0:1319 0:0529 0:0665 0:0574 0:0593 0:0323 0:0767 0:0845 0:1179 0:1092 0:1137 0:0402 0:0165 0:0119 0:0520 0:0216 0:0515 0:0325 0:0592 0:0521 0:0651

  table, with mid-classes u i and absolute frequencies n i :

	Table 5.
	i	1	2	3	4	5	6	7	8	9	10
	u i 18:18 24:16 30:13 36:11 42:09 48:07 54:04 60:02 66:00 71:97
	n i	5	16	31	45	69	62	53	51	36	11
	11	12	13
	77:95 83:93 89:91
	9	5	7
	In the following table we

  12 g 393 4:9575 4:7758 10 fu 1 ; u 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 ; u 12 ; u 13 g 400 4:9947 4:8856 11 fu 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 ; u 12 ; u 13 g 395 5:0047 4:9695 12 fu 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 ; u 12 g

	1 2 3 4 5 6	97 4:5856 4:5058 166 5:0985 4:8489 228 5:0488 4:7949 281 5:0174 4:7554 332 5:1145 4:7949 fu 388 4:9666 4:7577 fu 1 ; u 2 ; u 3 ; u 4 g fu 1 ; u 2 ; u 3 ; u 4 ; u 5 g fu 1 ; u 2 ; u 3 ; u 4 ; u 5 ; u 6 g fu 1 ; u 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 g fu 1 ; u 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 g
	13 14 15 16 17 18 19 20	fu 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 g fu 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 g fu 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 g fu 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 ; u 12 g fu 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 g fu 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 g fu 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 ; u 12 g fu	383 4:9901 4:7290 367 4:9973 4:7577 336 4:9665 4:7290 341 4:9373 4:7350 327 5:0200 4:7089 291 4:8478 4:7290 227 4:8117 4:7269

1 ; u 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 g 368 5:1386 4:8174 7 fu 1 ; u 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 g 379 5:0158 4:7949 8 fu 1 ; u 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 g 388 4:9790 4:7949 9 fu 1 ; u 2 ; u 3 ; u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 ; u 6 ; u 7 ; u 8 ; u 9 ; u 10 ; u 11 ; u 12 ; u 13 g 234 4:9247 5:3153 21 fu 7 ; u 8 ; u 9 ; u 10 ; u 11 ; u 12 ; u 13 g 172 4:8494 5:3153 10

Table 7 .

 7

	u i	1:5589	1:1294	0:6998	0:2703	0:1593 0:5888
	n i	1	3	6	17	24	41

Table 1

 1

	in Hartley

Concluding Remarks

In the foregoing study, we have presented a new statistical point estimation method which found be useful in truncated data situations. A new distance between probability distributions was introduced. It measures the di¤erence between the proportional variations between two given probability distributions. We introduced a new auxiliary distribution based on a truncation, from a chosen family of probability distributions. This new distribution will have the same parameters to estimate as the parent one.

We use then statistical methods to estimate the parameter of the variable under study using the empirical and the new auxiliary distribution in the region that captures the data, from which we determine the corresponding parent distribution. The later is the estimation by the new method. Using the new distance introduced we also estimate by the minimum distance approach. We have obtained a result which state that if we have to estimate the parameter of a probability distribution from the one parameter exponential family, then it su¢ ces to have two points with exact ratio of frequencies, that is equal to the theoretical one expressed by the ratio of the value of the density function on these two points, to obtain the true value of the parameter. We have conjectured that if we have in general r parameters, then it su¢ ces to have r + 1 points with exact ratios of their frequencies to obtain the r true parameters exactly. The later result need to be proved rigorously in a general setting for other distributions than the class considered. A large comparative study between the classical and new methods should also be investigated. We also shown the possible and interesting perspective of discrimination of merged populations and mixture of distributions by the new method by considering the case of two normal merged populations. It would be interesting to investigate this area for more than two distributions and to other standard families.


I Mancisidor
email: imancisidor@ideko.es

J Munoa

R Barcena

OPTIMAL CONTROL LAWS FOR CHATTER SUPPRESSION USING INERTIAL ACTUATOR IN MILLING PROCESSES

Keywords: Chatter, Active Control, Inertial Actuator

Actually, active control devices are considered as one of the most suitable chatter suppression methods for heavy duty machining operations. Generally, they are based on the generation of a reaction force over the main structure, controlled by a feedback control law, which depends on the real time vibration measurement. Several active control strategies have been proposed in the literature and due to the limited space offered by the machines for these actuators, the active control system has to be optimized, including the control law. In this work, the most employed strategies will be compared by means of a verified mechatronic model and the optimal law will be sought for each zone of stability lobes. Finally, experimental tests will be performed for verifying the results obtained by the mechatronic model simulations.

INTRODUCTION

The presence of self-excited vibrations, known as chatter, in machining processes is a classic problem that limits the productivity. These vibrations, caused principally by the regenerative effect [Tlusty 1957;[START_REF] Tobias | [END_REF], prevent obtaining required surface finishes and decrease the life of tools and mechanical components of the machine. The location of a tuned vibration absorber in the structure is a widely studied solution for increasing the stability margin [Den Hartog 1956;Sims 2007]. However, the employment of a passive absorber is not feasible in processes where the dynamics of the machine change according to the working position during machining. Therefore, the integration of active dampers has been widely analyzed for many years, as one of the most effective chatter suppression methods, due to their adaptability to variable conditions. In 1970, Cowley and Boyle proposed to use an electromagnetic inertial actuator with an accelerometer in order to introduce active damping into the structure of a machine tool [START_REF] Cowley | [END_REF]]. Some years later, Ehmann and Nordmann [START_REF] Ehmann | [END_REF]] recommended introducing inertial actuator in large machine tools with vibration problems. Among inertial actuators, electromagnetic devices are usually employed [Loix 2004;[START_REF] Sabalza | [END_REF][START_REF] Baur | [END_REF]], although some hydraulic actuators have been proposed in the literature [START_REF] Brecher | [END_REF]Brecher 2012]. However, the hydraulic actuators present several disadvantages, such as time delays to reach the required pressure, leak oil, hysteresis effects, maintenance problems or other nonlinearities. These issues make the electromagnetic to be the most employed technology to develop inertial actuators. When these active control systems are designed, one of the most important parameters is the ratio between the force capability and the occupied volume. The reason is that generally machine tools offer reduced spaces for the actuator location close to the tool and, hence, their volume has to be reduced. Therefore, the performance of the actuator has to be exploited to the fullest, which depends partially on the active control strategy employed for calculating the force. Therefore it is important the optimization of the performance of the control law. Several active control laws have been proposed in the literature, which can be classified under two general categories [START_REF] Preumont | Vibration Control of Active Structures: An Introduction[END_REF]]: feedback and feedforward. Generally, feedback active controls are applied in machine tools, since feedforward controllers require a reference disturbance signal, which is very complicated to achieve when it is referred to chatter vibrations. In a typical feedback control loop, a sensor will measure the total response of the mechanical structure, which is then fed via the controller to the secondary actuator, creating a mechatronic system (see Fig. 1). Among feedback control laws, some authors have proposed algorithms which require the knowledge of the system behavior based on a mathematical model. Some of them, such as Virtual Passive Absorber (VPA) and Virtual Passive/Active Absorber (VPAA) have been designed for chatter suppression [START_REF] Huyanan | [END_REF][START_REF] Bilbao-Guillerna | Control laws for chatter suppression in milling using an inertial actuator[END_REF]. These strategies can benefit the control since efforts of actuators can be focused on the modes or frequencies of interest, not wasting over all vibrations. However, in some cases, the elaboration of a mathematical model can be a complicated task and, at best, can only be a low-dimensional approach of the actual system. Moreover, for certain machine tools, the model should be modified for each position of the tool over the workspace. Model free algorithms have the advantage that the reduction of vibrations can be achieved without a model of the structure and, hence, they can adapt their behavior in cases where the dynamic properties are changing inside the workspace. Therefore, they are prepared to control a wide range of frequency, which is usually limited by a high-pass and a low-pass filter introduced by the user. The objective of typical feedback control laws found in the literature is to reduce the vibration of a mechanical system by automatic modification of the system's structural response [START_REF] Fuller | [END_REF]]. In this way, the modification of the main dynamic parameters has been sought using position, velocity and acceleration feedback [START_REF] Brecher | [END_REF]]. Recently, a new model free control strategy, called Delayed Position Feedback, has been proposed [START_REF] Munoa | [END_REF]]. Such control law is based on the reduction of the regenerative effect and theoretically it can improve the stability. One of the most important effects which limits the performance of feedback controller in mechanical systems is the unmodelled system delay. Such phase shift may arise because of the dynamic response of the sensors and actuators being used, or may be due to time delays in the controller. If either position or acceleration feedback are implemented in a control loop with a system delay, the effect of the delay changes considerably the damping of the system [START_REF] Fuller | [END_REF]]. This fact modifies significantly the behavior of the control and, in some cases, the system can become unstable. In contrast, when velocity feedback is employed, the delay has only a small effect on the effective mass and stiffness and, hence, it is thus seen to be a more robust control strategy. Meanwhile, when the novel strategy (Delayed Position Feedback) is employed, the force is introduced with the tooth passing period delay, which is usually higher than the existing delays due to the actuator and controller. Therefore, these undesired delays can be compensated, reducing the tooth passing period delay of the feedback control. The present work performs the comparison between these four model-free feedback control strategies, by means of the milling mechatronic model [START_REF] Munoa | [END_REF]]. In the first part, their mathematical explanation is provided briefly and then, they are compared theoretically. Finally, they are implemented in a real cutting process, showing which control strategy is the optimum for each zone on a stability lobe.

MATHEMATICAL FORMULATION

This section presents the mathematical development of the analyzed control laws. Both the frequency domain and time domain formulations can be employed in order to study the stability of a cutting process [START_REF] Munoa | [END_REF]]. The regenerative effect is clearly explained by frequency formulations, although non-linearities cannot be considered [Zatarain 2004;Stepan 2011]. Therefore, the present work has considered the formulation of the mechatronic model [START_REF] Munoa | [END_REF]]. Such mechatronic model simulates the cutting process, where the influence of inertial actuator control loop is integrated into the regenerative effect. It is based on time domain simulations, which permit the inclusion of non-linearities in its formulation. In this way, such model can be employed to simulate easily the possible limitations and saturations of the inertial actuators, such as force or stroke saturations.

Time domain formulation

The present subsection summarizes the formulation employed by the mechatronic model in order to simulate the milling cutting process with active inertial actuator effect. The general time domain equation is considered,

          t t t t t q p act, c, F F r K r C r M           (1)
where r={rp, rq, rm} T . The subscripts mean the location of each parameter, where p is the cutting point, q is the actuator location and m is the measurement point, for cases where the actuator and sensor are not collocated (see Fig. 2). Following the linear formulation [Altintas 2012], the cutting forces can be defined as,

                                ) ( sin ) sin( 1 T t f K K b K F F F p j z ac rc tc p a,j r,j t,j r u Δ     j ae re te g S K K K                  (2)
If the tangential, radial and axial forces (Ft, Fr, Fa) of the inserts are projected onto the xyz Cartesian axes:

                                 1 0 1 0 , , , , c , Z j p a,j r,j t,j j Z j p j z j y j x p F F F F F F t   β F (3) where                                sin cos 0 ) cos( cos ) cos( sin ) sin( ) sin( cos ) sin( sin ) cos( , j j j j j j j β
The actuator force Fact performed in point q, depends on the selected control law and is analyzed in the next subsection. Assuming proportional damping, if the modal transformation is performed, the equation ( 1) is as follows:

          t t t t t q q p p , act T , c T F Q F Q η K' η C' η M'           (4) being M'=Q T MQ, C'=Q T CQ and K'=Q
T KQ, the normalized and diagonalised modal mass, damping and stiffness, respectively and Q={Qp, Qq, Qm}. Then, it can be simplified for each mode,

        i i i i i i i i m t P t η t η t η     2 , n , n 2       (5)
where mi, ωn,i and ξi are the modal mass, the natural frequency and the relative damping ratio of ith mode of structure, while Pi is the modal force. Modal displacements are calculated numerically based on state space formulation [START_REF] Bediaga | [END_REF]]:

          i i i i i i i i i m t P t η t η t η t η                             1 0 2 1 0 , n 2 , n        (6)           i k i i i k i k i i i i i k i k i m t P t t t t η t η A A A A t η t η                                   2 , n 2 , 22 , 21 , 12 , 11 1 1 2       (7)
Then, the Cartesian displacements are calculated by the modal vector matrix and they will be used in the next step.

    1 T 1    k k t t η Q r (8)
In this way, the results can be showed in graphs. In particular, Fig. 3 shows the results for a chatter-free and for a chatter case.

Model-free control laws

This subsection analyses the most employed control laws when inertial actuators are used for chatter suppression. Four control strategies have been selected considering different terms of the general regenerative dynamic equation (4). In order to show clearly the effect obtained by each control algorithm, the system is analyzed in frequency domain. In such formulation, the zero order approximation has been considered. When the stability of the linear system is analyzed in frequency domain, only the dynamic part has influence, so equation ( 2) can be simplified as

    ) ( ) ( 1 , 1 0 T c,                            Z j j p ac rc tc j p g t K K b K t    r u β F Δ (9)
The different directional coefficients can be included in a matrix A(t) which depends on the position of the mill, the lead angle and the relation between the cutting coefficients [Altintas 2012;[START_REF] Munoa | [END_REF]]. Moreover, instead of using b, axial depth of cut ap will be used because it is usually the user data.

             t t t a K t p p tc p r r A F ) ( sin p c, ( 10
) If it is developed in frequency domain,          i i i , act T , c T 1 q q p p F Q F Q η    Φ (11)
where the frequency response function (FRF) of the system Φ is diagonal and defined as

    1 2 ' i ' ' i      C K M    Φ (12)
and the cutting force in frequency domain is defined as

   η Q A F     c i 0 p c, 1 sin i    e a K p tc p (13)
When Direct Velocity Feedback (DVF) is employed, the control forces appear as viscous damping and the stability lobes increase their limit depth of cut. Such strategy is based on the measurement of vibration velocity and its negative feedback, as shown in (6).

      η Q G F r G F         m q m q t t   i i : DVF , act act,  (14)     p p m q c, T T 2 i F Q η Q G Q C K M         (15)
where G=ga•gv, being ga and gv the actuator gain and the gain matrix introduced by the user, respectively. In a general case, the term related to the actuator gain (Qq T GQm), will have cross values and, hence, all the force will not be introduced directly as viscous damping of each mode. However, if the active damping term is a linear combination of the mass and stiffness matrices, the equation ( 15) can be decoupled.

K M G       (16)     i i i i i i P g c k m         i 2 (17)
These equations show, that in the most general case, a model of the system is required in order to diagonalize the active damping term.

However, in most of cases it can be considered that the frequency range of interest is dominated by a single mode i. Therefore, the actuator force term can be simplified as a scalar and it can be considered as proportional force

                         i m z i m y i m x vz vy vx a i q z i q y i q x m q g g g g , , , , , , T 0 0 0 0 0 0       Q G Q m z q z vz a m y q y vy a m x q x vx a i g g g g g g g , , , , , ,           (18)
When the active system is non-collocated, i.e. the m and q are not the same point, some instabilities can appear depending on the gain sign. It occurs when the signs of the modal vector values are different in both points. However, in general collocated systems (are usually recommended [START_REF] Preumont | Vibration Control of Active Structures: An Introduction[END_REF]] and then, the equation can be simplified:

  2 , 2 , 2 , T q z vz q y vy q x vx a i m q g g g g g         Q G Q (19)
If the same gain is introduced for the three directions (gvx=gvy=gvz) the gain is proportional to the modulus of the modal displacement. Moreover, the higher is the modal displacement on the actuation point, the higher is the introduced viscous damping. Currently, tridimensional actuator has not been developed in the literature, so the actuation in the three Cartesian directions is not possible when only one actuator is employed. When a biaxial actuator [START_REF] Munoa | [END_REF]] is used, all the force capability can be exploited if the modal displacement direction is mainly over a plane, as occurs usually in ram type machine tools. The single axis actuator can be also employed in these cases and if the actuation direction matches with the mode direction, the force can be also exploited totally. If these directions do not coincide, the actuator can also work properly, but some force will be wasted. Nevertheless, the complexity is highly increased in cases where modal interaction occurs in the frequency range of interest and a complete model of the system is required, since cross terms appear and the sign and values of the gains have to be deeply studied. The same problems appear with the other control strategies. Direct Acceleration Feedback (DAF) is similar to DVF, but the acceleration signal is fed back instead of the speed. As equation ( 22) shows, a proportional gain with DAF control law modifies the modal mass and, therefore, the natural frequency can be increased or reduced, depending on the gain sign. A similar effect is caused by Direct Position Feedback (DPF) algorithm, due to the modification of stiffness (see equation ( 25)).

      η Q G F r G F        m q m q t t 2 act, act, i : DAF     (20)     p p m q c, T T 2 i F Q η C K Q G Q M         (21)     i i i i i i P c k g m         i 2 (22)       η Q G F r G F        m q m q t t  i : DPF act, act, (23)     p p m q c, T T 2 i F Q η C Q G Q K M         (24)     i i i i i i P c g k m         i 2 (25)
The presented control laws are focused on the modification of dynamic parameters. However, the target is changed when the novel Delayed Position Feedback (DelPF) control strategy is employed. The objective of such control law is the distortion of the regenerative effect, by the reduction of the engagement between successive waves. It is based on the feedback of the previous tooth vibration, which is obtained by the introduction of tooth passing period delay feedback. The mathematical formulation is very simple for orthogonal cutting processes [START_REF] Munoa | [END_REF]]. However, the formulation is complicated for milling operations, since the directional matrix has to be considered when the commanded force is calculated. The force commanded by DelPF can be considered as

        η Q G F η Q G r G F      c i act, act, i : DelPF            e t t t m q m m q (26)
In order to see the influence of the direction matrix on the gain of the new control strategy, the effect is analyzed in frequency domain formulation.

  η GQ Q Q A Q η            m T q p T p tc e e a K      c c i i 0 p 1 sin 1 - Φ (27)
The force direction is very important in this control law since it is focused on the regenerative part distortion and, hence, the gain introduction is not as simple as in the other strategies. The directional matrix has to be taken into account inside the gain, while the same gain gv should be introduced for all Cartesian directions.

T 0 T p sin m p p q tc v a a K g g Q Q A Q Q G   (28) Then,    η Q A Q η v a p T p tc g g e a K     1 1 sin c i 0 p    1 - Φ (29)
Therefore, the perfect behavior of the control law could avoid chatter for any cutting process when the regenerative part is completely removed (gagv=1). Moreover, since the acceleration and position are counter phase signals, the acceleration signal can be directly employed in cases where the vibration measurement is performed by an accelerometer. It is named as Delayed Acceleration Feedback (DelAF). In this case, as equation ( 13) shows, the feedback has to be negative and the gain should be much smaller.

   η

Q A Q η v a p T p tc g g e a K 2 c i 0 p 1 1 sin c         1 - Φ (30)
Therefore, it is clear that the directional matrix have to be multiplied when the gain is introduced on the controller, in order to reduce properly the regenerative effect. In conclusion, the actuator force is calculated as follows:

         t a K g g t T m p T p q tc v a q r Q Q A Q Q F 0 p , act sin (31)
where A0 only depends on geometrical parameters:

  Z K K f e s ac rc , , , , 0    A (32)
Tab. 1: Dynamic parameters of the machine when Y=760mm and Z=550mm.

If the equation ( 31) is analyzed, it can be observed that when the actuator or measurement does not coincide with the cutting point, it is necessary the knowledge of their modal displacement relation. It means that a model of the system would be required. However, when these points are very close (p≈q≈m), it can be assumed that the modal matrices are the same and then, identity matrix is the result of their multiplication. Moreover, the terms, Ktc, ap and κ, can be considered constant, so the gain only would depend on the relation of directions and values of the directional matrix terms (33). The same mathematical development can be done for the DelAF.

        t g g t v a q r A F 0 , act (33)
However, when the cutting case is simplified for cases where the chatter is dominated by a single mode, the term Qp T A0Qp becomes into a scalar B0 and the gain is simplified as a scalar, as done in equation ( 18). In this case, the directional matrix would not have a decisive role.

    η g e B a K η i tc      1 1 sin Φ c i 0 p    -1 (34)
Nevertheless, a bidirectional actuator is recommended in order to improve properly the stability of a machine where the principal modes are defined in a plane. Moreover the proposed formulation can be inaccurate when interrupted cutting is performed, since zero order approximation is used. Therefore, a deeper study of mathematical formulation would be required for these cases.

EXPERIMENTAL SETUP

This section presents the machine and the actuator employed for simulations and the subsequent experimental verification. A prototype of a DANOBAT ram type machining center (Fig. 4) was selected, since its stability was previously studied in frequency domain with excellent results [START_REF] Munoa | [END_REF]].

Fig. 4: DANOBAT ram type machining center prototype.

The machine offers a variation of the ram overhang (Z) of 630 mm, height travel (Y) of 1100 mm and 1000 mm of X travel. In the present work, the ram overhang of 550 mm and the height of 760 mm remain constant in order to compare the results of different control laws. Tab 2. describes the cutting tool, the cutting conditions and the cutting coefficients used in all simulations and tests. Face milling operations with different spindle speeds were tested over the XY vertical plane, being X the feed direction.

Tab. 2: Description of cutting parameters. A complete modal analysis was performed in order to obtain the dynamic parameters of the machine, which are summarized in Tab. 1. It can be observed that the flexibility of the machine is mainly defined by XY directions, so the introduction of the control parameters is complicated. A single degree of freedom (DOF) small electromagnetic actuator [START_REF] Mancisidor | [END_REF]] has been located on the ram frontend (point q). Such device is characterized by its linearity (see Fig. 5) and provides a maximum force of 150 N. The ELMO TUBA servo amplifier has been used in order to convert the voltage commanded by the controller into current. A dSPACE DS-1005 board controller has been employed in order to implement the control laws, while the vibration signal is measured by a piezoelectric accelerometer.

Cutting tool

Reference

Natural frequency (ωn)

Damping ratio (ξ)

MECHATRONIC MODEL SIMULATIONS

The comparison of these four control strategies for orthogonal cutting process was reported in the literature [START_REF] Munoa | [END_REF], where DVF was showed as the best strategy at almost all spindle speeds. However, the complexity of milling process can distort the good behavior of the presented control laws, since than more than one mode appears and their directions do not match usually with the actuation direction [START_REF] Munoa | [END_REF]].

In general, single degree of freedom devices have been proposed in the literature. However, ram type machine tools present their flexibility on the two bending directions of the ram (see Tab. 1) and then, a biaxial inertial actuator can obtain better results [START_REF] Munoa | [END_REF]]. In this work, both a single direction actuator and a biaxial actuator have been simulated in the cutting process defined in Tab. 2. In the first case, the actuator described in the previous section has been considered. In this way, the mechatronic model has been employed in order to predict the stability improvement offered by each control law. The results shown in Fig. 6 present the DVF control strategy as the best one in almost all range of spindle speeds of the stability lobes. It is only overcome in speeds around the sweet spot (maximum stability zone), where DVF cannot improve the stability. However, as explained before, the machine flexibility is mainly defined by XY directions and more than one mode affects in the critical frequency range. Such issue does not affect negatively in strategies which are focused on the structural response modification. However, when DelPF or DelAF is employed, it affects drastically the results, since, at least, two directions should be controlled in order to reduce properly the regenerative effect, considering the directional matrix. In order to verify this statement, a second simulation where a biaxial actuator is simulated has been carried out.

The gain has been the same for both controllable directions and in this way, the real capability of each control law can be observed and compared in Fig. 7. It can be observed that DVF strategy remains the best in almost all zones of the stability lobes. However, the stability is improved considerably by the DelAF, principally in spindle speeds around the sweet spot and, hence, the maximum productivity (5 mm) is provided by this novel strategy. Therefore, it is proved that DelAF control law requires a biaxial actuator and the consideration of the directional matrix, since if not it cannot obtain a sizeable improvement.

EXPERIMENTAL RESULTS

Finally, the proposed control laws have been also compared experimentally in the DANOBAT machining centre prototype. The actuator proposed in section 3 is employed, so the fact that it is a single axis actuator has to be taken into account when the results are compared. Tab 3. shows numerically the results obtained for the selected position of the ram in the workspace (Y=760mm; Z=550mm). They can also be observed graphically in Fig. 8. The results are similar to those obtained by the mechatronic model with a single axis actuator. It can be clearly observed that almost all the other zones are clearly dominated by DVF and it is only overcome around the sweet spot. In such zone, DVF cannot improve the stability and then, DelAF achieve the best results. However, it should be commented that the DelAF results are poor, since a single axis actuator has been tested and hence, the directional matrix could not be considered properly. Therefore the theoretical results obtained for the biaxial actuator should be verified experimentally, in order to compare properly the control laws. The effect of the DVF control strategy in machine vibrations is clear if the actuator is switched on during the cutting process with chatter problems (see Fig. 9). In the one hand, the vibration acceleration signal is reduced more than a half when the actuator is working. On the other hand, the surface quality is considerably improved, which proves that the control suppress the chatter vibrations.

CONCLUSIONS

This work has presented a comparison of the most used model free feedback control strategies for suppressing chatter vibrations, in order to analyse the best control law for a certain milling process. A brief mathematical development of each control strategy is completed in the first part, where the complexity of the laws when they are used in milling operations is shown. The capital influence of the directional matrix consideration has been demonstrated when Delayed Position/Acceleration Feedback is employed in milling. If such matrix is not taken into account and the flexibility of the machine is defined by more than one mode, the delayed feedback strategies do not introduce the force correctly and in some cases, they can worsen the stability of the system. Then, the control laws have been tested theoretically by means of a validated mechatronic model, where all strategies can improve the stability of the system. The simulations of a single axis and a biaxial actuator have been carried out. The results show that DVF is the best control law for both actuators in minimum stability zone, although the results obtained by DelAF with biaxial actuator are also excellent, principally in spindle speeds close to the maximum stability zone. In this way, the maximum productivity is provided by the novel strategy, when a biaxial actuator is employed and the directional matrix is considered. Finally, these conclusions have been confirmed by experimental cutting tests with a single axis actuator, where the results obtained have been very similar to those obtained by the mechatronic model. It can be observed that DVF control strategy is the most suitable in most cases. However, the experimental tests with a biaxial actuator would be required in order to obtain the real difference between DVF and DelAF.

Spindle

Fig. 1 :

 1 Fig. 1: Chatter vibration scheme with feedback control.

Fig. 2 :

 2 Fig. 2: Points taken into account in the mechatronic model.

Fig. 3 :

 3 Fig. 3: Results obtained by the time domain model; a) Time domain signal of a chatter-free case; b) Spectrum of the speed of a chatter-free case; c) Time domain signal of a chatter case; d) Spectrum of the speed of a chatter case.

Fig. 4 :

 4 Fig. 4: DANOBAT ram type machining center prototype.

Fig. 6 :

 6 Fig. 6: Comparison of mechatronic model results of different control laws employed by a single axis actuator.

Fig. 7 :

 7 Fig. 7: Comparison of mechatronic model results of different control laws employed by a biaxial actuator.

Fig. 8 :

 8 Fig. 8: Comparison of experimental results of different control strategies.

  The minimum depth of cut obtained for each control strategy and spindle speed when Y=760mm and Z=550mm.

Fig. 9 :

 9 Fig. 9: Cutting tests (N=610rpm, ap=1.25mm) with DVF control strategy; a) Time domain acceleration signal of machine vibrations; b) surface finishing.

ACKNOWLEDGMENTS

This work has been partially supported by the PAINT ETORGAI (ER-2012/00019) project and REMEC project (IG-2013/0001136), both funded by the Basque Government.