A Reliability Modeling Method Based on Phase-Type Distribution Aiming at Shock Model

The traditional shock model generally describes the magnitude of the cumulative damage caused by a random shock sequence and compares the magnitude with a predetermined threshold to obtain the failure time of a component. There are two limitations in this kind of models in practice: First, the statistical characteristics of the damage due to a single shock may be difficult to obtain, which means the magnitude of the damage may not be described by an appropriate distribution; Second, the cumulative shock magnitude may be difficult to measure, or it may be difficult for a failure mode to be described by a threshold, meaning that the magnitude of the damage and the threshold may not be compared with each other. Considering both failure and censored samples, a reliability modeling method is proposed in this work to address the above problems. The shock model is first established by using both continuous and discrete phase-type (PH) distributions. Then the parameter estimation method of the shock model is derived based on EM method and the identifiability of the parameters in PH distributions is also given. Finally, the adaptability of the model is analyzed using three different types of simulation cases.


I. INTRODUCTION
In reliability theory, the shock model is often used to describe the probability distribution of a component's reliability under a random shock environment. The model was first proposed by Esary and Marshall [1], as shown in FIGURE 1.
When the cumulative magnitude (blue line) of the damage caused by the shock sequence (green line) exceeds the certain threshold (red line), the component fails. This moment, the reliability of the shock model satisfies: where N (t) represents the number of shocks that the component receives at moment t, and S represents the failure threshold. In the early literature [2], [3], in order to obtain the The associate editor coordinating the review of this manuscript and approving it for publication was Yu Liu . analytical expression of the parameters, the exponential distribution was frequently used to describe the magnitude of the damage due to a single shock, and the arrival time sequence was also characterized by a Poisson process. However, in many cases, the exponential distribution may be inconsistent with the actual situation. Some follow-up researchers VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ developed the modeling method with more types of distributions [4], [5], which greatly extended its adaptability. In addition to the use of various distributions, different types of shock magnitudes have been simultaneously considered and introduced into the model to characterize the multiple effects of various kinds of shocks on the life of a component. This kind of multi-shock model can be generally divided into two categories. One supposes that the magnitudes of damages has various forms, but only one failure threshold. For example, Eryilmaz and Tekin et al. [6] classified shocks into small type and extreme type and then compared the cumulative magnitude with a threshold to obtain the distribution of the failure time. The other supposes that different failure thresholds correspond to different failure modes, considering the competitive risk of cumulative magnitudes caused by multiple shocks in different failure modes [7], [8].
As the complexity of the modeling method increases, analytical methods usually fail to meet their demands. In many cases, only the approximate numerical result of the parameters can be obtained.
Phase-type distribution (also known as PH distribution) is a kind of probability distribution cluster based on continuous-time or discrete-time Markov process [9], which is widely used in shock model constructions. Neuts and Bhattacharjee [10] first used PH distribution for shock modeling. Literature [11], [12] used PH distribution to characterize the inter-arrival times between adjacent shocks and the magnitude of the damage due to a single shock, and then analyzed the reliability of the system.
The following works further extended the application scenario of PH distribution in shock modeling process. Eryilmaz and Kan [13] proposed an extreme shock model which has a change in the distribution of the shock magnitudes. In this article, the inter-arrival times follow PH distribution. Gong et al. [14] considered a system subject to random shocks from various sources with different probabilities, and used PH distribution to model the variables for discussing the reliability characteristics of the system. In [15], two random threshold shock models for a repairable deteriorating system were established using PH distribution and an ordering policy and a replacement policy ware then investigated to minimizing the long-run average cost rate. Zhao et al. [16] proposed a mixed shock model with a change point which is used to divide the failure process into two different stages to fit the situation of the changing of the failure mechanism during the system operation. They then proposed a multi-state shock model [17], in which the system failure patterns are mutative for different states. In this article, shocks were classified into three different types by the magnitudes of their damages. A type III shock causes a complete failure whenever it arrives, A type II shock triggers the system down to a lower state, and A type I shock has different effects in different states. They used continuous PH distributions to describe above assumptions. In [18], two stress-strength models with both single and multi-component were proposed and the systems with strength following the mixture of discrete PH distributions was considered. Montoro et al. [19] and [20] and Montoro and Rafael [21] considered two failure modes including shocks and continuous wear where the PH distribution is used to characterize the magnitude of the damage caused by a single shock.
Due to the different physical and geometric properties, the actual damage of a material caused by a single shock may be related to its cumulative magnitude of the damage [22]. In such cases, the statistical properties of a shock sequence are substantially unchanged, while the actual magnitude of the damage due to a shock in sequence may vary with increasing cumulative damage. For example, with a decline in performance, the actual shock magnitude tends to expand for electronic devices, which therefore does not satisfy the i.i.d. hypothesis. In response to this problem, Wang et al. [23] proposed a modeling method based on Weibull distribution to characterize the shock magnitude. In addition, models with different types of distributions related to different discrete degradation stages were introduced to characterize the evolution of cumulative magnitude caused by shocks [17], [24].
The above modeling methods have a common feature, that is, the cumulative magnitude of the damage due to shocks is characterized by a probability distribution and is required to be compared with a threshold to determine whether the component has failed. There are two limitations in this kind of modeling method in practical applications: First, the shock magnitude may not be characterized by a suitable distribution in some cases, or it may be difficult to grasp the changes law of the statistical characteristics of the damage. Second, the cumulative magnitude may not be compared with a failure threshold when it is difficult to measure the cumulative magnitude by direct or indirect method, or according to the failure mechanism, some failure modes are difficult to judge with one or more simple thresholds, such as fatigue. It is often difficult for traditional modeling methods to function well under such limitations.
Aiming at the issues mentioned above, this article proposes a modeling method based on continuous and discrete PH distributions, considering both failure and censored samples with the data of failure time, censored time and number of shocks. Then the EM method is used to derive the expressions of the parameter estimations, and the identifiability of the model parameters is also discussed. Finally, the adaptability of the proposed modeling method is tested and analyzed by three different types of simulation cases.
The paper is organized as follows: Section II briefly introduces the definitions and properties of both continuous and discrete PH distributions. Section III introduces the modeling process of the shock model based on the distributions mentioned above, and derives the reliability expression under shock impact. Section IV uses the EM method to derive the parameter estimation expressions based on the observed data of both failure and censored samples, and discusses the identifiability of the model parameters in such modeling method when using PH distributions. Section V introduces three different types of simulation cases, in which the corresponding sample data is generated, the parameter estimation results are obtained, the statistical characteristics of the model are compared with the assumptions in each case, and the adaptation of the method is also analyzed. Section VI summarizes the methods proposed in the paper and briefly introduces follow-up works.

II. PHASE-TYPE DISTRIBUTION
PH distribution is the basic mathematical distribution of the modeling method proposed in this article. This distribution is dense on the positive semi-axis; that is, it can fit any form of probability distribution with random variables located on the interval (0, +∞). In addition, PH distribution also has the characteristics of closed operation in the convolution, which can greatly simplify the calculation [25]. Concepts and properties of continuous and discrete PH distributions are briefly introduced here.

A. CONTINUOUS PH DISTRIBUTION
Consider a Markov process {X (t) , t ≥ 0} with continuous time and discrete states, where X (t) represents the state in which the random process is located at moment t. Assume that the states of the Markov process are {1, 2, . . . , k + 1}, where k + 1 is an absorbing state and the distribution vector of the initial state is σ = α 1 . . . α k 0 , which satisfies

The infinitesimal generator of the Markov process is
j = 0. Based on the Markov process described above, a k-th order continuous PH distribution (α, D 0 ) of time T taken from any initial state i = 1, 2, . . . , k to the absorbing state k + 1 can be defined [9], where α = α 1 . . . α k represents a probability distribution vector of the initial phase, D 0 represents a phase-type generator, and d 1 represents an exit vector.
The CDF of the PH distribution is: where e represents a column vector with appropriate dimensions whose elements are all 1. The PDF of the PH distribution is: The moments of the PH distribution are:

B. DISCRETE PH DISTRIBUTION
Similar to the continuous PH distribution, the discrete PH distribution represents the probability distribution of transition steps N = 0, 1, 2, . . . before entering the absorbing state. The s-th order discrete PH distribution can be expressed as (β, S 0 ) [9], where β = β 1 . . . β s represents a probability distribution vector of the initial phases, and β s+1 represents an absorption probability at the initial moment, which satisfies represents a transition probability matrix.
Before entering the absorbing state, the probability distribution of transition steps N is: where n = 0, or where n = 1, 2, . . ..

III. RELIABILITY MODEL BASED ON PH DISTRIBUTION
In practical situation, the number of shocks is often relatively easy to obtain, while the distribution of the cumulative shock magnitude and the failure threshold are usually difficult to determine. At this time, both continuous and discrete PH distributions can be combined to establish the reliability model.

A. MODEL ASSUMPTIONS AND ANALYSIS
Assume that the arrival times of shocks are independent from their magnitudes, and are also independent from the cumulative magnitude of the damage that has already been caused. The distribution and specific parameters of the arrival times and the number of shocks are described as follows.

1) DISTRIBUTION OF THE ARRIVAL TIME
Assume that the arrival time of the n-th (n = 1, 2, . . .) shock is T n > 0, while T 0 = 0, and inter-arrival time T n = T n −T n−1 > 0 is an i.i.d. random variable with CDF G (t) = P ( T n ≤ t) , t > 0, which satisfies a PH distribution (α, D 0 ) with m-th order. Wherein, α = α 1 . . . α m represents the probability distribution vector of the initial phase, which satisfies m i=1 α i = 1 and The corresponding PDF of T n is: Assume that the CDF of the n-th arrival time T n is G n (t) = P (T n ≤ t) , t > 0, which could be expressed by a new PH distribution µ n , M (0) n based on the combination of (α, D 0 ) [11]. Wherein, µ n = α 0 . . . 0 1,mn represents the probability distribution vector of the initial phase, and the subscripts respectively represent the number of rows and columns of the combined vector or matrix.
represents the phase-type generator where represents an absorption rate matrix, which is also a matrix describing the transition rate between phases when a shock arrives. The corresponding exit vector is m Then the CDF of T n is: n t e, t > 0 (9) The corresponding PDF of T n is:

2) DISTRIBUTION OF THE SHOCK NUMBER
Assume that the number of shocks N that components can withstand satisfies the discrete PH distribution (β, P 0 ) with order υ, where β = β 1 . . . β υ represents a probability distribution vector of the initial phase, β υ+1 represents an initial probability of the absorbing state, and also satisfies represents a transition probability matrix between phases.

B. RELIABILITY FUNCTION OF THE SHOCK MODEL
If the entire failure process caused by shocks is regarded as a Markov process with continuous time and discrete states, the infinitesimal generator Q of the process can be constructed based on the above mentioned PH distributions. The following discussions are divided into two cases: (1) Phase transfer when NO shock Arrives This means that the phase in distribution (α, D 0 ) has not transferred into the absorbing state. At this time, the transition-rate matrix between two adjacent shocks is D 0 .
(2) Phase transfer when A shock Arrives This means that the phase in distribution (α, D 0 ) has already transferred into the absorbing state, then has left and transferred to another phase immediately, and the corresponding transition-rate matrix is D 1 .
The shock arrives when the phase in distribution (α, D 0 ) transfers to the absorbing state, and the shock results depend on the phase transition in distribution (β, P 0 ). The discussions can be divided into the following cases: (1) No failure caused by a shock x The failure does not occur under both the (n-1)-th and the n-th shock, and the corresponding transition-rate matrix is where ⊗ represents the Kronecker product [26] of the two matrices.
y The failure does not occur under the first shock and the corresponding transition-rate matrix is β ⊗ D 1 . This means that the initial phase in distribution (β, P 0 ) is not an absorbing state.
(2) Failure caused by a shock x The failure does not occur under the (n-1)-th shock, and occurs under the n-th shock, then the corresponding transition-rate matrix is p 1 ⊗ D 1 .
y The failure occurs under the first shock, and the corresponding transition-rate matrix is β υ+1 ⊗D 1 . This means that the initial phase of the distribution (β, P 0 ) is the absorbing state. The phase transition-rate matrix in different cases can be obtained in TABLE 1. The infinitesimal generator Q for the entire shock process is: Then the life distribution satisfies a new continuous PH distribution (γ , E 0 ) with order (υ + 1) m. Wherein, the probability distribution vector of the initial phase is: and the transition-rate matrix E 0 can be found in TABLE 1. The corresponding exit vector is: The CDF of life at time t is: The PDF of life at time t is: Thus, the reliability function of the shock model can be expressed as:

IV. METHOD OF PARAMETER ESTIMATION
The form of the likelihood function is different with different types of observed sample data. Eryilmaz [27] considered shock models with two types of observed data and used PH distribution to obtain the survival probabilities and MTTF of a system. In [28], a shock model with only interval data was proposed and an approximate condition was assumed aiming at the parameter estimation of PH distribution.
Assuming that the samples are independent of each other, and the data can only obtain at censored time t (c) . The information of observed sample data is as follows: (1) The number of the samples is S > 0.
(2) The set of the failure modes of the samples is δ = δ 1 . . . δ S . Wherein, δ s = 0 represents the censored sample, and δ s = 1 represents failure sample. Based on the above sample data, parameters of the model can be estimated using the EM method, and the expressions of the parameter estimations can then be obtained.
According to the assumption, the arrival times of shocks are independent from their magnitudes, so the parameters in distribution (α, D 0 ) and distribution (β, P 0 ) can be estimated separately. Each of these are described in the following sections.  where e n, represents a column vector with appropriate dimensions, whose (n-1)m + 1-th to nm-th elements are 1 and the remaining elements are all 0. Then, the likelihood function for the total samples is: According to the EM method, the likelihood function with complete data needs to be constructed. Before that, three parameters and one function are required to be assumed.
(3) Let n (s) n represent the number of phase transitions of sample s = 1, 2, . . . , S in the n-th (n = 1, 2, . . . , n s ) inter-arrival time period, including the phase transition when the shock arrives.
(4) Let function Based on the above assumptions, four statistics can be constructed as follows. Notice that the value of n appearing in the following is based on the type of the sample which is shown in TABLE 3.
represents the set of statistics of sample s = 1, 2, . . . , S. The likelihood function with complete data is: wherein, the likelihood function of the censored sample (δ s = 0) is: and the likelihood function of the failure sample (δ s = 1) is: The posterior distribution of H for all the samples can be expressed as: The log-likelihood function can be obtained from function (19): It is assumed that the parameters θ 0 have been obtained through a certain step of iteration using EM method. Then, the expectation with the posterior distribution for H in the log-likelihood function can be expressed as: wherein, the expectation of the censored sample (δ s = 0) is function (25) and the expectation of the failure sample (δ s = 1) is: The parameters in θ to be estimated satisfy the following constraints: (1) The probabilities of the initial phases satisfy m j=1 α i = 1; i = 0 for any phase i = 1, 2, . . . , m in the phase-type generator.
Then parameters in θ can be expressed as follows (TABLE 4). The calculation method of the above expectations is based on improving the method mentioned in [29], [30], and the expressions are directly given here: n,i is: where i = 1, 2, . . . , m and n = 1, or function (28), as shown at the bottom of the next page, where i = 1, 2, . . . , m and n = 2, 3, . . . , n s + 1.
(2) The expectation of statistics Z (s) n,i is function (29), as shown at the bottom of the next page, where i = 1, 2, . . . , m, n = 1, 2, . . . , n s +1 and e n,i represents a column vector with appropriate dimensions whose (n − 1)m + i-th element is 1, and the remaining elements are all 0.
(   The EM method is still used here to estimate parameters in discrete PH distribution (β, P 0 ). A parameter and a function are required to be assumed first.
(1) Let p Based on the above assumptions, four statistics can also be constructed as follows.
wherein, the likelihood function of the censored sample (δ s = 0) is: and the likelihood function of the failure sample (δ s = 1) is: The log-likelihood function can be obtained from function (38): The posterior distribution of H for all the samples is: It is assumed that the parameters θ 0 = β (θ 0 ) , P have been obtained through a certain iteration using EM method. Then, the expectation with the posterior distribution for H in the log-likelihood function can be expressed as: wherein, the expectation of the censored sample (δ s = 0) is: and the expectation of the failure sample (δ s = 1) is: The parameters θ to be estimated satisfy the following constraints: (1) The probabilities of the initial phases satisfy  Then the parameters can be expressed as follows (TABLE 5) where i, j = 1, 2, . . . , υ.
Like the estimation method of the continuous PH distribution, the following statistics need to be bought into the expressions in TABLE 5 to obtain the new iteration values of the parameters. VOLUME 8, 2020 The expectations of the statistics with posterior distributions of the samples are: ij |n s , δ s , and E θ 0 Y (s) i |n s , δ s . The calculation method of above expectations is mainly based on improving the method proposed in [31], and the results are given directly here. is: where i = 1, 2, . . . , υ and n s = 1 or where i = 1, 2, . . . , υ and n s = 1.
In all of the above expressions in this section, parameters in θ 0 are expressed as follows.
The probability distribution vector of the initial phase is: The transition probability matrix is: The exit vector is:

C. IDENTIFIABILITY
When facing certain conditions, the parameters of the proposed model will be unidentifiable. This section will discuss these conditions.

1) CENSORED SAMPLE δ s = 0
It is not difficult to see that the statistics Y (s) i and Z (s) are not involved in the estimation of censored samples. Therefore, if a sample set only contains censored samples, the parameters cannot be identified.
2) FAILURE SAMPLE δ s = 1 (1) When n s = 1 does not exist in the set of failure samples.
It can be known from above analysis that if n s = 1 does not exist in the set of failure samples, then β υ+1 = 0. In fact, one shock is unlikely to cause failure. If it is possible to cause failure by only one shock in actual environment, since the failure sample data does not include this, the system will be unidentifiable based on the obtained samples. If the failure caused by one shock is an event with very small probabilities, then it can be assumed that υ i=1 β i = 1 and β υ+1 = 0, and the system will remain identifiable.
(2) Only n s = 1 exists in the set of failure samples. When all the failure samples are failed by one shock, according to the analysis above, it can be known that At this time, the system is also unidentifiable. Excluding the above conditions, the system will be identifiable.

V. CASE SIMULATION AND ANALYSIS
Generally speaking, a real data set is often used to prove the validity of a method. However, in this work, the original distribution of any real data set (if it really is a real data set) could not be found easily, meaning the lack of an accurate control group could not prove the validity of the proposed method. Besides, the characteristics of the PH distribution show that it can match any distribution, which means the PH distribution will not be limited by any specific distribution cluster. Thus, three different types of simulation cases are proposed and simulation results are used to analyze the adaptability of the proposed model.
Obviously, it is unlikely to obtain the original reliability function or its PDF based on function (1), especially with assumptions of general distributions, and the censored samples also limit the length of the time interval on X-axis. Considering realistic operability and computational complexity, two steps are considered to analyze the adaptability of the proposed modeling method.
(1) A relatively small group of the samples (include censored samples) is used to estimate the parameters; (2) Monte-Carlo simulation method is then used in which censored time is ignored and the sample size is up to a large number to verify the adaptability of the estimated results.
This section assumes that inter-arrival times between two adjacent shocks of all cases are i.i.d. random variables, and are independent from the single shock magnitude.

A. I.I.D. SHOCK MAGNITUDE
This case simulated a general shock situation, in which the single shock magnitude and the inter-arrival time were both i.i.d. random variables.

1) MODEL ASSUMPTIONS
(1) Assume that the sample size used to estimate parameters is S = 500, the sample size used for Monte-Carlo simulation is S MC = 50000, the failure threshold is fixed, satisfying S th = 50, and the censored time is t (c) = 70. (2) Assume that the single shock magnitude satisfies a gamma distribution with shape parameter 2 and scale parameter 2.
(3) Assume that the inter-arrival time between shocks also satisfies the gamma distribution with shape parameter 3 and scale parameter 2.

2) PARAMETER ESTIMATION RESULTS
The distribution of the inter-arrival time is fitted by a 4-th order Coxian distribution, and the number of shocks is fitted by a 5-th order discrete PH distribution.
(2) Phase-type generator:  The estimated PDF of the inter-arrival time based on PH distribution and its statistical characteristics compared with the assumed distribution are shown in FIGURE 2 and TABLE 6.  The initial value of the PH distribution (number of shocks): (1) Initial distribution vector: β = [ 1 0 0 0 0 ].

3) COMPARISON OF THE MODEL PARAMETERS
The estimated PDF of failure and the reliability function compared with the Monte-Carlo simulation samples are shown in

1) MODEL ASSUMPTIONS
Poisson process is a kind of stationary stochastic processes with a constant intensity λ > 0. When the intensity λ (t) = β η t η β−1 , t > 0 becomes a time-dependent function, the process obtained is called Weibull process [32]. Wherein, β > 0 represents the shape parameter and η > 0 represents the scale parameter.
The conditional CDF of X n is: When β = 1, it is obviously that the statistical characteristics of the magnitude X n change with the increasing of accumulated magnitude X * n−1 . Assume that the single shock magnitude satisfies a Weibull process with shape parameter 2 and scale parameter 2. The remaining parameters are assumed to be the same as those in Section V(A).

2) PARAMETER ESTIMATION RESULTS
The inter-arrival time is fitted by an Erlang-4 distribution, and the number of shocks is fitted by a 4-th order discrete PH distribution.
(2) Phase-type generator: (3) Exit vector: The estimated PDF of the inter-arrival time based on PH distribution and its statistical characteristics compared with the assumed distribution are shown in FIGURE 5 and TABLE 8.

3) COMPARISON OF MODEL PARAMETERS
The estimated PDF of failure and the reliability function compared with the Monte-Carlo simulation samples are shown in FIGURE 6 and FIGURE 7, respectively. The comparison of reliability at different time points are shown in TABLE 9.

C. RANDOM FAILURE THRESHOLD
This case simulated a situation in which the threshold was uncertain.

1) MODEL ASSUMPTIONS
Assume that the failure threshold is uniformly distributed in the interval [50, 60]. The remaining parameters are assumed to be the same as those in Section V(A).

2) PARAMETER ESTIMATION RESULTS
The inter-arrival time is fitted by a 3rd hypo-exponential distribution, and the number of shocks is fitted by a 6-th order discrete PH distribution; The initial value of the continuous PH distribution (interarrival time): (1) Initial distribution vector: α = [ 1 0 0 ].

3) COMPARISON OF MODEL PARAMETERS
The estimated PDF of failure and the reliability function compared with the Monte-Carlo simulation samples are shown in

D. CASES ANALYSIS
The above cases include different types of situations where the failure threshold and the single shock magnitude may be unknown. Based on the statistical characteristics, these cases can be divided into two categories: stationary and non-stationary. The stationary kind refers to the statistical characteristics of the sample, which does not change or changes slowly with time increasing; that is, the samples in Sections V(A) and V(C), while the samples in Section V(B) are non-stationary. In all cases, the same gamma distribution was used for simulation several times. However, in order to show the characteristics of parameter estimation based on PH distribution, three different types of PH distributions were used for fitting. Since the arrival time is independent from its magnitude, the differences in statistical characteristics among these three types of PH distributions can still be seen from the estimation results in the above tables from Section V(A) to V(C). From the above cases, it can also be clearly seen that the fitting condition of the 3rd order hypo-exponential distribution is the best. This is because a gamma distribution with shape parameter 3 and scale parameter 2 can be transformed into PH distribution as follows: (1) Initial distribution vector: α = [ 1 0 0 ].
(2) Phase-type generator: (3) Exit vector: The parameter estimation results with different types indicate that the selection of the initial form of PH distribution is not arbitrary. For general distributions, the estimated error is affected by the number and randomness of samples, the type of PH distribution and its initial value. The order (number of phases) should adequately represent the failure time of samples and the computational effort should be acceptable. Suggestions for selecting a PH distribution with suitable form and order are provided in [33]- [35], and PH distribution with different forms has been described in detail in [36].
In practice, the estimated reliability is usually conservative, which means the estimated reliability function under stationary conditions has a longer effective interval than non-stationary condition. Due to the stationary and non-stationary conditions of the single shock magnitude, it can be seen that when fitting through a general form of discrete PH distribution, the fitting errors under non-stationary condition increases over time. However, the fitting results are not unacceptable. Thus, the model proposed in this article has good adaptability.
A key to increase the length of the effective interval and reduce the error of the estimated results under non-stationary condition is to construct discrete PH distribution with a specific form. For the situation mentioned in Section V(B), the degradation process can be divided into different stages along with different discrete PH distributions to describe them respectively. One of the specific forms is as follows: Supposing that β (i) , P is the transition probability matrix between phases, and P (i) 1 = p (i) 1 β (i+1) . The above expansion makes the model more complex for fitting the actual situation, which will lead to a reduction in the reliability function error compared with the general form distributions. The expansion also reduces the scope of application of the model. This will be discussed in another paper. It can be clearly observed that the extension still does not cast off the modeling method proposed in this article, which means the method has good expansibility. VOLUME 8, 2020

VI. CONCLUSION
In view of the limitations of the traditional shock model in practical applications. A shock modeling method was proposed with both continuous PH distribution and the discrete PH distribution. The simulation cases clearly illustrate that the established model has good adaptability and expansibility.
It is relatively easy to use PH distribution to solve the convolution problem that general distribution may require a large amount of computation. Subsequent work will continue to be based on PH distribution, and will consider degradation modeling methods using more complex environments, including various shocks and wears.
XU LUO received the M.S. and Ph.D. degrees in mechanical engineering from the National University of Defense Technology (NUDT), Changsha, China, in 2008 and 2014, respectively.
He is currently a Lecturer with the Laboratory of Science and Technology on Integrated Logistics Support, NUDT. His current research interests include maintenance engineering, maintainability design, and maintainability evaluation. She is currently an Equipment Tester. Her current research interests include random modeling and probability analysis. VOLUME 8, 2020