domain: math stat

text:
Gabriela Ciuperca

Matúš Maciak

Detection of Similar Successive Groups in a Model with Diverging Number of Variable Groups

Keywords: Adaptive penalty, Different successive groups, Diverging-dimensional group model, Fused group. Subject Classifications: 62F12, 62F35, 62J07

In this paper, a linear model with grouped explanatory variables is considered. The idea is to perform an automatic detection of different successive groups of the unknown coefficients under the assumption that the number of groups is of the same order as the sample size. The standard least squares loss function and the quantile loss function are both used together with the fused and adaptive fused penalty to simultaneously estimate and group the unknown parameters. The proper convergence rate is given for the obtained estimators and the upper bound for the number of different successive group is derived. A simulation study is used to compare the empirical performance of the proposed fused and adaptive fused estimators and a real application on the air quality data demonstrates the practical applicability of the proposed methods.

INTRODUCTION

The idea of this paper is to automatically detect different successive groups of unknown coefficients of some explanatory variables in a multivariate linear model. The number of groups is supposed to be of the same order as the number of observations. For a given loss function, the fused type penalties allow this automatic detection of these successive groups of the unknown coefficients. Depending on the assumptions imposed on the model errors, two modeling frameworks are considered: either the standard least squares loss function is used or the robust quantile loss function is considered instead. Moreover, for each framework, two fused group penalties are proposed: firstly, the fused-type penalty which is later used to construct the adaptive fused penalty leading to a more accurate selection of different successive groups. For each of the two estimators the convergence rates are provided and the upper bound for the number of the successive groups is derived.

Address correspondence to G. Ciuperca, Université de Lyon, Université Claude Bernard Lyon 1, CNRS, UMR 5208, Institut Camille Jordan, Bat. Braconnier, 43, blvd du 11 novembre 1918, F -69622 Villeurbanne Cedex, France; E-mail: Gabriela.Ciuperca@univ-lyon1.fr 1 In order to highlight the novelty of our work, we firstly make the state of the art regarding the proposed fused method with the automatic detection of the grouped explanatory variables in the multivariate linear model. Let g ∈ N denote the number of variable groups and let n ∈ N be the total number of the available observations. The fused quantile method for a particular case of nongrouped variables with the quantile level τ = 0.5 was already considered in [START_REF] Liu | Fused LASSO penalized least absolute deviation estimator for high dimensional linear regression[END_REF] where the fused LASSO penalized least absolute deviation (LAD) estimator in a high-dimensional linear model is discussed and the proper convergence rate of the obtained estimator is derived together with a linearized alternating directional method for finding the numerical solution.

The quantile linear model with a finite number of non-grouped explanatory variables is investigated by [START_REF] Jiang | Interquantile shrinkage in regression models[END_REF] and [START_REF] Jiang | Interquantile shrinkage and variable selection in quantile regression[END_REF] by utilizing the adaptive fused penalization. In [START_REF] Jiang | Interquantile shrinkage in regression models[END_REF], the oracle property for the difference in the estimated coefficients for two different quantile levels is proved. More precisely, an automatically detection of the unchanged quantile slope coefficients across various quantile levels is discussed. In [START_REF] Jiang | Interquantile shrinkage and variable selection in quantile regression[END_REF], the adapted fused method is used to automatically select the explanatory variables and to identify their successive differences at the neighhoring quantile levels. For a linear quantile regression with g groups of explanatory variables, [START_REF] Ciuperca | Adaptive Fused LASSO in Grouped Quantile Regression[END_REF] shows the oracle property for the adaptive fused estimator when g = O(n c ), for 0 ≤ c < 1.

If the model errors satisfy the classical conditions (i.e., zero mean and bounded variance) then the least squares (LS) loss function is more appropriate: in such case, the high-dimensional linear model with the automatic selection of the corresponding groups of the explanatory variables with the adaptive LASSO penalty is considered by [START_REF] Wei | Consistent group selection in high-dimensional linear model[END_REF] for the Gaussian errors when the number of groups is much larger than the sample size (g n) and by [START_REF] Zhang | On the oracle property of adaptive group lasso in high-dimensional linear models[END_REF] for non Gaussian errors. These results are further elaborated in [START_REF] Wang | Variable selection in quantile regression[END_REF] for a generalized linear model when g = n c , with 0 < c < 1. The automatic selection of the grouped variables is also considered in [START_REF] Guo | Model selection and estimation in high dimensional regression models with group SCAD[END_REF] where the SCAD penalty is utilized under the assumption that the number of groups can grow at a certain polynomial rate with n. A combination of the L 1 and L 2 norms under the Gaussian model errors is investigated in [START_REF] Campbell | Within group variable selection through the Exclusive Lasso[END_REF], where the authors propose a structured variable selection in order to select at least one variable from each group. To our best knowledge, the only papers considering the fused penalty with the main focus on the selection of variable groups, there is a paper of [START_REF] Li | Linearized alternating direction method of multipliers for sparse group and fused LASSO models[END_REF] where the LS loss function is penalized with the fused LASSO penalty, where the L 1 norm si considered for the magnitudes of the parameters and also for the successive differences of between the estimated coefficients.

In the present paper, the penalty is of the fused type, that is, it is built against the L q,1 norm (with q ≥ 1) of the difference between two successive groups of parameters, while in the mentioned just before papers, the norm in the penalty is L 2,1 or L 1,1 of each parameter group, the goal being to automatically select the significant coefficient groups and not the identical successive coefficient groups. In a model, it can have successive vectors of non-zero parameters that are not different. A practical example is given in Section 4 of the present paper on the influence of the groups of meteorological variables measured every hour, on the daily maximum benzen concentration. It is this type of automatic detection that interests us in the present work. Whether for the quantile or LS methods, particular cases of L q,1 penalties of the difference between two successive groups of parameter vectors were considered within the change-points automatic detection framework in linear model. Except that, in the literature, for linear models with change-points the statistical model is different from that considered in this paper, because the parameter number of the model was constant. For the LS loss function, we have in [START_REF] Zhang | Multiple change-points estimation in linear regression models via sparse group lasso[END_REF] the sum of the penalties L 2,1 and L 1,1 , while in [START_REF] Qian | Shrinkage estimation of regression models with multiple structural changes[END_REF] the penalty is L 2,1 . For the quantile loss function, [START_REF] Ciuperca | Change-point detection in a linear model by adaptive fused quantile method[END_REF] consider the L 2,1 penalty. The paper is organized as follows. In Section 2 we introduce the model, assumptions and general notation. In Section 3, fused and adaptive fused group estimators for LS and quantile loss function are defined and asymptotically studied. In Section 4 we present a simulation study on the proposed estimators ans an application on real data. The proofs of the results in Section 3 are given in Section 5.

MODEL

In this section we state the model definition and some general assumptions imposed on the model design. Let us start, however, with some notation which will be used throughout the paper: All limits in are taken with respect to n → ∞; All vectors are columns and matrices and vectors are denoted with a bold face; For some matrix A we denote its transpose as A and for a set A, we denote by |A| its cardinality and by A its complement; Expressions µ min (.) and µ max (.) are used to refer to the smallest and largest eigenvalue of some positive definite matrix and for

x = (x 1 , • • • , x p ) ∈ R p being some p dimensional vector x q = |x 1 | q + • • • + |x p | q 1/q denotes its L q norm while x ∞ = max(|x 1 |, • • • , |x p |) stands for the maximum norm. If, in addition, x = (x 1 , • • • , x g
) is a vector split into g subvectors, then g j=1 x j q defines for the L q,1 norm of x.

Moreover, β j = (β j,1 , • • • , β j,p ) ∈ R p stands for the corresponding group specific vector of the dimension p ∈ N, for any j ∈ {1, . . . , g}, where g ∈ N is the number of the successive groups. Last but not least, C denotes some positive generic constant not depending on n which may take different values in different formulas throughout the paper.

In the present paper, we consider a multivariate linear model with g groups of explanatory variables. The number of groups g ∈ N depends on the sample size n ∈ N, g being known, such that g ≤ n/p, while the number of explanatory variables in each group is fixed and does not depend on n. Without reducing any generality, it is assumed that each group of the explanatory variables contains the same number of variables, p ∈ N. Thus, the overall number of all parameters in the regression model is r n ≤ n.

Let us consider the following linear model with the grouped explanatory variables

Y i = g j=1 X i,j β j + ε i = X i β g + ε i , i = 1, • • • , n, (2.1) with β g ≡ (β 1 , • • • , β g ) ∈ R rn
, where β j ∈ R p is the vector of parameters for the group j ∈ {1, . . . , g}. For each observation i ∈ {1, . . . , n}, the vector X i ∈ R rn contains the explanatory variables X i,j ∈ R p from all groups. These group specific explanatory variables are assumed to be deterministic, for any j = 1, . . . , g and i = 1, . . . , n. The error terms {ε i } i are assumed to be independent and the response variable is denoted as Y i ., ε i . The true (unknown) vector of parameters is β 0 = (β 0 1 , . . . , β 0 g ) . For p = 1, the model with ungrouped explanatory variables is obtained. Note that the order of appearance of the groups in the model in (2.1) is important and some natural ordering is required.

Given the data {(Y i , X i ); i = 1, . . . , n} we would like to automatically determine, using the fused method, whether two successive groups of the explanatory variables have the same influence on the response or not while, at the same time, quantifying the corresponding effect magnitudes. In addition to the example on the air pollution in Section 4, a nice demonstration of the practical applicability of the proposed estimation method can be also seen in the very recent work of [START_REF] Zhou | Modeling Disease Progression via Fused Sparse Group Lasso[END_REF], where the fused group method allows for capturing the temporal smoothness of the predictive biomarkers on the cognitive scores in the progression of the Alzheimer's disease. To achieve the sparsity property between two successive groups of the explanatory variables (in a sense that the corresponding vectors of estimated parameters for two successive groups are mostly the same), the fused and adaptive fused group estimators are proposed and studied with two loss functions: the standard least squares and the quantile check function.

The asymptotic behavior of the group specific estimators for the fused and the adaptive fused method with n ≥ gp are investigated for n → ∞ where, in addition, a deterministic sequence

(b n ) n∈N is needed, such that b n → 0, n 1/2 b n → ∞. (2.2) Example of such sequence which satisfies (2.2) is b n = n -1 log n 1/2 .
Unlike [START_REF] Ciuperca | Adaptive Fused LASSO in Grouped Quantile Regression[END_REF], where the number of groups is either fixed or it is of the order n c , with 0 < c < 1, the model in (2.1) assumes that the number of the successive groups may be of the same order as the sample size. A similar model is also considered in [START_REF] Ciuperca | Change-point detection in a linear model by adaptive fused quantile method[END_REF] where the change-point detection and estimation is performed in the quantile model with fused type penalty, however, for the unknown vector of parameters with the dimension p, not depending on n. The same model is also considered in [START_REF] Leonardi | Computationally efficient change point detection for highdimensional regression[END_REF] where the change-point locations are detected by utilizing the LS loss function with the LASSO type penalty.

Assumptions

The following regularity assumptions imposed on the model design are needed. The assumptions required for the model errors will be presented in Subsection 3.1 for the quantile framework and in Subsection 3.2 for the LS framework.

(A1) max 1 i n X i ∞ ≤ C 0 , for some constant C 0 > 0. (A2) There exist two positive constants, 0 < m 0 ≤ M 0 < ∞, such that m 0 ≤ µ min (n -1 n i=1 X i X i ) ≤ µ max (n -1 n i=1 X i X i ) ≤ M 0 .
Assumption (A1) is considered, for instance, in [START_REF] Leonardi | Computationally efficient change point detection for highdimensional regression[END_REF] for the highdimensional regression model and, also, by [START_REF] He | Regularized quantile regression under heterogeneous sparsity with application to quantitative genetic traits[END_REF] for the penalized quantile regression. Assumption (A2) is standard in the linear model to ensure the parameter identifiability (see, for example, [START_REF] Zhang | On the oracle property of adaptive group lasso in high-dimensional linear models[END_REF][START_REF] Ciuperca | Adaptive group LASSO selection in quantile models[END_REF][START_REF] Ciuperca | Adaptive Fused LASSO in Grouped Quantile Regression[END_REF][START_REF] Wu | Variable selection in quantile regression[END_REF].

ESTIMATION METHODS

In this section, two estimation frameworks are presented: the automatic detection and estimation of the successive groups of the explanatory variables is considered under two different model error assumptions. For each framework, the asymptotic properties are investigated. Firstly, the fused group estimator is proposed and, afterwards, the adaptive version of the fused group estimator is defined.

If the model errors {ε i } i i n in (2.1) do not meet the standard conditions for the existence of the first two moments then the robust version needs to be employed, therefore, the quantile estimation technique is appropriate. On the other hand, if the conditions E [ε i ] = 0 and Var [ε i ] < ∞ are satisfied, the penalized LS method is considered. The main results are presented for both scenarios in next two subsections while the proofs are all postponed to Section 5.

Quantile Loss Function

Let the model errors in (2.1) satisfy the following:

(A3) Random errors ε i , for i = 1, . . . , n, are independent and identically distributed (i.i.d.) with the continuous distribution function F , such that F (0) = P[ε ≤ 0] = τ , for some known τ ∈ (0, 1). The corresponding density function f with the nonzero compact support G ⊆ R is supposed to be continuous and strictly positive in a neighborhood of zero. Moreover, the first derivative of f is bounded in a neighborhood of zero.

Assumption (A3) on the errors is standard for the quantile regression models when the number of parameters depends on the sample size n ∈ N (see, for instance, [START_REF] Ciuperca | Adaptive group LASSO selection in quantile models[END_REF]Wu and[START_REF] Wu | Variable selection in quantile regression[END_REF]. The standard assumptions E [ε] = 0 and E [ε 2 ] < ∞ are not considered and, therefore, the least squares method is not appropriate. Since P[ε < 0] = τ , we can consider the quantile method with the fixed quantile level τ ∈ (0, 1), with the corresponding quantile check function ρ τ (u) = u(τ -1 1 {u<0} ), for u ∈ R. Thus, for the model in (2.1) the following quantile random process is obtained

G n (β g ) ≡ n i=1 ρ τ (Y i -X i β g ), (3.1)
with the group quantile estimator defined as

β g ≡ arg min β g ∈R rn G n (β g ). (3.2)
For the particular case of τ = 0.5 we obtain the median regression, for which the quantile process and the associated estimator (3.2) are reduced to the absolute deviation process and the least absolute deviation estimator respectively. The following Lemma gives the appropriate convergence rate of the group quantile estimator β g .

Lemma 3.1. Under Assumptions (A1), (A2), and (A3) it holds that

β g -β 0 1 = O P (b n ),
where

(b n ) n∈N is the sequence defined in (2.2).
The convergence rate of the group quantile estimator for the number of groups [START_REF] Ciuperca | Adaptive group LASSO selection in quantile models[END_REF] and the convergence rate of β g from (3.2) can not be obtained as a straightforward extension of the situation where g = O(n c ) for 0 ≤ c < 1, when c → 1.

g = O(n) is different from that obtained when g = O(n c ), with 0 ≤ c < 1. Indeed, for 0 ≤ c < 1 the convergence rate of β g is of the order O P (gn -1 ) 1/2 = O P (n (c-1)/2 ) (see Lemma 1 of
In order to preserve the group effect of the explanatory variables and to simultaneously detect the successive groups of identical parameter vectors the L q,1 norm of the consecutive differences β j -β j-1 , for j = 2, • • • , g, is used as a penalty with some q ≥ 1 fixed. Thus, the following quantile process is considered

Q n (β g ) ≡ G n (β g ) + nλ n g j=2 β j -β j-1 q .
(3.3)

For q = 1 the relation is (3.3) gives the process penalized with the standard L 1 norm while for q = 2 the process is penalized by the L 2,1 norm. The positive sequence (λ n ) n∈N plays a role of a tuning parameter, such that it converges to zero as the sample size tends to infinity. An additional condition on (λ n ) n∈N will be given later when formulating the theorems with the main results.

Based on the penalized process in (3.3), the corresponding fused group quantile estimator is obtained as

β g ≡ arg min β g ∈R rn Q n (β g ), (3.4)
where β g = β 1 , . . . , β g . The estimator β g depends on the norm considered in the penalty term of random process in (3.3) and, also, the tuning parameter λ n > 0.

Let us define the set of indexes which form the true different successive groups

B 0 = j ∈ {2, • • • , g}; β 0 j = β 0 j-1 . (3.5)
Since the values of the true parameter vector β 0 are unknown the set B 0 is left unknown too. Therefore, an analogous set is considered with respect to the differences of the estimated parameters of two successive groups as

B n = j ∈ {2, • • • , g}; β j = β j-1 .
It is obvious, that this set is used to provide a reasonable estimate for B 0 .

Remark 3.1. The results obtained in this section are also valid for p = 1, which is, to the authors' best knowledge, the case which has not been previously considered with in any literature. The number of the groups in Ciuperca ( 2017) is of order n c , with 0 ≤ c < 1 and, moreover, in Ciuperca ( 2017), the goal is to select the groups of significant variables simultaneously with the group's inheritance.

The following theorem provides the convergence rate of the fused group quantile estimator defined in (3.4), under the additional assumption that there is only a finite number of the successive groups with different coefficients. For a suitable choice of the tuning parameter this convergence rate is of the same order as the sequence (b n ) and, moreover, it is the same as the one obtained for the group quantile estimator in Lemma 3.1. The convergence rate of β g does not depend on the L q norm considered in the penalty term in (3.3).

Theorem 3.1. Under Assumptions (A1), (A2), and (A3), the condition in

(2.2), if, moreover, |B 0 | < ∞ and λ n b -1 n -→ n→∞ 0, then β g -β 0 1 = O P (b n ).
Examples of such sequences

(λ n ) n∈N , (b n ) n∈N which satisfy (2.2) and λ n b -1 n -→ n→∞ 0 are λ n = n -1 (log n) 1/2 and b n = n -1 log n 1/2 .
Similarly as for the standard LASSO type penalties, the consistent selection of the different successive groups does not occur with the probability converging to 1 and some overfitting is observed. The missclassification error | B n \ B 0 | is used to assess the number of the different successive groups being mistakenly detected as different. The following theorem provides the upper bound for this missclassification error.

Theorem 3.2. Under the same assumptions as in Theorem 3.1, there exists a positive constant C 1 > 0, such that

lim n→∞ P | B n \ B 0 | ≤ C 1 max b n λ n , 1 b n = 1.
Note, that the upper bound in Theorem 3.2 depends on the tuning parameter λ n > 0 and the sequence (b n ) n∈N abd thus, it can be hypothetically unbounded from above. Nevertheless, this result provides the upper bound for the number of elements in B n , more specifically, it gives the upper bound for the number of successive groups of explanatory variables which have different estimated effect on the response variable.

Corollary 3.1. Since |B 0 | < ∞ and B n \ B 0 ≥ | B n | -|B 0 |
with probability one, we can deduce by Theorem 3.2, that lim

n→∞ P | B n | ≤ C max b n λ -1 n , b -1 n = 1. Remark 3.2. For instance, if λ n = n -1 (log n) 5/2 and b n = n -1 log n 1/2
, then the upper bound given by Theorem 3.2 is

| B n | ≤ C max n 1/2 (log n) -2 , n 1/2 (log n) -1/2 = Cn 1/2 (log n) -1/2 ,
which implies that the number of elements contained by B n is much smaller than n 1/2 , however, it can converge to infinity for n → ∞.

To improve the estimation accuracy of B 0 we consider an adaptive penalty constructed on the basis of the estimator in (3.4). Let us consider the random process

∨ Q n (β g ) ≡ G n (β g ) + nλ n g j=2 ω n,j β j -β j-1 q , (3.6)
with the adaptive weights

ω n,j ≡ 1/ max n -1/2 , p k=1 | β j,k -β j-1,k | γ ,
for a fixed constant γ > 0, where β j = ( β j,1 , . . . , β j,p ) . Let us remark that for j ∈ B n we have β j -β j-1 = 0 p . The tuning parameter sequences in relations (3.3) and (3.6) may be different, both with a convergence rate faster than the sequence (b n ) n∈N . Therefore, the choice of n -1/2 in ω n,j is used as deterministic sequence that converges to 0 when β j = β j-1 , however, with the rate faster than b n because of the condition

n 1/2 b n → ∞ in (2.2). Notice that ∨ Q n (β 0 ) ≡ G n (β 0 ) + nλ n g j=2 ω n,j β 0 j -β 0 j-1 q
and the adaptive fused group quantile estimator for β 0 is defined as

∨ β g ≡ arg min β g ∈R rn ∨ Q n (β g ), where ∨ β g = ∨ β 1 , • • • , ∨ β g
. By Theorem 3.1, we have that for all j ∈ B 0 there exists a constant c > 0 such that

lim n→∞ P ω n,j > c | j ∈ B 0 = 1. (3.7)
Therefore, taking into account the relation in (3.7) and the fact that γ > 0 a similar proof to that of Theorem 3.1 can be used to derive the convergence rate of ∨ β g which is the same as for β g .

Theorem 3.3. Under Assumptions (A1), (A2), and (A3), the condition in

(2.2), if |B 0 | < ∞ the for any sequence (λ n ) n∈N such that λ n b -1 n -→ n→∞ 0 it holds that ∨ β g -β 0 1 = O P (b n ).
Considering the adaptive fused group quantile estimator ∨ β g we can also define an updated estimator for the set B 0 as

∨ B n ≡ j ∈ {2, • • • , g}; ∨ β j = ∨ β j-1 ,
which is indeed more appropriate as shown by the next theorem where the upper bound for the cardinality of ∨ B n \ B 0 is proved to be much smaller than the one for B n \ B 0 in Theorem 3.2.

Theorem 3.4. Under the same assumptions as in Theorem 3.3, there exist a positive constant C 2 such that,

lim n→∞ P | ∨ B n \ B 0 | ≤ C 2 max n -1/2 , b γ n max b n λ n , 1 b n = 1.
Remark 3.3. (i) For γ > 1 and the tuning parameter

(λ n ) n∈N such that n -1/2 b n λ -1 n → 0 and b γ+1 n λ -1 n → 0, we obtain that max n -1/2 , b γ n max b n λ -1 n , b -1 n → 0, as n → ∞. The examples of sequences (λ n ) and (b n ) from Remark 3.2 satisfy these conditions. (ii) If 0 < γ ≤ 1 then, max n -1/2 , b γ n = b γ n . In this case we have, b γ n max b n λ -1 n , b -1 n ≥ b γ-1 n
and the sequence on the right-hand side of this inequality converges to infinity for γ < 1 and it is bounded for γ = 1. Thus, in this case, it seems like we should take the value γ = 1 and the same sequences (b n ), (λ n ) as in Remark 3.2.

Comparing now Theorem 3.2 and Theorem 3.4, we can deduce that the adaptive weights ω n,j are responsible for a strong reduce the number of elements in ∨ B n ∩ B 0 , e.i, the false discoveries of different successive groups. This is also later confirmed by the simulation study performed in Section 4.

Least Squares Loss Function

In a standard linear regression model the least squares (LS) objective function is standardly used under the following assumptions imposed on the model errors:

(A4) The error terms (ε i ) 1 i n are i.i.d., such that E [ε] = 0 and Var [ε] < ∞;
We will now focus on the fused and adaptive fused group estimator based on the least squares objective function. In this case, instead of (3.1), an analogous empirical process is considered

L n (β g ) ≡ n i=1 (Y i -X i β g ) 2 ,
(3.8)

with the corresponding estimator given as β g (LS) ≡ arg min

β g ∈R rn L n (β g ),
and the penalized process analogous to (3.3) is

U n (β g ) ≡ L n (β g ) + nλ n g j=2 β j -β j-1 q , (3.9)
with the corresponding fused group LS estimator β g (LS) ≡ arg min

β g ∈R rn U n (β g ).
A similar linear model with non-grouped explanatory variables (p = 1) with the penalty of the form αν

(1) n g j=1 |β j | + (1 -α)ν (1) n j<k |β j -β k |,
for some α ∈ (0, 1], with the LS objective function is also considered in [START_REF] Jang | Some properties of generalized fused lasso and its applications to high dimensional data[END_REF], however, for the situation where g ∈ N is fixed. The corresponding fused estimator allows the selection of groups of predictors that are positively correlated.

As already pointed out in Remark 3.1, the results derived in this section for the LS framework are also novel for a model containing non-grouped variables (p = 1) as the number of groups is allowed to increase with the sample size.

The convergence rates of the proposed estimators β g (LS) and β g (LS) are the same as those obtained for the analogous estimators obtained for the quantile framework in Subsection 3.1. Lemma 3.2. Under Assumptions (A1), (A2), and (A4), and the sequence (b n ) n∈N as in (2.2), it holds that

β g (LS) -β 0 1 = O P (b n ).
Following the lines of the proof of Theorem 3.1 we also obtain the proof of the following theorem.

Theorem 3.5. Under Assumptions (A1), (A2), and (A4), the condition in

(2.2), if |B 0 | < ∞ and λ n b -1 n -→ n→∞ 0, then β g (LS) -β 0 1 = O P (b n ).
The estimator of B 0 based on

β g (LS) = β 1,(LS) , • • • , β g,(LS) is given in a straightforward way as B n,(LS) = j ∈ {2, • • • , g}; β j,(LS) = β j-1,(LS) ,
and a similar result to the one in Theorem 3.2 can be derived again.

Theorem 3.6. Under the same assumptions as in Theorem 3.5, there exists a positive constant C 1 such that,

lim n→∞ P | B n,(LS) \ B 0 | ≤ C 1 max b n λ n , 1 b n = 1.
Similarly as for the quantile framework before, one can again improve the estimation accuracy of B 0 by taking the advantage of β g (LS) and defining the adaptive fused penalty with the corresponding empirical process

∨ U n (β g ) ≡ L n (β g ) + nλ n g j=2 ω n,j,(LS) β j -β j-1 q , (3.10)
where the weights ω n,j,(LS) are again constructed on the basis of fused group LS estimator as ω n,j,(LS) ≡ 1/ max n -1/2 , p k=1 | β j,k,(LS) -β j-1,k,(LS) | γ , for some fixed γ > 0 and β j,k,(LS) being the kth component of β j,(LS) . Thus, the adaptive fused group LS estimator is ∨ β g (aLS) ≡ arg min

β g ∈R rn ∨ U n (β g ),
and the corresponding estimator for B 0 is defined as

∨ B n,(aLS) ≡ j ∈ {2, • • • , g}; ∨ β j,(aLS) = ∨ β j-1,(aLS) .
As for the quantile framework, the sequence (λ n ) n∈N , in relations (3.9) and (3.10), can be different. Finally, using now the same arguments as in Theorem 3.3 and following the same lines of the proof, we obtain an analogous results also for ∨ β g (aLS) . Theorem 3.7. Under Assumptions (A1), (A2), and (A4), the condition in (2.2), if |B 0 | < ∞, the for any sequence

(λ n ) n∈N such that λ n b -1 n -→ n→∞ 0, it holds that ∨ β g (aLS) -β 0 1 = O P (b n ).
The results presented in Subsection 3.1 and Subsection 3.2 show that the estimated number of different successive groups is of the same order for both estimation frameworks with the adaptive fused approach and, moreover, the convergence rates of the corresponding estimators for the model parameters are also of the same order, all under the assumption that the true number of groups is bounded. The finite sample performance is investigated in the next section.

NUMERICAL STUDY AND APPLICATION

In this section we firstly present a Monte Carlo simulation study to show some numerical properties of the proposed fused methods for the varying number of groups, different sample sizes and error distributions. Later, the application on the air quality data is presented. The goal is to detect daily moments when the temperature and humidity contribution change their effect with respect to the maximum daily benzene concentration.

Numerical Study

The fused group quantile estimator β g defined in terms of the minimization in (3.3) and the adaptive fused group quantile estimator ∨ β g in (3.6) are both compared with the fused group LS estimator β g (LS) in (3.9) and its adaptive version ∨ β g (aLS) in (3.10) with respect to a wide range of different simulation settings. In order to make the comparison meaningful the quantile level of τ = 0.5 is considered. The dimension of the unknown group specific vector of parameters β ∈ R p is either p = 1 or p = 3 and three options are used for number of groups, g ∈ {20, 100, 200}. The sample is given as n = pg. The model covariates are randomly generated from the normal distribution and two distributions are used for the error terms (standard normal and Cauchy). The true number of different successive groups in the model is either 2, 5, 10, or n/5, where the last option (20 % of the sample size) clearly does not satisfy the model assumptions but it is still included in the simulation setup for the comparison purposes. Obviously, if there are two change points in the group parameter then there are three successive groups. Analogously for 5 changes in the group specific vector parameter-there are 6 successive groups. The corresponding locations of changes between successive groups are determined randomly and the jump magnitudes are also assigned randomly on the scale from 0.5 to 2 to allow for various signal-to-noise ratio. The regularization parameter equals λ n = n -1 (log(n)) 1/2 for the fused method and λ n = n -1 (log(n)) 5/2 for the adaptive fused method with the adaptive weights defined in (3.6) and (3.10) for γ = 1.

All four methods are compared with respect to the quality of the final fit and, mainly, the different successive coefficient group detection performance. The median (MED) of (Y i -Y i ) 1 i n and the L 1 norm of the difference between the true vector of parameters and its estimate are used to evaluate the estimation performance while the true recovery rate (the proportion of truly detected different successive coefficient groups with respect to all unknown changes) and the overestimation rate (proportion of the number of detected different successive coefficient groups with respect to the number of true changes) are used to assess the detection performance. The results are summarized in Table 1 (for p = 1) and Table 2 (for p = 3).

For M independent Monte Carlo replications let β (m) denote the estimate of β g by one of the four estimation methods for the m-th Monte Carlo run, with m = 1,

• • • , M . The corresponding forecast for Y i is Y i,(m) = X i β (m) . For each Monte Carlo replication the median error med (m) = median(Y i -Y i,(m) ; i = 1, • • • , n
) is obtained and the reported results are averaged over all M simulation runs MED= M -1 M m=1 med (m) . For the parameter estimation the value of MAD=

M -1 M m=1 1 pg pg j=1 |β 0 j -β j,(m) | is reported.
For some illustration of the model there is an example in Figure 1: the number of true different successive groups is two (out of g = 20 in total) and the number of the explanatory variables within each group is three (p = 3). The true vector of the group specific parameters for the first group (group indexes j ∈ {1, . . . , 15}) is β 0 1 = (1, 2, 3) and the true vector of the group specific parameters for the second group (group indexes j ∈ {16, . . . , 20}) is β 0 1 = (1.5, 1, 5) . The sample size is n = gp (n = 80). All four proposed estimation methods are applied and the and three explanatory variables in each group (p = 3). The first group specific vector parameter is the same for the groups j ∈ {1, . . . , 15} and it differs from the second group specific vector parameter, which is the same for the groups j ∈ {16 . . . , 20}. The Cauchy error terms are considered to visualize the robust favor of the quantile estimation approach for τ = 0.5 (right panel) when compared with the standard least squares (left panel).

corresponding estimates are given in Figure 1 1. Simulation results for the situation where the dimension of the unknown parameter is p = 1.

Two goodness-of-fit quantities are provided: the median (MED) of (Yi -Yi) 1 i n and the mean absolute difference (MAD) between the true parameter vector β 0 and the corresponding empirical estimate Two goodness-of-fit quantities are provided: the median (MED) of (Yi -Yi) 1 i n and the mean absolute difference (MAD) between the true parameter vector β 0 and the corresponding empirical estimate β.

The Recovery column is given in terms of two values: the proportion of truly discovered different successive coefficient groups (value 1 stands for all true changes being discovered) and the proportion between the number of estimated different successive coefficient groups and true changes (value 1 stands for the situation where the number of estimated changes equals the number of true changes). An ideal situation is 1.00/1.00

which means that all true changes are discovered with no other detections in addition. The results are averaged over 1000 Monte Carlo simulation runs.

From Tables 1 and2 for the Gaussian errors, we deduce that for |B 0 | ∈ {2, 5, 10}, the fused estimations for the least squares and the quantile methods have the same properties and the same also applies for the adaptive frameworks which, moreover, have the recovery detection rates close to one. On the other hand, if the assumption |B 0 | < ∞ does not hold, as for the models with 20% different successive groups, then not all of the different successive groups are detected and the performance is worse.

The robustness of the quantile methods is obvious when the Cauchy errors are used instead: while the LS based methods fail in both, the estimation and the group detection, the quantile approaches perform comparably well as in the situations with the Gaussian errors.

Application to Air Quality Data

In order to demonstrate the practical applicability of the proposed model we use the air quality data from De [START_REF] De Vito | CO, NO2 and NOx urban pollution monitoring with on-field calibrated electronic nose by automatic bayesian regularization[END_REF] which can be downloaded from the Machine Learning Repository site http://archive.ics.uci.edu/ml/datasets/Air+Quality#. The hourly meteorological and air quality data were recorded from March 2004 to February 2005. The idea is to use the daily temperature and humidity profiles (recorded every hours) to predict the maximum benzene concentration level for the given day. Optimally, it would be appropriate to use the temperature and humidity information only from some few instant moments during the day instead of recording both continuous profiles over the whole day. Given the data, there are g = 24 hourly groups and for each group there is the corresponding vector parameter β j = (β T j , β H j ) ∈ R 2 , for j = 1, . . . , 24, where β T j is responsible for the contribution of the temperature at 'j' o'clock and β H j models the effect of the humidity, again at 'j' o'clock. Using the model formulation from Section 2 and the estimation in terms of (3.3) it can be achieved that most of the corresponding parameter vector estimates are the same. If otherwise, then the existing changes in the vector estimates identify some specific daily segments with the same temperature and humidity contribution with respect to the maximum daily benzene concentration. The corresponding magnitudes for both effects in each daily segment are all estimated simultaneously.

Similarly as in the simulation section, four different models are fitted: the fused group LS approach and its adaptive version both presented in Figures 3(a) and 3(a) and the proposed fused group quantile and the adaptive fused group quantile in Figures 3(c) and 3(d). The temperature data and the humidity data are heavily skewed and, therefore, it can be assumed that robust approaches are more appropriate for this situation.

Indeed, while the fused group LS and also its adaptive version can not identify any specific daily moments which should be used to determine the maximum daily benzene concentration, the fused group quantile and its adaptive version in particular clearly identify some segments during the day when the contribution of the temperature and humidity is obvious. days out of 357 available days with full profiles in total. In addition, the maximum benzene concentration is recorded for each day and the corresponding time of the maximum occurrence (in hours) is given in terms of the frequency histograms in both panels.

PROOFS

Throughout the proofs, the following identity for the quantile check function ρ τ is be used: for any x, y ∈ R it holds that

ρ τ (x -y) -ρ τ (x) = y(1 1 x<0 -τ ) + τ 0 (1 1 x≤v -u x≤0 )dv.
(5.1)

Proof of Lemma 3.1. We will show that for all > 0, there exists a constant C > 0, such that for n large enough, we have

P inf u∈R rn , u 1 =1 G n β 0 + C b n u > G n (β 0 ) ≥ 1 -. (5.2)
Then, for any constant c 1 > 0, we can write the difference

G n β 0 + c 1 b n u -G n (β 0 ) using the form G n β 0 +c 1 b n u -G n (β 0 ) = E G n β 0 +c 1 b n u -G n (β 0 ) +W n u+ n i=1 (R i -E [R i ]), (5.3) with the r n -dimensional random vector W n ≡ c 1 b n n i=1 D i X i , the random variables D i ≡ (1 - τ )1 1 {ε i <0} -τ 1 1 {ε i ≥0} , and R i ≡ ρ τ (ε i -c 1 b n X i u) -ρ τ (ε i ) -c 1 b n D i X i u.
Using the Holder's inequality, we have that |X i u| ≤ X i ∞ u 1 . Then, for all u ∈ R rn such that u 1 = 1, by Assumption (A1), we have that |X i u| ≤ C.

Firstly, we study the first term on the right-hand side of relation (5.3). Using the identity in , β H j ) ∈ R 2 , for j = 1, . . . , 24, for four different estimation techniques: fused group LS, adaptive fused group LS, fused group quantile and adaptive fused group quantile. The adaptive fused group quantile estimator in panel (d) clearly identifies some instant moments during a day when the temperature and humidity information is relevant for the maximum benzene concentration. In other words, it seems enough to record the temperature and humidity information at 2 pm and, also, after 6 pm.

(5.1), we obtain

G n β 0 + c 1 b n u -G n (β 0 ) = -c 1 b n n i=1 X i uD i + n i=1 c 1 bnX i u 0 [1 1 {ε i <v} -1 1 {ε i <0} ]dv.
Applying now the mean value theorem, taking into account the fact that the derivative of f is bounded in a neighborhood of zero by Assumption (A3), and the fact that u 1 = 1, E [D i ] = 0 and b n → 0, we obtain

E G n β 0 + c 1 b n u -G n (β 0 ) = n i=1 E c 1 bnX i u 0 1 1 {0<ε i <v} dv = n i=1 c 1 bnX i u 0 [F (v) -F (0)]dv = f (0) 2 c 2 1 b 2 n n i=1 (X i u) 2 + o b 2 n n i=1 u (X i X i )u .
Using Assumption (A2) together with f (0) > 0, we get that

n -1 E G n β 0 + c 1 b n u -G n (β 0 ) = Cf (0)b 2 n 1 n n i=1
u X i X i u(1 + o(1)) > 0.

(5.4)

Next, we study the last two terms on the right-hand side of relation (5.3). For the last term we have, with probability one, for any i = 1, . . . , n, that |R i | ≤ c 1 b n |X i u|1 1 {|ε i |≤c 1 bn|X i u|} . Since (ε i ) 1 i n are independent, then the random variables (R i ) 1 i n are independent as well and, therefore

E n i=1 (R i -E [R i ]) 2 = n i=1 E [R i -E [R i ]] 2
(5.5)

≤ n i=1 E [R 2 i ] ≤ Cb 2 n n i=1 |X i u| 2 E 1 1 {|ε i |≤c 1 bn|X i u|} .
Using the fact that the density f is bounded in a neighborhood of 0 by assumption (A3), adopting the Taylor's expansion, Cauchy-Schwarz and Holder inequalities, Assumption (A1), and the fact that u 1 = 1, we obtain β j -β j-1 q ≥ nλ n p -1+1/q j∈B∩B 0 β j -β j-1 1

E 1 1 {|ε i |≤c 1 bn|X i u|} = 2c 1 b n |X i u|f (d i,n ) ≤ Cb n max 1 i n X i ∞ = Cb n , (
= O P nλ n (|B ∩ B 0 |)b n > 0.
We have also S 4n = Cnλ n ≥ 0, therefore, taking into account the fact that the difference Q n (β g ) -Q n (β 0 ) must be negative for the minimizer β g in (5.12), using the relations in (5.13)) and (5.16), we deduce that nb 2 n + nλ n ≥ nλ n (|B ∩ B 0 |)b n , which also implies that

|B \ B 0 | ≤ C max b n λ -1
n , b -1 n . This finishes the proof.

Proof of Theorem 3.3. In this case, for a positive constant c 2 > 0, a vector u ∈ R rn such that u 1 = 1, we study the difference

∨ Q n (β 0 + c 2 b n u) - ∨ Q n (β 0
). The penalty related to this difference, similarly as in (5.11), becomes nλ n g j=2 ω n,j β 0 j + c 2 b n u j -(β 0 j-1 + c 2 b n u j-1 ) q -β 0 j -β 0 j-1 q ≥ -c 2 nλ n b n j∈B 0 ω n,j u j -u j-1 q .

Taking into account the relation in (3.7) and using similar arguments as in the proof of Theorem 3.1, we obtain that

∨ Q n (β 0 + c 2 b n u) > ∨ Q n (β 0
), which holds with probability converging to 1, as n → ∞.

Proof of Theorem 3.4. The proof is very similar to that of Theorem 3.2. We only give the main results, using the same notation as in the proof of Theorem 3.2. For β g ∈ V n (β 0 ), the difference between the adaptive processes can be expressed as

∨ Q n (β g ) - ∨ Q n (β 0 ) ≡ ∨ S 1n + ∨ S 2n + ∨ S 3n - ∨ S 4n ,

Figure 1 .

 1 Figure 1. An illustration of the model in (2.1) for two truly different successive groups (out of g = 20 in total)

  Model with 2 different groups Model with 5 different groups Model with 10 different groups Model with 20 %

  in terms of two values: the proportion of truly discovered different successive coefficients (value 1 stands for all true changes being discovered) and the proportion between the number of estimated different successive coefficients and true changes (value 1 stands for the situation where the number of estimated changes equals the number of true changes). An ideal situation is 1.00/1.00 which means that all true changes are discovered with no other detections in addition. The results are averaged over 1000 Monte Carlo simulation runs. D p = 3 Model with 2 different groups Model with 5 different groups Model with 10 different groups Model with 20 %

Figure 2 .

 2 Figure 2. Daily temperature profiles (left panel) and daily humidity profiles (right panel) for 50 randomly selected

Figure 3 .

 3 Figure3. The estimated parameter vectors β j = ( β T j , β H j ) ∈ R 2 , for j = 1, . . . , 24, for four different estimation techniques: fused group LS, adaptive fused group LS, fused group quantile and adaptive fused group quantile. The adaptive fused group quantile estimator in panel (d) clearly identifies some instant moments during a day when the temperature and humidity information is relevant for the maximum benzene concentration. In other words, it seems enough to record the temperature and humidity information at 2 pm and, also, after 6 pm.

  5.6)with d i,n between c 1 b n |X i u| and (-c 1 b n |X i u|). Then, using Assumption (A1) together with the relations in (5.5) and (5.6), and the fact that |X i u| ≤ C, we have also T 1n = o P n β gβ 0 2 1 , we have for the relation in (5.14) thatS 1n = O P (nb 2 n ).(5.16) For (5.13) it remains to study the sums S 2n , S 3n , and S 4n . Sinceβ g ∈ V n (β 0 ),together with the fact that the cardinality |B 0 | is bounded and λ n b -1 n -→ n→∞ 0, we obtain S 2n = O P (nλ n b n ) = o P (nb 2 n ) and also S 3n ≡ nλ n j∈B∩B 0

Table 2 .

 2 Simulation results for the situation where the dimension of the unknown parameter is p = 3.

	0.76/3.37	1.00/3.03	1.00/2.50	0.41/1.65	0.97/1.12	0.99/1.02	0.77/3.48	1.00/4.28	0.82/1.76	0.22/0.89	0.70/0.72	0.51/0.69	0.65/3.67	0.81/3.55	0.79/3.24	0.42/2.20	0.58/1.99	0.62/2.01	0.62/3.40	0.97/4.89	0.77/1.87	0.23/1.00	0.37/0.71	0.43/0.66
	-0.01 0.34	0.00 0.14	0.00 0.10	-0.05 0.37	0.04 0.08	0.02 0.05	-0.01 0.35	0.05 0.17	-0.06 0.36	-0.07 0.61	-0.07 0.38	-0.15 0.51	0.00 3.05	-0.01 3.16	-0.16 2.42	1.71 3.37	-1.53 3.42	-1.16 2.68	0.00 0.73	-0.63 0.67	-1.12 0.57	0.02 0.75	-0.24 0.60	-0.32 0.60
	0.65/1.11	1.00/6.72	1.00/11.03	0.39/0.64	1.00/1.04	1.00/1.17	0.57/1.02	1.00/10.08	1.00/4.73	0.29/0.38	0.88/1.31	1.00/1.00	0.65/1.21	0.82/7.64	0.80/14.52	0.45/0.78	0.62/4.23	0.65/8.49	0.56/1.09	0.98/10.38	1.00/6.18	0.22/0.35	0.49/1.61	0.95/1.48
	0.02 1.00	0.00 0.08	-0.07 0.07	-0.02 1.09	0.00 0.03	0.00 0.02	0.08 1.01	-4.86 0.91	-0.02 0.03	-0.16 1.15	-0.21 0.15	0.00 0.02	0.01 3.54	-0.22 3.12	-0.12 2.30	1.74 3.91	-1.58 3.29	-1.17 2.54	0.05 1.17	-7.85 1.62	-1.27 0.22	-0.17 1.28	-0.58 0.39	-0.41 0.11
	0.92/2.63	1.00/15.61	1.00/20.40	0.63/1.31	1.00/1.13	1.00/1.02	0.92/2.70	1.00/24.75	1.00/8.44	0.33/0.63	0.80/3.06	1.00/1.00	0.70/2.73	0.83/17.35	0.82/31.83	0.45/1.64	0.62/9.17	0.67/18.97	0.65/2.53	1.00/24.61	1.00/27.59	0.22/0.71	0.67/5.14	0.89/2.41
	0.02 0.40	0.00 0.06	0.00 0.03	0.07 0.44	0.00 0.02	0.00 0.01	0.00 0.41	-8.29 2.47	0.00 0.02	0.00 0.53	-0.09 0.13	-0.01 0.01	0.01 3.26	0.05 3.18	-0.10 2.29	1.76 3.45	-1.65 3.32	-1.11 2.49	0.02 0.71	-9.66 2.30	-2.45 0.40	0.12 0.71	-0.56 0.26	-0.26 0.07
	1.00/9.87	1.00/62.67	1.00/87.79	0.99/1.82	1.00/2.00	1.00/1.26	1.00/9.92	1.00/99.00	1.00/29.92	0.85/0.90	0.07/2.57	1.00/1.00	0.73/10.85	0.86/71.60	0.79/128	0.57/6.26	0.68/35.92	0.67/75.66	0.76/9.84	1.00/97.88	1.00/109	0.32/1.34	0.48/9.37	0.99/8.28
	17	0.00 0.05	0.00 0.02	-0.02 0.08	0.00 0.01	0.00 0.01	-0.04 0.18	-5.79 2.05	-0.01 0.02	0.03 0.24	0.06 0.07	0.00 0.01	0.00 3.24	-0.08 3.36	-0.03 2.29	1.69 3.54	-2.04 3.60	-1.22 2.49	-0.02 0.55	-8.20 2.03	-2.66 0.39	-0.02 0.48	-0.11 0.12	-0.52 0.06
		100	200	20	100	200	20	100	200	20	100	200	20	100	200	20	100	200	20	100	200	20	100	200
		β g (LS)			∨ β g (aLS)			β g			∨ β g			β g (LS)			∨ β g (aLS)			β g			∨ β g
													C

ACKNOWLEDGEMENTS

The authors would like to sincerely thank the Editor-In-Chief, Professor Nitis Mukhopadhyay, an Associate Editor and the anonymous referee for the time taken to review this paper. The work of Matus Maciak was partially supported by the Czech Science Foundation project GA ČR (No. 18-0522Y).

We consider a deterministic sequence (a n ) n∈N such that: a n → ∞ and nb 3 n a n n 2 b 4 n . An example of such sequence is a n = (log n) 3/2 if b n = n -1 log n 1/2 .

Considering the relation in (5.7) and since a n nb 3 n , then also

which implies, by the Bienaymé-Tchebychev inequality, that the last term of the right-hand side of the relation in (5.3) equals to

(5.8)

Finally, we study the second term of the right-hand side in (5.3). By the Central Limit Theorem (CLT) for the independent random variables ), and by the condition in (2.2) where n 1/2 b n → ∞, since a n n 2 b 4 n , together with the relations in (5.4) and (5.8), we have for (5.3) the following:

Therefore, the relation in (5.2) is proved. Moreover, it implies that β gβ 0 1 = O P (b n ) and, therefore, the lemma is proved.

Proof of Theorem 3.1. In order to prove the assertion of the theorem, let us consider a vector u ∈ R rn , such that u 1 = 1 and a constant c 2 > 0. Then the following holds

(5.9)

On the other hand, since u 2 ≤ u 1 = 1, by the proof of Lemma 3.1, we have with the probability converging to 1, that

(5.10)

If the components of u are denoted as u 1 , • • • , u g , then, using the triangular inequality, for the penalty in (5.9), we have (5.11) where for the last equality in (5.11) we have used the fact that

) and taking into account the relation in (5.10), we obtain for (5.9) and (5.11) that

which holds with the probability converging to 1, as n → ∞.

Proof of Theorem 3.2. By Theorem 3.1 we have

(5.12)

with the neighborhood V n (β 0 ) of β 0 with the radius c 2 b n defined as

for some constant c 2 > 0. Then, in order to prove the assertion of the theorem we consider the parameter vector β g = (β 1 , . . . , β g ) ∈ V n (β 0 ) and the index set B ≡ j ∈ {2, • • • , g}; β j = β j-1 . Note, that B and β g both depend on n and the vector of true unknown parameters β g is not random. Therefore, we consider only B ∩ B 0 = ∅, otherwise the theorem trivially holds.

Let us concentrate on the following decomposition:

(5.13)

Using the identity in (5.1) we can write the sum S 1n as

(5.14)

For T 1n , we have E [T 1n ] = n i=1 (β gβ 0 ) X i F (0) -F (0) = 0 and using Assumptions (A1), (A2), and (A3), we obtain for the variance that

Then, by the Law of Large Numbers, we also have

dv, we can apply the Taylor expansion

for some ṽ between 0 and v. Since the derivative f is bounded in some neighborhood of zero, taking into account Assumption (A1), we obtain

(5.15)

On the other hand, since the error terms (ε i ) 1 i n are independent, we have

where for brevity, we used the notation where

Hence, taking into account this last relation together with (5.15), since β g ∈ V n (β 0 ), b n → 0 as n → ∞, and applying the Bienaymé-Tchebychev inequality, we obtain

where S 1n is defined in (5.13) and the other sums are

For ∨ S 2n , taking also into account the relation in (3.7), similarly as for S 2n in (5.13), we ob- Therefore, for the vector parameter β g ∈ B which minimizes

, which holds with the probability converging to one as n → ∞. This also implies

Proof of Lemma 3.2.

For any constant c 1 > 0 and some r n -vector u, such that u 1 = 1, we have

(5.17)

By Assumption (A1), we have |X i u| ≤ C. Therefore, using Assumption (A4) and CLT we get n i=1 X i ε i u = O P (n 1/2 ). By Assumption (A2), we also get n i=1 (X i u) 2 = O(n) and taking into account the condition in (2.2), we get that (5.17) is

Thus, for any > 0, there exists a positive constant C > 0, such that,

Proof of Theorem 3.6.

The proof is similar to that of Theorem 3.2 with the only difference that for (5.13) the sum S 1n equals n i=1

which is, using the same arguments as in the proof of Lemma 3.2, of the order O P (nb 2 n ). The rest of the proof is omitted because it follows the same lines as the proof of Theorem 3.2.

Proof of Theorem 3.7. The proof follows the same lines as the proof of Theorem 3.3 and, therefore, it is omitted.


Jean Auriol

Islam Boussaada

Hugues Mounier

Silviu-Iulian Niculescu

Torsional-vibrations Damping in Drilling Systems: Multiplicity-Induced-Dominancy based design

The interest in using delayed design is emphasized in recent research on the stabilization of finite/infinite dimensional dynamical systems. In particular a property called multiplicity-induced-dominancy (or MID for short) allowing for a reduced-complexity delayedcontroller design showed its efficiency in fast damping of harmful oscillations. This contribution is concerned with oil-well drilling systems torsional vibrations, which constitute an important source of economic losses; drill bit wear, pipes disconnection, borehole disruption and prolonged drilling time, among other consequences. Such torsional vibrations are assumed to be governed by a wave equation with weak damping term. The MID-based design is further exploited to quench the torsional vibrations along the rotary drilling system. The proposed control law guarantees the existence of robustness margins with respect to delays and parameters uncertainties. Numerical investigation of the performance as well as the robustness of the corresponding control law will be presented.

INTRODUCTION

In this paper we consider the stabilization of torsional oscillations for a drilling system. Our approach is based on a multiplicity-induced-dominancy (MID) property allowing the design of a reduced-complexity delayed-controller (adjusted from Boussaada et al. (2019a)). The resulting closed-loop system is proved to be robust to uncertainties on the parameters and to delays in the loop.

The dynamical behavior of drill-strings is complex as many dynamic phenomena are involved, as vibrations, bending and twisting quasi-static motion, bit-rock interactions [START_REF] Kapitaniak | Unveiling complexity of drill-string vibrations: Experiments and modelling[END_REF]; [START_REF] Spanos | Oil and gas well drilling: a vibrations perspective[END_REF]. In particular, the drill-string interaction with the borehole gives rise to a wide variety of non-desired oscillations [START_REF] Dunayevsky | Application of stability approach to bit dynamics[END_REF]; [START_REF] Jansen | Nonlinear dynamics of oil well drillstrings[END_REF]; [START_REF] Saldivar | Stick-slip oscillations in oillwell drilstrings: distributed parameter and neutral type retarded model approaches[END_REF] which can be classified depending on the direction they appear. Among them, the longitudinal (or axial) oscillations may cause a bit-bouncing effect that is characterized by a repetitive loss of contact of the bit with the bottom of the well [START_REF] Aarsnes | Linear stability analysis of self-excited vibrations in drilling using an infinite dimensional model[END_REF], [START_REF] Depouhon | Instability regimes and self-excited vibrations in deep drilling systems[END_REF], [START_REF] Germay | Multiple mode analysis of the self-excited vibrations of rotary drilling systems[END_REF]; [START_REF] Zhou | Adaptive predictor control for stabilizing pressure in a managed pressure drilling system under time-delay[END_REF]. Torsional vibrations can appear due to downhole conditions (such as significant drag, tight or formation characteristics [START_REF] Nandakumar | Stability analysis of a state dependent delayed, coupled two dof model of drill-stringvibration[END_REF] for instance) or due to side forces induced by Coulomb friction Aarsnes and Shor (2018). They are known as stick-slip and are characterized by a series of stick (the bit stop rotating) and slip (a sudden release of energy) [START_REF] Kamel | Modeling and analysis of stick-slip and bit bounce in oil well drillstrings equipped with drag bits[END_REF]. These self-excited unde-sirable vibrations may cause fatigue of the equipment or a deterioration of the performance of the process implying a reduction of the Rate of Penetration (ROP). They can also lead to premature failing of the bit or provoke premature wear and tear of drilling equipment resulting in fatigue and induced failures such as pipe wash-out and twist-off [START_REF] Mason | Addressing bha whirl-the culprit in mobile bay[END_REF]. Thus, they may cause catastrophic damages and at least wear to expensive components of the drill-string [START_REF] Kriesels | Cost savings through an integrated approach to drillstring vibration control[END_REF]. In this context, a clear understanding of drill-string dynamics appears to be crucial to control these vibrations and consequently improve the performance of drilling systems (ROP), prevent any eventual damage and reduce safety risks. This explains why an extensive research effort has been conducted in the last five decades to suppress these undesirable oscillations and in so far reduce the costs of failures and increase the ROP (see the references cited in the survey [START_REF] Leine | Literature survey on torsional drillstring vibrations[END_REF] and in [START_REF] Saldivar | A control oriented guided tour in oilwell drilling vibration modeling[END_REF]).

During the drilling process, the operator wants to control the downhole behavior of the drill-string (e.g. reach a given rpm, a given orientation...) and optimize the ROP, while avoiding undesired oscillations. To do so, it is usually possible to impose (using the rotary table) the weight on the drill-string and the torque at the surface. Automated control laws have been designed to solve such control problems [START_REF] Serrarens | H-infinity control for suppressing stick-slip in oil well drillstrings[END_REF]. Recent control law (which are not simple PID controllers) are usually based on distributed parameter models. More precisely, the drill-string dynamics can be modeled by a set of hyperbolic partial differential equations (namely wave equations) coupled with ODEs through their boundaries Di Meglio and [START_REF] Di Meglio | A distributed parameter systems view of control problems in drilling[END_REF]. Such wave equations can then be rewritten as neutral systems. In general, it may not be possible to efficiently stabilize such systems using classical control laws. Specific techniques have to be developed as an approximation of delay systems by a finite-dimensional representation (early-lumping) may lead to the same degree of complexity in the control design which can be potentially harmful for the stability of the system. Consequently, different approaches have been developed to provide realistic and effective control laws. A flatness-based approach has for instance been developed in Aarsnes et al. (2018). More recently, a low complexity control law with a reduced number of tuning gains has been proposed in Boussaada et al. (2019a) in the case of a vertical well. The controllers gains are chosen such that the closed-loop system has a particular structure and the corresponding spectrum lies in the complex left-half plane.

In this paper, we modify the control law introduced in Boussaada et al. (2019a), which is based on the MID property, by adding a derivative action, thus introducing an additional degree of freedom that facilitates a targeted placement of the poles of the closed-loop system. The proposed control law only requires the measurement of the states at the boundaries of the system and consequently does not need the design of a state observer. This is of specific interest as the design of an observer for this class of system may be difficult or computationally expensive (see [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF]; [START_REF] Auriol | Combining formation seismic velocities while drilling and a pde-ode obvserver to improve the drill-string dynamics estimation[END_REF] in which stateobservers are designed usign the backstepping approach.) Moreover, we analyze the robustness properties of the resulting closed-loop system. More precisely, we show that it is robust with respect to uncertainties on the parameters and delays in the actuation. These results are stated in a general framework of quasi-polynomials functions since, for the proposed control law, the closed-loop system can be expressed using such a structure. It is worth mentioning that the rightmost root for quasipolynomial function corresponding to stable retarded time-delay systems (also in the neutral case under some assumptions) is actually the exponential decay rate of its time-domain solution, (see [START_REF] Mori | On an estimate of the decay rate for stable linear delay systems[END_REF] for an estimate of the decay rate for stable linear delay systems or [START_REF] Boussaada | Multiplicity and stable varieties of time-delay systems: A missing link[END_REF] for the dominancy properties). This dominancy property presented in [START_REF] Hayes | Roots of the transcendental equation associated with a certain difference-differential equation[END_REF] was used by the authors to design the control law in Boussaada et al. (2019a). The robustness aspects are considered by means of an algebraic analysis in the Laplace domain.

The paper is organized as follows. In Section II, we give the PDE model describing the torsional vibration of an oilwell rotary drilling system. This model is then expressed as a neutral equation. A stabilizing control law that guarantees the damping of the different undesired vibrations is given in Section III. Section IV is devoted to the robustness analysis of the proposed control law. The different results are expressed in a general framework of quasi-polynomials functions. An illustrative example completes the paper in Section V. Some concluding remarks are given in Section VI.

PROBLEM FORMULATION

In this section we start by providing an hyperbolic partial differential equations model to describe the torsional and axial motions of the drill-string vibrations. Then, we reduce such a model to a functional differential equation of neutral type.

A Coupled Axial/Torsional Vibrations Model

Let us consider a vertical well for which the rotary table (located at the top of the drill-string) sets the drillstring into a rotary motion around its main axis. The dynamics of interest can be obtained by assuming elastic deformations and using equations of continuity and state. More precisely, the torsional and axial vibrations can be modeled by a set of hyperbolic PDEs (namely wave equations) coupled through their boundaries. In what follows, we denote L the length of the drill-string, Φ(z, t) the rotary angle and U (z, t) the longitudinal position. They are are functions of (z, t) evolving in {(z, t)| z ∈ [0, L], 0 < t < T }, the position z = 0 corresponding to the top of the drill-string and the position z = L to the bottom. We have [START_REF] Boussaada | Analysis of drilling vibrations: A timedelay system approach[END_REF]:

∂ 2 z Φ(z, t) = c 2 ∂ 2 t Φ(z, t), (1a) ∂ 2 z U (z, t) = c 2 ∂ 2 t U (z, t), (1b)
with the boundary conditions

GJ∂ z Φ(0, t) = β∂ t Φ(0, t) -u T (t) (2a) GJ∂ z Φ(L, t) = -I B ∂ 2 t Φ(L, t) -T (∂ t Φ(L, t)) , (2b) and EΓ∂ z U (0, t) = α∂ t U (0, t) -u H (t) (3a) EΓ∂ z U (L, t) = -M B ∂ 2 t U (L, t) -T (∂ t Φ(L, t)) . (3b)
The axial and torsional propagation velocities v U , v Φ , are defined as: v U = c -1 and v Φ = c-1 . They depend on the system physical parameters (Young modulus E, shear modulus G and density ρ a ), by means of :c = ρa E and c = ρa G . The two actuators are denoted u H and u T . They respectively correspond to the brake motor control (upward hook force) and to the torque produced by the rotary table motor. The functions α∂U t (0, t) corresponds to a friction force of viscous type (where α is the viscous friction coefficient). Finally, the term β∂ t Φ(0, t)u T (t) corresponds to the difference between the motor speed and the rotational speed of the first pipe.

The different geometrical and physical parameters of the drill-string are assumed to be spatially constant. Note that this assumption may be idealistic since the Bottom Hole Assembly (lower part of the drillstring) is usually made of bigger pipes with different physical properties. However, the length of the BHA is usually much lower and its effect can consequently be lumped into the ODE. We denote Γ the shear modulus, J the polar moment of inertia, M B the mass of the drill-string, I B the inertia moment 1 of the drill bit. Finally the function T corresponds to the frictional torque resulting from the interaction between the drill bit and the rock. We model it by the following odd function Boussaada et al. (2012a): T (x) =2 kx k2 +x 2 , for which we have the following local approximation:

T (x) = 2 k x + O x 3 ,
For a further discussion on friction models in rotary drilling systems, the reader is invited to see [START_REF] Marquez | Analysis and control of oilwell drilling vibrations[END_REF]. This frictional torque is one of the cause of the torsional oscillations that creates the stick-slip phenomenon. However, it is worth mentioning that the stick-slip oscillations may also be caused by a negative difference between static and kinematic along-string Coulomb-type friction Aarsnes and Shor (2018). It emphasizes the action of nonlinear forces along the drill-string (which are combined with the bit rock interaction) in the torsional oscillatory behavior of drilling systems. Although the effect of these side-forces is not really important for vertical wells, they have a significant importance when the well is not vertical. The numerical values of the different parameters are given in Table 1.

Table 1. Numerical values of the system parameters.

Sym

Neutral-type model for the torsional oscillations

In this section we only consider the evolution of the rotary angle Φ(z, t) that corresponds to the torsional oscillations (1a) with the B.C. (2a)-(2b). Using the method of the characteristics, it is possible to rewrite the state Φ(L, t) as the solution of a neutral-type equation, thus relating the system variables at both ends of the drilling rod. We can show that there is a one-to-one correspondence between the solutions of the mixed problem for hyperbolic PDE and the initial value problem for the associated system of functional equations [START_REF] Rasvan | A method for distributed parameter control systems and electrical networks analysis[END_REF]. This is of interest as it implies that the techniques developed for time delay systems can be used for hyperbolic PDEs.

Namely, as described in Boussaada et al. (2019a), using the linearization of T given above, the system (1a) with the boundary conditions (2a)-(2b) reduces to a to an I/O system which is a neutral system of order 3 :

                                     (- βI B α L 2GJc - I B 2α L ) ... ΦL(t)+( βI B α L 2GJc - I B α L 2 ) ... ΦL(t -2τ) -( β 2α L + GJc 2α L ) ΦL (t)+( GJc 2 α L - βα L 2 ) ΦL (t -2τ) -(2 GJc 2α L ζ + β 2α L ) ΦL (t)+(2 GJcζα L 2 - βζα L 2 ) ΦL (t -2τ) + ζ 2 GJc 2α L Φ L (t) + ζ 2 GJc 2α L Φ L (t -2τ ) = 1 k(α L ΦL (t -2τ)+α -1 L ΦL (t)) - 2 k(α L ΦL (t -2τ) -α -1 L ΦL (t)) -uT (t -τ ) -ζu T (t -τ ), (4)
The delay τ = cL represents the total transport time for the torsional wave to travel from one to the other extremity of the drillstring.

LOW COMPLEXITY CONTROLLER DESIGN

In this section, we design a stabilizing control law that stabilizes (4). This control law is inspired from Boussaada et al. (2019a) and is based on the multiplicity induceddominancy (MID) property which was underlined in [START_REF] Boussaada | Multiplicity and stable varieties of time-delay systems: A missing link[END_REF][START_REF] Boussaada | Further remarks on the effect of multiple spectral values on the dynamics of timedelay systems. application to the control of a mechanical system[END_REF] for retarded time-delay systems and exploited in several applications [START_REF] Boussaada | On the coalescence of spectral values and its effect on the stability of time-delay systems: Application to active vibration control[END_REF]Boussaada et al. ( , 2018a)); Boussaada and Niculescu (2018). However, due to the specific structure of the system (neutral), the MID properties for time-delay systems of neutral type had to be adjusted. This has been done in Boussaada et al. (2019a) for instance. Compared to Boussaada et al. (2019a), the control law we propose in this paper includes a action on the derivative term. This introduces a supplementary degree of freedom that can be used to enforce the characteristic equation of the closed-loop system to have multiple negative spectral values. Under appropriate conditions, such a multiple root may be the dominant one among the spectrum, guaranteeing the stability of the closed-loop system. Adding this additional degrees of freedom (by means of this derivative term) means that the multiplicity of the dominant root can be increased compared to Boussaada et al. (2019a), thus improving the closed-loop performance. More precisely the control law we consider is defined by u

T (t) = κ 1 Φ L (t -τ ) + κ 2 ΦL (t -τ ),
where the delay τ is the same as the intrinsic delay of the nominal system. Consequently, we have uT (t) = κ 1 ΦL (tτ ) + κ 2 ΦL (t -τ ) and the closed-loop system (4) rewrites

         a 3 ... ΦL(t) + b 3 ... ΦL(t -2τ) a 2 ΦL (t) + b 2 ΦL (t -2τ ) a 1 ΦL (t) + b 1 ΦL (t -2τ) + a 0 Φ L (t) + b 0 Φ L (t -2τ ) = 0, (5)
where:

                                 a 3 = - βI B α L 2GJc - I B 2α L , b 3 = βI B α L 2GJc - I B α L 2 , a 2 = - β 2α L - GJc 2α L - 3 kα L , b 2 = κ 2 + GJc 2 α L - βα L 2 + α L k , a 1 = - GJc α L ζ - β 2α L , b 1 = κ 1 + κ 2 ζ + GJcζα L - βζα L 2 , a 0 = ζ 2 GJc 2α L , b 0 = κ 1 ζ + ζ 2 GJc 2α L where c = ρ a G , ζ = γ 2c 2 , α L = e -cLζ , τ = 2τ = 2cL, (6)
Such a system is of neutral type (since a 3 = 0 and b 3 = 0). Using the Laplace transform, we can rewrite it as a quasipolynomial including one delay.

     L(s) = P (s) + Q(s)e -sτ , P (s) = 3 i=0 a k s k , Q(s) = 3 i=0 b k s k , (7)
The main results can be resumed as follows.

Proposition 1. Boussaada et al. (2019a) The linearized system (5) subject to the controller u T (t) = κ 1 Φ L (t-τ )+ κ 2 ΦL (t -τ ) has the following properties in closed-loop :

• The maximal multiplicity of any root of the quasipolynomial L given in (7) with arbitrary coefficients (a k , b k ) 0≤k≤3 is bounded by 7.

• If | b k a3 | < 1 ∀ 0 ≤ k ≤ 3
and such a maximal multiplicity is reached, then the corresponding spectral value is real and dominant.

The proposed control law u T introduces two degrees of freedom by means of κ 1 and κ 2 . These two tuning parameters can be used to modify the dominant root of L. This dominant root can somehow be seen as an additional tuning parameter (under the constraint that the coefficients (b 0 , b 1 , b 2 ) must verify the constraint ( 5)). Having three degrees of freedom, it may become possible to find a dominant root of order 3. This is a considerable improvement compared to Boussaada et al. (2019a), where the dominant root was simply of order 2. More precisely, we have the following proposition. Proposition 2. If the system parameters are chosen in Table 1 and if b 0 , b 1 , b 2 verify ( 8)-( 10), then s 0 = -3 is a triple root for (7). In particular, if the parameter values are the ones given in Table 1 then the controller gains satisfy κ 1 = -16929630.76, κ 2 = -302236.4149 and the triple root at s 0 = -3 is dominant. In such a case the trivial solution is asymptotically stable.

Sketch of the proof

The proof is based on the argument principle applied on the translation of standard Bromwich contour at s 0 = -3, which allows to count the roots of the quasipolynomial with (s) > -3. by defining a specific integration contour γ. Then an integral over γ is defined as the sum of the integrals over the directed smooth curves that make γ up. See [START_REF] Boussaada | Multiplicity-induceddominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF] for details.

It is worth mentioning that the proposed control law only requires the value of the state at the boundary x = L. If a sensor located at x = L is available, this implies that we do not need to design a state-observer. This is of specific interest as it means that the proposed approach is not computationally expensive. More precisely, the design of a state observer may require solving a set of PDEs in real time (see [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF]; [START_REF] Auriol | Combining formation seismic velocities while drilling and a pde-ode obvserver to improve the drill-string dynamics estimation[END_REF] for the design of state-observers using the backstepping approach). However, from a physical point of view it may be hard and expensive to have a sensor available at the bottom of the well. Moreover, the frequency of such (noisy) measurements would be extremely low. Thus, for the considered drilling example, the design of an observer may be necessary. We believe that this can be done following a similar approach, which would be much easier to implement compared to [START_REF] Aarsnes | Estimating friction factors while drilling[END_REF] (even if in this reference, the authors consider a more general case of a non-vertical well subject to non-linear side-forces friction terms). However, for other physical applications for which such a sensor would be available, this is an crucial asset of our approach.

ROBUSTNESS ASPECTS

In this section, we consider the robustness properties of the control law u T (t) = κ 1 Φ L (t -τ ) + κ 2 ΦL (t -τ ) we have designed in the previous section. These results will be stated in a general framework. More precisely, let us consider two quasi-polynomial functions G(s) and K(s) respectively defined by

G(s) = n k=0 a k s k + n k=0 b k s k e -sτ , (11)
K(s) = p k=0 c k s k + p k=0 d k s k e -sτ , (12)
where s denotes the Laplace variable, and where the coefficients a k , b k , c k , d k and τ are real and constant. We assume that a n b n = 0. This means that the equation G(s) = 0 corresponds to the characteristic equation associated to a neutral system of order n. Let us consider the system defined in the Laplace domain by G(s)X(s) = u(s), (13) where X denotes the state of the system and where u denotes the actuation. Let us define the control law u as a feedback with the following structure u(s) = K(s)X(s).

(14) Consequently, the closed-loop system rewrites as a neutral equation (G(s) -K(s))X(s) = 0.

(15) In this section, we will show that if the control operator K(s) guarantees the exponential stability of the closedloop system (15), then, under some assumptions on the structure of K, this closed-loop system is robust to small uncertainties in the parameters and to delays in the actuation. The robustness with respect to delays in the actuation (delay-robustness) requires a nice behavior of the open-loop system at high frequencies. More precisely, we need to guarantee that the principal part of the operator associated to G(s) generates an exponentially stable semigroup. This leads to the following assumption Assumption 3. The coefficients b n and a n verify

|b n | < |a n |.
If this assumption does not hold, then the operator G(s) has an infinite number of zeros in the complex right half plane (RHP) which means that any linear control will lead to a zero delay margin [START_REF] Logemann | Conditions for robustness and nonrobustness of the stability of feedback systems with respect to small delays in the feedback loop[END_REF], Theorem 1.2) i.e. for any control law u(t) the closed-loop system becomes unstable when there is an (arbitrarily small) delay δ > 0 in the loop. Finally, the approach we propose to prove the existence of robustness margins requires the following assumption Assumption 4. The degree of the the quasi-polynomial function K is chosen such that n > p (i.e. the order of the open-loop transfer function is larger than the order of the control transfer function).

Let us now assume that the different parameters of the system (i.e. the coefficients that appear in the definition of the function G) are subject to small uncertainties and that there is a delay in the loop. More precisely, let us assume that the real plant operator now rewrites We also consider a delay δ > 0 acting on the actuation. Consequently, the closed-loop system now rewrites (G u (s) -K(s)e -sδ )X(s) = 0.

(17) The next proposition guarantees the existence of robustness margins for the control law u(t) defined in ( 14). Proposition 5. Suppose that Assumption 3 and Assumption 4 are satisfied. Suppose that the control law u(s) = K(s)X(s) stabilizes the nominal system (13). Then, there exist δ max > 0 and κ 0 > 0 such that for every delay δ < δ max and for every set of uncertainties δ a k , δ b k , δ c k , δ d k and δ τ such that κ < κ 0 , the state X(s) solution of the uncertain neutral system (17) exponentially converges to zero.

Proof. The characteristic equation associated to the system ( 17) is given by G u (s) -K(s)e -sδ = 0.

(18) We need to prove that the solutions of this equation are all located in the complex left half-plane if the delays and uncertainties are small enough. Equation (18) rewrites for s = 0 ān + bn e -sτ = -

1 s n n-1 k=0 (ā k + bk e -sτ )s k + 1 s n e -sδ ( p k=0 c k s k + p k=0 d k s k e -sτ ), (19)
Since Assumption 3 is satisfied, there exist κ > 0 such that if κ < κ, then | bn | < |ā n |. Thus (see [START_REF] Hale | Introduction to functional differential equations[END_REF]), the left part of ( 19) is lower-bounded by a constant η > 0 for s ∈ C + = {s ∈ C, (s) ≥ 0}. Due to Assumption 4 the right part of ( 19) converges to zero when |s| goes to infinity. Thus, there exists M 0 > 0 such that for all s ∈ C + with |s| > M 0 , the right part of ( 19) is smaller than η 2 . This implies that (18) does not have a solution on 20) Since the control law u(s) = K(s)x(s) stabilizes the nominal system, there exists η 1 > 0 such that for all s ∈ C + we have |G(s) -K(s)| > η 1 . Since |s| ≤ M 0 , the function K is bounded on the considered domain. Consequently, there exists δ max > 0 such that |K(s)(1e -sδ )| < η1 4 . Similarly, we have the existence of κ > 0 such that for κ < κ, |G(s) -G u (s)| < η1 4 . This implies that (20) Fig. 1. The spectral values' distribution for equation ( 5) is illustrated using QPmR toolbox from [START_REF] Vyhlidal | Mapping based algorithm for large-scale computation of quasi-polynomial zeros[END_REF]. We have a triple root at s = -3.

C + if |s| > M 0 . Let us now consider s ∈ C + such that |s| ≤ M 0 . Equation (18) rewrites G(s) -K(s) = (G(s) -G u (s)) -K(s)(1 -e -sδ ). (
does not have any solution on C + if |s| ≤ M 0 . Choosing κ 0 = min κ, κ concludes the proof.

This proposition can now be applied for the stabilization of torsional oscillations. Proposition 6. Consider the neutral system (5) with the control law u T (t) = κ 1 Φ L (t -τ ) + κ 2 ΦL (t -τ ). It this control law stabilizes the system (namely if all the roots of the quasipolynomial L(s) defined in (7) have negative real part), then it is robust with respect to uncertainties in the parameters and delays in the loop.

Proof. The control law u T (t) rewrites (using the Laplace transform) as u T (s) = κ 1 e -sτ + κ 2 se -sτ . The rest of the proof is a straightforward application of Proposition 5.

ILLUSTRATIVE EXAMPLE

In this section, we numerically illustrate our results. The values of the different parameters are defined in Table 1. In Figure 1, we have pictured the spectral values' distribution of the characteristic function ( 7) for the gains defined in Proposition 2. It is important to notice that the multiple root at s 0 = -3 is dominant, which guarantees that the closed loop-system is asymptotically stable. Finally, we have pictured in Figure 2, the distribution of the roots for the characteristic function (7) in presence of an uncertainty acting on the length of the well L. The red branches correspond to a positive uncertainty whose value vary between 0 and 30% of the real value of L. The blue branches correspond to a negative uncertainty. It is important noticing that such an uncertainty modifies all the coefficients a i (since α L directly depends on L) and the delay τ . As stated in Proposition 6, the system is robust with respect to this uncertainty, if it is small enough. measurement is available (which may not be the case for the considered drilling application), then our control law does not require an observer for practical implementation. Finally, this control law has been proved to be robust with respect to uncertainties and delays. This has been done by means of an analysis in the Laplace domain and validated through simulations.

  0 -27 a 1 + 81 a 2 -243 a 3 ) τ 2 2e 3 τ -(6 a 0 -54 a 2 + 324 a 3 ) τ 2e 3 τ --54 e 3 τ b 3 + 2 a 0 -54 a 3 2e 3 τ , (8)b 1 = -(3 a 0 -9 a 1 + 27 a 2 -81 a 3 ) τ 2 e 3 τ -(a 0 + 3 a 1 -27 a 2 + 135 a 3 ) τ e 3 τ --27 e 3 τ b 3 + a 1 -27 a 3 e 3 τ , (9) b 2 = -1/2 (a 0 -3 a 1 + 9 a 2 -27 a 3 ) τ 2 e 3 τ -1/2 (2 a 1 -12 a 2 + 54 a 3 ) τ e 3 τ -1/2 -18 e 3 τ b 3 + 2 a 2 -18 a 3 k we have āk = a k + δ a k , bk = b k + δ b k , ck = c k + δ c k , dk = d k + δ d k and τ = τ + δ τ , the terms δ a k , δ b k , δ c k , δ d kand δ τ representing small (unknown) constant uncertainties acting on the different parameters of the systems. We denote κ their common upper-bound.

Fig. 2 .

 2 Fig. 2. Sensibility of the spectral values for equation (5)with respect to an uncertainty acting on the length L (positive uncertainty in blue, negative in red)6. CONCLUDING REMARKSIn this paper, we have used a PDE model of the torsional oscillations of a drilling system. Following the approach given inBoussaada et al. (2019a), this system has been rewritten as a third-order neutral equation with a single delay. We then have designed a reduced complexity delayed stabilizing controller that enables an appropriate poleplacement of the closed-loop system. Compared toBoussaada et al. (2019a), this control law now include a derivative action, then introducing a new degree of freedom.The obtained controller acts for damping oscillation of the wave equation. The proposed control law only requires the measurement of the state at the boundary. If such a measurement is available (which may not be the case for the considered drilling application), then our control law does not require an observer for practical implementation. Finally, this control law has been proved to be robust with respect to uncertainties and delays. This has been done by means of an analysis in the Laplace domain and validated through simulations.

This research was partially supported by the iCODE Institute, research project of the IDEX Paris-Saclay in the "Programme des Investissements d'Avenir".

The inertia moment is such thatI B = M B r

, where r is taken as the averaged radius of drillpipe.