Jonathan Jaglin
email: jonathan.jaglin@lag.ensieg.fr

Carlos Canudas-De-Wit
email: carlos.canudas-de-wit@inpg.fr

Cyrille Siclet

Delta Modulation for Multivariable Centralized Linear Networked Controlled Systems

Keywords: Delta modulation, Networked controlled systems, NCS, quantized systems

 to any type of linear multivariable systems. In addition, we provide an estimation of the stability attraction domain, and we give some simulation results validating the proposed approach.

I. INTRODUCTION

T His paper deals with the stabilization problem of a linear multivariable system through a communication network where information is transmitted via a particular coding algorithm. Coding algorithms seeking to transmit a minimum number of information bits are appealing in wireless networks since they allows a substantial channel bandwidth reduction. Many of such types of control architecture using that type of codes have been studied in the past. See [START_REF] Ishii | Remote control of lti systems over networks with state quantization[END_REF], [START_REF] Elia | Stabilization of linear systems with limited information[END_REF], [START_REF] Liberzon | On stabilization of linear systems with limited information[END_REF], [START_REF] Lemmon | Control system performance under dynamic quantization: the scalar case[END_REF], [START_REF] Tan | Numerical study of joint quantization and control under block-coding[END_REF], [START_REF] Li | Robust quantization for diginal finite communication bandwidth (dfcb) control[END_REF], [START_REF] Hespanha | Towards the control of linear systems with minimum bit-rate[END_REF], [1], [START_REF] Tatikonda | Control under communications constaints[END_REF] among others.

Delta modulation (∆-M ) is one alternative to minimize the numbers of bits to be coded. The reason is that innovation increments (with a granularity depending on a quantization factor ∆) are coded rather than the absolute value of the signal. Recent works in [START_REF] Canudas-De-Wit | Differential coding in networked controlled linear systems[END_REF] have re-adapted the standard form of the delta modulation structure to their use in a feedback setup. One advantage of this type of strategy is that the coding algorithm can be built in a methodological and simple manner. A limitation is that re-synchronization may be needed, if the signal track is lost. Inspired by this approach several variants of [START_REF] Canudas-De-Wit | Differential coding in networked controlled linear systems[END_REF] have been studied: asynchronous entropy coding [START_REF] Canudas-De-Wit | Entropy coding in networked controlled systems[END_REF], energy-aware coding [START_REF] Canudas-De-Wit | Energy-aware 3-level coding and control co-design for sensor network systems[END_REF], adaptive delta modulation [START_REF] Rubio | Adaptive delta-modulation coding in networked controlled systems[END_REF], and gain scheduling multi-bit coding [START_REF] Abdallah | Compensation schemes for a delta-modulation-based ncs[END_REF]. Except for the trivial case of diagonalizable multivariable system that can be reformulated as a set of n-scalar ones, all these works deal exclusively with scalar system.

In this paper, we present a generalization of the deltamodulation coding presented in [START_REF] Canudas-De-Wit | Differential coding in networked controlled linear systems[END_REF], to MIMO systems. In particular we introduce a vector coding structure for multivariable centralized linear systems. The notion of centralization 1 MIMO : Multiple Input Multiple Output. refers here to the fact that both the encoder-decoder and the control law use the full available information from all sensors. The idea is shown in Fig. 1, where we can see that all the sensed system outputs are collected in a central point, then transformed into a different coordinate-basis (using the transform matrix T k ) before they are coded using a vectorcoding algorithm. At the receiver side, it is similarly assumed that the transmitted information arrives to a central receiver, then decoded, and finally the control is computed using this centralized information. It is worth to notice that decentralized case is clearly much more constrained, even in absence of a coding process. A recent work [START_REF] Basar | Decentralized multi controller stability over communication channels[END_REF] dealing with the case of decentralized multi-controller stability over communication channel illustrates well the fundamental difficulties, and provides an interesting preview on how to handle these problems when information is not centralized. The paper is organized as follows. After formalizing the problem in section II, we introduce in section III the general vector coding algorithm that can be adapted for all different forms of Jordan blocks resulting from the change of coordinate basis. Then, vector coding is performed in the transform domain. Vector coding here refers to the fact that a specific codeword is assigned for specific combinations between states. These are new information that are not present in standard scalar coding as they result from the combination of individual signals. Closed-loop stability properties resulting from this approach are also exposed here. Section V characterizes the attraction set associated to the previous local stability conditions. This allows a finer estimation of quantization values to be used in the coding process. Finally simulation results are shown in section VI.

II. PROBLEM FORMULATION AND ASSUMPTIONS

The problem considered here is the stabilization of a multivariable system in which sensor signals are centralized, and then transmitted through a digital communication link to the controller. At the controller side, the information is received in a unique point, and then decoding process provides the system n-dimensional estimated state, to be used for feedback. The coding design aims to achieve stability with a minimal information rate, thanks to a judicious coding strategy selection during the quantization step.

Let us assume the following:

• the coding process is centralized : a single encoder can be used to encode all the sensed states of the system, • the encoded information is transmitted through a noiseless perfect transmission channel. Hence possible impairments (delay, errors) due to the transmission are not considered, • information flow is unidirectional; the information is only transmitted from the encoder to the decoder, • the encoder and decoder clocks are assumed to be synchronized, and samples are assumed to occur at each T s . The following notations will be used:

• n is the state dimension that corresponds to the number of sensors, • m is the number of control inputs,

• x k = [x 1 k , . . . , x n k ] T ∈ R (n×1)
is the n-dimensional sensed state vector at instant kT s (each x i k corresponds to the i -th sensor) ; m×1) , is m-dimensional control input vector at instant kT s . The discretized system is described by:

• u k = [u 1 k , . . . , u m k ] T ∈ R (
x k+1 = Ax k + Bu k (1)
where A ∈ R (n×n) , and B ∈ R (n×m) . Moreover, the control law is given by

u k = -K xk (2)
with K ∈ R (m×n) such that the eigenvalues of A -BK are strictly lower in magnitude than 1. xk is an estimation of x k , and xk denotes the estimation error :

xk = x k -xk , (3)
and, more generally, for a given signal s k , ŝk represent an estimated value of s k and sk represent the error s k -ŝk . Without loss of generality, we suppose that system (1) is already expressed in its Jordan's form, such that A is of the form,

A =   J λ1 0 0 0 J λ l 0 0 0 J λγ   (4)
where we assume that there are α Jordan's blocks, of dimension µ l × µ l , with multi-valued real eigenvalue, and γ -α Jordan's blocks, of dimension 2µ l × 2µ l , with multi-valued complex conjugated eigenvalues.

For the multi-valued real eigenvalue case, the J λ l , for 1 l α, are of the form,

J λ l =     λ l 1 0 0 λ l 1 0 λ l 1 0 λ l     (5)
and, for the multi-valued complex conjugated eigenvalues, the J λ l , are, for all α + 1 l γ, of the form,

J λ l =     |λ l |R(θ l ) I 2 0 0 |λ l |R(θ l ) I 2 0 |λ l |R(θ l ) I 2 0 |λ l |R(θ l )     (6)
where λ l = |λ l |(cos(θ l ) + j sin(θ l )) describes the complex eigenvalues, with magnitude |λ l |, and angle θ l . R(θ l ) is the rotation matrix associated to the polar form adopted above, i.e.

R(θ

l ) = cos(θ l ) sin(θ l ) -sin(θ l ) cos(θ l ) (7)
Remark 1: It is worthwhile underlining the fact that µ l is not necessarily the multiplicity order of λ l since the eigenvalues λ l are not necessarily different from each other. Nevertheless, we have µ 1 + . . . + µ γ = n, n being the size of A.

The case where A is diagonal (µ l = 1), with real-valued eigenvalues has been treated in [START_REF] Canudas-De-Wit | Differential coding in networked controlled linear systems[END_REF]. In this paper we extend these results to the general case of multiple-valued, real and complex eigenvalues, with µ l possibly different from 1.

III. MULTIVARIABLE ∆-MODULATION CODING STRATEGY

In this section, we present the multivariable coding strategy. This strategy is inspired from the ∆-modulation algorithm studied previously in [START_REF] Canudas-De-Wit | Differential coding in networked controlled linear systems[END_REF] for the one-dimensional case. The n-dimensional case considered here does not result from the simple extension of the one-dimensional case, but requires a new vector coding strategy, and a particular change of coordinates (matrix T k ) for the multi-valued complex conjugated eigenvalue case. The role of the rotation matrix T k is to align the direction of the eigenvector (signal oscillation) to the vector quantizer block.

A. Principle of multivariable coding and decoding process

Figure 1 shows the architecture of the proposed differential coding algorithm. It is composed of three main components:

• The vector quantizer block transforms the error zk , into a finite codeword set, which is latter transformed into bits and sent through the communication channel, • The predictor, that transforms back the codeword into a system state prediction xk • The rotation matrix T k transforms the estimation error xk between the signal x k and its estimated (reconstructed) value xk into a new set of coordinates zk , i.e.

zk = T -1 k xk (8)
As shown in the Figure 1, the encoding algorithm has the 3 components described above, while the decoding algorithm is just the predictor whose inputs are the received information codewords at the decoding side.

Each of these components are explained in detail next. 1) vector quantizer: it maps the transformed vector zk into the quantized vector ẑk . The multi-level quantizer is constructed as follows:

• we consider M i (odd or even) subdivisions for each zi with respective quantization step ∆ i . The partition is centered at the origin, • This partition generates an hypercube of dimension n with a total of n C = i=n i=1 M i quantized volumes (see example in Figure 2), • To each quantized volume is associated a value for the quantized vector ẑk (see example in the Table I). The formula used to compute ẑk is the following:

If M i is odd, then ẑi k is given as:

ẑi k =    (M i -1)∆ i /2 if C 1 N ∆ i if C 2 -(M i -1)∆ i /2 if C 3
where the conditions C i are:

C 1 : zi k (M i -1)∆ i /2 C 2 : zi k ∈ [(N -1/2)∆ i , (N + 1/2)∆ i ], (N ∈ {-(M i -1)/2, . . . , (M i -1)/2}) C 3 : zi k < -(M i -1)∆ i /2 If M i is even, then ẑi k is given as: ẑi k =    M i /2∆ i if C 1 (N + 1/2)∆ i if C 2 -M i /2∆ i if C 3
where the conditions C i are:

C 1 : zi k (M i -1)/2∆ i C 2 : zi k ∈ [N ∆ i , (N + 1)∆ i ], (N ∈ {-(M i -1)/2, . . . , (M i -1)/2}) C 3 : zi k < -(M i -1)∆ i /2
Remark 2: Before transmission, note that quantizer vector ẑi k is associated to a codeword of dimension n C that can be coded directly into R = ⌈log 2 (n C )⌉ bits, where ⌈.⌉ denotes the ceil function.

2) Predictor: The estimation of the signal xk is computed thanks to a model-based predictor:

xk+1 = (A -BK)x k + A xk = (A -BK)x k + AT k ẑk (9)
where the last expression results from the use of the inverse transformation matrix, i.e.

xk = T k zk (10)
Due to the particular nature of this transformation (rotation matrix) its inverse always exists. Thus, using equations ( 8), ( 10) and ( 9), we get :

zk+1 = T -1 k+1 AT k (z k -ẑk ) (11)
Note that, as this predictor is used at both the encoder and the decoder side, their respective initial conditions x0 and ẑ0 are assumed to be the same.

3) Transformation matrix T k : The selection of this matrix for the general case is quite involved. In what follows we present two examples: one with a trivial choice of T k = I, and an other where its choice depends on the eigenvalues position in the complex plane. The general case will be treated in detail in section IV.

B. Example 1: two-dimensional system with a real eigenvalue

Consider a system of the form (1), with

A = λ 1 0 λ
and some B such that (A, B) is controllable. Then, as the system does not contains oscillatory modes, we can take T k = I 2 , where I n denotes the n-entry identity matrix, which leads, with xk = zk , to

zk+1 = λ 1 0 λ (z k -ẑk ) ∆ 1 ∆ 2 z1 z2 d 2 d 1 Ω ext Ω int 1 2 3 4 5 6 7 8
9

Fig. 2. Evolution of zk where z0 begins in

Ω ext = {[-1.5∆ 1 , 1.5∆ 1 ] × [-1.5∆ 2 , 1.5∆ 2 ]} and zk ∈ Ω int = {[-d 1 , d 1 ] × [-d 2 , d 2 ]
} and the dots delimit the nine subdivisions of the space.

Let us choose M i = 3 subdivisions per signal, with a different step for each one; a quantization step of ∆ 1 > 0 for z1 k , and ∆ 2 > 0 for z2 k . This partition is shown in Figure 2, and the associated coding strategy in Table III-B Now if we assume that |λ| < 3, and that the quantization steps are chosen such that

∆ 2 < ∆ 1 (3 -|λ|) ( 12
)
then it is easy to show that if the error signal z0 is initiated inside the centered rectangle set Ω ext , then the evolution of

k = 0 Ω ext (1) k = 1 k = 1 Ω ext (2) Ω int (2) T k = I θ = Π/4 Ω ext (2) Ω int wasted space T k = R(θ) a ) b ) c )
Fig. 3. Evolution of zk , in the first figure, we choose that zk ∈ Ω ext (1) and in the second figure we see that zk+1 ∈ Ω int and we see that if we code the signal zk+1 we loose some space and, to ensure that Ω int ⊂ Ω ext , the maximal possible eigenvalue is |λ| < 3/ √ 2. The third figure shows a forced rotation of the coder which permits to have better performances

TABLE I CODING STRUCTURE RELATED TO FIGURE 2 Bits Codeword Value of ẑk 0000 1 (∆ 1 , 0) 0001 2 (∆ 1 , ∆ 2 ) 0010 3 (0, ∆ 2 ) 0011 4 (-∆ 1 , ∆ 2 ) 0100 5 (-∆ 1 , 0) 0101 6 (-∆ 1 , -∆ 2 ) 0110 7 (0, -∆ 2 ) 0111 8 (∆ 1 , -∆ 2 ) 1000 9 (0, 0)
zk will enter (in one step) inside the set Ω int as defined in Figure 2.

To see that, note that if zk ∈ Ω ext , then we have

|z i k -ẑi k | ∆i 2
, ∀i ∈ {1, 2}. Now, from error equation in zk , we have that

|z 1 k+1 | < |λ| ∆1 2 + ∆2 2 = d 1 , and that |z 2 k+1 | < |λ| ∆2 2 = d 2 .
This defines the set Ω int . From here it is obviously needed that Ω int ⊂ Ω ext , which lead to the condition [START_REF] Lemmon | Control system performance under dynamic quantization: the scalar case[END_REF].

C. Example 2: two-dimensional system with complex conjugate eigenvalues

Consider a system of the form (1), with

A = |λ|R(π/4)
with R(π/4) is defined in [START_REF] Hespanha | Towards the control of linear systems with minimum bit-rate[END_REF], and B such that the pair (A, B) is controllable. Suppose that we take T k = I 2 , which gives xk = zk and from (3) we get zk+1 = |λ|R(π/4)(z k -ẑk ) As in the former example, let us choose M i = 3 subdivisions per signal, with a quantization step ∆ 1 > 0 for z1 k , and ∆ 2 > 0 for z2 k .

We suppose that the initial condition at k = 0 z0 ∈ Ω ext defined in the Figure 3 a), thus at k = 1 we obtain z1 ∈ Ω int (Figure 3

b)). It can be proved following similar steps as in

Example 1 that Ω ext is an invariant set if |λ| < M1 √ 2 with ∆ 1 = ∆ 2 .
This condition is more conservative than the one obtained in Example III-B, where we only require that |λ| < M 1 . It is also possible to retrieve the same result by redefining the transform matrix T k as shown below

Let us choose T k such that

T k = R(kπ/4)
Then zk = R(-kπ/4)x k with R(π/4) -1 = R(-π/4). Equation (3) becomes

zk+1 = R(-(k + 1)π/4)|λ|R(π/4)R(-kπ/4) -1 (z k -ẑk ) = |λ|R(-(k + 1)π/4)R(π/4)R(kπ/4)(z k -ẑk ) = |λ|I 2 (z k -ẑk )
Hence, we obtain a fully decoupled system and it is straight forward to show that if z0 begins in the set Ω ext , it is necessary that Ω int ⊂ Ω ext to ensure that Ω ext is an invariant set, this condition leads to |λ| < 3 and a independent choice of ∆ 1 and ∆ 2 . In this case, we see that we can find the same properties as in the real eigenvalues system. The generalization of this result needs an other transformation.

IV. CONSTRUCTION OF THE TRANSFORM MATRIX T k :

GENERAL CASE Consider a system of the form (1), with A defined in (4) and B such that (A, B) is controllable. The error equation:

xk+1 = A(x k -xk )
As we have assumed that A is a block diagonal matrix, the associated stability properties can be analyzed separately for J λ l . In the following paragraph, we will first deal with the case of real eigenvalues 1 l α and latter we will focus on the complex conjugate case α + 1 l γ.

To simplify the notation, we only note xk instead of xk (l) ∈ R µ l , J λ = J λ l and µ = µ l .

A. Case of multiple-valued real eigenvalues

Lemma 1: Case of multiple real eigenvalues. Assuming that ẑk is computed thanks to the quantization procedure given in section III-A1, and suppose that

z0 ∈ Ω ext = {z ∈ R µ : |z i | M i ∆ i 2 , 1 i µ}
and the quantization steps satisfy the equations

|λ| + ∆ i+1 ∆ i M i , 1 i µ -1 (13) Then i) Ω ext is an invariant set ii) zk ∈ Ω int , ∀k 1 where Ω int = {z ∈ R µ : |z i | |λ|∆ i /2 + ∆ i+1 /2 ∀i : 1 i µ -1 and |z µ | λ∆ µ /2}
Proof: According to (5):

zi k+1 = λ(z i k -ẑi k ) + (z i+1 k -ẑi+1 k ) (14)
zµ k+1 = λ(z µ k -ẑµ k ) (15)
Given that ẑi+1 k is quantized by the procedure given in section III-A1, we have

|z i+1 k -ẑi+1 k | ∆i+1 2 .
Then using ( 13), for 1 l µ -1, we get

|z i k+1 | |λ||z i k -ẑi k | + |z i+1 k -ẑi+1 k | |λ| ∆ i 2 + ∆ i+1 2
M i ∆ i 2 (16)
Finally, (13) implies that |λ| < M i , so that

|z µ k+1 | M i ∆ µ 2 (17)

B. Case of complex conjugate eigenvalues.

We now consider the case where λ ∈ C for α + 1 l γ. So, let us introduce the matrices W (θ) and Q(θ) defined by

W (θ) = R(θ) 0 0 R(θ) (18) Q(θ) = R(-θ) 0 0 R(-µθ) . ( 19
)
It can be shown after a few calculations that

Q -1 (θ)W -1 ((k + 1)θ)J λ W (kθ)Q(θ) =      |λ|I 2 I 2 0 0 0 |λ|I 2 I 2 0 . . . . . . |λ|I 2 0 . . . 0 |λ|I 2      = J λ Let us choose T k = W (kθ)Q(θ).
Then, as in the case of real-valued eigenvalues, we have

zk+1 = J λ (z k -ẑk ) (20)
and J λ is a block diagonal matrix, so that we can consider separately each block again. Then, considering separately even indices and odd indices, we exactly recover the results of the case of real-valued eigenvalues. Indeed, if we denote ze

k = [z 2 k , z4 k , . . . , z2µ k ] and zo k = [z 1 k , z3 k , . . . , z2µ-1 k ],
we have

ze k+1 =     |λ| 1 0 0 |λ| 1 0 |λ| 1 0 |λ|     (z e k -ẑe k ) ( 21
)
zo k+1 =     |λ| 1 0 0 |λ| 1 0 |λ| 1 0 |λ|     (z o k -ẑo k ) (22)
Lemma 2: Case of multiple complex eigenvalues. Assuming that ẑk is computed thanks to the quantization procedure given in section III-A1, and suppose that

z0 ∈ Ω ext = {z ∈ R 2µ : |z i | M i ∆ i 2 } |λ| + ∆ i+2 /∆ i M i , ∀i : 1 i 2µ -2 (23)
Then we ensure that i) Ω ext is an invariant set ii) zk ∈ Ω int , ∀k 1 where

Ω int = {z ∈ R 2µ : |z i | |λ|∆ i /2 + ∆ i+2 /2 1 i 2µ -2 and else |z i | |λ|∆ i /2}
Proof: The proof is identical to the one derived for the demonstration of Lemma 1 in the case of real-valued eigenvalues.

C. General case: combined real and complex eigenvalues

Theorem 1: Suppose the system (2)

x k+1 = Ax k + Bu k with the pair (A, B) controllable. And a channel rate R bounded by n l=1,|λ l |>1 ⌈|λ l |⌉ < 2 R
Then, the coding structure that ensure that x k is bounded, is realized with the Delta-modulation coding explained in section III where zk = T -1 k xk with T k defined as

T k =    I µ 1 0 0 0... 0 Iµ α 0 0... 0 0... W ι (kθ ι )Q ι (kθ ι ) 0 0 0... 0 W γ (kθγ )Q γ (θγ )    (24)
with α + 1 = ι. Then zk+1 = J (z k -ẑk ) and where

A = T k+1 J T -1 k J =       J 1 0 . . . 0 . . . J α . . . 0 . . . 0 J ι 0 0 . . . 0 J µ      
with the properties for M i and ∆ i given in lemma 1 for real eigenvalues and lemma 2 for complex eigenvalues.

Proof: For each signal with instable open loop, one of the condition is

|λ l | < M i , it is sufficient that ⌈|λ l |⌉ < M i with R = log 2 n i=1 M i .
If we multiply for all the coefficients, the result becomes

n l=1,|λ l |>1 ⌈|λ l |⌉ < 2 R
Using the previous lemmas, we ensure that x is bounded.

x k+1 = (A -BK)x k + Ax k
With the following system where A-BK has its eigenvalues strictly inferior than 1, the authors of [START_REF] Canudas-De-Wit | Differential coding in networked controlled linear systems[END_REF] have shown that the cascade system ensures that x k is bounded.

V. DOMAIN OF ATTRACTION AND NEW TUNING POLICIES

FOR ∆ i The aim of this section is twofold. First assuming the use of the tuning rule (13), we provide a less conservative method to estimate the attraction domain (named B ⊃ Ω ext ). Second, assuming the same attraction domain Ω ext , we provide a new tuning rule for the ∆ i that,compared to previous rule given in [START_REF] Basar | Decentralized multi controller stability over communication channels[END_REF], results in smaller values for ∆ i . As a consequence, the system precision can be improved. Specific simulation results concerning this last case, will be presented at the end of the paper.

A. Characterization of B

Let assume that the ∆ i are tuned following the rule in [START_REF] Basar | Decentralized multi controller stability over communication channels[END_REF], and denote B the new estimation of the attraction domain with Ω ext ⊂ B ⊂ R n . Let B be defined as the compositions of the sub-sets

B λ l , B = B λ1 × ... × B λw (25)
where the B λ l describes the attraction domain for the l-th Jordan's block, J λ l , under consideration,

zk+1 =   |λ l | 1 0 0 |λ l | 1 0 0 |λ l |   Jλ l (z k -ẑk )
This decomposition simplifies the analysis by looking at each block separately instead of considering the whole system together. Therefore, we only need to focus on a single block B λ l , and repeat the same analysis for other block when needed. Inspired by the Jordan block structure, assume in turn that B λ l = H λ l ,1 × ... × H λ l ,µ l where each subset, H λ l ,1 , correspond to a domain associated to each of the Jordan block components. For simplicity reasons, we omit the subindex λ l in the sequel. Hence, we simply note B = H 1 × ... × H µ l .

Theorem 2: Assume that ẑk is computed thanks to the quantization procedure given in section III-A1, and that ∆ i are tuned following the rule in [START_REF] Basar | Decentralized multi controller stability over communication channels[END_REF], and suppose that

z0 ∈ B = {z ∈ R µ : |z i | γ i } with, for 1 i µ -1, γ i = min (M -1)∆ i /2 + ε i max , (|λ|| ẑi k | -ε i+1 max )/(|λ| -1) ε i+1 max min (M -|λ|)∆ i /2, (M -1 -|λ|)∆ i /2 + ε i max then: i) B is an invariant set, i.e. zk ∈ B ∀k 0. ii) ∃k 1 > 0, such that, zk ∈ Ω int , ∀k k 1 .
where Ω int is the same set as defined in Lemma 1-ii).

Proof: Details of the proof are given in Appendix. Note that this analysis allows us to obtain a bigger attraction domain than the one obtained in section IV. To see this, note that ε i max ∆ i /2, which implies that γ i M ∆ i /2, and therefore we have that

B ⊃ Ω ext

B. Tuning policies for ∆ i

Assume now that the attraction domain Ωext , is given by

Ωext = {z ∈ R µ : |z i | δ i , 1 i µ}
where δ i are arbitrary values specified by the user. Note that the specification above imposes, in the previous tuning method, that M i ∆i 2 = δ i , whereas theorem 3 below will show that the new values ∆i < ∆ i = 2δi

Mi leading to a smaller convergence set Ωint ⊂ Ω int , where Ω int is the same set as defined in Lemma 1-ii).

Theorem 3: Suppose that z0 ∈ Ωext , and let the following rule to be applied to select the coding levels,for

1 i µ-1, ∆i = 2 |λ| -1 |λ|(M -1) δ i + 2 δ i+1 -(M -1) ∆i+1 /2 |λ|(M -1) ∆µ = δ µ (2(|λ| -1))/|λ| Then:
i) Ωext is an invariant set, and ii) ∃k 1 > 0, such that, zk ∈ Ωint , ∀k k 1 , where Ωint ⊂ Ω int is given as:

Ωint = z ∈ R µ : |z i | |λ| ∆i /2 + ∆i+1 /2 1 i µ -1 |z µ | |λ| ∆µ /2 i = µ
Proof: Property i) can be shown following the same proof as in part i) of Theorem 2, given in the Appendix. For the Property ii) the convergence of zk towards the set Ωint in finite time also follows the same lines as the proof of Theorem 2 and is omitted here.

Finally the fact that Ωint ⊂ Ω int follows by first observing that both sets Ωint , and Ω int have the same upper bound structure, and hence it is sufficient to prove that ∆i < ∆ i . This last inequality follows from inspection comparing the definition of the ∆i given in the theorem with the ones resulting from the imposed constraints to the previous tuning method, i.e. ∆ i = 2δi

Mi .

VI. SIMULATION RESULTS The aim of this section, is to compare by simulations, the precision improvements that the second tuning method derived in previous section can provide. For this, we consider a second order system, already in it Jordan form:

A = 1.1 1 0 1.1 , B = 0.2 0.3
The controller is designed on the basis of a full static state feedback with the desired closed-loop eigenvalues located at (0.5, 0.6). The control objective is to regulate the output states to a fix value;

x ref 1 = 1, x ref 2 = 1.
The desired attraction domain for the estimation error is specified as (δ 1 , δ 2 ) = (0.62, 0.52), and the initial error state are taken inside Ωext ; x 0 = (0.6, 0.5), and x0 = (0, 0). We choose 2 word-code by signal namely 2 bits per unit of sampling time;

M 1 = M 2 = 2.
Under this conditions on the quantization step ∆ i are computed according to the conditions given in theorem 1; ∆ 1 = 0, 62 and ∆ 2 = 0, 52. The ∆i are now computed following the procedure in Section V; ∆1 = 0, 35 and ∆2 = 0, 057. k (upper) and x 1 k (lower) using two different tuning methods discussed in this paper. The impact of quantization on the first state is less effective than on the second state.

Figure 4 shows the time-evolution of the resulting closedloop signals (coding including). In both runs, the initial condition are the same, and as it was expected the second methods provides smaller values for the coding gains, resulting in a better signal reconstruction quality, and hence better regulation precision.

VII. CONCLUSION

In this paper, we have investigated the closed-loop properties of multivariable (MIMO) linear systems where the sensed information is centralized and coded on the basis of a ∆modulation algorithm intended to be used for minimizing the number of transmitted bits.

In particular we had proposed a new centralized vector coding algorithm that allows us to extend our previous results in [START_REF] Canudas-De-Wit | Differential coding in networked controlled linear systems[END_REF] to linear multivariable systems of arbitrarily dimension and arbitrarily structure (any canonical form with arbitrarily eigenvalues). The key feature allowing this results was based on the idea of performing the differential coding in a timevarying rotation coordinates associated to the well known canonical Jordan forms.

We have also shown that this fixed-gain simple and methodic coding strategy results in a ultimately uniformly (local) stability. We have also provided an estimation of the attraction domain, and a new method to tune the coding gains, resulting in closed-loop precision improvements. Simulation results have also been presented validating the proposed approach.

Future extensions of this work envision to devise adaptation rules for the coding gains, in order to generalize these results to global stability with an arbitrarily small convergence set precision.

Fig. 1 .

 1 Fig. 1. NCS System with γ representing conversion from ẑ to codeword.

Fig. 4 .

 4 Fig. 4. Time-evolution of the closed-loop state x 2k (upper) and x 1 k (lower) using two different tuning methods discussed in this paper. The impact of quantization on the first state is less effective than on the second state.

APPENDIX

The principle of this demonstration is based on a cascade argument. We assure that z1 k remains in H 1 under conditions on z2

k and so on. In section IV, we have studied the comportment of zk in Ω ext . Here we only interest us in the case where zk is outside Ω ext , so we specify all the demonstration on the last quantization namely

Here we determine conditions on z1 k and z2 k which permits to assure that the upper bound of H 1 is bigger than (or equal to) M ∆ 1 /2. We characterize the maximal value

, which yet permits that z1 k remains in H 1 . To this purpose, let us introduce the function

max )/(|λ| -1). In the following analysis, we use a more restrictive condition to characterize

On this part, we interest ourselves in the zone where the rate of the function is positive and then negative. We search to find the maximal value of ϕ(z

This function is piecewise described on each interval of step ∆ 1 /2 and we easily obtain that the function ϕ 1 is totally increasing (decreasing) on each sector, hence the maximal values are at each extremity. So the value is

To ensure that z1 k remains bounded, when the rate ∆V 1 k is positive, the maximal value of z1 takes at time k + 1 has to be lower than α 2 . So we obtain:

Moreover, if we fix the worst case for the size of

This choice permits a bigger size for H 2 due to the inequality between the two quantization steps in (13), we have ε 2 max ∆ 2 /2, so we lose in conservatism and we obtained |z 2 k -ẑ2 k | ε 2 max . Here we conclude that

The analysis on H i is almost the same as the analysis on H 1 except that another constraint is imposed ∀k k > 0 |z i k | | ẑi k | + ε i max and we obtain:

), we can conclude that γ i M ∆ i /2 and we have

We have determined conditions that zµ k needs to fill in order to characterize H µ-1 . Let us introduce the function

There exists a set where

In the following, we obtain a restrictive condition. If that zone exists, we must interest ourselves in the zone

where the rate is positive. To understand its impact let us introduce ϕ µ .

With the analysis of ϕ µ we can show that on the zone where ∆V µ k is positive, we ensure that the worst case of zµ k+1 is inferior than min(ε µ max , |λ|(M -1) 2(|λ|-1) )∆ µ . So we obtain that γ µ = min(ε µ max , |λ|(M -1)∆ µ /(2(|λ| -1))). We conclude that


Stefan Ankirchner
email: s.ankirchner@uni-jena.de

Christophette Blanchet-Scalliet
email: christophette.blanchet@ec-lyon.fr

Diana Dorobantu
email: diana.dorobantu@univ-lyon1.fr

Laura Gay
email: laura.gay@univ-lyon3.fr

First passage time density of an Ornstein-Uhlenbeck process with broken drift

Keywords: Ornstein-Uhlenbeck process, piecewise-linear drift, first passage time density, Laplace transform, joint law of the running supremum and endpoint, Fokker-Planck equation

We consider an Ornstein-Uhlenbeck process with different drift rates below and above zero. We derive an analytic expression for the density of the first time, where the process hits a given level. The passage time density is linked to the joint law of the process and its running supremum, and we also provide an analytic formula of the joint density / distribution function. Results from a numerical experiment reveal that our formulas allow to numerically evaluate the joint law and the density of the first passage time faster than a simulation based method.

Introduction

In this article we consider an Ornstein-Uhlenbeck (OU) process with a drift coefficient that is linear above and below zero, but with different rates. The drift coefficient is thus a piecewise linear function with a break point at zero. The aim is to determine the law of the first time the process hits a given level, and the joint law of the process and its running supremum. For the classical Ornstein-Uhlenbeck process with a globally linear drift these laws are already known: the distribution of first hitting times is computed in [START_REF] Alili | Representations of the first hitting time density of an Ornstein-Uhlenbeck process[END_REF] and [START_REF] Linetsky | Computing hitting time densities for CIR and OU diffusions: Applications to mean-reverting models[END_REF], and the joint law of the process and its running supremum in [START_REF] Blanchet-Scalliet | Joint law of an Ornstein-Uhlenbeck process and its supremum[END_REF].

Ornstein-Uhlenbeck processes with piecewise linear drift coefficients have hardly been examined yet. If the slope is always positive, then the OU process has a stationary density, explicitly computed, e.g. in [START_REF] Browne | Piecewise-linear diffusion processes[END_REF]. Apart from that, little is known, in particular about first hitting times. We show that if there is only one kink in the drift, then one can derive an analytic formula for the density of the first time the process attains a fixed level. Moreover, one can explicitly describe the joint law of the process and its running supremum.

First hitting times are intensively studied within applied probability theory for various reasons. They can be used, for example, for a probabilistic representation of parabolic PDEs with finite boundary conditions. In economic applications first hitting times are frequently used for modeling the arrival time of events, e.g. death times or bankruptcy times. In queuing theory diffusion hitting times can be used for approximating the first time a queue attains a certain length.

The recent years have seen an increased interest in diffusion processes with piecewise linear drift and piecewise constant diffusion coefficients. Such processes naturally arise, for example, when considering diffusion approximations of some queues (see [START_REF] Browne | Piecewise-linear diffusion processes[END_REF]). Moreover, in many control problems the optimally controlled state is a diffusion with a piecewise linear drift or piecewise constant diffusion coefficient, respectively. In particular, the OU process with a single kink in its drift can be interpreted as an optimally controlled state process within a model where an agent aims at maximizing the process's occupation time above zero (see Section 4 in [START_REF] Gay | Processus d'Ornstein-Uhlenbeck et son supremum: quelques résultats théoriques et application au risque climatique[END_REF]). A Brownian motion with piecewise constant diffusion coefficients appears in optimal diffusion control problems with piecewise convex and concave payoff functions (see [START_REF] Mcnamara | Optimal control of the diffusion coefficient of a simple diffusion process[END_REF]).

Several authors have considered some optimization problems within models using diffusions with piecewise linear drift and piecewise constant diffusion coefficients. For example in [START_REF] Mordecki | Optimal stopping of oscillating Brownian motion[END_REF] and [START_REF] Mordecki | Optimal stopping of brownian motion with broken drift[END_REF], the authors solve optimal stopping problems for a Brownian motion with positive piecewise constant volatility ( [START_REF] Mordecki | Optimal stopping of oscillating Brownian motion[END_REF]) or a positive piecewise constant drift ( [START_REF] Mordecki | Optimal stopping of brownian motion with broken drift[END_REF]) changing at zero. Another optimal stopping problem related to the pricing of an American option in an extension of the Black-Merton-Scholes model with piecewise constant dividend and volatility rates appears in [START_REF] Gapeev | Perpetual American options in a diffusion model with piecewise-linear coefficients[END_REF].

A statistical perspective is taken in [START_REF] Lejay | Maximum likelihood drift estimation for a threshold diffusion[END_REF], analyzing a maximum likelihood estimator for diffusions with piecewise constant drift and diffusion coefficients. The benefit of diffusions with piecewise constant parameters in financial applications is demonstrated in [START_REF] Guterding | The Heston stochastic volatility model with piecewise constant parameters-efficient calibration and pricing of window barrier options[END_REF]. The authors consider a version of the Heston model with piecewise constant parameters and discuss its calibration to the foreign exchange options market. We obtain our results by studying the Fokker-Planck equation associated to the joint law of the process and its running supremum. We explicitly solve the Fokker-Planck equation with the spectral method, and thus obtain an infinite series representation of the joint density/ distribution function of the process and its supremum at a given time t ∈ [0, ∞). A rather short calculation then allows to derive, from the joint law, a density formula for first passage times.

We remark that an alternative way for determining hitting time densities is to compute first their Laplace transform and then to use a Laplace inversion theorem. Indeed, the Laplace transform of the first passage time of the OU with a single kink in its drift can be calculated explicitly (see Section 4 in [START_REF] Gay | Processus d'Ornstein-Uhlenbeck et son supremum: quelques résultats théoriques et application au risque climatique[END_REF]). From this one can, formally, derive an analytic formula for the density of the first passage time. The Laplace transform, however, has infinitely many poles and thus the inversion requires to determine the residues along these poles. It seems to us that is is difficult to check whether the conditions of the residue theorem are satisfied, and therefore we do not follow this path.

The paper is organized as follows. In Section 2 we present the expressions of the density for the first passage time of an Ornstein-Uhlenbeck process with broken drift. We also describe the joint density / distribution function of the process and its running supremum. Moreover, we present some results from numerical experiments. Section 3 is devoted to the proof of our main results. For the reader's convenience, in Section 4 we provide some properties and auxiliary results of Hermite functions and parabolic cylinder functions.

Results

Ornstein-Uhlenbeck process with broken drift

Let (B t ) t≥0 be a standard Brownian motion defined on a probability space (Ω, F, P). Moreover let k -, k + ∈ R * + and let k : R → R be the function defined by

k(x) = k -, if x ≤ 0, k + , if x > 0. (1)
Let x 0 ∈ R and denote by X = (X t ) t≥0 the strong solution of the stochastic differential equation (SDE)

dX t = -k(X t )X t dt + dB t , X 0 = x 0 . (2)
Note that the drift coefficient of X is a piecewise linear function with a kink in zero. In accordance with [START_REF] Mordecki | Optimal stopping of brownian motion with broken drift[END_REF] we refer to X as an Ornstein-Uhlenbeck process with a broken drift.

Let m ∈ R be fixed. We introduce the first passage time of (X t ) t≥0 of the level m:

τ = inf{t > 0, X t = m}. (3)
Note that if x 0 < m ≤ 0 or x 0 > m ≥ 0, then τ is the first passage time of a "classical" OU process with drift function x → k -x and x → k + x, respectively. The only cases requiring new computations are (x 0 ≤ m, m > 0) and (x 0 ≥ m, m < 0). Symmetry allows to reduce the second case to the first one. Therefore, in the following we consider only the case (x 0 ≤ m, m > 0).

Main theorems

Our first aim is to prove that τ has a density p τ and to derive an analytic formula for it. We introduce the following notations: for any µ ∈ R let D µ be the parabolic cylinder function of the parameter µ (see Appendix 4.1 for the definition) and

α µ = 1 2 D µ k - (0)
D µ k + (0) + k + k - D µ k - -1 (0) D µ k + -1 (0) , β µ = 1 2 D µ k - (0)
D µ k + (0) - k + k - D µ k - -1 (0) D µ k + -1 (0) . Theorem 1. We assume that x 0 ∈] -∞, m], m ∈ R * + . Then 1.
The probability P(τ ∈ dt) admits a density p τ , 2. The density p τ satisfies, for all t ∈ R * + ,

p τ (t) =              -e k - 2 x 2 0 - k + 2 m 2 ∞ n=1 e -µnt D µn k - (-x 0 √ 2k -) ∂µ αµD µ k + (-m √ 2k + )+βµD µ k + (m √ 2k + ) µ=µn if x 0 ∈ R -, -e k + (x 2 0 -m 2 ) 2 ∞ n=1 e -µnt αnD µn k + (-x 0 √ 2k + )+βnD µn k + (x 0 √ 2k + ) ∂µ αµD µ k + (-m √ 2k + )+βµD µ k + (m √ 2k + ) µ=µn if x 0 ∈]0, m], (4)
where α n = α µn , β n = β µn , and (µ n ) n≥1 is the ordered sequence of positive zeros of

µ → α µ D µ k + (-m 2k + ) + β µ D µ k + (m 2k + ).
We prove Theorem 1 by using an explicit expression for the joint density/distribution of the endpoint and the running supremum, more precisely for

P(X t ∈ dx, τ > t) = P(X t ∈ dx, sup u≤t X u < m).
The analytic formula for the joint density/distribution is given in the following theorem.

Theorem 2. We assume that

x 0 ∈] -∞, m], m ∈ R * + . Then 1. For all t ∈ R * + the probability P X t ∈ dx, sup 0≤u≤t X u ≤ m admits a density p with support in ] -∞, m]: P X t ∈ dx, sup 0≤u≤t X u ≤ m = p(m, x, t) dx.
2. The density p solves the following Fokker-Planck equation:

               ∂ t p(m, x, t) = k -∂ x (xp(m, x, t)) + 1 2 ∂ 2 x 2 p(m, x, t), (x, t) ∈] -∞, 0[×R * + , ∂ t p(m, x, t) = k + ∂ x (xp(m, x, t)) + 1 2 ∂ 2 x 2 p(m, x, t), (x, t) ∈]0, m] × R * + , p(m, 0 -, t) = p(m, 0 + , t), t ∈ R * + , ∂ x p(m, 0 -, t) = ∂ x p(m, 0 + , t), t ∈ R * + , p(m, m, t) = 0, t ∈ R * + , lim t↓0 p(m, x, t)dx = δ x 0 (dx).
(5)

With the notations of Theorem 1 and for

(m, x) ∈ {(z, y) ∈ R * + ×R | z ≥ y}, t ∈ R *
+ , the density p can be expressed as:

• If x 0 ∈ R -, then p(m, x, t) =      e k - 2 (x 2 0 -x 2 ) ∞ n=1 e -µnt N (µn) D µn k - (-x 2k -)D µn k - (-x 0 2k -) , x ≤ 0 e k - 2 x 2 0 - k + 2 x 2 ∞ n=1 e -µnt N (µn) D µn k - (-x 0 2k -) α n D µn k + (-x 2k + ) + β n D µn k + (x 2k + ) , 0 < x ≤ m (6) • If x 0 ∈]0, m], then p(m, x, t) =            e k + 2 x 2 0 - k - 2 x 2 ∞ n=1 e -µnt N (µn) D µn k - (-x 2k -) α n D µn k + (-x 0 2k + ) + β n D µn k + (x 0 2k + ) , x ≤ 0 e k + 2 (x 2 0 -x 2 ) ∞ n=1 e -µnt N (µn) α n D µn k + (-x 2k + ) + β n D µn k + (x 2k + ) × α n D µn k + (-x 0 2k + ) + β n D µn k + (x 0 2k + ) , 0 < x ≤ m (7)
where

N (µ n ) = 1 2 ∂ x α n D µn k + (-x 2k + ) + β n D µn k + (x 2k + ) x=m × ∂ µ α µ D µ k + (-m 2k + ) + β µ D µ k + (m 2k + ) µ=µn . (8)
From Theorem 1 and 2 we can derive the law of the running supremum alone.

Corollary 3. Let x 0 ∈] -∞, m], m ∈ R * + . Then P sup 0≤u≤t X u ≤ m =              -e k - 2 x 2 0 - k + 2 m 2 ∞ n=1 e -µnt µn D µn k - (-x 0 √ 2k -) ∂µ αµD µ k + (-m √ 2k + )+βµD µ k + (m √ 2k + ) µ=µn if x 0 ∈ R -, -e k + (x 2 0 -m 2 ) 2 ∞ n=1 e -µnt µn αnD µn k + (-x 0 √ 2k + )+βnD µn k + (x 0 √ 2k + ) ∂µ αµD µ k + (-m √ 2k + )+βµD µ k + (m √ 2k + ) µ=µn if x 0 ∈]0, m]. (9)

Numerical results

Even though Formulas (4), (6) (or [START_REF] Gay | Processus d'Ornstein-Uhlenbeck et son supremum: quelques résultats théoriques et application au risque climatique[END_REF]) are infinite sums of special functions, it is easy to evaluate them numerically.

In a numerical experiment we have evaluated the functions t → p τ (t) and x → p(m, x, t) first with expression ( 4) and (6) (or ( 7)) and, second, with a Monte Carlo method using 1 × 10 6 simulations and ∆ t = 10 -6 in the Euler Scheme. Figures 1 and2 plot, respectively, p τ (t) and p(m, x, t) for some parameters and initial laws with T = 1.

The experiments reveal that the computation with our formula is faster than with a Monte Carlo method. Indeed, for the density depicted in Figure 1, the evaluation with our formula only lasts 4.09 seconds time, whereas the Monte-Carlo method requires more than 16 hours, using a processor 3.6 GHz with 4 cores.

Proofs

This section is devoted to the proofs of Theorem 2, Theorem 1 and Corollary 3.

Proof of Theorem 2

1. The existence is given by Lemma 2.1. of [START_REF] Delarue | First hitting times for general non-homogeneous 1d diffusion processes: density estimates in small time[END_REF] .

2. Directly from Proposition 5.4.3.1. of [START_REF] Jeanblanc | Mathematical methods for financial markets[END_REF] .

3. We use a spectral method to solve the Fokker-Planck equation.

Let us note D(L

F P ) the set {f ∈ L 2 (]-∞, m], e k(.). 2 /2 )|f ∈ L 2 (]-∞, m], e k(.). 2 /2 ), f (m) = 0, f (0 + ) = f (0 -), f (0 + ) = f (0 -), x → -k -(xf (x)) - 1 2 f (x) ∈ L 2 (]-∞, 0[, e k -. 2 /2 ), x → -k + (xf (x)) - 1 2 f (x) ∈ L 2 (]0, m[, e k + . 2 /2 )}
The Fokker-Plank operator L F P is defined by

L F P : D(L F P ) -→L 2 (] -∞, m], e k(.). 2 /2 ) f -→ L F P f : x → -(k(x)xf (x)) - 1 2 f (x)
Our goal is to find a basis (e n ) n≥1 of eigenvectors of L F P in order to write

p(m, x, t) = ∞ n=1 c n (t)e n (x).
Hilbertian basis of L 2 (] -∞, m], e k(.). 2 /2 ) Let remark that L F P is Hermitian. We prove that any function of L 2 (]-∞, m], e k(.). 2 /2 ) can be decomposed on an Hilbertian basis of eigenvectors of L F P . Since, it is not easy to prove it directly, we introduce a new operator T defined as

T : L 2 (] -∞, m], e k(.). 2 /2 ) -→L 2 (] -∞, m], e k(.). 2 /2 ) g -→f g
where f g is solution of

( * ) =                -1 2 f g -k -(xf g ) = g on ] -∞, 0[ -1 2 f g -k + (xf g ) = g on ]0, m[ f g (m) = 0 f g (0 -) = f g (0 + ) f g (0 -) = f g (0 + ) f g , f g ∈ L 2 (] -∞, m], e k(.). 2 /2 )
The eigenvectors of T are the same as those of L F P with inverse eigenvalues. Since T and L F P are Hermitian, we can use classical results on Hermitian operators decomposition. To prove that T is well defined we write it as an integral operator.

Lemma 1. The solution of ( * ) is

f g (x) = 2 m -∞ K(x, z)g(z)dz, where K(x, z) = 1 x<0 m x∨z e k -(y 2 -x 2 ) dy + 1 0≤x≤m m x∨z e k + (y 2 -x 2 ) dy.
Proof. We remark that f g (m) = 0. We easily check that f g satisfies ( * ):

f g (x) = 2 m -∞ ∂ ∂x K(x, z)g(z)dz with ∂ x K(x, z) = -2k -xK(x, z)1 x<0 -2k + xK(x, z)1 0≤x≤m -1 x>z 1 x<0 -1 x>z 1 0≤x≤m -1x>z
Hence

f g (x) = -2 x -∞ g(z)dz -2k -xf g (x)1 x<0 -2k + xf g (x)1 0≤x≤m and f g is con- tinuous in 0. Therefore f g is also continuous in 0. The second derivative of f g is f g (x) = -2g(x)-2k -(xf g (x)) 1 x<0 -2k + (xf g (x)) 1 0≤x≤m , which concludes the proof.
To obtain that T is well defined, continuous and compact, we first prove that the kernel K is square-integrable. More precisely, one has

Lemma 2. The kernel K satisfies K ∈ L 2 (] -∞, m]×] -∞, m], e k(x)x 2 -k(z)z 2 2
).

Proof. We have that .

m -∞ m -∞ e k(x)x 2 -k(z)z 2 K 2 (x, z)dxdz = 0 -∞ 0 -∞ e -k -x 2 -k -z 2 m x∨z e k -
All we need to show now is that

I 1 = 0 -∞ 0 -∞ e -k -x 2 -k -z 2 m x∨z e k -y 2 dy 2 dzdx is finite. Note that I 1 = 2 0 -∞ e -k -x 2 x -∞ e -k -z 2 dz m x e k -y 2 dy 2 dx = C 0 -∞ e -k -x 2 Φ( 2k -x) m x e k -y 2 dy 2 dx,
where C is a constant. For some constants C 1 , C 2 , we have that Gaussian cdf

Φ( 2k -x) ∼ x→-∞ C 1 e -k -x 2 x and m x e k -y 2 dy 2 ∼ x→-∞ C 2 -e k -x 2 2k -x 2
, hence I 1 < ∞ and the conclusion holds.

Classical results on operators defined by a kernel imply

2 (] -∞, m], e k(.). 2 /2 ) → L 2 (] -∞, m], e k(.). 2 /2 )
is well defined, continuous and compact.

(b) L 2 (] -∞, m], e k(.). 2 /2 ) admits an Hilbertian (countable) basis of eigenvectors of T for the scalar product associated to the L 2 (] -∞, m], e k(.). 2 /2 )-norm.

Eigenvectors of L F P

Let φ be an eigenvector of L F P with eigenvalue µ. Then, by the definition of L F P , we have for x ≤ m

L F P φ(x) = µφ(x) ⇐⇒ -(k(x)xφ(x)) - 1 2 φ (x) = µφ(x). If φ ∈ L 2 (] -∞, m], e k(.). 2 /2 ), then e k -x 2 2 φ(x)1 x≤0 + e k + x 2 2 φ(x)1 0<x≤m = φ(x) ∈ L 2 (] -∞, m]).
Hence for x ∈ R -:

- 1 2 (e -k -x 2 2 φ(x)) -k -(xe -k -x 2 2 φ(x)) = µe -k -x 2 2 φ(x) ⇐⇒ - 1 2 φ (x) + 1 2 k 2 -x 2 φ(x) = µ + k - 2 φ(x).
With the variable change z = x 2k -and h(z) = φ z √ 2k - = φ(x), we obtain the following equation satisfied by h:

h (z) + - 1 4 z 2 + µ k - + 1 2 h(z) = 0, z ≤ 0. Since h is integrable on -∞, then the solution is h(.) = D µ k - (-.). We deduce that for x ∈ R -, φ(x) = e -k 2 -x 2 2 D µ k - (-x 2k -).
For x ∈]0, m], it is easy to adapt the previous reasoning and there exist α, β ∈ R such that

φ(x) = e -k 2 + x 2 2 αD µ k + (-x 2k + ) + βD µ k + (x 2k + ) on ]0, m]. Since φ ∈ D(L F P ), it is a C 1 -function in 0. Therefore α = α µ = 1 2 D µ k - (0) D µ k + (0) + k + k - D µ k - -1 (0) D µ k + -1 (0) β = β µ = 1 2 D µ k - (0) D µ k + (0) - k + k - D µ k - -1 (0) D µ k + -1 (0) . Moreover φ(m) = 0, i.e. α µ D µ k + (-m 2k + ) + β µ D µ k + (m 2k + ) = 0.
In the following we denote by

(µ n ) n≥1 the ordered sequence of zeros of µ → α µ D µ k + (-m 2k + )+ β µ D µ k + (m 2k + ). Proposition 5. The roots (µ n ) n≥1 of α µ D µ k + (-m 2k + ) + β µ D µ k + (m 2k + ) = 0 are strictly positive. Proof. Since µ n is a zero of µ → α µ D µ k + (-m 2k + ) + β µ D µ k + (m 2k + ), this is equivalent to a zero of µ → α µ H µ k + (-m 2k + ) + β µ H µ k + (m 2k + ). So, using (13) 1 Γ(-µn k + ) ∞ i=0 1 i! Γ i -µn k + 2 (2m 2k + ) i α n + (-2m 2k + ) i β n :=A i = 0.
Remark that

A 2i = (2m 2k + ) 2i (α n + β n ) = (2m 2k + ) 2i × D µn k - (0) D µn k + (0) , A 2i+1 = (2m 2k + ) 2i+1 (α n -β n ) = (2m 2k + ) 2i+1 × k + k - D µn k - -1 (0) D µn k + -1 (0)
.

Using a proof by contradiction, let us assume that

µ n < 0 (hence µn k + , µn k -, µn k + -1, µn k -- 1 < 0), then Γ(-µn k + ) > 0, Γ i-µn k + 2 > 0, D µn k - (0) D µn k + (0) > 0 and D µn k - -1 (0) D µn k + -1 (0) > 0.
Therefore the previous sum is strictly positive (all the terms are strictly positive). We deduce that µ n ≥ 0.

Suppose now that µ n = 0. Since H 0 (x) = 1, then µ n will be solution of α n + β n = 0. This is absurde because when µ n = 0, α n + β n = 1. We deduce that µ n > 0.

Expression of p

We have found the eigenvectors of L F P namely the functions

x → e -k 2 -x 2 2 D µn k - (-x 2k -)1 x≤0 +e -k 2 + x 2 2 α n D µn k + (-x 2k + ) + β n D µn k + (x 2k + ) 1 0<x<m with (µ n ) n≥1 being the ordered positive roots of µ → α µ D µ k + (-m 2k + )+β µ D µ k + (m 2k + ).
With a normalization, we finally have an Hilbertian basis (e n ) n≥1 of L 2 (]-∞, m], e k(.). 2 /2 ):

e n : x → e -k 2 -x 2 2 D µn k - (-x 2k -)1 x<0 + e -k 2 + x 2 2 α n D µn k + (-x 2k + ) + β n D µn k + (x 2k + ) 1 0<x<m 0 -∞ D 2 µn k - (-y 2k -)dy + m 0 α n D µn k + (-y 2k + ) + β n D µn k + (y 2k + ) 2 dy .
The following lemma gives another equivalent form of the denominator of e n . Lemma 3. Under the assumptions of Theorem 2, we have:

0 -∞ D 2 µn k - (-y 2k -)dy + m 0 α n D µn k + (-y 2k + ) + β n D µn k + (y 2k + ) 2 dy = 1 2 ∂ x α n D µn k + (-x 2k + ) + β n D µn k + (x 2k + ) x=m ∂ µ α µ D µ k + (-m 2k + ) + β µ D µ k + (m 2k + ) µ=µn .
This term is denoted N (µ n ).

Proof. Using the same reasoning as [START_REF] Blanchet-Scalliet | Joint law of an Ornstein-Uhlenbeck process and its supremum[END_REF] (see the proof of Proposition 6), we can prove that if f : R 2 → R is solution of the following partial differential equation

∂ 2 x 2 f (x, µ) + (-x 2 k 2 + k + 2µ)f (x, µ) = 0, then for all x, µ, a, b ∈ R, b a f 2 (x, µ)dx = - 1 2 [f (b, µ)∂ x ∂ µ f (b, µ) -∂ x f (b, µ)∂ µ f (b, µ) -f (a, µ)∂ x ∂ µ f (a, µ) + ∂ x f (a, µ)∂ µ f (a, µ)] .
We apply this result for the two integrals. In the first case, we take

f (x, µ) = D µ k - (-x 2k + ) and for the second integral f (x, µ) = α µ D µ k + (-x 2k + )+β µ D µ k +
(x 2k + ). We conclude using the relations between α n and β n .

We look for a solution p of (5) that belongs to L 2 (] -∞, m], e k(.). 2 /2 ). We can decompose p on the basis (e n ) n≥1 :

p(m, x, t) = ∞ n=1 c n (t)e n (x). However p(m, x, 0) = ∞ n=1 c n (0)e n (x), then for n ∈ N * : c n (0) =< p(m, ., 0), e n > L 2 (]-∞,m],e k(.). 2 /2 ) = 0 -∞ e k -z 2 p(m, z, 0)e n (z)dz + m 0 e k + z 2 p(m, z, 0)e n (z)dz = e k -x 2 0 2 D µn k - (-x 0 2k -)1 x 0 ≤0 + e k + x 2 0 2 [α n D µn k + (-x 0 2k + ) + β n D µn k + (x 0 2k + )]1 0<x≤m N (µ n )
It remains to find c n (t) for t > 0. On the one hand

∂ t p(m, x, t) = ∞ n=1 c n (t)e n (x)
and on the other hand

∂ t p(m, x, t) = -L F P p(m, x, t) = - ∞ n=1 c n (t)L F P e n (x) = - ∞ n=1 c n (t)µ n e n (x).
The decomposition is unique, thus c n (t) = -c n (t)µ n and c n (t) = c n (0)e -µnt .

Proof of Theorem 1

First of all recall that if (Y t ) t≥0 is a Brownian Motion (BM), y ∈ R , Py the law of (Y t ) t≥0 when Y 0 = y and (F t ) t≥0 the filtration generated by the BM, then Girsanov Theorem gives: 1. The proof of the existence of p τ is similar to [START_REF] Pauwels | Smooth first-passage densities for one-dimensional diffusions[END_REF]. Let h : R → R a function defined by h(x) = -x 0 k(s)sds. The derivative of h is absolutely continuous. Ito formula gives (we can apply it, see Problem 7.3 , Section 3.7 [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF])

∀t > 0, dP y d Py Ft = exp - t 0 k(Y s )Y s dY s - 1 2 t 0 (k(Y s ))
h(X t ) = h(X 0 ) - t 0 k(X s )X s dX s + 1 2 t 0 h (X s )ds, t ∈ R + where h = -k -1 R * - + k + 1 R * + exists a.e.
Let T > 0 and f a function C(R + , R) such that supp (f ) ⊂ [0, T ]:

E x [f (τ )] = Ẽx f (τ ∧ T ) exp - T 0 k(X s )X s dX s - 1 2 T 0 (k(X s )) 2 X 2 s ds = Ẽx f (τ ) exp h(X T ) -h(X 0 ) - 1 2 T 0 h (X s ) + (k(X s )) 2 X 2 s ds = T 0 f (s) +∞ -∞ e h(ξ)-h(x) Ẽx [e -1 2 T 0 g(Xu)du |τ = s, X T = ξ] Px (τ ∈ ds, X T ∈ dξ) .
where g(x) = h (x) + k(x) 2 x 2 , x ∈ R * . Since (X t ) t≥0 is a BM under Px :

Px (τ ∈ ds, X T ∈ dξ) = Px (τ ∈ ds) Px (X T -s + m ∈ dξ) = |m -x| √ 2πs 3 e -(m-x) 2 2s ds 1 2π(T -s) e -(ξ-m-x) 2 2(T -s)
dξ and we easily conclude that the density of τ exists.

2. We use the link between p τ and p, namely:

p τ (t) = -1 2 ∂ x p(m, x, t)| x=m . Indeed, p τ (t) = -∂ t P(τ > t) = -∂ t m -∞ p(m, x, t)dx = - m -∞ ∂ x (k(x)xp(m, x, t)) + 1 2 ∂ 2 x p(m, x, t)dx = -k(x)xp(m, x, t) + 1 2 ∂ x p(m, x, t) m -∞ = - 1 2 ∂ x p(m, x, t)| x=m
We consider here only the case x 0 ∈ R -(for x 0 ∈]0, m] we use similar arguments). We use Formula (6), case 2, for computing the derivative ∂ ∂x p(m, x, t) at x = m. Since we are given an infinite sum, we have to prove that

∞ n=1 e -µnt N (µ n ) D µn k - (-x 0 2k -)∂ x α n D µn k + (-x 2k + ) + β n D µn k + (x 2k + ) (10)
converges uniformly in a neighborhood of m. Consider a neighborhood of the form [m -ε, m] with m -ε ≥ 0. Using the expression (8) of N (µ n ) and Corollary 6 in the appendix, we obtain for all

x ∈ [m -ε, m] ∂ x α n D µn k + (-x 2k + ) + β n D µn k + (x 2k + ) N (µ n ) ≤ 1 ∂ µ α µ D µ k + (-m 2k + ) + β µ D µ k + (m 2k + ) µ=µn . (11)
Thus, for proving that the sum ( 10) is uniformly convergent it is sufficient to show that

∞ n=1 e -µnt D µn k - (-x 0 2k -) ∂ µ α µ D µ k + (-m 2k + ) + β µ D µ k + (m 2k + ) µ=µn < ∞.
To complete the proof it suffices to apply (21) and ( 22) below. Thus for large n we have

e -µnt D µn k - (-x 0 2k -) ∂ µ α µ D µ k + (-m 2k + ) + β µ D µ k + (m 2k + ) µ=µn n→∞ ≤ Ce -k -(2n-1)t , (12)
with C constant.

Proof of Corollary 3

Notice that

P ( sup 0≤u≤t X u < m) = P (τ > t) = ∞ t p τ (s)ds,
and hence we obtain the formula of Corollary 3 by integrating the density (4) over [t, ∞).

Observe that we can interchange the integral and the infinite sum, because the series on RHS of Equation ( 4) is absolutely convergent. For the case x 0 ∈ R -this follows from [START_REF] Lejay | Maximum likelihood drift estimation for a threshold diffusion[END_REF] and the fact that the sequence (µ n ) n converges to infinity (see (21) below).

Appendix

Definitions and some properties on Hermite functions and parabolic cylinder functions

We recall here the definitions and some known results about parabolic cylinder functions and Hermite functions. These results can be found in Chapters 4 and 10 of [START_REF] Lebedev | Special functions and their applications[END_REF] and in Section 7 of [START_REF] Alili | Representations of the first hitting time density of an Ornstein-Uhlenbeck process[END_REF].

For ν ∈ R, the Hermite function H ν is defined, for all x ∈ R, by

H ν (x) = 1 2Γ(-ν) +∞ i=0 (-1) i i! Γ i -ν 2 (2x) i , (13)
or, equivalently, by

H ν (x) = 2 ν Γ 1 2 Γ 1-ν 2 1 F 1 -ν 2 , 1 2 , x 2 + x Γ -1 2 Γ -ν 2 1 F 1 1 -ν 2 , 3 2 , x 2 , ( 14
)
where 1 F 1 is the Kummer confluent hypergeometric function. One can show that

H ν (0) = 2 ν Γ 1 2 Γ 1-ν 2 , H ν (0) = 2 ν Γ -1 2 Γ -ν 2 . ( 15
)
The Hermite function satisfies the following relations:

∂ x H ν (x) = 2νH ν-1 (x), (16)
and

∂ 2 x H ν (x) = -2νH ν (x) + 2x∂ x H ν (x). (17)
For ν ∈ R we define the parabolic cylinder function as the function

D ν (x) = 2 -ν 2 e -x 2 4 H ν x √ 2 , x ∈ R. (18)
One can show that D ν is one of the two solutions of the differential equation

y (x) + - x 2 4 + 1 2 + ν y(x) = 0, x ∈ R. (19)
Notice that the second solution of ( 19) is given by x → D ν (-x). We have the following relation:

∂ x D ν (x) = νD ν-1 (x) - x 2 D ν (x). (20)
Finally, we have the following asymptotic expansion: 4 ) .

H ν (x) = 2 ν Γ 1+ν 2 e x 2 2 √ π cos x √ 2ν + 1 - πν 2 + O(ν - 1

Asymptotic expansions for large µ

For the proof of Theorem 1 we need some asymptotic properties of the function defined by

G(x, µ) = α µ D µ k + (-x 2k + ) + β µ D µ k + (x 2k + ).
We start by deriving the asymptotic behaviour of G(x, µ) as µ → ∞ (cf. the method given in [START_REF] Slater | Confluent hypergeometric functions[END_REF] p. 69 for the function φ).

Lemma 4. When µ → ∞, the following asymptotic formula holds: 4 ) .

G(x, µ) = 2 µ 2k - Γ 1 2 + µ 2k - √ π cos x k + + 2µ - πµ 2k - + O(µ - 1
where F (x, µ) = cos µπ 2k -cos x k + + 2µ +sin µπ 2k -sin x k + + 2µ

√ k -Γ 1+ µ 2k - 1 4 + µ 2k + -1 2 √ k + Γ 1 2 + µ 2k - .
Since for large n, µ n is solution of F (m, µ n ) = 0, we have

∂ µ G(m, µ n ) =∂ µ   2 µ 2k - Γ 1 2 + µ 2k - √ π 1 + O(µ -1 4 )   µ=µn F (m, µ n ) =0 + 2 µn 2k - Γ 1 2 + µn 2k - √ π 1 + O(µ -1 4 n ) ∂ µ F (m, µ) µ=µn =(-1) n π 2k - +O(µ -1 2
n ) =2 µn 2k - √ πΓ 1 2 + µn 2k - 2k - (-1) n 1 + O(µ -1 4 n ) .
Hence : Proof. We apply the method given in [START_REF] Alexandrov | An inequality of Duffin-Schaeffer type for Hermite polynomials[END_REF]. Since (-αH ν (-x) + βH ν (x)) = 2ν (αH ν-1 (-x) + βH ν-1 (x)) , we have f (x) = 1 4ν 2 [(-αH ν (-x) + βH ν (x)) ] 2 + 2ν [-αH ν (-x) + βH ν (x)] 2 . Using the fact that -αH ν (-x) + βH ν (x) is a solution of y (x) + 2νy(x) -2xy (x) = 0, we see that

D µn k - (-x 0 2k -) ∂ µ G(m, µ n ) = (-1) n e -x 2 0 k -2k - π cos x 0 k -+ 2µ n + πµ n 2k - + O(µ -1 4 n ). (22)

Auxiliary results

f (x) =
x (-αH ν (-x) + βH ν (x))

2 ν 2 .
Hence f is a decreasing function on R -and increasing on R + . We deduce that for all x ∈ [m -ε, m] and m -ε ≥ 0, [αH ν-1 (-x) + βH ν-1 (x)] 2 ≤ f (x) ≤ f (m), and the conclusion hold.

A direct consequence of this lemma is the following result :

Figure 1 :

 1 Figure 1: Comparison of expression (4) and a Monte Carlo method when m = 1, X 0 = -1 and k + = 2, k -= 1.

Figure 2 :

 2 Figure 2: Comparison of expression (6) and a Monte Carlo method when m = 2, X 0 = -1 and k + = 0.5, k -= 2.

e -k + x 2 -k + z 2 m x∨z e k + y 2 dy 2 dzdxe k + y 2 dy 2 dx

 222 <∞ on[0,m]

Corollary 4 .

 4 (a) The operator T : L

Lemma 5 .

 5 Let ν, m, ε ∈ R + such that m -ε ≥ 0. For all α, β ∈ R and x ∈ [m -ε, m],we have :| αH ν-1 (-x) + βH ν-1 (x) |≤ f (m), where f (x) = [αH ν-1 (-x) + βH ν-1 (x)] 2 + 1 2ν [-αH ν (-x) + βH ν (x)] 2 .

  )Y s ds is a BM under P y and (Y t ) t≥0 is an OU process with broken driftdY t = -k(Y t )Y t dt + dW t .

		2 Y 2 s ds ,
	and W t := Y t +	t 0 k(Y s

Sketch of the proof. Using ( 14), [START_REF] Slater | Confluent hypergeometric functions[END_REF], and the definition of α µ and β µ , we find:

Tricomi's formula (see for ex. [START_REF] Slater | Confluent hypergeometric functions[END_REF], p.57) gives

where J is the Bessel function and the coefficients A n are recursively defined by

Hence, we can write

where

From the equalities J -1 2 4 ).

We finally find 4 ) .

For large µ, the zeros of G(x, µ) are the zeros of U 0 which are very close of the zeros of cos x k + + 2µ -πµ 2k -. Hence, for large n, µ n is close of the positive solution m k + + 2µ -πµ 2k -= (-n + 1 2 )π and we have the asymptotics :

For µ in a neighborhood of a zero µ n , we can prove that