C Sarotte

J Marzat

H Piet-Lahanier

G Ordonneau

M Galeotta

Model-based active fault-tolerant control for a cryogenic combustion test bench

Keywords: Health management system, Fault-tolerant control, Fault detection, Fault isolation, Rocket engine

In this paper a method is proposed to design a fault detection and isolation scheme based on quantitative physics-based models, as well as fault-tolerant control strategy to improve the reliability of a cryogenic combustion bench operation. The detection and isolation scheme is composed of an extended observer, a cumulative sum algorithm and an exponentially weighted moving average chart. In the case of interdependent parts, a dynamic parity space approach is proposed to isolate one or two simultaneous faults with constraints based on the mass flow rate continuity and the energy conservation for the overall system.

The method allows settling adaptive thresholds that avoid pessimistic decision about the continuation of tests while detecting and isolating faults in the system. Then a fault-tolerant system reconfiguration mechanism is provided with a control law which compensates for an estimated actuator additive fault to maintain the overall system stability or allows converging to a reference state in order to overcome instabilities and take into account actuator saturations. For that purpose, the fault-tolerant control algorithm comprises a linear quadratic regulator, an unknown input observer to estimate the fault and an anti-windup scheme. The model and the estimation part were validated on real data from the ONERA/CNES MASCOTTE test bench, and the reconfiguration control law was validated in realistic simulations of the same system.

Introduction

Monitoring engines and test benches is a major challenge in the development and integration of new propulsion systems for rockets, including reusable ones [START_REF] Wu | Liquid-propellant rocket engines health-monitoring: A survey[END_REF]. The increase of launchers operation safety, economical efficiency and launcher's engines reliability [START_REF] Duyar | Fault diagnosis for the space shuttle main engine[END_REF][START_REF] Fekih | Fault diagnosis and fault tolerant control design for aerospace systems: A bibliographical review[END_REF] are major problematic [START_REF] Nie | Liquid-propellant rocket engine online health condition monitoring base on multi-algorithm parallel integrated decisionmaking[END_REF], [START_REF] Hayton | Static and dynamic novelty detection methods for jet engine health monitoring[END_REF], [START_REF] Figueroa | Rocket testing and integrated system health management[END_REF]. For those reasons Health Management Systems (HMS) have to be improved to perform technical and economical optimization of launchers systems, and in particular of engine systems. The use of HMS has been initiated in the early 70's. The objectives are to design efficient, fast and reliable approaches to detect faults of various magnitudes. Those approaches can be divided in two different categories, data-based and model-based ones. In the case of rocket engines, since it is not realistic to collect enough data to use data-based methods, qualitative or quantitative model-based methods are essentially used. However, the use of model-based methods implies the description of complex physical phenomenon as well as the compliance with sensors sensitivity and thermo-mechanical positioning constraints. Moreover, since the developed algorithms have to allow fault detection in real time [START_REF] Betta | Instrument fault detection and isolation: State of the art and new research trends[END_REF] the methods used have to be fast and robust. Then, if a minor component and/or instrument fault is detected by the model-based Fault Detection and Isolation (FDI) approaches [START_REF] Hwang | A survey of fault detection, isolation, and reconfiguration methods[END_REF][START_REF] Yin | A review on recent development of spacecraft attitude fault tolerant control system[END_REF], non-shutdown actions have to be defined to maintain the overall system current performances close to the desirable ones and preserve stability conditions ( [START_REF] Ye | Adaptive fault-tolerant tracking control against actuator faults with application to flight control[END_REF], [START_REF] Yang | Adaptive fault-tolerant H inf control via state feedback for linear systems against actuator faults[END_REF], [START_REF] Yang | Fault tolerant control and hybrid systems[END_REF]). For that reason, it is required to perform a reconfiguration [START_REF] Ioannou | Robust adaptive control[END_REF] of the engine using Fault-Tolerant Control Systems (FTCS). Active FTC Systems are characterized by on-line FDI processes as described by Zhang et al. [START_REF] Zhang | Bibliographical review on reconfigurable fault-tolerant control systems[END_REF], [START_REF] Zhang | Active fault-tolerant control system against partial actuator failures[END_REF]. This system firstly detects and estimates faults, the second step is to achieve a steady-state tracking of the reference input by compensating the fault [START_REF] Theilliol | Actuator fault tolerant control design based on a reconfigurable reference input[END_REF]. For that purpose, FDI methods have been developed to evaluate failures and take a decision using all available information with the help of explicit or implicit models [START_REF] Zhong | A survey on model-based fault diagnosis for linear discrete time-varying systems[END_REF]. The most common model-based approach for FDI makes use of observers to generate residuals as presented by Ding et al. and Hwang et al. [START_REF] Ding | Model based fault diagnosis techniques: design schemes, algorithms, and tools[END_REF], [START_REF] Hwang | A survey of fault detection, isolation, and reconfiguration methods[END_REF]. Faults are then detected by setting a fixed or variable threshold on each residual signals as described by Basseville [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. Those FDI methods assume that the mathematical model used is representative of the system dynamics [START_REF] Frank | Handling modelling uncertainty in fault detection and isolation systems[END_REF][START_REF] Gertler | Fault detection and diagnosis[END_REF]. The methods commonly used nowadays for engine HMS [START_REF] Gubanov | USSR main engines for heavy-lift launch vehicles-status and direction[END_REF][START_REF] Huang | Key reliability drivers of liquid propulsion engines and a reliability model for sensitivity analysis[END_REF] are a basic engine redline system as well as advanced sensors and algorithms including multiple engine parameters that infer an engine anomaly condition from sensor data and take mitigation action accordingly. Basic redlines are straightforward in that they usually act on a single operating parameter anomaly [START_REF] Feng | Research on health evaluation system of liquidpropellant rocket engine ground-testing bed based on fuzzy theory[END_REF]. Those methods can induce false alarms or undetected failures that can be critical for the operation safety and reliability. Moreover, designing representative mathematical models is challenging in practice because of the presence of modeling uncertainties and unknown disturbances [START_REF] Chen | Design of unknown input observers and robust fault detection filters[END_REF], [START_REF] Yang | Observers for linear systems with unknown inputs[END_REF], [START_REF] Bittner | An integrated process for FDIR design in aerospace[END_REF] to which the developed FTCS should be robust. Finally, due to physical actuators characteristics or performances, unlimited control signals are not available and saturations should be taken into account in the control law design.

The development and the evaluation of the performances of FDI methods should rely on real measurements [START_REF] Johnson | A review of fault management techniques used in safetycritical avionic systems[END_REF]. For that purpose, this work is partly based on the exploitation of real data from the MASCOTTE test bench developed at ONERA (see [START_REF] Vingert | Dossier de Définition et de réalisation de Mascotte V05[END_REF]). This test facility is dedicated to the experimental study of cryogenic rocket engines fueled with oxygen and hydrogen or methane.

The obtained measurements will allow updating and adapting the simulation models as well as validating by identification the engine characteristics on offline tests. The different types of faults were simulated with CARINS simulation software (CNES/ONERA). CARINS is a software developed for simulation and modeling with a system based approach (see [START_REF] Ordonneau | CARINS: A Future Versatile and Flexible Tool for Engine Transient Prediction[END_REF]).

The threefold objective consists first in (Fig. 1). The first objective is the modeling of the bench. A first difficulty is to propose a model of evolution of the phenomenon whose states can be identified online and which makes it possible to highlight a change of behavior in a robust and fast way. On the basis of the previous works of Iannetti et al. [START_REF] Iannetti | Development of model-based fault diagnosis algorithms for MASCOTTE cryo-60 genic test bench[END_REF], approaches have been developed to allow the comparison between the evolution of the complete state (pressure, temperature, mass flow) and a prediction under nominal operating hypothesis.

The second objective is to design fault detection and isolation filters based on this modeling [START_REF] Emami-Naeini | Effect of model uncertainty on failure detection: the threshold selector[END_REF], [START_REF] Isaksson | An on-line threshold selector for failure detection[END_REF]. The use of methods based on an adaptive threshold value (ACUSUM) allows to detect a fault regardless of the component of the affected state instead of comparing residuals with predetermined thresholds that are assigned to each system output used in the FDI algorithm (for example see [START_REF] Palmer | Active fault diagnosis with sensor selection in a diesel engine air handling system[END_REF]). In the cooling system of the test bench which is under-monitored, methods to isolate faults based on parity-space approaches have been developed. Those methods are based on fluid mechanical constraints to determine the parity space matrix instead of defining robustness/sensitivity criteria as presented by Zhong et al. [START_REF] Zhong | Parity space-based fault detection for linear discrete time-varying systems with unknown input[END_REF]. The third objective is to define a real-time control system, which would make it possible to overcome certain failures. On the basis of the dynamical modeling of the bench, the first step is to model the link between the input set-points (flow rates, pressures, ...) and the desired operating point, then in a second step to determine the commands to be applied to maintain this operat-ing point when a fault occurs on the bench. A control law allowing to recover the characteristics of the nominal operation after detecting a fault affecting actuators (valves) and to take into account input saturations has been developed.

The nominal control law is a Linear Quadratic Regulator (LQR) type error feedback. Those kind of approaches allow to ensure the system stability around an operating trajectory and to compensate for an additive actuator failure which was not taken into account with usual engine control methods based for example on Proportional-Integral-Derivative (PID) controller (see [START_REF] Fang | Closed-loop combustion phase control for multiple combustion modes by multiple injections in a compression ignition engine fueled by gasoline-diesel mixture[END_REF]). Moreover, the error feedback allows to take into account the state estimation error directly in the control design. An anti-windup scheme has been proposed to account for actuator saturations. In this approach, the set of admissible initial states and its associated domain of stability are determined to take into account the compensation of additive actuator faults. The developed Fault Detection, Isolation and Reconfiguration (FDIR) scheme has been validated on real data of MASCOTTE test bench and simulations. In Section 2, a model of the test bench system is proposed. In Section 3 an active fault-tolerant control system is designed, this system is composed of an FDI part based on extended observers and adaptive cumulative sum algorithms. Section 4, describes the reconfiguration part of the Active FTCS composed of an error feedback and an Extended Unknown Input Observer (EUIO) to compensate the fault.

Description of MASCOTTE Test Bench Subsystems

The MASCOTTE test facility was developed by ONERA to study elementary processes (atomization, droplets vaporization, turbulent combustion...) which are involved in the combustion of cryogenic propellants [START_REF] Ordonneau | PLIF investigation of reactive flows in the separation region of an over-expanded twodimensional nozzle[END_REF][START_REF] Ordonneau | First results of heat transfer measurements in a new water-cooled combustor on the mascotte facility[END_REF]. Those studies in well-controlled and representative operating conditions are needed to optimize the design of high performance LPREs. For this purpose, MAS-COTTE is aimed at feeding a single element combustor with actual propellants [START_REF] Vingert | Evolution of the MASCOTTE test bench to high pressure operation and related combustor technology issues[END_REF]. Five successive versions of this test facility were built up. The thrust chamber body sub-assembly (Figure 2) consists of:

• a cylindrical section in which the combustion occurs;

155

• a section narrowing toward a throat;

• an expanding nozzle section through which the combustion gases are expelled.

This chamber is composed of three ferrules (Figure 3):

• Two heat measurement ferrules.

160

• An upstream ferrule (slightly more complex); it is equipped with the ignition torch and a larger number of thermocouples. The igniter is located in the upstream flange of the first ferrule. MASCOTTE test bench operates with oxygen (liquid or gaseous) and hydrogen or methane (gaseous) propellants. The combustion is initiated by ignition devices such as chemical pyrotechnic igniters (ignition torch). The propellants flow through the coaxial injector orifices into the thrust chamber combustion zone. The flow of liquid or gaseous propellants (see Figure 2), brought into the injection plane (see Figure 4), is calibrated by means of a cavitating Venturi.

The injection head of MASCOTTE has two modes of operation, gas / gas and liquid / gas.

MASCOTTE configurations

The test bench exists in 3 different versions (see Figure 5, and Table 1):

• Thermal measurements configuration: it consists of a cooled assembly composed of a multi-injector injection head, a modular combustion chamber and an axisymmetric nozzle. This configuration is compatible with operation at high pressure and high mixing ratio to conduct research on wall heat transfers in both the combustion chambers and nozzles (CON-FORTH) under the same conditions to those encountered in rocket engine combustion chambers. The thermal measurement chamber consists of two ferrules equipped with thermocouples. These ferrules are water-cooled structures. A set of parts, placed between the core and the body, forms the cooling part and allows the installation of thermocouples. The ground water supplies are provided by four tubes connected to a torus. Inside, a distribution grid composed of 36 holes distributes water in the cooling system.

• ATAC configuration: to study the flow detachment in a nozzle more or less over-expanded under hydrogen / oxygen combustion operating conditions representative of the conditions of a Vulcain 2 engine it is desirable to have an operating time of about 60 seconds. This objective cannot be reached with only a cooled cuff whose structure quickly reaches thermal equilibrium. For ATAC, the two-dimensional (or planar) nozzle can also be used. It consists of five main elements: three flat walls (left, right and floor), the main nozzle (convergent-divergent) and the helium throat which includes the instrumented divergent and upstream of it, the injection of parietal film simulating the re-injection of turbine gases into the Vulcain 2 nozzle extension. The nozzle is essentially equipped with wall temperature measurements. Only the whole ensemble "convergent-divergent" is concerned by this equipment. For heat flux estimation, thermocouples are used, located near the gas-side wall and the cooling-side wall. No thermocouple is placed at the throat because the local wall thickness is too small to receive thermocouples. Similarly, thermocouples directly upstream and downstream of the nozzle are located at the same thickness.

• Visualization module configuration: it has two identical flanges so that it can be turned over in order to move away the optical measurement area from the injection location. In addition, it can be mounted directly behind the injection head (need for a suitable cuff) or after a section of thermal measurements. Depending on the configuration, one or two cross-section transformation parts are required to switch from the cylindrical section to the section with the four plates. The eight cooling channels have special forms to bring water to all areas to be cooled.

MASCOTTE cooling systems

Because of the high combustion temperatures and the high heat transfer rates from the hot gases to the chamber wall, thrust chamber cooling is a ma-jor design consideration. For short duration operation (up to a few seconds), uncooled chamber walls can be used. In this case, the heat can be absorbed by the sufficiently heavy chamber wall material which acts as a heat sink, before the wall temperature rises to the failure level. For most of longer duration applications, a steady-state chamber cooling system has to be employed. The following chamber cooling techniques are used in MASCOTTE:

• Water cooling via a tubular heat exchanger: with this principle, the water is fed through passages in the thrust chamber wall for cooling and a part of it is dumped through an opening at the rear end of the nozzle skirt.

• Film cooling: here, exposed chamber wall surfaces are protected from excessive heat with a thin film of coolant (Helium) which is introduced through manifold orifices in the chamber wall near the injector and toward the throat. So, we can consider the HP air pressure sensor downstream of the regulator as part of the water circuit.

Along the water cooling system are (see Figures 3 and6):

• Six "inlet" pressure sensors: located respectively at the top of the water sphere, after a valve, at the inlet of the cuff, at the inlet of the first ferrule, at the inlet of the nozzle and on the water torus.

• Six "output" pressure sensors: located respectively at the bottom of the sphere, at the outlet of the sleeve, at the exit of the first ferrule, at the exit of the ferrule 2, at the exit of the ferrule 3 and at the neck of the nozzle.

• Three flow-meters: one at the inlet of the nozzle, one at the exit of the ferrules and one at the exit of the sleeve.

The diameters of the cavitating Venturis fixing mass flow rates were determined at the end of the development tests of the water cooling system.

The ATAC nozzle cooling part is designed as follow, see Figure 7:

• The total pressure at the outlet of the spherical tank, called "sphere" is of 39 bars,

• The part before the visualization window composed of three lines, • The part cooling the walls before the visualization window,

• The part cooling the bottom before the visualization window,

• The part after the visualization window composed of four lines,

• The part cooling the walls after the visualization window,

• The part cooling the bottom after the visualization window,

260

• The line cooling the helium throat.

MASCOTTE modeling

The Failure Modes, Effects and Criticality Analysis (FMECA) points out the necessity to perform a monitoring of the lines and cooling system (see Fig- The modeling and state estimation results were validated with the real measurements of 11 ATAC configuration trials with different mixture ratios and / or valves opening profiles. For those trials, it was assumed that the combustion chamber is long enough so that the combustion is finished and the mixture is homogeneous before entering the convergent. To simplify the implementation it has been decided to operate under a gas / gas configuration. The combustion chamber is fed with gaseous oxygen and hydrogen at ambient temperature. The models were validated by integration with measured inputs in the same initial operating conditions of the concerned MASCOTTE campaigns.

In this part we first consider a non-viscous ideal fluid system with heat exchanges. For the modeling representation, two main elements are considered with respect to the balance equations for non-viscous compressible unsteady flows:

• Cavities: defined in pressure and temperature, described by the energy conservation and mass continuity equations.

• Pipes: defined in mass flow rate, described by momentum conservation.

Continuity equation:

The total mass can be represented by the sum of the densities over the total volume and does not change over time. According to the Leibniz-Reynolds theorem and the Gauss theorem, we find an equation of mass balance (continuity equation).

dM dt = d dt V (t) ρdV (1)
∂ ∂t V (t) ρdV + S ρu.ndS = 0 (2)
The time evolution of the mass is equal to the sum of the input and output flows.

Momentum balance equation:

The Euler momentum equation is an extension of Newton's law M Γ = F ext to fluids. According to the Leibniz-Reynolds, the Gauss and Green-Ostrogradsky theorems, we find an equation for the momentum conservation.

Γ = du dt (3)
M du dt = F ext (4)
M u = V (t) ρudV (5)
∂ ∂t V (t) ρudV + S ρu(u.n)dS = F ext (6)
V (t) ∂ρu ∂t dV + S ρu(u.n)dS + S P ndS = 0 (7)
The difference between the input and output momentum over a period ∆t causes an increase in the momentum contained in the control volume. The pressure gradient creates a movement.

For a moderate turbulent flow in a smooth pipe the momentum balance equation considering friction forces is given by:

V (t) ∂ρu ∂t dV + S ρu(u.n)dS + S P ndS = -F f (8)
The friction forces can be expressed using the Blasius relation and the Darcy-Weisbach friction factor [START_REF] Nakayama | Introduction to fluid mechanics[END_REF]:

F f = λ f ρ Lu 2 2D h (9)
with

λ f := 0.316R -1 4 e (10)
Energy balance equation:

From the first law of thermodynamics, considering the total amount of energy in the entire control volume E t = V (t) ρEdV , we obtain the following equation:

d dt V (t) ρEdV = - S ρEu.ndS - dW dt + dQ dt (11)
d dt V (t) ρEdV = - S ρEu.ndS - S P u.ndS + S q.ndS (12)
For heat exchanges, it can be written, taking into account that the wall-fluid system tends towards thermal equilibrium:

d dt V (t) ρEdV = - S ρEu.ndS - S P u.ndS + S λ∇T.ndS (13)
Then the global heat transfer coefficient λ can be calculated by taking into account the thermal conduction in the walls and the convection over a heat transfer surface.

In the case of internal forced convection for short pipes with laminar flow, an initial simple approach is to utilize the dimensional analysis to obtain important parameters and dimensionless numbers. For the coolant side flow, considering a steady laminar flow of an incompressible fluid in a convectional tube. The local heat transfer coefficient can then be determined from the Nusselt number as a function of the fluid properties, geometry, temperature, and flow velocity:

N u := hL c λ ( 14
)
N u := 1.86 R e P r D L 0.33 µ µ wall 0.14 (15)
The Reynolds number R e is given by:

R e := ρD h u µ = D h ṁ µS = 4 ṁ πD h µ (16)
for a fully established flow in a circular pipe. The Prandtl number P r is defined as:

P r := µC p k (17)
The global heat transfer coefficient is given by:

λ := h 1 1 + he wall k wall S exc (18)
Here, he wall k wall is the Biot number characterizing the impact of the internal flux and external flux via the ratio of the heat transfer resistances.

Water Cooling System model

It is of prime importance to be able to detect and overcome an anomaly in this part of the engine since the chamber pressure value for a high-performance engine system is largely limited by the capacity and efficiency of the chamber cooling system. In turn, chamber pressure will affect other design parameters such as nozzle expansion area ratio, propellant feed pressure, and rockets weight (propellants consumption).

The different parts can each be modeled by two cavities defined in pressure and temperature linked by a pipe where friction forces and heat flux exchanges are taken into account, see the work of Iannetti et al. [START_REF] Iannetti | Development of model-based fault diagnosis algorithms for MASCOTTE cryo-60 genic test bench[END_REF]. For the joint CNES-ONERA ATAC research program, the two-dimensional nozzle cooling part of MASCOTTE cooling system is also modeled by a succession of cavities and pipes in parallel.

A first model with a constant mass flow rate, of the cooling system (ferrules) has been proposed by Iannetti et al. [31]. This model presented approximations in the transient assuming that the mass flow rate was constant (see Fig. 9, model 1). In the established model 2, the flow is assumed to stay single-phased, is ideal (no force due to viscosity acts) and incompressible subjects to the balance and continuity equations. The cavities volumes and the velocity of sound in the pipes are assumed constant.

Cavities equations:

We assume that the fluid flow velocity is small in comparison to the velocity of sound. The velocity in cavities is disregarded. The flow crossing cavities respects the continuity equation, after integrating (2) over the cavity volume, we obtain:

∂P ∂t = c 2 V ( ṁe -ṁs ) (19)
The temperature gradient can be determined via the energy balance [START_REF] Ioannou | Robust adaptive control[END_REF] in the cavities assuming that the temperature transfers are faster than the mechanical ones for an incompressible flow. The heat flux is written:

∆Q = h 1 1 + hew kw (T w -T av )S exc (20)
We denote ∆T := T s -T e . To obtain the water convection coefficient we use the Colburn correlation [START_REF] Kakaç | Handbook of single-phase convective heat transfer[END_REF]:

h = k D 0.023 ṁL µ 0.8 µC v k 1/3 (21)
The temperature model is given by:

∂T av ∂t = S exc θ 4 ṁ0.8 (1 + θ 1 ṁ0.8 θ 5 ) -1 ρC v V (T w -T av ) - ṁ ρV ∆T ( 22
)
with

θ 4 := k D 0.023 L µ 0.8 µCv k 1/3
, θ 5 := ew kw , and T av := 1 2 (T s + T e ).

Pipes equations:

The flow crossing the pipe between the two cavities respects the momentum balance equation with friction forces, expressed with the Darcy-Weisbach and Blasius equations for moderate turbulent flows in a smooth pipe [START_REF] Nakayama | Introduction to fluid mechanics[END_REF]. After integrating [START_REF] Hwang | A survey of fault detection, isolation, and reconfiguration methods[END_REF] along the pipe volume we obtain:

1 S 2 ∂ ṁ ∂t + ∆P V = -0.316 4 ṁ πDµ -1 4 L D h ṁ2 2ρV S 2 (23)
with ∆P := P 2 -P 1 , where 1 is for the input cavity and 2 for the output cavity. Cooling system models:

The model of the cooling system is then:

         ∂ ṁe ∂t = θ 1 ṁ 7 4 e -θ 2 ∆P ∂P2 ∂t = -θ 3 ∆ ṁ ∂Tav ∂t = Sexcθ4 ṁ0.8 (1+θ1 ṁ0.8 θ5) -1 ρCvV (T w -T av ) -ṁ ρV ∆T (24)
with ∆ ṁ := ṁ2s -ṁ2e , θ 1 := -0.316

4 πDµ -1 4 L D h 1 2ρV , θ 2 := S 2 V , θ 3 := c 2 V , 340 θ 4 := k D 0.023 L µ 0.8 µCv k 1/3
, θ 5 := ew kw , and T av := 1 2 (T s + T e ).

Parameter identification:

The parameter θ 1 must be identified since the distance L is unknown. We can assume here that the density and the viscosity remain constants for the considered pressures and temperatures ranges. One way to identify θ 1 is to use recursive least-squares by selecting one steady-state equilibrium point for the mass flow rate and the pressures. An alternative used here is the Hagen-Poiseuille formula [START_REF] Nakayama | Introduction to fluid mechanics[END_REF] in one steady-state equilibrium point for the mass flow rate and the pressures to express the unknown length as a function of the average mass flow rate ṁav : The model presented here permits to determine the pressure but also the mass flow rate and it is now possible to model their dynamics during the engine transients. The model was tested off-line with real measurements of MAS-COTTE as inputs. The final evolution of the pressure dynamics is well reconstituted (Fig. 9 and Table 2). This model has been compared to the previous one and shows a more accurate estimation of the pressure dynamics, the deviations for the first model was of 10.58% while the new one reduces the error to 5.44%.

L = - ρS 32µ ∆P ṁav D 2 (25)
The set of parameters is chosen in order to fit with the measurements and in accordance with the known properties of the test bench.

Description of GOX/GH2 Lines Downstream the Heat Exchanger

The propellants for MASCOTTE are liquid cryogenic oxygen (O2) or gaseous oxygen and gaseous hydrogen (H2). The combustion is initiated by ignition devices such as chemical pyrotechnic igniter (ignition torch). Since the use of gaseous propellants is easier to implement for studies in which the cryogenic nature of propellant is not a key point (most of the studies focus on the flow in the nozzle and not in the combustion chamber), we will model gaseous lines.

The propellant feeding lines are composed of a dome-loaded pressure regulator, a Coriolis flow-meter, a valve and cavitating Venturis to set injection mass flow rates. Pressure sensors are located upstream and downstream of the dome and the nozzle. The portion of the gaseous oxygen (GOX) / gaseous hydrogen (GH2) lines modeled is located between the outlet of the pressure dome-loaded 365 regulators and the measurement upstream of the nozzle fixing the injection rates.

Distributing lines model (pipes equation):

Using the momentum balance equation, taking into account regular pressure drops for perfect gases and assuming that the temperature remains constant along this section of the line (the speed of sound is also assumed to be constant); then after integrating along the orifice volume and assuming that the pressure evolution along the pipe is linear, we obtain:

∂ ṁ ∂t = - c 2 λL γ2DV ∆P ṁ2 ln P (L) P (0) - S L ∆P - c 2 ṁ2 γV 1 P (L) - 1 P (0) (26)
with ∆P := P (L) -P (0), where L and 0 are respectively the pressure measurements at the end and the beginning of the pipe. The model has been tested

Description of the GOX/GH2 Injections in the Combustion Chamber

The detection of failures in the injection is then an important part of the process performance and safety, a certain ratio of oxidizer (mixing ratio) to fuel in a bi-propellant combustion chamber will usually yield a maximum performance value. The propellants flow through an injector into the radial injector passages, and finally through the injector orifices into the thrust chamber combustion zone. The flow of liquid or gaseous oxygen, brought into the injection plane by Pitot tubes, is calibrated by means of a diaphragm. It is the same for the fuel (hydrogen or methane) whose distribution is then ensured by the cuff.

Injections model (cavities equation):

The flow after the cavitating Venturi of the propellant feeding lines is given by the compressible equations. The characteristic speed is assumed to be given for a nominal operation, the mixing ratio can be calculated from the flow measurements or assumed to be constant in nominal operation (these values are predetermined before a test and must remain constant in order to maintain the engine performance, see Fig. 11). The continuity equation at the injection plus the expression of the mass flow rate after the cavitating Venturi is given by:

ṁl = γP v S v,l c 2 γ + 1 γ+1 2(γ-1) (27)
The injected propellant flow rate approximated for the fuel is given by (for the oxidizer one replaces R m with 1/R m ):

ṁinj = P c,d S v,d c (R m + 1) (28)
which gives after integration, the evolution of the injection pressure over the time:

∂P inj ∂t = - c 2 V γP v S v,l c 2 γ + 1 γ+1 2(γ-1) - P c,d S v,d c (R m + 1) (29)
The average deviations for the GOX and GH2 injections are respectively of 7.603% and 12.018%. From the figure and the deviations we can see a deviation of the GOX injection pressure model from the measured output. This can be An observer is used to estimate the state of the system and to generate residuals for detection purposes. In the case of non-linear systems one of the developed techniques is to linearize and design an Extended Kalman filter (EKF) 395 used to generate the residuals (cooling system temperature, lines mass flow rates and pressures) or an EUIO as described by Witczak [START_REF] Witczak | Modelling and estimation strategies for fault diagnosis of nonlinear systems: from analytical to soft computing approaches[END_REF] in the case of model uncertainties that can be described as unknown inputs.

Models from Section 2 can be rewritten in the more general form:

   X k+1 = f (X k , U k ) + ED k + w k Y k+1 = CX k+1 + v k+1 (30)
Then they are written as a linear time-varying system with an unknown input by linearizing around a steady-state equilibrium trajectory. The system can be transformed into an equivalent discrete-time state space system as follows:

   X k+1 = A k ( X)X k + BU k + ED k + w k Y k+1 = CX k+1 + v k+1 (31)
where The first objective is to design an observer depending only on known input and output measurements. We propose to use an EUIO with the following structure [START_REF] Witczak | Modelling and estimation strategies for fault diagnosis of nonlinear systems: from analytical to soft computing approaches[END_REF]:

X k ∈ X ⊂ R n + is the state vector, Y k ⊂ Y ⊂ R m + is the measured output, U k ∈ U ⊂ R l + is the known input, D k ∈ U ⊂ R l + is
   Z k+1 = N k+1 Z k + K k+1 Y k + GU k Xk+1 = Z k+1 + HY k+1 (32)
The gain matrix is chosen to be:

K k+1 = T A k+1 P k C T (CP k C T -R k ) -1 + N k+1 H (33)
in order to minimize the variance of the state estimation error. The matrices are designed in such a way as to ensure unknown input decoupling as well as the variance of the state estimation error minimization .

e k = Xk -X k = Z k -X k + HY k ( 34
)
e k+1 = (T A k -K k+1 1 C)e k + (G k+1 -T B k )U k (35)
-(T A k -N k+1 -K k+1 C)Z k + (K k+1 2 -(T A k -K k+1 1 C)H)Y k -T ED k with K k+1 = K k+1 1 +K k+1 2 .
To reduce its expression to a homogeneous equation 415 we impose:

G = T B (36)
T A k -N k+1 -K k+1 1 C = 0 (37)
T E = 0 (38)

K k+1 2 = N k+1 H (39) with:
T = I n -HC and n the dimension of the state, N k+1 Hurwitz to ensure the asymptotic convergence of the state estimation.

A necessary condition for the existence of a solution is given under the assumption 1. A particular solution is then:

H = E((CE) T (CE)) -1 (CE) T N k+1 = T A k -K k+1 1 C (40)
In order to obtain the gain matrix K k 1 which minimizes the variance of the state estimation error, it is chosen to be:

K k+1 1 = T A k+1 P k C T (CP k C T -R k ) -1 (41)
The covariance matrix is then obtained as:

P k+1 = T A k+1 P k T A T k+1 -K k+1 1 CP k T A T k+1 + HR k+1 H T + T Q k T T ( 42
)
The residual is given by:

e k+1 = Xk+1 -CX k+1 (43)
For control purpose it is useful to dispose of all the system information. In the works of Zhu et al. [START_REF] Zhu | State estimation and unknown input reconstruction via both reduced-order and high-order sliding mode observers[END_REF] and Kalsi et al. [START_REF] Kalsi | Sliding-mode observers for systems with unknown inputs: a high-gain approach[END_REF], an auxiliary output vector is introduced so that the observer matching condition is satisfied and is used as the new system output to asymptotically estimate the system state without suffering the influence of the unknown inputs. From this result, it is possible to build an unknown input reconstruction method based on both the state and the auxiliary output derivative estimates. The auxiliary output is defined as:

Y i a,k := C i a,k X k with i = 1, .
.., p and p is the number of rows of Y k . The auxiliary output vector contains the output information of the original system.

If we denote: C a,k := C 1 ... C 1 A γ1-1 k ... C p ... C p A γp-1 k T with 1 ≤ γ i ≤ n i i = 1, ..., p
where n i is defined as the smallest integer such that:

   c i A γi k E = 0 γ i = 0, 1, ..., n i -2 c i A ni-1 k E = 0 (44)
and

C i the i th row of C then, we denote C i a,k := C i ... C i A γi-1 k T
. Since the auxiliary output vector depends on unmeasured variables, we can design a high-order sliding mode observer to get the estimates of both the auxiliary output vector and its derivative as presented by Kalsi et al. [START_REF] Kalsi | Sliding-mode observers for systems with unknown inputs: a high-gain approach[END_REF]. An estimation of the unknown input is then given by:

Dk = (M T k M k ) -1 M T k ( ξk+1 -C k (A k Xk + BU k )) (45)
with

M k := C k E, C k :=         (C 1 A γ1-1 k ) T (C 2 A γ2-1 k ) T . . . (C p A γp-1 k ) T        
, and ξk+1 :=

        C 1 A γ1+1 k Xk + C 1 A γ1-1 k BU k C 2 A γ2+1 k Xk + C 2 A γ2-1 k BU k . . . C p A γp+1 k Xk + C p A γp-1 k BU k         . Since rank(M k ) = rank(C a,k D k ) = rank(D k ) = q, M T k M k is invertible because M k
has full column rank for all the systems considered in this paper.

State Estimation and Unknown Input Reconstruction Application

The estimation cadence used on real measurements of a campaign dedicated to the determination of the combustion chamber wall temperature dispersion for different gas/gas operating points is fixed at 0.03 second, the acquisition machine acquires one point each 0.01 seconds and delivers information to the surveillance machine at a rate of one point each 0.03 seconds. The state and output noises are assumed to be Gaussian, zero-mean, white (the standard deviation is denoted σ). The state estimation error ( 43) is taken as a residual. 3 report the estimation results of the UIO for the cooling system ferrules model (the state is composed of the output pressure and input mass flow rate, the unknown input is considered to be the output mass flow rate, the known input is the input pressure), which are very satisfactory. Moreover in the case where it is not possible to measure the mass flow rates we can obtain an accurate estimate, in the permanent regime of the engine. The first peak in Fig. 4 corresponds to the transient.

To validate the unknown input reconstruction method, the results are compared to the cavity 2 output mass flow rate measurements available for these trials. Results are reported in Fig. 13. The deviation is of 17.6% and show a correct convergence after the transient phase (35.2% in the transient phase due to the linearization and 1.19e -2 % in the steady-state phase). This method can The fault detection mechanism is supposed to detect and diagnose any relevant failure and shall react sufficiently early to set up timely safe recovery.

To complete the FDIR system one needs to define residual analysis algorithms.

The objective is then to be able to detect a residual mean shift from a nomi-460 nal behavior, see [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. The observers from the previous subsection permits to estimate outputs and generate the residual defined as the state estimator error,

r k := Y k -C
Xk , in order to detect a failure and estimate its amplitude. Under some conditions, it is possible to simultaneously perform the fault estimation and unknown input decoupling in order to detect, isolate and identify faults [START_REF] Witczak | An LMI approach to robust fault estimation for a class of nonlinear systems[END_REF] [46]. In the present case, due to the model design (E = B) it is not possible, then we consider the two following hypotheses: H0: The mean value of the residual is nominal µ = µ 0 .

H1: The mean value of the residual has a shift µ = µ 1 .

For an unknown mean shift, in most common practical cases, µ 1 is unknown.

It is then possible to use an Adaptive CUSUM algorithm (ACUSUM) which estimates µ 1 as presented by Jiang et al. [START_REF] Jiang | Adaptive CUSUM procedures with EWMAbased shift estimators[END_REF]. To estimate the unknown mean shift δ, a generalization of the Exponentially Weighted Moving Average (EWMAC) chart has then be proposed allowing for a same set of parameters to improve the algorithm detection performances in the case of failures of various amplitudes and dynamics. The fault estimate is defined as:

δk = δk-1 + Φ(e p,k ) (46)
with e p,k = r k -δk-1 the prediction error, Φ is defined as a Huber score function.

Φ γ :=          e p + (1 -λ)γ , e p < -γ λe p , |e p | ≤ γ e p -(1 -λ)γ , e p > γ
with γ ≥ 0, usually fixed constant. γ is defined here at each step by γ :=| r k -δk-1 | /2 to improve the algorithm efficiency for the detection of small shifts. This leads to the following ACUSUM Statistic:

s k = ± δ± σ 2 (r k -µ 0 ± δ± 2 ) ( 47
)
where for a mean shift increase or decrease: δ+ := max (δ +,min , δk ), and δ-:= min (δ -,min , δk ). δ +,min and δ -,min are here the minimum mean shifts amplitudes to detect. The threshold is chosen to be a security coefficient times δ+ . The objective of this fault detection system is to be able to detect abrupt changes and to differentiate state perturbations and speed transients characterized by slower variations from a failure. After eliminating the effect of process input signals, filtering the effect of disturbances and model uncertainties on the residual, a residual evaluator has been designed by choosing an evaluation function and determining the threshold. To evaluate the effectiveness of the designed algorithm, the good detection (GDR) and false detection rates (FDR)

have been calculated for a simulated obstruction in the cooling system.

Fault Detection application

To choose the coefficients values and evaluate the algorithm performances, three sets of faults, composed of ten trials with different noises, have been simulated using CARINS. Each set has been simulated with various closure and opening profiles of the cooling system inflow (see Table 4, Fig. 14). The algorithm parameters are the following: δ +,min and δ -,min are fixed at ±4e -2 , the threshold security coefficient is chosen to be equal to 4.5e 4 and λ is set to 1.1055. The total time of the simulation is 60 seconds with a time step of 1 4). The cadence of the estimation and the detection is 1 time step per 30 milliseconds. The settings have been chosen to optimize the Results on Fault 2 are satisfactory since it is mandatory not to detect slow variations that can be confused with transients. Good results are obtained for Faults 1 and 3. The last case permits to evaluate the algorithm performance for successive faults of different sizes. In some rare cases the system behavior between two faults can be considered to be faulty but in most cases the two faults in 3 are well detected.

Parity Check Fault Isolation Algorithm

Once failures are detected with the ACUSUM algorithm it is necessary to be able to isolate one or several failures. The objective of this part is to isolate a fault in one or two branches (simultaneously) of the cooling system. We still consider an additive actuator failure on the system. Once the fault has been detected by an on-line and real-time first FDI mechanism the goal is to isolate the fault by a parity check (Fig. 15). The parity space-based fault detection approach is also one of the most common approaches to residual generation by using parity relations [START_REF] Chow | Analytical redundancy and the design of robust failure detection systems[END_REF]. Those relations are rearranged direct input-output model equations subjected to a linear dynamic transformation. The design freedom obtained through the transformation can be used to decouple disturbances and improve fault isolation [START_REF] Ding | A characterization of parity space and its application to robust fault detection[END_REF]. The parity space methodology using the temporal redundancy may allow to overcome time delays is to use recursion over a sliding window (see [START_REF] Wang | Fault Detection and Isolation Scheme Based on Parity Space Method for Discrete Time-Delay System[END_REF], [START_REF] Kratz | Fault detection for time-delay systems: a parity space approach[END_REF]) especially for discrete-time systems [START_REF] Naik | Recursive identification algorithms to design fault detection systems[END_REF]. In most existing works, the projection matrix for a parity check is chosen arbitrarily as described by Gao et al. [START_REF] Gao | A survey of fault diagnosis and fault-tolerant techniques-part i: Fault diagnosis with model-based and signal-based approaches[END_REF] or by Schneider et al. [START_REF] Schneider | Parity space based FDIscheme for vehicle lateral dynamics[END_REF]. The introduced approach assumes that the fault is constant in time which rules out recent works proposing new parity space approach including methods to design the projection matrix for realistic situations considering the general system with both system and measurement noises and both actuator and sensor faults, simultaneously [START_REF] Kim | A new parity space approach to fault detection for general systems[END_REF]. In our case, the fault has its own known dynamics which allow us to use direct fluid mechanics constraints.

An obstruction has been simulated on the part before the visualization window of the cooling system (surface reduction) for fault isolation. The faults have been simulated for each case on one or two different parallel lines (1, 2 or 3).

For our model of this part, we consider 3 input cavities (1, 2, 3), giving input pressures, linked by orifices (4e, 5e, 6e), giving the mass flow rates, to 3 output cavities [START_REF] Nie | Liquid-propellant rocket engine online health condition monitoring base on multi-algorithm parallel integrated decisionmaking[END_REF][START_REF] Hayton | Static and dynamic novelty detection methods for jet engine health monitoring[END_REF][START_REF] Figueroa | Rocket testing and integrated system health management[END_REF], giving the output pressures (Fig. 16).

To perform a parity check, we define the faulty system as:

   X k+1 = A k X k + BU k + ED k + F f k Y k+1 = CX k+1 (48)
The fault distribution matrix F could be different from the unknown input distribution matrix E. In this case, the projection matrix for the parity test will remain of the same form but its coefficients will change. In the studied system (the cooling system) and for the type of simulated fault (an obstruction), those

Y L,k = A L X k-L + B L U L + E L (D L + f L ) (49)
Assuming

A L := C T (CA) T . . . (CA L ) T T , B L :=         0 0 . . . 0 0 CB 0 . . . 0 0 . . . . . . . . . . . . . . . CA L-1 B CA L-2 B . . . CB 0        
, and

E L :=         0 0 . . . 0 0 CE 0 . . . 0 0 . . . . . . . . . . . . . . . CA L-1 E CA L-2 E . . . CE 0         .
The aim is to design a residual signal which is close to zero in fault free case and non-zero when a fault occurs in the monitored system. Then, for the parity check we search H L in such that:

H L (Y L -B L U L -E L D L ) = H L A L X k-L + H L E L f L = H L E L f L ( 50
)
With H L the projection matrix. The projection matrix for the parity check can then be chosen by augmenting our previous system of equations with the following relations ( 51), ( 52), ( 53), [START_REF] Schneider | Parity space based FDIscheme for vehicle lateral dynamics[END_REF]. The parallel lines have to respect the mass flow rate continuity and the energy conservation. An obstruction in a line induces an increase of the mass flow rate in the other lines and a pressure drop in a line induces a pressure increase in the other lines. The mass flow rate continuity gives:

ṁ0,k = ṁ1,k + ṁ2,k + ṁ3,k (51)
We can then use Euler conservation equations for an incompressible fluid.

P i,k+1 -P i,k = -d t c 2 V i ( ṁi,k,e -ṁi,k,s ) (52)
ṁj,k+1,e -ṁj,k,e = -d t S 2 i (P j,k -P i,k )

V i + k p d t ṁ2 j,k,e 2ρV i ( 53
)
P j,k+1 -P j,k = -d t c 2 V i ( ṁj,k,s -ṁj,k,e ) (54)
We denote ∆P q,k+1,k := P q,k+1 -P q,k for q = 1, ..., 6. This yields

ṁ0,k = ṁ4,k,e + V 1 ∆P 1,k+1,k d t c 2 + V 2 ∆P 2,k+1,k d t c 2 + V 3 ∆P 3,k+1,k d t c 2 + ṁ5,k,e + ṁ6,k,e
The detection algorithm is then triggered after the transient to not consider With the help of those expressions we can then find the projection matrices.

We have:

Y k+1 -CBU k -CED k = CEf k + CA k X k (55)
Since CB = 0, we have:

Y k+1 -CBU k -CED k = Y k+1 -CED k ( 56
)
CED k = -c 2 dt V1 ṁ4,k,s -c 2 dt V2 ṁ5,k,s -c 2 dt V3 ṁ6,k,s T ( 57
)
and ṁj,k,s = ṁj,k,e -

Vi(∆P j,k+1,k ) c 2 dt
for i = 1, ..., 3, j = 4, ..., 6. Then:

Y k+1 -CED k =                     c 2 dt V1 ( ṁ0,k -ṁ6,k,e -ṁ5,k,e ) -∆P 1,k+1,k - V2(∆P 2,k+1,k ) V1 - V3(∆P 3,k+1,k ) V1 + P 4,k c 2 dt V2 ( ṁ0,k -ṁ4,k,e -ṁ6,k,e ) - V1(∆P 1,k+1,k ) V2 -∆P 2,k+1,k - V3(∆P 3,k+1,k ) V2 + P 5,k c 2 dt V3 ( ṁ0,k -ṁ5,k,e -ṁ4,k,e ) - V1(∆P 1,k+1,k ) V3 - V2(∆P 2,k+1,k ) V3 -∆P 3,k+1,k + P 6,k                     (58)
with: ṁj,k,e = The projection matrix H has to verify:

HCA k X k = 0 (59)
Using ( 58), H is then equal to:

H :=      h 1 h 2 h 3 h 1 h 2 h 3 h 1 h 2 h 3      with: h i := 3ω k 3 d t c 2 V i ṁj,k,e +3P j,k -ω k ,i = 1...3, j = 4...6,

545

and

ω k := ( ṁ0,k - V1∆P 1,k+1,k dtc 2 - V2∆P 2,k+1,k dtc 2 - V3∆P 3,k+1,k dtc 2
). Since for i = 1, ..., 3, j = 4, ..., 6 we have:

P j,k = P j,k+1 + d t c 2 V i ( ṁj,k,s -ṁj,k,e ) (60)
P j,k = Y i,k+1 -(CED k ) i - d t c 2 V i ṁj,k,e (61)
Then:

H L :=         H 0 ... 0 0 H ... 0 ... ... ... ... 0 0 ... H        
. The estimate of faults f L is then obtained from (see Fig. 17 for results on the example):

f L = (H L E L ) -1 H L (Y L -E L D L ) (62)
To isolate the fault we can compare the variation of faults. • If the variation is negative for two pipes (↓/↑) then the fault occurs in the other pipe (+/-, one fault case): an obstruction implies a mass flow rate decrease in the impacted line so that the mass flow rate continuity for the overall system allow us to conclude that the mass flow rate increases in the other lines.

• As long as the sign of variations remain the same, faults are persisting

(+/-).
This analysis is summarized in Tables 5 and6. The termsmax , + max , min , + min indicates if a fault in a line is of greater amplitude than the fault in the other line. The arrows indicate if the mean values of the residuals increase or decrease. fault in line 1 fault in line 2 fault in line 3

r 1 +/- ↓/↑ ↓/↑ r 2 ↓/↑ +/- ↓/↑ r 3 ↓/↑ ↓/↑ +/-
To evaluate the effectiveness of the designed algorithm, the good detection and false detection rates (GDR, FDR) have been calculated for ten runs. For simultaneous faults we consider to be a good detection the simultaneous detection and isolation of the faults in the two impacted lines, if at least one detection is false then we consider it to be a false detection. Those rates, which are satisfying for the considered application, (Table 7) have been calculated from ten runs for each simulation and the settings have been chosen to optimize the good isolation rate and minimize the false isolation rate of abrupt mean shifts.

Reconfiguration Algorithms for Non-Shutdown Actions

Once the fault has been detected and isolated by an on-line and real-time FDI mechanism, in the case of non-shutdown actions, the goal is to maintain the overall system stability and an acceptable performance in spite of the occurrence of faults and saturations by reconfiguring the nominal control law.

Residuals

Case 4 fault in lines 1 and 2

r 1 + - -max + max + min -min r 2 + - + min -max -max + min r 3 ↓ ↑ / ↓ ↑ /

Residuals

Case 5 fault in lines 1 and 3

r 1 + - -max + max + min -min r 2 ↓ ↑ / ↓ ↑ / r 3 + - + min -max -max + min

Residuals

Case 6 fault in lines 2 and 3

r 1 ↓ ↑ / ↓ ↑ / r 2 + - -max + max + min -min r 3 + - + min -max -max + min

Design of a Closed-Loop Fault-Tolerant Control System

For systems such as MASCOTTE ones, which relies on main propellant flow variation, the closed-loop control system has to operate on the principle of variable fluid resistances (pressure dome-loaded regulators) in the main oxidizer and fuel feed lines to achieve propellant flow-rate modulation or in the cooling system lines to overcome performances losses. The principal reasons for mixtureratio control are recalled:

• Optimum engine performance (important)

• Complete propellant utilization; i. e., minimum residuals (most important)

Based on the FMECA, in a first approach one can see that an obstruction or a leakage in the propellant manifolds may be critical and imply a shutdown action. For that reason we will validate our active fault-tolerant control algorithm by considering in this part an additive actuator failure on the system, which may correspond to an obstruction or a leakage. We also study the possibility of a reconfiguration of the propellants mass flow rates in order to maintain a suitable mixing ratio. The method proposed here consists in the design of a controller based on an unknown input observer by considering the fault to be the unknown input similar to the work of Hamayun et al. [START_REF] Hamayun | A fault tolerant control allocation scheme with output integral sliding modes[END_REF] and the design of an anti-windup strategy, extending the results of Da Silva et al. [START_REF] Da Silva Jr | Anti-windup design with guaranteed regions of stability for discrete-time linear systems[END_REF], in order to ensure the asymptotic stability of the saturated system for a given set of initial conditions and determine the stability domain. This FTC strategy permits to compensate the fault and maintain current performances in the presence of actuator saturations but also to converge if necessary to another state reference.

The system considered is the following:

(Σ)    X k+1 = A c X k + B c U k + B c f k Y k+1 = CX k+1 (63)
where X k ∈ n is the state vector, Y k ∈ m the measured output, U k ∈ l the known input, f k ∈ l the unknown actuator failure, A c ∈ n×n the state matrix, B c ∈ n×l the known input distribution matrix and C ∈ m×n the output distribution matrix, with m ≤ n. In the case of models with unknown inputs, the linear quadratic regulator is defined from an equivalent system where the unknown input is expressed as a function of the known state and known input vectors in order to decouple the fault effect on the system from the unknown input dynamics (see Appendix C).

In the unsaturated case, we can use a control law of the form:

U k := -B + c B c fk + W c ( Xc,k -X k )
, where we denote B + c as the pseudo-inverse of B c , -B + c B c fk is the fault compensation part and W c ( Xc,k -X k ) is the reconfiguration part. The fault magnitude estimate fk is estimated with an EUIO and the gain W c is calculated with a LQR approach [START_REF] Theilliol | Actuator fault tolerant control design based on a reconfigurable reference input[END_REF]. For a global state configuration the computational burden might be too high to calculate the gain at each time step. To overcome this issue, one can use the result on polytopes in Appendix A. This result gives a global gain matrix based on a Lyapunov stability demonstration, considering that the matrices A k are bounded and belong to a polytopic set (see [START_REF] Ramos | A less conservative LMI condition for the robust stability of discrete-time uncertain systems[END_REF], [START_REF] Ramos | An LMI approach to compute robust stability domains for uncertain linear systems[END_REF]). The state reference trajectory X k is predetermined and its dynamics is given by:

X k+1 = A c Xk + B c Ūk (64)
with Ūk the nominal input. Since the fault tolerant control is activated once a fault has been detected, the nominal input can be chosen as the mean input over a sliding window during nominal performances.

The control law can be alternatively written as U k := -B + c B c fk +W c e c,k +W c η k . The dynamics of the augmented state is expressed as: stability of the saturated system as inputs. Unlike conventional anti-windup methods based on the resolution of bilinear matrix inequalities, this method is relatively simple and proposes an iterative algorithm of linear matrix inequalities in the same spirit as presented by Da Silva et al. [START_REF] Da Silva Jr | Anti-windup design with guaranteed regions of stability for discrete-time linear systems[END_REF]. In this approach, the set 640 of admissible initial states and its associated domain of stability are determined to take into account the compensation of additive actuator faults.

ζ k+1 =   A c + B c W c B c W c 0 N c   ζ k (
The aim of this part is to design an anti-windup control law in order to ensure the asymptotic stability of the system with a saturated input for a given set of initial conditions and determine the associated stability domain. We want to determine the anti-windup gain matrix E c such that for a set S of admissible initial states (ζ 0 ∈ S), the corresponding trajectory converges asymptotically to the origin of the subset E ⊂ S. Then, E is a region of asymptotic stability.

For that, we want to determine a new control law of the form

U k + = U k -Gζ k
when the control law U k reaches its bounds with G ∈ l×2n . The reference state dynamics for the anti-windup strategy is given as:

X k+1 := A c X k + B c U k + E c (sat(U k ) -U k ), and
U k := U k -B + c B c fk + W c ( Xc,k -X k ). If the control law is saturated then U k = ±U sat : X k+1 = A c X k + B c U k + E c (±U sat -U k + B + c B c fk -W c e c,k -W c η k ) (67)
Which gives

η k+1 = A c η k + B c W c e c,k + B c W c η k + E c (±U sat -U k + B + c B c fk -W c e c,k -W c η k ) (68)
we then have:

ζ k+1 =   A c + B c W c B c W c 0 N c   ζ k -(RE c )Ψ(Kζ k ) (69) Ψ(u) :=          u i -U i,sat + U i,k -(B + c B c fk ) i , if u i + U i,k -(B + c B c fk ) i > U i,sat 0, if -U i,sat ≤ u i + U i,k -(B + c B c fk ) i ≤ U i,sat u i + U i,sat + U i,k -(B + c B c fk ) i , if u i + U i,k -(B + c B c fk ) i < -U i,sat          where R =   I n 0   , K = W c W c , ∀i = 1, . . . , l.
The set of admissible initial states S considered will be defined as a polyhedral set and the domain of stability E will be designed as an ellipsoid (see proofs in Appendix D).

Determination of the set of admissible initial states

Lemma 1. Consider a matrix G ∈ l×2n and define the following polyhedral set:

S = [ζ k ∈ 2n ; -U i,sat ≤ ((K -G)ζ k ) i + U i,k -(B + c B c fk ) i ≤ U i,sat ; ∀i = 1, . . . , l] (70)
For the function Ψ(u) defined in (70), if ζ k ∈ S then:

Ψ(Kζ k ) T T [Ψ(Kζ k ) -Gζ k ] ≤ 0 (71)
for any matrix T ∈ l×l diagonal and positive definite.

This property will be used in the proof of Theorem 2 (D. [START_REF] Theilliol | Actuator fault tolerant control design based on a reconfigurable reference input[END_REF]) to find the gain E c depending on the choice of G to ensure the exponential asymptotic stability of the system.

E(P ) ⊂ S.        W 0 2n,1 0 1,2n -1   WK T i -(GW)i T (B + c Bc fk )i-U i,k KiW-(GW)i (B + c Bc fk )i-U i,k Ui,sat 2 (B + c Bc fk )i-U i,k 2      ≥ 0 (72) ∀i = 1, . . . , l
Remark 1. The case where fi = 0 is not considered since the reconfiguration is triggered only in the case of a fault detection.

Remark 2. The case where U i,k = 0 is not considered since in our application case the input reference is a pressure consign. In other cases, a state feedback as presented by Da Silva et al. [START_REF] Da Silva Jr | Anti-windup design with guaranteed regions of stability for discrete-time linear systems[END_REF] might be used due to the reference state dynamics form.

Remark 3. We assume that (B + c B c fk ) i -U i,k = 0 since the case where the fault amplitude is equal to the reference input is considered as a critical failure inducing the system shutoff.

Determination of the associated domain of stability

In this part, we denote:

A :=   A c + B c W c B c W c 0 N c   . Z ∈ n×l and
∆ ∈ l×l a diagonal positive definite matrix are parameters which will be chosen in order to maximize the size of the set of admissible initial states and ensure the exponential asymptotic stability of the augmented system (69).

Theorem 2.

The ellipse

E(P ) = [ζ k ∈ 2n , ∀i = 1, . . . , l; ζ k T P ζ k ≤ 2] with P = W -1 is
a region of exponential asymptotic stability for the augmented system, if for

E c = Z∆ -1 :      W -(GW) T -WA T -(GW) 2∆ Z T R T -AW RZ W      > 0 (73)
for the considered Lyapunov candidate quadratic function:

V (ζ k ) := ζ k T P ζ k , P = P T > 0, P ∈ 2n×2n V (ζ k ) is a Lyapunov function since: 1. δV (ζ k ) < 0, ∀ζ k ∈ E(P ), ζ k = 0 2. ∃α ∈ + , δV (ζ k ) ≤ -αV (ζ k ) 4.3.3. Application
When a fault is detected the system switches to the FTCS in the case of the cooling system. In the case of the propellant injection control, the system switches to the FTCS at a prefixed switch-time. The desired transient behavior depends on the gain choice; we have to limit the overshoots to maintain the system performances. The aim of those simulations is to see if the controller is able to stabilize the closed-loop system after the detection or when the switchtime is imposed, see (Tables 8,9, and Fig. 18, and 21). The water cooling system is regulated with a pressure dome-loaded regulator (sphere) and valves.

The actuator is saturated since the pressure is limited by thermo-mechanical

constraints. An obstruction at the input of the ferrules part has been simulated by computing a closure profile of the valves. The closure profile is computed as a modification of the cross-sectional area of the actuator. The faults were implemented as in the previous section. The simulated cases for the cooling system are the following:

• Case 1: constant valve closing profile, no actuator saturations.

• Case 2, 3, and 4: time varying closing profile (successive faults of different magnitudes) no actuator saturations.

• Case 5: constant valve closing profile, with actuator saturations and a new state reference.

For the first four trials, faults are compensated and it can be seen that the control law for the rewritten system permits to stabilize the system around the reference steady-state equilibrium with sufficient accuracy. The FTCS with anti-windup trial (Case 5) aims at compensating the fault and at converging to This trial shows that the fault is well compensated and the convergence to the nominal value is faster than in the case of a controller with a fixed limit value.

We can also see that since the reference state dynamics is modified by the antiwindup scheme in order to ensure the exponentially asymptotic convergence, the trajectory is more stable in this case than in the case of an fixed imposed limit (Figures 19 and20). The new reference state dynamics is consistent with the established model; we can see that the dynamics relations between the state and the input are respected.

The simulated cases for the propellants injection are the following: In this case, the performances (Table 9 and Fig. 21) are satisfying however, 710 the simulations consider isolated systems and do not take into account the impact on the combustion chamber pressure which is linked to the mass flow rate injection. Future work will address its modeling. The closed-loop performances for the cooling system are lower than the performances for the propellants injection regulation due to the fault compensation error; however, the results are satisfying. Table 8 shows that even in the case of successive faults, once the system has switched to closed-loop FTCS the performances are maintained. It is then not required to switch back to the open-loop system once the stability around a nominal value is obtained.

Conclusion

A new model for the purposes of FDIR was proposed for the evolution of pressure, temperature and mass flow rates in the cooling system of a cryogenic test bench and for the evolution of pressure and mass flow rates in the propellants feeding lines and injection for gas/gas operations. An Active FTCS based on a FDI system, a fault estimator and an EUIO permits to compensate an actuator failure and to converge if necessary to a chosen steady state.

The FDI method consists in a detection part with an observer-based residual generation, using an EKF or an EUIO. Residuals are analyzed with the means of an ACUSUM algorithm. Then in the case of the cooling system where the lines are interdependent, a parity-space based fault isolation method has been proposed to isolate faults, using a projection matrix defined by dynamics fluid mechanic relations for the overall system. The Active FTCS is based on a second EUIO where the unknown input is then considered to be the fault and a linear quadratic regulator. A method to design an anti-windup scheme has also been proposed in order to compute another steady point which may be reachable in the case where the previous nominal steady point cannot be reached because of the actuator failure and the effect of the saturation. This method is based on the resolution of linear matrix inequalities and ensures exponential asymptotic stability in an ellipsoidal domain for a polyhedral set of admissible initial states.

This method was tested on a model proposed for the evolution of pressure and mass flow rates in the cooling system, the propellants feeding lines and injection of MASCOTTE test bench. Future work will address the design of a method to improve the performances during transients by fixing predetermined nominal trajectories and the extension of this work to non-linear systems. The rest of the work will consist in validating the performance of FDI on-line, in order to evaluate the quality of the global detection and reconfiguration loop. for T diagonal and positive definite.

V
(2) We consider the case where: (Kζ k ) i + U i,k -(B + c B c fk ) i < U i,sat then,

(Kζ k ) i + U i,sat + U i,k -(B + c B c fk ) i < 0 (D.4)
We have:

Ψ(Kζ k ) = (Kζ k ) i + U i,sat + U i,k -(B + c B c fk ) i (D.5) If ζ k ∈ S, ((K -G)ζ k ) i + U i,k -(B + c B c fk ) i ≥ -U i,sat
, then:

[(Kζ k ) i + U i,sat + U i,k -(B + c B c fk ) i ] T T i,i [((K -G)ζ k ) i + U i,sat + U i,k -(B + c B c fk ) i ] ≤ 0 (D.6)
for T diagonal and positive definite. Then we have:

  P 0 2n,1 0 1,2n -1   -   K T i -G T i Yi -1   U i,sat -2 Y i -2
Ki-Gi Yi

  ζ k Yi Yi   T   K T i -G T i Yi -1   U i,sat -2 Y i -2 Ki-Gi Yi -1   ζ k Yi Yi   then ζ k T P ζ k - Y 2 i Y i 2 ≥ ζ k T K T i -G T i Yi -Yi Yi U i,sat -2 Y i -2 Ki-Gi Yi ζ k -Yi Yi (D.11) So ζ k ∈ S since ζ k T P ζ k - Y 2 i
Yi 2 ≤ 1:

ζ k T P ζ k - Y 2 i Y i 2 U i,sat 2
Y i 2 ≥ ζ k T K T i -G T i Yi -Yi Yi Ki-Gi Yi ζ k -Yi Yi (D.12) U i,sat 2
Y i 2 ≥ ζ k T K T i -G T i Yi -Yi Yi Ki-Gi Yi ζ k -Yi Yi (D.13)
We then have: We can write this inequality under the form: Then we have δV (ζ k ) < 0 for all ζ k ∈ E(P ), ζ k = 0 , so V (ζ k ) is strictly decreasing along the system trajectories. Then E(P ) is a stability region for the system. We can see that there always exists a positive scalar δ such that:

-U i,sat ≤ ((K i -G i )ζ k + U i,k -(B + c B c fk ) i ) ≤ U i,sat (
δV (ζ k ) ≤ -ζ k T Ψ T   X 1 X 2 X 2 T X 3     ζ Ψ   (D.
δV (ζ k ) ≤ -δ ζ k 2 -δ Ψ 2 ≤ -δ ζ k 2 ≤ -δζ k T P ζ k (D.20)
which ensures the exponential convergence with δ := δ λmax(P ) and λ max (P ) the maximum eigenvalue of P .

µ

  Dynamic viscosity, P a.s ρ Density, kg/m 3

Figure 1 :

 1 Figure 1: Diagram of the complete loop

Figure 2 :

 2 Figure 2: MASCOTTE test bench -Simplified Synoptic -gas / gas operations

Figure 3 :

 3 Figure 3: MASCOTTE test bench -Ferrules

Figure 4 :

 4 Figure 4: MASCOTTE test bench -Injector

Figure 5 :

 5 Figure 5: MASCOTTE test bench -Configurations

Figure 6 :

 6 Figure 6: MASCOTTE test bench -Water cooling circuit

Figure 7 :

 7 Figure 7: MASCOTTE test bench -Cooling system -ATAC + visualization configuration

ure 8 )

 8 pressures or mass flow rates and temperatures, as well as the injection 265 pressure-drops. Such failures can result in critical events or impact the engine performances.

Figure 8 :

 8 Figure 8: Flow diagram of MASCOTTE test bench -Propellant distributing manifolds and cooling system models

Figure 9 :

 9 Figure 9: Pressure model -Ferrules

Figure 10 :

 10 Figure 10: GH2 -Mass flow rate model

Figure 11 : 3 .

 113 Figure 11: GOX injection pressure model

  the unknown input vector, and X the equilibrium state. The considered mass flow rates and pressures are bounded by thermomechanical constraints, X∈ X = [X inf ; X sup ], Y ∈ Y = [Y inf ; Y sup ] and U ∈ U = [U inf ; U sup ].With A ∈ R n×n the state matrix, B ∈ R n×l the known input distribution matrix, E ∈ R n×l the unknown input distribution matrix and C ∈ R m×n the output distribution matrix, with m ≤ n, w k and v k are respectively the state noise and the measurement bounded Gaussian white noise of covariance matrices Q k and R k . Assumption 1. rank(CE) = rank(E) = rank(D k )

Figure 12 :

 12 Figure 12: Pressure residual -Ferrules

Figure 12 ,

 12 Figure 12, and Table 3 report the estimation results of the UIO for the

Figure 13 :

 13 Figure 13: Unknown input reconstruction

Figure 14 :

 14 Figure 14: Fault 3 residual -CARINS simulation

Figure 15 :

 15 Figure 15: Nozzle cooling system FDI

Figure 16 :

 16 Figure 16: Nozzle cooling system synoptic -upstream part

540

  them as failures in a first time. A failure is assumed to impact proportionally the mass flow rate: ṁj,k,e := (f r,i,k + 1) ṁj,k,e,nominal or again:ṁj,k,e := (f r,i,k + 1) 2S 2 (∆P nominal ) kp -2ρV (∆ ṁe,nominal ) kpdt .We obtain the expression of faults in each line f r,i,k in the case of a single fault and two simultaneous faults (full expressions are given in Appendix B).

2S 2 i

 2 ρ(P j,k -P i,k ) kp -2ρVi( ṁj,k+1,e -ṁj,k,e ) kpdt for i = 1, ..., 3, j = 4, ..., 6.

Figure 17 :

 17 Figure 17: Case 1 reconstructed faults

  65) where ζ k := η k e c,k T , with e c,k = Xc,k -X k the estimation error, η k = X k -X k the reconfiguration error and X k the state reference. N c is the gain of an observer ensuring the estimation error convergence so that its dynamics reduces to e c,k+1 = N c e c,k .

Theorem 1 .

 1 Define E(P ) = [ζ k ∈ 2n , ∀i = 1, . . . , l; ζ k T P ζ k ≤ 2] with P ∈ 2n×2n a positive definite matrix and W := P -1 . If W satisfies (72) for each input value then

Figure 18 :

 18 Figure 18: Cooling system pressure and mass flow rate control -CARINS simulator -Case 1

Figure 19 :

 19 Figure 19: Fault compensation & reconfiguration -Case 5 -Control law & system input

Figure 20 :

 20 Figure 20: Fault compensation & reconfiguration -Case 5 -System state

Figure 21 :

 21 Figure 21: GH2 mass flow rate control and estimation

( 3 ) 8 )

 38 Ψ(Kζ k ) = 0, then:Ψ(Kζ k ) T T [Ψ(Kζ k ) -Gζ k ] = 0 (D.7)for T diagonal and positive definite.Proof of Theorem 1:770 By Schur's complement, (72) gives ∀i = 1, . . . , l:With P = W -1 , Y i = (B + c B c fk ) i -U i,k, K i and G i are the i th lines of K and G.

D. 14 )Proof of Theorem 2 :

 142 so that E(P ) ⊂ S775We calculate ∆V (ζ k ):δV (ζ k ) = V (ζ k+1 ) -V (ζ k ) (D.15) = ζ k T A T P Aζ k -2ζ k T A T P (RE c )Ψ(Kζ) + Ψ(Kζ k ) T (RE c ) T P (RE c )Ψ(Kζ k ) -ζ k T P ζ kUsing Lemma 1, we have:δV (ζ k ) ≤ -(ζ k T A T P Aζ k + 2ζ k T A T P (RE c )Ψ(Kζ k ) -Ψ(Kζ k ) T (RE c ) T P (RE c )Ψ(Kζ k ) + ζ k T P ζ k ) -2Ψ(Kζ k ) T T [Ψ(Kζ k ) -Gζ k ] (D.[START_REF] Theilliol | Actuator fault tolerant control design based on a reconfigurable reference input[END_REF]

  )with:X 1 := P -A T P A, X 2 := A T P (RE c )-G T T , X 3 := 2T -(RE c ) T P (RE c ).By Schur's complement, (73) gives: T := ∆ -1 and P := W -1 we have:

Table 1 :

 1 MASCOTTE configurations

		Thermal measurements	ATAC
	Composition	Two water-cooled ferrules	Two dimensional nozzle
		Axisymetric nozzle
	Operations	High pressure	60 second operating time
		High mixing ratio	Vulcain 2 representative
	Study	Heat transfers in combustion	Flow detachment in an
		chambers and nozzles	over-expanded nozzle (O2/H2)
		Visualization
	Composition	Visualisation window and 8 cooling channels
	Operations	Mountable behind injection head or thermal measurements part

Table 2 :

 2 Deviations of the pressure model 1 and 2 -Ferrules

	Model	Total	Transient Permanent
		(%)	(%)	(%)
	Pressure (1)	10.58	13.35	5.04
	Pressure (2)	5.44	8.01	0.31
	Input mass flow rate (2) 3.31e-5	4.97e-5	6.17e-8

Table 3 :

 3 Deviations of the pressure and input mass flow rate estimations -Ferrules

	Model	Total	Transient Permanent
		(%)	(%)	(%)
	Pressure, Pa	9.92e-2	7.00e-2	1.16e-2
	Mass flow rate, kg/s	6.27	31.4	1.18e-2

also be useful in the case of Vulcain 2 engine during an Ariane flight, where it is difficult or expensive to measure the mass flow rate.

Table 4 :

 4 Simulated faults

	Fault	Type	GDR F DR N begin N end
			(%)	(%)
	Fault 1	Abrupt	98.8	0.0	1367	1540
		large mean shift
	Fault 2	Slow	27.4	0.0	1032	1252
		large mean shift
	Fault 3	Slow			1310	1368
	(1)	small mean shift
		abrupt recovery	98.5	14.5
	Fault 3	Abrupt			1532	2000
	(2)	large mean shift
		slow recovery
	millisecond (Table

Table 5 :

 5 Residuals variations -1 failure cases

	Failure	Case 1	Case 2	Case 3

Table 6 :

 6 Residuals variations -2 failures cases

Table 7 :

 7 Isolation rates

		Faults Fault 1 Fault 2 Fault 3 Fault 4 Fault 5 Fault 6
		GDR 96.06% 91.07% 90.40% 80.29% 93.06% 90.21%
		FDR	0.00%	4.51%	4.51%	4.51%	3.75%	3.75%
		4.1. MASCOTTE Test Bench Open-Loop Control System
		With this system, control is accomplished by preset control means, such as
	585	orifices, and on off command devices as it is currently done for most existing
		rocket engine systems. The extent of correction is determined from systems
		pretest calibration test data. Open-loop control has the advantage of simplicity,
		however, it is limited to a specific set of operating parameters, and is unable to
		compensate for variable conditions during operation.

Table 8 :

 8 Deviations (%) of the simulated cooling system pressure and input mass flow rate

	from references
	Model (2) with control	Case 1 Case 2 Case 3 Case 4
	Pressure, Pa	3.02	2.67	2.67	3.94
	Input mass flow rate, kg/s	2.57	1.91	2.14	2.75

a different reference state than the nominal one (chosen arbitrarily). We fixed the saturated value at U sat = 3.782 • 10 6 P a in this case, the saturation value has been chosen in order to allow the convergence to the new state reference.

Table 9 :

 9 Deviations (%) of the simulated GH2 injection pressure and mass flow rate from

	references
	Model with control Propellant Case 1	Case 2
	Pressure, Pa	GH2	1.39e -4 2.68e -6
		GOX	7.31e -4 6.84e -4
	Mass flow rate, kg/s	GH2	9.47e -3 2.375e -3
		GOX	2.28e -2 1.18e -3

  2 ∆P 5,2,k + ( ṁ5,k+1,e -ṁ5,k,e )) -

								c 2 V 3 √ ρ	(P 3,k+1 -P 3,k )
	-	2V 3 k p d t	(d t	S 2 3 V 3	∆P 6,3,k + ( ṁ6,k+1,e -ṁ6,k,e ))
		2S 2 1 ∆P 4,1,n,k kp	+	1 2V1( ṁ4,k+1,e,n -ṁ4,k,e,n ) kpdt	-1	(B.1)
	f r,2,k =		ṁ0,k √ ρ	-	c 2 V 2 √ ρ	(P 2,k+1 -P 2,k ) -	c 2 V 1 √ ρ	(P 1,k+1 -P 1,k )
	-(f r,1 + 1)			2S 2 1 ∆P 4,1,n,k k p	+	2V 1 ( ṁ4,k+1,e,n -ṁ4,k,e,n ) k p d t
	-	2V 3 k p d t	(d t	S 2 3 V 3	∆P 6,3,k + ( ṁ6,k+1,e -ṁ6,k,e )) -	c 2 V 3 √ ρ	(P 3,k+1 -P 3,k )
	2S 2 2 ∆P 5,2,n,k kp	+	1 2V2( ṁ5,k+1,n,e -ṁ5,k,e,n ) kpdt	-1	(B.2)
	f r,3,k =		ṁ0,k √ ρ	-	c 2 V 2 √ ρ	(P 2,k+1 -P 2,k ) -	c 2 V 1 √ ρ	(P 1,k+1 -P 1,k )
	-(f r,1 + 1)			2S 2 1 ∆P 4,1,n,k k p	+	2V 1 ( ṁ4,k+1,e,n -ṁ4,k,e,n ) k p d t
	-	2V 2 k p d t	(d t	S 2 2 V 2	∆P 5,2,k + ( ṁ5,k+1,e -ṁ5,k,e )) -	c 2 V 3 √ ρ	(P 3,k+1 -P 3,k )
	2S 2 3 ∆P 6,3,n,k kp	+	1 2V3( ṁ6,k+1,e,n -ṁ6,k,e,n ) kpdt	-1	(B.3)

For the nominal system, the gain W c must stabilize (A c + B c W c ). Since the pair (A c , B c ) is assumed to be controllable, a LQR formulation can be adopted where W c is selected to minimize J

where Q c,k and R c,k are symmetric positive definite design matrices. When the input is assumed to be saturated the system considered becomes:

where U k ∈ l is the control law and U sat ∈ l + is the actuator limit, ∀i = 1, . . . , l.

Anti-windup control law

Multiple solutions have been studied to compensate for a decrease in system performance caused by the saturation of one or more actuators, one way is to add a so-called anti-windup command, another way is to use direct synthesis methods by taking into account the saturations in the control law. The idea of the anti-windup approach is to add a state, output or error feedback so that the actuator remains within its limits. This consists in neglecting the saturation in the first stage of the control design process, and then to add some problemspecific schemes to deal with the adverse effects caused by saturation. In the case of discrete systems, our interest is the development of control laws that provide a semi-global convergence on any arbitrarily large set of state space.

They usually have a simpler structure and the controller is less sensitive to model and disturbance uncertainties. The system performance one wants to achieve can range from the classic system stabilization problem to expanding the area of attraction, rejecting disturbances, and regulating the output of the system [START_REF] Hu | Control systems with actuator saturation: analysis and design[END_REF]. The advantage of the presented control method is that it studies the determination of the stability regions of a discrete linear system over time and allows determining an anti-windup control law which ensures the asymptotic Appendix A. Gain determination with polytopic sets

In the case where the full state is of greater size it may not be possible to solve the Riccati equation in real time. In this case it is possible to use results on polytopes in order to compute a global gain. We consider a polytope A(α)

A is bounded and closed with Ā the upper limit, A the lower limit and

If all A can be written as a convex combination of Ā and A we want to show that the Lyapunov function P assuring the Lyapunov stability of the system is the same convex combination of the limits corresponding Lyapunov functions. We denote P the Lyapunov function so that we have:

ĀT P + P Ā + P BR -1 B T P < -Q and P the Lyapunov function so that we have:

where Q and R are symmetric positive definite matrices. We calculate: A(α)

with A(α) = ᾱ Ā + αA and P (α) = ᾱ P + αP . We obtain:

if we assume that (A T P + P A+ P BR -1

= ᾱ2 ( ĀT P + P Ā + P BR -1 B T P ) + α 2 (A T P + P A + P BR -1 B T P ) + α ᾱ(A T P + P A + P BR -1 B T P + ĀT P + P Ā + P BR -1 B T P )

Appendix B. Faults dynamics calculations

In the case of an obstruction in the line 1:

In the case of 2 failures, for example in lines 1 and 2:

Cases

Since the system is controllable in pressure and a pressure failure induces a mass flow rate failure, it is necessary to be able to decouple the fault effect on the system dynamics from the unknown input one. The first step is then to decouple the reconstructed cavity 2 output mass flow rate from the system 760 dynamics in order to evaluate the actuator fault and compensate it.

Now that the unknown input expression is available (see the complete loop

Fig. 1), we can rewrite the first system without the mass flow rate as an unknown input when it is necessary. Then, in order to annihilate the actuator fault effect on the system, another EUIO is used to estimate the fault magnitude. A control law has then to compensate the fault and be computed such that the faulty system is as close as possible to the nominal one. We can use the previous result from the unknown input reconstruction part to rewrite the system under a second form only depending on known inputs for control purposes:

The system is linearized around a steady-state equilibrium, the nominal state to reach. This method requires matrix inversions, which may be numerically unstable due to possible ill-conditioning. In the problem considered, the matrices were invertible.

765

We obtain a new system under the form:

with A c the new state matrix and B c the new known input distribution matrix:

Appendix D. FTCS proofs Proof of Lemma 1:

(1) We consider the case where:

We have:

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐ The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

[Grab your reader's attention with a great quote from the document or use this space to emphasize a key point. To place this text box anywhere on the page, just drag it.]

x


Dominique Tournier

Philippe Godignon

José Millan

Dominique Planson

Jean Pierre Chante

Franck Sarrus

Jean François De Palma

On-Chip Temperature Monitoring of a SiC Current Limiter

Keywords: Current limiter, temperature sensor, reliability, packaging

. Carrier mobility dependence with temperature was extracted from cryogenic measurements. The temperature estimation is based on the measurement of the variation of the electrical resistance (caused by mobility variation) of the sensing device integrated with the current limiter. In this paper we will first describe the temperature estimation methodology using various technological solution (from metallic resistor solution to the SiC integrated sensor). Then experimental temperature measurements using an integrated SiC sensor within a packaged current limiting devices will be presented. Electro-thermal measurements on the fabricated devices show that the current limiter is able to work at 205°C under steady state conditions (320 V), without degrading their electrical performances. Finally, perspectives in terms of integration and reliability will be proposed.

Introduction

In order to get the full benefit of silicon carbide material properties; i. e. high voltage and high temperature working operation [START_REF] Matsunami | Progess in Wide Bandgap Semiconductor SiC for Power Devices[END_REF], the development of specific package is a great challenge on the roadmap to silicon carbide technology. Therefore a great effort has to be put on packaging development and reliability study. In this field and in order to monitor devices temperature, we have implemented a very simple integrated temperature sensor.

It is critical to reduce thermal interfaces between heating and sensing elements in order to reach a certain precision on temperature rising during transient operation. Conventional metallic sensor can be used (Pt, or Al made). Nevertheless, the drawback of such devices lies on their thermal response time. Moreover their thermal sensibility and the different thermal expansion coefficients in comparison to SiC leads to mechanical stress and finally makes the sensor subject to failure or derating.

In the case of a metal layer, the variation of the resistivity of the metal can be assumed to vary lineally with temperature, that makes simple the temperature rising estimation. However, we have preferred the integrated SiC sensor option to get a more reliable temperature monitoring.

Sensor and methodology description

Temperature sensor description

Unlike metallic sensor, an integrated sensor has the advantage to reduce interface layers effect and mechanical stress since the sensor is integrated with the device (figure 2). The sensor consists in a resistor that can be vertically or laterally integrated within the current limiting device. The lateral sensor has the advantage of being fully insulated by the Pburied layer.

In addition, the response time of the sensing element will be almost the same to that of the heating element. We have used a LVJFET current limiter with a high power density (140 kW/cm²), which has been previously checked [START_REF] Tournier | A 4H-SiC High power density VJFET as Controlled current limiter[END_REF], to demonstrate the possibility to measure the temperature evolution. The mobility variation with temperature is assumed to vary as describe in (1).

. Equation 1.

The dependence coefficient "a" was extracted through cryogenic measurements (a = 2.2 was found, close to literature reported results [START_REF] Ruff | SiC devices: Physics and Numerical Simulation[END_REF]).The temperature estimation technique, based on this law, will be analysed below and then illustrated by experimental measurements on a packaged current limiting device.

Temperature estimation methodology

The temperature estimation is based on the measurement of the resistance variation (caused by mobility variation) of the sensing device. Fig. 2 shows the integrated temperature sensor with the current limiter. The measured sensor resistance is inversely proportional to mobility. Measurements of the resistance variation, DR, with temperature Dq allow to simply compute the temperature, q, if one knows a reference resistance value at room temperature (R0, qamb):

. Equation 2.

(K includes sensor parameters: channel section and width). From (2), we can deduce the relation between q and Rq:

. Equation 3.

( ) 0 µ f a q q q ae ö ae ö ç ÷ = ç ÷ ç ÷ è ø è ø ( ) ( ) ( ) 0 0 1 1 L L R K R S S N q µ µ a r q q b b q q q é ù ae ö ae ö ´ê ú D = ´= ´´= ´= ḉ ÷ ç ÷ ç ÷ ´´ê ú è ø è ø ë û 1 0 0 0 R R R a q q q ae ö - = ´ç ÷ è ø

Experimental results

Method validation: IR measurements similar as presented in following figure .   IR-measurements were performed to check the temperature distribution between the heating device (H) and the sensing device (S). The figure 5 shows two cross sections along line L01 of figure 4 for two different VDS bias applied to the current limiter (VDS = 21V and VDS = 42 V). We can see from figure 5 that the contact metal temperature is similar in both the heating element (H, at 0.4 mm) and sensing element (S, at 1.2 mm) for both VDS applied bias. The IR measurements confirm the fact that the temperature is similar in both devices and corroborates the temperature measurement feasibility.

Temperature measurements

Once experimentally confirmed the homogenous distribution of temperature, electro-thermal measurements were carried out. Resistance measurements were done using two Keithley controlled sources. Different set-ups C1, C2 and C3 were used to get various dissipated power values. No heat sink was used for set-ups C1 and C2, while a cooling element was used in set-up C3. Another difference consists in the gate electrode polarization that was either leaved floating or connected to the source electrode to achieve different output current levels. For each set-up, a maximum case temperature Tmax=150°C was fixed. The measurement was stopped when Tmax was reached. IDS(VDS) curves where measured in steady state and the SiC temperature computed form the sensor. Current and temperature measurements are presented in figure 6. According to Mc Cluskey [START_REF] Mccluskey | Reliable use of commercial technology in high temperature environments[END_REF], plastics packaged devices are able to work at temperature around 125°C with the same reliability as ceramic or metallic package. Thermal ageing study shows that a stable behaviour is achieved at temperatures as high as 175°C. Thus, we have assumed a maximal case temperature of 150°C to perform our study. Back case temperature was measured with a K-type thermocouple. Measurements were done using a cellular current limiter as heater and a vertical VJFET as sensor. The die (4 mm × 4 mm) was soldered in a TO220 metallic package (figure 3). Thus, we have been able to measure the temperature in limiting operation, which is an important issue when dealing with current limiting devices since temperature is a feed-back parameter for power losses control. Furthermore, the use of an appropriate gate drive will permit to reduce power losses and temperature rising and, consequently to reduce the cooling system size.

Conclusion

We have proposed and illustrate a very simple way to perform on chip temperature monitoring of SiC current limiter with an integrated sensor. Electro-thermal measurements on the fabricated devices have shown that the current limiter is able to work at 205°C, under steady state conditions (320 V), without degrading their electrical characteristics. Different cooling set-ups have been used to illustrate both the temperature measurement and the high current density working operation of SiC devices. Next challenge will consist in the design of a fully integrated current limiter including temperature sensor to be able to prevent any thermal runaway of the packaged device. Such a device will be based on a combination of LVJFET, MESFET and a SiC thermal sensor.

For set-up C3, a maximum SiC temperature of 201°C is reached at a very high source power density (825 kW/cm²). These results show that the current limiter is able to work at 201°C under steady state conditions (320 V), without degrading neither the electrical performances nor the current limiting function. Results summarized in table 2 shows that high current density can be achieved as well as the temperature monitoring.

Fig. 2 .

 2 Fig. 2.-Cross-section of the integrated temperature sensor with the current limiter.

Fig. 4 .

 4 Fig. 4.-IR measurement at low drain bias (VDS = 42 V) Fig. 5.-Temperature cross-section along LI01 of the on heating (H) and sensing (S) chips left IR-picture (fig. 4).

Fig. 3 .

 3 Fig. 3.-Packaged current limiter (metallic TO-220 package). (die de 4mm ´ 4 mm) (limiter 460µm ´ 460 µm )

Table . 2. Measurements and estimated temperature for the various experiments (C1,C2, C3)

 .

	Experiment	C1 (VGS floating)	C2 (VGS = 0V)	C3 (VGS floating)
	COOLING	NO	NO	YES
	VDS max (V)	70 V	151 V	320 V
	IDS max (mA)	31.4 mA	17.3 mA	52.4 mA
	J MAX (A/cm²)	Total Area Source Area	38 A/cm² 1558 A/cm²	21 A/cm² 858 A/cm²	63 A/cm² 2580 A/cm²
	P max (W)	2.2 W	2.6 W	16.8 W
	TSiC MAX (°C)	204 °C	215 °C	201 °C
	TPackage MAX	156 °C	156 °C	151 °C