Bernard Chalmond
email: bernard.chalmond@cmla.ens-cachan.fr

A macro-DAG structure based mixture model

Keywords: Mixture model, DAG structure, Bayesian network, EM algorithm

In the context of unsupervised classification of multidimensional data, we revisit the classical mixture model in the case where the dependencies among the random variables are described by a DAG structure. The structure is considered at two levels, the original DAG and its macro-representation. This twolevel representation is the main base of the proposed mixture model. To perform unsupervised classification, we propose a dedicated algorithm called EM-mDAG, which extends the classical EM algorithm. In the Gaussian case, we show that this algorithm can be efficiently implemented. The experiments reveal that this method favors the selection of a small number of classes.

Introduction

Let X be a random vector with values in R n for which we have a N -sample X = {x 1 , ..., x N } with n<N. Our goal is the clustering of X . This task is approached through a mixture model but with a particular constraint that makes the specificity of our contribution.

The dependency structure among the n components X j of X is subject to a structure represented by a DAG, in other words X is a Bayesian network. This structure induces a partition of X into M +1random vectors called macro-variables : X = ⊎ M m=0 X Jm ,w h e r eX Jm = (X j1 , ..., X jm ) when J m = {j 1 , ..., j m }.Fig. 1 depictes an example with M =3and J 0 = {1}, J 1 = {2, 3}, J 2 = {4, 5}, J 3 = {6, 7, 8}.

Each macro-variable X Jm is dependent on a hidden class variable C m with values in K m = {1, 2, ..., ν m }. Each occurrence in K m is the number of a class called elementary class. Therefore X is dependent on the hidden multi-class variable C =( C 1 , ..., C M ) whose values are in K = ⊗ M m=0 K m . Each (M +1)-tuple of K refers to a set of elementary classes called composite class. The (M +1)-tuples can be interpreted as pathways connecting the elementary classes through the macro-variables as it is illustrated in Table 1. The objective is to find the most probable pathways. We considerer the mixture model

pθ(x)= k∈K α k pθ k (x | k) ,
where the probability distribution pθ k (x | k) is that of the Bayesian network conditional on the composite class k and θk denotes the set of parameters defining this distribution. In this paper we describe this mixture model and we give a version of the EM algorithm, called EM-mDAG, for performing unsupervised classification. One of the main role of the EM-mDAG algorithm is to reveal probabilistic relationships among hidden elementary classes. Its implementation is done in the Gaussian case. Simulations illustrate the method and reveal a specific property. The EM-mDAG algorithm can select a small number of significant composite classes in K.

Models and Method

Basic knowledge

• Conventional mixture model for non supervised classification.

Let a random vector X =(X1 , ..., X j , ..., X n ) with values in R n . We assume that its probability distribution p φ (x) is a mixture of ν distributions {p θ k (x)} as follows :

p φ (x)= ν k=1 α k p θ k (x) with ν k=1 α k =1.
(2.1)

p θ k (x)
is defined by a parametric law of parameters θ k , as for instance the Gaussian law. The parameter set is denoted φ = {α, θ} where α = {α k } and θ = {θ k }. This mixture model can be interpreted in the context of unsupervised data classification. Let C be the hidden variable, which is an indicator variable of classes with values in {1, ..., ν}. Then, (2.1) is rewritten as

p φ (x)= ν k=1 P (C = k) p θ k (x | C = k). (2.
2)

The classification is to assign a class to every observation x, 1 .Whenφ is given, the MAP decision rule consists to choose the class

k(x) = arg max k P φ (C = k | x). (2.3)
Otherwise, things are more complicated because k(x) and φ have to be simultaneously estimated.

On the basis of maximum likelihood, the EM algorithm allows this estimation from a sample X = {x 1 , ..., x N } of X.

X4

J 1 =(X 2 ,X 3 ), X J 2 =( X 4 ,X 5
) and X J 3 =(X 6 ,X 7 ,X 8 ) ; the small circles depict the elementary classes.

The general formulation of the EM algorithm, which is also valid for our particular case, reads as follows. If φ(ℓ) is an estimation of φ, then an updated estimation is :

φ(ℓ + 1) = arg max φ Q(φ|φ(ℓ)) , Q(φ|φ(ℓ)) = IE C|X [log p φ (X , C) | φ(ℓ)] , (2.4)
where C = {C 1 , ..., C N } is a series of i.i.d. variables related to C. Q is an expected loglikelihood with respect to p φ(ℓ) (C |X ). The EM algorithm is an iterative procedure. From an initial estimate φ(0), it computes successively φ(0) → ... → φ(ℓ) → .... The marginal likelihood series {p φ(ℓ) (X ),ℓ =0, 1, ...} is non-decreasing.

• Bayesian network.

The previous classical formalism is the primal version for mixture modeling in the context of classification [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. The EM algorithm also applies to more complex situations such as those where the X i are not i.i.d. variables, but are dependent through hidden variables C i governed by a Markov chain [START_REF] Baum | A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains[END_REF] or a Markov random field [START_REF] Chalmond | An iterative Gibbsian technique for the reconstruction of m-ary images[END_REF]. In this article, we remain in the case where X is a sample from i.i.d. variables, but we consider a Markov structure for the dependence of the components X j . This Markovian structure is based on a DAG denoted G =( V, E). V = {1, ..., j, ..., n} denotes the variable numbers. The edges E ⊂ V × V are directed : (j ′ ,j) ∈ E is denoted j ′ → j.Theset j = {j ′ : j ′ → j} denotes the parents of the node j. The DAG structure has a fundamental property due to its acyclic nature : there is a numbering of the nodes such that j ⊂{1, 2, ..., j -1}. We assume that the nodes have been ordered in this way. With this property and that of Markov, we get the factorization p(x)= j p(x j | x j ).

(2.5)

The set B =( X,G,{p(x j | x j }) is called Bayesian network. When the distribution p(x) is non homogeneous, a mixture model as (2.2) can be considered in which p θ k (x | C = k) denotes a Bayesian network conditional on the hidden class C. This mixture model has been investigated in [START_REF] Thiesson | Learning mixtures of DAG models[END_REF] with a particular interest for DAG structure estimation.

Mixture model, composite class and Bayesian network

Composite class model

Let a partition of V composed of M +1macro-nodes : V = J 0 ⊎ J 1 ⊎ ... ⊎ J M , built from the DAG structure : J m is a macro-node if all its nodes have the same parents (In Fig. 1, M =3,and

J 0 = {1}, J 1 = {2, 3}, J 2 = {4, 5}, J 3 = {6, 7, 8}
). J 0 is the root of the tree and most often is a single node2 .LetJ 1 , ..., J M be the parents of J 1 , ..., J M , respectively. Given the definition of macro-nodes, each J m is composed of a single macro-node (In Fig. 1, J1 = J 0 , J2 = J 1 , J3 = J 1 ). The macro-nodes V = {J m } and their connexions E induced by {J m } defineane w directed acyclic graph G =(V, E) called macro-DAG.

Given a set of specifications {p(x Jm | x Jm } for B, a Bayesian network B =(X, G, {p(x Jm | x Jm }) can be defined for the macro-variables {X J m } M m=0 . The difference with B is essentially that B is a vectorial process whose factorization formula is written as

p(x)= m p(x Jm | x Jm ).
(2.6)

The factorization (2.6) assumes that the probability distribution is homogeneous, whereas it is not the case in our context. The distribution is depending on a hidden class variable C,which implies that p(x) is a mixture of distributions.

Firstly we assume that each macro-variable X J m is characterized by ν m classes, called elementary classes, whose parameters are denoted θ m = {θ m 1 , ..., θ m νm }. If we forget for a while the DAG structure, then each variable taken independently of the others, is defined by a mixture model for which (2.1) is rewritten as 1. The classification is to assign a composite class to each observation x. This involves selecting an elementary class k m for each macro-variable. An immediate solution would be to perform M independent classifications, based on (2.7) but this approach would have the disadvantage of not considering the DAG structure. Therefore we must address the classification as a whole.

p(x Jm )= νm k=1 P (C m = k) p θ m k (x Jm | C m = k). (2.7) Secondly, we consider the indicator variable of composite classes C =( C 1 , ..., C M ) with values in the set of M -tuples K = {k =( k 1 , ..., k M )} where k m ∈{ 1, ..., ν m }, as represented in Table
Considering the DAG structure, a composite class k is not only defined by the parameters θ k = {θ m km } M m=1 of its elementary classes, but also by the dependency parameters θk = { θm km } M m=1 that define the specifications of the Bayesian network X conditionally to C = k3 . These parameters are related to the parameters θ k . For each composite class, the factorization formula (2.6) based on the macro-DAG is written as

pθ k (x | k)=p(x J0 ) M m=1 pθm k (x Jm | x J m ,k m ,k m ) , (2.8)
where k m denotes the class number associated to x J m and appearing in k. In the notation pθm k , only the classes k m and k m of k are active. Finally the mixture model is written as

pθ(x)= k∈K α k pθ k (x | k) . (2.9)
Initially in (2.7) the definition of elementary classes has been made independently for each macro-variables. Now, the Markov dependence (2.8) introduces dependencies among these classes.

The parameter setting of the mixture model (2.9) differs from the classical mixture model (2.1).

Two M -tuples may have common components. For example, all components of (1,k 2 , ..., k M ) and (2,k 2 , ..., k M ) are identical, except the first. Thus, since two M -tuples may have common components, two components of the mixture may have common parameters 4 . In fact, there is one parameter setting per class, totaling m ν m settings, while there are |K| = m ν m composite classes.

EM-mDAG algorithm

The ultimate objective is to assign a composite class to every observation x :

x → k(x) = arg max k∈K P φ (C = k | x) .
Therefore, it is necessary to estimate φ =(α, θ). In an equivalent manner to (2.4), the estimation of φ is based on the log-likelihood by maximizing the Lagrangian function

L(α, θ)= N i=1 log k∈K α k pθ k (x i | k) + λ k∈K α k -1 , = N i=1 log k∈K α k M m=1 pθm k (x Jm i | x J m i ,k m ,k m ) + λ k∈K α k -1 , (2.10)
where λ denotes the Lagrangian parameter associated to the constraint k∈K α k =1 .A tt h e iteration ℓ of the EM algorithm, the re-estimation formula of α is written as in the classical case :

α k (ℓ +1)= 1 N i p φ(ℓ) (k | x i ) , (2.11)
where the a posteriori probability of the composite class k is defined by

p φ(ℓ) (k | x i )= α k (ℓ) pθ k (ℓ) (x i | k) p φ(ℓ) (x i ) = α k (ℓ) pθ k (ℓ) (x i | k) k∈K α k (ℓ)pθ k (ℓ) (x i | k)
.

(2.12)

As we said above, the peculiarity of this variant of the EM algorithm is the fact that a same parameter θm km can be present in several composite classes. In the classical case (2.1), the gradient of the Lagrangian function with respect to θ k concerns only p θ k while in (2.10), the gradient with respect to θm km relates to several components pθ k .

Proposition 1. The re-estimation formula of θ is given by the solution θ(ℓ +1)of the linear system N i=1 k=(k1,...,kM ):

km=τm p φ(ℓ) (k | x i ) ∂ ∂ θm τm log pθm τm (x Jm i | x J m i ,τ m ,k m ) θ= θ(ℓ+1) =0,
τ m =1, ..., ν m ,m=1, ..., M.

(2.13)

Proof. Let's focus on θm τm where τ m is a particular class number in {1, ..., ν m }.

∂L(α, θ) ∂ θm τm = N i=1 1 p φ (x i ) ⎡ ⎢ ⎢ ⎣ k=(k1,...,kM ): km=τm α k ∂ ∂ θm τm pθ k (x i | k) ⎤ ⎥ ⎥ ⎦ , = N i=1 1 p φ (x i ) ⎡ ⎢ ⎢ ⎣ k=(k1,...,kM ): km=τm α k pθ k (x i | k) pθ k (x i | k) ∂ ∂ θm τm pθ k (x i | k) ⎤ ⎥ ⎥ ⎦ , = N i=1 k:km=τm p φ (k | x i ) ∂ ∂ θm τm log pθ k (x i | k) .
(2.14)

Recalling the factorization formula (2.5), the gradient can be written as

∂L(α, θ) ∂ θm τm = N i=1 k:km=τm p φ (k | x i ) ∂ ∂ θm τm M a=1 log pθa k (x Ja i | x J a i ,k a ,k a ) , = N i=1 k:km=τm p φ (k | x i ) ∂ ∂ θm τm log pθm τm (x Jm i | x J m i ,τ m ,k m ) , (2.15)
which leads after a shortcut, to the system (2.13).

Gaussian case, linear dependency model and DAG

• Linear dependency model and DAG.

Under the Gaussian assumption, conditionally on the elementary classes, the law of the macrovariables are

X Jm|km =[X Jm | k m ] ∼N(µ m km , Γ m km ) , (2.16)
and with respect to the DAG, the transition laws among these variables are

[X Jm|km | x J m ,k m ]=[ X Jm | x J m ,k m ,k m ] ∼N(µ m,x km|k m , Γ m,x km|k m
) .

(2.17)

We assume the linear regression model

µ m,x km|k m = A m km,k m x J m + b m km,k m , Γ m,x km|k m =Γ m km|k m .
(2.18) Therefore, the respective parameter settings of (2.16)and(2.17) are respectively

θ m km = {µ m km , Γ m km } , θm km = {A m km,k m ,b m km,k m , Γ m km|k m } .
Note that the linear regression model (2.18) depends on the direction of the DAG. A m km,k m characterizes the regression of X Jm on x J m and not the reverse. Note also that the regression model is multidimensional in output since X Jm|km is a random vector in R |Jm| .

• Re-estimation formulas.

The first approach to obtain these formulas would be to use (2.13) taking into account the Gaussian log-density that for any observation x i can be written as :

log pθm km (x Jm i | x J m i ,k m ,k m )=cst - 1 2 | log Γ m km|k m | - 1 2 (x Jm i -µ m,xi km|k m ) ′ (Γ m km|k m ) -1 (x Jm i -µ m,xi km|k m
) .

Establish equation (2.13) requires to derivate with respect to all the components of θm km .F o r instance, the partial derivative with respect to A m τm,k m is written as

∂ A m τm,k m log pθm τm (x Jm i | x J m i ,τ m ,k m )=(Γ m τm,k m ) -1 x J m i (x Jm i -µ m,xi τm,k m ) ′ .
We could go on, but there is a more direct way to operate.

Proposition 2. Since (2.18) is a linear regression model, it is easier to use the classical formulas of the maximum likelihood estimation of this model. In our context, these formulas [5]a r e written as

b m τm,k m (ℓ +1)= IE(X Jm|τm ) -A m τm,k m (ℓ +1) IE( X J m |k m ) , A m τm,k m (ℓ +1)= Cov(X Jm|τm ,X J m |k m ) Var(X J m |k m ) -1
, where . denotes an empirical estimate.

However, the empirical estimate of the moments (expectations, covariance matrix and variancecovariance matrix) must be weighted by weights w derived from the DAG :

b m τm,k m (ℓ +1)= N i=1 k:km=τm w τm i,k (ℓ) x Jm i -A m τm,k m (ℓ +1) N i=1 k:km=τm w τm i,k (ℓ) x J m i , (2.19)
where the weights at iteration ℓ are

w τm i,k (ℓ)= p φ(ℓ) (k | x i ) N i=1 k:km=τm p φ(ℓ) (k | x i )
.

Similarly, by denoting μJm|τm = IE (X Jm|τm ),wehave

A m τm,k m (ℓ +1)= N i=1 k:km=τm w τm i,k (ℓ)(x J m i -μJ m |τ m )(x Jm i -μJm|τm ) ′ × N i=1 k:km=τm w τm i,k (ℓ)(x J m i -μJ m |τ m )(x J m i -μJ m |τ m ) ′ -1
.

(2.20)

Finally, we have also

Γ m τm,k m (ℓ +1)= N i=1 k:km=τm w τm i,k (ℓ) x Jm i -µ m τm,k m (ℓ +1) x Jm i -µ m τm,k m (ℓ +1) ′ .
(2.21)

Note that the programming of the re-estimation formulas (2.19, 2.20, 2.21) is relatively difficult because two data structures interfere : the dependency structure derived from the DAG, and the list structure of the composite classes (cf. (1) for instance).

• Back to the elementary class parameters.

For all x i , the estimated composite class k(x i )=( k 1 (x i ), ..., k M (x i )) has been computed. The user is also interested by the parameters θ k of the elementary classes { k m (x i )} M m=1 in order to interpret the leaves of the tree. The law (2.16) ignores the DAG, contrary to the law (2.17). However, in the Gaussian case, the parameters θ k are related to the parameters θk as follows [START_REF] Hastie | The Elements of Statistical Learning[END_REF]:

µ m,x km|k m = µ m km +Γ m km,k m (Γ m k m ) -1 (x J m -µ m k m ) , Γ m km|k m =Γ m km +Γ m km,k m (Γ m k m ) -1 Γ m k m ,km .
(

2.22)

To avoid the difficulty of solving the system with respect to µ m km and Γ m km , we consider

μm km = 1 N i x Jm i 1 km(xi)=km , Γm km = 1 N i (x Jm i -μm km )(x Jm i -μm km ) ′ 1 km(xi)=km .
• Initial solution.

The solution at the first step of the EM-mDAG algorithm is obtained by performing M independent classifications using the conventional EM algorithm. Therefore, for each macro-variable, we have ν m clusters in R |Jm| whose labels are { k m (x Jm i ),i=1 , ..., n}. From there, the initial solution θm km,k m (0) at iteration ℓ =0is computed using ordinary linear regression for every pair of clusters (k m ,k m ) for which there are observations : {i :

k m (x Jm i )=k m , k m (x J m i )=k m } = ∅.
Starting from this initial solution, the role of the EM-mDAG algorithm is to re-organize the clusters in order to extract from K a set of composite classes of high likelihood.

=3, ν 0 =1 ,ν 1 =2 , ν 2 =3 , ν 3 =4 .Herec =1.5.
µ k :

X 1 : 0 ,0 ,0 ,0 ,0

X 2 :
-c, -c, +c, +c, +c X 3 : +c, +c, -c, -c, -c

X 4 :
-c, -c, +c, +c, +c X 5 :

-c, -c, -c, +c, +c

X 6 :
-c, +c, -c, +c, +c X 7 :

-c, +c, +c, -c, -c X 8 :

-c, +c, -c, +c, +c

Experiments on simulated data

The random vector X of dimension n =8consists of M =3macro-variables with ν 1 =2 , ν 2 =3 , ν 3 =4 . Among the |K| =2 4potential composite classes, only 5 significant composite classes K 0 were considered. It implies that for k ∈K 0 no data has been generated, or in another words

α k =0, ∀ k ∈K 0 . (3.1)
Table 2 gives the expectation of these 5 composite classes and therefore the expectation of the elementary classes within the 3 macro-variables. Fig. 4 shows the observations of the macro-variables with their labels. This simulation was inspired by the cytometry data analysis domain (see Appendix) but with much more overlapping of the elementary classes. The first macro-variable X J1 =( X 2 ,X 3 ) shows two groups that it is possible to manually split, giving rise to two elementary classes denoted X 2 + and X 3 +. Each group is a mixture that the other two macro-variables help to identify. The macro-variable X J2 =( X 4 ,X 5 ) highlights the components of the group X 2 +, while the macro-variable X J3 =( X 6 ,X 7 ,X 8 ) highlights the components of the group X 3 +. However the overlapping of the mixture components in the groups X 2 + and X 3 + does not allow a partitionning of these groups as easy as for X J1 . Therefore we must address the classification as a whole.

• Data processing.

Fig. 5 shows the initial solution of the EM-mDAG algorithm at step ℓ =0. This solution results from M =3independent classifications by applying the classical EM algorithm on each macrovariable. The class labels k(x Jm i ) definedin(2.3) are gathered for making composite class labels by using a table similar to Table 1. This initial solution is unsatisfactory. The macro-variables X J2 and X J3 are strongly blurred by several small composite classes that are artefacts. The final solution of the EM-mDAG algorithm is shown in Fig. 6. The representation in terms of mixture components is close to the original in Fig. 4. The macro-variables X J2 and X J2 respectively highlight the components of the group X 2 + and X 3 +, despite the fifth class that is divided into two neighbour classes. Fig. 7 and Fig. 8 show respectively the classifications obtained with the usual EM algorithm successively performed on the basis of 24 classes and 5 classes. With 24 classes, the number of non-empty classes is large and therefore the classification is greatly erroneous. With 5 classes, the classification provided by the EM algorithm does not meet the specificity of the data. The macro-variable X J2 =(X 4 ,X 5 ) does not highlight the mixture components in the group X 2 +. This is a major problem that hinders significantly the interpretation of composite classes in terms of macro-variables.

• Property of the EM-mDAG algorithm.

The experiments show that the EM-mDAG algorithm has the property to keep a small number of α k different from zero when there is a limited number of significant composite classes K 0 ⊂ K. This selection ability is not so surprising. Firstly, X is not observable along k when k ∈K 0 , which means that its conditional distribution is not defined for these k. There exists at least one couple (k m ,k m ) in k whose observability of [X Jm |X J m ,k m ,k m ] is undefined. At every step ℓ of the algorithm, there are several couples (k m ,k m ) such that no observation x i is simultaneously present in the clusters k m and k m : {i :

k m (x Jm i )=k m , k m (x J m i )=k m } = ∅.
Secondly, the Markovian dependency introduced by the specifications p(x Jm i |x

J m i ,k m ,k m )
has for effect to reorganize the initial clustering while maintaining a well-contrasted partitioning. This is a wellknown property of the Markovian approach. These two remarks should be useful for establishing an analytical proof of the selection property.

Discussion

In this paper, we have presented a mixture model dedicated to the case where the dependencies among the components of the multidimensional random vector are governed by a DAG structure. The mixture model takes advantage of a two-level structure, which is composed by the DAG itself and its macro-representation. A dedicated EM algorithm has been efficiently implemented for the Gaussian case. The experiments show that this algorithm is able to select a small number of composite classes. This selection ability is important because it allows to circumvent the difficulty of choosing the exact number of elementary classes for each macro-variable. In fact, one of the main role of the EM-mDAG algorithm is to reveal significant relationships among the hidden elementary classes, some of them becoming empty during the procedure. The analysis, which is based on a dependency tree as illustrated in Fig. 1-a, is usually accomplished by sequential manual partitioning (called "Gating") of the sample X from the top to the bottom of the tree. Rather than watching simultaneously the n dimensions, that is to say the cloud X in the space R n , the biologist works in subspaces of smaller dimensions, 1, 2 or 3, according to associations of variables X j , here called macro-variables, as shown in Fig. 1-b.

At the top of the tree, only one coordinate of the cloud X is analyzed. This is the variable X 1 corresponding to high values CD45+ of the biological variable CD45. In this example, to simplify, the tree height was reduced by starting the tree with X 1 = CD45+ instead of (CD45-,CD45+). To determine the two groups CD45et CD45+, a threshold τ CD45 is manually selected for separating the small and large values of CD45. Conditionally on the elementary class X 1 = CD45+ , the procedure continues along the tree structure, as follows.

Three elementary classes are extracted from the 2-D distribution of the sample {(x 2 i ,x

i )} N i=1 3
and denoted (X 2 +), (X 2 -,X 3 -), (X 3 +) as illustrated in Fig. 2 and Fig. 3-a, On each group, this operation is repeated on the following macro-variables in dimension 2 for (X 4 ,X 5 ) conditionally on (X 2 +) as illustrated in Fig. 3-b, and in dimension 3 for (X 6 ,X 7 ,X 8 ) conditionally on (X 2 -,X 3 -). This conditional and sequential procedure can be represented by a DAG and then modeled by a Bayesian network. The main advantage of using the EM-mDAG is its ability to global classification while keeping the biological dependency structure, which is necessary for identifying the cell types. for a sample of size N = 1000. The graphic "macro-variable (X2,X3)" depicts the coordinates {(x 2 i ,x 3 i )} N i=1 ,t h e "macro-variable (X4,X5)" depicts {(x 4 i ,x 5 i )} N i=1 , and the "macro-variable (X6,X7,X8)" depicts {(x 6 i ,x 7 i ,x 8 i )} N i=1 . The histogram gives the empirical distribution of the composite classes. 5) and (X 6 ,X 7 ,X 8 ) respectively highlight the components of the group X 2 + and X 3 + of the macro-variable (X 2 ,X 3 ), despite the fifth class that is divided into two neighbour classes, numbered k =4and k =5in the histogram. The macro-variable (X 4 ,X 5 ) does not highlight the mixture components of the group X 2 + of (X 2 ,X 3 ).

Figure 1 .

 1 Figure 1. Two-level structure. (a) DAG structure. (b) Macro-DAG structure with its macro-variables XJ 1 =(X 2 ,X 3 ), X J 2 =( X 4 ,X 5) and X J 3 =(X 6 ,X 7 ,X 8 ) ; the small circles depict the elementary classes.

Figure 2 .

 2 Figure 2. A partial decision tree with biological classes on leaves, (thanks to Xiaoyi Chen, Institut Pasteur).

Figure 3 .Figure 4 .

 34 Figure 3. Two steps of the sequential procedure leading to the biological classes "Th cells" and "Tc cells" in Fig.2, (thanks to Xiaoyi Chen, Institut Pasteur). From the distribution of the sample {(x 2 i ,x 3 i )} N i=1 shown in (a), 3 elementary classes (X 2 +), (X 2 -,X 3 -), (X 3 +) are manually extracted. (b) shows the distribution of the sample {(x 4 i ,x 5 i )} limited to the records i coming from the class (X 2 +). Conditionally to (X 2 +), 2 new elementary classes (X 4 +) and (X 5 +) are manually extracted.

Figure 5 .

 5 Figure 5. Initial solution of the EM-mDAG at step ℓ =0 . M =3independent classifications was achieved by applying the classical EM algorithm on each macro-variable. Compared with the ground true in Fig.4, this representation is strongly blurred.

Figure 6 .

 6 Figure 6. EM-mDAG based clustering at step ℓ =2 0 .AsinFig.4, the macro-variables (X 4 ,X5) and (X 6 ,X 7 ,X 8 ) respectively highlight the components of the group X 2 + and X 3 + of the macro-variable (X 2 ,X 3 ), despite the fifth class that is divided into two neighbour classes, numbered k =4and k =5in the histogram.

Figure 7 .

 7 Figure 7. Standard EM algorithm for 24 classes. The number of non-empty classes is large and therefore the classification is greatly erroneous.

Figure 8 .

 8 Figure8. Standard EM algorithm for 5 classes. The macro-variable (X 4 ,X 5 ) does not highlight the mixture components of the group X 2 + of (X 2 ,X 3 ).

Table 1 .

 1 Composite class numbering for M =3and ν 0 =1,ν 1 = ν 2 =2,ν 3 =4. This table gives the exhaustive list of the 16 composite classes K, where each column is an (M +1)-tuple (1, k) with k ∈K.

	m=0 :	1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
	m=1 :	1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2
	m=2 :	1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2
	m=3 :	1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4

Table 2 .

 2 Expectations µ k of the 5 composite classes K 0 for data simulation with M

A class is defined by its number and its parameters. More often, we confound "class" and "class number".

J 0 has only one class and therefore ν 0 =1 .

In this paper, the notation . is reserved to parameters associated to the DAG dependencies.

In the classical case, several components can also be concerned by a same parameter, for instance the same variance in the Gaussian case, but it is not a constraint contrary to our mixture model where many parameters are necessarily shared.

Appendix : a case study

This section presents a domain in which our method should be useful, as it is currently being undertaken by Xiaoyi Chen at Institut Pasteur (Systems Biology team). A N -sample X of tens and even hundreds of thousands of cells is observed by flow cytometry. For each cell i =1, ..., N , the instrument provides a measurement vector x i of dimension n. This sample is a mixture of several cell populations. The goal is to group these measurements so that each class corresponds to a well-identified cell type [START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF].


Moussab Bennehar

Ahmed Chemori
email: chemori@lirmm.fr

Sbastien Krut
email: krut@lirmm.fr

Franois Pierrot
email: pierrot@lirmm.fr

Continuous Closed Form Trajectories Generation and Control of Redundantly Actuated Parallel Kinematic Manipulators

This paper deals with continuous closed form trajectories generation and control of redundantly actuated parallel kinematic manipulators. Continuous closed form trajectories are usually given by means of analytical functions of time and show inherent discontinuities in velocity and/or acceleration profiles. In order to solve this problem, we propose an approach consisting of modifying the time function without changing the overall closed shape of the trajectory in the operational space. To deal with the actuation redundancy, an extended version of the PD controller with computed feedforward is proposed. The resulting control torques are projected using a regularization matrix in order to remove internal forces that may be dangerous for the manipulator. The proposed trajectory generator as well as the controller are implemented in real-time on Dual-V; a redundantly actuated parallel manipulator developed in our laboratory.

I. INTRODUCTION

Parallel Kinematic Manipulators (PKMs) are mostly known for their better performances over their serial counterparts [START_REF] Merlet | Parallel Robots[END_REF]. Indeed, in contrast with the serial ones, the actuators of PKMs are located on the fixed base and, hence, result in a very lightweight moving mass [START_REF] Corbel | Towards 100G with PKM. Is Actuation Redundancy a Good Solution for Pick and Place?[END_REF]. Furthermore, the closed kinematic chains yield to a stiffer and more accurate structure. Nevertheless, PKMs exhibits many drawbacks that moderate their expansion in industrial applications. The abundance of singularities and the limited workspace are the most noteworthy ones. While the latter is a matter of optimal design, the former can be solved by actuation redundancy. Hence, Redundantly Actuated PKMs (RA-PKMs) are increasingly attracting more interest of researchers in the last two decades.

Actuation redundancy significantly improves the performance of PKMs by homogenizing the physical properties throughout the workspace, allowing the manipulator to achieve higher accelerations. However, despite the aforementioned qualities, PKMs potential is not yet completely explored and is still an open research field. From one hand, mechanical design, identification and optimization can be further investigated in order to increase their acceptance in industrial applications. From the other hand, the close-loop constraints and the inherent nonlinearities give rise to more challenging tasks in terms of trajectory generation and control. This research was supported by the French National Research Agency, within the project ANR-ARROW Thanks to their high dynamic performances, PKMs are typically used for industrial high-speed applications such as pick-and-place in food industry [START_REF] Bonev | Delta parallel robot-the story of success[END_REF]. Hence, most researches in terms of trajectory generation focused only on classical Point-to-Point (PtP) trajectories [START_REF] Macfarlane | Jerk-bounded manipulator trajectory planning: design for real-time applications[END_REF], [START_REF] Liu | Time-Optimal Trajectory Generation of a Fast-Motion Planar Parallel Manipulator[END_REF], [START_REF] Su | Trajectory planning for constrained parallel manipulators[END_REF]. Nevertheless, the inherent stiffness and accuracy of PKMs awarded them to be used in more complex industrial tasks such as laser cutting [START_REF] Bruzzone | Modeling and design of a parallel robot for laser-cutting applications[END_REF] or machining [START_REF] Chanal | Etude de l'emploi des machines outils structure parallle en usinage[END_REF]. In this case, Continuous Closed Shape (CCS) trajectories draw more attention than traditional PtP ones. This class of trajectories, however, has not been sufficiently investigated in the literature.

It is conventional that PKMs share many physical properties with serial robots. Consequently, most of the control literature developed for serial manipulators has been straightforwardly implemented on PKMs [START_REF] Paccot | A Review on the Dynamic Control of Parallel Kinematic Machines: Theory and Experiments[END_REF], [START_REF] He | A survey on control of parallel manipulator[END_REF]. However, in the case of RA-PKMs, a particular property characterizing this class of manipulators needs a special attention. The nonuniqueness of the inverse dynamics solution yields to control forces with no effects on the moving platform of the robot [START_REF] Miiller | A projection method for the elimination of contradicting control forces in redundantly actuated PKM[END_REF]. These control forces produce undesired internal pre-stress in the kinematic chains that can harm the manipulator and cause a loss of energy. Consequently, the control scheme should take into consideration such a phenomenon.

In this work, we address the two aforementioned problems, namely; the CCS trajectories generation and the reduction of internal forces for RA-PKMs. We propose a trajectory generator that takes into consideration continuity constraints, both in task and joint spaces. As regards to the control solution, a joint space PD with computed feedforward controller [START_REF] Reyes | Experimental evaluation of model-based controllers on a direct-drive robot arm[END_REF] is proposed. This control scheme has the advantage of being lightweight in terms of computation time while compensating some of the nonlinearities. In order to deal with the internal forces issue, the proposed controller is extended with a projection of the generated control torques with the regularization. The proposed approach is then validated in real-time on a 3 dof RA-PKM, namely, the Dual-V. This paper is organized as follows. In section 2, the dynamic model of the Dual-V robot is presented. Section 3 is devoted to the proposed trajectory generator. The proposed control solution is detailed in Section 4. Experimental results are presented in Section 5. Finally, conclusions and perspectives are drawn in Section 6.

SSD'14 1569847595

II. DYNAMIC MODELING OF THE DUAL-V ROBOT Dual-V is a 3 dof planar redundantly actuated parallel manipulator belonging to the 4-RRR family. The arrangement of the four chains allows three independent movements for its end-effector (a traveling plate): two translations in the plane and one rotation about the z-axis. Hence, the Cartesian coordinates the traveling plate are described by the vector X = [x, y, θ ] T . For the dynamic modeling of the robot, the approach developed in [START_REF] Corbel | Towards 100G with PKM. Is Actuation Redundancy a Good Solution for Pick and Place?[END_REF] has been extended to take into account the rotational inertia of the forearms.

The required torque to move the robot's mechanical structure can be decomposed into three sub-torques, namely; τ 1 , τ 2 , τ 3 ∈ R 4 . Each sub-torque is responsible for the movement of a part of the structure. The torque τ 1 is the required torque to move the traveling plate and a part of the mass of the forearms, it can be calculated using the equations of power of the actuators as follows:

τ 1 = J T * m M I Ẍ (1)
where M I ∈ R 3×3 is the mass matrix of the concerned part of the links to be moved in Cartesian space, J T m * = J m (J T m J m ) -1 ∈ R 4×3 is the pseudo-inverse of the transpose of the inverse Jacobian matrix J m and Ẍ ∈ R 3 is the acceleration vector of the traveling plate. The component τ 2 is the required torque to move the arms, the counter-masses and the remaining part of the couplers, it can be expressed as:

τ 2 = M II q = M II ( Jm Ẋ + J m Ẍ) (2)
where M II ∈ R 4×4 is the mass matrix regrouping dynamic parameters of the involved parts of the structure in joint space, Ẋ ∈ R 4 is the traveling plate velocity vector and Jm is the time derivative of the inverse Jacobian matrix. Up to now, only the inertia of the two mass model [START_REF] Van Der Wijk | Design and exper-imental evaluation of a dynamically balanced redundant planar 4-RRR parallel manipulator[END_REF] of the couplers is considered. Though this is a righteous hypotheses when light material is used, it fails when the links are made with heavy material such as steal or aluminum. This is the case of the Dual-V robot and, hence, an additional term should be added. The torque term τ 3 is the additional torque component to compensate for the real inertia of the couplers and that of the equivalent mass model. Due to the number of pages limitation, the details of this additional term are not provided herein and the reader is referred to [START_REF] Van Der Wijk | Design and exper-imental evaluation of a dynamically balanced redundant planar 4-RRR parallel manipulator[END_REF] for a full description of the Dual-V robot dynamic modeling. It is worth noting though that τ 3 (X, Ẋ, Ẍ) is a highly nonlinear function with respect to its arguments.

Finally, the full dynamic model of the Dual-v is written by regrouping the three sub-torques as follows:

τ = J T * m M I Ẍ + M II ( Jm Ẋ + J m Ẍ) + τ 3 (X, Ẋ, Ẍ) (3)
The CAD model of the Dual-V with parameters definition is depicted in Fig. 1 and its geometric and dynamic parameters are summarized in Table I. Further details about the Dual-V robot are given in [START_REF] Wu | A comparison study on the dynamics of planar 3-DOF 4-RRR, 3-RRR and 2-RRR parallel manipulators[END_REF] and [START_REF] Wijk | Generic Method for Deriving the General Shaking Force Balance Conditions of Parallel Manipulators with Application to a Redundant Planar 4-RRR Parallel Manipulator[END_REF]. CCS trajectories are more complex than traditional PtP ones [START_REF] Liu | Time-Optimal Trajectory Generation of a Fast-Motion Planar Parallel Manipulator[END_REF], [START_REF] Macfarlane | Jerk-bounded manipulator trajectory planning: design for real-time applications[END_REF] (e.g pick-and-place trajectories) and many of them are generated by means of analytical functions parametrized by the time variable t. Therefore, generating trajectories with continuity constraints on velocities and accelerations may not be possible using their original analytical equations. Continuity of trajectories is a crucial constraint for the robot drives since discontinuities in the trajectories can generate discontinuous control torques and thus, lead to undesired behavior of the mechanical structure such as vibrations or instabilities.

In this paper we are interested in the class of trajectories described by a sum of weighted sin and cosine functions. This class of trajectories has been chosen because it covers a large amount of simple and complex geometrical shapes (e.g. circles, ellipses, deltoids, . . . ). The general analytical form of the trajectories in question is given by:

x(t) = n ∑ i=1 a i cos α i 2π n 1i T t + b i sin β i 2π n 2i T t (4)
where T ∈ R + is the trajectory duration and a i , α i , b i , β i , n ji ∈ R; i = 1, . . . , n; j = 1, 2 are the parameters defining the shape of the trajectory.

In the sequel, we will first give an illustrative example of direct application of (4). Then, its major drawback is highlighted and a solution to overcome it is proposed.

A. Classical method

Let x d (t) be the desired Cartesian position defined by ( 4) to be tracked by the traveling plate of the robot. These desired trajectories are fully included in the workspace of the robot away from singularities. The desired velocity and acceleration trajectories for this motion are obtained by differentiating x d (t) with respect to time, thus:

v d (t) = - 2π T n ∑ i=1 a i α i n 1i sin α i -1 2π n 1i T t - b i β i n 2i cos β i -1 2π n 2i T t (5) a d (t) = - 2π T 2 n ∑ i=1 a i α i (α i -1)n 2 1i cos α i -2 2π n 1i T t + b i β i (β i -1)n 2 2i sin β i -2 2π n 2i T t (6)
Without loss of generality, assume that the robot starts and finishes its movement with a null velocity and a null acceleration (i.e. v(0) = v(T ) = 0 and a(0) = a(T ) = 0). Using the original analytical functions ( 5) and ( 6), however, may lead to non-zero initial velocity and/or acceleration. Indeed, if we replace t = 0 or t = T in ( 5) and [START_REF] Miiller | A projection method for the elimination of contradicting control forces in redundantly actuated PKM[END_REF] we will obtain values dependent on initial and final conditions and thus, does not necessarily equal zero. For instance, the obtained initial conditions are given by:

v d (0) = - 2π T n ∑ i=1 b i β i n 2i = f (b i , β i , n 2i ) (7)
a d (0) = - 2π T 2 n ∑ i=1 a i α i (α i -1)n 2 1i = f (a i , α i , n 1i ) (8)
therefore, the resulting trajectories can be velocity and/or acceleration discontinuous leading to discontinuous control torques and therefore big tracking errors or even instability.

To further illustrate this inconvenience let's consider a simple circular trajectory with a radius r and a center whose coordinates are (x c , y c ) = (0, 0). The resulting trajectory needs to be C 2 continuous (in velocities and accelerations) both in operational and joint spaces in order to be appropriately tracked by the traveling plate. The analytical equations of positions, velocities and accelerations are given by:

x d (t) = r cos 2π T (t -t 0 ) y d (t) = r sin 2π T (t -t 0 ) (9)
ẋd (t) = -r 2π T sin 2π T (t -t 0 ) ẏd (t) = r 2π T cos 2π T (t -t 0 ) (10)
ẍd (t) = -r( 2π T ) 2 cos 2π T (t -t 0 ) ÿd (t) = -r( 2π T ) 2 sin 2π T (t -t 0 ) (11)
These trajectories satisfy equations ( 4), ( 5) and ( 6) being only shifted in time by t 0 . Fig. 2 illustrates the plots of the obtained trajectories using ( 9), ( 10) and (11) (t 0 = 0.5s, T = 1s). One can notice the presence of discontinuities on both velocity and acceleration profiles that have to be removed. In what follows, we propose a novel method to overcome this drawback and generate C 2 continuous reference trajectories. In order to overcome the discontinuity problem in the velocity and acceleration trajectories, we propose to revisit the general form of the desired operational trajectories. Consider again the previous illustrative example of a circular trajectory. If we check the previous equations of the trajectory, it can be seen that they can be rewritten as follows:

x d (t) = r cos(λ (t)) y d (t) = r sin(λ (t)) (12)
With λ (t) = 2π T (t -t 0 ) which is a simple affine function of time t and when implemented, does not satisfy any continuity constraint on the velocity and acceleration trajectories. The corresponding velocities and accelerations are obtained by differentiating [START_REF] Bruzzone | Modeling and design of a parallel robot for laser-cutting applications[END_REF] with respect to time, which gives:

ẋd (t) = -r λ (t) sin(λ (t)) ẏd (t) = r λ (t) cos(λ (t)) (13)
   ẍd (t) = -r λ (t)sin(λ (t)) + λ 2 (t)cos(λ (t)) ÿd (t) = r λ (t)cos(λ (t)) -λ 2 (t)sin(λ (t)) (14)
From ( 13) and ( 14) one can notice that if λ (t) is chosen such that it satisfies the following boundary constraints:

λ (t 0 ) = 0; λ (t 0 ) = 0; λ (t 0 ) = 0 λ (t 0 + T ) = 2π; λ (t 0 + T ) = 0; λ (t 0 + T ) = 0 ( 15
)
one can ensure the continuity on the velocity and acceleration trajectories. One possible solution would then be to choose λ (t) with a trapezoidal profile on velocity. Indeed, this choice allows to specify boundary conditions and hence, the constraints ( 15) can be satisfied. A comparison between trapezoidal and linear velocity profile is depicted in Fig. 3. The analytical equations of λ (t) are now expressed as:

       λ (t) = λ f -λ 0 2τ 3 (T -τ) [-t 4 + 2τt 3 ] + λ 0 ; t 0 ≤ t ≤ τ λ (t) = λ f -λ 0 (T -τ) [t -τ 2 ] + λ 0 ; τ ≤ t ≤ T -τ λ (t) = λ f -λ 0 2τ 3 (T -τ) [(t -T ) 4 + 2τ(t -T ) 3 ] + λ f ; else (16)
where τ is the duration of initial acceleration and final deceleration phases that has to be chosen according to the robot dynamic capabilities. λ 0 and λ f are the initial and final values respectively for λ (t). For instance, in the case of the above circular trajectory example, λ 0 = 0 and λ f = 2π.

The obtained trajectories using the proposed strategy are now continuous in velocity and acceleration and result in the same trajectory in operational space as the original one (cf. Fig. 4). Therefore, they can be safely used in the control scheme without risks on the actuators.

IV. PROPOSED CONTROL SOLUTION

It has been shown that PKMs exhibit many similarities with their serial counterparts [START_REF] Cheng | Dynamics and Control of Redundantly Actuated Parallel Manipulators[END_REF]. Consequently, many control strategies that were mainly developed for serial manipulators have been successfully applied on PKMs. However, in the particular case of RA-PKMs, we notice the existence of antagonistic control torques not affecting the motion of the manipulator. Thus, using conventional control strategies will certainly give occurrence to internal forces (incompatible with the robot's kinematics) [START_REF] Miiller | A projection method for the elimination of contradicting control forces in redundantly actuated PKM[END_REF] and hence, create internal prestress in the mechanism. The antagonistic forces can cause a multitude of undesirable phenomena such as the loss of energy, instability and vibrations. Consequently, the control architecture has to take in consideration this peculiarity.

A. Joint space PD controller

Consider a RA-PKM with m degrees of freedom and n > m actuators and let the joints position vector be denoted by q ∈ R n . The PD control is a decentralized uncoordinated control strategy relying on the measured error between the desired and the actual position. Let e(t) be the joints error vector, i.e : e(t) = q d (t)q(t), where q d (t) denotes the desired joint trajectory. Then, the PD control action is expressed by:

τ PD = K p e(t) + K d ė(t) (17)
where K p and K d are positive definite feedback gain matrices usually chosen diagonal. If the actuators are all identical, which is usually the case for PKMs, the same feedback gains could be selected for each axis, namely; K p = k p I n×n and K d = k d I n×n , where I n×n is an identity matrix and k p and k d are positive scalars denoting the feedback gains that should be tuned to achieve satisfactory closed-loop performance. Good accuracy on high accelerations is mostly required in parallel robots.

B. Desired compensation by nonlinear feedforward

If the inverse dynamics of the robots are well known, it would be interesting to exploit this knowledge to further improve the tracking performance. Indeed, the tracking errors can be significantly improved by compensating the inherent nonlinearities. One possible way would be the addition of a feedforward term. This strategy enables in partially compensating the nonlinear dynamics in the feedforward way, i.e. only desired values of positions, velocities and accelerations are fed to the inverse dynamics model. For the case of the Dual-V robot, the inverse dynamics are given by (3), hence, the feedforward torque is computed as:

τ f f = J T * m M I Ẍd + M II ( Jm Ẋd + J m Ẍd ) + τ 3 (X d , Ẋd , Ẍd ) (18)
where the subscript d refers to the desired quantities.

C. Projection method to reduce internal forces

For RA-PKM, the control torques may contain antagonistic forces that does not create any motion as they are in the null-space of J T m and, thus, need to be removed. This can be achieved by projecting the input torques onto the range of J T m by the following projector [START_REF] Miiller | A projection method for the elimination of contradicting control forces in redundantly actuated PKM[END_REF]:

X d q d + - e K p + + Inverse Kinematics X d , Ẋd , Ẍd J T m * J T m Inverse Dynamics + + RA-PKM q τ * τ τ PD τ f f
where N J T m = I -J T * m J T m is the projector to the null-space of J T m . Hence, the internal antagonistic forces are eliminated by projecting the control torques onto the range space of J T m . This is achieved by using the projection matrix R J T m as:

τ * = R J T m (τ PD + τ f f ) (20)
The block diagram of the proposed control scheme is depicted in Fig 5.

V. REAL-TIME EXPERIMENTAL RESULTS

A. The Dual-V robot and its hardware configuration

The Dual-V links are all made with Aluminum and mounted on four direct drive actuators provided by ETEL Motion Technology. The actuators are fixed on an Aluminum base and can supply up to 127 N.m torques. We use Simulink and Real-Time Workshop (both from Mathworks Inc.) to implement and build the control algorithm. The generated code is then uploaded to the target PC; an industrial computer cadenced at 10 kHz and running xPC Target in real-time. The experimental setup is depicted in Fig. 6.

B. Real-time experimental results

The proposed trajectory generator as well as the proposed controller presented in the previous section were implemented on the Dual-V robot. For comparison purpose, we also implement the classical approach using the original unmodified equations. In order to further investigate the benefits of the proposed approach, four different geometric shapes are implemented, namely; a circle, an ellipse, a tear-like and a deltoide (cf. Fig. 7a). The traveling plate of the manipulator has to perform every geometric trajectory in T = 0.5 s. Between the end point of a shape trajectory and the starting point of the next one, the end-effector has to perform a PtP trajectory defined by a 5 th order polynomial function.

The controller feedback gains have been experimentally tuned using the trial-and-error technique. Usually, industrial robots are not equipped with velocity sensors, thus velocities have to be estimated. In our case, a carefully designed low-pass filter is used to estimate the velocities.

Figures 7b and7c depict the Cartesian tracking performance of the end-effector on both x and y axes. We can clearly see that in the case of original trajectories, the discontinuities cause a poor tracking performance. This is more noticeable at the beginning and the end of the circle and the ellipse. If the robot was equipped with accelerometers, one could notice the high frequency vibrations. However, if we check the obtained performance using the proposed approach, we can clearly see the improvements in terms of tracking errors at the end and starting points. It can be clearly seen that the proposed trajectories are accurately tracked thanks to the removed discontinuities. The improvements are more noticeable on the starting and end points of the circle and ellipse trajectories (t = 5, 5.5, 7.5 and 8s). These results were expected since the main purpose of this paper is to present a technique to remove these discontinuities from the trajectories generated by means of classical analytical functions. These discontinuities are a major source of tracking loss.

The control inputs generated when using the proposed method as well as those generated by classical one are depicted in Fig. 8. As expected, the controller generates very high torques when it faces discontinuities on the reference trajectories, because the drives need to generate high torques in order step from zero to a certain velocity in a very short time, which is practically very delicate.

VI. CONCLUSIONS AND FUTURE WORK

In this paper we investigated the continuous closed shape trajectories generation and the reduction of internal forces by means of control for RA-PKMs. These trajectories are defined using analytical functions that show inherent discontinuities. For this reason, we proposed in this work an approach to overcome this inconvenience for a specific class of trajectories by modifying the time function without changing the overall closed shape of the trajectory in operational space. Furthermore, to overcome the phenomenon of internal forces, a PD with computed feedforward controller has been implemented. The resulting control torques have been projected using the Jacobian matrix to reduce internal forces that can harm the mechanical structure. The latter helps in significantly reducing the energy consumption and the vibrations on the robot's traveling plate. To validate the proposed approaches, realtime experiments have been conducted on a real 3 dof planar redundantly actuated parallel manipulator developed in our laboratory. The proposed trajectory generator exhibited its superiority over the original trajectory generator that utilizes original analytical equations.

Fig. 1 :

 1 Fig. 1: The CAD of the Dual-V with parameters definition

Fig. 2 :

 2 Fig. 2: View of the generated Cartesian trajectories using the classical method

Fig. 3 :

 3 Fig.3: Evaluation of λ (t), λ (t) and λ (t): original (dashed line) and modified (solid line)

Fig. 4 :

 4 Fig. 4: Obtained C 2 continuous operational trajectories

Fig. 5 :

 5 Fig. 5: Block diagram of the proposed control scheme

Fig. 7 :

 7 Fig. 7: Reference trajectory and evaluation of the Cartesian tracking errors versus time

Fig. 8 :

 8 Fig. 8: Evolution of control torques: (top) original discontinuous trajectories and (bottom) proposed continuous trajectories

TABLE I :

 I Dual-V geometric and dynamic parameters