domain: spi phys

text:
S Duzellier

T Nuns

C Inguimbert

C Pons

R Rey

T Le Cocq

C Cavel²

Radiation response of 3J solar cells to radiation at low temperature

Degradation of 3G28C 3J Solar cells from AZUR SPACE has been evaluated for the specific JUICE environment. Electrons and protons data on full scale cells from two different batches is presented with a focus on the influence of temperature (annealing mechanisms) and particles species. Further analysis of the adequacy of NIEL dose approach is also discussed.

INTRODUCTION

Space radiations strongly affect the long-term performance of most active semi-conductor devices on-board spacecraft. Solar Generator Subsystem (SGS) is particularly exposed to solar and trapped particles as not shielded by spacecraft structure.

JUpiter ICy moons Explorer (JUICE, 2022) aims at investigating Jupiter and its icy moons. The Solar Cell Assembly (SCA) is first degraded during the cruise phase (≈ 7.5 years) under the influence of solar protons then by the harsh radiation environment around Jupiter (trapped electrons and protons). The expected cumulated mission dose corresponds to 1MeV electronequivalent fluence of 2.5-3 10 15 e-/cm² (End-Of-Life -EOLfluence).

Additionally, very specific operating conditions are encountered: at Jupiter the temperature will be around -120°C with very low (photons) solar flux (between 46W/m² and 54W/m²).

The JUICE SGS design consists in two identical deployable wings with the baseline Photovoltaic Assembly (PVA) based on 40mm x 80mm 3G28C solar cells from AZUR Space (thickness of 80μm, external bypass diode, 300μm ITO coated and grounded cover glass). JUICE SGS design relies on accurate estimate of radiation constraints and representative ground testing of cells from flight lot. Ideally low temperature and very low solar flux level shall be applied at ground testing.

II. TEST PLAN AND FACILITIES

In the course of technological development activities, test campaigns have been conducted on small test cells (20mm x 20mm) to investigate the defect introduction rate dependence on temperature and annealing mechanisms [START_REF] Park | Origin of the degradation of triple junction solar cells at low temperature[END_REF][2] [START_REF] Baur | Development status of triple-junction solar cells optimized for low intensity low temperature applications[END_REF]. The main objective was to answer the question: can final qualification of solar cells be carried out at room temperature?

The purpose of this second phase study was thus to check if same conclusions can be drawn from full scale solar cells testing and to acquire "reference" EOL data for the design of JUICE SGS. It was also crucial to evaluate degradation and annealing with great accuracy and reproducibility.

The test plan included irradiations with 1MeV electrons (e-) and protons (p+) separately or combined with in situ dark/light(3.7% AM0) I(V) at different temperature and fluences.

III. MAIN RESULTS

Typical response is provided in Fig. 4 showing I sc and V oc parameters at each step of the test sequence (electron testing). First observations here are that 1) no annealing is observed at low temperature (steps -125IR1 to 3, <0.3%), 2) but can be observed after warming up at RT (step -125IR4, I sc here 1.22%). This limited recovery can amount up to ≈3% for P max (5-6% for the largest electron fluence) which is worst case (WC).

Fig. 4. I sc and V oc parameters at each step of "nominal e-" testing.

In [START_REF] Park | Origin of the degradation of triple junction solar cells at low temperature[END_REF] P max degradation level much greater than I sc and V oc ones is observed and explained by the radiation-induced tunneling current taking place in the top and bottom sub-cells. We do observe here again this significant higher P max degradation as disclosed in Table 1. The plotting of P max versus electron or proton fluence shows that degradation follows a Logdependence law (as reported in [START_REF] Baur | NIEL Dose dependence for solar cells irradiated with electrons and protons[END_REF] and [START_REF] Salzberger | Degradation fitting of irradiated solar cells using variable threshold energy for atomic displacement[END_REF]). It is proposed here to fit EOL/BOL parameter ratios (or Remaining Factor RF) with function expressed in terms of NIEL dose or displacement damage dose in middle cell (GaAs) (DDD method [START_REF] Messenger | SCREAM: a New code for solar cell degradation prediction using the displacement damage dose approach[END_REF], Figure 2). Therefore, the comparison of combined data with separate electron and proton contributions at low temperature can be done.

With this graph, it is shown that electron, proton and combined degradations at low temperature fit quite well with the displacement damage dose approach considering the middle cell is the main contributor to the 3J degradation. This large set of data will be completed with data of annealing after one year, degradation with a second electron energy beam in order to investigate fitting curves and displacement dose approach. All data will be provided in the final paper. well with the displacement damage dose approach considering the middle cell is the main contributor to the 3J degradation.

Annealing effects have been investigated and RT annealing coefficients measured with accuracy.

Additional testing is planned for further investigate annealing and displacement dose dependence of degradation at low temperature.

Figure 1 -

 1 Figure 1 -test plan, generic sequence and picture of 2 cells mounted into MIRAGE (AXEL lab.,) of ONERA DPHY department in Toulouse

Figure 2 -

 2 Figure 2 -electron, proton and combined data fitted with NIEL dose fitting curve. Data lines up with good agreement.

Table 1 -MEAN DEGRADATION AT -125°C FOR ALL CONFIGURATIONS

 1

		Nnm. e-double e-nom. p+	double p+ combined
	I sc	-12,3%	-17,9%	-6,1%	-10,7%	-12,0%
	V oc	-6,5%	-7,7%	-7,8%	-8,9%	-8,8%
	P max -21,9%	-34,8%	-18,2%	-24,3%	-27,2%
	FF	-4,8%	-13,9%	-5,6%	-7,0%	-9,3%


Maxime Juston

Nicolas Damay

Christophe Forgez

Guy Friedrich

Karrick Stephane Vivier

Mergo Mbeya

Bogdan Vulturescu

Direct determination of a single battery internal resistances distribution using a heterogeneous model

Keywords: battery, heterogeneous electrical model, internal resistance distribution, non-invasive diagnosis

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Lithium-ion batteries are becoming one of the best solution to store energy in a wide range of applications, in particular among systems where weight or volume are major constraints, such as drones, cars and trains. In the industry, this technology is not often well known and it is to our advantage to develop a plain, yet robust model, that can be used in a wide range of operating conditions. In order to anticipate the voltage response of a cell to a current profile, a model can be used. Three main types of models exist, with a first one being the electrochemical model. It takes into account the chemical reactions, diffusion reactions and energy equations at the atom scale. Example of this model can be found in many articles such as in Di [START_REF] Di Domenico | Lithium-Ion Battery State of Charge and Critical Surface Charge Estimation Using an Electrochemical Model-Based Extended Kalman Filter[END_REF] where the author use an electrochemical-based model coupled to a Kalman filter to estimate a cell state of charge during usage, or in [START_REF] Yuan | A transfer function type of simplified electrochemical model with modified boundary conditions and Padé approximation for Li-ion battery: Part 1. lithium concentration estimation[END_REF] where the authors estimate the lithium concentration within the electrode. In [START_REF] Dong | Numerical modeling and analysis of the thermal behavior of NCM lithium-ion batteries subjected to very high C-rate discharge/charge operations[END_REF] the authors use an electrochemical model to simulate the thermal behaviour of a cell under very high current charge or discharge, that is more than 8C. Finally, [START_REF] Delacourt | Analyse des Performances et du Vieillissement des Batteries à Ions Lithium par des Méthodes de Modélisation Mathématique[END_REF] presents the construction of an electrochemical model. The interested reader will find in this long document all the information regarding the electrochemical modelling of a cell. However, the complexity of those models and the number of parameters needed require invasive measurements that are not possible for most battery managers.

Another main type of model is the black box model that uses mathematical description of a cell using probability and large data set to be constructed. Example of this model can be found in [START_REF] Saha | Prognostics Methods for Battery Health Monitoring Using a Bayesian Framework. Instrumentation and Measurement[END_REF] where the authors estimate the remaining life of a cell using a Bayesian framework or in [START_REF] Nuhic | Health diagnosis and remaining useful life prognostics of lithium-ion batteries using data-driven methods[END_REF] where the author use machine learning to estimate the cell remaining life as well. Those models can be built with various methodologies (machine learning, fuzzy logic, genetic algorithms) but have often in common the fact that they are constructed thanks to a large set of data. The result is "black box" model with no or very few physical meaning, which cannot be used to diagnose the internal state of the battery.

One last option in term of model is to use electrical equivalent circuits, also known as behavioural models. Those model are generally done at the cell level and model it thanks to usual electrical components, that is resistances, capacitances, voltage sources, etc. Even though they do not take into account all the atom-level reactions, a good accuracy can be obtained. The interested reader can find more information in the article of [START_REF] Berrueta | A comprehensive model for lithium-ion batteries: From the physical principles to an electrical model[END_REF] where the modelling of electrochemical phenomena using an electrical framework is explained. Articles [START_REF] Andre | Characterization of high-power lithium-ion batteries by electrochemical impedance spectroscopy[END_REF] and [START_REF] Kuhn | Modelling Ni-mH battery using Cauer and Foster structures[END_REF] focuses on the modelling of unusual electrical component, that are constant phase element and Warburg impedance in RC circuits. Finally, in the article [START_REF] Merla | An easy-to-parameterise physics-informed battery model and its application towards lithium-ion battery cell design[END_REF] introduces a multi-particle, diagnosis oriented model. Although those behavioural models are easier to parametrize, their parameters have to depend on current, temperature and State of Charge (SoC) to bring accurate results. Thus, those models usually require lookup tables that need many measurements to be completed.

All those models are useful to predict the performances of a cell in many operating conditions in order to optimize the sizing of a battery pack, choose an appropriate cooling system or predict its capacity to fulfil a mission. A previous work [START_REF] Kuhn | Modelling Ni-mH battery using Cauer and Foster structures[END_REF] led to a homogeneous equivalent circuit of a Ni-MH cell, built on a physical basis. This model only takes into account the main phenomena that occur inside a cell, and associate to them a combination of resistances and capacitances to model those phenomena.

Recently, the homogeneous model developed in our laboratory was extended to a heterogeneous one, the socalled "multibunch model", by [START_REF] Damay | Heterogeneous behavior modeling of a LiFePO 4 -graphite cell using an equivalent electrical circuit[END_REF]. It aims to model the heterogeneity of a cell thanks to a distribution of one or more parameters. The origin of those distributed parameters can have many sources that will be listed. This model is more accurate than a homogeneous one and although it may appear more complicated at first sight, it is faster to characterize than a usual electrical equivalent circuit model. Damay et al. found that the addition of the heterogeneous behaviour representation in the model allows the removal of the parameters SoC-dependencies (except for the open circuit voltage). This led to model that is faster to characterize: around 10 times faster than the previous homogeneous equivalent circuit model, which had a SoC-dependency precision of 10%. Once characterized, this model is able to simulate the cell over its whole range of operation. The characterization of the heterogeneity of a cell, through the distribution of one or more parameters, also brings valuable insights into its behaviour and its internal electrical properties. We expect this distribution to change during the cell lifetime due to ageing and the tracking of this distribution may be used as a diagnosis tool.

In this paper, we demonstrate that the measurement of a single discharge, used to build our model, allows us to simulate the cell in 4 different operating conditions and to extract the distribution of one of our parameters. For this article, this parameter is a single resistance modelling fast-dynamic phenomena such as charge transfer, current-collector resistance and electrolyte resistance. The distribution of this resistance for two different cells, one aged and one new, will be computed and compared in order to support the idea that this tool can be used as a non-invasive diagnosis tool. The problems related to the determination of this distribution will be presented and our solutions explained. In section 2 we present the construction of the model, from a homogeneous one to a heterogeneous one, its parameters and variables. The simulation and optimization algorithms are presented in section 3 while section 4 describes our experimental setup and the used cells. Discussions about the obtained results are in section 5 and last, the conclusion and perspectives of this study are presented in section 6.

Model and parameter dependency

Construction of the multibunch model

Our model is based on the assumption that one cell can be divided into elementary volumes that behave homogeneously and have specific electrical properties. Those volumes can be seen as paths that are used by electrons and lithium ion to travel through the cell internal parts. Those volumes are connected in parallel, as showed on figure 1.a. We assume that certain volumes, not necessarily spatially bounded, have close electrical properties when compared to each other. Because their behaviours are comparable, we group those volumes into a "bunch" that is considered homogeneous. By repeating this operation, the cell is [START_REF] Park | Probing electrical degradation of cathode materials for lithium-ion batteries with nanoscale resolution[END_REF], who found that a resistance distribution is found by scanning the surface of an active material sample. The number of bunches of our model, n, is a choice from the user and its value is discussed later in this section.

The bunches have a proper current I i but the same voltage: V cell . The sum of the bunch currents I i is the cell current: I cell . Each bunch is modelled using an equivalent circuit model (ECM). The proposed representation of a bunch comes from a previous article from Damay et al. [START_REF] Damay | Thermal modeling of large prismatic LiFePO4/graphite battery. Coupled thermal and heat generation models for characterization and simulation[END_REF] in which the hypothesises are the following:

• The double layer capacity is neglected because of its fast dynamic and the fact that our tests are conducted at constant current during thirty to sixty minutes. We call "high frequency" resistance, noted R HF , the sum of the so called ohmic resistance and the charge transfer one.

• In order to reduce the computation time, the Warburg impedance representing the diffusion phenomenon is modelled by a single RC. More RC circuits can be added to better model the diffusion, as presented in the studies of [START_REF] Kuhn | Modelling Ni-mH battery using Cauer and Foster structures[END_REF] where the transformation from a Warburg element into an infinite sum of RC is explained, and of Andre et al. [START_REF] Andre | Characterization of high-power lithium-ion batteries by electrochemical impedance spectroscopy[END_REF] where the modelling of a constant phase element with RC circuits is explained. The subscript d under those elements stands for diffusion. This hypothesis is assumed reasonable as the errors in constant current discharge will be concentrated during the first tens of seconds for simulation of several tens of minutes.

Those two hypotheses lead us to the multibunch dynamic model used in this paper, presented on figure 1c: one voltage source, one high frequency resistance and one parallel RC. Because each bunch represent elementary parts with different electrical properties, the value of the resistance will be different from one bunch to another. Among those parameters and to simplify our model, all the heterogeneity of the electrical properties is considered to be concentrated in the values of the R HF,i . We expect this parameter to be the most heterogeneous one in a cell due to varying distance to the current collector and variation in the manufacturing processes. The determination of the distribution of the R HF,i is done either by searching a mathematical distribution (a Weibull one), or by searching directly values fitting the experiments. Differences between both methods will be discussed later

Parameters values, dependencies and number of bunches

The voltage source associated to bunch number i has a values U(Q i ) that varies with the bunch charge like an increasing and nonlinear function, following the open circuit voltage (OCV) of the cell. The bunch charge is defined by its initial state and the bunch current by

Q i,t = Q i,t=0 + t 0 I i (t) dt (1)
Thus the value of the bunch voltage source U(Q i ) only depends on the proper SoC of the bunch. All others parameters have values that vary with the temperature and the current. The lookup tables we use to represent those variations are obtained through the analysis of the voltage response of the cell to a current pulse. By using different current values at different temperatures, one can extract the variation of the parameters through the analysis of the response, as can be seen on figure 2. The interested reader will found detailed protocol in the article from Damay et al [START_REF] Damay | Heterogeneous behavior modeling of a LiFePO 4 -graphite cell using an equivalent electrical circuit[END_REF]. Values of the parameters at different temperatures and currents are then computed as follow, where P represents a parameter and β P(I,T ) is a coefficient which is calculated thanks to measurements:

P(I, T ) = P(50%SoC, 1C, 25 • C) • β P (I, T ) (2)
where

β P (I, T ) = P meas (50%SoC, I, T ) P meas (50%SoC, 1C, 25 • C)
Each parameter has its own non-linearity look-up tables regarding the temperature and the current. If no measurement for the particular I, T does exist, a linear interpolation between the closest measured values is done. Conversely to the R HF,i values, the values of R d and C d are the same for each bunch. Thus no subscript i is added to their names. The distribution of the parameters R HF,i is done by optimization and will be discussed in section 3.2 A summary of the differences between a usual equivalent electrical model and the multibunch model can be found in Table 1 where I refer to the cell current, as opposed to the local current.

U(Q cell ) U(Q i ), i = [1..n] R HF R e (SOC, I, T ) + R HF,i (I, T ) R ct (SOC, I, T ) and distributed R d R d (SOC, I, T ) R d (SOC 50%, I, T ) • n C d C d (SOC, I, T ) C d (SOC 50%, I, T )/n
Regarding the number of bunches n, several configurations have been tested, and the optimal number found for a single cell seems to be the higher one. The higher the number of bunches, the more precise the model is. Because the optimization time needed to find a set of resistance is roughly a first order function of the number of bunch, the chosen number here is a compromise between accuracy and computation time. The impact of n on accuracy and computation time for a given method can be found on figure 3. The result regarding the value of n on accuracy and computation time for other methods are similar and are not presented here. All used methods are presented in section 3.2. For this study, we choose n = 20.

On figure 4 are displayed the voltages simulated for a various number of bunches for a given distribution of parameter. The error decrease with the increase of bunch number for multiple reasons. First of all, due to the fact that a SoC modification does not affect the parameter's value, the over-voltage is only a function of the current, which is constant here, and of the temperature, which is monotonically increasing. So in a first approximation, one can estimate that the over-voltage is a function of the temperature only. As a temperature increase lead to a decrease in resistance value, the over-voltage for one bunch is a decreasing function of the temperature, thus of the time. That is why the shape of the OCV, and especially the transition between the plateaus are visible on the 1-bunch voltage. This lead to a great error in the end of the discharge. Usual homogeneous model uses a SoC-dependency lookup-table to tackle this problem, with an increase in the resistance at low and high SoC (the effect is the same during a charge, where the over-voltage is usually quickly increasing near the end of the discharge, this is why the constant voltage protocol is so widely used). In this article, the SoC-dependency is not modelled by choice. As one can see on figure 5, increasing the number of bunches give a satisfactory result regarding the voltage drop in the end of the discharge, as the least resistive bunches begin to be completely discharged, letting more current flow through high resistances one, effectively increasing their over-voltage, and lowering the equivalent voltage of the whole model.

Another effect of the increase of the number of bunches is the disappearance of the transitions between the plateaus. This lead to a smoother curve, which is in accordance with the measurement one can get from a high-current battery discharge. This is due to the fact that the transition between the OCV plateau does not occur at the same time for all the bunches, as the least resistive bunches are discharged quicker. This phenomenon "hide" the plateaus transition by smoothing them over a larger SoC range. As a result, for an increasing number of bunches, the error decrease, until a certain point where an augmentation of the number of bunches do not significantly decrease the error, and where the computation time becomes too large for the achieved precision compared to other choices of number of bunches, as depicted on figure 3.

Simulation and optimization algorithms

Simulation of the cell

Based on a homogeneous bunch model where the parameter and non-linearity tables are known, we build the multibunch model by determining the IRD. That is, we need to find a set of resistances R HF,i that allows the model to simulate the battery for any discharge inside the operating conditions. This is done by using an optimization in order to find the best set of R HF,i . The flow chart of our coupled simulation-optimization algorithm is represented on figure 6.

The voltage is simulated thanks to several inputs: the measured temperature and current during a discharge, the initial state of the cell and the IRD. The simulation algorithm has three steps.

• First the electrical parameters values (resistances, capacities, time constants) are updated with regard to current and temperature (Equation ( 2));

• Then, based on the previous bunch currents and the assumption that they remain constant during a time step, the local states of charge Q i and voltages U(Q i ) and the voltage across the RC element are computed (Equation ( 1));

• Finally, the new bunch currents and cell voltage are computed by solving a linear system using the cell current, the locals OCV U(Q i ) and the voltages across the electrical parameters.

Optimization algorithm and methods

In order to simplify the definition of the optimization function, we chose to use high frequency conductances G HF,i instead of high frequency resistances R HF,i . This section is dedicated to the determination of the R HF values, that are distributed among the bunches. To find the best set, three methods are compared:

• WD method: an indirect determination using a Weibull Distribution;

• FD-WI method: a direct determination using a Free Distribution with a Weibull Initialization;

• FD method: a direct determination using a Free Distribution initialized by a uniform one between two values.

A Weibull distribution is a mathematical distribution law that allows a great variety of shape, expected value and initial value. It's probability density function is

f (x; λ , k, θ ) = k λ ( x -θ λ ) k-1 exp -( x-θ λ ) k (3)
Using the WD method was already done in a previous article [START_REF] Damay | Heterogeneous behavior modeling of a LiFePO 4 -graphite cell using an equivalent electrical circuit[END_REF]. Here we propose, on the one hand, an optimization of the three Weibull parameters of the indirect determination and, on the other hand, the possibility to directly find the R HF,i values that were previously extracted from the WD. As the aim of the IRD is to be monitored in order to make a diagnosis of the cell, a shape modification of the IRD must remain possible. Thus a Weibull distribution is not the best suited one as we could miss an outsider in the IRD during the life of the cell. Using a free distribution allows us to be able to find those points, at the cost of increasing the number of parameters, therefore the computation time, to the number of bunches n=20. In order to decrease the computation time, we imposed an ordered solution, meaning that an initial conductance G 1 is searched for and then the difference between G n and G n+1 (equation ( 4))

∀i, j ∈ [2, n], G HF,i = G HF,1 + i ∑ j=2 ∆G j (4) subject to ∀ j ∈ [2, n], 0 ≤ ∆G j ≤ ∆G max
This way of ordering our solution was necessary as the algorithm used in the software Matlab includes the optimization function lsqnonlin which does not have constraint options in order to sort the solution.

Only lower and upper bounds can be defined, which goes well with a research of the differences, as the lower bound can be defined to be strictly superior to 0, forcing an increase from one value to another. The objective function for a free distribution is then:

min G HF,i f e (t, I, T, G HF,i ) = [U sim (t, I, T, G HF,i )-U meas (t)]
(5) where f e is the error vector for a given discharge.

The objective function to be minimized is the difference between the simulated voltage and a measured voltage. Thus, we need a test used as optimization test in order to determine our IRD. This test, the setup and the cells used are presented in the next section.

Experimental conditions

We used two 40Ah LiFePO 4 /graphite cells from the same batch for this study. The discharge current limit is 2C, the charge current limit is 1C, for a voltage varying between 2,5 and 3,7V. The operating-temperature range is -15 to 50°C. The first cell used is already aged and its actual capacity is around 37.9 Ah. This cell was kept at room temperature for 2 years and various tests have been made to acquire data for different studies. No particular aging plan was used to age this cell. The second cell used for this article is a cell from the same batch, that was kept at 5°C, and which capacity is around 39.4Ah.

The first cell was used to make several discharges in various conditions in order to determine its IRD and to run several validation simulations in order to validate the behaviour of the model when used in different operating conditions. The same method has been used to determine the IRD of the second cell in order to compare it with the IRD of the first cell and discuss the applicability of the proposed IRD determination method for diagnosis purpose.

Experimental setup

As explained in section 2.2, the parameters values are varying with regards to the temperature. It is measured on the centre of the largest face of the battery by using a thermocouple. But as we explained previously, our model is not a spatial one, as no assumptions are made regarding the position of the elementary volumes that are grouped inside a given bunch. Thus the temperature to be taken into account for a bunch may be different from the measured surface temperature. Thus, as the measured temperature will be used to compute the parameters values of all bunches, the cell needs to be as homogeneous in temperature as possible. We chose to thermally insulate the cell with a 10-cm thick polyurethane box during all our tests, as pictured on figure 7, in order to have a quasi-adiabatic setup. We used glass wool on the top of the cell to be able to connect the power wires. The latter are also insulated with polyurethane foam. The whole setup was then placed into a climatic chamber to control the initial temperature of the cell and to connect it to a Bio-Logic system with a 100A booster.

We then consider the temperature to be the same in the surface and inside the battery. Support to this assumption can be found in the article of [START_REF] Damay | Thermal modeling of large prismatic LiFePO4/graphite battery. Coupled thermal and heat generation models for characterization and simulation[END_REF], in which the thermal exchange between the core of the battery and the ambient air, through the largest face of the cell, has been modelled. The experimental setup was the same except that the bottom face was cooled by a thermally regulated plate with thermal glue. The thermal resistance between the core and the largest face of the battery R in has been estimated to be of 0.8K/W. The thermal resistance between the largest face and the ambient air R out has been estimated to be of 141 K/W. As R in is very small compared to R out , we assume that the temperature difference between the largest face and the core of the cell is very small compared to the temperature difference between the ambient air temperature and the largest face temperature.

In the same article, a 3 dimensional thermal model was used to simulate the core temperature of the cell during a full discharge at 2C with an initial temperature of 18°C. The difference between the estimated core temperature and the surface temperature is less than 3°C when the surface temperature reaches 40°C and less than 1°C when the surface temperature reaches 25°C. As the cell is cooled by the bottom, effectively creating gradient of temperature inside the cell, thus this can be viewed as "worst case" scenario. This small difference in temperature is due to the fact that the internal heat resistance is very low, thus ensuring a good homogeneity of the temperature inside the cell. Considering a thermal difference of 3°C between the core and the surface of the cell leads to an error of the estimated electrical parameters used in this article of 8 to 20% for the high

Tests descriptions

The IRD determination test and the validation tests are constant current discharge from the higher cut-off voltage to the lower cut-off voltage. The IRD determination test consists in a discharge at low initial temperature (15 °C) and high current (2C, that is 80A). The validation tests are similar but with an initial temperature of 25°C and 38°C with a current of 1C and 2C for each temperature, as well as a test at 15 °C, with a current The maximum temperature allowed by the constructor has been hit during the 2C discharge with initial temperature of 38 °C, resulting in a partial discharge only. This discharge is not exploited here and no simulation have been made using experimental condition.

The table 2 contains the different operating conditions and the associated tests.

Results and discussions

Choice of the IRD determination test

Our model aims to model the heterogeneity of a single cell. Thus a discharge that brings out the most heterogeneous behaviour of the cell should be the best suited one in order to determine the IRD of our model. Once this IRD is determined, it will be used to simulate other operating conditions with a comparison to experimental measurements to determine if it can model the cell in a wide range of operating conditions. Based on our observations, the discharge with the most heterogeneous behaviour is a low temperature -high current one. We tried to determine IRD with other operating conditions, such as higher temperature or lower current, but it gave greater error when the IRD is applied to the other conditions of simulation. Thus this particular discharge at low temperature and high current was kept as determination test for the IRD.

The low temperature brings more dispersed impedance values as temperature is a major factor for impedance variation. The high current brings state of charge heterogeneity among parts of the cell, because our cell has a flat OCV as a function of the SoC. This lead to no electrochemical constraint between a low SoC and a mid to high SoC local part of the cell. Moreover, as the setup only is a quasi-adiabatic setup, the temperature increases in the case of a high current should be greater than with a lower current, effectively covering a larger temperature range. This will ensure our model robustness with respect to large temperature variations.

We used our IRD determination algorithm (figure 6), tests (see section 4.2) and methods (see section 3.2) to determine all three IRDs. The result of the indirect WD method, that was developed in a previous study [START_REF] Damay | Heterogeneous behavior modeling of a LiFePO 4 -graphite cell using an equivalent electrical circuit[END_REF], are compared on figure 8 with the results of the direct FD-WI method that we propose in this paper. From an initial temperature of 15 °C, we observe an increase up to 43 °C. The linear shape of the temperature curve supports the quasi-adiabatic setup that we did, as a thermal cooling from the cell would be detected as a variation of the derivative of the temperature curve.

As the impedance is strongly linked to temperature in a cell, we observe a particular shape of the voltage curve. In a first part a voltage drop due to ohmic resis-tance can be seen. Then the voltage increases as the state of charge decrease. As the cell is an LFP/graphite cell, the OCV is nearly-flat at the beginning of the discharge, which should induce a flat voltage during discharge. However, the rising temperature of the cell leads to a decrease of its impedance and, thus, leads to an increase of the cell voltage. This happens until the OCV of the cell decreases at low SoC, effectively decreasing the cell voltage.

Relevance of the determined IRD to a diagnosis purpose

The aim of the optimization is to find a suitable set of resistances to model the cell. However, we also want to obtain a distribution that is physically meaningful. Because the goal of the IRD is to monitor the internal state of the cell, we have to exclude the results from the WD method from the acceptable results, as it is so far not proved that the resistances inside a cell vary with respect to a Weibull distribution. Besides, a Weibull distribution is constrained by equation (3) and cannot represent all possible IRD that a cell may have.

We found that the FD method may lead to a questionable solution. On figure 9 are presented the distributions that our algorithm finds for each method presented in section 3.2. The FD method brings a distribution with a few equals & subsequent values and one greater value at the end. Such an evolution is possible but is considered less likely because it doesn't suit the hypothesis of a continuum of electrical properties inside the cell. The FD-WI method brings a better regularity among the IRD, as expected from the model construction. As we find it hard to have a consistent criterion to evaluate the qualities of those solutions, we recommend the FD-WI method for the determination of the IRD as it gives satisfactory results regarding the hypothesis of an internal electrical-properties continuum. It has the advantage to be quickly initialized by the WD method in which only three parameters have to be optimized, and to avoid purely mathematical IRD.

Validation tests

To verify that the determined IRD is consistent, we simulated discharges for other operating conditions. The results are presented on figure 10. On each subfigure are displayed the measured voltage and the simulated voltage for the indicated operating condition. The results validate our IRD determination method as the error is kept below 25mV for an average voltage around 3.1V (that is to say, less than 1% relative error). The various operating points that are simulated strongly support the fact that our model has a physical meaning, that can be used to track the IRD during the cell lifetime.

The different simulations show a smaller error for a higher discharge current. We suppose that high currents are more likely to induce a heterogeneous behaviour inside the cell, which are to be modelled by this heterogeneous model. Although the determination test was also a high current test, which can lead to a bias in the error, IRD extracted from low-current test did not give better result, neither for high current, nor for low ones, as specified in section 5.1. Thus this IRD is considered to be the best one we could find. The overall shape of the simulation is satisfying compared to measured voltage, which validate our goal to simulate the cell voltage over a wide range of operating points.

As the cell voltage can be simulated during a constant current discharge, with only one IRD obtained by optimization on a single discharge at low temperature and high current, it becomes possible to track the evolution of this IRD. This open the way to a non-invasive diagnosis tool for a single cell or a group of parallel cells.

Direct IRD determination as a diagnosis tool

After having determined that the proposed IRD determination method is suitable for a diagnosis purpose, we tested this method in its usage as a diagnosis tool to compare the IRD of two cells that have different ageing conditions (those cells were described in section 4). We expect the distribution of the aged cell to show higher values of resistances, as its resistance should have increased with use. The IRD values are displayed on the Y axis and the number of the bunch on the X axis. Both cells show very similar shape of IRD, with the aged battery having a higher equivalent resistance value as expected. The values of the aged cell resistances are a bit higher by around 4mΩ for bunches 4 to 14 but the shape is preserved, as if the aging impacted bunches 1 to 18 homogeneously. The bunches 18 to 20 are interesting as we see a crossover between the IRD. This can be due to several causes, including • a difference between the two batteries initial IRD;

• an unexpected degradation process;

• a difficulty in finding a global optimum as local minima may exist (equation ( 5));

• an artefact due to the model construction hypothesis (very simple bunch model, no separation of the 2 electrodes behaviours...).

So far, no evidence has been found toward any of those causes, and further studies should help us to have a better understanding of these results. Tracking a single cell IRD during its lifetime would help to distinguish that, as its resistance should only grow.

The IRD analysis shows good consistency with the expected resistance values, as the aged cell has a higher average value of equivalent resistance. More tests on other types of battery are necessary to conclude about the usage of this criterion as a SoH diagnosis tool, but those first results are encouraging. Moreover, new tool can be used, based on the IRD shape or values, such as a heterogeneity level that is the difference between the min value of the IRD and the max value of the IRD. As the cell ages normally, we expect this ratio to decrease as less resistive parts are privileged paths for the current, inducing more degradations and a resistance increase. As this tool can be used on aged cell without a need to rebuild non-linearity tables, it makes a good candidate for a SoH diagnosis tool. Other criteria such as the standard deviation of the IRD could be used and will be investigated in other studies. This work may bring new methods in order to diagnose the internal behaviour of the cell, without opening them.

Conclusions and perspectives

In this paper we presented an equivalent circuit model called multibunch model, whose goal is to model the heterogeneity of a commercial LiFePO 4 /graphite cell. This model was initially developed in a previous study [START_REF] Damay | Heterogeneous behavior modeling of a LiFePO 4 -graphite cell using an equivalent electrical circuit[END_REF] and, in this article, we demonstrate that it can be used as a non-invasive diagnosis tool. This model gives access to an internal resistance distribution (IRD) which is assumed to represent the distribution of the cell internal electrical properties.

We propose a direct determination method of the IRD in this paper, as opposed to the indirect determination method of the previous study. This allows us to have a responsive distribution that is not bound to a mathematical law of a distribution. As the cell internal properties may vary in many different ways, the ability to spot outsiders values within the IRD is a strong improvement of the method.

An experimental protocol has been proposed for the IRD determination and three determination methods have been compared. The first one is the previous study indirect determination method using a Weibull distribution. The two others are a direct determination using a free optimization of the IRD values, after initializing them with either a Weibull distribution or a uniform distribution. We demonstrated that the best method is to first determine the best Weibull distribution, then to allow an optimization of the internal resistances around the Weibull values in order to fit the measurements obtained thanks to the IRD determination protocol. The found IRD was validated for different other operating conditions with an error of less than 25mV for an average voltage of 3.1V.

The comparison of the IRD of two different cells highlights the new insights that our method provides about a single cell internal state. We propose an "heterogeneity level" as a new criterion to determine a cell state of health, which allows to have a better accuracy on the cell internal resistance. Other criteria could be found to better represent the cell ageing with regard to the IRD information, but further study and tests on other batteries are necessary to have a motivated opinion on the subject. This is the next step of our work and will be presented in future publications.

The use of thermal model to estimate the core tem-perature, and the correct the estimation of the electrical parameter may be an interesting perspective of this study. Especially when low errors are achieved and where small variations of the parameters due to temperature may account for an unneglectable part of the remaining error.

Also, application of this method to a small battery pack which contains one aged cell and several new cells will be done in order to validate this method as a diagnosis tool for larger battery pack. The construction of the multibunch model allows it to model several batteries in parallel. As our heterogeneous model is made of parallel bunches, to model parallel cells is easy to parametrize, as we only need to increase the number of bunches.

Finally, we plan to apply this simulation method to other cell chemistries to verify that we are able to simulate them over their whole operating range, to get an insight of their internal states, and to follow the evolution of their IRD in order to have a non-invasive diagnosis tool for large battery packs.

Figure 2 :

 2 Figure 2: Expected voltage response of a cell

Figure 3 :

 3 Figure 3: Influence of the number of bunches

Figure 4 :

 4 Figure 4: Influence of the number of bunches on the simulated voltage

Figure 5 :

 5 Figure 5: Zoom on the end of the discharge

Figure 6 :

 6 Figure 6: flow chart of our algorithm

Figure 7 :

 7 Figure 7: Experimental setup

  (a) Overall measurements and simulations (b) Zoom on the end of the discharge

Figure 8 :

 8 Figure 8: Measured temperature, voltage and simulated voltages

Figure 9 :

 9 Figure 9: Comparison of the IRD

  Figure 10: Measured temperature, voltage and simulated voltage for different operating conditions

Figure 11 :

 11 Figure 11: Comparison of the IRD for an aged and a new cell

Table 1 :

 1 Differences between a homogeneous and the multibunch model

	Parameter	Homogeneous	Multibunch, order n
	OCV

Table 2 :

 2 Experimental conditions of the tests

	Current →	1C (40A)	2C (80A)
	Initial temp. ↓
	15 °C	Validation	IRD determination
		test	test
	25 °C	Validation test	Validation
		test	test
	38 °C	Validation test Aborted validation
		test	test
	of 1C.

Acknowledgements

This work was funded by the French railway company SNCF as part of a CIFRE Thesis between SNCF, Université de technologie de Compiègne and the ANRT (National Agency for Research and Technology).


Mengfei Fan

Rui Kang

Zhiguo Zeng

Modeling dependent competing failure processes based on stochastic hybrid systems

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Industrial components, systems and products might fail due to multiple failure processes, i.e. wear, corrosion, fracture, fatigue and etc. [START_REF] Jiang | Modeling zoned shock effects on stochastic degradation in dependent failure processes[END_REF]. In literature, failure processes are categorized into two types: degradation processes (or soft failures) and catastrophic failure processes (or hard failures) [START_REF] Li | Reliability modeling of multi-state degraded systems with multi-competing failures and random shocks[END_REF]. In general, soft failures are modeled by a continuous degradation process, e.g., diffusion process [START_REF] Lemoine | On failure modeling[END_REF], continuous-time semi Markov's process [START_REF] Lin | Integrating random shocks into multi-state physics models of degradation processes for component reliability assessment[END_REF], etc. Hard failures are modeled by a random shock process, e.g., a Poisson process [START_REF] Lemoine | On failure modeling[END_REF][START_REF] Lin | Integrating random shocks into multi-state physics models of degradation processes for component reliability assessment[END_REF][START_REF] Klutke | The availability of inspected systems subject to shocks and graceful degradation[END_REF]. Soft failures occur when the total degradation exceeds its threshold. Hard failures are caused by the random shock process in various shock patterns such as extreme shock, cumulative shocks,  -shocks, m -shocks, run shocks [START_REF] Rafiee | Reliability modeling for dependent competing failure processes with changing degradation rate[END_REF].

In practice, those processes are often dependent, making the reliability modeling more complex and therefore referred to as multiple Dependent Competing Failure Processes (DCFP). As a result of the dependence, the arrival shocks may change some attributes of the degradation process [START_REF] Rafiee | Reliability modeling for dependent competing failure processes with changing degradation rate[END_REF] or bring additional damages to the degradation process [START_REF] Peng | Reliability and maintenance modeling for systems subject to multiple dependent competing failure processes[END_REF][START_REF] Jiang | Reliability and maintenance modeling for dependent competing failure processes with shifting failure thresholds[END_REF]. The hard failure process may be effected by the degradation process in terms of shock intensity [START_REF] Caballe | A condition-based maintenance of a dependent degradation-threshold-shock model in a system with multiple degradation processes[END_REF][START_REF] Huynh | Modeling age-based maintenance strategies with minimal repairs for systems subject to competing failure modes due to degradation and shocks[END_REF][START_REF] Huynh | A periodic inspection and replacement policy for systems subject to competing failure modes due to degradation and traumatic events[END_REF][START_REF] Bagdonavicius | Statistical analysis of linear degradation and failure time data with multiple failure modes[END_REF], shock loads [START_REF] Fan | Multicomponent lifetime distributions in the presence of ageing[END_REF], hard failure threshold [START_REF] Jiang | Reliability and maintenance modeling for dependent competing failure processes with shifting failure thresholds[END_REF] and etc. Various models are developed to calculate the reliability of the DCFPs. For example, [START_REF] Peng | Reliability and maintenance modeling for systems subject to multiple dependent competing failure processes[END_REF] assume that the arrival of shocks would bring an abrupt increase to the normal degradation process and develop a model to calculate the reliability. [START_REF] Jiang | Reliability and maintenance modeling for dependent competing failure processes with shifting failure thresholds[END_REF] develop a reliability model to consider the dependency between degradation and shocks where the failure threshold is shifted by random shocks. [START_REF] Rafiee | Reliability modeling for dependent competing failure processes with changing degradation rate[END_REF] consider the circumstance where the degradation rate is shifted by different shock patterns. [START_REF] Jiang | Modeling zoned shock effects on stochastic degradation in dependent failure processes[END_REF] categorize shocks into different shock zones based on their magnitudes, and consider dependence where shocks in different zones have different effects on the degradation process.

In the existing reliability models concerning DCFP, systems fail when either a soft failure or a hard failure happens. Hence, in those models, system's reliability is regarded as the probability that neither soft failure nor hard failure occurs in a period of time, which is calculated by the product of the conditional probability of no soft failure and that of no hard failure conditioned on parameters shared by the two probability functions. Due to the introduction of dependence, the developed reliability functions are usually too complicated to be evaluated analytically. Rather, the reliability has to be calculated by Monte Carlo method, which needs large amount of computation to ensure the accuracy of the results.

Modeling dependent competing failure processes based on stochastic hybrid systems

Mengfei Fan & Rui Kang

Center for Resilience and Safety of Critical Infrastructures, Beihang University, Beijing, China

Zhiguo Zeng & Enrico Zio

Chair System Science and the Energy Challenge, Fondation Electricite de France (EDF), CentraleSuplec, Universit Paris-Saclay, Grande Voie des Vignes, 92290 Chatenay-Malabry, France Energy Department, Politecnico di Milano, Milano, Italy ABSTRACT: A novel approach is developed to calculate the reliability of dependent competing failure processes (DCFP) based on stochastic hybrid systems (SHS). The DCFP considers the dependency between a degradation process and a random shock process, where the arrival of a random shock will bring an abrupt change to the degradation process. A SHS model is developed to describe the evolution of system behaviors, where the degradation process is characterized by stochastic differential equations (SDEs). The SDEs are dependent on the transition of system modes, which represent the random shock process. The conditional moments of the degradation measure in each system mode is derived. Based on the conditional moments, the system reliability and its lower bound are estimated using First Order Second Moment (FOSM) method and Markov's inequality, respectively. Results show that the SHS-based method yields an accurate estimation of reliability with less computational costs compared to traditional Monte Carlo-based method.

Stochastic Hybrid System (SHS) models are capable of characterizing the system behaviors that involve both continuous and discrete stochastic processes [START_REF] Hespanha | Stochastic Hybrid Systems: Application to Communication Networks[END_REF]. Using Dynkin's formula, the conditional moments of each system state of an SHS can be determined by solving a set of differential equations [START_REF] Hespanha | Stochastic Hybrid Systems: Application to Communication Networks[END_REF]. In this paper, we try to exploit the SHS model to reduce the computational costs of calculating the reliability of DCFPs. As an initial attempt, we consider the system discussed by [START_REF] Peng | Reliability and maintenance modeling for systems subject to multiple dependent competing failure processes[END_REF] in this paper, where the degradations depend on random shocks, i.e., the arrival of a random shock would cause an abrupt increase on the cumulative degradation measures.

The reminder of this paper is organized as follows. Sect. 2 gives a brief introduction of stochastic hybrid systems and its solutions. Sect. 3 formulates the DCFP as an SHS and derives the conditional moments of the system states based on Dynkin's formula. The conditional moments are used in Sect. 4 to estimate the reliability and its lower bound for the DCFP. Sect. 5 presents a numerical example to demonstrate the proposed model and compares its solution with the results from Monte Carlo simulation. A brief conclusion of the this paper is made in Sect. 6.

BACKGROUND KNOWLEDGE ON SHS

Generally speaking, stochastic hybrid systems are suitable to characterize systems whose state evolution contains both continuous stochastic processes and discrete stochastic processes. A latest formally defined SHS model proposed by [START_REF] Hespanha | Stochastic Hybrid Systems: Application to Communication Networks[END_REF] is applied in this paper. This SHS model is a subset of a more general class of stochastic processes referred to as the piecewise-deterministic Markov's processes [START_REF] Davis | Markov models and optimization[END_REF]. For formal definitions and analytic solution methods for the SHS, one can refer to [START_REF] Hespanha | Stochastic Hybrid Systems: Application to Communication Networks[END_REF][START_REF] Hespanha | A model for stochastic hybrid systems with application to communication networks[END_REF][START_REF] Hespanha | Modelling and analysis of stochastic hybrid systems[END_REF].

The state space of a SHS is a combination of a discrete state

  q t Q  and a continuous state   n t  x
, where Q is a finite set containing all the possible discrete modes of the system. Generally, a SHS is defined by the following items [START_REF] Hespanha | Modelling and analysis of stochastic hybrid systems[END_REF]:

(1) A SDE that describes the evolution of the continuous state:

              , , , , , t d t f q t t t dt g q t t t d  x x x w (1)
where

: k t   w is a k-dimensional
        , , , j q t t t t o t     x (2)
and if transition j occurs at time T the combined state of the SHS will be reset instantaneously by:

            , , , . j q T T q T T T    xx (3)
According to [START_REF] Hespanha | Modelling and analysis of stochastic hybrid systems[END_REF], given a test function

      , , , q t t t  x :, n Q     
which is twice continuously differentiable with respect to x and once continuously differentiable with respect to t , we have:

              ,, ,, dE q t t t E L q t t t dt        x x (4)
where    ,, L q x t  is known as the extended generator of SHS, and for   ,,

n q x t Q                               2 2 1 , , , , , , : , , ,, 1 , , , , 2 , , , , , , , l jj j
q x t q x t L q x t f q x t xt q x t trace g q x t g q x t x q x t t

q x t q x t                      (5)
q x t 
with respect to x , respectively. Generally, test functions are defined as functions with respect to measures of the continuous state. Hereby, we define a family of test functions of the form:

    ,, 0 i m m i q i x q q qx qq          (6)
where m can be any natural number of interest. Let

    i m q t 
denote the m -order conditional moment of the continuous state x , which is defined by

                :, Pr . ii mm qq m ii t E q x E x t q t q q t q           (7)
Then, by (4)(5) we can obtain a group of differential equations with respect to arbitrary order condi-

tional moments    ,, i m qi t q Q m    :           ,. i i m m q q d t E L q x dt     (8)
By solving the differential equations in (8) , we can obtain knowledge about dynamics of the continuous state, based on which further analysis can be conducted.

3 SHS-BASED FAILURE BEHAVIOR MODELING

System description

As shown in Figure 1, we consider a system suffering two dependent failure processes: soft failures caused jointly by continuous degradation and additional abrupt degradation increment due to random shocks, and hard failures caused by sudden overload from the same shock process. Dependence exists among the two process, that is, the arrival of a shock would bring an additional damage to the degradation process. Failures occur whenever one of the following two events happens: i. the degradation process reaches its threshold, denoted by H ; ii. a shock whose magnitude exceeds a critical strength level, denoted by D , occurs.

Additional assumptions of this model include: (1) The continuous degradation process is modeled by a Stochastic Differential Equation (SDE).

 

, t dx t dt b dw     (9)
where  

xt is the degradation measure,  is a constant denoting the degradation rate, t w represents the Wiener process with

SHS modeling of the system

The SHS concerning the DCFP is described by the state-transition diagram in Figure 2. The system is a combination of (i) a continuous-time, discrete-state stochastic process     12 , q t q q  , where 1 q and 2 q are two possible states of the system. When 1 qq  , the system is degrading according to the degradation process, while 2 qq  indicates that the system fails due to hard failures; and (ii) a continuous-time continuous-state stochastic process   xt that describes the degradation of the system over time. As shown in Figure 2, the system's discrete state starts from 1 qq  , under which the dynamics of continuous state   xt is governed by ( 9). When a shock with magnitude less than D arrives, system's discrete state remains at 1 qq  , and the continuous state  

xt will be reset:

    , x t x t d
 where d is the damage caused by the random shock. When a shock with magnitude over D arrives, the system transfers to the hard failure mode 2 qq  and the continuous state   xt is reset to zero. The transition rates and transition reset maps of the above SHS are defined as follows:

    1 1 2 1 2 2 , ,: 0, 1, ,: 0. W W W W D qq qx qq D qx qq                                        (10)         11 22 , :
, , , : ,0 , 11) In this paper, we are interested in the evolution of degradation level x . So we define the test functions

q x q x d q x q     (
    , i m q qx  in each discrete state   , 1, 2 i qi 
to be:

        1 2 1 1 2 2 ,, 0
,. 0 m m q m m q x q q qx qq x q q qx qq                  (12)
where m can be any natural number of interest. We also define the conditional moments of the continuous state x in each discrete state by: m m q m m q t E x t q t q q t q t E x t q t q q t q

                       
              (13)
With the conditional moments defined above, we can obtain any order of conditional moments of interest by choosing m . It is obvious that     0 i q t  represents the probability that system's discrete state belongs to i qq  . Substituting (9)(10)(11)( 12) into (5), we have the extended generator of the SHS as:

                                          11 1 11 1 21 2 2 2 10 1 12 0 0 2 ,, 1 , 2 ,, ,, , , . mm
m qq q m qq m q qq q x q x L q x b xx q x d q x qx L q x q x                    (14)
According to ( 8) and ( 14) , we have differential equations with respect to the conditional moments of 𝑥 as follows:

                                1 1 1 11 21 12 2 1 1 2 0 00 2 1 1 2 , . m m m q q q m m k m k qq k qq d t m t b m m t dt m t d t k d tt dt                              (15)
Because of the fact that when 2 qq 

, the continuous state   xt always equals to zero, in (15) we only concern the zero-order conditional moment at state 2 q , i.e. the probability that the system fails due to a hard failure. It's clear that

    12 00 1, qq   (16)
since the system must belong to one of the two discrete modes. By (15), we can obtain the following set of differential equations:

              11
                                                 (17)
By solving ( 16) and ( 17 x q q t   and standard variance

  1 x q q t   . The   1 x q q t   and   1 x q q t   can be calculated by                                       1 1 1 1 11 11 1 1 0 1 2 2 11
2 21 00 ˆ, ˆ. q x q q q x q q qq qq t t E x q q t t D x q q E x q q E x q q tt tt

                     (18)
Therefore, the reliability of the DCFP is

                                  1 1 1 1 1 1 11 1 0 Pr , Pr Pr ˆ.
ˆx q q

x q q x q q q x q q R t q t q x t H q t q x t H q t q Ht P q t q t Ht t t

                                  (19)
The lower bound of reliability can also be estimated using Markov's inequality. According to Markov's inequality (Pishro-Nik, 2014), if X is an nonnegative random variable and Pr .

E x t x t H H 

This inequality could be easily extended to the following inequality with respect to the conditional probability case:

        1 1
Pr .

E x t q q x t H q q H     (22)
Therefore, the lower bound of reliability is estimated by . q q R t q t q x t H q t q x t H q t q q t q x t H q t q E x t q t q q t q H In this paper a stochastic hybrid system characterizing systems subject to dependent competing failure processes are developed. The dependence is modeled by assuming the arrival of a random shock would lead to an abrupt increment of the degradation process. By solving the SHS-based DCFP model, dynamics of the conditional moments of degradation are obtained. Based on that, the lower boundary of system's reliability is estimated by Markov's inequality, and the exact reliability is approximated by moment method. Results show that the developed method provides accurate predictions of system's evolution with better computational performance.

                                              1
t t H                            (

1

  , and b is a constant.(2) Random shocks arrive according to a Poisson process with intensity  .(3) The magnitude of each shock load, denoted by i W , is an i.i.d. random variable following normal dis-To consider the dependence, we assume the arrival of each shock would bring a degradation increment of a constant d .

Figure 1 .

 1 Figure 1. (a) Soft failure process (b) hard failure process

  ), we are able to find out (i) the probability of the SHS residing in each discrete mode, and (ii) the first-order and second-order moments of   xt conditioned on   qt .

Figure 2 .

 2 Figure 2. State-transition diagram for the SHS

  is carried out using the parameters listed in Table1. To demonstrate the developed model, analytic solutions obtained by solving the differential equations in (17) are compared with results of Monte Carlo simulation, as shown in Figures3-6. Results show that the dynamics of are very accurately predicted by the SHS model. The estimated reliability by moment method is consistent with results by Monte Carlo simulation. The estimated lower bound provide a relatively conservative reliability estimation. Besides, the running time of Monte Carlo simulation is 2796.3 times more than that of the developed analytic approach (By Monte Carlo simulation: 548.08s, by analytic approach: 0.196s).

  Figure 3. Results of

Table 1 .

 1 Parameter values ______________________________________________ Parameter Value ______________________________________________

ACKNOWLEDGEMENT

This work has been performed within the initiative of the Center for Resilience and Safety of Critical Infrastructures (CRESCI, http://cresci.buaa.edu.cn). This research is supported by the National Natural Science Foundation of China under grand number 61573043.