A Delga

L Doyennette

A Buffaz

V Berger

F R Jasnot

N Péré-Laperne

L A De Vaulchier

H C Liu

Space charge mediated negative differential resistance in terhertz quantum well detectors

Keywords:

In THz Quantum Well Infrared Photodetectors, a built-in-charge-mediated regime transition of the electronic transport is thoroughly investigated. The very strong current discontinuity and negative differential resistivity behavior are explained in terms of band structure reorganizations. The analysis of bias versus current measurements reveals that the transition occurs when the first two wells of the structure become partially drained, and that second well enters the ionized regime before the first one. Both many body effects and a careful model of the contact have to be considered to account for these features. The probable source of the built-in charge is identified as Intersubband Impact Ionization, of which this regime transition would be one of the few experimental evidence, and provide an original approach to investigate hot electrons kinetics in these structures.

I. INTRODUCTION

In the electromagnetic spectrum, the terahertz region (1-10 THz) is under investigation by both optics and electronics communities 1,2 . THz efficient sources have shown a fast development, thanks to the rise of THz quantum cascade lasers [3][4][5] . However there is still a need for fast and high sensitive detectors, in spite of the investigation of many different devices, such as quantum dots 6 , charge sensitive phototransistors 7 , superconducting bolometers 8 or blocked-impurity-band detectors 9 . Although traditionally used in the near or mid-infrared regions, Quantum Well Infrared Photodetectors structures (QWIPs) have been recently investigated as another valid approach 10 . Built-in charge effects are a key factor controlling electron kinetics in semiconductors. They have a crucial -positive or negative-impact in a broad range of devices, such as double-barrier heterostructures 11 , avalanche photodiodes 12,13 , or solar cells 14 . In QWIPs optical gain comes from such built-in charge effects, but it has also been shown in the THz spectrum that they trigger a regime transition where the structure switches from a 'down' resistive state to an 'up' conductive state with a current discontinuity up to five orders of magnitude 15 . In this work we focus on a THz QWIP structure and we thoroughly analyze the impact of the built-in positive charges in the quantum wells with the device macroscopic behavior.

II. EXPERIMENTAL FACTS

The studied sample is a THz QWIP structure (9.7 THz) as in Ref. 10. It contains 40 GaAs quantum wells (QWs) of thickness L w = 11.9 nm separated by wide Al 0.05 Ga 0.95 As barriers, with L b = 55.2 nm. Each QW is doped with a nominal value of n 2D = 10 11 cm -2 , and is designed to put in resonance both the second bound state and the bottom of the conduction band in the barriers [START_REF] Schneider | Quantum well infrared photodetectors: physics and applications[END_REF] . The structure is embedded between two contacts doped up to 10 17 cm -3 . The performances of this detector (referred as v265) were studied by Luo et al. 10 between 4 and 23 K.

Figure 1 shows I(v) and V (i) dark characteristics measured at a temperature of 4 K. On the whole I(v) curve, the bias is swept upwards then downwards, and the threshold values for the shifting up (0.84 V) is higher than the one for the shifting down (0.65 V). The hysteresis pattern is shown more clearly in the expanded scale plot. While increasing the bias from zero, the current is first dominated by interwell tunneling 15 . Temperature and bias are too low to allow electrons to flow through the continuum of states above the barriers. These channels are opened as the bias goes up, when the threshold value of 0.84 V is reached. At this stage, a sudden magnification of the process draining the QW occurs. Note that at this value, electrons injected at the contact Fermi level will have enough energy to ionize bound electrons from the wells by Coulomb interaction. A rough approximation will locate the hot electrons at an energy ∆E above the top of the barriers, where ∆E = E 2 -E 1 ∼ 29 meV is the transition energy between bound and resonant levels in the QW. Because the capture process is not en- hanced in the same way, some extracted electrons will leave unscreened donor impurities in the well, resulting in a built-in positive charge and an electric field discontinuity. From this point upwards the I(v) and V (i) curves show different behaviors. If bias is imposed, the energy difference between the two contacts is clamped. The only way for the structure to satisfy both constraints -bias and field discontinuity-is to breakdown the upstream barrier of the ionized well. This leads to a large shift of the current and the switching of the sample from the 'down' to the 'up' regime, where the current flows through the continuum of states over the barriers. On the contrary, if current is imposed, the injection field is set to a fixed value. The field discontinuity will simply lower the downstream field, and the whole 'S' pattern will be observed.

In reference 15, a qualitative interpretation of I(v) measurements shows that impact ionization is responsi-ble for the built-in charge , and that this process happens in the first wells of the structure. Two main microscopic pieces of evidence advocate for this statement. The wells in this structure are in under-populated regimes. There are fewer electrons than positive donors in the wells, which means the downstream electric field is lower than the upstream one. Hence the first wells of the structure sustain the highest fields (i.e. energy drops), and impact ionization is more likely to occur there. Moreover from the current level where the transition stands, an energy drop above the first period of the structure can be calculated: it is comparable to ∆E. Therefore the first wells after the contact are those involved in the regime transition. However a more quantitative theory of the built-in charge effects is needed to validate this hypothesis.

III. MODEL

The model presented here is inspired from the usual photoemissive models in QWIPs [START_REF] Ershov | [END_REF][18][19][20] . It is based on the fact that a macroscopic V (i) curve carries all the information about a phenomenon happening above a few wells. In order to get a better insight into the physics at stake, we translated this information to a microscopic scale. Two equations can describe the behavior of one well. The field discontinuity D between upstream and downstream fields is linked to the well population by the Poisson equation (1). The conservation equation (2) states that in steady state the current of particles J in that flows into the well equals J out that flows out of it:

D = F down -F up = -e(1 -n)N dop ε ( 1
)
dn dt = 0 = J in (F up , F down , n) -J out (F up , F down , n) (2)
N dop is the number of Si donors in a well, e the elementary charge, the dielectric permittivity, F up and F down the upstream and downstream electric fields, and n = N/N dop the relative well population (n = 1 when the well is neutral). These three variables D, F up and n linked by two equations can be merged into one law: D = f (F up ). This discontinuity law will be the object of interest in the rest of this paper. It is indeed a direct image of the built-in charge in the ionized QW.

For a given current I flowing through the structure, it is possible to calculate the injection field F 0 (I) as described in section IV. The D = f (F up ) law iteratively applied enables to calculate all the fields F k along the structure, and their sum gives the bias value. Let g be the unknown function that gives F down = g(F up ) (g(x) = f (x) + x), thus F k = g (k) (F 0 ). We solve the functional equation:

V exp (I) = L b × 40 k=0 g (k) (F 0 (I)) (3)
It is impossible to solve (3) analytically because V exp (I) is extracted from experimental data. Therefore section V will be dedicated to the discontinuity law D = f (F up ) through the analysis of its impact on the V exp (I) curve. In section VI we will present and discuss the solution obtained by numerical optimization.

IV. THE CONTACT

As shown in equation ( 2), an accurate model linking the current and the injection field is required. In THz QWIPs, the transition energy is around 10 meV, a value comparable to the corrections introduced by many-body interactions. Guo et al. showed that photoresponse peak positions cannot be explained without including their contribution 21 . Therefore it is critical to develop a good model of the band structure close to the contact. Another difficulty arises from the nature of the problem: while the contact is described by a 3D model, the quantum confinement in the wells is treated with a 2D one. This section is organized as follows. First we will determine the band structure at zero bias, by assuming that all the band structure rearrangements happen upon the first well. The criterion is that no total electric charge shall be present over the whole contact/first-well area. This will give the value F 0 (I)| 0V of the injection field when no current flows through the structure. Subsequently we shall calculate the band structure and the current for different values of F 0 , and use these points to interpolate the whole F 0 (I) law.

The effective mass m * (z) and parabolic band approximations lead to the following Schrödinger equation for the electron:

- 2 2 ∂ ∂z 1 m * (z) ∂ ∂z +V (z) ϕ l (z) = l ϕ l (z) (4)
where the total electron wave function is written as a product of a plane wave parallel to the layers, and a one-dimensional envelope function in the epitaxy growth direction z: ψ l,k / / r / / , z = e ik / / .r / / ϕ l (z). The energy of the electron is E l,k / / = 2 k 2 / / /2m * + l , where l is the subband index. V (z) is the bottom of the conduction band. The difference between 3D and 2D models lies in the way the electron density n(z) is calculated. Assuming a Fermi-Dirac distribution, we write for the contact:

n 3D (z) = 1 2π 2m * 2 3/2 V (z) ∞ E -V (z)dE 1 + exp E-E F k B T (5)
with T the temperature and E F the Fermi level. Beyond the middle of the first barrier, the density is considered two-dimensional:

n 2D (z) = m * k B T π 2 log 1 + exp E F -l k B T |ϕ l (z)| 2 (6)
Three contributions are taken for the potential term:

V (z) = V QW (z) + V H (z) + V XC (z). V QW (z) repre-
sents the potential of the quantum wells. V H (z) is the Hartree potential that accounts for the Coulomb interaction between electrons and ionized donors, and is calculated thanks to the Poisson equation:

∂ 2 ∂z 2 V H (z) = e 2 ε (N d (z) -n (z)) ( 7
)
where ε is the dielectric constant and N d (z) the doping density. V XC (z) accounts for the fermionic nature of the electrons, and the many-body wave function should be written as a Slater determinant: theoretically, it cannot be factorized in one-electron wave functions. We use the local density approximation (LDA) 22 to express it:

V XC (z) = e 2 4π 2 εa B (z) r s (z) 9π 4 1/3 × 1 + 0.0545r s (z) log 1 + 11.4 r s (z) (8)
with a B = ε 2 /e 2 m * (z) is the effective Bohr radius and r s (z) = (3/4π)(a 3 B n(z)) -1 1/3 . The charge neutrality condition gives the last equation of this close set (eq 4-9):

L 0 (N d (z) -n (z)) = 0 (9)
Because of the mixed-dimensionality of the problem, a numerical method was chosen to solve it and is presented in appendix. The main physical approximation is that we do not take into account the triangular quantum confinement at the end of the contact. This condition limits the validity of this model to small electric fields, relevant in the THz range. The band structure of the contact and the first period at zero bias is shown in figure 2. Figure 3 highlights the impact of many-body effects in THz QWIPs, as previously shown in 2D calculations 21 .

We use a similar method to determine the band structure for every given value of the injection field F 0 . The current is given by:

I (F 0 ) = S em * k B T 2π 3 ∞ 0 D (E, F 0 ) log 1 + e E F -E k B T dE (10
) with S the mesa area and D(E, F 0 ) the transparency of the injection barrier calculated with a transfer-matrix method 23 :

D (E, F 0 ) -1 = k in 4k out T 11 + k out k in T 22 2 + k out m * T 12 - m * k in T 21 2 (11)
with k in and k out the wavevectors on each side of the injection barrier, and the transfer matrix coefficients are Because of the triangular quantum confinement, we assume that the only electrons contributing to the current are the ones whose energies are above 0, for whom the 3D model is justified. For the lower electrons, their contribution to the current is exponentially decreasing with the height of the barrier, thus is expected to be small. This will lead to a slight under-evaluation of the current. For the sake of simplicity, we also only considered the current flowing from the contact to the first well, neglecting the current going back upstream. The latter is indeed negligible as soon as the effective Fermi levels of the contact and first well are no longer aligned (F 0 > 8 kV.cm -1 ), which is the range of interest. Figure 4 plots the current versus injection field for three different models: without any correction, with Hartree potential, and with both Hartree and exchangecorrelation contributions. Compared to the zeroth order calculation, the Hartree correction lowers the height of the barrier because of the electric field continuity condition. The current is then increased up to one order of magnitude (from 0.16 mA to 2.17 A at F 0 = 20 kV.cm -1 ). Moreover the 7.6 meV rise in the height of the barrier due to exchange-correlation effects makes the current to drop by one or two orders of magnitude. This clearly demonstrates how critical it is to consider those corrections in order to have a quantitative interpretation of such electronic transport regimes in THz QWIPs.

The inverse function of the high-field part (≥ 8 kV.cm -1 ) of the circle-marked curve on figure 4 gives the F 0 (I) law that will be implemented in the model exposed in section II.

V. ANALYZE OF THE DISCONTINUITY LAW

This section is dedicated to giving a practical insight into the shape of the discontinuity law. A sample D law is implemented in the model in order to calculate the corresponding V (i) curve. The macroscopic modifications induced by the modification of microscopic parameters are then analyzed. Figures 5 to 8 present different D laws and the corresponding current variations induced in the structure. The aim is to highlight the impact of each feature of the shape of the D law on the V (i) curve. This work is fundamental to parameterize the upcoming optimization in section V. The 'up' and 'down' regimes are separated by the threshold field F c . Below it, no electrons are present above the barriers, so the well is almost full, and the discontinuity is very small. For higher values of the upstream field, the well starts to be drained. If no capture mechanism is taken into account, the well will reach either of the following limit values. For low fields, the discontinuity is capped at a zero downstream field, in order for electrons to be extracted efficiently. For higher fields, the discontinuity cannot exceed the value given by a built-in charge equal to the doping amount in the well. These two values are labeled respectively F down = 0 and D max in figures 5 to 8. We will not consider the doping value as a varying parameter because its nominal amount has explained response peak in previous work 21 .

Figure 5 shows that the current at the transition I c goes up with the value of the threshold field F c . The injection field value controls both the current flowing through the structure and the amount of energy acquired by the electrons of the continuum over the first period. It shall be noticed that the usual threshold fields are far below the value where the F down = 0 and D max straight lines cross. Figure 6 shows that even in the down regime a positive charge in the wells has to be considered. Indeed, the bias at the transition V c decreases when this 'down' ionization increases. If the D law is capped at zero, the structure will be rigid and strong bias modulations will occur. On the contrary at higher ionization levels the structure will look bent like a chain of pearls, and bias effects will be attenuated. At a given F c this parameter directly sets the value of V c .

Another interesting feature is the height of the transition, or more precisely the evolution of the relative distance between the D law and the F down = 0 curves across and after the transition. Figure 7 proves that the 'S' pattern on the V (i) widens when the 'up' ionization gets more effective. As a matter of fact if the relative distance shrinks to zero, the bias on the structure will only be the one applied upstream of the ionizing well. On figure 1, the 'S' is thin (0.18 V wide), which can be explained in two ways. Either a well at the end of the structure is strongly drained, or a well at the beginning of the structure is lightly drained. We already mentioned that ionization occurs on one of the first wells of the structure. This clearly advocates for a partial ionization of the wells.

VI. OPTIMIZATION AND DISCUSSION

The arguments of section IV prove that the D law features have quite uncorrelated effects on the V (i) curve. This allows for a partially sequential resolution, much simpler than a raw multi-dimensional optimization, when trying to solve the functional equation ( 3).The fit of the experimental V (i) in order to output an optimized D law uses the quasi-Newton method BFGS scheme [START_REF] Press | Numerical recipes: the art of scientific computing[END_REF] . Firstly the threshold field F c is found by adjusting the transition current. Secondly the 'down' regime is calculated to fit the transition bias V c . The third calculated parameter is the height of the transition, given by the width of the 'S' pattern. To finish with, the shape of the 'S' grants access to the rest of the 'up' regime shape of the D law. The main result given by this first model is that the well entering its regime transition is only partially ionized. As will be shown in figure 9 even at the highest field, the wells are barely half drained. A refill mechanism attributed to capture by LO phonon emission has to balance the draining process.

A model with a single D law for all the wells (as considered in section IV) cannot explain the V (i) curve, and at the same time satisfy the impact ionization hypothesis. In this model the first well should be ionized first. The injection field should account for the amount of current at the transition (4.2 × 10 -6 A), which gives F 0 = 9.2 kV.cm -1 . However it should also grant enough energy to electrons over the first period for impact ionization, which gives F 0 ≈ 11 kV.cm -1 . This means that it is not the first well of the structure that gets ionized at this transition.

Another fact advocates for the ionization of several wells: the second cusp on the V (i) curve of figure 1, at 0.72 V and 9 × 10 -5 A. This cusp is the signature of another well entering its ionized regime. The band structure calculations of section III highlight the differences between the D laws of the first and the other wells. Figure 2 shows that at zero bias (a single Fermi level throughout the whole structure), there is already a non zero injection field (F 0 (I)| 0V = 2.061 kV.cm -1 ). Due to quantum confinement, the first well will be more ionized than the other wells in the down regime, and its threshold field will be higher.

D well 1 down = D other wells down -F 0 (I)| 0V F well 1 c = F other wells c + F 0 (I)| 0V + F ef f
where F 0 (I)| 0V = 2.061 kV.cm -1 is the value of the injection field calculated at zero bias, and F ef f stands for the fact that all the electrons that can potentially ionize the first well are injected from the contact quasi Fermi level. On the other hand the remaining wells see not only electrons injected at the upstream well quasi Fermi level, but also electrons already flowing through the continuum over the barriers. We complete our model by considering two different D laws: one for the first well after the contact, and one for all the other wells.

Figure 9 shows the optimized D laws and the comparison between experimental and calculated V (i). The latter shows an extremely good agreement, which should not be given too much credit: equation 3 has an infinite number of solutions. The sensitivity of the calculated V (i) towards variations of the D laws is very acute next to the regimes-switch zone. Consequently the relevance of the adjusted D laws reaches its maximum around F c . The good fit means that the objective of translating the information from the V (i) to the D law curves has been fulfilled.

It shall be noticed that both the D laws of the first and other wells have the same shape, modulo the D down and F c translations. This renforcer the hypothesis that the process at stake is local. This is consistent with the same mechanism being responsible for the 'up' regime part in both cases: balance between well ionization and electron capture.

Note that the calculation of the output D laws is purely numerical, and does not rely on any physical assumption about the mechanism building up the charge. However the returned threshold field values exactly match the energy condition for intersubband impact ionization. For the first well, F c = 10.94 ± 0.06 kV.cm -1 The ionizing electrons injected at the contact Fermi energy arrive with an energy of E c = 68 meV above the bound level. The energy condition for ionization E c > 2∆E = 64.5 meV is satisfied. For the second well at the transition, the threshold field value is F c = 6.83 ± 0.03 kV.cm -1 . In the hypothesis of a ballistic transport from the contact, this gives an incident energy of 103 meV above the bound level. Considering one LO phonon emission, the same energy condition E c = 67 meV arises. It is thus very likely that the mechanism responsible for the built-in charge effects that the regime transition in THz QWIPs is intersubband impact ionization. Further work will be published on these considerations about hot electrons kinetics.

In this sample, the second well of the structure will be ionized before the first one. In the latter case no NDR regime is observed. This counter-intuitive order of regime transitions highlights the critical importance of having a precise model of the contact. In Fig. 9, there is a small discrepancy in the shapes of the two D laws at the beginning of the 'up' regime part. When the first well begins to be ionized, the effective energy of the electrons arriving above the second well is increased, which downshifts the threshold field F c ) of the second well. In order to avoid complicating further the model by introducing another D law for the second well, we chose to add a small correction on the D law of the first well to keep a good fit of the V (i).

Focusing back to the I(v), it is possible to understand the asymmetry of the hysteresis pattern. When sweeping bias upwards, the breakdown of the second barrier drags along the first well into ionization. The third well is expected to be ionized at higher bias, out of the measurement range. When sweeping down the bias, the electric fields configuration allows the first well of the structure to shift to the 'down' regime without dragging along the second one.

VII. CONCLUSION

We have shown that in THz QWIPs, the apparition of a built-in charge in the quantum wells generates drastic band structure rearrangements, which have been thoroughly investigated and successfully explained a regime transition in electronic transport. In particular the first wells of the structure are proven to be only partially drained at high electric fields. In the studied sample the second well enters its ionized regime before the first one. A precise 3D-2D model of the contact considering many-body effects highlights the difference in behavior between the first well and the other ones in the structure. Intersubband impact ionization has been identified as the probable origin of the built-in charge. Calculations of the D laws down from a quantum description are expected to give a good insight into the electron distribution in the continuum above the barriers, and thus raises interest in wavelengths broader than the THz spectrum.

Using the regime switch caused by impact ionization allows to optimize the design of very high sensitivity fast THz detectors. When the structure is polarized just below the threshold bias V c , a few photons absorbed in the second well can trigger the shift to the 'up' regime. The main obstacle in the way of effective detectors is the fluctuations area near the transition, wherein the regime shift is not controllable. Understanding and deleting this area is the focus point for future work. Noise surge next to the instability is expected to be an issue, but the huge gain value (the current shifts by five orders of magnitude) soothes this concern. APPENDIX This appendix presents the numerical method used to solve the set of equations 5-9. The wave function of the well ground level is first calculated with a shooting method [START_REF] Harrison | Quantum wells, wires and dots: theoretical and computational physics of semiconductor nanostructures[END_REF] . If V XC is set to zero, equations 5-7 and 9 can be fused in the charge condition Θ(E F , V H ) = 0. Θ is a decreasing function of E F , so a simple way to solve this equation is to find two E F values where Θ is positive and negative, and use dichotomy. The key for the success of the method is the nature of equation 7 that allows to calculate the correction V H (z i ) thanks to the density at the previous step n(z i-1 ). This is not possible while introducing the exchange correlation term. After getting a band structure V (1) (z) obeying to the Poisson equation (see figure 3), equation 8 gives the V

(1) XC correction. It is then input along V (1) (z) to recalculate the wave function and resolve Θ = 0 to get V (2) (z). The result will be considered satisfactory as long as the change of the band structure in the second step is not significant:

V (2) (z) ≈ V (1) (z) + V (1)
XC , which is always the case. A similar method is used to calculate the band structure for any given value of F 0 . The new criterion for the Θ function is the value of the injection field. The 2D part of the problem is simplified: no calculations are made for the well area, except that its contribution to V XC at zero bias is added to the 3D exchange-correlation potential. This slightly modifies the shape of the end of the barrier, and has a noticeable effect only at low fields when the barrier is not yet triangular.

FIG. 1 :

 1 FIG.1: I(v) (dashed) and V (i) (plain) dark measurements for positive bias at T = 4 K, for upwards and downwards sweeps. Top: whole curves. Bottom: zoom onto the hysteresis pattern. The current discontinuity reaches here three orders of magnitude, from 2.4 × 10 -6 to 4.7 × 10 -3 A.

FIG. 2 : 1 FIG. 3 :

 213 FIG. 2: Band diagram of the contact and first period of v265 calculated including Hartree and exchange-correlation corrections. At zero bias the injection field is 2.061 kV.cm -1

FIG. 4 :

 4 FIG.4: I(F0) curves with three different levels of approximations. The cusps around 6 -7 kV.cm -1 appear when the barrier becomes triangular.

FIG. 5 :

 5 FIG. 5: Top: sample D laws with varying threshold fields Fc. The ionization degree of the well is given on the right axis. Bottom: corresponding V (i) curves. The value of Fc sets the level of the transition current.

FIG. 6 :

 6 FIG.6: Top: sample D laws with varying 'down' regime behavior. Bottom: corresponding V (i) curves. The value of the threshold bias decreases when the wells are more ionized in the 'down' regime.

FIG. 7 :

 7 FIG. 7: Top: sample D laws with varying amplitudes of the drain at the transition. Bottom: corresponding V (i) curves. A wide 'S' will follow an effective ionization of the wells.

Figure 8

 8 Figure 8 highlights the impact of the last critical parameter of the D law: the local slope. When it is inferior to -1 in the 'up' regime, an increase of the upstream field will cause a decrease of the downstream field. The structure is then in the negative differential resistance (NDR) regime. The rest of the 'S' pattern is shaped by the steepness of the D law.

FIG. 8 :

 8 FIG. 8: Top: sample D laws with varying negative differential resistance thresholds. These are indicated by the -1 slope tangents. Bottom: corresponding V (i) curves.

FIG. 9 :

 9 FIG. 9: Top: Adjusted D laws with the first well singled out. Bottom: Experimental and adjusted V (i) curves.


François-Régis Jasnot

Nicolas Péré-Laperne

Louis-Anne De Vaulchier

Yves Guldner

Francesca Carosella

Robson Ferreira

Amandine Buffaz

Laetitia Doyennette

Vincent Berger

Mathieu Carras

Xavier Marcadet

Magnetotransport in quantum cascade detectors : analyzing the current under illumination

Photocurrent measurements have been performed on a quantum cascade detector structure under strong magnetic field applied parallel to the growth axis. The photocurrent shows oscillations as a function of B. In order to describe that behavior, we have developed a rate equation model. The interpretation of the experimental data supports the idea that an elastic scattering contribution plays a central role in the behavior of those structures. We present a calculation of electron lifetime versus magnetic field which suggests that impurities scattering in the active region is the limiting factor. These experiments lead to a better understanding of these complex structures and give key parameters to optimize them further.

The quantum cascade detector [1] (QCD) recently proposed and realized in both the mid-infrared [2] and in the THz [3,4] range is a photovoltaic version of the quantum well infrared photodetector (QWIP) [5]. Their band structure are designed as quantum cascade laser (QCL) without any applied bias voltage [1,3]. QCD are totally passive systems and show a response only to photon excitation. As such, the QCD structure is designed to generate an electronic displacement under illumination through a cascade of quantum levels without the need of an applied bias voltage.

In a semiconductor quantum well structure, magnetic field applied along the growth direction breaks the 2D in-plane continuum into discrete Landau levels (LLs). This experimental technique has been used to evaluate the different contributions of scattering mechanism in complex quantum cascade structures [5][6][7][8][9].

We present in this paper experimental photocurrent measurements under magnetic field applied along growth direction. We develop a simple model of transport under illumination in a QCD. Through a comparison between experimental and calculation results, we evidence the mechanism limiting the response of the QCD.

The QCD under study is a GaAs/Al 0.34 Ga 0.66 As heterostructure with a detection wavelength of 8 µm as described in ref. [9]. It consists of 40 identical periods of 7 coupled GaAs quantum wells. Figure 1 recalls the principle of the device.

QCDs are mounted inside an insert at the centre of a superconducting coil where a magnetic field B up to 16 T can be applied parallel to the growth axis. Light is emitted by a globar source from a FTIR spectrometer and guided to the sample. The experiment consists in measuring the current under illumination (I light ) without any applied voltage at 80 K while the magnetic field is swept from zero to 16 T.

Experimental result is illustrated on figure 2(a). The photocurrent shows oscillations as a function of the magnetic field, superimposed on a continous decreasing background which is removed of the experimental data in figure 2 We propose a model of transport in one period based on a rate equation approach. We assume that electrons are in the upper detector state |up through absorption of a photon. Current as a function of lifetimes involved in this structure can be written :

J q = αN down τ up-down τ up-down + τ up-c = αN down QE . (1)
The parameters α and N down are respectively the absorption factor and sheet density of |down and are constant. The subscribe c stands for the whole cascade. We present in table I the calculated scattering rates of the different processes at B = 0 T. For interface roughness, we used a Gaussian autocorrelation of the roughness, with an TABLE I: Scattering rates in s -1 are calculated using different scattering processes for an electron in the |up subband at B = 0 T. average height of ∆ = 2.8 Å and a correlation length of Λ = 60 Å. LO phonon emission scattering rate has been calculated as in ref. [10]. In our structure impurities scattering is the most efficient process [11]. Usually in GaAs quantum cascade structures this mechanism is neglected because the doped layers are not in the active region. In order to take into account the main scattering process we calculate ionized-impurities scattering as a function of magnetic field. The details of the calculation are presented elsewhere [12].

Figure 3 represents a comparison between experimental data and electron-ionized impurities scattering time as a function of magnetic field. Figure 3(b) and 3(c) show the two lifetimes involved in Eq. 1 as a function of B calculated with electron-ionized impurities scattering. Figure 3(d) shows the calculation of the related quantum efficiency.

The oscillating behavior at high magnetic field (B > 9 T) is a result of the electronic transfer from |up to |down . This transfer leads to minima in the current which fit well with τ imp up-down and QE. The long period oscillating behavior of τ imp up-c as a function of B enhances the peak at B = 14 T in QE in agreement with experimental data. QE, which describes the performance of the detector, is oscillating between 74% and 85% under B. By extrapolating, at B = 0 T, QE is equal to 75%, a value that should be increased to improve the detector performance. An optimized structure should take these results into account by shifting the ionized impurities from the active region, where they are enhancing τ imp up-down , to a position where they would only enhance τ imp up-c . The series of peak at B < 9 T corresponds to a characteristic energy of 40 meV. This energy is attributed to transitions in the cascade involving an elastic scattering mechanism.

In conclusion, we observe oscillations of the photocurrent in a mid infrared QCD as a function of B. These socillations are due to electron-ionized impurities scattering. This mechanism is dominant in this structure because impurities are located in the active region. In order to improve further this efficiency, we suggest to shift the impurities in another location of the structure in order to minimize τ imp up-down . The Laboratoire Pierre Aigrain is a " Unité Mixte de Recherche " between École Normale Supérieure, the CNRS, the University Paris 6 and the University Paris 7.

This work has been supported by a grant of the Agence Nationale pour la Recherche (ANR).

  (b). Minima of current are located at B = 10.1 T, 11.4 T, 13.0 T and 15.3 T

FIGURE 1 :

 1 FIGURE 1: Conduction band diagram of one period of an 8 µm QCD showing the energy levels. Note that the ground state of the first QW belongs to the former period and is noted |down . The arrows illustrate the electronic path during a detection event. The layer sequence is as follows 67.8 / 56.5 / 19.8 / 39.6 / 22.6 / 31.1 / 28.3 / 31.1 / 33.9 / 31.1 / 39.6 / 31.1 / 45.2 / 50.8 (the barriers are represented in bold types). The n-doping of the large QW is 5 × 10 11 cm -2 .

FIGURE 2 :

 2 FIGURE 2: (a) Current under illumination as a function of B at 80K and at zero bias. (b) I light as a function of B where the decreasing backgroung as been substracted. (c) Fan chart of |up, 0 and |down, p as a function of B taking into account the band nonparabolicity.

FIGURE 3 :

 3 FIGURE 3: (a) I light as a function of the magnetic field where the backgroung has been substracted. (b) Ionized impurity scattering τ imp up-down under magnetic field between |up and |down levels. (c) Ionized impurity sattering τ imp up-c under magnetic field between |up and levels in the cascade. (d) Quantum efficiency (QE) calculated with Eq. 1.


Ahmed Khelassi
email: ahmed.khelassi@cran.uhp-nancy.fr

Didier Theilliol
email: didier.theilliol@cran.uhp-nancy.fr

Philippe Weber
email: philippe.weber@cran.uhp-nancy.fr

Jean-Christophe Ponsart

Fault-tolerant control design with respect to actuator health degradation: An LMI approach

The active fault-tolerant control approach relies heavily on the occurred faults. In order to improve the safety of the reconfigurable system, a methodology to incorporate actuator health in the fault-tolerant control design is proposed for a tracking problem. Indeed, information about actuator health degradation due to the applied control is considered in addition to fault estimation. The main objective is to design a fault-tolerant control system which guarantees a high overall system reliability and dependability both in nominal operation and in the presence of faults. Such an objective is achieved by a control performance index, which is proposed based on system reliability analysis. The fault-tolerant controller is synthesized by using a linear matrix inequality approach.

I. INTRODUCTION

The safety-critical nature of embedded systems used in air-craft, space, tactical, and automotive applications, mandates that the functionality for which they are designed be performed even in the presence of component faults. Thus, safety-critical systems should be able to adapt and compensate the component faults in systematic way. These types of adaptable systems are known as fault-tolerant control systems (FTCS). The aim of FTCS is to keep plant available by the ability to achieve the objectives assigned in the faulty situations and to accept reduced performance when critical faults occur [START_REF] Blanke | Diagnosis and Fault Tolerant Control[END_REF].

Various approaches for FTCS design have been proposed in the literature. Overviews on the development of FTCS have been provided in survey papers and books as for example in [START_REF] Zhang | Bibliographical review on reconfigurable tolerant-control system[END_REF] and [START_REF] Noura | Fault Tolerant Control Systems: Design and Practical Application[END_REF]. The FTCS design approaches can be classified under two types, the passive approaches and the active approaches [START_REF] Zhang | Bibliographical review on reconfigurable tolerant-control system[END_REF]. In the passive approaches, the same controller is used for the normal case as well as for the faulty cases where, a presumed set of process component faults are considered in the design stage. However, active FTC requires an on-line Fault Detection and Diagnosis (FDD) process and control reconfiguration mechanism [START_REF] Zhang | Bibliographical review on reconfigurable tolerant-control system[END_REF].

Therefore, it is important to enhance the system safety not only by improving reliability of individual components, or by designing control systems to compensate the effect of faults but also, by taking into consideration the degradation of components health in the controller design [START_REF] Khelassi | Reconfiguration of control inputs for over-actuated systems based on actuators health[END_REF]. Indeed, as illustrated in [START_REF] Blanke | Diagnosis and Fault Tolerant Control[END_REF], fault-tolerant control can improve the overall system reliability, because with an effective control input management, the overall system remains operational for a longer duration with an acceptable performance requirement. This problem is more considerable in practical safety-critical systems and in flight control systems where a catastrophic consequences from degraded functional operation can occur [START_REF] Zhang | Bibliographical review on reconfigurable tolerant-control system[END_REF].

Some research works have introduced reliability analysis for FTCS. [START_REF] Wu | Reliability and fault tolerant control system Part II[END_REF] and [START_REF] Wu | Reliability and supervisory control[END_REF] have used Markov models to dictate the system reliability where subsystems are supposed to reach two states: intact (available) of failed (unavailable) states. [START_REF] Staroswiecki | Sensor network design for fault tolerant estimation[END_REF] and [START_REF] Staroswiecki | General fault tolerant control design[END_REF] the authors proposed a sensor and actuator reconfiguration strategy based on physical redundancy. The reliability analysis provides some informations to select the optimal reconfigurable set. In a similar way, [START_REF] He | Robust H∞ filtering for timedelay systems with probabilistic sensor faults[END_REF] has considered the reliability of sensor faults in the filtering design issue. A reconfiguration mechanism of FTC strategy incorporating reliability analysis under a dynamic behavior constraints has been proposed in [START_REF] Guenab | Design of a fault tolerant control system incorporating reliability analysis under dynamic behaviour constraints[END_REF] .

In this paper, informations about actuator health are integrated in the fault-tolerant control design. The proposed design contribute to a fault-tolerant controller which achieves the control objective with a high overall system reliability. Indeed, such design improves the dependability of the system and keeps the set of the actuators available as long as possible by minimizing the use of the more sensitive actuators. The degradation of actuator health is evaluated by estimating the impact of the required load on the actuator life time. In this context, the solution of fault-tolerant tracking problem incorporating actuator health is proposed by using an LMI approach.

This paper is organized as follows. In section II, the faulty system model and the fault-tolerant tracking control problem are presented. In section III, a performance index for control design is proposed with respect to actuator health level. Section IV, a solution of fault-tolerant control tracking problem with respect to overall system reliability is proposed. In section V, a numerical example of an F-16 aircraft model and its simulation results are given and compared to a classical control performance index. Some concluding remarks are given in section VI.

II. PRELIMINARIES AND PROBLEM FORMULATION

Let us consider the linear time-invariant model given by:

ẋ(t) = Ax(t) + Bu(t) y(t) = Cx(t) z(t) = Sy(t) (1)
where A ∈ IR n×n , B ∈ IR n×m and C ∈ IR p×n are respectively, the state, the control and the output matrices.

x ∈ IR n is the system state, u ∈ IR m is the control input, y ∈ IR p is the system output. S ∈ R l×p is a known constant matrix used to select the output required to track the reference signal. z(t) ∈ IR l is the controlled output and (A, B) is stabilizable.

In this study, the considered faults are loss of control effectiveness where the system (1) can be written in degraded operating mode as follows:

ẋ(t) = Ax(t) + B f u(t) y(t) = Cx(t) z(t) = Sy(t) (2)
The matrix B f is written according to the nominal control input matrix B and the control effectiveness factors ω i ∈ [0 1], i = 1, . . . , m, as follows:

B f = Bω, ω = diag([ω 1 , ω 2 , . . . , ω m ]) (3)
Indeed, if ω i = 1, then the i th actuator is considered to be fault-free. Nevertheless, when 0 < ω i ≤ 1, the considered fault is a partial loss of control effectiveness. Moreover, when ω i = 0 failure is considered and the actuator is out of order. For the faulty system described by ( 2), the fault-tolerant control design problem is to determine a controller such that: i) During normal operating mode, the closed-loop system is stable and the output z(t) tracks the reference signal y r (t) without steady-state error:

lim t-→∞ e(t) = 0, e(t) = y r (t) -Sy(t) (4)
In fact, this study deals with active fault-tolerant control AFTC. ii) After faults diagnosis, the closed-loop system is still stable and the required output z(t) tracks the reference signal y r (t) with an acceptable lower level of tracking performance. Is is known that integrator control can be added in the closed-loop system to eliminate the steady state tracking error. The integral term can be constantly calculating the value of the control required at the set-point to force the error to go to zero. To accomplish the design of the control system for the integral and the original state vector, an augmented model is considered by combining (2) and (4) as follows:

ε(t) ẋ(t) = 0 -SC 0 A ε(t) x(t) + 0 B ωu(t) + I 0 y r (t) (5)
where ε(t) = t 0 e(τ )dτ is the integral error state. By considering x(t) = [ε T (t) x T (t)] T , the augmented system (5) can be described as follows:

ẋ(t) = Āx(t) + Bωu(t) + Gy r (t) (6)
where

Ā = 0 -SC 0 A , B = [ 0 B ] T , G = [ I 0 ] T .
Indeed, the performance index which the designed controller should minimize is written as follows:

J t = t 0 [ε(t) T Q 1 ε(t) + x(t) T Q 2 x(t) + u(t) T ωRωu(t)]dt (7) where Q 1 ∈ IR l×l , Q 2 ∈ IR n×n are semidefinite matrices and R ∈ IR m×m is a positive definite matrix.
Several fault-tolerant control approaches have been proposed in the literature but without interest in the improvement of the reconfigurable system life time. The main contribution of this work is to consider the overall system reliability as a principle objective which guides the design of the fault-tolerant controller. Thus, information about actuator health can be introduced in the performance index [START_REF] Liao | Reliable robust flight tracking control: An LMI approach[END_REF] in order to guarantee that the control efforts are distributed to the system by taking into account the reliability of the actuators. The aim is design a fault-tolerant controller that maximizes the overall system reliability. The stability and tracking problem will be solved, both in the nominal case and in the faulty conditions, by using an LMI approach.

III. CONTROL PERFORMANCE INDEX BASED ON RELIABILITY ANALYSIS

One motivation to integrate information about actuator health in the controller design is to improve the life time of the reconfigurable systems which is quantified by the overall system reliability estimation. Let us consider the following definitions:

Definition 1: Actuator health is a measure that can be affected by several parameters and factors. In this study, actuator reliability and loss of control effectiveness are both considered.

Definition 2: Reliability R(t) is defined as the probability that units, components, equipments and systems will accomplish their intended function for a specified period of time under some operating conditions and specific environments, [START_REF] Gertsbakh | Reliability Theory with Applications to Preventive Maintenance[END_REF].

In the useful period of life, the component can be characterized at a given time t by a baseline reliability measure R 0 (t). In following, R 0 i (t) will be assigned to the reliability of the i th actuator obtained under nominal operating conditions in the useful period of life such as:

R 0 i (t) = exp(-λ 0 i t), i = 1 . . . m (8)
where in the reliability engineering field, λ 0 i is a baseline failure rate of the i th component defined under a nominal load level called stress level.

However, a realistic health measurement should also include the trend of actuator degradation according to the variation of the operating conditions. Indeed, in many situations and especially in the considered study, failure rates are obtained from actuators under different levels of load depending on the applied control input. Several mathematical models have been developed to introduce the impact of the variable load level in the reliability estimation [START_REF] Gertsbakh | Reliability Theory with Applications to Preventive Maintenance[END_REF]. Proportional hazard model firstly proposed in [START_REF] Cox | Regression models and life tables[END_REF] is adopted in this paper.

Definition 3: For a variable functional condition, the actuator reliability can be estimated according to the baseline failure rate as follows:

R i (t) = exp(-λ i t) λ i = λ 0 i × g i (ℓ, ϑ) (9)
where λ 0 i represents the baseline failure rate (nominal failure rate) for the i th actuator and g i (ℓ, ϑ) is a load function (independent of time) which represents the effects of stress on the component failure rate ℓ represents an image of the applied load and ϑ is a component parameter.

Assumption 1: In this study, the exponential distribution for reliability estimation is considered. The load function (9) will be defined according to the load applied on the actuator until the end of the mission t = t M . The component reliability measure can be evaluated for t = t M and noted R(t M ).

Different definitions of the load function g i (ℓ, ϑ) exist in the literature. However, the exponential form is commonly used [START_REF] Guenab | Design of a fault tolerant control system incorporating reliability analysis under dynamic behaviour constraints[END_REF]. In this study, the load function is defined according to the root-mean-square of the applied control input as follows:

g i (ℓ, ϑ) = exp(β i tM 0 u 2 i (t)dt), i = 1, . . . , m (10)
where β i is an actuator parameter defined as follows:

β i = (t M (u i -u i )) -1 , i = 1, . . . , m (11)
where u i = and u i are the upper and lower saturation bound of u i respectively. Property 1: Based on ( 9) and ( 10), the actuator health at a given time t M is modeled as a function of the applied load where the actuator reliability degradation increases for a large control input.

Indeed, the actuator reliability R i (t M ) can be estimated depending on the baseline reliability level R 0 i (t M ) as follows:

R i (t M ) = α i (u i )R 0 i (t M ), i = 1, . . . m α i (u i ) = exp(-λ 0 i β i u i 2 ) (12)
where, . is the Euclidien norm. α i (u i ) ≤ 1 represents the rate of reliability degradation due to the applied load during the mission. In fact, the following relation is satisfied:

if u i (t M ) 2 → 0 : R(t M ) → R 0 (t M ).
Lemma 1: For a system composed by m redundant actuators, the overall system reliability Rg(t M ) -→ Rg max for an effective control inputs u * (t) = -Kx(t) ∈ U that stress/solicit less the more reliable actuators in the control efforts distribution and satisfy the following condition as close as possible:

U = {u * (t) | min α i (u i ) α j (u j ) if λ 0 i λ 0 j < 1}
where K is the nominal controller, i ∈ [1, . . . , m], j ∈ [1, . . . , m] and i = j. Rg max is the optimal value of overall system reliability that can be obtained at t = t M . Proof: Systems composed by m redundant actuators can be presented by a parallel scheme of Reliability Block Diagram (RBD). The overall system reliability under nominal conditions Rg 0 (t M ) can be obtained as follows:

Rg 0 (t M ) = 1 - m k=1 (1 -R 0 k (t M
)) The sensitivity of the overall system reliability versus the i th actuator s i can be evaluated as:

s i = dRg 0 (tM ) dR 0 i (tM ) = m k=1; k =i (1 -R 0 k (t M )) = m k=1 (1 -R 0 k (t M ))(1 -R 0 i (t M )) -1 (13)
It can be shown that, if R 0 i > R 0 j (or λ 0 i < λ 0 j ) then s i < s j . This result means that a small degradation α i (u i ) of the more sensitive actuator i causes a large degradation of the overall system reliability Rg(t M ). So, the optimal value of the overall system reliability under a variable functional condition can be obtained for control inputs as follows:

u * = arg min αi(ui) αj (uj ) : Rg(t M ) -→ Rg max s.t. λ 0 i λ 0 j < 1
which limit the use of the less reliable actuators in the control efforts distribution. The actuators failure rates can be ordered from the less sensitive actuator to the more sensitive as follows:

l 0 = {λ 0 ms , λ 0 ms+1 , λ 0 ms+2 , . . . , λ 0 ps } (14)
where ps and ms are respectively the indexes of the most sensitive and the less sensitive actuator with ps := ms+m-1. The associated degradation rates α i (u i ) can be re-arranged also as follows:

Θ = {α ms (u), α ms+1 (u), α ms+2 (u) . . . , α ps (u)} (15)

where to simplify the notation, α i (u i ) := α i (u) is adopted.

In order to design an effective fault-tolerant controller which guarantees a high overall system reliability taking into consideration the actuator health, the control inputs delivered by the controller should minimize the degradation rate of the most sensitive actuators and maximize the less sensitive for each subset (α ms+k-1 (u), α ms+k (u)) as close as possible. This objective can be reformulated as follows:

arg min u m-1 k=1 ( α ms+k (u) α ms+k-1 (u) ), k = 1, . . . , (m -1) (16)
Based on (12) associated with Taylor approximation (α i (u) < 1), the optimization objective (16) can be written as follows,

arg min u m-1 k=1 ( u ms+k-1 2 - λ 0 ms+k β ms+k λ 0 ms+k-1 β ms+k-1 u ms+k 2 ) (17) = arg min u ( u ms 2 -ρ ms+1 u ms+1 2 )+ ( u ms+1 2 -ρ ms+2 u ms+2 2 ) + . . . . . . + ( u ms+m-2 2 -ρ ps u ps 2 )
which is equivalent to,

arg min u m i=1 r i u i 2 = arg min u m i=1 r i t 0 u 2 i (t)dt (18) where    r i = 1 i = ms r i = -ρ i i = ps r i = 1 -ρ i ps < i < ms (19) and ρ i = λ 0 i β i (λ 0 i-1 β i-1 ) -1 < 1, (20)
The objective function (18) can be written as:

arg min u t 0 (r 1 u 2 1 (t) + r 2 u 2 2 (t) + . . . + r m u 2 m (t))dt = arg min u t 0 u(t) T Ru(t)dt (21)
where, R = diag([r 1 , r 2 , . . . , r m ])

However, in order to reformulate the function objective (21) such ( 19) and (20), as an index control performance with respect to the actuator health, the following performance index is proposed:

arg min u t 0 u(t) T (R + σI m )u(t)dt = t 0 u(t) T Ru(t)dt (22)
where σ > 1 is introduced in the criterion to ensure the convexity of the optimization problem with R ≻ 0, ∀ σ > 1.

In fact, σ is an adapted parameter used in the design stage to choose the desired dynamic of the closed-loop system defined by A, B, C, Q and R with respect to the actuators health levels.

To design an affective fault-tolerant controller which maximizes the overall system reliability, the classical performance index ( 7) is reformulated as in the following proposition: Proposition 1: To design an effective fault-tolerant control system by taking into account the actuator health, the controller can be synthesized in order to minimize the upper bound of the following performance index:

J * t = t 0 [ε(t) T Q 1 ε(t) + x(t) T Q 2 x(t) + u(t) T ω Rωu(t)]dt (23
) where R ≻ 0 is a diagonal matrix designed based on actuators failure rates as in ( 19) and ( 22).

IV. ACTIVE FAULT-TOLERANT SCHEME FOR A RELIABLE DESIGN

A solution to design an effective active fault-tolerant controller which minimizes the performance index (23) (both in the nominal case ω = I m and after fault occurrence) is presented for the tracking problem. In fact, the active faulttolerant methods involve the introduction of an FDI module in the control loop. We assume that the fault is detected at t = t f and the reconfigurable control inputs can be applied at time τ = t f + ∆t.

Let us consider the augmented system ( 6) with the following state feedback tracking controller:

u(t) = K x(t) = [K e , K x ] ε(t) x(t) (24)
The closed-loop augmented system is given by:

ẋ(t) = ( Ā + BK)x(t) + Gy r (t) (25)
The effective fault-tolerant controller which holds the tracking problem and minimizes the proposed control index (23) can be obtained by the following lemma.

Lemma 2: Consider the closed-loop augmented system (25) and the performance index (23). The following controller stabilizes the closed-loop augmented system for a given positive constant γ:

K = [K e , K x ] = XZ -1
(26) such the linear matrix inequality:

      ĀZ + Bf X +( ĀZ + Bf X) T G X T ω R1/2 ZQ 1/2 * -γI 0 0 * * -I 0 * * * -I       ≺ 0 ( 27
)
where

Z ∈ IR (n+l)×(n+l) , X ∈ IR m×(n+l) , Q = diag[Q 1 , Q 2 ] ≥ 0 and R ≻ 0.
Moreover, the upper bound of performance index J * t can be minimized by solving the following optimization problem min Trace(T ) s.t. (27) and

T I I Z (28)
where T ∈ IR (n+l)×(n+l) . In the nominal case, the problem is solved for ω = I m and B f = B. After fault occurrence, the controller K is calculated for the considering faulty actuator ω i = 1.

Proof: For the nominal case, the proof can be found in [START_REF] Ye | Adaptive fault tolerant tracking control against actuator faults with application to flight control[END_REF] by considering ω I m and B f = B in the constraint (27). For the degraded functional mode and after faults occurrence, the result follows by considering the value of loss control effectiveness ω.

Indeed, in nominal functional mode, the weighting matrix R is proposed based on the baseline actuators failure rates λ 0 i . After faults occurrence, a new values of actuators failure rates can be estimated on-line taken into account the applied load as illustrated in [START_REF] Staroswiecki | Sensor network design for fault tolerant estimation[END_REF]. Then, a new value of the weighting matrix R is obtained and the fault-tolerant controller is calculated subject to the constraint (27).

V. FLIGHT CONTROL APPLICATION

In this section, an example of tracking control for a linear F-16 aircraft model is given to demonstrate the proposed approach. The linearized aircraft model is described as in [START_REF] Blanke | Diagnosis and Fault Tolerant Control[END_REF], where, x(t) = [u x , u z , q, u y , p, r a ] T is the state, u(t) = [δ hr , δ hl , δ ar , δ al , δ r ] T is the control input, and y(t) = [a, μ, r s , α, β] T is the measured output, respectively. u x , u y and u z are components of aircraft velocity along X a , Y a and Z a body axes respectively. p, q and r a are the roll rate about X a body axis, the pitch rate about the Y a body axis, and the yaw rate about Z a body axis, respectively. δ hr , δ hl , δ ar , δ al , δ r are right horizontal stabilator, left horizontal stabilator, right aileron, left aileron, and rudder, respectively. μ is the stability-axis roll rate and r a is the stability-axis yaw rate. α is angle of attack and β is angle of sideslip. The matrices, A, B and C are given in [START_REF] Liao | Reliable robust flight tracking control: An LMI approach[END_REF] (Example 1). For the degraded situations, each of the five actuators may lose its effectiveness. For the desired system dynamic and as chosen in [START_REF] Ye | Adaptive fault tolerant tracking control against actuator faults with application to flight control[END_REF], we consider γ = 2 with:

S =   0 1 0 0 0 0 0 0 1 0 0 0 0 0 1   Q 1 = diag [0.09, 1, 0.36], Q 2 = diag [0.
01, 0.2, 0.01, 0, 0, 0] where the matrix S determines the output required to track, i.e., μ, α and β.

To illustrate the proposed approach in the short time window, the values of the actuators failure rates are adapted with the time of the considered scenario. The failure rates are considered with a very huge values and given in Table 1.

First, in the nominal case, two nominal controllers are calculated to guarantee the tracking requirement by solving the LMI problem (27). A nominal controller is obtained subject to the classical performance index [START_REF] Liao | Reliable robust flight tracking control: An LMI approach[END_REF] where the weighting matrix R = diag([0.49, 0.49, 0.004, 0.04, 0.04]) is chosen as in [START_REF] Ye | Adaptive fault tolerant tracking control against actuator faults with application to flight control[END_REF]. Moreover, an effective nominal controller is calculated by considering the proposed performance index (23) where the matrix R = diag([0.10, 0.54, 0.35, 0.60, 1.30]) is obtained based on the actuators failure rates λ 0 i . The parameter σ = 1.1 is chosen in order to obtain a closedloop dynamic close to the classical one. In fact, the nominal The effective nominal controller is designed in order to manage smartly the set of the actuators by considering the actuator health level. The actuator rate limits are not taken into account in the design. The Figure [START_REF] Blanke | Diagnosis and Fault Tolerant Control[END_REF] shows the trajectory aircraft in the nominal case for the proposed performance index and for the classical one. It can be seen that the controlled outputs track the desired outputs without steadystate error. In Figure (2) the generated control inputs are presented where as can be seen clearly, the sensitive actuators are less solicited compared to the classical controller. In the same way, the required efforts are distributed to the less sensitive actuators in order to guarantee the stability and the tracking requirement. In fact, as presented in the Table 1, the set of the ordered actuators ( 14) is obtained from the less sensitive to the most sensitive one as: l 0 = {λ 0 5 , λ 0 4 , λ 0 3 , λ 0 2 , λ 0 1 }. As consequence, the effective controller minimizes the load applied to the actuators 5, 4 and 3. The required efforts are compensated by the solicitation of the actuators 1 and 2. For the overall system reliability, it can be seen that the reliability of the system with the proposed controller R * (t) is more than in the classical one R(t) where In fact, this result is evident where the sensitive actuators which cause an important reliability degradation of the overall system are less used.

For the faulty situation, we consider a failure of the actuator 3 at t = 7sec. This degraded mode can be characterized by ω 3 = 0 from the time of fault occurrence. The controller is calculated in order to find a novel distribution of the desired efforts that tracks the reference and stabilizes the system taken into consideration the impact of the applied load on the health level associated to each actuator. For the proposed index performance, the control efforts are distributed with respect to the actuators health where a new estimation of the failure rates λ i using ( 9) is considered in the performance index (23). The aim is to obtain a high overall system reliability by minimizing the use of the sensitive actuators taking into account the health degradation caused by the applied load before reconfiguration.

The Figure [START_REF] Guenab | Design of a fault tolerant control system incorporating reliability analysis under dynamic behaviour constraints[END_REF], shows the outputs responses in the faulty case. It can be seen that even the presence of fault, the controlled outputs track the desired outputs without steadystate error. In Figure [START_REF] He | Robust H∞ filtering for timedelay systems with probabilistic sensor faults[END_REF], the control inputs are presented both for the classical controller compared to the proposed controller. Again, the use of the most sensitive actuators and the faulty one is minimized. This distribution maximizes the overall system reliability and improves the dependability of the reconfigurable system where, as shown the Figure (3.b), the life time of the proposed configuration is greater. This result means that the probability that the system remains operational is larger.

VI. CONCLUSION

In this work, a fault-tolerant control design based on actuators health analysis is proposed. The aim is to contribute to a reliable reconfigurable systems where the distribution of the desired efforts is performed by informations about the actuators reliability in addition to the occurred faults. The obtained results show the performance of the proposed approach in term of overall system reliability and the life time of the reconfigurable system. In fact, an effective management of the desired control inputs taking into account the actuators health improves the overall system reliability and the system dependability by keeping the system operational longer in time.

Fig. 1 .

 1 Fig. 1. Aircraft trajectory in the nominal case

Fig. 2 .

 2 Fig. 2. Control inputs in the nominal case

Fig. 3 .

 3 Fig. 3. Comparison between the reliabilities

Fig. 4 .

 4 Fig. 4. Aircraft trajectory in the faulty case

Fig. 5 .

 5 Fig. 5. Control inputs in the faulty case

This work was supported by the SAFFE project 2010-2011 (Ageing Management in fault-tolerant control system design project) from GIS 3SGS -France (https://www.gis-3sgs.fr/).


Alain Bensoussan

Laurent Mertz

An analytic approach to the ergodic theory of stochastic variational inequalities 1

In an earlier work made by the first author with J. Turi (Degenerate Dirichlet Problems Related to the Invariant Measure of Elasto-Plastic Oscillators, AMO, 2008), the solution of a stochastic variational inequality modeling an elasto-perfectly-plastic oscillator has been studied. The existence and uniqueness of an invariant measure have been proven. Nonlocal problems have been introduced in this context. In this work, we present a new characterization of the invariant measure. The key finding is the connection between nonlocal PDEs and local PDEs which can be interpreted with short cycles of the Markov process solution of the stochastic variational inequality.

Résumé

Une approche analytique de la théorie ergodique des inéquations variationnelles stochastiques. Dans un travail précédent du premier auteur en collaboration avec Janos Turi (Degenerate Dirichlet Problems Related to the Invariant Measure of Elasto-Plastic Oscillators, AMO, 2008), la solution d'une inéquation variationnelle stochastique modélisant un oscillateur élastique-parfaitement-plastique a été étudiée. L'existence et l'unicité d'une mesure invariante ont été prouvées. Des problèmes nonlocaux ont été introduits dans ce contexte. La conclusion importante est la connexion entre des EDPs nonlocales et des EDPs locales qui peuvent être interprétées comme les cycles courts du processus de Markov solution de l'inéquation variationnelle stochastique.

Version française abrégée

La dynamique de l'oscillateur élastique-parfaitement-plastique s'exprime à l'aide d'une équation à mémoire (voir (1)-( 2)). A. Bensoussan et J. Turi ont montré que la relation entre la vitesse et la composante élastique de l'oscillateur est un processus de Markov ergodique qui satisfait une inéquation variationnelle stochastique (voir (3)). La solution admet une mesure invariante caractérisée par dualité à l'aide d'une équation aux dérivées partielles avec des conditions de bord non-locales (voir (4)). Dans ce travail, une nouvelle preuve de la théorie ergodique est présentée ainsi qu'une nouvelle caractérisation de l'unique distribution invariante. Dans ce contexte, nous déduisons des nouvelles formules reliant des équations aux dérivées partielles avec des conditions de bord non-locales à des problèmes locaux (voir (10)).

Introduction

In the engineering literature, the dynamics of the elastic-perfectly-plastic (EPP) oscillator has been formulated as a process x(t) which stands for the displacement of the oscillator, evolving with hysteresis. The evolution is defined by the problem

ẍ + c 0 ẋ + F(x(s), 0 ≤ s ≤ t) = ẇ (1)
with initial conditions of displacement and velocity x(0) = x, ẋ(0) = y. Here c 0 > 0 is the viscous damping coefficient, k > 0 the stiffness, w is a Wiener process; F(x(s), 0 ≤ s ≤ t) is a nonlinear functional which depends on the entire trajectory {x(s), 0 ≤ s ≤ t} up to time t. The plastic deformation denoted by ∆(t) at time t can be recovered from the pair (x(t), F(x(s), 0 ≤ s ≤ t)) by the following relationship:

F(x(s), 0 ≤ s ≤ t) =    kY if x(t) = Y + ∆(t), k(x(t) -∆(t)) if x(t) ∈] -Y + ∆(t), Y + ∆(t)[, -kY if x(t) = -Y + ∆(t).
(

) 2
where Y is an elasto-plastic bound. Such elasto-plastic oscillator is simple and representative of the elastoplastic behavior of a class of structure dominated by their first mode of vibration, they are employed to estimate prediction of failure of mechanical structures. Karnopp & Scharton [START_REF] Karnopp | Plastic deformation in random vibration[END_REF] proposed a separation between elastic states and plastic states and introduced a fictitious variable z(t

) := x(t) -∆(t).
Recently, the right mathematical framework of stochastic variational inequalities (SVI) modeling an EPP oscillator with noise has been introduced by one of the authors in [START_REF] Bensoussan | Degenerate Dirichlet Problems Related to the Invariant Measure of Elasto-Plastic Oscillators[END_REF]. Although SVI have been already studied in [START_REF] Bensoussan | Contrôle impulsionnel et inéquations quasi variationnelles[END_REF] to represent reflection-diffusion processes in convex sets, no connection with random vibration had been made so far. The inequality governs the relationship between the velocity y(t) and the variable z(t):

dy(t) = -(c 0 y(t) + kz(t))dt + dw(t), (dz(t) -y(t)dt)(φ -z(t)) ≥ 0, ∀|φ| ≤ Y, |z(t)| ≤ Y. (3)
Let us introduce some notations. In [START_REF] Bensoussan | Degenerate Dirichlet Problems Related to the Invariant Measure of Elasto-Plastic Oscillators[END_REF], it has been shown that the probability distribution of (y(t), z(t)) converges to an asymptotic probability measure on D ∪ D + ∪ D -namely ν. Moreover, ν is the unique invariant distribution of (y(t), z(t)). In addition, from [START_REF] Bensoussan | On a class of partial differential equations with nonlocal Dirichlet boundary conditions[END_REF] we know also that there exists a unique solution u λ to the following partial differential equation (PDE)

λu λ + Au λ = f in D, λu λ + B + u λ = f in D + , λu λ + B -u λ = f in D - (4)
with the nonlocal boundary conditions given by the fact that u λ (y, Y ) and u λ (y, -Y ) are continuous, where λ > 0 and f is a bounded measurable function. The function

u λ satisfies u λ ∞ ≤ f ∞ λ , u λ is continuous and for all (η, ζ) ∈ D, we have lim λ→0 λu λ (η, ζ) = ν(f ).
We use the notation u λ (y, z; f ). Now, we introduce the short cycles to provide a new proof of the ergodic theory for [START_REF] Bensoussan | On a class of partial differential equations with nonlocal Dirichlet boundary conditions[END_REF]. In this context, we derive new formulas linking PDEs with nonlocal boundary conditions to local problems.

Short cycles

Let λ > 0, consider v λ (y, z) the solution of

λv λ + Av λ = f in D, λv λ + B + v λ = f in D + , λv λ + B -v λ = f in D - (5)
with the local boundary conditions v λ (0 + , Y ) = 0 and v λ (0 -, -Y ) = 0. Also, if f is symmetric (resp. antisymmetric) then v λ is symmetric (resp. antisymmetric). We use the notation v λ (y, z; f ). As λ → 0,

v λ → v with Av = f in D, B + v = f in D + , B -v = f in D - (P v )
with the local boundary conditions v(0 + , Y ) = 0 and v(0 -, -Y ) = 0. We use the notation v(y, z; f ). We call v(y, z; f ) a short cycle. We detail the solution of (P v ) in the next section. We introduce next π + (y, z) and π -(y, z) such that

Aπ + = 0 in D, π + = 1 in D + , π + = 0 in D - (6)
and

Aπ -= 0 in D, π -= 0 in D + , π -= 1 in D -. (7)
We have π + + π -= 1, so the existence and uniqueness of a bounded solution to (6) and (7) are clear. A new formulation of the invariant distribution is given by the following theorem.

Theorem 1.1 (New formulation of the invariant distribution ν). Let f be a bounded measurable function on D, we have the following analytical characterization of the invariant distribution:

ν(f ) = v(0 -, Y ; f ) + v(0 + , -Y ; f ) 2v(0 + , Y ; 1) . Denote ν λ (f ) := v λ (0 -,Y ;f )+v λ (0 + ,-Y ;f ) 2v λ (0 -,Y ;1)
. As λ → 0,

u λ (y, z; f ) - ν λ (f ) λ → u(y, z; f ), ν λ (f ) → ν(f ) (8)
where u satisfies

Au = f -ν(f ) in D, B + u = f -ν(f ) in D + , B -u = f -ν(f ) in D - (9)
with the nonlocal boundary conditions given by the fact that u(y, Y ) and u(y, -Y ) are continuous.

Then, we obtain also the representation formula

u(y, z; f ) = v(y, z; f ) -ν(f )v(y, z; 1) + π + (y, z) -π -(y, z) 4π -(0 -, Y ) (v(0 -, Y ; f ) -v(0 + , -Y ; f )) (10)

Analysis of the short cycles

We describe the solution of (P v ). We can write v(y, z; f ) = v e (y, z; f

) + v + (y, z; f ) + v -(y, z; f ) with v e , v + , v -satisfying Av e = f (y, z) in D, v e = 0 in D + , v e = 0 in D -, ( 11
)
Av + = 0 in D, v + (y, Y ) = ϕ + (y; f ) in D + , v + = 0 in D -, (12)
and

Av -= 0 in D, v -= 0 in D + , v -(y, -Y ) = ϕ -(y; f ) in D -. ( 13
)
where ϕ + (y; f ) and ϕ -(y; f ) are defined by

- 1 2 ϕ + yy + (c 0 y + kY )ϕ + y = f (y, Y ), y > 0, ϕ + (0 + ; f ) = 0 (14) and - 1 2 ϕ - yy + (c 0 y -kY )ϕ - y = f (y, -Y ), y < 0, ϕ -(0 -; f ) = 0. (15)
We check easily the formula ϕ

+ (y; f ) = 2 ∞ 0 dξ exp(-(c 0 ξ 2 + 2kY ξ)) ξ+y ξ f (ζ; Y ) exp(-2c 0 ξ(ζ -ξ))dζ, if y ≥ 0 and also ϕ -(y; f ) = 2 ∞ 0 dξ exp(-(c 0 ξ 2 -2kY ξ)) -ξ y-ξ f (ζ; -Y ) exp(-2c 0 ξ(ζ -ξ))dζ, if y ≤ 0.

Ergodic Theorem

Proof (Proof of Theorem 1.1). We first prove the result when f is symmetric. In that case, we can write

u λ (y, z; f ) = v λ (y, z; f ) + v λ (0 -, Y ; f ) v λ (0 -, Y ; 1) 1 λ -v λ (y, z; 1) (22)
Indeed, we know that u λ (y, z; f ) and v λ (y, z; f ) are symmetric. Setting ũλ (y, z; f ) = u λ (y, z; f )v λ (y, z; f ), we obtain

λũ λ + Aũ λ = 0 in D, λũ λ + B + ũλ = 0 in D + , λũ λ + B -ũλ = 0 in D - (23)
with the boundary conditions ũλ (0

+ , Y ; f )-ũ λ (0 -, Y ; f ) = v λ (0 -, Y ; f ) and ũλ (0 + , -Y ; f )-ũ λ (0 -, -Y ; f ) = -v λ (0 + , -Y ; f ).
This last condition is automatically satisfied, thanks to the previous one and the symmetry. The function 1 λv λ (y, z; 1) satisfies the three partial differential equations on D, D + and D -. So, ũλ = C 1 λv λ (y, z; 1) and writing the first boundary condition, we have ũλ (0 1) and formula ( 22) has been obtained. Now, we have

+ , Y ; f ) -ũλ (0 -, Y ; f ) = -C (v λ (0 + , Y ; 1) -v λ (0 -, Y ; 1)) = Cv λ (0 -, Y ; 1). Hence, C = v λ (0 -,Y ;f ) v λ (0 -,Y ;
ν λ (f ) → ν(f ) = v(0 -,Y ;f ) v(0 -,Y ;1) , as λ → 0. If we define u ⋆ λ (y, z; f ) = u λ (y, z; f ) -ν λ (f ) λ = v λ (y, z; f )-ν λ (f )v λ (y, z; 1). The function u ⋆ λ (y, z; f ) → v(y, z; f )-ν(f )v(y, z; 1) = v(y, z; f -ν(f )), λ → 0. Also from its definition the function u ⋆ λ (y, Y ; f ) and u ⋆ λ (y, -Y ; f ) are continuous. From the choice of ν(f ) the function v(y, Y ; f -ν(f )) is continuous. Now, since f -ν(f ) is symmetric v(0 + , -Y ; f -ν(f )) - v(0 -, -Y ; f -ν(f )) = v(0 + , Y ; f -ν(f )) -v(0 -, Y ; f -ν(f )) = 0.
So the result is proven when f is symmetric. We now consider the situation when f is antisymmetric. We know that u λ (y, z; f ) is antisymmetric. Similarly v λ (y, z; f ) is antisymmetric. Consider π - λ and π + λ defined by

λπ + λ + Aπ + λ = 0 in D, λπ + λ + B + π + λ = 0 in D + , π + λ = 0 in D - (24)
with the boundary condition π + λ (0 + , Y ) = 1 and

λπ - λ + Aπ - λ = 0 in D, π - λ = 0 in D + , λπ - λ + B + π - λ = 0 in D - (25)
with the boundary condition π - λ (0 -, -Y ) = 1. We have π - λ (y, z) = π - λ (-y, -z), we then state the formula u λ (y, z

; f ) = v λ (y, z; f ) - (π + λ (y,z)-π - λ (y,z))v λ (0 + ,-Y ;f ) 1-π + λ (0 -,Y )+π + λ (0 + ,-Y )
. So we see that u λ (y, z; f ) converges as λ → 0,without substracting a number ν λ (f ) λ . The function u λ (y, z; f ) converges pointwise to u(y, z; f ) = v(y, z; f ) -(π + (y,z)-π -(y,z))v(0 + ,-Y ;f ) 2π -(0 -,Y )

. So when f is antisymmetric, the results (8)-(9) hold with ν λ (f ) = 0 and ν(f ) = 0.

For the general case, we can write f = f sym + f asym with f sym (y, z) = f (y,z)+f (-y,-z) 2 , f asym (y, z) = f (y,z)-f (-y,-z) 2

. We have ν(f sym ) = . We obtain also the representation formula u(y, z; f ) = v(y, z; f )ν(f )v(y, z; 1) + π + (y,z)-π -(y,z)

4π -(0 -,Y )

(v(0 -, Y ; f )v(0 + , -Y ; f )) and the result is obtained.

Notation 1 .

 1 D := R × (-Y, +Y ), D + := (0, ∞) × {Y }, D -:= (-∞, 0) × {-Y }, and the differential operators Aζ := -1 2 ζ yy +(c 0 y +kz)ζ y -yζ z , B + ζ := -1 2 ζ yy +(c 0 y +kY )ζ y , B -ζ := -1 2 ζ yy +(c 0 y -kY )ζ y . where ζ is a regular function on D.

v( 0

 0 -,Y ;fsym) v(0 -,Y ;1) and thus ν(f sym ) = v(0 -,Y ;f )+v(0 + ,-Y ;f ) 2v(0 + ,-Y ;1) Since ν(f asym ) = 0, we deduce ν(f ) = ν(f sym ) = v(0 -,Y ;f )+v(0 + ,-Y ;f ) 2v(0 + ,-Y ;1)

This research was partially supported by a grant from CEA, Commissariat à l'énergie atomique and by the National Science Foundation under grant DMS-0705247. A large part of this work was completed while one of the authors was visiting the University of Texas at Dallas and the Hong-Kong Polytechnic University. We wish to thank warmly these institutions for the hospitality and support.Email addresses: alain.bensoussan@utdallas.edu (Alain Bensoussan), mertz@ann.jussieu.fr (Laurent Mertz)

This research in the paper was supported by WCU (World Class University) program through the National Research Foundation of Korea funded by the Ministry of Education, Science and Technology (R31 -20007).

Solution to Problem (11)

The proof will be based on solving a sequence of Interior Exterior Dirichlet problems and a fixed point argument. Thus, we need to state the two following lemmas as preliminary results. It is sufficient to consider f = 1, with no loss of generality.

Interior Dirichlet problem

We begin with the interior problem, let

We consider the problem

Lemma 2.1. There exists a unique bounded solution to the equation (16).

Proof. It is sufficient to prove an a priori bound. For that we can assume ϕ + , ϕ -= 0. Consider λ > 0 and the function θ(y, z) = exp(λc 0 (y 2 + kz 2 )) then -

If we pick λ > max(1, 1 c0 ) the right hand side of (17) is positive. Therefore the minimum of H can occur only on the boundary y = ȳ1 and z = Y with y > 0 or z = -Y with y < 0. It follows that H(y, z) ≥ -exp(λc 0 (ȳ 2 1 + Y 2 )) and thus also 0

Exterior Dirichlet problems

Now, we proceed by considering two exterior Dirichlet problems. Let 0 < ȳ < ȳ1 , we define

We use the same notation η(y, z) for the two problems (18),(19) for the convenience of the reader. We have Lemma 2.2. For any ȳ > 0 there exists a unique bounded solution of (18),(19).

Proof. It is sufficient to prove the bound, we claim that Proof. Uniqueness comes from maximum principle. Setting Φ = (ϕ + (z), ϕ -(z)) and using the notation

We thus have defined a map Γ from C + 1 , C - 1 into itself. If Γ has a fixed point, then it is clear that the function

and the required regularity is available at boundary points ȳ, ȳ1 , -ȳ, -ȳ 1 . The result will follow from the property : Γ is a contraction mapping. This property will be an easy consequence of the following result.

Consider the exterior problem

where

, then the maximum is attained on the line y = ȳ1 , and this is impossible because it cannot be at z = Y , nor at z = -Y , nor at the interior, by maximum principle considerations.

Solution to Problems (12) and (13)

We now consider the function ϕ + and ϕ -solution of ( 14) and ( 15). Note that if y < 0, we have ϕ -(y; 1) = ϕ + (-y; 1). So it is sufficient to consider (14) and we easily see that

), if y > 0. We next want to solve the problem (12). We proceed as follow. We extend ϕ + for y < 0, by a function which is C 2 on R and with compact support on y < 0. It is convenient to call ϕ(y) the C 2 function on R, with compact support for y < 0 and ϕ(y) = ϕ + (y; 1) for y > 0. We set w + (y, z) = v + (y, z)ϕ(y) then we obtain the problem

with g(y, z) = --1 2 ϕ yy + (c 0 y + kz)ϕ y .

But, g(y, z) = 1 {y>0} (-1 + k(Yz)ϕ y (y)) + 1 {y<0} -(-1 2 ϕ yy + (c 0 y + kz)ϕ y ) and thus, taking into account the definition of ϕ when y < 0, we can assert that g(y, z) is a bounded function. Again, from the definition of ϕ(y) when y < 0, we obtain that on the boundary, w + is bounded. It follows from what was done for Problem (11) that (21) has a unique solution. So we can state the following proposition. Proposition 2.2. There exists a unique solution to (12) of the form v + (y, z) = ϕ + (y)1 {y>0} + ṽ+ (y, z) where ṽ+ (y, z) is bounded. Similarly, there exists a unique solution to (13) of the form v -(y, z) = ϕ -(y)1 {y<0} + ṽ-(y, z) where ṽ-(y, z) is bounded.

Proof. We just define ϕ(y) extension of ϕ + (y) for y < 0 as explained before and consider w + (y, z) solution of (21). We know that w + (y, z) is bounded and we have v + (y, z) = ϕ(y) + w + (y, z) = ϕ + (y)1 {y>0} + ϕ(y)1 {y<0} + w + (y, z) which is of the form (12) with ṽ+ (y, z) = ϕ(y)1 {y<0} + w + (y, z).

The complete Problem (P v )

Finally, we consider the complete Problem (P v ), we can state Theorem 2.1. There exists a unique solution of (P v ) of the form v(y, z; f ) = ϕ + (y; f )1 {y>0} +ϕ -(y; f )1 {y<0} + w(y, z) where w(y, z) is a bounded function which can be written as w = v e + w + + w -.

Proof. We just collect the results of Propositions 2.1 and 2.2.


Alain Bensoussan

Laurent Mertz

Phillip Yam

Behavior of the plastic deformation of an elasto-perfectly-plastic oscillator with noise

Earlier works in engineering, partly experimental, partly computational have revealed that asymptotically, when the excitation is a white noise, plastic deformation and total deformation for an elasto-perfectly-plastic oscillator have a variance which increases linearly with time with the same coefficient. In this work, we prove this result and we characterize the corresponding drift coefficient. Our study relies on a stochastic variational inequality governing the evolution between the velocity of the oscillator and the non-linear restoring force. We then define long cycles behavior of the Markov process solution of the stochastic variational inequality which is the key concept to obtain the result. An important question in engineering is to compute this coefficient. Also, we provide numerical simulations which show successful agreement with our theoretical prediction and previous empirical studies made by engineers.

Résumé

Le comportement de la déformation plastique pour un oscillateur élastique-parfaitementplastique excité par un bruit blanc. Des résultats expérimentaux en sciences de l'ingénieur ont montré que, pour un oscillateur élasto-plastique-parfait excité par un bruit blanc, la déformation plastique et la déformation totale ont une variance, qui asymptotiquement, croît linéairement avec le temps avec le même coefficient. Dans ce travail, nous prouvons ce résultat et nous caractérisons le coefficient de dérive. Notre étude repose sur une inéquation variationnelle stochastique gouvernant l'évolution entre la vitesse de l'oscillateur et la force de rappel non-linéaire. Nous définissons alors le comportement en cycles longs du processus de Markov solution de l'inéquation variationnelle stochastique qui est le concept clé pour obtenir le résultat. Une question importante en sciences de l'ingénieur est de calculer ce coefficient. Les résultats numériques confirment avec succès notre prédiction théorique et les études empiriques faites par les ingénieurs.

Version française abrégée

Dans cet article, nous étudions la variance de l'oscillateur élastique-parfaitement-plastique (EPP) excité par un bruit blanc. La dynamique de l'oscillateur s'exprime à l'aide d'une équation à mémoire. (voir

Introduction

In civil engineering, an elasto-perfectly-plastic (EPP) oscillator with one single degree of freedom is employed to estimate prediction of failure of mechanical structures under random vibrations. This elastoplastic oscillator consists in a one dimensional model simple and representative of the elasto-plastic behavior of a class of structure dominated by their first mode of vibration. The main difficulty to study these systems comes from a frequent occurrence of plastic phases on small intervals of time. A plastic deformation is produced when the stress of the structure crosses over an elastic limit. The dynamics of the EPP-oscillator has memory, so it has been formulated in the engineering literature as a process with hysteresis x(t), which stands for the displacement of the oscillator. We study the problem

ẍ + c 0 ẋ + F t = ẇ (1)
with initial conditions of displacement and velocity x(0) = 0, ẋ(0) = 0. Here c 0 > 0 is the viscous damping coefficient, k > 0 the stiffness, w is a Wiener process; F t := F(x(s), 0 ≤ s ≤ t) is a non-linear functional which depends on the entire trajectory {x(s), 0 ≤ s ≤ t} up to time t. The plastic deformation denoted by ∆(t) at time t can be recovered from the pair (x(t), F t ) by the following relationship:

F t =    kY if x(t) = Y + ∆(t), k(x(t) -∆(t)) if x(t) ∈] -Y + ∆(t), Y + ∆(t)[, -kY if x(t) = -Y + ∆(t). (2)
where ∆(t) = t 0 y(s)1 {|Ft|=kY } ds and Y is an elasto-plastic bound. Karnopp & Scharton [START_REF] Karnopp | Plastic deformation in random vibration[END_REF] proposed a separation between elastic states and plastic states. They introduced a fictitious variable z(t) := x(t) -∆(t) and noticed the simple fact that between two plastic phases z(t) behaves like a linear oscillator. In addition, other works made by engineers [START_REF] Borsoi | Approche probabiliste de la ruine d'un oscillateur élasto-plastique sous séisme. 2ème colloque national de l'AFPS[END_REF], partly experimental, partly computational, have revealed that total deformation has a variance which increases linearly with time:

lim t→∞ σ 2 (x(t)) t = σ 2 . ( 3
)
where σ 2 is a positive real number.

Recently, the right mathematical framework of stochastic variational inequalities (SVI) modeling an elasto-plastic oscillator with noise (presented below) has been introduced by the first author with J. Turi in [START_REF] Bensoussan | Degenerate Dirichlet Problems Related to the Invariant Measure of Elasto-Plastic Oscillators[END_REF]. The inequality governs the relationship between the velocity y(t) and the elastic component z(t):

dy(t) = -(c 0 y(t) + kz(t))dt + dw(t), (dz(t) -y(t)dt)(φ -z(t)) ≥ 0, ∀|φ| ≤ Y, |z(t)| ≤ Y.
(SVI)

The plastic deformation ∆(t) can be recovered by t 0 y(s)1 {|z(s)|=Y } ds.

Throughout the paper, the objective is to prove (3) and to provide a exact characterization of (3) based on (SVI). Now, we introduce long cycles in view to the fact that we are interested in identifying independent sequences in the trajectory.

Long cycles

Denote τ 0 := inf{t > 0, y(t) = 0 and |z(t)| = Y } and s := sign(z(τ 0 )) which labels the first boundary hit by the process (y(t), z(t)). Define θ 0 := inf{t > τ 0 , y(t) = 0 and z(t) = -sY }. In a recurrent manner, knowing θ n , we can define for n ≥ 0 τ n+1 := inf{t > θ n , y(t) = 0 and z(t) = sY }, θ n+1 := inf{t > τ n+1 , y(t) = 0 and z(t) = -sY }. Now, according to the previous setting we can define the n-th long cycle (resp. first part of the cycle, second part of the cycle) as the piece of trajectory delimited by the interval [τ n , τ n+1 ), (resp. [τ n , θ n+1 ) and [θ n+1 , τ n+1 )). Indeed, at the instants {τ n , n ≥ 1}, the process (y(t), z(t)) is in the same state at the instant τ 0 . In addition, there are two types of cycles depending on the sign of s. The set of stopping times {τ n , n ≥ 0} represents the times of occurence of long cycles. Let us remark that the plastic deformation and the total deformation are the same on a long cycle since

τ1 τ0 y(t)dt = τ1 τ0 y(t)1 {|z(s)|=Y } dt + τ1 τ0 y(t)1 {|z(s)|<Y } dt and that τ1 τ0 y(t)1 {|z(s)|<Y } dt = z(τ 1 ) -z(τ 0 ) = 0.
As main result, we have obtained the following theorem. Theorem 1.1 (Characterization of the variance related to the plastic/total deformation). In the previously defined context, we have shown

lim t→∞ σ 2 (x(t)) t = E τ1 τ0 y(t)dt 2 E(τ 1 -τ 0 ) . (4)
Our proofs are based on solving nonlocal partial differential equations related to long cycles. Simpler formula will be given below at equation (13).

The issue of long cycles and plastic deformations

Let us introduce notations. Let f be a bounded measurable function, we want to solve

Av + = f (y, z) in D, B + v + = f (y, Y ) in D + , B -v + = f (y, -Y ) in D - (P v + )
with the nonlocal boundary conditions v + (y, Y ) continuous and v + (0 -, -Y ) = 0, and

Av -= f (y, z) in D, B + v -= f (y, Y ) in D + , B -v -= f (y, -Y ) in D - (P v -)
with the nonlocal boundary conditions v -(0 + , Y ) = 0 and v -(y, -Y ) continuous. The functions v + (y, z; f ) and v -(y, z; f ) are called long cycles. In addition, we define π + (y, z) and π -(y, z) by

Aπ + = 0 in D, π + (y, Y ) = 1 in D + , π + (y, -Y ) = 0 in D - (P π + )
and

Aπ -= 0 in D, π -(y, Y ) = 0 in D + , π -(y, -Y ) = 1 in D - (P π -)
Note that π + (y, z) + π -(y, z) = 1, so the existence and uniqueness of a bounded solution (P π + ) and (P π -) is clear.

Proposition 2.1. We have the properties π -(0 + , -Y ) > 0 and π + (0 -, Y ) > 0.

Proof. We check only the first property. We consider the elastic process (y yz (t), z yz (t)):

z yz (t) = e -c 0 t 2 {z cos (ωt) + 1 ω (y + c 0 2 z) sin (ωt)} + 1 ω t 0 e -c 0 2 (t-s) sin (ω(t -s))dw(s), y yz (t) = - c 0 2 z yz (t) + e -c 0 t 2 {-ωz sin (ωt) + (y + c 0 2 z) cos (ωt)} + t 0 e -c 0 2 (t-s) cos (ω(t -s))dw(s)
where

ω := √ 4k-c 2 0 2
(we assume 4k > c 2 0 ). Note that the condition 4k > c 2 0 is needed so that (y yz (t), z yz (t)) have real valued solutions. Set

τ yz := inf{t > 0, |z yz (t)| ≥ Y } then we have the probabilistic interpretation π + (y, z) = P(z yz (τ yz ) = Y ), π -(y, z) = P(z yz (τ yz ) = -Y ). We can state π -(y, z) → 1 as y → -∞, z ∈ [-Y, Y ]. (5)
Indeed, ∀t with 0 < t < π ω we have z yz (t) → -∞, as y → -∞ a.s. Therefore ∀t with 0 < t < π ω , a.s. z yz (t) < -Y for y sufficiently large. Hence, a.s. τ yz < t for y sufficiently large. Therefore a.s. lim sup y→-∞ τ yz < t. Since t is arbitrary, necessarily a.s. τ yz → 0, as y → -∞ which implies [START_REF] Karnopp | Plastic deformation in random vibration[END_REF]. Moreover the function π + (y, Y ) cannot have a minimum or a maximum at any finite y < 0. It is then monotone and since

π -(-∞, Y ) = 1, it is monotone decreasing. It follows that π -(0 -, Y ) < 1. It cannot be 0. Otherwise, π -(y, Y
) is continuous at y = 0, and (0, Y ) is a point of minimum of π -(y, z). Since for y < 0, π - y (y, Y ) < 0 from the equation of π -we get lim sup y→0 π - yy (y, Y ) ≤ 0 which is not possible since (0, Y ) is a minimum. We next define η(y, z) by

- 1 2 η yy + (c 0 y + kz)η y -yη z = f (y, z) in D, η(y, Y ) = 0 in D + , η(y, -Y ) = 0 in D -(P η )
with the local boundary conditions η(0 + , Y ) = 0 and η(0 -, -Y ) = 0. For f bounded (P η ) has a unique bounded solution. Then define ϕ + (y; f ) by solving

- 1 2 ϕ +,yy + (c 0 y + kY )ϕ +,y = f (y, Y ), y > 0, ϕ + (0; f ) = 0. (6)
We check easily the formula

ϕ + (y; f ) = 2 ∞ 0 dξ exp(-(c 0 ξ 2 + 2kY ξ)) ξ+y ξ f (ζ; Y ) exp(-2c 0 ξ(ζ -ξ))dζ, y ≥ 0. (7)
We consider ψ + (y, z; f ) defined by

   -1 2 ψ +,yy + (c 0 y + kz)ψ +,y -yψ +,z = 0 in D, ψ + (y, Y ) = ϕ + (y; f ) in D + , ψ + (y, -Y ) = 0 in D - (P ψ+ ) and similarly ϕ -(y; f ), ψ -(y, z; f ) defined by - 1 2 ϕ -,yy + (c 0 y -kY )ϕ -,y = f (y, -Y ), y < 0, ϕ -(0; f ) = 0. (8)
which leads to

ϕ -(y; f ) = 2 ∞ 0 dξ exp(-(c 0 ξ 2 -2kY ξ)) -ξ y-ξ f (ζ; -Y ) exp(-2c 0 ξ(ζ -ξ))dζ, y ≤ 0. and    -1 2 ψ -,yy + (c 0 y + kz)ψ -,y -yψ -,z = 0 in D, ψ -(y, Y ) = 0 in D + , ψ -(y, -Y ) = ϕ -(y; f ) in D - (P ψ-)
We can state Proposition 2.2. The solution of (P v + ) is given by

v + (y, z; f ) =η(y, z; f ) + ψ + (y, z; f ) + ψ -(y, z; f ) (9) + η(0 -, Y ; f ) + ψ + (0 -, Y ; f ) + ψ -(0 -, Y ; f ) π -(0 -, Y ) π + (y, z).
and the solution of (P v -) is given by

v -(y, z; f ) =η(y, z; f ) + ψ + (y, z; f ) + ψ -(y, z; f ) (10) + η(0 + , -Y ; f ) + ψ + (0 + , -Y ; f ) + ψ -(0 + , -Y ; f ) π + (0 -, Y ) π -(y, z).
Proof. Direct checking.

Complete cycle

First, let us check that E[x(t)] = 0 and then σ 2 (x(t)) = E[x 2 (t)]. Indeed, by symmetry of the inequality (SVI) and by the choice of initial conditions y(0) = 0, z(0) = 0, the processes (y(t), z(t)) and -(y(t), z(t)) have same law. Then, E[x(t)] = E t 0 y(s)ds = 0. In addition, (y(τ 1 ), z(τ 1 )) is equal to (0, -Y ) or (0, Y ) with probability 1/2 for both, therefore with no loss of generality we can suppose that (y(0), z(0)) = (0, -Y ) or (0, Y ) with probability 1/2 for both. Let us treat the case y(0

) = 0, z(0) = Y . So, τ 0 = 0 and θ 1 = inf{t > 0, z(t) = -Y, y(t) = 0}. We can assert that Eθ 1 = v + (0, Y ; 1) hence θ 1 < ∞ a.s.. Next we define τ 1 = inf{t > θ 1 , z(t) = Y, y(t) = 0} then Eτ 1 = v + (0, Y ; 1) + v -(0, -Y ; 1) = 2v + (0, Y ; 1)
. At time τ 1 the state of the system is again (0, Y ). So the sequence {τ n , n ≥ 0} is such that τ n < τ n+1 and in the interval (τ n , τ n+1 ) we have a cycle identical to (0, τ 1 ). Consider the variable τ1 0 y(t)dt. We have

E τ1 0 y(t)dt = E θ1 0 y(t)dt + E τ1 θ1 y(t)dt = v + (0, Y ; y) + v -(0, -Y ; y). However if f is antisymmetric f (-y, -z) = -f (y, z) we have v + (0, Y ; f ) = -v -(0, -Y ; f ), therefore E τ1 0 y(t)dt = 0. Hence, E τn 0 y(t)dt 2 = E n-1 j=0 τj+1 τj y(t)dt 2 = nE τ1 0 y(t)dt 2 . Then E τn 0 y(t)dt 2 Eτ n = E τ1 0 y(t)dt 2 Eτ 1 .
Let T > 0 and N T with τ NT ≤ T < τ NT +1 , N T = 0 if τ 1 > T, τ 0 = 0 and then calculations lead to the following

E τN T +1 0 y(t)dt 2 Eτ NT +1 = E τ1 0 y(t)dt 2 Eτ 1 .
Next, we can justify that lim

T →∞ 1 T E T 0 y(t)dt 2 = lim T →∞ 1 T E τN T +1 0 y(t)dt 2 (11)
and that we have a lower bound and an upper bound for (11), that is

E τ1 0 y(t)dt 2 Eτ 1 ≤ 1 T E τN T +1 0 y(t)dt 2 ≤ 1 + Eτ 1 T E τ1 0 y(t)dt 2 Eτ 1 . Therefore lim T →∞ 1 T E T 0 y(t)dt 2 = E τ1 0 y(t)dt 2 Eτ 1 . (12)
Moreover, we can simplify (12), indeed we have Since Eτ 1 = 2v + (0, Y ; 1) we obtain lim

T →∞ 1 T E T 0 y(t)dt 2 = E θ1 0 y(t)dt 2 -(v + (0, Y ; y)) 2 v + (0, Y ; 1) . ( 13
)
4. Numerical evidence in support of our result

In this section, we provide computational results which confirm our theoretical results. A C code has been written to simulate (y(t), z(t)). See [START_REF] Bensoussan | An Ultra Weak Finite Element Method as an Alternative to a Monte Carlo Method for an Elasto-Plastic Problem with Noise[END_REF] for the numerical scheme considered to do direct simulation. Let T > 0, N ∈ N and {t n = nδt, 0 ≤ n ≤ N } where δt = T N . Then, to compute the left hand side of (4), we consider M C ∈ N and we generate M C numerical solutions of (SVI) {y i (t), 0 ≤ t ≤ T, 1 ≤ i ≤ M C} up to the time T . By the law of large numbers, we can approximate

1 T E T 0 y(s)ds 2 by X MC := 1 T 1 MC MC i=1 N i=1 y i (t n )δt 2 and also 1 T 2 E T 0 y(s)ds 4 by X 2 MC := 1 T 2 1 MC MC i=1 N i=n y i (t n )δt 4 . De- note s X := X 2 MC -(X MC ) 2
, we also know by the central limit theorem that 1

T E T 0 y(s)ds 2 ∈ [X MC - 1.96sX √ MC , X MC + 1.96sX √ MC ]
with 95% of confidence. Similarly, to compute the right hand side of (4), we generate M C numerical long cycles. For each trajectory {y i (t), t ≥ 0}, we consider N i c the required number of time step to obtain a completed cycle. Denote δ

MC := 1 MC MC i=1 N i c n=0 y i (t n )δt 2 , δ 2 MC := 1 MC MC i=1 N i c n=0 y i (t n )δt 4 , τ MC := 1 MC MC i=1 N i c δt , τ 2 MC := 1 MC MC i=1 (N i c δt) 2 , s δ := δ 2 MC -(δ MC ) 2
and s τ := τ 2 MC -(τ MC ) 2 . We also know that δMC τMC ∈ [ δMC -

1.96s δ √ M C τMC + 1.96sτ √ M C , δMC + 1.96s δ √ M C τMC-1.96sτ √ M C
] with 95% of confidence. In table 1, is shown a comparison of the results obtained for X MC and δMC τMC where T = 500 , δt = 10 -4 and M C = 5000. Remark 4.1. From a numerical point of view, the behavior in long cycles is very relevant. Indeed, computing the left hand side of (4) is much more expensive in time compared to the right hand side (see the τ MC -column of table 1).

c 0 = 1, k = 1 Y X MC , T =

Notation 1 .

 1 D := R × (-Y, +Y ), D + := (0, ∞) × {Y }, D -:= (-∞, 0) × {-Y }, and the differential operators Aζ := -1 2 ζ yy +(c 0 y +kz)ζ y -yζ z , B + ζ := -1 2 ζ yy +(c 0 y +kY )ζ y , B -ζ := -1 2 ζ yy +(c 0 y -kY )ζ y . where ζ is a regular function on D.

2 -

 2 (v + (0, Y ; y)) 2   .

Table 1 :

 1 500 δMC τMC τ MC Relative error % 0.1 0.807 ±0.031 0.834 ±0.069 6.61 ±0.11 3.2 0.2 0.649 ±0.026 0.624 ±0.047 8.74 ±0.13 3.8 0.3 0.493 ±0.020 0.464 ±0.034 10.45 ±0.16 5.8 0.4 0.361 ±0.0140.355 ±0.026 12.12 ±0.18 1.7 0.5 0.266±0.011 0.257 ±0.019 13.80±0.21 Monte-Carlo simulations to compare numerical solution of the left and right hand sides of (4), T = 500 , δt = 10 -4 and M C = 5000.

	3.3

This research was partially supported by a grant from CEA, Commissariat à l'énergie atomique and by the National Science Foundation under grant DMS-0705247. A large part of this work was completed while one of the authors was visiting the University of Texas at Dallas and the Hong-Kong Polytechnic University. We wish to thank warmly these institutions for the hospitality and support.Email addresses: alain.bensoussan@utdallas.edu (Alain Bensoussan), mertz@ann.jussieu.fr (Laurent Mertz)Preprint submitted to Elsevier December 18,

This research in the paper was supported by WCU (World Class University) program through the National Research Foundation of Korea funded by the Ministry of Education, Science and Technology (R31 -20007).