Dynamic Risk-Adjusted Monitoring of Time Between Events: Applications of NHPP in Pipeline Accident Surveillance

An important focus of Statistical Process Monitoring methods is to develop efficient control charts to monitor the quality of a process over time, and to improve performance by providing valuable insights into process stability and variability (Montgomery, D. C. (2019)). Though the ideas were first introduced to monitor industrial processes (e.g., a manufacturing process), they have now gained a strong foothold in many nonindustrial sectors, including medicine, healthcare and public health, surveillance, financial markets, climate and environment, chemical analysis, network monitoring, and change-point problems, see for example (You and Qiu 2020; Stevens et al. 2021; Bisiotis et al. 2022; Qiu and Xie 2022). In an industrial context, SPM techniques can noticeably contribute to improving the quality of manufacturing processes by monitoring the output of the process over time. As another example in a financial context, SPM tools can be used for various financial applications such as portfolio monitoring and stock trading (Kovarik et al. (2015)).

One important application of SPM is in the context of monitoring defectives or failures, where attribute control charts (e.g. $p$ chart) have been traditionally used. However, for monitoring high-quality processes, where the number of defectives (failures) is often most likely to be zero, a more efficient monitoring procedure is to focus on the elapsed time between two successive defective items (failures), and not on monitoring the number of defectives. Time-between-events (TBE) control charts are widely used for these processes. Note that in this context, although traditionally, the event of interest has been “producing a defective item”, the idea is not limited to this case, and the event of interest can include a broad range of applications, even in non-manufacturing contexts, such as the time between disease outbreaks in a health care context, or the time between failures of operating systems in the industrial sector. To be consistent throughout the paper, “system” may refer to a “manufacture item”, a “production process”, or a “repairable system”. Likewise, a “failure or defective” will be used to refer to any event of interest.

Most of the studies on TBE charts assume that failures occur according to a homogeneous Poisson process (HPP), which results in the distribution of TBE being the well-known exponential distribution (Castagliola et.al 2021; Kumar et al. 2022; Ahmad et al. 2023a, b). This means that the deterioration mechanism of the production equipment, say its hazard rate function, is assumed to be constant over time. While this assumption may apply in some cases, it may not be valid even for non-complex production equipment where deterioration is expected with time. This argument could also be true in other monitoring settings. For example, if the goal is to monitor the real-time condition of a repairable system via monitoring the time between its failures, then the HPP may not be an efficient choice. A repairable system is a system that can be repaired after breaking down, thus its failure rate is most likely to change with time and accordingly the non-homogeneous Poisson process (NHPP) is a more realistic choice to model the number of failures of such a system (Guler et al. 2022). In practice, many systems exhibit complexities and dynamic behaviours, making it challenging to assume a constant failure rate. In such cases, the NHPP offers a more accurate representation of such systems by accounting for varying rates and capturing different failure patterns (Pradhan et al. 2022; Hong et al. 2021). That is why the NHPP plays a central role in modelling the failure counts of complex systems in reliability analysis and maintenance literature (Xu et al. 2017; Cha and Finkelstein 2018; Safaei et al. 2019; Safaei and Taghipour 2022). Despite this, the attention to NHPP has been comparatively low in the SPM literature.

It is essential to recognize that any failure of a system may come with costs to bring it back to operation. These costs may include expenses for repairs and replacements, parts and labour, downtime, warranty claims, as well as potential penalties and fines, either individually or in combination. Since to total cost (TC) associated with each failure may be affected by various factors, it makes more sense to treat it as a random quantity. In the TBE literature, the “amplitude” variable is used as a representative of the random TC. Many techniques have been developed for monitoring the TBE and amplitude/TC variables simultaneously that are also known as the TBEA control charts (Sanusi et al. 2020; Qu et al. 2018; Rahali et al. 2019).

The failure time datasets often include information about some risk factors (also known as covariates or confounding variables) in addition to the time to failures and costs. These factors can noticeably affect the performance of the system. They encompass a wide range of information and ignoring them in the analysis can cause a high degree of variation and overall deterioration in the system’s performance (Zheng et al. 2021; Zhang et al. 2023). There are several models in the literature to incorporate the relevant risk factor information. These include the well-known Cox proportional hazards model, which is more likely to be applied to account for system-reliability-related risk factors in an industrial context, and the accelerated failure time model which is more commonly used to account for health-related risk factors in a clinical context (Steiner and Jones (2010)). Applying a proper risk model to account for the effect of risk factors provides a more informative framework (Cockeran et al. 2021). In the healthcare-related SPM domain, several control charts have been developed that account for the effect of patient’s preoperative risk factors. Among the pioneering works in this area, Steiner et al. (2000, 2001) studied risk-adjusted control charts to track patient death rates. Since then, a huge body of risk-adjusted monitoring methods has been proposed by different scientists (Steiner and Mackay 2014; Woodall et al. 2015; Sachlas et al. 2019). Also, Paynabar et al. (2012) introduced a comprehensive phase I risk-adjusted control chart for monitoring binary surgical outcomes by considering categorical covariates. Steiner (2014) provides a list of justifications for risk adjustment in the health context. For example, they stated that patients may have different conditions before treatment and thus are not expected to be homogeneous (unlike manufactured parts). One can argue the same is true in the case of monitoring the performance of systems in the presence of different risk factors such as material, labour, weather conditions, and so forth. Risk adjustment in this case takes care to distinguish the failure mechanism of the system itself from the risk model. By doing this, the monitoring method would be more sensitive to the changes of system failures and associated costs (Tian et al. 2021; Mun et al. 2021).

The main objective of this paper is to develop a risk-adjusted control chart for monitoring the average cost per time unit, defined as the ratio TC/TBE, using the NHPP as the underlying stochastic process for the failure counts. The adoption of the NHPP model offers a crucial advantage, as it allows the failure intensity of the system to change with time. Additionally, the incorporation of the total cost and risk factors directly into the control chart analysis represents a novel approach to monitoring procedures. The rest of the paper is structured as follows: Section 2 presents a motivational example based on pipeline accident data that forms the background of this work. Section 1 briefly introduces the NHPP with two widely used intensity models, as well as discusses the risk adjustment approach. Section 4 offers an overview of the basic principles and background of copula theory. Moving on, Section 5 focuses on the construction of the proposed control chart, while Section 6 presents a numerical study that evaluates the performance of the charts. In Section 7, a practical application of the proposed charts based on real data examples is discussed. Finally, the study is concluded in Section 8 with a concise summary of the findings, and conclusions drawn.

The Pipeline and Hazardous Materials Safety Administration (PHMSA) of the United States Department of Transportation has released a comprehensive report, including a dataset on oil pipeline accidents between January 1, 2010, and September 1, 2017. The accident/failure of oil pipelines refers to an incident in which the pipeline leaks or ruptures, releasing the contents of the pipeline. Such failures can be hazardous and can potentially cause damage to the environment and human health as well as economic losses (Lu et al. 2023; Xi et al. 2023; Wang et al. 2010; Li et al. 2016). The financial consequences of pipeline accidents are also substantial (White et al. 2003). The total cost of oil spill response and cleanup operations globally ranges from tens of millions to billions of dollars annually (Etkin. 2004). That is why it is crucial to provide efficient strategies to inspect the pipelines, monitor their potential failures, and develop emergency response plans in case of a failure. All these kinds of efforts could be highly beneficial to minimize the risk and cost of failures, as well as to prevent any unexpected failures. By closely monitoring pipeline accidents, operators can detect potential issues at an early stage and take corrective and preventive actions before they either escalate or happen again in the future.

The database comprises details of 2,789 accidents and is publicly available on the Kaggle website (https://www.kaggle.com/datasets/usdot/pipeline-accidents). Various costs (e.g., property damage costs, lost commodity costs, emergency damage costs, etc.) were incurred in order to repair the pipeline and bring it back into the network. In addition, the dataset provides the TC associated with each accident, which is the aggregate sum of all incurred costs. Furthermore, at the time of each accident, information regarding the pipeline type and location, the type of hazardous liquid involved, the accident location, and the accident cause has also been recorded. All these risk factors could potentially affect the intensity or pattern of the accidents.

A number of articles on the analysis of pipeline failures show that an NHPP provides a satisfactory model for pipeline incidents, for example (Cobanoglu et al. 2016; Kermanshachi et al. 2020). Accordingly, an NHPP that accounts for a time-dependent accident rate and adjusts for risk factors is expected to provide a better statistical model for the accident mechanism of pipelines. In this context, there is a need to develop a monitoring technique based on a risk-adjusted NHPP model. Thus, the proposed monitoring approach can serve as a real-time surveillance system to detect statistically and practically significant (unusual) deviations (increases and decreases) in the times between successive accidents and the associated total cost from the baseline of (expected) accident patterns. Accordingly, the proposed methodology can help in dealing with and mitigating these undesirable, and unsafe situations and manage the associated risks and costs.

In what follows, we present additional information about the dataset to provide a deeper understanding. Note that the majority of these incidents are concentrated in the central and southern regions of the United States (US). This geographical distribution of accidents across the US is depicted in Figure 1. For instance, the state of Texas has reported 1004 accidents, whereas the number of accidents is considerably lower in the northwest of the US, such as in the state of Oregon state, where only four accidents have been recorded (see Figure 2). According to Figure 2, approximately 50$\%$ of accidents occur in pipelines that are used for transporting crude oil. However, when considering the location of the pipelines, 99$\%$ of the accidents happen with onshore pipelines, with 53$\%$ of them being above ground. Additionally, about 51$\%$ of accidents are attributed to material, welding, or equipment failures.

Figure 3 (a) and (b) show the letter-value plots for the TBE data (observations) (ranging from 0 to 250 hours) and the TC data (ranging from 0 to 50 million USD). The letter-value plot provides a more accurate visualization compared to the traditional box plots for large and highly skewed data sets (Hofmann et al. (2017)), which is indeed the case in this pipeline accidents dataset. Based on these figures, it’s evident that the majority of TBE data points fall within the interval of $(0,25)$ hours, indicating that most failures occur in close succession. Additionally, a significant portion of the TC data is concentrated in the range of 0 to one million dollars, suggesting that most incidents are associated with relatively lower costs. Both these figures suggest right-skewed distributions for both TBE and TC variables while showing a small number of outliers in the dataset.

The NHPP is a counting process that allows the failure rate/intensity to change with time. Denoting an NHPP with $N(t)$ for $t>0$, where $N(t)$ represents the number of failures that occurred in the time interval $(0,t]$. The process $N(t)$ can be fully characterized by the intensity function $\lambda(t)$ that quantifies the failure rate at time $t$. The NHPP satisfies the following conditions:

• •

  $N(0)=0$, which means that no failures occur at time $t=0$.

• •

  The process has independent increments which means that the number of events in disjoint intervals $\left(s_{1},s_{2}\right]$ and $\left(t_{1},t_{2}\right]$ are independent.

• •

  At most one failure can occur in any infinitesimal time interval, i.e., $P\left(N(t+\Delta t)-N(t)\geq 2\right)=0$.

• •

  The probability that a failure occurs in the time interval $[t,t+\Delta t]$ is given by $P(N(t+\Delta t)-N(t)=1)=\lambda(t)\Delta t+o(\Delta t)$.

It follows that in an NHPP, the number of failures in the time interval $(s,t]$ has a Poisson distribution with parameter $\Lambda(t)-\Lambda(s)$ where $\Lambda(t)=E(N(t))=\int_{0}^{t}\lambda(u)du$ is called the mean function (also known as cumulative intensity function) of the process. Thus, we have:

According to equation (1), it can be inferred that $N(t)$ follows a Poisson distribution with the parameter $\Lambda(t)$. When the intensity function is constant, i.e., $\lambda(t)=\gamma$, say, the NHPP reduces to HPP with the constant rate $\gamma$. The intensity function $\lambda(t)$ can take different mathematical forms. However, in many applications, the log-linear and the power-law models are the most widely used. The following two subsections provide a brief introduction to these models.

Copula modelling has emerged as a widely adopted technique for capturing dependencies in various domains of application. It offers the ability to disentangle the dependence structure from the joint distribution function of a set of variables while preserving their individual univariate marginals. Sklar’s theorem plays a pivotal role in this process, allowing us to establish the joint distribution $H(x,y)$ of two random variables, namely $X$ and $Y$, with marginal cumulative distribution functions $F_{X}(.|\bm{\theta}_{X})$ and $F_{Y}(.|\bm{\theta}_{Y})$, respectively. The equation for $H(x,y|\bm{\theta}_{X},\bm{\theta}_{Y},C)$ involves the copula model $C$ once its parametric form is determined, where:

where $\bm{\theta}_{X}$ ($\bm{\theta}_{Y}$) is the vector of $X$ ($Y$) distribution’s parameters and $\theta_{c}$ is the parameter of the copula function $C$. A key advantage of Sklar’s theorem lies in its guarantee that a unique copula $C$ exists when dealing with continuous marginal distributions, as is the case in our study. Before delving into specific copula models, it is essential to introduce the dependence measure known as Kendall’s tau $\tau$.

Recall that the main aim of this study is to develop a control chart to monitor the TBE variable and its associated TC while accounting for the risk factors recorded at the time of each failure. In this setting, one reasonable quantity for monitoring the stability of such a process, that incorporates the effects of these two factors, is the average cost per unit time (AC) so that AC=TC/TBE. The average cost per unit time metric also has a meaningful interpretation, making the control chart interpretation more practically meaningful and providing valuable insights into the financial implications. Both of these are relevant aspects of the decision-making process.

We assume that the system starts to operate at time $T_{0}=0$. As time passes, failures may occur at random times $T_{1},T_{2},...$ resulting in TBEs $X_{1}=T_{1},X_{2}=T_{2}-T_{1},...$. In addition, at time $t_{i}$ of failure $i=1,2,\ldots$, data on the risk factors and the TC variable denoted by $\bm{z}_{i1},...,\bm{z}_{ip}$ and $Y_{i}$, respectively, are collected. Figure 4 shows a schematic of this data collection process.

To calculate the control limits of the proposed monitoring scheme, we first derive the conditional cumulative distribution function (CDF) of the $i$th failure time $T_{i}$ given the vector of risk factors $\textbf{{z}}_{i}$ and the fact that $T_{i-1}=t_{i-1}$ as:

where $\bm{\theta}_{X}=(\gamma,\eta,\bm{\beta})$ is the vector of parameters regarding the intensity function and risk model and $\Lambda\left(.|\bm{z}_{i}\right)$ is given in (7). The second equality in (5) is obtained using the independent increments of NHPP and the third one is obtained using the fact that:

Eventually, the conditional CDF of $X_{i}$ given $T_{i-1}=t_{i-1}$ and $\bm{z}_{i}$ can be obtained using equation (5) as follows:

As a result, the conditional probability density function (PDF) of $X_{i}$ can be calculated as:

Let $Y_{i}$ be the corresponding random total cost of failure $i=1,2,\ldots$ having the CDF $F_{Y_{i}}(.)$ and PDF $f_{Y_{i}}(.)$ with the vector of parameters $\bm{\theta}_{Y}$. Furthermore, let us assume $X_{i}$ and $Y_{i}$ to be dependent and their joint CDF is given by:

where $C(u,v|\theta_{c})$ is a Copula model containing all information on the dependence structure between $Y_{i}$ and $X_{i}$ and $\bm{\theta}=\left(\bm{\theta}_{X},\bm{\theta}_{Y},\theta_{c}\right)$. Accordingly, the joint PDF of $Y_{i}$ and $X_{i}$ can be derived as:

where $f_{\left(X_{i}|t_{i-1},\bm{z_{i}}\right)}(.|\bm{\theta}_{X})$ is given in (15), $f_{Y_{i}}(.|\bm{\theta}_{Y})$ is the PDF of $Y_{i}$, and $c(u,v\mid\theta_{c})=\frac{\partial C(u,v\mid\theta_{c})}{\partial u\partial v}$ is the copula density.

In order to assess the performance of the proposed control chart and investigate its sensitivity with respect to various factors, this section aims to conduct a numerical study based on the $ARL$ metric. First of all, let us set $ARL_{0}=200$. In addition, consider a system where its failures come from an NHPP with power law or log-linear intensity function. The IC parameters are assumed to be $\eta_{0}=1.50$ for the power law model, $\eta_{0}=2.00$ for the log-linear model, and $\gamma_{0}=0.05$. These choices of intensity parameters represent increasing failure rates, a scenario more probable in deteriorating systems. Furthermore, we assume the associated TC represented by $Y$ follows an Exponential distribution with IC rate $\mu_{0}=1.00$. The joint distribution of $Y$ and $X$ is also characterized by a Gumbel copula model. Furthermore, we introduce a risk factor $z$ ($p=1$), which impacts the TBE. The risk factor will be a factor ranging from 1 to 3, generated using a truncated Poisson distribution with rate $1$, which means the numbers 1, 2 and 3 will be generated with probabilities $0.36,0.18$ and $0.06$, respectively. In addition, the regression coefficient $\beta$ takes values from the set $\{-2.00,-0.50,0.00,0.50,2.00\}$ to cover different association levels between the risk factor and TBE. To capture different levels of dependence between the TBE and the TC when they are positively associated, we consider the dependence parameter $\tau$ from the set $\{0.3,0.8\}$. This will enable us to explore moderate and strong positive dependencies.

To conduct the numerical analysis, we assume the parameters $\gamma$ and $\eta$ regarding the TBE variable and the parameter $\mu$ of the TC variable are subject to shift, while $\beta$ and $\theta_{c}/\tau$ are assumed to remain unchanged. Let us donete the vector of IC parameters by $\bm{\theta}_{0}=(\gamma_{0},\eta_{0},\mu_{0})$. To calculate the $ARL_{1}$ values, the vector of OC parameters is also defined as $\bm{\theta}_{1}=(\gamma_{1},\eta_{1},\mu_{1})=(\delta_{\gamma}\gamma_{0},\delta_{\eta}\eta_{0},\delta_{\mu}\mu_{0})$ where $\delta>0$ quantifies the shift magnitude in each parameter. We set $\delta_{k}=0.25,0.75,1.25,2.00$ for $k=\gamma,\eta$ and $\mu$ to account for both upward and downward shifts in the associated parameters. By examining the $ARL_{1}$ values in these cases, we can gain insights into the performance and effectiveness of the control charts under various scenarios and deviations from the stable condition.

$ARL$ values are computed through Monte Carlo simulations described in Algorithm 1. Though $ARL_{0}$ can be calculated theoretically by $ARL_{0}=\frac{1}{\alpha}$ from (5.1), it could also be calculated from Algorithm 1 by generating $X_{i}$ and $Y_{i}$ from $\bm{\theta}_{0}$. Moreover, $ARL_{1}$ can be computed from Algorithm 1 by generating $X_{i}$ and $Y_{i}$ using $\bm{\theta}_{1}$. Tables 1 and 2 show the results of the ARL analysis for power-law and log-linear intensity models, respectively. The primary result shown by the tables is that the proposed chart is $ARL$-unbiased across all shift combinations and both intensity models, i.e., we always have $ARL_{0}>ARL_{1}$.