Reliability Analysis and Optimal Replacement Policy for Systems with Generalized Pólya Censored δ Shock Model

: A fresh censored δ shock model is investigated. The arrival of random shocks follows a generalized Pólya process, and the failure mechanism of the system occurs based on the censored δ shock model. The generalized Pólya process is used for modeling because the generalized Pólya process has excellent properties, including the homogeneous Poisson process, the non-homogeneous Poisson process


Introduction
Shock models are widely used in practical life and have significant implications in the study of engineering reliability.For example, earthquakes, sudden epidemics, etc., can be considered external shocks.Esary et al. [1] pioneered classical shock models, in which the lifetime distributions of systems suffering from shocks were derived and some lifetime distribution properties were investigated.Shortly afterward, two conventional shock models, the cumulative shock model [2] and the extreme shock model [3] were put forward.Subsequently, Mallor et al. [4] extended the models of [3] by proposing a relaxed run shock model.The shock effects produced by the three shock models described above all depend on the shock magnitude.These three shock models have also been intensively studied in the last three decades.For example, Shanthikumar et al. [5] investigated further general shock models by assuming that the pairs of shock interval and shock magnitude follow a correlated renewal process.More extensions can be seen in [6][7][8].
In addition, two special shock models have also been put forward successively, i.e., the δ shock model and the censored δ shock model.The failure mechanism in both models depends on the inter-arrival time between successive shocks and not on the shock magnitude produced by the shocks.The δ shock model is based on the traffic problems proposed by Li et al. [9].For the δ shock model, Wang et al. [10] considered a mixed shock model, in which the system fails when the magnitude of a shock surpasses a defined level γ or the time inter-arrival between consecutive shocks is no greater than the limit δ, whichever occurs first.Eryilmaz [11] studied a generalized δ shock model by considering the run shock model and δ shock model simultaneously.Bian et al. [12] studied reliability analysis for systems subject to mutually dependent competing failure processes under the δ shock model.Additional relevant research works focusing on the systems subjected to the shock model can be found in [13][14][15][16].
Later, the δ shock model was extended, and the censored δ shock model proposed by Ma et al. [17] was based on the customer lifetime value.The censored δ shock model has crucial practical implications in customer relationship management (CRM).In [17], the lifetime of customers was characterized utilizing the censored δ shock model; the customer's lifetime was defined as the sum of all previous trade times and δ, and the probabilistic properties of the lifetime were also derived.Based on the above background, the failure mechanism of the censored δ shock model can be defined, i.e., the system fails when the time interval between two successive shocks does not reach the given threshold value δ.Next, a model comparison between the δ shock model and the censored δ shock model is given in Table 1, which contains the parameters of both models, the failure mechanism, and the lifetime of the system.The censored δ shock model can also be applied to medical treatment, insurance, and other aspects.There are relatively few studies on the censored δ shock model.In the past decade, Eryilmaz et al. [18] investigated the lifetime distribution of the censored δ shock model assuming that the external shocks arrive according to a renewal process.Bai et al. [19] studied the parameter estimations of the censored δ shock model in uniform intervals.Bian et al. [20,21] considered the lifetime distribution of two types of discrete censored δ shock models under the Markov shock arrival process.Ma et al. [22] exhibited a different way to find a system's reliability with the Poisson censored δ shock model by directly calculating probabilities.Lately, Chadjiconstantinidis et al. [23] considered the survival function of a system subjected to the censored δ shock model, in which the inter-arrival times have well-known discrete and continuous distributions.There is no other literature except the above, so the censored δ shock model is a hot topic worth studying.
In the existing methods of reliability modeling, under a random shock environment, the arrival of random shocks is mostly characterized by homogeneous Poisson processes [24] (HPP) and the non-homogeneous Poisson processes [24] (NHPP).It is well-known that the HPP is a typical renewal process where the inter-arrival time of the random shocks follows an exponential distribution, while the NHPP is a stochastic process with a non-constant hazard rate.Although most studies assume that the shock arrival process follows HPP and NHPP only for the sake of mathematical formulation simplicity and mathematical tractability, this is extremely restrictive due to their independent increments.At the same time, this assumption can lead to significant errors in assessing the impact of shocks on the system.Therefore, to avoid this error, Teugels et al. [25] studied the occurrence of external shocks according to a Pólya process.The Pólya process is introduced as a mixed Poisson process in [26]; it has a general dependence increment.It is characterized by a negative binomial distribution of shock arrivals.Moreover, another essential important property is that the shock arrival times are mutually dependent.Based on these advantages, the Pólya process is increasingly used to model the point process of real life.Eryilmaz [27] considered the modeling of a shock model using the Pólya process and investigated the system's reliability.To the authors' knowledge, there is not enough research on shock model modeling using the Pólya process.This is mainly because the Pólya process is considered mutually dependent on the inter-arrival time between the successive shock arrival times.
Later on, Konno [28] further extended the Pólya process and proposed the generalized Pólya process (GPP) within the framework of a non-stationary-type master equation approach.The GPP can be viewed as a further generalization of the NHPP, which possesses neither the independent increments nor the stationary increment property.Only the marginal distribution of the number of events in (0, t] was obtained in Konno [28] based on a differential equation.The GPP processes the Markovian property.The GPP is a different and more general counting process, and the HPP, the NHPP, and the Pólya process can be regarded as special situations of the GPP.Although the GPP has neither independent increments nor stationary increments, it has pleasant properties and can be used in various applications to obtain closed-form results.In addition, one of the most valuable contributions to the field of reliability is that the GPP allows the definition of a new type of maintenance and a different type of failure process.It has been a great inspiration for the development of various different maintenance modeling techniques.Moreover, the GPP can simulate shocks triggered by a correlated system failure.The current literature on shock models with GPP is as follows.Cha [29] expounded the GPP at great length by deriving a variety of properties that can be used in numerous applications.Cha et al. [30] derived and analyzed the corresponding survival and failure rate functions of the extreme shock model under the GPP.Goyal et al. [31] studied the survival function and correlation properties of a history-dependent mixed shock model under the GPP.Goyal et al. [32] investigated the correlation properties of the time-dependent δ shock model under the GPP and studied the optimal replacement strategy of the established model and the associated random properties.Goyal et al. [33] studied a general δ shock model when the recovery time depends on both the arrival times and the magnitudes of shocks.In view of the above discussion, the extension of the existing literature results to the GPP characterization shock model is of great significance for practical application.
Currently, in the literature, we can conclude that all the published research on the censored δ shock model has been carried out under the assumption of the HPP and the NHPP.However, this assumption is overly restrictive and unrealistic to describe various real-life scenarios.Therefore, we aim to study some lifetime distribution properties of the censored δ shock model based on the generalized Pólya shock process, which is highly rare among existing research.Thus, the present paper aims to consider the extension of the shock arrival process from the renewal process to the GPP.This assumption is efficient and mathematically tractable for modeling point processes with correlated increments.Thus, it is necessary to consider the lifetime properties of the censored δ shock model system under the generalized Pólya shock process.Finally, an optimal replacement policy N under the generalized Pólya censored δ shock process is developed.
The rest of the paper is organized as follows.In Section 2, the definition of the censored δ shock model and the generalized Pólya process are introduced.There are some preliminaries, followed by a description of the model.In Section 3, the closed-form reliability function of the system and its corresponding properties under the generalized Pólya censored δ shock model are derived.In Section 4, under the generalized Pólya shock process, a replacement policy N based on the failure number of the system using geometric processes is developed.In Section 5, numerical examples are presented to validate the rationality and effectiveness of the proposed model.Some conclusions are given in Section 6.

Preliminaries and Model Description Model Description
Consider a system operating in a stochastic environment.The system is subjected to external shocks, and the failure mechanism of external shocks occurs according to the censored δ shock model.In this section, some basic assumptions are given, and the definition of the generalized Pólya process, the generalized Pólya censored δ shock model, the failure rate, the increasing (decreasing) mean failure rate type, and some degradation models for the generalized Pólya process are given.
Let T denote the lifetime of a system that has started operation at t = 0.The arrival of external shocks follows a generalized Pólya process {N(t), t > 0} with a set of parameters {λ(t), α, β}, α ≥ 0, β > 0. Let 0 = T 0 < T 1 < T 2 < • • • be a sequence of the arrival times and X m = T m − T m−1 , m = 1, 2, • • • represent the inter-arrival time between the mth and the (m − 1)th shocks.Let X 0 = 0; the system fails if no shock occurs within a δ length time period from the last shock, where δ is the system failure threshold.The lifetime of the system is defined as T = M ∑ i=0 X i + δ, where M represents the number of shocks until the system fails.Based on the above assumptions, the definition of the censored δ shock model is given, i.e., the system fails when the inter-arrival time between two consecutive shocks is more than the time interval threshold δ. Figure 1 shows the trajectory diagram for the censored δ shock model.A system survives in [0, t) if all shocks' inter-arrival times are less than δ.In the following, some definitions of the corresponding models are given.

Definition 3 ([34]
).Let the distribution and the density function of system lifetime T be F(t) and f (t), respectively.Then, the failure rate is Definition 4 ([34]).If − 1 t lnR(t) is an increasing (decreasing) function of t ≥ 0, then it belongs to the increasing (decreasing) mean failure rate type.
According to the above assumptions and definitions, under the generalized Pólya shock process, the relevant reliability indices for the built model are studied in Section 3.

Reliability of a System with a Generalized Pólya Censored δ Shock Model
The reliability function R(t) of the system lifetime is one of the major reliability indices in the study of the system.Therefore, in this section, the reliability of the system subjected to the generalized Pólya process censored δ shock model is analyzed.Based on the established model, we will firstly obtain some reliability indices of the system, such as the system reliability, and the mean lifetime of the system.Then, the failure rate of the system, the upper bound of reliability of the system, and distribution properties are also derived.Before giving the main results above, the definition of the reliability of the system is given as which is very important for assessing the reliability of the lifetime of a system.To obtain the reliability function for the system under SM{[GPP(λ(t), α, β)], CD(δ)}, the following lemma is needed.
(1) The distribution of N(t) is given by In the following theorem, the survival function of a system for the generalized Pólya censored δ shock model is obtained. where The following two special situations of the system's lifetime should also be of concern.
According to the characteristics of the survival function curve in [17], the distribution of lifetime T = δ can be calculated as and the distribution of lifetime T > δ can be derived as It is worth noting that these two special cases include the results of the non-homogeneous Poisson censored δ shock model (SM{[NHPP(λ(t), 1)], CD(δ)}) and Poisson censored δ shock model SM{[HPP(λ)], CD(δ)}, which are given in [17].That is, when λ t = λ(t) and λ t = λ, respectively, the above two results degenerate to where X is an exponential random variable with the parameter λ t = λ(t) and where X is an exponential random variable with the parameter λ t = λ.
Moreover, the following corollary follows from Theorem 1 using Remark 1(b).

Corollary 2. The reliability of SM{[HPP
where From the second equation of Equation (A4), one has That is, In Theorem 1, the calculation of the function f m,δ (t) is very complicated, so the accurate reliability of SM{[GPP(λ(t), α, β)], CD(δ)} is not easy to obtain.However, sometimes, only the upper bound of the reliability is needed in practice.
Next, the upper bound of reliability is given by Theorem 2.
where t δ represents the largest positive integer not greater than t δ .
In the study of shock models, there is another very important quantity to be studied, which is the number of shocks before the system fails, i.e., M. Due to the generality of the intensity function in GPP and the interdependence of the arrival times of shocks, it is quite difficult to obtain explicit expression for the distribution of M. Therefore, for convenience, we assume that the shock arrival process arrives according to SM{[GPP(λ(t) = 1/(b + αt), α, β)], CD(δ)}.Under this assumption, the probability mass function, the survival function, and the mean of M are obtained by Proposition 1.

Proposition 1.
The following results hold.
(a) The probability mass function of M is given by Proof.

Note that
where Λ is a structure random variable with the probability density given by Furthermore, by conditioning on Λ = λ, the GPP reduces to the HPP with the rate λ.Thus, the conditional probability density function of (X Thus, where the last equality follows from the fact that the binomial expansion and (b) The survival function of M is given by (c) The mean of M is given by Next, the remaining reliability indices of the system under the generalized Pólya process are addressed.
Firstly, the mean lifetime of the system with SM{[GPP(λ(t), α, β)], CD(δ)} is considered.Let ET be the mean lifetime of the system.Based on ET = ∞ 0 R(t)dt and Equation (7), the mean lifetime of the system is given by the following theorem.Secondly, when studying a system, one often wants to know how likely it is that the system will fail at some point.Therefore, according to Definition 3 and the established model, the failure rate of the system is studied as follows.
Theorem 4. The failure rate of the system lifetime is The failure rate of the system is derived by Theorem 4. Now, the lifetime distribution type of the system is investigated.Theorem 5. Let F T (t) be the distribution function of system lifetime; then, where w(t) = −βΛ(t) It should be emphasized that the detailed proofs of Theorem 1, Theorem 2, Theorem 4, and Theorem 5 are shown in Appendix A, Appendix B, Appendix C, and Appendix D.

Optimal Replacement Policy under SM{[GPP(λ(t), α, β)], CD(δ)}
In many practical problems, a repaired system is more likely to fail again than the original one.For the censored δ shock model, the shock interval for the system to fail again should be smaller than the threshold δ value of the original system when the repaired system is subjected to continuous shocks.Therefore, in this section, based on the number of failures of the system, a replacement strategy N of the generalized Pólya censored δ shock model is considered.The problem is to choose the optimal replacement policy N * such that the long-run excepted cost per unit time is minimized.
The following assumptions are made for the repairable system under SM{[GPP(λ(t), α, β)], CD(δ)}: (1) A new system is installed at time t = 0.The system is repaired immediately upon its failure.The system is immediately replaced with a new system when the system is observed to fail for the Nth time.(2) The system suffers from external shocks.The arrival of the external shocks follows a generalized Pólya process {N(t), t > 0} with a set of parameters {λ(t), α, β}, α ≥ 0, β > 0. The system fails when the interval time of successive shocks is larger than δ.After the nth repair, the system failure threshold decreases as a n−1 δ, 0 ≤ a ≤ 1.
Denote T 1 as the first failure time of the system and let T n be the operating time of the system after the (n − 1)th repair to the nth failure, where n = 2, 3, • • • .(3) Let Y n represent the repair time of the system after the nth failure, where In particular, when µ = 0, the maintenance time is ignored.(4) The system repair cost per unit time is c 1 , the operating reward rate per unit time is c 2 , the replacement cost is c 3 , and the replacement time is negligible.(5) The GPP, T n , and Y n are independent.
Denote W as a random length of a cycle under the replacement policy N.Then, According to Theorem 3 and Assumption (3), the expected length of the renewal cycle is obtained as Under the replacement policy N, according to the renewal reward theorem, the longrun average cost per unit time C(N) can be computed as

Case Study
Case studies of customer relationship management (CRM) are developed in this section to illustrate the theory obtained in Sections 3 and 4. In CRM, the customer lifetime value (CLV) plays an important role.A customer relationship management tool can be important for almost any kind of business.In [17], the authors utilized the censored δ shock model to characterize the lifetime of customers and define the customers' lifetimes as the sum of all previous trade times and δ.Based on the CRM background, the failure mechanism of the censored δ shock model was defined, i.e., the system was said to fail when the time interval between two successive shocks did not reach the given threshold value δ.Therefore, the customer's lifetime in CRM is modeled using the generalized Pólya censored δ shock.Some research works focusing on the system with Poisson censored δ shock based on CRM can be found elsewhere [17,22].The parameter settings are provided in Table 1, where some parameter values are derived from existing studies and others are assumed.
We first consider CRM to illustrate the established model.The customer operates in a random environment with customs arriving as the generalized Pólya process.Assume that the intensity function of the generalized Pólya process is λ(t) = λ(2 + sin(t)).Then, the mean function of the generalized Pólya process becomes Λ(t) = λ(2t − cos(t) + 1), and the distribution law of N(t) when the parameters α and β are constant is However, for the convenience of simulation, we assume that the generalized Pólya process with the intensity function of λ(t) = λ.Under this assumption, some numerical examples are developed to illustrate the theoretical results obtained in Sections 3 and 4.
Under SM{[GPP(λ(t), α, β)], CD(δ)}, based on Equation (7), the reliability curves of the system under different parameter settings of λ and δ are plotted in Figures 2 and 3. When the values of the parameter λ are fixed as λ = 4, λ = 6, λ = 8, and λ = 10, respectively, the variation curves of the reliability of δ with different values, i.e., δ = 0.1, δ = 0.2, δ = 0.3, and δ = 0.5, are shown in Figure 2. From Figure 2, one can see that the failure parameter δ has a significant influence on the system reliability R(t).In each subgraph of Figure 2, when the value of λ is fixed, the larger the failure threshold δ and the quicker the increase in the system's reliability.This result indicates that a small shock arrival threshold δ ensures a better reliability performance.That is, the larger the value of δ, the less likely the system is to fail.
In addition, the shock failure rate λ also has a great influence on system reliability.When the values of the failure threshold δ are fixed as δ = 0.1, δ = 0.2, δ = 0.3, and δ = 0.5, respectively, the variation curves of the reliability of λ with different values, i.e., λ = 4, λ = 6, λ = 8, and λ = 10, are shown in Figure 3. From Figure 3, one can observe that with the increase in the failure rate λ, the system's reliability decreases under shock.The result is reasonable.This result indicates that a small shock arrival rate λ ensures a better reliability performance.For the replacement model discussed in Section 4, in Table 2, we compute the values of C(N) when the parameters are c1 = 2, c2 = 0.9, c3 = 6, µ = 3, α = 6, β = 10, a = 0.85, and δ = 2. From the results of the numerical simulation, part of the resulting data is selected and made into Table 1.From Table 2, C(29) = −0.9944 is the minimum value of the average cost of long-term operation.That is, the optimal policy is N * = 29.The plot of C(N) as a function of N is plotted in Figure 4.It can be seen from Figure 4 that C(N) first decreases as N increases, and when N = 29, C(29) = −0.9944 is the minimum; when N > 29, C(N) increases again as N increases.Thus, C(29) = −0.9944 is the minimum of the long-run average cost.That is, the optimal policy is N * = 29, and moreover, it can be seen from Figure 4 that the optimal policy uniquely exists.

Conclusions
In this paper, a repairable system with the generalized Pólya censored δ shock model was studied.We have not only given some reliability indices of the system, such as the corresponding reliability function, the upper bound of the reliability function, and the mean lifetime, but also the failure rate.The lifetime distribution class of the system was also proved.The replacement policy for the system using a geometric process and determining the optimal replacement policy by objective function minimization was also considered.Finally, a numerical example was presented to verify the plausibility of the model using a Monte Carlo simulation in Matlab2016.Since the generalized Pólya process contains the homogeneous Poisson process, the non-homogeneous Poisson, and the Pólya process, compared with the existing censored δ shock model, the generalized Pólya censored δ shock model proposed in this paper is an extension of the censored δ shock model.Data Availability Statement: Not applicable.

Acknowledgments:
The authors would like to thank the reviewers for their helpful comments, which improved the presentation of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.The arrival of external shocks M

Abbreviations
The number of shocks until the system fails X i The inter-arrival time between the ith and the (i − 1)th shocks δ The failure threshold of the censored δ shock model r(t) The failure rate of the system f (t) The density function of system lifetime T F(t) The distribution function of system lifetime T R(t) The reliability function of system lifetime T ET The mean lifetime of the system J The Jacobian determinant The repair time of the system after the nth failure, n = Taking Equations ( 5), (A2), and (A4) into Equation (A1), one has Thus, the theorem follows.

Appendix B. Proof of Theorem 2
Proof.Obviously, P(T > t|N(t) < t−δ δ ) = 0.Then, where where By a linear transformation and the density transformation formula, one has From Equations ( 6) and (A8), one has where J is the Jacobian determinant and Taking Equation (A9) into Equation (A7), one has Substituting Equations ( 5) and (A10) into Equation (A6), it can be seen that This completes the proof.

Appendix C. Proof of Theorem 4
Proof.Due to P(T ≤ δ) = 0, we only consider the case of t > δ.
When t > δ, the probability density function of the system lifetime is From the result of Theorem 1, the derivative of the reliability function of the system is obtained as According to Definition 3, substituting Equations ( 7) and (A13), the failure rate of the system is formulated as .

Appendix D. Proof of Theorem 5
Proof.According to Definition 4, let l(t) = − 1 t lnR(t).From Theorem 1, when t ≤ δ, l(t) = 0; when t > δ, one has Now, for the derivative of l(t) with respect to t, we have

Figure 4 .
Figure 4. Average cost C(N) curve with respect to maintenance policy N.

Author Contributions:
Conceptualization, B.P. and Y.Y.; methodology, L.B.; software, B.P.; validation, L.B., B.P. and Y.Y.; formal analysis, L.B.; writing-original draft preparation, L.B. and B.P.; writing-review and editing, Y.Y.; supervision, Y.Y.; project administration, L.B.; funding acquisition, L.B.All authors have read and agreed to the published version of the manuscript.Funding: The research was supported by the Doctoral Scientific Research Foundation of Northwest Normal University (202203101203).

Table 1 .
Model comparison between δ shock model and censored δ shock model.

Table 2 .
The values of average cost C(N).