Souad Bezzaoucha

Benoît Marx

Didier Maquin

José Ragot

On the Unknown Input Observer Design : a Decoupling Class Approach with Application to Sensor Fault Diagnosis

This paper addresses fault diagnosis for observerbased residual generators for linear discrete-time systems subject to unknown input. The proposed approach is a new method allowing to characterize a class of unknown inputs from which the estimation error is decoupled. This contribution is divided into two parts. The first one concerns the design of the UIO satisfying either an exact decoupling or an L 2 -attenuation of the unknown input to the state estimation error. The second part is dedicated to the implementation of a bank of such observers for sensor fault detection and isolation.

I. INTRODUCTION

Due to an increasing demand for higher performances, safety and reliability, fault diagnosis for uncertain systems with Unknown Input (UI) has received considerable interest. Since in many cases a part of the system input is inaccessible (e.g. plant disturbance or actuator failure), a conventional observer that requires the knowledge of all inputs cannot be used directly ; then Unknown Input Observers (UIOs) were developed to estimate the state of uncertain systems despite the existence of UIs or disturbances [START_REF] Chang | Design of General Structured Observers for Linear Systems with Unknown Inputs[END_REF], [START_REF] Basile | On the observability of linear, time-invariant systems with unknown inputs[END_REF], [START_REF] Darouach | Full-order observers for linear systems with unknown inputs[END_REF], [START_REF] Guidorzi | On Wonham stabilizability condition in the synthesis of observers for unknown input systems[END_REF], [START_REF] Hostetter | Observing systems with unmeasurable inputs[END_REF], [START_REF] Hou | Design of observers for linear systems with unknown inputs[END_REF].

Classically, the state estimation of a system subject to UI can be obtained by means of the so-called UIO. The goal of the UIO is to provide state reconstruction of the system with some robustness with respect to possible UI. Design of UIO has been extensively investigated in the literature and is based either on the decoupling such that the estimation error do not depend on the UI [START_REF] Yang | Observers for linear systems with unknown inputs[END_REF] [START_REF] Hou | Design of observers for linear systems with unknown inputs[END_REF] [START_REF] Darouach | Full-order observers for linear systems with unknown inputs[END_REF], or on the synthesis of an Integral Observer for the estimation of disturbances [START_REF] Shumsky | Algebraic Approach to the Problem of Fault Accommodation in Nonlinear Systems[END_REF] [START_REF] Ichalal | State and unknown input estimation for nonlinear systems described by Takagi-Sugeno models with unmeasurable premise variables[END_REF] [START_REF] Ichalal | Simultaneous state and unknown inputs estimation with PI and PMI observers for Takagi Sugeno model with unmeasurable premise variables[END_REF]. These strategies frequently require structural and rank constraints on the system matrices.

In this paper, the proposed strategy consists to decompose any UI into two terms. The first one is a sum of exponential functions from which the state estimates can be exactly decoupled. For a given system, the class of the UI satisfying that property is clearly established. The effect of the remaining part of the UI on the state estimates is then attenuated in an L 2 framework.

Then, the proposed UIO design will be applied to investigate a sensor fault detection and isolation problem. A method based on the design of an observer bank will be used. This approach uses the proposed observer in a bank of observers. Residues corresponding to each observer are generated and defined in such a way to detect the fault occurrence ; coupled with a residual analysis methods, the faulty instrument (sensor in the considered case) is identified.

Different schemes of bank observers can be used for fault diagnosis (Dedicated Observer Scheme, Generalized Observer Scheme,..), see [START_REF] Chen | Robust detection of faulty actuators via unknown input observers[END_REF] [14] [START_REF] Ding | Model-Based Fault Diagnosis Techniques Design Schemes, Algorithms and Tools[END_REF]. In this paper, a Generalized Observer Scheme (GOS) is applied. The bank consists of N + 1 observers that include N fault observers and one normal condition observer. The i th observer is driven by all inputs and N -1 outputs of the system and generates the corresponding residual vector r i k . Then, these residual vectors are evaluated for fault detection and localisation.

This paper is organised as follows : Section II presents a second order system to introduce the decoupling strategy and the Uknown Input Class for exact decoupling notion and how to generate this class. Section III is a generalization of the second section. In section IV, we introduce the notion of partial decoupling and the linear matrix inequalities conditions to ensure the L 2 attenuation of the UI effect on state estimates. In order to improve the obtained results, a pole assignment is also implemented. However, the usual linearization approaches are not suitable to the present problem since BMIs (Bilinear Matrix Inequalities) are to be dealt with. A gain adjustment technique is then applied. This synthesis linearize the inequalities by fixing one of the unknown variable [START_REF] Losser | Ajustement de lois de commande : Application en aéronautique[END_REF]. This kind of procedure can be found in the centrage-XY procedure [START_REF] Iwasaki | The XY-centring algorithm for the dual LMI problem : a new approach to fixed-order control design[END_REF], the D-K iteration mentionned in [START_REF] Rotea | An alternative to the D-K iteration[END_REF] or Yamada's approach [START_REF] Yamada | An LMI approach to local optimization for constantly scaled H ∞ control problems[END_REF]. In section V, simulations are presented to show the efficiency of the proposed approach. And finally, in the last section, as an application, a residual generator design in a case of a sensor fault is addressed.

II. ILLUSTRATIVE EXAMPLE

To begin with, the procedure is introduced with the help of a simple example being a second order system. The different steps leading to the Uknown Input Class for exact decoupling are detailed. Consider a second order system described by :

x k+1 = Ax k + Bu k + Dη k-1 y k = Cx k + eη k-1 (1)
where x k ∈ R 2 , u k ∈ R, η k ∈ R and y k ∈ R are the system state, input, unknown input and the output vector respectively. The system matrices are real valued, constant and of appropriate dimensions :

x k = x 1k x 2k A = a 11 a 12 a 21 a 22 B = b 1 b 2 C = c 1 c 2 D = d 1 d 2 (2)
The proposed Proportional Integral Observer (PIO) of gain K and the UIO depending on an auxiliary variable z k ∈ R are respectivelly given by the following equations :

   xk+1 = A xk + Bu k + D ηk-1 + K ỹk ŷk = C xk + e ηk-1 ỹk = y k -ŷk (3)
z k+1 = γ 1 ỹk + λ 1 z k ηk+1 = γ 2 z k + λ 2 ηk (4)
with :

K = k 1 k 2
It can be noted that the UIO has a filter structure with as an input the output reconstruction error ỹk . The filter parameters γ 1 , γ 2 , λ 1 and λ 2 allow to modify the gain and time constants of the UIO. Depending on the value of the UIO parameters, we can have either a proportional observer, an integral or a multiple integral observer. In this particular case, the choice λ 1 = 1 or λ 2 = 1 introduces the two integrators in this filter structure.

In the following, the state and UI estimation errors are expressed as function of the UI. Since the system and its observer are linear, the time shift operator q (q f k = f k+1 ), is adequate to express the reconstruction errors :

xk = x k -xk ηk = η k -ηk (5)
From equations ( 1), ( 3) and ( 4), the state reconstruction error of the UI is given by :

x1k = N 1 (q) D(q) η k x2k = N 2 (q) D(q) η k (6)        N 1 (q) = (q -λ 1 )(q -λ 2 )( ã12 d2 -ã22 d1 + q d1 ) N 2 (q) = (q -λ 1 )(q -λ 2 )( ã21 d1 -ã11 d2 + q d2 ) D(q) = ((q -λ 1 )(q -λ 2 ) + γe)((q -ã11 )(q -ã22 ) -ã12 ã21 ) +γc 1 ( ã12 d2 -ã22 d1 + q d1 ) + γc 2 ( ã21 d1 -ã11 d2 + q d2 ) (7) with :                ã11 = a 11 -k 1 c 1 ã12 = a 12 -k 1 c 2 ã21 = a 21 -k 2 c 1 ã22 = a 22 -k 2 c 2 d1 = d 1 -k 1 e d2 = d 2 -k 2 e (8)
From ( 6), conditions for the estimation errors to be independent from the UI can easily be derived. Then, the UI familly satisfying an exact decoupling is solution of :

N 1 (q) D(q) η k = 0 N 2 (q) D(q) η k = 0 (9)
In order to find the solution η k assuring the previous conditions, it is imposed that polynomials N 1 (q) and N 2 (q) have the same roots. However, before that, it should also be checked if some solutions are common to D(q) and N 1 (q) (or N 2 (q)).

That leads to :

q 0 = ã21 d1 -ã11 d2 d2 (10)
which is a common root between N 1 (q), N 2 (q) and D(q). Thus condition ( 9) is reduced to :

(q -λ 1 )(q -λ 2 )η k = 0 ( 11
)
The solution is given by an UI being the sum of two exponential functions :

η k = A 1 λ k 1 + A 2 λ k 2 ( 12
)
where coefficients A 1 and A 2 are arbitrarily set. Finally, the choice of the observer values λ 1 and λ 2 (4) gives the UI class assuring the exact decoupling of the state error from η k for any values of the coefficients A 1 and A 2 .

III. RECONSTRUCTION ERRORS : DISTURBANCES

DECOUPLING

Let us now return to the general case by using the following system equations :

x k+1 = Ax k + Bu k + Dη k-1 y k = Cx k + eη k-1 (13)
Vectors

x k ∈ R n , u k ∈ R m , η k ∈ R
and y k ∈ R p are the system state, input, unknown input and the output vectors respectively. The system matrices

A ∈ R n×n , B ∈ R n×m , D ∈ R n×1 , C ∈ R p×n et e ∈ R p×1
are known real values.

The proposed system observer of gain K and the UIO are respectivelly given by the following equations :

   xk+1 = A xk + Bu k + D ηk-1 + K ỹk ŷk = C xk + e ηk-1 ỹk = y k -ŷk (14)
z k+1 = Γ ỹk + Λz k ηk+1 = γz k + λ ηk (15)
with appropriate dimensions :

z k ∈ R q , K ∈ R n×p , Γ ∈ R q×p , γ ∈ R 1×q , Λ ∈ R q×q and λ ∈ R.
By following the same steps as in the previous section, the state and UI reconstruction errors are expressed ; we get from equation [START_REF] Iwasaki | The XY-centring algorithm for the dual LMI problem : a new approach to fixed-order control design[END_REF] with the time operator q : ηk = (qλ ) -1 γz k (16)

(qI q -Λ)z k = Γ ỹk ( 17
)
which leads to :

[(qI q -Λ) + Γe(q -λ ) -1 q -1 γ]z k = ΓC xk + Γeq -1 η k (18)
The state error dynamics is obtained from ( 13) and ( 14) :

xk+1 = A xk + D ηk-1 A = A -KC D = D -Ke (19)
That gives the state estimation error :

xk = (qI n -A) -1 Dq -1 η k -(qI n -A) -1 Dq -1 (q -λ ) -1 γz k (20)
By replacing this expression in [START_REF] Peaucelle | A new robust D-stability condition for real convex polytopic uncertainty[END_REF], we have :

z k = Z -1 Λη k (21)
with :

Λ = ΓC(qI n -A) -1 D + Γe Z = q(qI q -Λ) + Λ(q -λ ) -1 γ (22)
Finally, replacing ( 21) in ( 16) and ( 20) leads to :

ηk = (q -λ ) -1 γZ -1 Λη k xk = (qI n -A) -1 Dq -1 1 -(q -λ ) -1 γZ -1 Λ η k (23)
The UI estimation error becomes :

ηk = 1 -(q -λ ) -1 γZ -1 Λ η k
From ( 23) the state estimation error decoupling condition from the UI can be written as :

(qI n -A) -1 Dq -1 1 -(q -λ ) -1 γZ -1 Λ η k = 0 (24)
In order to decouple the state from the UI and assure its exact estimation, the following condition has to be verified :

1

-(q -λ ) -1 γZ -1 Λ η k = 0 (25)
Equation ( 25) may be extended as N(q) D(q) η k = 0. Solving this last equation gives roots defining the UI class that ensure an exact decoupling of the estimation error from the UI. This class is written as : ∑ i A i λ i where the λ i correspond to the roots of (25) and A i are totally free parameters.

IV. PARTIAL DECOUPLING OBSERVER

In the previous section, was detailed how to find the class of UI ensuring an exact decoupling of the UI in respect to the state estimation error. In the following section, a general case with an UI that does not satisfy the decoupling condition is considered. In this case, the problem is solved by attenuating the effect (transfer) of the UI to the estimation error and propose linear matrix inequalities to determinate the observer gain so that the estimated state asymptotically tends to the real one. In addition to the two previous cases (exact and partial decoupling), we also have a third one, which is a mix between the two solutions. In fact, any UI may be decomposed into a sum of two terms η k = η d k + η a k . The first term corresponds to the exact decoupling term obtained as explain in section III, and the second one is the approximaion term onto L 2 attenuation is applied. In subsection A, we only present the attenuation approach ; but, in the simulation section the combined approach will be illustrated.

A. L 2 Attenuation

System and observer equations are given by :

   xk+1 = A xk + D ηk-1 ηk = η k -λ η k-1 -γz k-1 + λ ηk-1 z k+1 = ΓC xk + Γe ηk-1 + Λz k (26)
The corresponding matrix form is given by :

e k+1 = A 1 e k + B 1 η a k (27)
with :

A 1 =     A D 0 0 0 λ 0 -γ ΓC Γe Λ 0 0 0 1 0     B 1 =     0 0 1 -λ 0 0 0 0     , e k =     xk ηk-1 z k z k-1     η a k = η k η k-1 (28)
In particular, (27) gives the UI influence on the estimation errors. To focus on the impact of the UI on the state estimation xk , a new observer output is considered :

g k = C 1 e k (29)
with : C 1 = I 0 0 0 . Considering the Real Bounded Lemma [START_REF] Boyd | Matrix Inequalities in System and Control Theory[END_REF], the system ( 27) is stable and the L 2 gain from η a k to g k is bounded by

g k 2 η a k 2
< µ if there exists a positive symmetric matrice P and a positive scalar µ such that the following condition holds :

  A T 1 PA 1 -P A T 1 PB 1 C T 1 B T 1 PA 1 B T 1 PB 1 -µ 2 I 0 C 1 0 -µ 2 I   < 0 (30)
According to [START_REF] De Oliveira | A new discrete-time robust stability condition[END_REF] and [START_REF] Peaucelle | A new robust D-stability condition for real convex polytopic uncertainty[END_REF], the previous problem can be reformulated by searching a positive symmetric definite matrice P, gains K and G such that :

    -P A T 1 PB 1 C T 1 A T 1 G T B T 1 PA 1 B T 1 PB 1 -µ 2 I 0 0 C 1 0 -µ 2 I 0 GA 1 0 0 -G -G T + P     < 0
(31) where A 1 defined in (28) with the help of (19), depends on K. Due to this dependence, let us remark that inequality (31) is not linear. For that reason some transformations are needed to obtain LMIs. Let us write the matrix A 1 such that :

A 1 = A 1 -RKB 1 ( 32
)
with :

A 1 =     A D 0 0 0 λ 0 -γ ΓC Γe Λ 0 0 0 1 0     R =     I 0 0 0     B 1 = C e 0 0 (33)
Replacing A 1 by (32) in (31), we have :

    -P A T 1 PB 1 C T 1 A T 1 G T B T 1 PA 1 B T 1 PB 1 -µI 0 0 C 1 0 -µI 0 GA 1 0 0 P -G -G T     +M T N +N T M < 0 (34) with M =     -B T 1 K T 0 0 0     T , N =     0 B T 1 PR 0 GR     T and µ = µ 2
Let us recall the following lemma [START_REF] Zhou | Robust stabilization of linear systems with norm-bounded time-varying uncertainty[END_REF]. Consider two real matrices Π and Λ with appropriate dimensions, for any positive matrix Σ the following inequality holds :

Π T Λ + Λ T Π ≤ Π T ΣΠ + Λ T Σ -1 Λ (35)
Applying this lemma, (34) becomes :

    -P A T 1 PB 1 C T 1 A T 1 G T B T 1 PA 1 B T 1 PB 1 -µI 0 0 C 1 0 -µI 0 GA 1 0 0 P -G -G T     + M T ΣM + N T Σ -1 N < 0 (36)
Applying Schur's complement, we get :

        -P A T 1 PB 1 C T 1 A T 1 G T B T 1 K T 0 B T 1 P B T 1 PB 1 -µI 0 0 0 B T 1 PR C 1 0 -µ 0 0 0 GA 1 0 0 -G -G T + P 0 GR KB 1 0 0 0 -Σ -1 0 0 R T PB 1 0 R T G T 0 -Σ         <0
(37) At last, by congruence, (37) becomes :

        -P A T 1 PB 1 C T 1 A T 1 G T B T 1 F T 0 B T 1 P B T 1 PB 1 -µI 0 0 0 B T 1 PR C 1 0 -µ 0 0 0 GA 1 0 0 -G -G T + P 0 GR FB 1 0 0 0 -Σ T 0 0 R T PB 1 0 R T G T 0 -Σ         < 0
(38) with F = ΣK. The LMI must be solved in respect to P, G, F and the gain K is obtained by K = Σ -1 F.

B. Pole Assignment

The minimization of the attenuation factor µ may result in slow dynamics of the state estimation error. This problem can be solved by pole assignment of the closed loop system in a specified region. The considered region is a disk centred at (q, 0) with radius α. Thus, the condition to answer this constraint is given by the following : find P = P T > 0 and Q = Q T > 0 such that the following LMI [START_REF] Chilali | H ∞ design with pole placement constraints : an LMI approach[END_REF] holds :

-αQ -qQ + QA 1 -GC (-qQ + QA 1 -GC) T -αP < 0 (39)
with G = QK. We have to solve this LMI regarding to Q and G then we deduce K. Thus, to ensure the stabilty and pole assignment, the conditions (38) and (39) must be fulfilled simultaneously.

C. Gain Ajustement

From matrices F and G definitions, there is a dependence between the two LMIs (38) and (39). Then we have to solve simultaneously these two LMIs which can be noted LMI1(P, K) and LMI2(Q, K). The proposed method is based on an ajustment technique allowing to set some variables and calculate others in an iterative way. More precisely, if the gain K is fixed, we solve LMI1(P, K) regarding to P. Then we solve LMI2(Q, K) regarding to Q and K and use the obtained result K for the next iteration (see table 1). This procedure was chosen in reason of its simplicity, but one should be aware that no optimality or convergence guarantee is given. However, since our study goal is to find a solution to the given conditions, an optimal solution is not a necessity.

Iterative optimisation for gain K :

1) Set i = 0. Choose a stabilisable value K 0 . Put

K (i) = K 0 . 2) L 2 attenuation : Find P (i+1) > 0 solution of LMI1(P, K (i) ). 3) Pole assignment : Find Q (i+1) and K (i+1) solution of LMI2(Q, K) . 4) Stopping condition : -If ||K (i+1) -K (i) || < ε stop the algorithm : K f inal = K (i+1) . -Else, set i = i +
Λ = 0.33   1 1 1 0 1 1 0 1 1   λ = 0.7 Γ = 0.2I 3 γ = -0.4 1 1 1
At the first step, let us determine the observer gain K with the proposed iterative algorithm. The obtained gain K and attenuation µ for a pole assignment in a disk centred at (0.3,0) with radius 0.2 are : The second step consists of finding the UI class for an exact decoupling. Let us recall that the state decoupling condition from the UI with its exact estimation of the UI is given by (25). In this example, it corresponds to an UI composed of a linear combination of seven exponential functions : two roots of (25) are complex conjugate and the others have real values given by : λ 1 = 1 ; λ 2 = 0.7 ; λ 3 = 0.83 ; λ 4 = 0.15 ; λ 5 = 0.5 λ 6,7 = 0.085 ± 0.23i

K =       0.
Then, the class of UI for an exact decoupling is given by :

η d k = A 1 λ k 1 +A 2 λ k 2 +A 3 λ k 3 +A 4 λ k 4 +A 5 λ k 5 +A 6 a k cos(φ k +ψ) (

VI. RESIDUAL GENERATOR DESIGN

The residual generator design is addressed in this section. Based on the system structure, a bank of observers is then designed using the developed UIO in order to detect and isolate a sensor fault through the estimation of system outputs using measurable signals and the model of the system. The procedure is performed by analysing the time-evolution of the residual signals obtained by the comparison between the measured outputs and the estimated outputs [17] [7]. A GOS structure for the observer bank is adopted (Fig. 4). In theory, the residual signals (i.e. the output estimation error) are null under normal operating conditions of the system. The residual signal structuring, in order to generate appropriated fault indicators, can be obtained by replacing the use of only one observer by the use of a bank of observers where each observer is driven by a partial set of the available signals.

Let us consider the case of a sensor fault occuring at the first sensor of magnitude 1 in the time interval t ∈ [START_REF] De Oliveira | A new discrete-time robust stability condition[END_REF][START_REF] Iwasaki | The XY-centring algorithm for the dual LMI problem : a new approach to fixed-order control design[END_REF]. For the simulation, a normally distributed noise of standard deviation equal to 0.1 is added to all the outputs whose magnitudes are varying between -5 and 8. The first row of Fig. 5 represents the residuals under normal operating conditions with measurement noise. A fault signature localization method had been considered. By comparing the theoretical study (truth table) [START_REF] Patton | Issues of Fault Diagnosis for Dynamic systems[END_REF] [7] and the obtained residues (simulations), the faults susceptible to Summurising this section, an observer bank for residual generator was considered. From the obtained results we showed that the previously developped UIO can be used for the detection and isolation of faults when the system is subject to UIs, fault sensors and measurement noise.

VII. CONCLUSION AND PERSPECTIVES

This paper addresses new method to design observers with unknown inputs. The proposed approach is based on a partial decoupling of the state estimation from the UI without any rank constraints on the system matrix. The main result is about the way to find the UI class ensuring an exact decoupling. The proposed work can be extended to the nonlinear case, in particular, systems with Takagi-Sugeno representation.

  1 and go back to step 2.

2 √

 2 41)with : a = Re(λ 6 ) 2 + Im(λ 6 ) 2 and cos(φ) = Re(λ 6 ) 2 Re(λ 6 ) 2 +Im(λ 6 ) 2The UI is defined by :η k = 0.1 -0.1(0.7) k -0.5(0.83) k + 0.3(0.15) k -0.4(0.5) k + 0.4(0.47) k cos(1.78k) (42) Finally, the considered UI η k can be written as η k = η d k + η a k where η d k corresponds to the UI for exact decoupling and η a k to the approximation error. The following figures are obtained for the initial conditions x 0 = 0.5 0.1 0.2 -0.1 0 T and x0 = -0.5 0.5 -0.4 0.2 0.2 T . Fig.1 shows the system inputs. Fig.2 represents the UIs (for the exact η k = η d k and partial decoupling cases η k = η d k + η a k ) and their estimates and Fig 3. represents the system state and their estimate for both situations of exact and partial decoupling. In both situations, the state estimation is satisfactory.

Fig. 1 .

 1 Fig. 1. System inputs

Fig. 2 .Fig. 3 .

 23 Fig. 2. UI and its estimate

Fig. 5 .

 5 Fig. 4. GOS Structure


Jean-Jacques Sinou

Béatrice Faverjon

The vibration signature of chordal cracks in a rotor system including uncertainties

The aim of this paper is to investigate the effects of the presence of a transverse crack in a rotating shaft under uncertain physical parameters in order to obtain some indications that might be useful in detecting the presence of a crack in rotating system. The random dynamic response of the cracked rotor is evaluated by expanding the changing stiffness of the crack (i.e. the breathing mechanism) as a random truncated Fourier series. To avoid the use of the Monte Carlo simulations (MCS), an alternative procedure that is based on a combination of the Harmonic Balance Method and the Stochastic Finite Element Method (SFEM) using the Polynomial Chaos Expansion (PCE) is proposed. So the response of the Fourier components of the cracked rotor is expanded in the polynomial chaoses. The random dynamic response obtained by applying this procedure is compared with that evaluated through numerical integration based on the Harmonic Balance Method and the Monte Carlo simulations.

Introduction

The effect of crack in rotors is important to ensure the integrity of structural components for a wide range of engineering applications. In recent years a lot of effort has been devoted to the detection of transverse cracks in shafts [START_REF] Gasch | Dynamic behaviour of a simple rotor with a cross-sectional crack[END_REF][START_REF] Wauer | Dynamics of cracked rotors: Literature survey[END_REF][START_REF] Dimarogonas | Vibration of cracked structures: a state of the art review[END_REF][START_REF] Davies | The vibrational behaviour of a multi-shaft, multibearing system in the presence of a propagating transverse crack[END_REF][START_REF] Mayes | Analysis of the response of a multi-rotor-bearing system containing a transverse crack in a rotor[END_REF]. In most of the studies for crack detection in mechanical or rotor systems, researchers used changes in natural frequencies [START_REF] Salawu | Detection of structural damage through changes in frequency: a review[END_REF] and evolution of the non-linear behaviour of the system at the super-harmonics components as the diagnostic tools. For example, Gash [START_REF] Gasch | Dynamic behaviour of a simple rotor with a cross-sectional crack[END_REF][START_REF] Gasch | A survey of the dynamic behaviour of a simple rotating shaft with a transverse crack[END_REF] demonstrated that a slight decrease in the natural frequencies and the appearance of new resonances when the rotational speeds of the rotor reach 1 2 , 1 3 and 1 4 of the resonant frequencies of the rotor are key indicators for the detection of transverse cracks in a rotating shaft. He also suggested that the non-linear behavior of the damage rotor and the opening and closing of the crack during its rotation are due mainly to the shaft's self-weight. Then, Sinou [START_REF] Sinou | Detection of cracks in rotor based on the 2x and 3x super-harmonic frequency components and the crack-unbalance interactions[END_REF] indicated that the vibration amplitudes in the 1 n sub-critical resonances (with n ≥ 2) depend not only on the rotor damping, unbalance, position and depth of the crack, but also on the combinations of the unbalance and the crack parameters. The sensitivity of the magnitudes of 1 2 and 1 3 sub-critical resonances with respect to the unbalance angle 1 and the unbalance-crack interactions are different in the vertical and horizontal directions. However, the author suggested that the magnitudes of the sub-critical resonance peaks do not greatly change if the crack effect is predominant. Adewusi and Al-Bedoor [START_REF] Adewusi | Wavelet analysis of vibration signals of an overhang rotor with a propagating transverse crack[END_REF][START_REF] Adewusi | Experimental study on the vibration of an overhung rotor with a propagating transverse crack[END_REF] investigated an experimental study on the non-linear dynamic response of an overhung rotor with a propagating transverse crack. They observed that the appearance of second harmonic component is the first signature of the presence of damage. It was demonstrated that the third harmonic components are excited just before fracture, and so can not be used as an efficient and robust factor for the detection of damage in early stages. Some researchers [START_REF] Adewusi | Wavelet analysis of vibration signals of an overhang rotor with a propagating transverse crack[END_REF][START_REF] Sekhar | Crack detection and monitoring in a rotor supported on fluid film bearings: start-up vs run-down[END_REF][START_REF] Sekhar | Crack detection through wavelet transform for a run-up rotor[END_REF][START_REF] Zou | A comparative study on time-frequency feature of cracked rotor by wigner-ville distribution and wavelet transform[END_REF][START_REF] Sinou | An experimental investigation of condition monitoring for notched rotors through transient signals and wavelet transform[END_REF] proposed to introduce time-scale signal processing tools based on wavelet transform for the detection of damage in rotating machinery due to the limitation of the conventional Fourier analysis that is suitable for steady state vibration signals, but provides a poor representation of transient signals due to a strong increase of the shaft's rotating speed.

Due to the fact that crack breathes when the system rotates, the cracked rotor has a time-dependent coefficient. The amount of open part of the crack varies with the rotation of the shaft, thereby changing the stiffness of the cracked rotor. Gash [START_REF] Gasch | Dynamic behaviour of a simple rotor with a cross-sectional crack[END_REF][START_REF] Gasch | A survey of the dynamic behaviour of a simple rotating shaft with a transverse crack[END_REF] also suggested that the opening and closing of the crack during its rotation are due mainly to the shaft's self-weight. Researchers highlight the possibility of crack detection through the observation of the non-linear dynamic behaviour of rotor systems due to the this breathing mechanism [START_REF] Gómez-Mancilla | The influence of crack-imbalance orientation and orbital evolution for an extended cracked jeffcott rotor[END_REF][START_REF] Gómez-Mancilla | The influence of crack-imbalance orientation and orbital evolution for an extended cracked jeffcott rotor[END_REF][START_REF] Bachschmid | Some remarks on breathing mechanism, on non-linear effects and on slant and helicoidal cracks[END_REF][START_REF] Sinou | Influence of cracks in rotating shafts[END_REF][START_REF] Sinou | A non-linear study of a cracked rotor[END_REF][START_REF] Machorro-López | Identification of damaged shafts using active sensing-simulation and experimentation[END_REF][START_REF] Sinou | Experimental study on the nonlinear vibrations and nx amplitudes of a rotor with a transverse crack[END_REF]. They concluded that the primary response characteristic resulting from the changes in the non-linear dynamical behaviour of the rotor system through half resonance speeds appears to be the characteristic signature for detecting the presence of a crack rotor. They also indicated that emerging of the 3X super-harmonic frequency components and the distortion of the orbit, and formation of a double loop and inside loop in the orbit [START_REF] Sinou | Detection of cracks in rotor based on the 2x and 3x super-harmonic frequency components and the crack-unbalance interactions[END_REF], may provide useful information on the presence of a crack and may be used on an on-line crack monitoring rotor system. For example, Darpe [START_REF] Darpe | Transient response and breathing behaviour of a cracked jeffcott rotor[END_REF] studied both numerically and experimentally the evolution of orbital patterns during the passage through sub-critical resonances. The authors explained that the inner loop is present at one-half of the first resonant speed due to dominant second harmonic component. When the rotor is passing through the 1 2 sub-critical speed, the inner loop changes its orientation by almost 180 degrees. Moreover, Darpe et al. [START_REF] Darpe | Analysis of the response of a cracked jeffcott rotor to axial excitation[END_REF] used impulse axial excitation to a rotating cracked shaft to detect damage. They indicated that the presence of damage induces not only the coupling mechanism in lateral and longitudinal vibrations but also the combination of harmonics due to the interaction of rotational frequency and its harmonics with the constant excitation frequency and its harmonics. Then the same results of coupling measurements and appearance of combination harmonics have been observed for slant crack in rotor: coupled longitudinal, bending and torsional vibrations [START_REF] Darpe | Coupled bending, longitudinal and torsional vibrations of a cracked rotor[END_REF][START_REF] Darpe | Coupled vibrations of a rotor with slant crack[END_REF].

In practice the material parameters and the unbalance excitation force in the rotor systems are subjected to inherent physical variability and so are not exactly known. The dynamic properties of the rotor system can be deeply affected by small variations of some of these parameters and so they can be within a band of possible values due to uncertainties of some parameters such as stiffness and unbalance. However, the efficiency of the most practical indicators for health monitoring purposes is not analyzed in the presence of uncertainties present on rotating systems. The paper addresses this important issue of the changes in response of a cracked rotor under uncertain parameters and the robustness of indicators for monitoring rotating machinery.

Numerous methods have been developed to quantify physical uncertainties in a variety of computational problems e.g. the Monte Carlo approach, the perturbation method and the polynomial chaos expansions. The Monte Carlo simulations consists in including the uncertainties by calculating n times the deterministic problem, that is, for n independent samples of the random parameter. Therefore, the statistical response is easily obtained. With its easy use, the Monte Carlo method is one of the most frequently used. However, the convergence rate of this method is very slow which brings computations to be intensive really fast and not adapted to high-dimensional deterministic problems or high CPU time deterministic computations. In particular, rotordynamic problems are quite complex in a deterministic sense. The second type of methods includes either the perturbation method [START_REF] Loève | Perturbation Methods[END_REF] random quantities in expansion of Taylor series or the Neumann method based on Neumann series expansion [START_REF] Benaroya | Finite element methods in probabilistic structural analysis: A selective review[END_REF][START_REF] Yamazaki | Neumann expansion for stochastic finite element analysis[END_REF]. These methods obtain acceptable results only if the random fluctuations are small, therefore if we consider a problem with one excitation close to the resonance frequency, they can not be efficiently used. Finally, the polynomial chaos expansion associated with a Galerkin projection, named the Stochastic Finite Element Method [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] showed its effectiveness for quantifying uncertainties through the value of stochastic processes or variables in orthogonal bases of random variables. It is based on the homogeneous chaos theory of Wiener [START_REF] Wiener | The homogeneous chaos[END_REF] and the original polynomial chaos expansion [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] that consists in a mean-square convergent expansion of multidimensional Hermite polynomials of normalized Gaussian variables. The Hermite polynomials are orthogonal with respect to the Gaussian measure then optimal exponential convergence is obtained for Gaussian inputs [START_REF] Ghanem | Stochastic finite elements with multiple random non-gaussian properties[END_REF]. Besides, many extensions were made on the nature of the expansion (e.g. non Gaussian variables) to solve very different types of problems, especially in structural dynamics [START_REF] Soize | Physical systems with random uncertainties: Chaos representations with arbitrary probability measure[END_REF].

In this paper, due to the fact that determining all the Fourier components of the rotor system for the deterministic problem is already computationally intensive, the computational costs of Monte Carlo simulations can become impractical. Consequently, we use the polynomial chaos expansion method in the form of the stochastic finite element model that seems to be well adapted for a combination with the Harmonic Balance Method: the response of the Fourier components of the rotor system returns as random and is then expanded in the polynomial chaos basis. To validate the method, all the calculations are compared to those obtained from the combination of the Harmonic Balance Method and the Monte Carlo simulations.

First of all, a brief introduction of the rotor system and the modeling of the crack based on the notion of stiffness reduction are described. Then, the paper presents the numerical procedure for calculating the response of a cracked rotor under uncertain parameters. To quantify the effects of uncertainties on response variability in cracked rotor systems the stochastic response via the combination of the Harmonic Balance Method and the Polynomial Chaos Expansion procedure are explained and studied. Finally, the results obtained by applying the alternative approach via the Harmonic Balance Method and the Polynomial Chaos Expansion are compared with those evaluated by the Harmonic Balance Method with the Monte Carlo Simulations. The efficiency and the robustness of the combination of the Harmonic Balance Method with the Polynomial Chaos Expansion is demonstrated through different numerical simulations in order to analyze the random non-linear response of a cracked rotor against uncertain parameters and random excitation to assess its accuracy and the computation time.

The model of the cracked rotor

The layout of the cracked rotor system under consideration is shown in Figure 1. The rotor is composed of a shaft of length 0.5m with one disc at the mid-span. All the values of the material properties and dimensions of the rotor are given in Table 1. The uncraked rotor is discretized into Timoshenko beam finite elements, with four degrees of freedoms at each node, the axial and torsional degrees of freedom being not considered. After assembling the various shaft elements, the rigid disc and the bearing supports, the equations of the uncracked rotor can be written as

Mẍ + D ẋ + Kx = f + q (1)
where f and q contain the gravitational and balance forces, respectively. The mass matrix M includes the mass matrices of the shaft and the rigid disc. The matrix D combines the effects of the shaft's internal damping, damping of the supports, and gyroscopic moments. The matrix K includes the stiffness matrices of the shaft and the supports, together with the circulatory matrix which accounts for the shaft's internal damping. Table 1: Value of the physical parameters Then, due to strain energy concentration in the vicinity of the tip of the crack under load, the presence of a transverse crack introduces local flexibility. For a comprehensive literature survey of various crack modeling techniques, see [START_REF] Wauer | Dynamics of cracked rotors: Literature survey[END_REF][START_REF] Dimarogonas | Vibration of cracked structures: a state of the art review[END_REF]. In this paper, the model proposed by Mayes and Davies [START_REF] Davies | The vibrational behaviour of a multi-shaft, multibearing system in the presence of a propagating transverse crack[END_REF][START_REF] Mayes | Analysis of the response of a multi-rotor-bearing system containing a transverse crack in a rotor[END_REF] is used in order to locally represent the stiffness properties of the crack cross section. This model considers the reduction of the second moment of area ∆I of the element at the location of the crack that may be defined by

∆I = I 0    R l 1 -ν 2 F (µ) 1 + R l 1 -ν 2 F (µ)    , (2)
where I 0 , R , l, and ν are the second moments of area, beam radius, length of the section and Poisson's ratio, respectively. µ is the non-dimensional crack depth and is given by µ = h R where h defines the crack depth of the beam, as shown in Figure 1. F (µ) defines the non-linear compliance as a function of variations in non-dimensional crack depth µ, which can be derived from a series of experiments using chordal cracks (see Mayes and Davies [START_REF] Davies | The vibrational behaviour of a multi-shaft, multibearing system in the presence of a propagating transverse crack[END_REF][START_REF] Mayes | Analysis of the response of a multi-rotor-bearing system containing a transverse crack in a rotor[END_REF]). For the interested reader, the complete expressions of the stiffness matrix K crack at the transverse crack location are given in [START_REF] Sinou | Influence of cracks in rotating shafts[END_REF]. In this study, the crack model is restricted to a chordal crack. Then, one of the models of a crack is that of Mayes and Davies [START_REF] Mayes | Analysis of the response of a multi-rotor-bearing system containing a transverse crack in a rotor[END_REF] where the opening and closing of the crack were described by a cosine function by assuming that the gravity force is much greater than the imbalance force: when a cracked rotor rotates slowly under the load of its own weight, the crack will open and close once per revolution. This periodic opening and closing of the crack is called "breathing" phenomenon [START_REF] Gasch | Dynamic behaviour of a simple rotor with a cross-sectional crack[END_REF]. The function describing a simple crack breathing phenomenon such as discussed by Sinou and Lees [START_REF] Sinou | Influence of cracks in rotating shafts[END_REF][START_REF] Sinou | Effects of a crack on the stability of a non-linear rotor system[END_REF], typical in weight dominated systems, may be approximated by a cosine function g(t)

g (t) = 1 -cosωt 2 (3)
where ω defines the rotational speed of the rotor. If g(t) = 0, the crack is closed and has no effect on the dynamic behaviour of the rotor (i.e. the rotor may be treated as uncracked). If g(t) = 1, the crack is fully open. Finally, the equations of the cracked rotor can be written as

Mẍ + D ẋ + (K -g (t) K c ) x = f + q (4)
where ẍ, ẋ and x are the acceleration, velocity and displacement vectors. M and K are the mass and stiffness matrices of the complete uncracked rotor. f and q are the gravitational and balance forces, respectively. The global stiffness matrix K c due to the presence of the crack is given by

diag (K c ) = ( 0 • • • 0 K crack 0 • • • 0 ) ↑ i th element (5)
where K crack defines the stiffness matrix of the crack element. It is situated at the i th beam location. 0 defines the 8 × 8 null matrix.

3 Response of the cracked rotor under uncertain parameters via the Harmonic Balance Method and the Polynomial Chaos Expansion (HBM-PCE)

Deterministic response via the Harmonic balance Method

As previously indicated, the above equations of the cracked rotor have a time-dependent coefficient due to the fact that the crack breathes when the system rotates. The amount of open part of the crack constantly varies with the rotation of the shaft, thereby changing the stiffness of the cracked rotor. The global stiffness matrix of the rotor consists of a constant part K and a time dependent part g (t) K c . In order to determine the periodic solutions of the cracked rotor system, we assume that the dynamical responses of the cracked rotor system can be approximated by a truncated Fourier series of order m with a fundamental frequency

x (t) = A 0 + m k=1 (A k cos (kωt) + B k sin (kωt)) (6)
where ω defines the rotational speed of the cracked rotor system. A 0 , A k and B k (with k = 1, • • • , m) define the unknown coefficients of the finite Fourier series. The number of harmonic coefficients m is selected on the basis of the number of significant harmonics expected in the dynamical response of the crack rotor system. Then, it may be observed that the gravitational force f and the unbalance force q can also be exactly defined by finite Fourier series (with constant components and first-order periodic components in the frequency domain, respectively). Expressions are given by

f = C f 0 (7) q = C q 1 cos (ωt) + S q 1 sin (ωt) (8)
where C f 0 is the vector of the constant Fourier components of the gravitational force. C q 1 and S q 1 are first-order periodic components in the frequency domain of the unbalance force. By combining the approximated dynamical response of the rotor x (see equation 6) and the the global mechanism of the time dependent part (see equation 3), the expression g (t) K c x can be approximated by considering truncated Fourier series of order m

g(t)K c x = m k=0 K c - 1 2 A k-1 + 1 2 A k - 1 4 A k+1 cos(kωt)+ m k=1 K c - 1 2 B k-1 + 1 2 B k - 1 4 B k+1 sin(kωt) (9)
with B 0 = A -1 = 0 and k + 1 < m. Finally, substituting Equations 6, 7, 8, 9 into Equations 4 yields a set of (2m + 1) × n linear equations (where n is the number of degrees of freedom for the complete cracked rotor)

(Λ + Λ c ) Θ = Γ (10) with Θ = [A 0 A 1 B 1 A 2 B 2 • • • A k B k • • • A m B m ] T (11)
Γ = C f 0 C q 1 S q 1 0 0 • • • 0 0 • • • 0 0 T (12) Λ =             K Λ 1 Λ 2 . . . Λ k . . . Λ m             with Λ k = K -k 2 ω 2 M kωD -kωD K -k 2 ω 2 M ( 13
)
Λ c = 1 4                     2K c -K c -2K c 2K c 0 -K c 0 0 0 2K c 0 -K c -K c 0 2K c 0 -K c 0 0 -K c 0 2K c 0 -K c . . . -K c 0 2K c 0 -K c 0 0 -K c 0 2K c -0 -K c . . . -K c 0 2K c 0 0 -K c 0 2K c                     (14)
The unknown Fourier coefficients Θ can be found solving Equation 10.

Stochastic response via the Harmonic Balance Method-Polynomial Chaos Expansion (HBM-PCE)

Uncertainties on the stiffness parameters and on the unbalance excitation force have to be modeled. Λ, Θ and Γ u are random processes with argument τ denoting the random character, Λ c is considered as deterministic. Equation 10 can then be rewritten in a random way such that

(Λ(τ ) + Λ c ) Θ(τ ) = Γ(τ ) (15)
In fact, excitation force Γ consists of a deterministic part that is the gravity component Γ g and of a random part that is the unbalance component Γ u (τ ) such as

Γ = Γ g + Γ u (τ ) (16)
From Equation 12, it yields :

Γ g = C f 0 0 • • • 0 T and Γ u (τ ) = [ 0 C q 1 (τ ) S q 1 (τ ) 0 • • • 0] T (17)
Λ, defined by its deterministic expression by Equation 13, can be randomized using the Karhunen Loeve expansion with a Gaussian law on the stiffness of the rotor. The use of Karhunen-Loeve expansion allows to work with random variables instead of random processes [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Nayfeh | Probability Theory[END_REF]. It should be noted that the support boundary conditions are taken as deterministic so the stiffness associated and denoted K b is not taken as random. Then, considering that K = K b + K, Equation 13 becomes

Λ = Λ d + Λ r (τ ) (18)
where Λ d is the deterministic part of Λ defined by

Λ d =      K b Λ d 1 . . . Λ d m      with Λ d k = K b -k 2 ω 2 M kωD -kωD K b -k 2 ω 2 M ( 19
)
and Λ r (τ ) is the random one given by

Λ r (τ ) =    K(τ ) . . . K(τ )    (20)
K(τ ) is expanded using the Karhunen-Loeve expression implemented in the Galerkin formulation of the finite element method [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] such as

K(τ ) = L i=1 K i (1 + δ K i ξ K i (τ )) , (21)
in which {ξ i }, i = 1 to L, is a set of orthonormal random variables, Z is the mean of random quantity Z and δ Z is its variation coefficient.

Considering the unbalance excitation force Γ u (τ ), several random ways of modeling can be possible. One simple one but sufficiently representative of the reality is done through the Polynomial Chaos Expansion with a Gaussian law. The associated basis is an orthogonal basis set of random variables, represented in a mean-square convergent expansion in terms of multidimensional Hermite polynomials of normalized Gaussian variables [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. We have

Γ u (τ ) = ∞ j=0 F u j Ψ j (ξ(τ )) (22)
The unknown solution Θ(τ ) is also then expanded in chaos polynomials such as

Θ(τ ) = ∞ j=0 Y j Ψ j (ξ(τ )) (23)
where Ψ j (ξ(τ )) refers to a rearrangement of the p-order finite dimensional orthogonal polynomials with respect to the Gaussian function which constitute a complete basis in the space of second-order random variables ; Y j is the unknown deterministic j th vector associated with Ψ j (ξ(τ )) and ξ = {ξ r }, set of orthonormal random variables [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. Finally, after projection on the chaos polynomial basis, the stochastic system to solve is written as

AY = F (24)
where each component

[A] ij of A is [A] ij = Λ c + Λ d + L k=1 Λ rk Ψ 2 i δ ij + L k=1 Λ rk δ K k c ijk K (25)
Z corresponds to the mathematical expectation of quantity Z, δ ij is the Kronecker symbol and c ijk K is defined by [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]

c ijk K = ξ K k Ψ i Ψ j (26)
Finally, vector F(τ ) of forces projected on the polynomial chaos basis has polynomial chaos vector terms given by

F j (τ ) = Γ g δ 0j Ψ 2 i δ ij + F u j Ψ 2 i δ ij (27)
4 Numerical results

In this section, efficiency of the proposed methodology based on the combination of the Harmonic Balance Method with the Polynomial Chaos Expansion will be firstly demonstrated through numerical simulations of the rotor system with a breathing crack. Then, the random non-linear response of a cracked rotor against uncertain parameters and the effects of uncertainties from stiffness properties of the rotor and external forces will be analyzed in order to propose a robust identification of damage in rotors. The random aspect will be presented through the mean and the standard deviation of the response. The standard deviation σ r is defined by

σ r = (Θ(τ ) -Θ 0 ) 2 (28)
Given the orthogonality of the polynomial chaoses, Equation 28 can be rewritten as

σ r = ∞ j=1 Θ 2 i (τ ) Ψ 2 i (29)

Choice of the terms in the expansions

Clearly, in a computational context, the summations in the previous expressions are truncated after some term P . In particular, response Θ(τ ) given by Equation 23can be written as

Θ(τ ) = P j=0 Y j Ψ j (ξ(τ )) (30)
where P , the total number of polynomial chaoses used in the expansion (excluding the 0 th -order term), can be determined by P = 1 + p s=1 1 s! s-1 r=0 (L + r), p being the order of the homogeneous chaos being used [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. {Ψ j (θ)} j=0 to P refers to a rearrangement of the p th -order finite-dimension orthogonal polynomials with respect to the Gaussian function, which constitute a complete basis of the space of the second-order random variables. Besides, the uncertainty is quantified identically along the stiffness of the rotor, consequently the random modeling of K uses only one random variable. Thus, Equation 21 can be rewritten as

K(τ ) = K(1 + δ K ξ K (τ )) . (31)
It should be noted that the normal distribution assumed has been truncated to avoid negative values of the parameter, since it is physically strictly positive. Moreover, the expansion of unbalance force Γ u can reasonably be truncated to order 1 which is equivalent to consider the Karhunen-Loeve expansion. Then Equation 22 can be reduced to

Γ u (τ ) = Γ u (1 + δ Γu ξ Γu (τ )) , (32)
considering one random variable. Integrating these hypotheses in Equation 24, the system of equations to solve becomes

A Y = F (33)
where each component

[ A] ij of A is [ A] ij = Λ c + Λ d + Λ r Ψ 2 i δ ij + Λ r δ K c ij1 K (34)
and polynomial chaos vectors F j are given by

F j (τ ) = Γ g δ 0j Ψ 2 i δ ij + F u j Ψ 2 i δ 0j δ ij + δ Γu Ψ i (ξ) ξ Γu (35)
where ξ = {ξ K , ξ Γu }.

Finally, the standard deviation is also built from a truncated expansion that is

σ r = P j=1 Θ 2 i (τ ) Ψ 2 i /Θ 0 (36)
It should be noted that the size of the system of equations depends on P and is obviously higher than those of the deterministic problem.

Comparison between Monte-Carlo Simulations and Polynomial Chaos Expansion

The example of the non-linear vertical response at the middle of the cracked rotor, for one nondimensional crack depth µ equal to 1 and with variations for both shaft stiffness of the rotor system and the excitation forces of 2% is taken to show the efficiency and the robustness of the HBM-PCE method versus the HBM-MCS method. The reference solution is obtained with the HBM-MCS method using 5000 samples. We present the results from the HBM-PCE method using several orders in the expansion : p = 5, 10 and 20 (for listing of coefficients, see [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]). Here, for a better legibility, the mean (Figure 2) and the standard deviation (Figure 3) only for harmonic 2 of the random response are presented. It should be noted that similar conclusions were observed for the other harmonics of the response. Figures 2 and3 present each : (a) the response on the whole frequency range studied and (b) a zoom for example on the resonance second peak. Clearly, for both mean and standard deviation, we can observe that, outside the resonances, the HBM-PCE solution is identical to that obtained by the MCS whatever the order p. The more important area where differences between the HBM-PCE solution and the the HBM-MCS is located at each sub-critical speed. Figures 4 and5 show the mean and the standard deviation respectively for the four harmonics of the non-linear vertical response computed from both methods. Here again, the results for the random response are very close between both methods even around the resonances. These variations are only due the very small damping taken into account in the present study. For the interested reader, Figures 6 illustrated the non-linear vertical response with a higher damping (i.e. η = 10 -5 ). It clearly demonstrated that the difference between HBM-PCE and HBM-MCS solutions decreases with the value of the damping as also mentioned in PhD thesis of Dessombz [START_REF] Dessombz | Analyse dynamique de structures comportant des paramètres incertains[END_REF]. However, it is clearly seen that the HBM-PCE gets better with p, even if the damping is very small (i.e. η = 0.1 10 -5 ). Besides, we can stop to order 20 since it seems to be sufficient to represent correctly the random behavior of the response of the cracked rotor due to the uncertain inputs.

Efficiency of the crack detection based on the non-linear responses with uncertainties

In this first part of the parametric studies, the non-dimensional crack depth µ is equal to 1 (corresponding to the loss of half the shaft's area). We consider random variations for the stiffness properties of shaft of the rotor system (applied actually on the Young modulus) and for the excitation forces. Three cases are investigated in which the coefficients of variation of the shaft's stiffness and of the excitation forces are taken equal to 1%, 2% and 5% respectively. We recall that in this case the cost of calculation may be high since it is directly linked to the number of polynomials and consequently to the order of polynomial chaoses and the number of random parameters. Figures 7 show the mean value and the standard deviation of the non-linear vertical responses at the middle of the cracked rotor for these three cases. First of all, it appears that variations of the mean value and of the standard deviation increase when the rotor is passing through 1 n sub-critical speeds. For example, the 1× vertical response of the cracked rotor generates highly prominent amplitudes when the running speed is at the first critical speed and approximately 1 2 of the first critical speed. An increase of the 2× vertical response at the critical speed, 1 2 and 1 3 of the first critical speed can also be observed. Finally, the 3× and 4× vertical responses increase at the critical speed, 1 2 , 1 3 and 1 4 of the first critical speed. However it is clear that only 1× component is dominant at first critical speed and 2× component is very small. Similarly, 2× component is very small while 2× component is very dominant at 1 2 of critical speed. So, increases of the n× amplitudes at 1 m of the first critical speed (with m ≤ n) are clues for the detection of a transverse breathing crack in a rotating shaft, even if uncertainties are introduced in the rotor system. The same conclusions can be observed for the evolutions of the the mean value and the standard deviation of the non-linear horizontal response of the cracked rotor for the three cases (see Figures 8).

Moreover, it appears that increasing uncertainties of the rotor's stiffness and of the excitation forces affects the standard deviation of all the n× amplitudes when the rotor is passing through the 1 n subcritical speed. The range of speed where the non-linear amplitudes of the rotor system are predominant increase drastically with increasing uncertainties. This last observation clearly demonstrates that basic methods based on linear condition monitoring techniques (such as the changes in frequencies and modes shapes) can be rendered ineffective in the case of specific configurations taking into account uncertainties for the detection of cracks in rotor systems. However, whatever the levels and different kinds of uncertainties (such as material and/or loading characteristics) presented here, these numerical results demonstrate that emergences and variations in n× harmonic components of the system in the frequency domain are key indicators for the detection of transverse cracks in a rotating shaft. It is reminded for the reader comprehension that even if there is uncertainties in the value of stiffness property and unbalance, the evolution of damage will take a predictable path and the expected value of the n× amplitudes will be within a band of possible values due to uncertainties.

Finally Figures 9, 10,11 and 12 illustrate the mean values and the standard deviation of the nonlinear n× amplitudes (with n = 1, 2, 3 and 4) for a rotor system containing a small crack (i.e. µ = 0.5). The previous results (with µ = 1) are also given for comparison. In this last part of the study, the random stiffness properties of shaft (i.e. for the Young modulus) and the excitation forces have each a variation coefficient equal to 2%. As previously seen for µ = 1, the emergence of the n× amplitudes are clearly detected or µ = 0.5. So it may be concluded that using condition monitoring techniques based on non-linear approaches can be considered to be the first step for a rapid estimation of the presence of damage in rotors with structural uncertainties even if the crack size is small. Nevertheless, it is observed that the evolutions of the mean value and the standard deviation of the 1× harmonic components (2× and 3× harmonic components, respectively) are not significant at 1 2 of the first critical speed ( 13 and 1 4 respectively). Showing Figures 9, it is observed a large difference in the mean values of the 1× amplitudes for two different crack depths in horizontal rotor response while almost no difference in the vertical rotor responses. Then, considering Figures 11 or 12, variations of the mean values of the 3× and 4× amplitudes are significant for both the vertical and horizontal directions. These results illustrate interaction between the contributions of the unbalance and the breathing of the crack.

Finally, it is shown that increasing the crack size induces a small decrease of the frequency values of the critical and sub-critical speeds that may be observed even if uncertainties of the rotor's stiffness and of the excitation forces are introduced. However, these changes in frequencies are quite small and can not be used for a robust detection of the presence of damage in rotors.

Conclusion

This paper presents an extension of the use of Harmonic Balance Method in the stochastic domain using the Polynomial Chaos Expansion to calculate the random non-linear dynamical response of mechanical systems. Here, the proposed methodology is applied to investigate the effects of the presence of a transverse crack in a rotating shaft under uncertain physical parameters in order to outline some robust indicators for detecting damage in rotating system. We considered the stiffness properties of shaft of the rotor system and the excitation forces random by Gaussian law. The random non-linear response consists in the mean and the standard deviation of the non-linear response and is obtained by finding the deterministic coefficients associated to the Polynomial Chaos Expansion applied on the n× harmonics which are given from the Harmonic Balance Method. The comparison of the results from this method and the reference one, combining the simulations of Monte Carlo and the Harmonic Balance Method, for one case in particular, validates the proposed methodology. Then a complete study is done for different cases of uncertainties (from 0% to 5%) appearing on the stiffness properties of shaft of the rotor system and on the excitation forces to carry out the effect on the non-linear response. Besides, two values of depths were studied. For a practical point of view, using condition monitoring techniques based on linear approaches (i.e. frequency shifts and changes in modal parameters) can be considered to be unpractical and not efficient for an estimation of the presence of damage if uncertainties are present in the mechanical systems. So alternative indicators based on the non-linear signature of dynamical systems can be more effective for a robust detection of the presence of damage in complex mechanical engineering structures. Even if uncertainties are introduced in real engineering structures (more particularly in rotating applications), an important observation is the recognition of the appearances of the n× harmonic components when the running speed reaches 1 n of the critical speeds. However due to interaction between the unbalance and the breathing of the crack, the variations of the mean and standard deviation of amplitudes for two different crack depths can be more significant in the horizontal direction than in the vertical direction. Moreover, in some cases, no difference for two different crack depths should be observed in the mean and standard deviation of 1× amplitudes if the unbalance is predominant versus the breathing of the crack.

Figure 1 :

 1 Figure 1: Finite-element model of the rotor and the cracked-beam section

Figure 2 :Figure 3 :ωFigure 4 :ωFigure 5 :Figure 6 :

 23456 Figure 2: Mean of the second harmonic order computed by the MCS (solid line) and by the PCE of several orders : order 5 (dotted lines), order 10 (dotted-dashed lines), order 20 (dashed lines). (a) : whole frequency range, (b) : zoom on the second peak of resonance

Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 789101112 Figure 7: Evolution of the vertical n× displacements at the middle of the shaft for a non-dimensional crack depth µ = 1 with variations of the Young modulus and the excitation forces ( black=order 1, --blue=order 2 , -.red=order 3 , • • • green=order 4) (a) mean -1% variations (b) standard deviation -1% variations (c) mean -2% variations (d) standard deviation -2% variations (e) mean -5% variations (f) standard deviation -5% variations 18


Hippolyte D'albis

Emmanuel Thibault

Johanna Etner

Firouz Gahvari

Pierre Pestieau

Sébastien Pouget

Annuities, Bequests and Portfolio Diversification *

Keywords: Annuities, uncertain longevity, risk aversion. JEL codes: D11, D81, G11, G22

In this article, the diversification motives of the demand for annuities is analyzed. Using a model allowing for the uncertainty of both the human life length and the interest rate, the Decision Maker is supposed to choose an optimal portfolio to maximize a bequest. Conditions under which an increase in the risk of bond returns increase the demand for annuities are proposed and discussed. Moreover, it is shown that, contrary to previous claims, more risk adversion is associated with a lower demand for annuities.

Introduction

Since [START_REF] Yaari | Uncertain lifetime, life insurance, and the theory of the consumer[END_REF], the demand for annuities has become the cornerstone of the theory of consumption under uncertain lifetime. Provided that annuities are fairly priced, it has been shown that their demand should be relatively high to finance the last period's consumption. The remainder of the portfolio may then be composed of regular bonds, which, if the Decision Maker (DM) exhibits some bequest motives, are intended for her heirs. As shown in [START_REF] Davidoff | Annuities and individual welfare[END_REF], this result holds under quite general specifications and is therefore at odds with most empirical studies. This is known to be the annuity puzzle.

In most of the literature, a rather specific financial environment with no other uncertainty than life duration is assumed. The annuity is therefore a risky asset, whose stochastic yield is compared to a risk free interest rate. We claim, that introducing a stochastic interest rate strongly modifies the demand for annuities, as this becomes an instrument for portfolio diversification. The problem we consider is the following. At each period of life, annuities provide deterministic returns if the DM is alive and zero returns if she is dead. Conversely, bonds (or stocks) yield an uncertain return, which is independent of their owner's survival probability. We propose the simplest possible model to analyze the diversification problem of a DM that faces these two types of risk. The simplicity lies on the assumption that only the bequest yields some utility, or equivalently, that there is no consumption in the case of survival. This eases the comparison with the literature, since the optimal demand for annuities is zero if bond returns are deterministic. In such a setting, a positive demand for annuities is only due to the uncertainty of the interest rate.

We show that, in this stochastic financial environment, the demand for annuities may be positive. Annuities are purchased to diversify the portfolio, but their demand does not necessarily increase with the risk of bond returns. We show that a meanpreserving spread on bond returns increases the optimal demand for annuities if the DM is prudent, but not too prudent. We therefore extend the initial result proposed by [START_REF] Hadar | The effects of shifts in a return distribution on optimal portfolios[END_REF] to demand for annuities. Within this framework, we also study the impact of risk aversion on the demand for annuities. We provide a sufficient condition, under which the demand for annuities is reduced when risk aversion is increased.

The condition, based only on prices, simply ensures that the optimal utility if alive is greater than the utility if not alive. This is not true for all sets of parameters, since the bequest yields some utility. Under that condition, a more risk adverse DM chooses to purchase fewer annuities to reduce the gap between the utilities in the two states of nature, life and death. This result puts into perspective the previous studies results [START_REF] Blake | Annuity markets: Problems and solutions[END_REF][START_REF] Friedman | The cost of annuities: Implications for saving behavior and and bequest[END_REF][START_REF] Milevsky | Optimal asset allocation toward the end of the life cycle[END_REF][START_REF] Warshawsky | Private annuity markets in the US[END_REF]) that claim that demand for annuities increases with risk aversion. Section 2 presents the basic static framework, while the results lie in Section 3. In Section 4, the model is extended to a dynamic setting with savings.

The basic framework

We consider a portfolio choice model under uncertain lifetime in which the Decision Maker (DM) can invest in two assets to maximize her bequest utility. To focus on the diversification issue, the problem is static and there is no utility derived from consumption. In Section 4, this basic framework is extended to a dynamic problem with endogenous savings.

The DM faces two independent risks: the survival and the return of one of the two assets. The length of life is at most, two periods, but only the second one is uncertain.

Life uncertainty is characterized using the random variable x which follows a Bernoulli law whose expectancy is E x (x) = p ∈ (0, 1). The bequest might therefore happen at the end of Periods 1 and 2.

At the first period, the DM is endowed with a positive initial income ω 1 that can be shared between bonds1 and annuities. The returns in Period 2 of bonds purchased Period 1 are given by the random variable r with support on [r -, r + ], where 0 < r -≤ r + < ∞, and whose expectancy is E r (r) = r. Bond returns are paid to the DM if alive and to her heir if not alive. Conversely, the annuities return is r/p in Period 2 if the DM is alive and nothing if she is not alive.

Following [START_REF] Yaari | Uncertain lifetime, life insurance, and the theory of the consumer[END_REF], annuities returns are fair. However, we suppose that annuity sellers can eliminate risk at no cost. It follows that annuities and bonds are two risky assets that have the same expected return:

E x x r p = E r (r) = r.
Conditional to survival, investments in annuities on average are more profitable that investments in bonds. This is nevertheless not true when the realizations of r are larger than r/p. As a consequence, holding annuities may permit an increase in wealth that will be bequeathed: annuities can thus be purchased even if the DM does not derive utility from her consumption. Remark that it could be possible to extend our framework by assuming that bonds have a higher expected return than annuities, e.g., due to diversification costs supported by annuity sellers. If the difference in returns is small enough, our results are not modified.

If alive in Period 2, the DM receives a non-negative income ω 2 and bequeaths her entire wealth to her heir. Since death is certain at the end of Period 2, the bequest is exclusively a demand for bonds. More simply, it is supposed that, between Periods 2 and 3, there is no uncertainty on bond returns, which are normalized to 1. In addition, whatever the length of the DM's life, bequests are received in Period 3. We denote respectively by a and ω 1 -a the demand for annuities and the demand for bonds in Period 1. As of Period 1, the wealth that will be bequeathed in Period 3 hence satisfies:

x r p a + ω 2 + r (ω 1 -a) .

sumptions in Periods 1 and 2 are fixed. Then, ω 1 and ω 2 represent the exogenous difference between income and consumption in each period. Note also that the case of endogenous savings is considered in Section 4.

Due to uncertain life spans, it is not possible to borrow by selling bonds. It is also assumed that short sale positions on annuities are not possible. Thus, the demand for annuities and the demand for bonds are supposed to be non negative:

0 ≤ a ≤ ω 1 . (1)
Finally, the utility derived from the bequest is computed using the function u (.), which is C 3 , satisfies u ′ > 0, u ′′ < 0, as well as the usual limit conditions: lim y→0 u ′ (y) = +∞ and lim y→+∞ u ′ (y) = 0. The DM is an expected utility maximizer and her objective is therefore given by:

max a E x E r u x r p a + ω 2 + r (ω 1 -a) . (2)
The DM's problem is to solve (2) subject to (1).

Annuities and portfolio diversification

This section is devoted to the analysis of the optimal portfolio. The first proposition establishes the condition under which the DM invests in annuities. After which, the impact of risk and risk aversion are studied.

Proposition 1 -There exists a unique optimal demand for annuities.

If ω 2 = 0, the optimal demand for annuities is positive.

If ω 2 > 0, the optimal demand for annuities decreases with ω 2 when the DM is prudent. Moreover, there exists ω2 > 0 such that the optimal portfolio is composed of bonds and annuities, if ω 2 ≤ ω2 and only of bonds if ω 2 > ω2 .

Proof -As a preliminary, we define the function φ (a) as the first derivative of (2) with respect to a, such that:

φ (a) = E x E r x r p -r u ′ x r p a + ω 2 + r (ω 1 -a) . (3)
This function is well defined for a, such that, whatever the state of nature, the realization x (ra/p + ω 2 ) + r (ω 1 -a) is non-negative. As a result, a is lower than the positive threshold a max = ω 1 which is necessary and sufficient to guarantee the positivity of r(ω 1 -a). This is larger than the negative threshold a min = -(ω 2 + r -ω 1 )/(r -pr -), which is necessary and sufficient to guarantee the positivity of ra/p + ω 2 + r (ω 1 -a). Subsequently, φ(.) is well defined for a ∈ (a min , a max ).

By construction of these thresholds and using the Inada conditions we have:

lim a→a min φ(a) = lim ̺→0 (r -pr -)u ′ (̺) = +∞, lim a→amax φ(a) = lim ̺→0 E r [-ru ′ (r̺)] = -∞.
Consequently, φ(.) has a real root â ∈ (a min , a max ). Since

φ ′ (a) = E x E r x r p -r 2 u ′′ x r p a + ω 2 + r (ω 1 -a)
is negative, â is the unique real root of φ(.).

The solution of the DM's problem is thus denoted as a * and satisfies:

a * = â, if â ≥ 0, 0 otherwise. ( 4
)
We now establish the positivity of â, when ω 2 = 0. Since u(.) is concave and lim ̺→0 u ′ (̺) = +∞, we have:

lim a→0 φ(a)| ω 2 =0 = E r [(r -r)u ′ (rω 1 )] = cov[r -r, u ′ (rω 1 )] > 0, lim a→amax φ(a)| ω 2 =0 = lim ̺→0 E r [-ru ′ (r̺)] = -∞.
The continuity of φ(.) implies that the unique root â belongs to (0, ω 1 ) when ω 2 = 0.

When ω 2 is positive, by applying the implicit function theorem to (3), we can establish that:

dâ dω 2 = E r [(r -pr) u ′′ (râ/p + ω 2 + r (ω 1 -â))] -φ ′ (â) .
The numerator of the RHS can now be rewritten as:

p cov r -r, u ′′ r p â + ω 2 + r (ω 1 -â) + (1 -p) rE r u ′′ r p â + ω 2 + r (ω 1 -â) .
This numerator is negative if u ′′′ (.) > 0, or equivalently, if the DM is prudent. The denominator is positive. Hence, the optimal demand for annuities decreases with ω 2 (i.e., dâ/dω 2 < 0).

Finally, using Inada conditions, we can prove the existence of the positive threshold ω2 by computing the following limit:

lim ω 2 →+∞ φ (a) = -(1 -p) E r [ru ′ (r (ω 1 -a))] .
Since this limit is negative, there exists, by continuity, a threshold ω2 above which the optimal portfolio is only composed of bonds (i.e., a * = 0).

The optimal portfolio always includes bonds, as it is the only possibility to bequeath in case of death after the first period. The assumption of an infinite marginal utility when the bequest goes to zero indeed ensures a positive demand for bonds.

Conversely, investing in annuities may increase the bequest value in case of survival, but this is not necessary. Proposition 1 states that diversification using annuities arise if ω 2 , the income received if alive in the second period, is low enough.

The intuition is the following. For a given realization of r, an increase in ω 2 reduces the marginal utility if alive. At the optimum, this is compensated by an increase in the demand for bonds, meaning a reduction of the marginal utility if not alive. This additional demand for bonds corresponds to a reduction in the demand for annuities.

Nevertheless, the concavity of the utility function is not sufficient to obtain our result. Indeed, for large realizations of r, i.e. those which are greater than the annuities return r/p, a reduction in the demand for annuities increases the utility if alive. A prudent behavior, as characterized by [START_REF] Kimball | Precautionary savings in the small and in the large[END_REF], is then a sufficient condition for this latter effect to be dominated. It indeed implies the concavity of the marginal utility if alive with respect to r.

Further assumptions on the utility function permit us to derive additional results.

In the following corollary, this is notably the case, when the preferences of the DM are homothetic. Such preferences can be represented by a homogeneous utility function, such as the CRRA function.

Corollary 1 -When preferences are homothetic, the demand for annuities is linearly increasing in ω 1 and linearly decreasing in ω 2 .

Proof -Homothetic preferences can be represented by a concave utility function u(y) = y λ with 0 < λ < 1. According to Mas-Colell et al (1995, p. 50) this is a necessary and sufficient condition for representing homothetic and continuous preferences.

Moreover, for homothetic preferences, the DM's problem can be simplified by a change in variable. Hence, we define k such that:

a = kω 1 -ω 2 r/p + k ,
with k = -r/p. Then, (3) can be rewritten as a function of k as follows:

φ kω 1 -ω 2 r/p + k = E x E r x r p -r u ′ (xk + r) rω 1 + pω 2 r + pk .
As the function u ′ (.) is homogeneous of a given degree, the real roots of φ (.) are those of ϕ (.), where:

ϕ (k) = E x E r x r p -r u ′ (xk + r) . ( 5
)
2 When the indifference curves are homothetic with respect to the origin, we say that the preferences are homothetic. Such preferences can be represented by the composition of a function homogeneous of degree 1 with an increasing function.

Using the Inada conditions we have:

lim k→+∞ ϕ(k) = (1 -p)E r [-ru ′ (r)] < 0.
Importantly, we also have:

ϕ(0) = E r [(r -r)u ′ (r)] = cov[r -r, u ′ (r)] > 0.
Concerning the derivative of ϕ(.) we find that:

ϕ ′ (k) = E r [(r -pr)u ′′ (k + r)] = p cov (r -r, u ′′ (k + r)) + (1 -p) rE r [u ′′ (k + r)] . Since u ′′′ (y) = λ(λ -1)(λ -2)y λ-3 > 0, ϕ(.) is a decreasing function of k. As ϕ(0) > 0, ϕ ′ (k)
< 0 and lim k→+∞ ϕ(k) < 0, the function ϕ(.) has a unique real positive root, k. Then, the optimal demand for annuities, denoted a * , satisfies:

a * = kω 1 -ω 2 r/p + k , if k ≥ ω 2 /ω 1 , 0 otherwise. ( 6
)
Since k is independent of ω 1 and ω 2 , one can conclude that a * is linearly increasing in ω 1 and linearly decreasing in ω 2 .

Proposition 1 and Corollary 1 reveal that the demand for annuities may be positive, even if there is no utility derived from the consumption in Period 2. As shown in the two following propositions, this demand can be explained by the uncertainty on bond returns.

Proposition 2 -If the DM is prudent, the demand for annuities increases with the survival probability.

Proof -Suppose that the demand for annuities a * is positive. According to the proof of Proposition 1, we have a * = â. To perform comparative statics, it is useful to rewrite (3) as follows:

φ (a) = E r (r -pr) u ′ r p a + ω 2 + r (ω 1 -a) -(1 -p) ru ′ (r (ω 1 -a)) , (7)
which will be denoted E r [Γ(r)] for convenience.

From [START_REF] Hadar | The effects of shifts in a return distribution on optimal portfolios[END_REF], and using the implicit function theorem, we have:

dâ dp = -E r r u ′ (r (ω 1 -â)) -u ′ r p â + ω 2 + r (ω 1 -â) + râ p 2 E r (r -pr) u ′′ r p â + ω 2 + r (ω 1 -â) φ ′ (â).
The first term of the numerator is negative by the concavity of u (.), while the second one is negative if u ′′′ (.) > 0. Thus, if the DM is prudent, the numerator is negative and we have da * /dp > 0.

Proposition 2 is rather intuitive. An increase in the survival probability reduces the annuities return if alive but not its expected return, which is still equal to the return of the bonds. However, increasing the survival probability provides more weight to the marginal utilities if alive which is compensated by an increase in the demand for annuities.

As in Proposition 1, the assumption of prudence is sufficient to eliminate counterintuitive behaviors. Moreover, the limit cases of determinisitic survival and death are useful to understand why annuities are purchased.

In the limit case, such that p = 1, (7) becomes:

φ(a)| p=1 = cov r -r, u ′ (ra + ω 2 + r(ω 1 -a)) > 0.
As the concavity of u(.) ensures its positivity, φ (.) | p=1 has no real root and the optimal solution is a corner solution, a * = a max = ω 1 .

In the limit case, such that p = 0, there is no annuity market and the objective function (2) can be rewritten as: E r u (r (ω 1 -a)) . In this case, the optimal solution is a corner solution: a * = 0.

In both limit cases, there is by definition no annuity market, but they permit to understand the logic behind Proposition 2. If survival is certain, the "annuity" appears to be a standard non-risky asset whose return equals the expected return of the bonds.

Because of risk aversion, the optimal portfolio is then composed of non-risky assets only. Conversely, if death is certain, there is no non-risky asset and all the initial endowment are invested in bonds, whose minimal return, r -, is greater than zero.

The next proposition studies the effect of uncertainty on bond returns. To begin, we define the relative prudence as:

RP = - yu ′′′ (y) u ′′ (y) .
Proposition 3 -If r + = r -, the optimal demand for annuities is zero.

If r + > r -, a mean-preserving spread on bond returns increases the optimal demand for annuities if the relative prudence is positive and less than 2.

Proof -When there is no risk on bond returns, equation (3) can be rewritten as:

φ (a)| r=r = (1 -p) r u ′ r p a + ω 2 + r (ω 1 -a) -u ′ (r (ω 1 -a)) .
The unique root of φ (.) | r=r , â = -pω 2 /r, is non-positive and, according to (4), the optimal solution is then a * = 0.

Suppose now that â > 0. Consider a mean-preserving spread on r. Following Rothschild and Stiglitz (1971) and using [START_REF] Hadar | The effects of shifts in a return distribution on optimal portfolios[END_REF], it increases â if Γ(r) is convex. Differentiating twice Γ(.) and rearranging the equation yields:

Γ ′′ (r) = -2 (ω 1 -a) E x u ′′ x r p a + ω 2 + r (ω 1 -a) -(ω 1 -a) E x x r p a + ω 2 + r (ω 1 -a) u ′′′ x r p a + ω 2 + r (ω 1 -a) + (ω 1 -a) (ω 1 r + pω 2 ) u ′′′ r p a + ω 2 + r (ω 1 -a) .
Conclude that Γ(r) is convex when, for all y > 0, u ′′′ (y) ≥ 0 and RP ≤ 2.

Whatever the uncertainty of the length of life, there is no demand for annuities if there is no risk on bond returns. In this case, the optimal behavior simply aims at equalizing the utilities in the two states of nature: life and death. Since there is no consumption if alive, the optimal demand for annuities is zero. Uncertainty on bond returns is, in this framework, a necessary condition for annuitization. This motive complements the one traditionally studied in the annuity literature, which relies on consumption.

The impact of a change in risk on the demand for annuities was also discussed in Proposition 3. As bonds become more risky, the DM optimally diversifies her portfolio by increasing the share of annuities. To obtain such a behavior, sufficient conditions on preferences are exhibited. The DM has to be prudent, but not relatively too prudent, with a benchmark value at 2. This result is similar to Hadar and Seo (1990) who studied an optimal portfolio problem with two risky assets. It is also generalized in [START_REF] Gollier | The Economics of Risk and Time[END_REF]. Interestingly, [START_REF] Eeckhoudt | A benchmark value for relative prudence[END_REF] proposed an interpretation of the benchmark value in terms of preferences for disaggregating harms.

Finally, note that the conditions exhibited in Proposition 3 are satisfied if the relative risk aversion is both lower than one and non-decreasing. The following example reveals that the demand for annuities may decrease with the risk on bonds if the relative risk aversion is greater than one.

Assume that the random variable r has only two realizations: r + ≡ r + ε/q with probability q and r -≡ r -ε/ (1 -q) with probability 1 -q. Hence, r is the expected return and ε ≥ 0 is a mean-preserving spread measure. Assume that the utility function is CRRA: u (y) = y 1-α / (1 -α), where α > 0 stands for the relative risk aversion coefficient. As u(.) is homothetic, the optimal demand for annuities a * is defined by [START_REF] Gollier | The Economics of Risk and Time[END_REF] where k is given by (5) (see Corollary 1).

Obviously we have da * /dε = da * /d k × d k/dε. Since the first term of the RHS is positive for all a * > 0, the sign of the effect of an increase in risk on the demand for annuities can be evaluated with d k/dε. According to (5):

d k dε = - ϕ ′ ε ϕ ′ k = E x u ′ (xk + r + ) - x r p -r + u ′′ (xk + r + ) -E x u ′ (xk + r -) - x r p -r -u ′′ (xk + r -) E r (r -pr)u ′′ (k + r) .
If the individual is prudent and exhibits a relative risk aversion greater than one, the demand for annuities increases with the realization of r. The effect of a larger r + is then the opposite of the effect of a larger r -. If relative prudence is strong enough, it is not clear which effect dominates. This is illustrated in Figure 1, which plots a * as a function of ε for the following parameter values: α = 4, p = 0.9, q = 0.8, r = 1.5, ω 1 = 1 and ω 2 = 0. Figure 1 shows that the demand for annuities is reduced by an increase in risk when the meanpreserving spread is sufficiently high. This example has been computed for a relative risk aversion of 4, which is standard. Nevertheless, it is a crucial parameter. In Figure 2, a * as a function of ε is plotted for various values of α.

The impact of risk aversion on the demand for annuities will be discussed in Propo- Proof -Consider two DMs, namely A and B. We assume that the utility function of A is u A (.), whereas the utility function of B is u B (.) ≡ T • u A (.) where T (.) is an increasing and concave function. In this case, the DM B is more risk averse than the DM A. 3 Each DM i (i = A, B) maximizes their expected utility

E x E r [u i (x (ra/p + ω 2 ) + r (ω 1 -a))
] subject to (1). We denote by a i ⋆ (i = A, B) their optimal choices, which are the real roots of φ i (a) if this latter is positive, or zero otherwise. One has:

φ i (a) = E x E r x r p -r u i ′ x r p a + ω 2 + r (ω 1 -a) .
Assume that a A ⋆ and a B ⋆ are both positive. We now exhibit a condition such that

φ B a A ⋆ < 0. Such a condition implies that a B ⋆ < a A ⋆ because φ B (.) is a decreasing function satisfying φ B a B ⋆ = 0.
Observe first that φ B a A ⋆ can be rewritten as follows:

φ B a A ⋆ = E r T ′ u A r p a A ⋆ + ω 2 + r ω 1 -a A ⋆ ψ A 1 (r) -E r T ′ u A (r (ω 1 -a)) ψ A 2 (r) ,
where:

ψ A 1 (r) = (r -pr) u A ′ ra A ⋆ /p + ω 2 + r ω 1 -a A ⋆ , ψ A 2 (r) = (1 -p) ru A ′ r ω 1 -a A ⋆ .
As ψ A 2 (r) is positive we have:

φ B a A ⋆ < E r T ′ u A r p a A ⋆ + ω 2 + r ω 1 -a A ⋆ ψ A 1 (r) -T ′ u A (r + (ω 1 -a)) E r [ψ A 2 (r)]. (8)
Importantly, φ A (a A ⋆ ) = 0 implies that:

E r (ψ A 1 (r)|r < r/p) Pr(r < r/p) = E r ψ A 2 (r) -E r (ψ A 1 (r)|r ≥ r/p) Pr(r ≥ r/p).
Moreover, the concavity of T (.

) implies that E r [T ′ (u A (ra A ⋆ /p+ω 2 +r(ω 1 -a A ⋆ )))ψ A 1 (r)]
is lower than:

T ′ (u A (r -a A ⋆ /p + ω 2 + r -(ω 1 -a A ⋆ )))E r [ψ A 1 (r)|r < r/p] Pr(r < r/p) +T ′ (u A (r + a A ⋆ /p + ω 2 + r + (ω 1 -a A ⋆ )))E r [ψ A 1 (r)|r ≥ r/p] Pr(r ≥ r/p).
Then, according to (8), we have:

φ B a A ⋆ < T ′ u A r p a A ⋆ + ω 2 + r -ω 1 -a A ⋆ -T ′ u A r + ω 1 -a A ⋆ E r ψ A 2 (r) + T ′ u A r p a A ⋆ + ω 2 + r + ω 1 -a A ⋆ -T ′ u A r p a A ⋆ + ω 2 + r -ω 1 -a A ⋆ ×E r ψ A 1 (r) r ≥ r/p Pr (r ≥ r/p) .
Since T (.) is concave, the second term of the RHS of this inequality is negative. A sufficient condition for φ B (a A ⋆ ) < 0 is therefore:

T ′ u A r p a A ⋆ + ω 2 + r -ω 1 -a A ⋆ < T ′ u A r + ω 1 -a A ⋆ ,
which can be rewritten as:

r p a A ⋆ + ω 2 -(r + -r -) ω 1 -a A ⋆ > 0.
Since the RHS of this inequality is increasing with a A ⋆ , condition ω 2 ≥ (r + -r -) ω 1 is sufficient.

Proposition 4 states that more risk aversion may increase the optimal share of the portfolio invested in bonds. The sufficient condition given only relies on parameters and is therefore true for any increasing and concave utility function.

The condition is simply that the endowment growth factor should not be lower than the support of the distribution of the interest factor. Were bonds supposed to be riskless, the condition would be satisfied. More precisely, this condition ensures that, when computed at the optimal point, the utility if alive in Period 2 is always greater than the utility if not alive. And consequently, living with the lowest possible interest rate bears more utility than being dead with the highest bond returns. Note that this is not necessarily the case since it is the bequest which yields some utility.

If our condition is satisfied, more risk aversion induces to reduce the difference in utility for each realization of r. A more risk adverse DM therefore lowers her demand for annuities, which increases the utility if dead and decreases the utility if alive. Note finally that our condition is sufficient and, as shown in the following example, that a positive income in Period 2 is not necessary. We now consider the example used previously and plot the demand for annuities a * as a function of the relative risk aversion coefficient α, for the following parameter values: p = 0.9, q = 0.8, r = 1.5, ω 1 = 1, ω 2 = 0, and an interest rate spread of ε = 0.1.

Despite the assumption of ω 2 = 0, the demand for annuities typically monotonically decreases with relative risk aversion. However, this is not always true. For instance, by reducing the interest rate spread, the relationship may be reversed for some value of α. As shown in Proposition 3, a reduction of the spread ε may indeed lower the demand for annuities, and therefore imply, for some realizations of r, that the utility if not alive is higher than the utility if alive. As shown in Figure 4, the relationship between the demand for annuities and the risk aversion may become positive (case where ε = 0.05).

Extension to savings

This section introduces a dynamic behavior allowing for a first period consumption.

The assumption of no consumption if alive in the second period (or, equivalently, of a fixed consumption) is kept to show that the demand for annuities is still positive, provided that bond returns are risky.

The framework is the one described Section 2 extended to a first period consumption. Before the lotteries, the DM may allocate her initial income ω 1 > 0 between consumption and savings, which are respectively denoted ω 1 -s and s. The utility derived from consumption is described by the function v (.) that satisfies the same assumptions as those of function u (.). The DM's program can now be rewritten as:

max s,a v (ω 1 -s) + E x E r u x r p a + ω 2 + r (s -a) (9)
subject to 0 ≤ s ≤ ω 1 and 0 ≤ a ≤ s.

Proposition 5 -If r + = r -, there exists a unique optimum (s * , a * ) where s * ∈ (0, ω 1 )

and a * = 0. If r + > r -, there exists ω2 > 0 such that if the DM is prudent and if ω 2 < ω2 , there exists a unique optimum (s * , a * ) ∈ (0, ω 1 ) 2 .

Proof -First, we establish the existence and uniqueness of the optimum (s * , a * ), when the bond return is random. The first order conditions of problem ( 9) are:

-v ′ (ω 1 -s * ) + E x E r ru ′ x r p a * + ω 2 + r (s * -a * ) = 0, (10)
E x E r x r p -r u ′ x r p a * + ω 2 + r (s * -a * ) = 0. (11)
From [START_REF] Milevsky | Optimal asset allocation toward the end of the life cycle[END_REF], we can define the following function ξ a (s):

ξ a (s) = -v ′ (ω 1 -s) + E x E r ru ′ x r p a + ω 2 + r (s -a) .
The function ξ a (.) is well defined for all s, such that x (ra/p + ω 2 ) + r (s -a) is non-negative. ξ a (.) is well defined for s ∈ (a, ω 1 ) and using the Inada conditions:

lim s→a ξ a (s) = lim ̺→0 (1 -p)E r [ru ′ (r̺)] = +∞, lim s→ω 1 ξ a (s) = lim ̺→0 -v ′ (̺) = -∞.
As a result, ξ a (.) has at least one root s * that belongs to (a, ω 1 ). According to Proposition 1, it is straightforward that for any root s * , there exists a positive threshold ω2 such that (11) has one real root a * that belongs to (0, s * ), if ω 2 < ω2 .

To prove the uniqueness of the pair (s * , a * ), we first combine [START_REF] Milevsky | Optimal asset allocation toward the end of the life cycle[END_REF] and [START_REF] Pratt | Aversion to one risk in the presence of others[END_REF] to obtain:

γ(s * , a * ) = -v ′ (ω 1 -s * ) + rE r u ′ r p a * + ω 2 + r (s * -a * ) = 0.
From γ(s * , a * ) = 0, we can use the implicit function theorem to define a continuous application a * = η (s * ) where:

η ′ (s * ) = - γ ′ s (s * , a * ) γ ′ a (s * , a * ) = - v ′′ (ω 1 -s * ) + rE r [ru ′′ (ra * /p + ω 2 + r (s * -a * ))] rE r [(r -pr) u ′′ (ra * /p + ω 2 + r (s * -a * ))] /p , (12)
As γ ′ a (.) = r cov[r -r, u ′′ (ra * /p + ω 2 + r(ω 1 -a * ))] + (1/p -1)r 2 E r u ′′ (ra * /p + ω 2 + r(ω 1 -a * )) , the denominator of ( 12) is negative when u ′′′ (.) > 0. Since γ ′ s (.) < 0, η(.)

is a decreasing function when the DM is prudent. Now, replacing η (s) in [START_REF] Milevsky | Optimal asset allocation toward the end of the life cycle[END_REF], we can define the function ξ (s) as follows:

ξ (s) = ξ η(s) (s) = -v ′ (ω 1 -s) + E x E r ru ′ x r p η (s) + ω 2 + r s -η(s) .
We have:

ξ ′ (s) = v ′′ (ω 1 -s) + E x E r r2 u ′′ x r p -r η (s) + xω 2 + rs +η ′ (s) E x E r r x r p -r u ′′ x r p -r η (s) + xω 2 + rs .
Next, a sufficient condition for ξ ′ (s) < 0 is:

-η ′ (s)E r [r(r -pr)u ′′ (z)] > v ′′ (ω 1 -s) + pE r [r 2 u ′′ (z)],
or, equivalently, using ( 12):

E r [r (r -pr) u ′′ (z)] rE r [(r -pr) u ′′ (z)] /p < v ′′ (ω 1 -s) + pE r [r 2 u ′′ (z)] v ′′ (ω 1 -s) + rE r [ru ′′ (z)] ,

where z = rη/p + ω 2 + r (s -η).

The condition can be rewritten as:

- 1 p v ′′ (ω 1 -s)E r [(r/p -r) 2 u ′′ (z)] < (r/p) 2 {E r [r 2 u ′′ (z)]E r [u ′′ (z)] -(E r [ru ′′ (z)]) 2 }.
According to the Cauchy-Schwartz inequality, the RHS of this inequality is positive.

Subsequently, the condition is satisfied and ξ ′ (s) < 0. Consequently, ξ(s) has at most one real root s * and the pair (s * , a * ) is unique.

Finally, we focus on the riskless case. When there is no risk of bond returns (i.e.

for r = r), [START_REF] Pratt | Aversion to one risk in the presence of others[END_REF] can be rewritten as follows:

(1 -p) u ′ r p a + ω 2 + r(s -a) -u ′ r(s -a) = 0.

Then, one may explicitly compute the solution, â = -pω 2 /r, which is negative.

The optimal demand for annuities is therefore: a * = 0.

Proposition 5 shows that the results contained in Propositions 1 and 3 still hold when a saving motive is introduced. The following numerical simulations also show that, as claimed in Proposition 4, the demand for annuities may still decrease with risk aversion.

We extend the basic example studied in Section 3 with the following utility function for the first period consumption: v (y) = y 0.6 /0.6. Using the parameter values of Figure 3, Figure 5 plots the demand for annuities a * and the savings as functions of the relative risk aversion coefficient α. The demand for annuities still monotonically decreases with relative risk aversion, despite an increase in saving.

However, considering a lower spread ε (0.05 rather than 0.1) is sufficient to obtain a non-monotonic relationship. Figure 6 plots the demand for annuities a * and the savings as functions of the relative risk aversion coefficient α, with a lower spread. In this article, we have studied the demand for annuities when the alternative asset has a risky return. We have provided conditions under which the demand for annuities increases with the financial risk and reduces with risk aversion. An important assumption of this work is the DM's perfect knowledge of the probability distributions. It would be interesting to instead suppose some ambiguity in the probabilities and study the DM's attitude when facing those two different sources of uncertainty.

  Figure 1

Figure 2 Proposition 4 -

 24 Figure 2

  Figure 3

  Figure 4

  Figure 5

  Figure 6

For simplicity, we call the "bond" the risky asset.

This result was stated in[START_REF] Pratt | Aversion to one risk in the presence of others[END_REF] who showed that concavity is preserved under mixture of independent risks (see also[START_REF] Finkelshtain | On risk aversion with two risks[END_REF]).

support of the Europlace Institute of Finance is gratefully acknowledged.