Pierre Druilhet
email: druilhet@ensai.fr

Walter Tinsson
email: walter.tinsson@univ-pau.fr

Optimal repeated measurement designs for a model with partial interactions

Keywords: AMS 2000 subject classifications. Primary 62K05, 62K10 Crossover design, optimal design, universal optimality, total effects

 methods to find optimal designs when the parameters of interest are total effects, which corresponds to the use of a single treatment, once selected. We give a lower bound of Φ p criteria and propose efficient designs of small sizes.

Introduction

In traditionnal crossover designs, the statistical model consider that each treatment has a carryover effect, which, for parsimony, is assumed not interacting with direct treatment effects. To enrich this model, Hedayat and [START_REF] Afsarinejad | Repeated measurements designs for a model with self and simple mixed carryover effects[END_REF] propose a model with partial interactions: the carryover effect of a treatment is different if the treatment is preceded by itself.

For such model, [START_REF] Kunert | Optimal crossover designs in a model with self and mixed carryover effects[END_REF] find that universally optimal designs for direct treatment effects have no consecutive pairs of identical treatments. This feature is mainly due to the fact that partial interactions between carryover and direct treatment effects are considered as nuisance parameter. However, the aim of an experiment is to select a single treatment which will be used alone and then will be preceded by itself. In that case, the parameter of interest, called total effects, is the sum of direct treatments and self carryover effects. [START_REF] Bailey | Optimality of neighbor-balanced designs for total effects[END_REF] give a review of situations where total effects are considered in the literature. They show that neighbour balanced designs having no treatment preceded by itself are efficient when the number of periods is small. However, for models with partial interactions, such designs are inefficient because total effects are not estimable.

In this paper, we propose a general method to construct optimal circular crossover designs for total effects under models with partial interactions.

Upper bounds of the information matrix are obtained by a new method that generalizes the approaches of [START_REF] Kushner | Optimal repeated measurements designs: the linear optimality equations[END_REF], [START_REF] Kunert | On the determination of optimal designs for an interference model[END_REF] and [START_REF] Bailey | Optimality of neighbor-balanced designs for total effects[END_REF]. Then we show that, in most cases, optimal designs are generated by one sequence or a mixture of two sequences that belongs to a set of k-1 possible sequences, where k is the number of periods.

Finally, we give some examples of optimal designs for several values of k. We also propose efficient designs generated by only one sequence of treatments.

To simplify our presentation, we first construct optimal design for a model without period effects. Then we obtain optimal designs for a model with period effects by generalizing [START_REF] Kunert | Optimal design and refinement of the linear model with applications to repeated measurements designs[END_REF] results on orthogonal effects.

The designs and the models

Let d(i, j) ∈ {1, ..., t} be the treatment assigned to subject i in period j, 1 ≤ i ≤ b, 1 ≤ j ≤ k. We assume that the response y ij follows the model:

y ij = β i + τ d(i,j) + λ d(i,j-1) + χ d(i,j-1)d(i,j) + ε ij (1)
where β i is the effect of subject i, τ i is the effect of treatment i, λ i is the general carryover effect of treatment i and χ ij is the additional specific carryover effect when treatment i is followed by itself (χ ij = 0 if i = j). The errors ε ij are assumed i.i.d. with expectation 0 and variance σ 2 . The model with period effects is considered in Section 5. In vector notation, we have:

Y = B β + T d τ + L d λ + S d χ + ε (2)
where B, T d , L d and S d are incidence matrices of subjects, direct treatments, carryover and specific self-carryover effects. Note that the parametrization, although equivalent, is slightly different from that of Hedayat and [START_REF] Afsarinejad | Repeated measurements designs for a model with self and simple mixed carryover effects[END_REF] or [START_REF] Kunert | Optimal crossover designs in a model with self and mixed carryover effects[END_REF]: in our model, the self carryover effect is equal to ρ d(i,j-1) + χ d(i,j-1)d(i,j) . We define the vector φ of total effects by φ = τ + λ + χ. It corresponds to the effect of a treatment when it is preceded by itself. Note that if

θ ′ = (τ ′ , λ ′ , χ ′ ) and K ′ = (I t |I t |I t ), then φ = K ′ θ.
Because the effect of having no treatment differs from the carryover effect of any treatment, we consider only designs with pre-periods, i.e. designs with one period, called pre-period, added at the beginning. On this preperiod, each subject receive a treatment but the response is not used in the analysis.

We assume that the designs are circular, i.e. the treatment assigned to a subject in the pre-period is the same as the treatment assigned in the last period. The circularity condition may be written d(i, 0) = d(i, k). We denote

Some technical tools

In this section, we present a general method to derive optimal crossover designs. We generalize and unify the techniques developed by [START_REF] Kushner | Optimal repeated measurements designs: the linear optimality equations[END_REF], [START_REF] Kunert | On the determination of optimal designs for an interference model[END_REF] and [START_REF] Bailey | Optimality of neighbor-balanced designs for total effects[END_REF].

Let denote by 1l k , I k and J k respectively the vector of ones of length k, the (k, k) identity matrix and the (k, k) matrix of ones. For any matrix A, denote by A + the Moore-Penrose inverse of A. The projection matrix onto the column span of matrix A is denoted by pr

(A) = A(A ′ A) + A ′ . We denote pr ⊥ (A) = I -pr (A) and Q k = pr ⊥ (1l k ) = I k -k -1 J k .

Information matrices and its extremal representation

Consider a generic partitioned linear model:

Y = Aα + Bβ + ε with E(ε) = 0 and Var(ε) = σ 2 I,
where α is a vector of length q. It is well known (see e.g. [START_REF] Kunert | Optimal design and refinement of the linear model with applications to repeated measurements designs[END_REF] that the information matrice C[α] of the parameter α is

C[α] = A ′ pr ⊥ (B) A. (3)
Consider now a subsystem K ′ α where K is a (q, s) matrix. The information matrix C[K ′ α] of K ′ α may be defined by the extremal representation [START_REF] Gaffke | Further characterizations of design optimality and admissibility for partial parameter estimation in linear regression[END_REF][START_REF] Pukelsheim | Optimal design of experiments[END_REF]:

C[K ′ α] = min L∈R q×s :L ′ K=Is L ′ C[α]L. (4)
where the minimum, taken relative to the Loewner ordering, exists and is unique. If L * is a (q, s) matrix that minimizes L ′ C[α]L under the constraint

L ′ K = I s , then C[K ′ α] = L * ′ C[α]L * .
(5) Formula (5) will be central in the following mainly because it provides a linear relationship between C[α] and

C[K ′ α].
The main issue to construct optimal designs will be to find L * for designs candidate to optimality. Because a global minimum in (4) exists,

L * minimizes L ′ C[α]L iff L * minimizes tr(L ′ C[α]L). So, we have tr(C[K ′ α]) = tr(L * ′ C[α]L * ) = min L∈R q×s :L ′ K=Is tr(L ′ C[α]L). (6)
In our problem, the matrix C[α] will have a natural block structure with completely symmetric blocks. In that case, the following results show that the matrice L * have also completely symmetric blocks.

Proposition 1 Let C[α] = ((C ij )) i,j=1,.
..,q be a block matrix, where the blocks C ij are (t, t) completely symmetric matrices for i = 1, ..., q and j = 1, ..., q. Let K ′ = (K ′ 1 , ..., K ′ q ) where for all i, K i are (t, t) completely symmetric matrices. Then C[K ′ α] is completely symmetric. Moreover, if

L * ′ = (L * 1 ′ , ..., L * q ′
) satisfies both Eq. (5) and the constraint L * ′ K = I t , then L * i , i = 1, ..., q, can be chosen completely symmetric.

Proof : The proof is given in Appendix A. ✷

The interest of this result is that

L * i = a * i I t + b * i J t can be found out by minimizing q(a * 1 , b * 1 , ..., a * q , b * q ) = tr(L * ′ C[α]L * ) which is a polynomial of second degree in a * 1 , b * 1 , ..., a * q , b * q .
Corollary 2 Under the notations and assumptions of Proposition 1, if moreover C ij 1l t = 0 for i = 1, ..., q and j = 1, ..., q, then L * i can be chosen equal to a * i I t for some scalar a * i .

Proof : By Proposition 1, L * i can be chosen equal to a * i I t + b * i J t . Since, The following lemma, although straightforward, is useful to calculate explicitly L * .

C ij J t = J t C ij = 0, L ′ i C ij L j = a * i a * j C ij
Lemma 3 Let C be a nonnegative (u, u) matrix, γ be a k-vector, and γ → L γ be a linear mapping with L γ a (u, v) matrix. Then, q(γ) = tr

L ′ γ CL γ is a multivariate convex quadratic polynomial in γ 1 , ..., γ k . Moreover, q(γ * )
is a minimum of q(γ) if and only if ∂q ∂γ (γ * ) = 0.

Some examples

Under the notations and assumptions of Proposition 1 and Corollary 2, i.e.

assuming that C ij are completely symmetric and that C ij 1l t = 0, we show in some cases how the matrices L * can be obtained. The three first examples give new presentations of known results. The last one will be used in this paper. We denote c ij = tr(C ij ).

Example 1 : let q = 2 and K ′ = (I t |0). By Corollary 2 and because L * ′ K = I t , L * ′ can be chosen to be equal to (I t |x * I t ) for some x * ∈ R. By

Eq. ( 6), x * can be found out by minimizing

q(x) = tr(L ′ C[α]L) = c 11 + 2x c 12 + x 2 c 22 .
This quadratic function was used by Kushner (1997, Eq. 4.1). The minimum is obtained for x * = c 12 c 22 and we found

C[K ′ α] = C 11 -C 12 C + 22 C 21 , the Schur- complement of C 22 .
Example 2 : let q = 3 and K ′ = (I t |0|0) . As in example 1, L * ′ can be chosen equal to (I t |x * I t |y * I t ), where x * and y * minimize the quadratic function

q(x, y) = tr(L ′ C[α]L) = c 11 + x 2 c 22 + y 2 c 33 + 2xc 12 + 2yc 13 + 2xyc 23 .
This quadratic function was used by Kunert and Martin (2000, Proposition 3).

Example 3 : let q = 2 and K ′ = (I t |I t ). We have

L * ′ = (x * I t |(1 -x * ) I t ), where x * minimizes q(x) = x 2 c 11 + (1 -x) 2 c 22 + 2x(1 -x)c 12 . If C 11 = C 22 , then x * = 1 2 and L * = 1 2 K. Therefore, C[K ′ α] = 1 4 K ′ C[α]K.
This equation was obtained in a different way by [START_REF] Bailey | Optimality of neighbor-balanced designs for total effects[END_REF] in order to construct optimal designs for total effects under models without interaction.

Example 4: let q = 3 and K ′ = (I t |I t |I t ). We assume that C 11 = C 22 and that C 13 = C 23 . We no longer assume that C 33 1l t = 0. We have

L * ′ = (x * 1 I t |x * 2 I t |(1 -x * 1 -x * 2 )I t + y * J t ), where x * 1 , x * 2 and y * minimize q(x 1 , x 2 , y) = (x 2 1 + x 2 2 )c 11 + 2(x 1 + x 2 )c 11 -(x 1 + x 2 ) 2 (2c 13 + c 33 ) +2(x 1 + x 2 )(c 13 -c 33 -c 33 ) +2y c 33 + ty 2 c 33 + c 33 (7) with c 33 = tr(J t C 33 ) = t -1 tr(J t C 33 J t ). Note that L * does not satisfy L * ′ K = I t . This is in fact a simplified form of L * ′ = (x * 1 I t + u * J t |x * 2 I t -(u * + v * )J t |(1 -x * 1 -x * 2 )I t + y * J t ) ,
noting that the terms u * J t and (u * + v * )J t vanish in the expression of

L * ′ C[α] L * .
By lemma 3 and by symmetry of q(x 1 , x 2 , y) in x 1 , x 2 , it is easy to see that x * 1 and x * 2 can be chosen equal. Denote

x * = x * 1 = x * 2 . From
Eq. ( 7), x * and y * can be found out by minimizing

q (x, y) = c 33 + 4 (c 13 -c 33 ) x + 2 c 33 y + 2 (c 11 + c 12 -4c 13 + 2c 33 ) x 2 + t c 33 y 2 -4 c 33 xy, (8)
and the minimum q(x * , y * ) is equal to tr(C[K ′ α]).

4 Optimal circular crossover designs for total effects

From [START_REF] Kiefer | Construction and optimality of generalized Youden designs[END_REF], a design d * for which the information matrix

C d * [φ] is
completely symmetric and that maximizes the trace of C d [φ] over all the designs d in Ω t,b,k is universally optimal. In this section, we propose a method to construct universal optimal design for total effects.

Upper bound of tr C d [φ]

For a design d, the information matrix for the whole parameter 1) is:

θ ′ = (τ ′ , λ ′ , χ ′ ) in Model (
C d [θ] =         T ′ d pr ⊥ (B) T d T ′ d pr ⊥ (B) L d T ′ d pr ⊥ (B) S d L ′ d pr ⊥ (B) T d L ′ d pr ⊥ (B) L d L ′ d pr ⊥ (B) S d S ′ d pr ⊥ (B) T d S ′ d pr ⊥ (B) L d S ′ d pr ⊥ (B) S d         =         C d11 C d12 C d13 C ′ d12 C d22 C d23 C ′ d13 C ′ d23 C d33        
.

By circularity of the designs,

B ′ T d = B ′ L d , T ′ d T d = L ′ d L d and T ′ d S d = L ′ d S d = S ′ d S d = S ′ d L d = S ′ d T d and therefore C d22 = C d11 and C d23 = C d13 . As in Section 3.2, we define c dij = tr(C dij ) and c dij = tr(JC dij ).
A design is said to be symmetric if all the blocks C dij are completely symmetric, or equivalently, if C d [θ] is invariant by any permutations of the treatment labels.

Proposition 4 For any design d, the information matrix for total effects φ satisfies:

tr C d [φ] ≤ min x,y q d (x, y)
where

q d (x, y) is defined by (8). Equality holds if C dij [θ] are completely symmetric for 1 ≤ i ≤ 3 and 1 ≤ j ≤ 3. Proof : For a design d, we denote C d [θ] = 1 t! σ∈St (I 3 ⊗ P σ ) C d [θ] (I 3 ⊗ P ′ σ )
the symmetrized information matrix for θ. By conctruction, C dij are completely symmetric for all i and j. Moreover,

c dij = tr(C dij ) = tr(C dij ) and c dij = tr(J t C d33 ) = tr(J t C d33 ). Denote by C d [φ] the information matrix of φ associated to C d [θ]. By concavity of C d [φ] w.r.t. C d [θ], (see Pukelsheim, 1993 p. 77), C d [φ] ≤ C d [φ] and tr C d [φ] ≤ min x,y q d (x, y). ✷ A design is called degenerate if c d33 = tr(J t C d33 [θ]) = 0 or equivalently if
each block either contains only one treatment or has no treatment preceded by itself. Note that the information matrix of a degenerate design is null and therefore such a design is not considered. The following lemma shows that the minimization of q d (x, y) may be reduce to the minimization of a one variable quadratic function.

Lemma 5 Let d be a non-degenerate design. The values x * and y * that minimize q d (x, y) satisfy 2x * -ty * -1 = 0. Moreover, x * minimizes q d (x),

where

q d (x) = c d33 - 1 t c d33 +4(c d13 -c d33 + 1 t c d33 )x+2(c d11 +c d12 -4c d13 +2c d33 - 2 t c d33 )x 2 , (9)
and q d (x * , y * ) = q d (x * ).

Proof : by lemma 3, ∂q d ∂y (x * , y * ) = 0 at the minimum. For a non-degenerate design, this is equivalent to 2x * -ty * -1 = 0. Replacing y by 2(x -1)/t, we found q d (x, y) = q d (x). Therefore, x * necessarily minimizes q d (x). ✷

Decomposition over the blocks

It is well known that C d [θ] is the sum of the information matrices C du corresponding to Block u. Denote by T du , L du and S du the incidence matrices for Block u. Thus,

T ′ d = (T ′ d1 | ... | T ′ db ) , L ′ d = (L ′ d1 | ... | L ′ db ), S ′ d = (S ′ d1 | ... | S ′ db ) and: C d [θ] = b u=1 C du [θ] = b u=1         T ′ du Q k T du T ′ du Q k L du T ′ du Q k S du L ′ du Q k T du L ′ du Q k L du L ′ du Q k S du S ′ du Q k T du S ′ du Q k L du S ′ du Q k S du        
We decompose in the same way tr (C dij ) and tr (C dij J t ):

c dij = tr (C dij ) = b u=1 c (u) dij and c dij = tr (C dij J t ) = b u=1 c (u) dij
denoting by c

c d11 = tr (C d11 ) = tr b u=1 T ′ du Q k T du = b u=1 c (u) d11 with c (u) d11 = tr(T ′ du Q k T du ).
Simplifications of these forms give:

c (u) d11 = k - n s u k , c (u) d12 = m u - n s u k , c (u) d13 = m u - l u k , c (u) d33 = m u - m s u k , c (u) d33 = m u - m 2 u k , with n s u = t i=1 n 2 ui , m u = t i=1 m ui , m s u = t i=1 m 2 ui , l u = t i=1 n ui m ui ,
denoting by n ui the number of plots in Block u which receive Treatment i and by m ui the number of times Treatment i is preceded by itself in the Block u. It follows that:

q d (x) = b u=1 h (u) d (x)
where

h (u) d (x) = c (u) d33 - 1 t c (u) d33 +4(c (u) d13 -c (u) d33 + 1 t c (u) d33 )x+2(c (u) d11 +c (u) d12 -4c (u) d13 +2c (u) d33 - 2 t c (u) d33 )x 2 .
Two sequences of treatments in two blocks u 1 and u 2 are said to be equiv- the same for all u from one equivalent class ℓ, we change the notation and write n s ℓ , m ℓ , m s ℓ , l ℓ and c dij (ℓ) instead. We define:

alent if n s u 1 , m u 1 , m s u 1 , l u 1 = n s u 2 , m u 2 , m s u 2 , l u 2 ,
h ℓ (x) = m ℓ - m s ℓ k - δ ℓ t + 4 k m s ℓ -l ℓ + kδ ℓ t x +2 (k -m ℓ ) + 2 k (2l ℓ -n s ℓ -m s ℓ ) - 2δ ℓ t x 2
where δ ℓ = m ℓ (1 -m ℓ /k) . For a design d, we denote by π dℓ the proportion of blocks assigned to the class ℓ (1 ≤ ℓ ≤ K). So, we have:

q d (x) = b K ℓ=1 π dℓ h ℓ (x) . (10)

Optimality conditions

Our goal, in order to find an optimal configuration d * , is to obtain a value

x * and proportions π ′ d * = (π d * 1 , ..., π d * K ) such that:

q d * (x * ) = max π d q * d with q * d = min x q d (x).
The following propositions characterizes universally optimal designs.

Proposition 6 [START_REF] Kunert | On the determination of optimal designs for an interference model[END_REF] Consider a symmetric design d * ∈ Ω t,b,k and a point x * such that the first derivative of q d * is zero.

If we have also:

∀ ℓ = 1, ..., K , bh ℓ (x * ) ≤ q * d * ,
then d * is universally optimal over Ω t,b,k .

The following proposition reduce to k-1 the number of equivalent classes that possibly appear in an optimal design. We denote by ⌊x⌋ the integer part of x.

Proposition 7 If a symmetric design d * ∈ Ω t,b,k is universally optimal, then the treatment sequences present in the design necessarily satisfies:

1. all the periods receiving the same treatment are side by side in the sequence.

2. each treatment present in the sequence occurs ⌊k/v⌋ or ⌊k/v⌋+1 times, where v is the number of different treatments present in the sequence.

The number of treatments that occur ⌊k/v⌋ + 1 times is k -v⌊k/v⌋ and the number of treatments that occurs ⌊k/v⌋ times is v(⌊k/v⌋ + 1) -k.

Proof : The proof is given in Appendix B. ✷

It is worth noting that for each value of v, there is only one equivalent class of treatment sequences that may appear in the optimal design.

Construction of optimal designs

Consider a sequence ℓ present in an optimal design. From Proposition 7, a treatment i in ℓ satisfies m ℓi = n ℓi -1. Moreover we have v ℓ = k -m ℓ . So

n s ℓ = m s ℓ + m ℓ + k, l ℓ = m s ℓ + m ℓ and then: h ℓ (x) = m ℓ - m s ℓ k - δ ℓ t + 4 k kδ ℓ t -m ℓ x+2 (k -m ℓ ) k -2 k - 2δ ℓ t x 2 , with δ ℓ = m ℓ (1 -m ℓ /k) .
Therefore an optimal design for this one-dimensional problem is obtained using one class or a mixture of two different classes (see [START_REF] Kushner | Optimal repeated measurements designs: the linear optimality equations[END_REF]).

The following method (for given k and t) can be used in order to prove that a design d * is optimal.

• If the optimal design d * is constituted by one sequence of treatment ℓ 1 (i.e. q d * (x) = bh ℓ 1 (x)):

-find x * that minimize h ℓ 1 and then the minimum q * d * of q d * ,

-check that: ∀ ℓ = 1, ..., K , bh ℓ (x * ) ≤ q * d * (cf. Prop. 6).

• If the design d * is constituted by two different classes of treatments ℓ 1 and ℓ 2 (i.e.

q d * (x) = b (π d * ℓ 1 h ℓ 1 (x) + π d * ℓ 2 h ℓ 2 (x))):
-find an admissible intersection point x * according to the definition of [START_REF] Kushner | Optimal repeated measurements designs: the linear optimality equations[END_REF], i.e. h ℓ 1 (x * ) = h ℓ 2 (x * ) and

∂h ℓ 1 ∂x (x * ) ∂h ℓ 2 ∂x (x * ) ≤ 0,
-find the proportions π * d * ℓ 1 and π * d * ℓ 2 in order to obtain: ∂q d * ∂x (x * ) = 0, -use these proportions for finding the minimum q * d * = q d * (x * ) of q d * ,

-check that: ∀ ℓ = 1, ..., K , bh ℓ (x * ) ≤ q * d * (cf. Prop. 6).

Note that the optimal proportions can be found by the following method.

∂q d * ∂x (x * ) = 0 ⇐⇒ b π d * ℓ 1 ∂h ℓ 1 ∂x (x * ) + π d * ℓ 2 ∂h ℓ 2 ∂x (x * ) = 0.
Denote a i = (∂h ℓ i /∂x) (x * ), i = 1, 2. The optimal proportions are then:

π * d * ℓ 1 = a 2 a 2 -a 1 and π * d * ℓ 2 = a 1 a 1 -a 2 .
Now, we give optimal designs for several values of k. For k = 3, 4 we give explicit formulae. For k = 5, ..., 10, we present numerical results.

1) The case k = 3. From proposition 7, an optimal design is constituted by one or two blocks in the following set of sequences:

Sequence n s ℓ m ℓ m s ℓ l ℓ [ 1 2 3 ] 3 0 0 0 [ 1 1 2 ] 5 1 1 2
So the functions h ℓ are given by (we identify each class with its value of m ℓ ):

h 0 (x) = 2x 2 and h 1 (x) = 2 3 1 - 1 t - 4 3 1 - 2 t x + 4 3 1 - 2 t x 2 .
It is impossible in that case to satisfy Proposition 6 using only one of these two sequences. So, we must find an admissible intersection point x * . Some algebra shows that:

x * = (t -2) 2 + (t -1) (t + 4) -(t -2) t + 4 .
The proportions in the optimal design d * are then: -1) .

π * d * 0 = 1 -π * d * 1 and π * d * 1 = 3t t + 4 1 - (t -2) t (2t
The following table gives the optimal proportions for several values of t. The sequence [ 1 1 2 ] is predominating in this mixture. So it can be interesting in practice to use designs generated by only this sequence. The quality of such designs can be quantified by the classical Φ p criteria. We know (see e.g. [START_REF] Druilhet | Conditions for optimality in experimental designs[END_REF]) that when the information matrix is completely symmetric Φ p does not depend on p. Thus we can derive the efficiency factor of a design d ∈ Ω t,b,3 generated by one sequence ℓ:

Eff (d) = b h * ℓ q d * (x * ) with h * ℓ = min x h ℓ (x) .
Numerical applications are given in the following table:

t 3 4 5 10 20 ∞ Eff.
[ 1 1 2 ] 0.989 0.985 0.982 0.976 0.974 0.971

2) The case k = 4. An optimal configuration is constituted by one or two blocks belonging to the following set of sequences:

Sequence n s ℓ m ℓ m s ℓ l ℓ [ 1 2 3 4 ] 4 0 0 0 [ 1 1 2 3 ] 6 1 1 2 [ 1 1 2 2 ] 8 2 2 4
The functions h ℓ are then given by:

h 0 (x) = 4x 2 , h 1 (x) = 3 4 1 - 1 t -1 - 3 t x + 3 1 - 1 t x 2 , h 2 (x) = 3 2 - 1 t -2 1 - 2 t x + 2 1 - 2 t x 2 .
For any t, the minimum of h 2 is obtained for x * = 0.5 and satisfies the conditions of proposition 6. So the optimal design is generated by the sequence [ 1 1 2 2 ]. As an example we can consider, for t = 4, the optimal design such that:

D =             1 2 1 3 1 4 2 3 2 4 3 4 1 2 1 3 1 4 2 3 2 4 3 4 2 1 3 1 4 1 3 2 4 2 4 3 2 1 3 1 4 1 3 2 4 2 4 3             ∈ Ω 4,12,4 .
Note that, by circularity, the design obtained by taking away one out of every two columns is universally optimal over all the designs in Ω 4,6,4 .

3) The case k = 5. The optimal design is generated by the following mixtures:

t 3 4 5 10 20 ∞ Prop.
[ 1 1 2 2 3 ] 0.167 0.250 0.300 0.400 0.450 0.500

Prop.

[ 1 1 1 2 2 ] 0.833 0.750 0.700 0.600 0.550 0.500

The efficiencies of designs generated by only one sequence are:

t 3 4 5 10 20 ∞ Eff.
[ 1 1 2 2 3 ] 0.844 0.895 0.918 0.949 0.961 0.970

Eff.

[ 1 1 1 2 2 ] 0.984 0.977 0.972 0.963 0.959 0.955 It can be observed that the design generated by [ 1 1 1 2 2 ] is more efficient than the one generated by [ 1 1 1 2 2 ] for t > 18.

4) The case k = 6. The optimal design is generated by the following mixtures:

t 3 4 5 6 7 ≥ 8
Prop.

[ 1 1 2 2 3 3 ] 0.400 0.628 0.775 0.878 0.954 1.000

Prop.

[ 1 1 1 2 2 2 ] 0.600 0.372 0.225 0.122 0.046

The efficiencies of designs generated by only one sequence are:

t 3 4 5 6 7 ≥ 8
Eff.

[ 1 1 2 2 3 3 ] 0.962 0.989 0.997 0.999 0.999 1.000

Eff.

[ 1 1 1 2 2 2 ] 0.962 0.942 0.930 0.922 0.917

We note that the sequence [ 1 1 2 2 3 3 ] is always more efficient than the sequence [ 1 1 1 2 2 2 ]. It also generates an optimal design when t ≥ 8.

5)

The case k = 7. The optimal design is generated by the following mixtures:

t 3 ≥ 4
Prop.

[ 1 1 1 2 2 3 3 ] 0.682 1.000

Prop.

[ 1 1 1 1 2 2 2 ] 0.318

The efficiencies of designs generated by only one sequence are:

t 3 ≥ 4
Eff.

[ 1 1 1 2 2 3 3 ] 0.994 1.000

Eff.

[ 1 1 1 2 2 3 3 ] 0.938

We note that he sequence [ 1 1 1 2 2 3 3 ] is always better than the sequence [ 1 1 1 1 2 2 2 ] and generates the optimal design when t ≥ 4.

6)

The case k = 8, 9, 10. We find that the optimal design is generated, for every t, by only one sequence. These optimal sequences are given by: t 8 9 10

Sequence [ 1 1 1 2 2 2 3 3 ] [ 1 1 1 2 2 2 3 3 3 ] [ 1 1 1 1 2 2 2 3 3 3 ] 5 Models with period effects

We consider the model:

y ij = α j + β i + τ d(i,j) + λ d(i,j-1) + χ d(i,j-1)d(i,j) + ε ij , (11)
where α j is the effects of the period j. We denote by A the corresponding incidence matrix. A balanced design is a design such that all the sequences belonging to the same equivalent class appear equally often. This notion is slightly more restrictive than the notion of symmetric design defined in Section 4.1 but correspond to the design generated by one sequence or a mixture of sequences, considering all the treatment permutations.

Proposition 8 A balanced design which is universally optimal for total effects under Model (1) is also universally optimal under Model (11).

Proof : The proof is a direct consequence of the following lemma. ✷ Lemma 9 For a balanced design, the information matrix for the total effects is the same under model ( 1) and ( 11)

Proof : Let d be a balanced design. The difficulty is that interaction and period effects are not orthogonal. Write θ′ = (τ |λ|χ|α) and Cd [ θ] the corresponding information matrix. We want to prove that Cd

[ K′ θ] = C d [K ′ θ]
where K′ = (I t |I t |I t |0). We write:

Cd [ θ] =     C d [θ] C d12 C d21 C d22    
where

C ′ d12 = (0|0|D ′ d ) with D d = (p d1 1l t |p d2 1l t |...|p db 1l t ).
The real p dj depends only on the number of times a treatment is preceded by itself on period j. Note that the usual orthogonality condition (see [START_REF] Kunert | Optimal design and refinement of the linear model with applications to repeated measurements designs[END_REF] between interaction and period effects is D d = 0 which is not the case here.

The key point of the proof is that

Q t D d = 0. ( 12
that C d [K ′ θ]1l t = 0 and, since Cd [ K′ θ] ≤ C d [K ′ θ], we have Cd [ K′ θ]1l t = 0. So, C d [ K′ θ] = Q t Cd [ K′ θ]Q t . Therefore, L * can be chosen equal to (M * d Q t |N * d Q t ).
For any permutation σ, we have (I 3 ⊗ P σ )C d12 = C d12 and therefore, similarly to Proposition 1, it can be shown that the three (t × t) blocks of M * d can be chosen completely symmetric.

So, L * d ′ can be chosen equal to (x * 1 Q t |x * 2 Q t |(1 -(x * 1 + x * 2 ))Q t |N * ). Put L′ = (M ′ |N ′ ) with M ′ = (x 1 Q t |x 2 Q t |(1 -(x 1 + x 2 )Q t )
. By (12), we have

M ′ C d12 = 0 and then Cd [ K′ θ] = min x 1 ,x 1 ,N L′ C[ θ] L = min x 1 ,x 1 ,N (M ′ C d [θ]M + N ′ C d12 N ).
Since x 1 , x 2 and N vary freely, N can be chosen equal to 0 and therefore, from Example 4 of Section 3.2 and the constraint on y obtained in Lemma

5, M * ′ C[ θ]M * = L * ′ C d [θ]L * . The result follows. ✷ Appendix A
We give here the proof of Proposition 1.

Step 1: We show that the set E of matrices L that minimize L ′ C[α]L under the constraint L ′ K = I t is an affine subspace: consider L a and L b in E with

L a = L b . It is sufficient to prove that L γ = γL a + (1 -γ)L b belongs to E
for any real γ. We have L ′ γ K = I t and

L ′ γ C[α]L γ = γ 2 L ′ a C[α]L a +(1-γ) 2 L ′ b C[α]L b +γ(1-γ)(L ′ a C[α]L b +L ′ b C[α]L a )
This quadratic function w.r.t. γ admits two minima relative to the Loewner ordering at γ = 0 and γ = 1 and therefore it is constant. So, L γ minimizes

L ′ C[α]L for any γ ∈ R.
Step 2: Consider a permutation σ on {1, ..., t} and denote by P σ the corresponding permutation matrix. We want to prove that if L * ∈ E, so does L * σ = (I 3 ⊗ P ′ σ ) L * P σ . Since C ij and K i are completely symmetric:

(I 3 ⊗ P σ )C[α](I 3 ⊗ P ′ σ ) = C[α], (13)
and

(I 3 ⊗ P ′ σ )K = KP σ . (14)
Moreover,

C[K ′ α] = min L∈R 3t×t :L ′ K=It L ′ C[α]L, = min L∈R 3t×t :L ′ K=It L ′ (I t ⊗ P σ )C[α](I t ⊗ P ′ σ )L (by (13)), = P σ min L∈R 3t×t :L ′ K=I t P ′ σ L ′ (I t ⊗ P σ )C[α](I t ⊗ P ′ σ )LP σ P ′ σ .
Put L σ = (I 3 ⊗ P ′ σ )LP σ . By ( 14), L ′ σ K = I t . Since, L ↔ L σ is a one to one mapping, we have:

C[K ′ α] = P σ min Lσ∈R 3t×t :L ′ σ K=It L ′ σ C[α]L σ P ′ σ = P σ C[K ′ α]P ′ σ . (15)
Since Eq. ( 15) holds for any permutation σ, C

[K ′ θ] is completely symmetric. So, C[K ′ α] = min Lσ∈R 3t×t :L ′ σ K=It L ′ σ C[α]L σ and therefore C[K ′ α] = L * σ ′ C[α]L * σ .
Step 3: If L * ∈ E then, by step 1 and 2, L * = 1 t! σ L * σ also belongs to E.

By construction, L *

i is completely symmetric and the proof is complete.

that h ℓ (x) admits a minimum. In that case, any x is a minimum and can be chosen positive. We assume now that C ℓ > 0. By Eq. 16, we have The proof of the Proposition 7 is given below.

Step 1: Let ℓ be a sequence containing v ℓ different treatments numbered 1, ..., v ℓ and ℓ the associated sequence defined above. We want to prove that .

Since n e ℓi = m e ℓi + 1 we have:

n s ℓ -2l ℓ +m s ℓ = v ℓ i=1 (n ℓi -m ℓi ) 2 = v ℓ i=1 γ 2 i ≥ v ℓ i=1 γ i = k-m ℓ = v e ℓ i=1
n e ℓi -m e ℓi 2 with equality if and only if ℓ = ℓ. The result follows.

• From Lemma 10, we have B ℓ < B e ℓ .

• We show that A ℓ < A e ℓ . It is equivalent to show that m s ℓ > m s e ℓ

. This follows from the fact that i m ℓi = i m e ℓi and, for all i, m ℓi ≥ m e ℓi with equality if and only if ℓ = ℓ.

Step 2: Consider the set L of sequences ℓ having the same number v ℓ = v of distinct treatments and such that ℓ = ℓ. For a sequence ℓ ∈ L, n ℓi = m ℓi -1 and m ℓ = k -v. So we have for every ℓ 1 , ℓ 2 ∈ L (see Section 5.2):

h ℓ 1 (x) -h ℓ 2 (x) = 1 k m s ℓ 2 -m s ℓ 1 .
If there exists a sequence ℓ * ∈ L such that h ℓ * (x) -h ℓ (x) > 0 (for every ℓ = ℓ * ) then only this sequence can belong to the optimal design. Such a sequence minimizes m s ℓ = v ℓ i=1 m 2 ℓi under the constraint v ℓ i=1 m ℓi = m ℓ .

So, it is well known that the m ℓ * i or equivalently the n ℓ * i (i = 1, ..., v ℓ * ) must be as equal as possible: if k/v ℓ * is an integer, a treatment present

  For a square matrix A, we denote by tr(A) the trace of A. For two symmetric matrices M and N , M ≤ N means that N -M is a nonnegative definite matrix (Loewner ordering). A matrix M is completely symmetric if M = aI + bJ for some scalars a and b.

  does not depends on b * i and b * j and the result follows. ✷

  dij the contributions of Block u. For example, for i = j = 1 we have:

  which is the case if one sequence is obtained from the other one by relabelling the treatments. So, for given k and t, we can divide the set of all possible treatment sequences into K equivalence classes of treatments. Since n s u , m u , m s u , l u and c (u) dij are

′

  ) where M d is a (3t×t) matrix and N * d is a (b×b) matrix. It is easy to see

l

  ℓ = i m ℓi n ℓi ≥ m ℓ + m s ℓ and thenB ℓ ≤ 4 kt m ℓ (k -m ℓ -t)with equality if and only if the occurrences of the treatments present in the sequence are placed side by side. In that case a sequence ℓ satisfies k -m e ℓ = v e ℓ and, since t ≥ v e ℓ , B e ℓ ≤ 0. Consider now a general sequence ℓ and denote by ℓ its associated sequence with side-by-side treatments. Since m e ℓ = m ℓ we have also δ ℓ = δ e ℓ and then B ℓ ≤ B e ℓ if and only if (l ℓ -m s ℓ ) ≥ l e ℓand only if ℓ = ℓ. So we have for every sequence ℓ: B ℓ ≤ B e ℓ ≤ 0

  h e ℓ (x) > h ℓ (x) for all x ≥ 0 and ℓ = ℓ. The idea is to show that A ℓ < A e ℓ , B ℓ < B e ℓ and C ℓ < C e ℓ . Note that m ℓ = m e ℓ so δ ℓ = δ e ℓ . Then:• We show that C ℓ < C e ℓ .It is equivalent to show that (n s ℓ -2l ℓ + m s ℓ )

by Ω t,b,k the set of all circular designs with t treatments, b subjects and k periods.

in ℓ * occurs k/v ℓ * times; otherwise, a treatment present in ℓ * occurs either ⌊k/v ℓ * ⌋ or ⌊k/v ℓ * ⌋ + 1 times in ℓ *

Appendix B

We prove here Proposition 7. From Eq. ( 10), we have h ℓ (x) = A ℓ + B ℓ x + C ℓ x 2 where:

. We denote by v ℓ the number of treatments present in the sequence ℓ. For any treatment i in ℓ, we have:

with equality if and only if treatment i occurs in consecutive periods. In general, each treatment appears in the sequence ℓ in several groups of consecutive periods. The number of such groups is

we construct a new sequence ℓ as follows: in each group of consecutive periods receiving the same treatment, we replace the treatment that occurs in this group by a new treatment, such that in the new sequence ℓ, each treatment have only one group. For example ℓ = (1, 1, 2, 2, 2, 1, 3, 2, 2) gives ℓ = (1, 1, 2, 2, 2, 4, 3, 5, 5). We have m e ℓ = m ℓ and v e ℓ = γ = k -m ℓ , where γ = i γ i is the number of different groups in ℓ.

Lemma 10 If d * is optimal, then x * that minimizes q d * (x) is positive.

Proof. From Lemma 3 we know that h ℓ (x) is convex, then C ℓ ≥ 0. So it is sufficient to show that B ℓ ≤ 0. If C ℓ = 0 then B ℓ = 0 because we know


Mohamed Daoudi

Zakia Hammal
email: zakiahammal@yahoo.fr

Anis Kacem

Jeffrey F Cohn

Gram Matrices Formulation of Body Shape Motion: An Application for Depression Severity Assessment

Keywords: Gram matrices, body movement, dynamics, depression severity

We propose an automatic method to measure depression severity from body movement dynamics in participants undergoing treatment for depression. Participants in a clinical trial for treatment of depression were interviewed on up to four occasions at 7-week intervals with the clinicianadministered Hamilton Rating Scale for Depression. Body movement was tracked using OpenPose from full-body video recordings of the interviews. Gram matrices formulation was used for body shape and trajectory representations from each video interview. Kinematic features were extracted and encoded for video based representation using Gaussian Mixture Models (GMM) and Fisher vector encoding. A multi-class SVM was used to classify the encoded body movement dynamics into three levels of depression severity: severe, mild, and remission. Accuracy was high for severe depression (68.57%) followed by mild depression (56%), and then remission (37.93%). The obtained results suggest that automatic detection of depression severity from body movement is feasible.

I. INTRODUCTION

Depression is one of the most common psychiatric disorders worldwide. According to the Diagnostic and Statistical Manual of Mental Disorders (DSM-5) behavioral indicators of depression include reduced facial expressiveness, decreased vocal intensity and speech, inability to sit still, pacing, handwringing, and slowed body movements . Yet, often these indicators are not taken into account in depression severity assessment that relies almost entirely on patients' verbally reported symptoms in clinical interviews. To improve depression severity assessment, automatic and objective assessment of depression from behavioral signals has emerged as a powerful option.

Most previous efforts in automatic assessment of depression severity have focused on facial expression [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF], [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF], voice quality and timing [START_REF] Scherer | Dyadic behavior analysis in depression severity assessment interviews[END_REF], and head pose [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF], [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF], [START_REF] Alghowinem | Head pose and movement analysis as an indicator of depression[END_REF]. For instance, Hdibeklioglu and colleagues [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF] proposed a state of the art multimodal (face, head, and voice) deep learning based approach to detect depression severity in participants undergoing treatment for depression. The dynamics of facial, head, and vocal prosody was important for the measurement * Co-first authorship: M. Daoudi and Z. Hammal contributed equally to this paper. of depression severity, but one could not say whether the dynamics of those measures was increasing, decreasing, or varying in some non-linear way. To overcome that limitation, Kacem and colleagues [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF] proposed a method to measure depression severity from facial and head movement dynamics using affine-invariant barycentric and Lie algebra representation, respectively. Consistent with clinical data, the extracted kinematic features revealed that the velocity and acceleration of facial movement strongly mapped onto depression severity symptoms [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF].

Although body movement also changes with depression severity (e.g., inability to sit still, slowed body movement [DSM-5]), the communicative functions of body movement dynamics in relation to depression severity have been fairly unexplored using automatic approaches. One notable exception is [START_REF] Joshi | Can body expressions contribute to automatic depression analysis?[END_REF]. In [START_REF] Joshi | Can body expressions contribute to automatic depression analysis?[END_REF], the authors investigated body movement for the detection of depression. To do so, the authors extracted the relative orientation and radius of body parts as well as the holistic body motion measured using space-time interest points [START_REF] Laptev | On space-time interest points[END_REF]. The two set of descriptors were then fused to train a support vector machine (SVM) classifier for the detection of presence from the absence of depression. From a clinical perspective, it is also critical to measure change over time in depression severity. In the continuity with the initial efforts listed above, we investigated the discriminating power of body movement dynamics for depression severity assessment.

Compared with previous efforts, including ours [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF], [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF], we focused on the unique contribution of body movement dynamics in clinically relevant participants in treatment of moderate to severe depression. The combination of body movement with other modalities (e.g., face and head [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF], [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF]) are subject to future work. To capture changes in the dynamics of body movement that would reflect the psychomotor retardation and agitation of depressed participants, Gram matrices formulation is used for body shape and trajectories representation. Relevant kinematic features are then extracted from body shape trajectories (i.e., velocities amics and accelerations). Gaussian Mixture Models (GMM) combined with an improved fisher vector encoding are then used to obtain a single vector representation for each sequence (i.e., clinical interview). Finally, a multi-class SVM with a Gaussian kernel is used to classify the encoded body movement dynamics into three depression severity levels: severe, mild, and remission. To the best of our knowledge, this is the first time that depression severity rather than only presence from absence of depression has been measured automatically from body movement dynamics. The proposed approach is illustrated in Fig. 1.

II. PARTICIPANTS

Fifty-seven participants with moderate to severe depression (34 women, 23 men) were recruited from a clinical trial for treatment of depression [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF]. At the time of the study, all met DSM-4 criteria [START_REF] First | Structured clinical interview for DSM-IV axis I disorders -Patient edition (SCID-I/P, Version 2.0)[END_REF] for Major Depressive Disorder (MDD). Data from 49 participants was available for analysis. Participant loss was due to change in original diagnosis (e.g., severe suicidal ideation) and methodological reasons (e.g., missing or corrupted video recordings). Symptom severity was evaluated on up to four occasions at 1, 7, 13, and 21 weeks post diagnosis and intake by four clinical interviewers (the number of interviews per interviewer varied) [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF].

Interviews were conducted using the Hamilton Rating Scale for Depression (HRSD) [START_REF] Hamilton | A rating scale for depression[END_REF]. HRSD is a clinician-rated multiple item questionnaire to measure depression severity and response to treatment. HRSD scores of 15 or higher are generally considered to indicate moderate to severe depression; scores between 8 and 14 indicate mild depression; and scores of 7 or lower indicate remission [START_REF] Fournier | Antidepressant drug effects and depression severity: A patient-level meta-analysis[END_REF] [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF].

Participant and interviewer were sitting face-to-face during the entire interview (see Fig. 2). Interviews were recorded using three hardware synchronized analogue cameras and two unidirectional microphones (see Fig. 2). Audio data was not used in the current paper. For video data, one camera recorded the interviewer's shoulders and face from approximately 15 degrees to the participant's right. Two additional cameras were positioned approximately 15 degrees to the interviewers left and right (see Fig. 2). One camera recorded the participants face and one camera recorded a full body view (see Fig. 2). In the current work, we used the data recorded from the full body view camera digitized at a frame rate of 29.97 fps (see Fig. 2). Using these data and the cut-off scores of depression severity described above, we defined three ordinal depression severity classes: moderate to severe depression, mild depression, and remission (i.e., recovery from depression). The final sample from the body camera view was 89 available full body recorded sessions from 49 participants: 35 moderate to severely depressed, 25 mildly depressed, and 29 remitted.

III. AUTOMATIC BODY TRACKING

We used a real-time fully-automatic multi-person 2D body tracker (OpenPose, [START_REF] Cao | Realtime multi-person 2d pose estimation using part affinity fields[END_REF]) to track body parts' locations. Open-Pose uses a non-parametric representation (key-points association that encodes both position and orientation of human limbs) to learn the association between body parts' locations (i.e., Head, Shoulders, Elbows, Wrists, Hips, Knees, and Ankles, see Fig. 3) and individuals present in a given image. The tracker encodes global context, allowing a greedy bottomup parsing step that maintains high accuracy of body parts tracking while achieving real time performance, irrespective of the number of people in the image [START_REF] Cao | Realtime multi-person 2d pose estimation using part affinity fields[END_REF]. OpenPose reliability has been evaluated with state-of-the-art accuracy on multiple public benchmarks [START_REF] Cao | Realtime multi-person 2d pose estimation using part affinity fields[END_REF]. Table I reports the proportions of tracked video frames for each body point for all data. Valid tracking was lower for points 10 and 13 (left ad right ankles, respectively). Self-occlusion (occluded ankles) contributed to tracking failure of points 10 and 13.

Because we are interested in body points only, tracked points 14, 15, 16, 17 (i.e., eyes, nose, and ears, respectively) were not used for analyses. Because points 10 and 13 are usually not detected, they were not used for analyses. This resulted in skeletons of 12 body joints used for the analysis of body movement dynamics in participants in treatment for moderate to severe depression (for a more detailed description of the data please see section Participants).

IV. MODELING BODY MOVEMENT DYNAMICS USING GRAM MATRICES FORMULATION

Evolutionary theories of depression highlights symptoms of psychomotor retardation and agitation [START_REF] Fowles | A motivational theory of psychopathology[END_REF], [START_REF] Nesse | Is depression an adaptation?[END_REF]. We used dynamic measures of body movement to capture aspects of psychomotor retardation and agitation in participants in treatment for moderate to severe depression.

A. Frame Based Body Shape Representation

Body movement dynamics is represented using the time series of the coordinates of the 12 tracked body points (i.e., p 1 = (x 1 , y 1 ), . . . , p 12 = (x 12 , y 12 )) during each video sequence (see section Automatic Body Tracking). Each video sequence is thus characterized by a set of landmark configurations {Z 0 , . . . , Z τ }, where τ is the duration in frames of the video sequence (interview). Each configuration Z i (1 ≤ i ≤ τ ) is an n × d, where d = 2 and n = 12, matrix of rank d encoding the positions of n distinct landmark points in d dimensions. We seek to measure the dynamic changes of the curves formed by the landmark configurations invariant to rotation and translation. To do so, we computed the matrix of pairwise distances between the landmarks of the same body shape (at a given frame). The measured distances are then centered by subtracting the landmarks center of mass referred to as p 0 . We then compute the Gram matrices as the squares of all measured centered distances such as:

G = ZZ T = p i , p j , 1 ≤ i, j ≤ n , (1)
The Gram matrices of the form ZZ T , where Z is an n × d matrix of rank d, are characterized as n × n positive semidefinite matrices of rank d. Conveniently, the Riemannian geometry of the space of these matrices, called the positive semidefinite cone S + (d, n), was well studied in [START_REF] Bonnabel | Riemannian metric and geometric mean for positive semidefinite matrices of fixed rank[END_REF]- [START_REF] Massart | Quotient geometry with simple geodesics for the manifold of fixed-rank positive-semidefinite matrices[END_REF].

B. Video Based Body Shape Trajectories

The dynamic changes of body joints movement are characterized by parametric trajectories on the Riemannian manifold of positive semi-definite matrices of fixed-rank (see Fig. the data are well acquired (i.e., without tracking errors or missing data) [START_REF] Kacem | A novel geometric framework on gram matrix trajectories for human behavior understanding[END_REF], [START_REF] Kacem | A novel space-time representation on the positive semidefinite cone for facial expression recognition[END_REF], [START_REF] Otberdout | Deep covariance descriptors for facial expression recognition[END_REF]. To account for both missing data and tracking errors we apply a curve fitting algorithm (as "smoothing" on Bezier curves [START_REF] Gousenbourger | Piecewise-bézier c1 smoothing on manifolds with application to wind field estimation[END_REF]) on the trajectories of the Gram matrices in the manifold. More precisely, the trajectories in the manifold are considered as Bezier curves and can be generalized with the De Casteljau algorithm (see Fig. 1

.b).

To do so, let's consider the positive, semi-definite, of fixed rank matrices S + (d, n). Our goal is to approximate a Bezier curve B to data points d 0 , ..., d t , with t the number of points composing the trajectories. The control points estimated for the curve are chosen such that B is close the the data points and the mean squared acceleration is small. Therefore, the Bezier curves are composed of t Bezier functions. Let p 0 , ..., p t ∈ S + (d, n) denotes the endpoints of each Bezier curve and let

(b - i , b + i ) ∈ S + (d, n) × S + (d, n
) denotes the right and left control points of p i . Fig. 1.b illustrates Bezier curve estimation.

V. KINEMATIC FEATURES AND FISHER VECTOR ENCODING

A. Kinematic Features

Because videos of interviews varied in length, the extracted body curves (of different videos) vary in length. This variation in the obtained curves' lengths could introduce distortions in feature extraction. To overcome this problem, we apply a cubic spline interpolation to the obtained C(t) curves. This step results in smoother and shorter curves of fixed length. We set empirically the new length of the curve given by spline interpolation to 5000 samples for body curves. Then, we compute the velocities by computing the tangent vectors between successive frames. These can be obtained using the logarithm map with respect to the Riemannian metric of S + (d, n) introduced in [START_REF] Massart | Quotient geometry with simple geodesics for the manifold of fixed-rank positive-semidefinite matrices[END_REF]. More specifically, given two successive Gram matrices

C(t) = G t = Z t Z T t and C(t+1) = G t+1 = Z t+1 Z T
t+1 at time t and t+1, respectively, the velocity is a vector at the tangent space of C(t) and is given by,

V C (t) = log C(t) (C(t + 1)) = Z t+1 Q * -Z t , (2)
where Q * is the optimal 2 × 2 rotation aligning Z t+1 and Z t which can be efficiently computed using SVD. By doing so, the velocities will lie on different tangent spaces, i.e., each of them belongs to the tangent space of the current frame.

As a solution, we use the parallel transport provided in the Manopt Toolbox [START_REF] Boumal | Manopt, a matlab toolbox for optimization on manifolds[END_REF] in order to transport all the velocities to the tangent space of a common reference frame. After that, we compute the acceleration A C (t) = ∂ 2 C(t) ∂t 2 from the transported velocities after reducing their dimensions (see Fig. 1.c). Finally, body shapes, velocities, and accelerations are concatenated as:

K C (t) = [C(t); V C (t); A C (t)] , (3)

B. Fisher Vector Encoding

Our goal is to obtain a single vector representation from the kinematic curves K C (t) for depression severity assessment. Following [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF], we used the Fisher Vector representation using a Gaussian mixture model (GMM) distributions [START_REF] Zivkovic | Improved adaptive Gaussian mixture model for background subtraction[END_REF] (see Fig. 1.d). Assuming that the observations of a single kinematic curve are statistically independent, a GMM with c components is computed for each kinematic curve by optimizing the maximum likelihood (ML) criterion of the observations to the c Gaussian distributions. To encode the estimated Gaussian distributions in a single vector representation, we use the

VI. RESULTS

We seek to discriminate three levels of depression severity (moderate to severe, mild, and remission) from body movement dynamics. To do so, we used leave-One-Subject-Out cross validation scheme. Performance was quantified using the mean accuracy over the three levels of depression severity. A multi-class SVM was used to classify the encoded body movement dynamics into three levels of depression severity.

To evaluate the reliability and effectiveness of curve fitting, we first compared the overall classification results with and without curve fitting. As shown in Table II, Gram with curve fitting performed better than without curve fitting. The obtained results can be explained by the fact that curve fitting smooths the curves without modifying the velocity and the acceleration of body movement.

Table . III shows the confusion matrix for depression severity assessment using body movement dynamics. Misclassification was more common between adjacent categories (e.g., Mild and Severe) than between distant categories (e.g., Remission and Severe). Consistent with previous findings, using the dynamics of facial and head movements [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF], [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF], highest accuracy was found for Severe depression (68.57%) compared to Mild (56%) and Remission (37.93%). The assessment of the unique contribution of different body joints movement for depression severity assessment is work in progress.

We evaluated the possible interpretability of the proposed kinematic features (i.e., K C (t) in Eq. 2) for depression severity assessment. To do so, we compute the l2-norm of velocity and acceleration intensities for body (i.e., V C (t) and A C (t)) curves for each video. Since each video is analyzed independently, we compute the histograms of the velocity and acceleration intensities over 10 samples (videos) from each level of depression severity. This results in histograms of 50000 velocity and acceleration intensities for each depression level. Fig. 4 shows the histograms of body velocity and acceleration intensities. The level of depression severity is proportional to the acceleration intensity of body movement. Acceleration decreases as participants improve from severe to mild and remission. Velocity on the other hand failed to vary systematically with change in depression severity. Previous automatic based studies have found that depression severity is marked by reductions in general facial expressiveness [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF]. Velocity and acceleration of facial expressiveness both increased as participants improved from severe to mild and then to remission [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF]. On the other hand, our findings suggest that body movement is highest during severe depression and decrease during mild to remission. These results highlights symptoms of psychomotor agitation (e.g., fast or repetitive body movements) and are consistent with the observation that psychomotor agitation in body movements (e.g., unable to sit still, fidgeting, hand-wringing) in depression lessens as severity decreases (i.e., participants become less agitated when remitted). While encouraging, our results should be interpreted with caution given the unbalanced distribution of depression severity scores. Future research, with a larger sample, is needed to replicate our findings on the relation between depression severity and the dynamics of body movement. Finally, the next step will be to compare the relative contribution of each body joint for depression severity assessment and related dynamics changes during the course of treatment.

Overall the obtained results show that depression severity can be assessed from body movement dynamics. Because other behavioral markers (e.g., face and head movement dynamics) also showed discriminating power for depression severity assessment [START_REF] Dibeklioglu | Dynamic multimodal measurement of depression severity using deep autoencoding[END_REF], [START_REF] Kacem | Detecting depression severity by interpretable representations of motion dynamics[END_REF], the combination of different modalities (i.e., face and head) with body movement dynamics is work under progress.

VII. CONCLUSION

We proposed a space-time geometric representation of body movement dynamics to measure depression severity. The proposed approach model the temporal evolution of joint points configurations as parameterized trajectories on the positive semidefinite matrices of fixed-rank. The results of proposed approach obtained in adults participants with history of chronic depression show promising results.

  1.a). More specifically, let β G : I → S + (d, n) defines a curve (I denotes the time domain, e.g., [0, 1]) to model the spatiotemporal evolution of elements on S + (d, n). Given a sequence

Fig. 1 .

 1 Fig. 1. Overview of the proposed approach, (a) Automatic landmark detection for each frame of the video and Gram matrix computation, (b) Building trajectories on S + (2, n) and Bezier curve smoothing, (c) Kinematic features extraction, (d) Depression severity classification.

Fig. 2 .

 2 Fig. 2. Face-to-face clinical interview setup.

Fig. 3 .

 3 Fig. 3. Overview of OpenPose Targeted Body Points [10].

Fig. 4 .

 4 Fig. 4. Histograms of velocity and acceleration intensities for body movements.

  C (t) separately. The number of Gaussian distributions c are chosen by a a leave-one-subject-out crossvalidation and are set to 9 and the fisher vectors dimension is12 × 2 × 2 × 2 × 9 = 864.The set of encoded dynamic features are then fed to a multiclass SVM with a Gaussian kernel to classify the extracted body movement dynamics into different depression severity levels (see Fig.1.d).

		TABLE II
	OVERALL DEPRESSION SEVERITY ASSESSMENT WITH AND WITHOUT
		CURVE FITTING
		Method		Accuracy
	Gram without curve fitting (89 sequences)	42.89
	Gram with curve fitting (89 sequences)		54.17
		TABLE III
	DEPRESSION SEVERITY ASSESSMENT RESULTS
		Remission Mild	Severe
	Remission	37.93	12	11.42
	Mild	17.24	56	20
	Severe	44.82	32	68.57
	improved fisher vector encoding, which is suitable for large-
	scale classification problems [25]. This step is performed for
	kinematic curves K

VIII. ACKNOWLEDGEMENTS Research reported in this publication was supported in part by the National Institute of Nursing Research of the National Institutes of Health under Award Number R21NR016510, and the National Institutes of Health under awards MH096951 and MH51435. The content is solely the responsibility of the authors and does not necessarily represent the official views of the sponsors.


A Klimchik
email: a.klimchik@innopolis.ru

Anatol Pashkevich
email: anatol.pashkevich@imt-atlantique.fr

Damien Chablat
email: damien.chablat@cnrs.fr

Stiffness modeling of NAVARO

Keywords: parallel robots, variable actuation transmission, stiffness analysis, MSA technique

published or not. The documents may come

INTRODUCTION

For most industrial robots, the actuator transmission system is one of the principal components defining the stiffness properties of an entire manipulator. For conventional serial manipulators, the elasticity of the transmission can be easily taken into account by introducing equivalent one-dimensional springs in actuated joints [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF], Klimchik and Pashkevich, 2017, Klimchik et al., 2017b, Klimchik et al., 2017a). A similar approach is also used for many parallel manipulators (GS platform, etc.) [START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF], where the actuator elasticity is modeled by onedimensional springs integrated into the legs connected moving platform and the base. However, recent developments in parallel robotics produced new architectures for which traditional approach operating with onedimensional equivalent springs can be hardly applied. In such manipulators, the moving platform is actuated via multi-link chains whose elasticity is multi-dimensional and configuration-dependent. This paper presents a stiffness analysis of one of such systems used in parallel manipulator NAVARO II and gives a general framework for elastostatic modeling of similar robots.

At present, there exist three main techniques for the stiffness modelling (Pashkevich et al., 2011b), which are the Finite Element Analysis (FEA) [START_REF] Yan | Stiffness analysis of parallelogram-type parallel manipulators using a strain energy method[END_REF][START_REF] Wang | Finite Element Analysis and Comparison of Two Hybrid Robots-the Tricept and the TriVariant[END_REF][START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF], the Matrix Structural Analysis (MSA) [START_REF] Nagai | A systematic approach to stiffness analysis of parallel mechanisms[END_REF][START_REF] Marie | Elasto-geometrical modeling and calibration of robot manipulators: Application to machining and forming applications[END_REF][START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF] and the Virtual Joint Modelling (VJM) [START_REF] Klimchik | Stiffness matrix of manipulators with passive joints: computational aspects[END_REF], Quennouelle and Gosselin, 2008[START_REF] Gosselin | Stiffness analysis of parallel mechanisms using a lumped model[END_REF], Pashkevich et al., 2011a). The most accurate but computationally expansive is the FEA [START_REF] Yan | Stiffness analysis of parallelogram-type parallel manipulators using a strain energy method[END_REF], while the MSA used in this paper is considered as a compromise technique since it operates with rather large flexible elements connected by the actuated and passive (Klimchik et al., 2019b).

The NAVARO II robot, whose transmission system is studied here, is a novel variable actuation mechanism based on modified 3-RPR parallel kinematics where the second revolute joint of each leg is replaced by a scissor [START_REF] Chablat | NaVARo II: A New Parallel Robot With Eight Actuation Modes[END_REF]. Such modification gives a number of benefits and allows excluding kinematic parallelograms that are traditionally used in similar manipulators. To obtain better spatial rigidity, the leg mechanism is constructed by placing the scissors in an orthogonal plane to the displacement.

The actuation system of NAVARO II can be considered as a double clutch and contains an electric motor and two electromechanical clutches connecting the motor shaft either to the first revolute joint or to the shaft of the prismatic joint. Two position sensors give the leg angular position relative to the base and the current length of the prismatic joint. Depending on the clutches' state, such kinematics allows to implement four actuation schemes where: (#0) the leg is not controlled and may freely change its length and orientation; (#1) the leg length may be arbitrary but its orientation is controlled; (#2) the leg length is controlled while its orientation may be arbitrary; (#3) both the leg length and orientation are controlled in order to execute synchronized translation/rotation. However, only the first and second schemes are studied in this paper because of their practical significance. It is clear that actuation scheme #0 is passive and cannot be used to control the mobile platform location.

FUNDAMENTALS OF MSA MODELLING FOR ROBOTIC MANIPULATORS

For the stiffness modeling of NAVARO II manipulator, which includes numerous closed-loops and passive connections, it is reasonable to use the MSA-based technique. Main theoretical results in this area suitable for the analysis of under-constrained, fully actuated and over-constrained serial and parallel manipulators were presented in our previous works (Klimchik et al., 2018, Klimchik et al., 2019a, Klimchik et al., 2019b). Here, let us present a summary of this approach that will be further applied to the stiffness analysis of the NAVARO II parallel manipulator.

Modeling of manipulator links and platform

In the frame of the MSA technique, the manipulator is presented as a set of nodes connected by flexible or rigid links. The flexible link is described by the following linear matrix equation relating displacements at the end-nodes

( ) ( ) 11 12 ( ) ( ) 21 22 1 12 2 ij ij ii ij ij jj        =            Wt KK Wt KK (
W K 0 0 K t W t 0 K 0 K W t 0 0 K K t W K K K K K K (2)
where ,,

i j k
   t t t are the deflections at the leg clamping points, e t is the deflection at the end-effector reference point, ,

i j k
W W W are the wrenches at the leg clamping points, e W is the total wrench applied to the end-effector. It should be mentioned that all stiffness matrices in eqs. ( 1) and (2) should be presented in a global coordinate system. Relevant numerical values of the stiffness matrix elements can be obtained either from the approximation of the links by regular beams or from the FEA-based experiments in the virtual environment using link's CAD models [START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF].

If the link flexibility is negligible, the stiffness model (1) should be replaced by the "rigidity constraint"

3 3 3 ) 3 6 ( 6 [ ] T j i i j      =   d I D 0I 33
[] ij  d denotes the 3 3  skew-symmetric matrix derived from the vector () ij d describing link geometry. In addition, it is necessary to take into account the force equilibrium that here can be presented in the following form

() 6 6 6 1 ij T i j    =    W 0 W ID (5)
Totally, this yields 12 scalar equations describing relations between the displacements/wrenches, similarly to (1).

MSA models of manipulator's joints

To take into account connections between the adjacent links by means of joints, the stiffness model must include relations describing two principal rules of structural mechanics: (a) displacement compatibility, (b) force equilibrium. For the rigid joint, the displacement compatibility is expressed as

  6 6 6 6 6 1 6 12 i j       -  =    t 0 t II (6)
and the force equilibrium is presented as

  6 6 6 6 6 1 6 12 i j      =   W I I 0 W (7)
For the passive joint, the displacement compatibility condition contains only 6 r  equations

1 12 ** rr i r r j ij ij      -  =     Λ 0 t Λ t (8)
where

  1 * 6 ,..., r r T r  = Λ u u is a rectangular matrix of the rank r that is formed the vectors 1 2 6
,. , .., u u u composing the orthonormal basis associated with the orientation of passive joint(s) in such way that the unit vectors 1 ,..., r uu describe the directions of the rigid connection and the unit vector(s) 16 ,..., r + uu correspond to the passive connection allowing free relative motions. In this case, the force equilibrium gives the following relations For the elastic joint, the deflection compatibility condition (8) remains the same. However, the force equilibrium condition must be replaced by

              =      -        0 0 0 Λ Λ Λ 0 t II t W KK Λ W W ( 10
)
where the matrix  stiffness matrix describing elastic properties of the joint, 0 ij W is preloading in the corresponding springs.

Connection to the base and loadings

The rigid connection of the link to the robot base can be presented as a special case of the rigid joint with zero deflection j = t0 . This simplifies the deflection compatibility constraint (6) down to

  6 6 6 1 j     =  I0 t (11)
The passive connection of the link to the robot base can be presented as a special case of the passive joint with zero deflection in the non-passive directions * r ij j   = Λ t 0 . This reduces the deflection compatibility constraint (8) and force equilibrium to

        =         =   -       Λ 0 W W t t K Λ Λ Λ (13)
To take into account the external loading e W , the global stiffness model must be extended by the linear equation ensuring the force equilibrium, i.e.

... A similar technique is used to all nodes where the external wrenches are applied.

Aggregation of MSA model components

After aggregation of all equation describing the manipulator links and their connections, the global stiffness model of the manipulator can be presented in the form

(1) (1) (1) W0 (2) (2) ( t 2 W Δ 0 Δt ) {} {} ... ... ... i i          =                t A A b W A A b ( 15
)
where

           =                   S K K W b E F C W W t ED t F (16)
where the matrices are obtained using relevant equations describing the links and connections, while the vectors where the matrix inversion usually exists (if the manipulator does not include redundant passive joints).

NAVARO II MANIPULATOR ARCHITECTURE

Particularities of NAVARO II architecture

NAVARO II is a special case of the 3-RPR manipulator (Fig. 2) where the translational joint is implemented by means of double clutch pantograph-type transmission presented in Fig. 3. In more details, the architecture of this manipulator is presented in Fig. 1. Due to the double-clutch system, the transmission allows to implement three following actuation modes:

Mode #1: Clutch 1 is active while clutch 2 is passive. The first leg axis (green) is driven by the rotation of the motor shaft. In this case, the angle  is active while  is passive.

Mode #2: Clutch 2 is active while clutch 1 is passive. The first leg joint is free, but the rotation of the motor leads to a displacement of the slider, which activates the scissors.

In this case, the  is passive and  is active.

Mode #3: Both clutches 1 and 2 are active. Both joints cause a synchronized rotation and translation motion. The end of the leg will make a spiral motion. In this paper, we will focus on the stiffness analysis of the NAVARO II manipulator and its transmission for the actuation modes #1 and #2, which are the most interesting for practical application.

Kinematics of NAVARO II manipulator

To estimate the manipulator stiffness behavior within its workspace it is required to obtain relations allowing to find the transmission variables i  and i  corresponding to the given end-effector location, i.e. to solve the inverse kinematic problem. For this manipulator, for any given end-effector position O and platform orientation  , the positions of the platform connection points i B can be computed as 2 ( 2) , 1, 2 0 6 3 0 , 3

iz i r i    -+ -+ =    = -         B O R (18)
where r is the distance from the platform center to the connection points (the same for all points), and

( )

.

i yi yi xi xi i i i i   -- = - = = AA B A B B (19)
It should be noted that the inverse kinematic solution exists for reachable mobile platform location only, which is specified by minimum and maximum pantograph lengths, i.e.

min max i     ( 20
)
The values i  will be used further to compute orientations of the pantograph internal links.

Kinematics of NAVARO II transmission system

As it was mentioned before, the NAVARO II controller may utilize two actuation mode (types of transmission): with active and passive pantographs (Fig. 4). Clearly, this effects on the transmission stiffness properties while from the kinematical point of view these cases are identical.

The latter allows us to compute the orientation of each link of the pantograph transmission and to transform the link stiffness matrices into the global coordinate system.

MSA MODELS MANIPULATOR COMPONENTS

To apply the MSA technique for the stiffness modeling of NAVARO II, let us split the manipulator into four parts: three kinematically identical legs and the mobile platform, which is actuated by means of either passive or active pantographs. The MSA representations for these principal components are shown in Fig. 5. Relevant models are composed of a number of flexible and rigid links (i,j) connected in a special way by passive, rigid and elastic joints <i,j>. The model components are also listed in Table 1.

MSA models of pantograph-type transmission system

In order to obtain stiffness model of entire NAVARO II mechanism let us first derive the stiffness model for separate legs represented via passive or active pantographs. AS follows from Figure 5a,b the leg is composed of a set of eight flexible links (3,5), (4,6), (7,9), (8,10), (11,13), (12,14), (15,17), (16,18) with rigid connections <5,8>, <6,7>, <13,16>, <14,15>, elastic connections <1,1'> or <0,1> for active and passive pantographs and passive connections between the nodes <1,3>, <2,4>, <5,6>, <9,11>, <10,12>, <13,14>, <17,e>, <18,19>, <e,19> complemented by the connections <0,1> or <1,2> for active and passive pantographs. Active pantograph also contains a rigid link (1',2). Such presentation allows us to describe the elastic properties of pantographs by the following matrix equation

Structural components

Passive joints Elastic joints Rigid joint Elastic links Rigid links

Active pantograph <0,1>, <1,3>, <2,4>, <5,6>, <9,11>, <10,12>, <13,14>, <17,e>, <18,19>, <e,19> <1,1'> <5,8>, <6,7>, <13,16>, <14,15>

(3,5), (4,6), (7,9), (8,10), (11,13), (12,14), (15,17), (16,18) (1',2)

Passive pantograph <1,2>, <1,3>, <2,4>, <5,6>, <9,11>, <10,12>, <13,14>, <17,e>, <18,19>, <e,19> <0,1> <5,8>, <6,7>, <13,16>, <14,15>

(3,5), (4,6), (7,9), (8,10), (11,13), (12,14), (15,17)

agr leg agr     =   -      W 0 K t I (23)
where agr t and agr W are the aggregated displacement and wrenches defined at the nodes from 3 to 18, leg K aggregates stiffness matrices of links (3,5), (4,6), (7,9), (8,10), (11,13), (12,14), (15,17), (16,18) presented in the form (1) into a single matrix of a size 96 96  allowing to describe the manipulator leg elasticity.

For the considered manipulator, the MSA based model includes three types of passive joints and two types of elastic joints with passive directions around z and y-axis and translational joint in z directions, which are described by the following matrices

    * * *
1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 ; 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 ; 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1

rp it    ==     ΛΛ ( 26
)
where the matrix * r yi Λ depends on the orientation of the ith leg of the manipulator and are different for all chains in nonsymmetrical configurations. Using this notation, the constraints imposed by the passive joints <1,3>, <2,4>, <5,6>, <9,11>, <10,12>, <13,14>, <17,e>, <18,19> can be presented as

              -      =     -      0 0 0 0 0 0 0 t ΛΛ t II W K W Λ Λ Λ 0 K (31)
where the coefficient e  K defines rotation joint stiffness. For both pantographs, it is required to take into account that pantograph contains only 4 component by introducing rigid joint constraints via equations ( 6) and ( 7) for the joints <5,8>, <6,7>, <13,16>, <14,15>. Due to particularities of the abovedescribed assumptions, for both active and passive pantographs, boundary constraints are introduced in the system for the node <0> via rigid constraint (11) and external force introduces via eq.( 14) for the node <e>. Combining all above-mentioned equations it is possible to write stiffness model for active and passive pantographs in the form (15), where for passive pantograph the size of the main matrix has the size of 252 288  and for active pantograph has the size of 264 294  . It is obvious that the system of equation is over-constrained, so to find the desired Cartesian stiffness matrix using expression(17) pseudo inversion should be used.

MSA models of NAVARO II parallel manipulator

Using the stiffness model of the transmission system, it is possible to obtain the model of the entire manipulator can be presented as a non-rigid mobile platform connected via passive joints to equivalent models of pantograph-type transmissions fixed at another end. In this case to describe elasticity of mobile platform one can use eq. ( 2) and for legs flexibility reduced form of eq. ( 1), where the first node is eliminated because of the rigid connection to the base. The latter allows us to use only 6 6  stiffness matrices to describe elasticities of the pantograph-type transmission system. Passive joints for the connections <2a,3a>, <2b,3b>, <2c,3c> describes using expression (29). External loading applied to the platform end-effector <ep> is introduced via eq. ( 14) Combining all above-mentioned equations into a single system of equations one can get a system of order 84 (with the matrix of the size 84 84  ) in the form ( 15) and ( 16) that can be solved with respect to the reaction of the end effector node <ep> on the external wrench applied to it using eq. ( 17). This stiffness matrix can be computed for any reachable end-effector location and used for stiffness analysis of manipulator behavior within its workspace.

CONCLUSIONS

This paper deals with the stiffness modeling of a new type of parallel manipulator NAVARO II with a pantograph-type transmission system. The desired models were obtained using the enhanced matrix structural analysis (MSA) approach that is able to analyze the under-actuated and over-constrained structures with numerous passive joints. The MSA technique was applied to active and passive pantograph-type mechanisms for which the stiffness models operate with the matrices of sizes 252 288  and 264 294  respectively. MSA-based model for the entire manipulator with equivalent models for the pantograph-type transmission legs stiffness models operates with the matrix of sizes 84 84  . These models will be used for manipulator parametric optimization and stiffness analysis of NAVARO II manipulator with different types of actuation modes.

Figure 1 .

 1 Figure 1. Parallel planar manipulator NaVaRo.

  be noted that all passive connections should be treated separately (in contrast to rigid ones).

  elastic connection of the link to the robot base can be presented as a special case of the elastic joint with zero deflection in the non-elastic directions *

  can be rewritten in the form (for three adjacent links)

  wrench and displacement of the i th node, the matrices () above using expressions presented in the previous sections. After rearranging the matrix rows and introducing relevant definitions for the blocks, the global stiffness model can be rewritten as

  internal nodes and the end-effector node, where the external wrench e W is applied. The latter allows us to obtain, the desired Cartesian stiffness matrix

Figure 2 .

 2 Figure 2. Kinematic model of NAVARO II manipulator.

Figure 3 .

 3 Figure 3. Kinematic model of NAVARO II transmission.

Figure 4 .

 4 Figure 4. Equivalent models of NAVARO II transmission. For both passive and active pantographs, the inverse kinematics problem (giving the orientation of the transmission links) can be solved using the relation between the length of the leg  and the length of pantograph links l 2 cos l   = (21) That gives the following expressions for the links orientations acos( / 2 ), /2 up down up l     == -(22)

  Figure 5. MSA-based representation of the NAVARO II manipulator.

Table 1

 1 Components of the MSA-based model of NAVARO II manipulator and its legs

	Mechanical
	structure

ACKNOWLEDGMENTS

The work presented in this paper was supported by the grant of the Russian Science Foundation №17-19-01740.