An application of survival analysis to disruption prediction via Random Forests

TL;DR

This paper integrates survival analysis with machine learning, specifically Random Forests, to improve disruption prediction in tokamak plasmas, providing better warning times and reducing false alarms.

Contribution

It introduces a novel application of survival analysis to disruption prediction, combining it with Random Forest classifiers for enhanced warning time estimation.

Findings

Survival analysis improves disruption warning times.

The method reduces false alarms compared to threshold-based schemes.

Application to Alcator C-Mod data demonstrates practical effectiveness.

Abstract

One of the most pressing challenges facing the fusion community is adequately mitigating or, even better, avoiding disruptions of tokamak plasmas. However, before this can be done, disruptions must first be predicted with sufficient warning time to actuate a response. The established field of survival analysis provides a convenient statistical framework for time-to-event (i.e. time-to-disruption) studies. This paper demonstrates the integration of an existing disruption prediction machine learning algorithm with the Kaplan-Meier estimator of survival probability. Specifically discussed are the implied warning times from binary classification of disruption databases and the interpretation of output signals from Random Forest algorithms trained and tested on these databases. This survival analysis approach is applied to both smooth and noisy test data to highlight important features of the survival and hazard functions. In addition, this method is applied to three Alcator C-Mod plasma discharges and compared to a threshold-based scheme for triggering alarms. In one case, both techniques successfully predict the disruption; although, in another, neither warns of the impending disruption with enough time to mitigate. For the final discharge, the survival analysis approach could avoid the false alarm triggered by the threshold method. Limitations of this analysis and opportunities for future work are also presented.

Figures

Equations

$$0,$$

$$\displaystyle\mathrm{P}(t_{\mathrm{D}}-t<\Delta t|\mathrm{D})=\Bigg{\{}\;$$

$$\mathrm{P}((t+\Delta t>t_{\mathrm{D}})\cap \mathrm{D})=\mathrm{P}(t+\Delta t>t_{\mathrm{D}}|\mathrm{D})\times \mathrm{P}(\mathrm{D}).$$

$$S(t)=\mathrm{P}(T>t).$$

$$h(t)=-\frac{1}{S(t)}\frac{\mathrm{d}S(t)}{\mathrm{d}t}.$$

$$S(t_b)=S(t_a)\exp \left(-{\int }_{t_a}^{t_b}h(t)\mathrm{d}t\right).$$

$$S(t_n)=\prod _{i=0}^n\mathrm{P}(T>t_i|T\ge t_i).$$

$$S(t_n)=\prod _{i=0}^n\left(1-{\mathrm{P}}_{i\to i+1}\right).$$

$$S(t_n|t_i)=\prod _{j=i}^n\left[1-\mathrm{P}(t_{j+1}>t_{\mathrm{D}}|\mathrm{D}){\mathrm{P}}_{\mathrm{D}}(t_j)\right].$$

$$S(t+\Delta t|t)=\prod _{j=0}^n\left[1-{\mathrm{P}}_{\mathrm{D}}(t+j\Delta t_{\mathrm{s}})\frac{\Delta t_{\mathrm{s}}}{\Delta t_{\mathrm{c}}}\right].$$

$${\mathrm{P}}_{\mathrm{D}}(t+j\Delta t_{\mathrm{s}})\approx {\mathrm{P}}_{\mathrm{D}}(t)+\frac{\mathrm{d}{\mathrm{P}}_{\mathrm{D}}}{\mathrm{d}t}\times j\Delta t_{\mathrm{s}},$$

$$S(t+{\tau }_{50}|t)=\mathrm{0.5.}$$

$$\tau ={\int }_0^{\infty }S(t+t^{\prime }|t)\mathrm{d}{t}^{\prime }.$$

$${\tau }_{\mathrm{eff}}={\int }_0^{1-t}S(t+t^{\prime }|t)\mathrm{d}{t}^{\prime },$$

$$\tilde{\tau }={\int }_0^{1}S(t+t^{\prime }|t)\mathrm{d}{t}^{\prime }.$$

$$S(t+t^{\prime }|t)=\prod _{j=0}^n\left[1-{\mathrm{P}}_{\mathrm{D}}\left(t+\frac{j{t}^{\prime }}{n}\right)\frac{t^{\prime }}{n\Delta t_{\mathrm{c}}}\right],$$

$$\tau ={\int }_0^{\infty }S(t+t^{\prime }|t)\mathrm{d}{t}^{\prime }\ge {\int }_0^{T}S(t+t^{\prime }|t)\mathrm{d}{t}^{\prime },$$

$$\tau \ge {\int }_0^T\prod _{j=0}^n\left[1-\frac{t^{\prime }}{n\Delta t_{\mathrm{c}}}\right]\mathrm{d}{t}^{\prime }={\int }_0^T{\left(1-\frac{t^{\prime }}{n\Delta t_{\mathrm{c}}}\right)}^{n+1}\mathrm{d}{t}^{\prime }.$$

$$\tau \ge -\frac{n}{n+2}\Delta t_{\mathrm{c}}{{\left(1-\frac{t^{\prime }}{n\Delta t_{\mathrm{c}}}\right)}^{n+2}|}_0^{n\Delta t_{\mathrm{c}}}=\frac{n}{n+2}\Delta t_{\mathrm{c}}.$$

$$h(t)=-\frac{\mathrm{d}\ln S(t)}{\mathrm{d}t}.$$

$$S(t)=\prod _{j=0}^n\left[1-{\mathrm{P}}_{\mathrm{D}}\left(\frac{jt}{n}\right)\frac{t}{n\Delta t_{\mathrm{c}}}\right],$$

$$\ln S(t)=\sum _{j=0}^n\ln \left[1-{\mathrm{P}}_{\mathrm{D}}\left(\frac{jt}{n}\right)\frac{t}{n\Delta t_{\mathrm{c}}}\right].$$

$${\mathrm{P}}_{\mathrm{D}}\left(\frac{jt}{n}\right)\approx {\mathrm{P}}_{\mathrm{D}0}+\frac{jt}{n}\frac{\mathrm{d}{\mathrm{P}}_{\mathrm{D}0}}{\mathrm{d}t},$$

$$h(t)=\sum _{j=0}^n\left[\frac{{\mathrm{P}}_{\mathrm{D}0}}{n\Delta t_{\mathrm{c}}}+\frac{2jt}{n^2\Delta t_{\mathrm{c}}}\frac{\mathrm{d}{\mathrm{P}}_{\mathrm{D}0}}{\mathrm{d}t}\right]{\left[1-{\mathrm{P}}_{\mathrm{D}}\left(\frac{jt}{n}\right)\frac{t}{n\Delta t_{\mathrm{c}}}\right]}^{-1}.$$

$$h(t)\le \frac{1}{\Delta t_{\mathrm{c}}}\frac{n+1}{n}\left({\mathrm{P}}_{\mathrm{D}0}+t\frac{\mathrm{d}{\mathrm{P}}_{\mathrm{D}0}}{\mathrm{d}t}\right){\left(1-\frac{t}{n\Delta t_{\mathrm{c}}}\right)}^{-1}.$$

$$h(t)\le \frac{n+1}{n}{\left(\Delta t_{\mathrm{c}}-\frac{t}{n}\right)}^{-1}.$$

$$h(t)\lesssim \frac{1}{\Delta t_{\mathrm{c}}-\Delta t_{\mathrm{s}}}.$$

Full text

Keywords: tokamak plasma, disruption prediction, survival analysis, machine learning, Random Forest, binary classification

1 Introduction

Plasma disruptions in tokamaks pose a serious risk to current experiments and future fusion devices. During a disruption, the total thermal and magnetic energies—upwards of tens or hundreds of megajoules in future devices like SPARC [1] or ITER [2]—can be dissipated in tens of milliseconds or fewer, leading to (i) high heat fluxes on plasma-facing components, (ii) large induced eddy and halo currents in the surrounding vacuum vessel, and (iii) generation of highly-relativistic “runaway” electrons. Ideally, plasma disruptions should be avoided altogether, thereby preventing damage and continuing operation. Although in reality, many disruptions will need to be mitigated; that is, the plasma will be terminated in such a way—e.g. through massive gas or shattered pellet injection—so as to minimize both damage and delay in operation.

However, for both avoidance and mitigation, an impending disruption must first be predicted with enough time to actuate an appropriate response. There are many past and ongoing efforts to develop disruption prediction algorithms, including work on tokamaks ADITYA [3, 4], Alcator C-Mod [5, 6], ASDEX-U [7, 8, 9], DIII-D [10, 11, 5, 12, 6], EAST [6], JET [13, 8, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23], JT60-U [24, 25], J-TEXT [26, 27], and NSTX [28], among others. At present, most predictors estimate the plasma state at each moment in time, usually identified as either non-disruptive or disruptive. While this is an important and necessary step toward accurate disruption prediction, it is not sufficient to know, even with 100% accuracy, that a disruption is currently occurring. Instead, the ultimate predictive capability is to estimate the future probability of a disruption. The goal is to answer the question, what is the probability that the plasma will survive a time $\Delta t$ into the future? From this prediction, combined with knowledge of plasma actuators and mitigation systems, an attempt can be made to navigate the plasma state away from disruptive territory, or mitigation would be employed.

A framework for predicting events and estimating survival times already exists in the area of statistics called survival analysis. Decades of research in such fields as medicine, engineering, and sociology have utilized and refined these tools and techniques. The fusion community has an opportunity to leverage this established knowledge base to confront the challenge of disruptions. The present work is not the first application of survival analysis for event prediction in the realm of fusion. Only recently, “direct hazard modeling” was used to study and predict the onset of neoclassical tearing modes [29, 30] which can often lead to disruptions of tokamak plasmas. The authors of [29] actually note the potential of the survival analysis approach for disruption prediction.

The aim of this work is to demonstrate the integration of existing disruption prediction techniques into the survival analysis framework. This work focuses on the Disruption Predictor using Random Forests (DPRF) [6] algorithm with real-time prediction capabilities, trained and tested on binary-classified data from the Alcator C-Mod tokamak [31], a high magnetic field ($B_{0}$ = 2-8 T), compact ($R_{0}$ = 68 cm, $a$ = 22 cm) device with plasma currents and densities of order 1 MA and $10^{20}$ m*-3*, respectively. A database of disruption parameters has been compiled for over four thousand C-Mod discharges, containing over $10^{5}$ data points in total [5, 6].

This integration of an existing algorithm within a new framework is important for two reasons: First, as mentioned, most disruption predictors only consider a binary-classification of the current plasma state; thus, the relationship between a class label and the time until the disruption must be carefully considered when making time-to-event predictions with survival analysis. Second, these existing disruption prediction algorithms are the current standard against which survival analysis must be compared; this comparison is performed, in part, in the present work. The authors note here (and in section 7) that there is certainly more work to be done beyond the present analysis; the ultimate goal is a framework incorporating the evolution models of many plasma and operational parameters in a time-to-disruption prediction. However, this is outside the scope of the current paper, and optimization of these methods is left for future work.

The outline of the rest of the paper is as follows: Section 2 discusses binary classification of disruption databases and RF algorithms trained on those data sets. In section 3, survival analysis is introduced and explored as an approach to disruption prediction. This methodology is then applied to both test and real experimental data in sections 4 and 5, respectively. In section 6, a summary of results is given. Finally, limitations of this analysis and opportunities for future work are presented in section 7.

2 Disruption database classification and prediction

This section discusses the binary classification scheme commonly used to classify data in disruption databases. The framework described here, in its current form, is only applicable to predictors using binary classification. Future work could expand these methods to multi-class data. We focus on the supervised Random Forest (RF) [32] machine learning approach to disruption prediction. Most importantly, we discuss the implications of binary classification and interpretations of RF algorithm output signals. Regarding notation in this and following sections, we use $t$ for time, $\Delta t$ for time intervals or windows, and $\tau$ for times of interest in survival analysis.

2.1 Implication of binary classification

Many efforts to develop disruption prediction algorithms follow a similar methodology, summarized here: First, a set of plasma parameters—e.g. Greenwald fraction and internal inductance—and operational parameters—e.g. differences between programmed and measured values of plasma current or vertical position—are chosen based on their relevance to disruption physics as well as their real-time measurement capability. Second, databases of these parameters are assembled for many machines and many times throughout plasma discharges, some terminated by disruptions; for examples, see [10, 7, 13, 17, 28, 11].

These data are then typically split into two classes: disruptive (D) and non-disruptive (ND). All times from non-disrupting plasma discharges are considered ND, but only a subset of times from plasma discharges ending in disruptions fall in class D. Here, it is assumed that times far from the disruption reside in a ND region of parameter space; thus, only times close to the disruption are classified as D. The time $\Delta t_{\rm c}$ dividing the two classes can be physically-motivated or optimized for best performance of the disruption prediction algorithm. In addition, due to the finite response times of the plasma control system (PCS), actuators, and mitigation system, data within $t_{\mathrm{D}}-t<\Delta t_{\rm min}$ of the disruption time $t_{\mathrm{D}}$ are often excluded from analyses, where $\Delta t_{\rm min}$ is the minimum time required to avoid or mitigate a disruption. See figure 1 for an illustration of the timeline of a disruptive discharge.

Note that in much of the previous work on this topic, e.g. [14, 28, 21, 5, 6], the time evolution of parameters is not taken into account in the algorithm development; that is, each time-slice of each class is treated as independent of all other times. While time-independent studies have set up a good framework for disruption prediction, a temporal analysis is absolutely necessary for future prediction algorithms. Ultimately, it would be extremely beneficial for the PCS to know not only the current plasma class (here D or ND), but also the probability of a disruption in the next time $\Delta t$, as discussed below, in order to avoid or mitigate the disruption.

Consider a database populated with plasma and operational parameters from many plasma discharges, which has been split into classes ND and D by a threshold time $\Delta t_{\rm c}$ before the disruption time $t_{\mathrm{D}}$, excluding data $t_{\mathrm{D}}-t<\Delta t_{\rm min}$. If a data point is randomly selected from class D, the only information known regarding the time until the impending disruption is that $\Delta t_{\rm min}\leq t_{\mathrm{D}}-t\leq\Delta t_{\rm c}$, where $t_{\mathrm{D}}-t$ is the unknown time interval between the time $t$ at which the data point was collected and the disruption time.111Of course, both $t$ and $t_{\mathrm{D}}$ can be determined from the database, but time dependence is not considered in the development of these RF algorithms. In other words, we know that each data point in class D comes from a disruptive discharge; that is how we assembled the database, classifying data based on $\Delta t_{\rm c}$. However, because our algorithm neglects the time dependence of data within class D, it is ignorant of the time at which each data point occurred, only knowing that $t_{\mathrm{D}}-\Delta t_{\rm c}\leq t\leq t_{\mathrm{D}}-\Delta t_{\rm min}$ (refer to figure 1). The probability that the random variable $t$ lies within a time interval $\Delta t>0$ before the disruption time is given by

[TABLE]

Here, the probability is conditional on the data point being in class D, and a uniform probability density is chosen, as a conservative assumption, due to our ignorance of the actual time until the disruption. Other probability density functions could be explored, as discussed in section 7, but since no one option seems particularly advantageous or physically-motivated, uniform is adopted for simplicity. It is important to note that this formulation is not ad hoc; rather, it is a consequence of the binary classification scheme. Furthermore, equation (2.1) is just another way of writing the probability that a disruption occurs in the interval $\Delta t$ following any given data point in class D, i.e. $\mathrm{P}(t+\Delta t>t_{\mathrm{D}}|\mathrm{D})$. This is the probability of interest for real-time disruption prediction, where $t$ is known but $t_{\mathrm{D}}$ is not.

From (2.1), we see that our choice of $\Delta t_{\rm c}$ and the value of $\Delta t_{\rm min}$ (typically machine-dependent) affect the interpretation of data within each class. Nonetheless, this formulation provides a physical interpretation of data within class D beyond the label “disruptive.” That is, given any data point in class D, the probability of a disruption occurring in the next time interval $\Delta t$ increases linearly from 0 to 1 over $\Delta t_{\rm min}\leq\Delta t\leq\Delta t_{\rm c}$. For simplicity in the upcoming analyses, $\Delta t_{\rm min}$ has been set to zero; this is an allowable assumption since $\Delta t_{\rm min}$ is oftentimes much less than the chosen or optimized $\Delta t_{\rm c}$.

It follows from the product rule that, given a random data point from either class ND or D, the joint probability that the data point is in class D and the plasma will disrupt in the next time $\Delta t$ is given as

[TABLE]

where P(D) is the probability of the data point belonging to class D. The calculation of P(D) is discussed in the following section.

2.2 Interpretation of Random Forest predictions

Given a disruption database split into two classes (ND and D, as described in the previous section) and an array of input parameters corresponding to one plasma state, most disruption predictors have been designed to evaluate the state’s membership in either class or the distance to the “boundary” between classes. There are many approaches to create such algorithms. Perhaps the simplest is a scheme in which thresholds (of physical quantities in parameter space) between classes ND and D are determined by minimizing overlap in histograms of data from each class, as discussed in [28, 11, 5, 23] among other works. Machine learning methods have also been implemented, trained and tested on these data sets (see the references listed in section 1). Generally, algorithms are optimized to maximize correct classifications of states in class D—i.e. to “catch” disruptions (true positives)—and minimize incorrect classifications of states in class ND—i.e. to reduce false positives.

This work focuses on the RF machine learning approach to disruption prediction, which has been applied on several tokamaks; see, for example, [15, 16, 11]. In the RF framework, the parameter space of the disruption database is randomly and iteratively subdivided until the final nodes (i.e. “leaves” at the ends of “branches” of one decision “tree”) contain instances of only one class from the training set. Many trees comprise the total “forest.” As new data is fed through all trees, each tree results in an evaluation of 0 (ND) or 1 (D). Thus, an average over the entire forest results in an effective probability that plasma state is in class D.

This output signal, called the disruptivity in this work, is exactly the quantity of interest P(D) in (2). At each time step, the disruptivity is calculated based on the current vector of plasma and operational parameters; therefore, the time-evolving disruptivity signal is denoted $\mathrm{P}_{\mathrm{D}}(t)$. Real-time calculations of this signal have been demonstrated by the DPRF algorithm implemented in the PCS system of the DIII-D tokamak [12]. Efforts are underway to perform full-discharge analyses of disruptivity signals, incorporating time evolution and optimizing both disruptivity thresholds and time-windows to trigger warnings [6]. These are discussed further in section 5.

3 Survival analysis

The aim of survival analysis is to estimate the time until a specified event—oftentimes failure or death—occurs. For instance, in medicine, survival analysis techniques are used to study patients subject to different medical treatments, importantly assessing their expected lifetime (i.e. time until death). In such a study, some data will be censored222Specifically, this data would be right-censored., whereby the event (death) will not be observed for all patients during the time of the study. This is similar, in some ways, to disruption studies as researchers often only consider data from part of the plasma discharge, such as during the flattop portion of the plasma current, before the discharge ends. However, unlike the individuals in the medical study presented above, plasmas can (effectively) be sustained indefinitely, as could be the case in a future fusion power plant. Thus, the event (here, a disruption) might never occur at all for some plasmas, no matter their discharge length.

There are many approaches within the field of survival analysis which could be applied to disruption prediction. Some attempt to determine the parametric dependence of future survival. This would be, of course, very useful from the perspective of plasma control and disruption avoidance, especially if event predictions are early enough for plasma actuation. However, disruption mitigation is perhaps the more pressing issue since, as of yet, many disruptions cannot be predicted with adequate warning time for mitigation. Therefore, in this paper, we focus on ascertaining plasma survival probabilities and times. The survival function $S(t)$ computes the probability that the survival time $T$ surpasses a given time $t$,

[TABLE]

The survival function is defined such that $S(t_{0})\equiv 1$ at the starting time $t_{0}$ and $S(t)$ decreases monotonically to zero as $t\to\infty$. The corresponding hazard function $h(t)$ gives the instantaneous rate of failure, assuming survival until time $t$,

[TABLE]

If we know the hazard function over a time interval $t_{a}\leq t\leq t_{b}$, the survival function can be computed from (4),

[TABLE]

There are many different approaches to compute $S(t)$ for time-evolving data sets. One of the earliest and most straightforward methods is the Kaplan-Meier formalism [33], which is adopted in this work. Given a discretized timebase, the probability of survival beyond time $t_{n}$ is given by the “product-limit” formula,

[TABLE]

That is, the probability that the survival time $T$ exceeds time $t_{n}$ is the product of the probabilities of $T$ exceeding each time step $t_{i}$ before $t_{n}$ given that $T$ actually reaches each $t_{i}$. If the probability of failure between consecutive times, i.e. $\mathrm{P}_{i\rightarrow i+1}$, is known, the survival function can then be written as

[TABLE]

This is essentially equation (7) in [29].

Applying the Kaplan-Meier formalism to disruption prediction is relatively straightforward: The event or failure is simply the disruption itself. For any time $t_{i}$, the probability of disruption in the incremental step between $t_{i}$ and $t_{i+1}$ is given by (2), where $\Delta t=t_{i+1}-t_{i}$ and $\mathrm{P}(\mathrm{D})=\mathrm{P}_{\mathrm{D}}(t_{i})$ is the evaluation of the disruption prediction algorithm at time $t_{i}$. Therefore, if the current time is $t_{i}$, the survival probability beyond future time $t_{n}$ is calculated

[TABLE]

This is written as a conditional probability because it implicitly assumes that the plasma has survived until time $t_{i}$. Because such a predictor would ultimately be coupled with the PCS, we assume that the disruptivity is calculated with a sampling time smaller than the class time, i.e. $\Delta t_{\rm s}<\Delta t_{\rm c}$. Then, the probability of survival $\Delta t$ into the future from current time $t$ is

[TABLE]

Here, the future time $\Delta t$ is split into $n$ sampling times (i.e. $\Delta t=n\Delta t_{\rm s}$) with the constraint $\Delta t_{\rm s}\leq\Delta t_{\rm c}$. The disruptivity $\mathrm{P}_{\mathrm{D}}$ is evaluated at each time step $j\Delta t_{\rm s}$ from the current time $t$, and the final term $\Delta t_{\rm s}/\Delta t_{\rm c}$ follows from the assumption of uniform probability density in (2.1). In reality, the future disruptivity signal is unknown, so a simple linear extrapolation is used in this work,

[TABLE]

where $\mathrm{P}_{\mathrm{D}}(t)$ and $\mathrm{d}\mathrm{P}_{\mathrm{D}}/\mathrm{d}t$ are obtained from the DPRF output. For the purposes of this paper, $\mathrm{d}\mathrm{P}_{\mathrm{D}}/\mathrm{d}t$ is calculated from a linear least-squares fit over a predetermined time window of length $\Delta t_{\rm fit}$ before the current time, and the extrapolated disruptivity is restricted to the interval $\mathrm{P}_{\mathrm{D}}\in\left[0,1\right]$. While this is certainly not an optimized form of extrapolation, it at least provides a starting point from which to evaluate the survival analysis approach. Possible alternatives will be discussed in section 7.

Equation (9) is powerful in that its evaluation provides an actual estimate for the probability of plasma survival (i.e. no disruption) a time $\Delta t$ into the future, and it requires only a RF algorithm (well-)trained on binary-classified disruption data. Still, this approach does not actually trigger an alarm itself; we need to determine the threshold for mitigation or specific avoidance techniques. One possible alarm could be a threshold for the median time $\tau_{50}$, defined by

[TABLE]

This is the time into the future beyond which the probability of survival is less than 50%. Another metric could be the expected future lifetime $\tau$, calculated as

[TABLE]

It is worth noting that the expected future lifetime is usually (approximately) greater than or equal to the class time, i.e. $\tau\gtrsim\Delta t_{\rm c}$, for cases when $\Delta t_{\rm s}\ll\Delta t_{\rm c}$; this is due to the ignorance introduced from binary classification and a fast sampling rate which frequently updates the state of the plasma. (See A for a derivation of this inequality and further discussion.) For either case, if $\tau_{50}$ or $\tau$ would reach some threshold time, say $\Delta t_{\rm min}$ or $\Delta t_{\rm c}$, then the disruption mitigation system could be triggered.333Note that for some specific cases, like $\mathrm{P}_{\mathrm{D}}(t)=0$, $\tau_{50}$ does not exist and $\tau$ diverges. In this paper, we assume that $\tau_{50}$ exists and $\tau$ is finite. The median time and lifetime are used as examples throughout the rest of this paper, but choices of the most reliable metrics and optimization of thresholds are left for future work.

Here, it is prudent to distinguish the present application of survival analysis from that in [29, 30]. In [29, 30], a direct hazard model was used to calculate multivariable hazard functions from machine learning of experimental data, which could then be related to survival probabilities via (5). This is, in a way, opposite to the approach adopted in this paper; nevertheless, their model has benefits including the incorporation of time-dependent covariates. It is unclear, though, if these hazard functions can be calculated in real-time, which would be necessary for integration with a PCS.

4 Application to test data

In this section, the Kaplan-Meier formalism for survival analysis is applied to test data in two cases: First, smooth test data is used to demonstrate the capabilities and limitations of this approach, including calculations of the survival and hazard functions. Second, an application to noisy test data highlights the effects of linear extrapolation, noise, and fitting windows. Later, in section 5, the survival analysis approach is applied to experimental data from the Alcator C-Mod tokamak to show its potential use in “real time.”

4.1 Smooth data

In this test case, a simple, ad hoc disruptivity signal was created and is shown in figure 2a. As is seen, there is an initial rapid increase in $\mathrm{P}_{\mathrm{D}}$, peaking at $t$ = 0.2 s and then decaying to a steady-state value of $\mathrm{P}_{\mathrm{D}}=0.3$ from $t$ = 0.3-0.6 s. At $t$ = 0.6 s, the disruptivity signal increases gradually until reaching $\mathrm{P}_{\mathrm{D}}=1$ at $t$ = 1 s. Note that this does not necessarily imply that a disruption would occur at $t$ = 1 s in this test case; instead, a disruption would be likely to occur within the next class time, i.e. $t\in[1,1+\Delta t_{\rm c}]$ s, if the disruptivity remained at $\mathrm{P}_{\mathrm{D}}=1$.

Contours of survival probability, calculated from (9), are shown in figures 2b-c, with current time $t$ and time into the future $\Delta t$ plotted on the horizontal and vertical axes, respectively. Here, a class time of $\Delta t_{\rm c}$ = 100 ms and sampling time of $\Delta t_{\rm s}$ = 1 ms were assumed. Each color corresponds to one decade within the range 0-100%. In figure 2b, the future disruptivity signal is known for all times; that is, this calculation predicts the future. Thus, data beyond $t+\Delta t>$ 1 s is unknown (grey region). Overlaid are contours of $\tau_{50}$ (solid) and the effective future lifetime (dotted), calculated as

[TABLE]

since data are only available up to $t$ = 1 s.

It is important to note several features of figure 2b that may seem counter-intuitive. First, even though $\mathrm{P}_{\mathrm{D}}(t)$ is known for all $t$ (and $\Delta t$), the contours vary in time. This is because the calculation “domain” changes at each $t$, and the probability of survival at the current time ($\Delta t=0$) is always 1. Also, one may have hoped that the effective future lifetime would always satisfy $t+\tau_{\rm eff}$ = 1 s since the calculation “knows” that $\mathrm{P}_{\mathrm{D}}$($t$ = 1 s) = 1. However, recall that a disruptivity value of 1 only conveys that a disruption is likely to occur within the next time $\Delta t_{\rm c}$. In fact, the value of survival probability at $t+\Delta t$ = 1 s is not always zero. What is more, the data of figure 2b represents the best possible predictive capability of survival analysis applied to this particular test case. This is the “ideal” scenario against which we should compare more realistic scenarios.

Figure 2c shows contours of the survival function calculated using a linear extrapolation of $\mathrm{P}_{\mathrm{D}}(t)$, as in (10), with fitting window $\Delta t_{\rm fit}$ = 100 ms. Also overlaid are contours of $\tau_{50}$ and the approximate future lifetime $\widetilde{\tau}$, calculated as

[TABLE]

Here, predictions are only made $\Delta t$ = 1 s into the future, a limitation which could be varied. This example highlights some drawbacks of the linear extrapolation method. The initial rise in disruptivity extrapolates to $\mathrm{P}_{\mathrm{D}}$ = 1 around $t\approx$ 0.3 s, causing both $\tau_{50}$ and $\widetilde{\tau}$ to dip close to $\Delta t_{\rm c}$ around $t\approx$ 0.2 s. Then, the decrease in disruptivity from $t\approx$ 0.2-0.3 s extrapolates to $\mathrm{P}_{\mathrm{D}}$ = 0 around $t\approx$ 0.4 s, at which time the plasma is predicted to survive far ($\Delta t\geq$ 1 s) into the future.

The absolute difference between survival functions in figures 2b and 2c is given in figure 2d, with magnitude in the range $[-1,1]$. While there is significant disagreement between the linearly-extrapolated and “ideal” cases around $t\approx$ 0.3 s, most other times agree within $\pm$270 ms and some within $\pm$90 ms. Future work should explore how to best minimize these differences; some suggestions are made in section 7.

Finally, contours of the hazard function are plotted in figure 2e, calculated using the data in figure 2c and (4). As expected, the hazard is highest from $t$ = 0-0.2 s and 0.7-1.0 s and lowest from $t$ = 0.2-0.4 s. As derived in B, the hazard function can take values approximately in the range $h\in[0,(\Delta t_{\rm c}-\Delta t_{\rm s})^{-1}]$ s*-1*; these failure rates correspond to predicted future lifetimes ranging from approximately $\Delta t_{\rm c}$ to $\infty$ when $\Delta t_{\rm s}\ll\Delta t_{\rm c}$. Because these ranges are dependent on the chosen class time, we must rely on comparing relative values. Future work should explore “calibrating” $h(t)$ and/or finding an appropriate combination of thresholds for both $S(t)$ and $h(t)$.

4.2 Noisy data

The disruptivity signal used in the above analysis is smooth, unlike real data which will be noisy. Artificial noise, with random amplitude in the range $[-0.1,0.1]$ and with period (peak-to-peak) of 50 ms, was added to the smooth signal and is plotted in figure 3a. These variations in $\mathrm{P}_{\mathrm{D}}$ are similar to those observed in the disruptivity calculated for real C-Mod data, as seen in figures 4-6. The survival probability, calculated assuming the future disruptivity is known, is shown in figure 3b; for low amplitude noise, it is (unsurprisingly) quite similar to that in figure 2b.

Figures 3c-e show contours of the survival function computed using linear extrapolation, from (10), with fitting windows of $\Delta t_{\rm fit}$ = 50, 100, and 200 ms, respectively. The median times and approximate future lifetimes, from (14), are also overlaid. Note, as seen in figure 3c, how using a time window on the order of the period of noise oscillations greatly affects $S(t)$ , with $\tau_{50}$ and $\widetilde{\tau}$ varying widely on this timescale. As the time window increases, the contours smooth out in time $t$, and variations in $\tau_{50}$ and $\widetilde{\tau}$ decrease. However, the “responsiveness” to changes on the order of $\Delta t_{\rm fit}$ is also delayed; this is observed as both a rightward shift and broadening of the light-colored feature ($\mathrm{P}_{\mathrm{D}}\geq 0.7$) in figures 3c-e as $\Delta t_{\rm fit}$ increases from 50 to 200 ms. Future work should consider optimization of the fitting time window as well as proper uncertainty quantification. Here, it appears that a value of $\Delta t_{\rm fit}$ between 100-200 ms, i.e. 2-4 times the noise period, would best approximate the smooth data in figure 2c.

5 Application to Alcator C-Mod data

As mentioned in section 1, the Random Forest algorithm DPRF has been trained and tested on C-Mod data [5, 6], with most recent results optimizing an alarm window to trigger mitigation, as reported in [6]: If the disruptivity increases above a high threshold without decreasing below a low threshold over a time window444Note that this alarm time window $\Delta t_{\rm W}$ is different from the fitting time window $\Delta t_{\rm fit}$ discussed elsewhere in the present paper. $\Delta t_{\rm W}$, an alarm is triggered. In this work, the optimized DPRF algorithm had high and low threshold values of $\mathrm{P}_{\mathrm{D}}$ = 0.35 and 0.05, respectively, an alarm window $\Delta t_{\rm W}$ = 5 ms, and a class time $\Delta t_{\rm c}$ = 325 ms. See [6] for more details on the DPRF algorithm optimization.

This section investigates three C-Mod discharges on which this threshold-based alarm system has been tested in “real time.” The first, in section 5.1, is a good prediction, or true positive, meaning that an alarm would have been triggered with sufficient time to mitigate (or avoid) the impending disruption. The second, in section 5.2, is an example of a late warning, or false negative, for which an alarm is triggered with too little time to mitigate. Finally, section 5.3 presents a case of a false positive; i.e. an alarm would have been triggered for a non-disrupting discharge. An example of a true negative—when a non-disrupting discharge is correctly identified—is not included in the present work because disruptivity values $\mathrm{P}_{\mathrm{D}}\approx 0$ do not illuminate any additional capabilities of the survival analysis approach.

5.1 A good prediction

Figures 4a-b show plasma parameters for C-Mod discharge #1140226013, during which the density increased steadily and toroidal magnetic field (not shown) decreased in time from $B_{0}$ = 5.3 T to 3.9 T over $t\approx$ 0-0.7 s. Since the plasma current was held constant at $I_{\rm p}$ = 0.8 MA, the edge safety factor decreased to $q_{95}\approx$ 3. A locked mode, identified by a reduction in sawteeth observed in the plasma temperature (not shown), began around $t\approx$ 1 s and was likely the cause of the disruption at $t\approx$ 1.4 s. The disruptivity signal in figure 4a seems to accurately predict the disruption: $\mathrm{P}_{\mathrm{D}}$ remains low during the first part of the flattop current ($t\approx$ 0.5-1.0 s) with $\mathrm{P}_{\mathrm{D}}$ increasing after $t\approx$ 1 s and reaching $\sim$1 approximately 200 ms before the disruption. The optimized DPRF predictor, discussed above, performs well on this discharge. The calculated disruptivity $\mathrm{P}_{\mathrm{D}}$ crosses the high threshold $\mathrm{P}_{\mathrm{D}}=0.35$ and remains above the low threshold $\mathrm{P}_{\mathrm{D}}=0.05$ (shown as solid horizontal lines in figure 4a) over the required time window $\Delta t_{\rm W}$ = 5 ms before the class time $\Delta t_{\rm c}$ = 325 ms, indicated by the vertical dashed line .

The application of survival analysis to this experimental data gives another perspective of the disruptivity’s predictive capability. Figures 4c-d show the survival and hazard functions, using linear extrapolation from (10) and a fitting time window $\Delta t_{\rm fit}$ = 100 ms. Note that the disruptivity was actually calculated with two sampling times distinguished by proximity to the disruption time $t_{\mathrm{D}}$: $\Delta t_{\rm s}$ = 20 ms when $t_{\mathrm{D}}-t>$ 20 ms, and $\Delta t_{\rm s}$ = 1 ms when $t_{\mathrm{D}}-t\leq$ 20 ms.666These sampling rates were chosen to reduce the total database size and were based on considerations unrelated to the present study. The database can be updated to vary/increase the sampling rate. Here, it is assumed that a real-time sampling rate of $\Delta t_{\rm s}$ = 1 ms is achievable, so $\mathrm{P}_{\mathrm{D}}(t)$ was interpolated appropriately. The median time and approximate future lifetime overlay the survival probability in figure 4c. Note that $\tau_{50}$ falls below $\Delta t=\Delta t_{\rm c}$, plotted also as a horizontal dashed line in figure 4c, and $\widetilde{\tau}$ approaches $\Delta t_{\rm c}$ around the time that $t_{\mathrm{D}}-t=\Delta t_{\rm c}$. In fact, $\widetilde{\tau}\approx\Delta t_{\rm c}$ during almost all times $t_{\mathrm{D}}-t<\Delta t_{\rm c}$. This is due to the linear extrapolation of $\mathrm{P}_{\mathrm{D}}$ to 1 around this time. Thus, we conclude that the survival function has “accurately” predicted the disruption in this discharge. In this scenario, we could imagine setting an alarm with a requirement that $\tau_{50}$ remains below $\Delta t_{\rm c}$ (or $\widetilde{\tau}\approx\Delta t_{\rm c}$) for a certain time interval, perhaps a few confinement times ($\tau_{\rm E}\sim$ 20-30 ms in C-Mod). The hazard function is shown for completeness in figure 4d. As expected, $h(t)$ increases in amplitude as $\tau_{50}$ and $\widetilde{\tau}$ decrease, reaching a maximum at $h=(\Delta t_{\rm c}-\Delta t_{\rm s})^{-1}\approx$ 3.1 s*-1*.

5.2 A late warning

Plasma parameters for C-Mod discharge #1150722006, which ended in a disruption at $t\approx$ 1.15 s, are shown in figures 5a-b. The disruptivity signal $\mathrm{P}_{\mathrm{D}}$ output from the DPRF algorithm remained low during much of the flattop current, i.e. $t\approx$ 0.2-0.7 s. During this time, the plasma $\beta$ was steadily increasing, while the edge safety factor (not shown) decreased. As the plasma density decreased from $t\approx$ 0.7-0.9 s, $\mathrm{P}_{\mathrm{D}}$ increased slightly, although not enough to pass the high threshold $\mathrm{P}_{\mathrm{D}}=0.35$ and trigger an alarm. Finally, only $\sim$20 ms before the disruption, $\mathrm{P}_{\mathrm{D}}$ rose rapidly above the threshold; this likely resulted from an impurity injection, which was observed as a fast increase in radiated power (not shown). Unfortunately, the disruptivity passed the threshold $\mathrm{P}_{\mathrm{D}}>0.35$ too close to the disruption, qualifying this as a false negative; that is, the DPRF failed to predict the disruption with enough time to mitigate or avoid it.

For this discharge, the survival and hazard functions, calculated using a fitting window $\Delta t_{\rm fit}$ = 100 ms and sampling time $\Delta t_{\rm s}$ = 1 ms, are plotted in figure 5c-d. As is seen, neither the median time $\tau_{50}$ nor approximate future lifetime $\widetilde{\tau}$ fall below the class time $\Delta t=\Delta t_{\rm c}$ at any point during the discharge. Thus, the predictive capabilities of survival analysis would fail to adequately warn of an impending disruption in this case. In fact, unlike a disruptivity-threshold scheme, the survival analysis framework may have not triggered an alarm at all since both $\tau_{50}$ and $\widetilde{\tau}$ are greater than $\Delta t_{\rm c}$, even near the disruption time. More work must be done to optimize this framework.

5.3 A false alarm

Figures 6a-b show plasma parameters for C-Mod discharge #1140227018. During this discharge, tokamak operators applied error fields to purposefully achieve a locked mode. However, the plasma did not lock and was successfully ramped down without disrupting, as seen in the time traces of plasma current and density. The DPRF algorithm generated a disruptivity signal which has a large spike during the current ramp-up and later a steady increase from $\mathrm{P}_{\mathrm{D}}\approx$ 0-0.5 over $t\approx$ 0.75-1.5 s. Considering only the flattop portion of the discharge, the rise in $\mathrm{P}_{\mathrm{D}}$ was likely due to the peaking of the current profile, seen as an increase in internal inductance, as well as magnetic pickup from the error field coils. Because $\mathrm{P}_{\mathrm{D}}$ crossed the high disruptivity threshold of $\mathrm{P}_{\mathrm{D}}$ = 0.35 at $t\approx$ 1.2 s, this discharge qualifies as a false positive; that is, in the proposed alarm-threshold scheme, this discharge would be identified as disruptive even though it does not actually disrupt.

However, using the estimated future lifetime threshold suggested above, the survival analysis framework correctly classifies this discharge. The survival and hazard functions, computed using the linear extrapolation method with fitting window $\Delta t_{\rm fit}$ = 100 ms and sampling time $\Delta t_{\rm s}$ = 1 ms, are shown in figures 6c-d. The data are quite noisy, with large variations in median time $\tau_{50}$ and approximate future lifetime $\widetilde{\tau}$. A larger fitting time window, as discussed in section 4.2, could improve this. Note that $\tau_{50}$ falls below the class time $\Delta t_{\rm c}$ at $t\approx$ 1.2 and 1.35 s, the same times as the peaks in $\mathrm{P}_{\mathrm{D}}$ seen in figure 6a. However, the estimated future lifetime $\widetilde{\tau}$ never decreases to $\Delta t_{\rm c}$ during the flattop portion of the discharge. Of course, a wide variety of thresholds for $\tau_{50}$ and $\widetilde{\tau}$ could be chosen such that a false alarm would have been triggered by this survival analysis approach. Nevertheless, physics intuition might tell us that if the expected future lifetime exceeds the class time chosen for disruptive data, then an alarm should not be sounded.

6 Summary

In this paper, we proposed a survival analysis approach to disruption prediction, applying the Kaplan-Meier formalism to existing Random Forest (RF) algorithms trained on binary-classified data. While it is unknown whether this method could predict disruptions with higher accuracy or more warning time than traditional machine learning (ML) methods, we argued that survival analysis provides a more physically-meaningful and more easily-interpretable framework for disruption prediction.

This work required a database of binary-classified data, split into non-disruptive and disruptive classes by a class time $\Delta t_{\rm c}$ before the disruption. Given that a plasma state was disruptive, we assumed that the time of the impending disruption was uniformly distributed. In addition, the output “disruptivity” signal $\mathrm{P}_{\mathrm{D}}(t)$ from a Random Forest (RF) algorithm was interpreted as the probability that the plasma state at time $t$ belongs to the disruptive class. Therefore, at time $t_{i}$, knowledge of $\mathrm{P}_{\mathrm{D}}(t_{i})$ and $\Delta t_{\rm c}$ allows evaluation of the probability of a disruption in a future time interval $\Delta t$.

The Kaplan-Meier formalism, from survival analysis, was then used to calculate the survival probability $S(t)$ from a time-evolving disruptivity signal. At any time $t$, the survival function $S(t+\Delta t|t)$ actually gives the probability of survival (i.e. no disruption) a time $\Delta t$ into the future. From this, the median time $\tau_{50}$ and (approximate) future lifetime $\widetilde{\tau}$ were used as indicators of plasma “health,” where low values of either $\tau_{50}$ or $\widetilde{\tau}$ could be used to trigger an alarm of an impending disruption.

Application of this model to test data illuminated some non-intuitive, probabilistic features of the survival function—including $\tau_{50},\widetilde{\tau}>0$ for $\mathrm{P}_{\mathrm{D}}(t)=1$—and drawbacks of linearly-extrapolating $\mathrm{P}_{\mathrm{D}}(t)$ to predict future disruptivity values. Moreover, the hazard function $h(t)$, which indicates the instantaneous rate-of-failure, was found to be consistent with expectations: low values of $h(t)$ for high values of $\tau_{50}$ and $\widetilde{\tau}$, and vice versa. However, $h(t)$ is dependent on the choice of $\Delta t_{\rm c}$, so comparison of relative values may only be useful. When noise was added to the smooth test disruptivity signal, $S(t)$ proved to be quite sensitive to the length of the “fitting” time window $\Delta t_{\rm fit}$ used for linear extrapolation. If such a method is used in future analyses, this time window should be several ($\sim$2-4) times longer than the longest period of noise oscillations.

This survival analysis methodology was also applied to real disruption data from three Alcator C-Mod discharges. The disruptivity signal was generated from the DPRF algorithm trained on the C-Mod disruption database and optimized with $\Delta t_{\rm c}$ = 325 ms. In one discharge, the median time $\tau_{50}$ fell below $\Delta t_{\rm c}$ and the expected future lifetime $\widetilde{\tau}$ approached $\Delta t_{\rm c}$ for times within $\Delta t_{\rm c}$ of the disruption time $t_{\mathrm{D}}$; that is, $\tau_{50},\widetilde{\tau}\lesssim\Delta t_{\rm c}$ for $t_{\mathrm{D}}-t\leq\Delta t_{\rm c}$, as hoped. For this specific scenario, the application of survival analysis to disruption prediction should be considered successful. However, in another case, both the disruptivity threshold scheme and survival analysis approach were not able to adequately warn, with sufficient time, of an impending radiative collapse caused by an impurity injection. For the final discharge, the optimized DPRF disruptivity signal $\mathrm{P}_{\mathrm{D}}$ and thresholds falsely triggered an alarm for a non-disrupting plasma; the approximate future lifetime $\widetilde{\tau}$ from survival analysis, however, remained above the class time $\Delta t_{\rm c}$, indicating that such a false alarm could have been avoided.

7 Opportunities for future work

The survival analysis approach presented in this paper has several limitations, and there is much work to be done before this framework could be implemented in real time in a plasma control system. Future work must explore a variety of metrics (i.e. beyond $\tau_{50}$ and $\widetilde{\tau}$) and optimize thresholds and durations required to trigger mitigation or pursue avoidance strategies. For instance, perhaps a predictor would be more successful by monitoring the third-quartile time, i.e. $S(t+t_{75}|t)=0.75$, and triggering an alarm if $t_{75}$ decreases below the minimum time $\Delta t_{\rm min}$ required for mitigation .

Here, the authors would like to comment on some specific limitations of this analysis and suggest several opportunities for improvement:

In this paper, we assumed the existence of an RF disruption prediction algorithm already trained on disruption-relevant covariates. However, the field of survival analysis has many robust statistical methods for calculating survival probabilities, among other quantities, which incorporate time-dependent covariates. See, for instance, the Cox proportional hazards model [34] or a more general “direct” hazard model [29, 30]. These methods should be implemented for disruption prediction and compared to current ML algorithm performance. In addition, other established fields, like probabilistic risk assessment, should be explored as alternative methods for disruption prediction and avoidance. 2. 2.

The analysis in this work relied on the assumption of a well-calibrated disruptivity signal; that is, a value of $\mathrm{P}_{\mathrm{D}}=0$ or 1 was assumed to truly indicate that the plasma state resides in the non-disruptive or disruptive class, respectively. Future work with ML algorithms should take care to properly calibrate output signals, for instance as described in [29, 30, 35]. This approach might also be applied to other ML methods trained on two classes and with output signals in the range $[0,1]$, such as those described in [24, 13, 8, 15, 9, 22]. 3. 3.

Additionally, the physical meaning of ML algorithm output signals must be well understood. Fortunately for RF algorithms trained on two-class data, this interpretation is relatively straightforward. However, the “distance to a boundary” as calculated by a Support Vector Machine model, for example, might not be as clearly interpreted. 4. 4.

Another assumption of this analysis was a uniform probability density of times of data in the disruptive class. However, in principle, the probability density could be learned from the database itself. One extension of this work would be to consider more carefully the sampling rate of data in the database and/or to adjust the ML algorithm to remove any related biases. 5. 5.

Finally, the linear extrapolation of disruptivity employed in this paper is quite simple (even crude), especially considering that plasma parameters can be actuated in real time to navigate in “disruptivity space.” A smarter method, such as that described in [36], could be used to calculate gradients in parameter space and map possible trajectories away from the disruptive boundary ($\mathrm{P}_{\mathrm{D}}=1$, in this case). Furthermore, dynamical models of the the plasma state vector $\mathbf{x}(t)$—which is input into the disruption predictor, i.e. $\mathrm{P}_{\mathrm{D}}[\mathbf{x}(t)]$—could be implemented to more realistically extrapolate the future disruptivity.

Acknowledgements

The authors thank E. Olofsson and M. Parsons for enlightening discussions, as well as the entire Alcator C-Mod team. Comments from the reviewers also strengthened this paper. This work was supported by US DOE Grant DE-FC02-99ER54512, using Alcator C-Mod, a DOE Office of Science User Facility.

R. Sweeney is supported by the U.S. Department of Energy Fusion Energy Sciences Postdoctoral Research Program administered by the Oak Ridge Institute for Science and Education (ORISE) for the DOE. ORISE is managed by Oak Ridge Associated Universities (ORAU) under DOE contract number DE-SC0014664. All opinions expressed in this paper are the authors’ and do not necessarily reflect the policies and views of DOE, ORAU, or ORISE.

Appendix A Calculation of the minimum expected future lifetime

In this section, it is shown that the expected future lifetime $\tau$, as calculated from (12), does not fall far below the class time $\Delta t_{\rm c}$, chosen to bifurcate data into non-disruptive and disruptive classes, for most cases of interest. Recall that the probability of survival beyond time $t+t^{\prime}$, assuming survival until time $t$, is given by (9)

[TABLE]

where the future time $t^{\prime}$ has been split into $n$ steps. As described in section 3, equation (15) is the Kaplan-Meier estimator of survival probability, where the probability of failure (i.e. a disruption) is intuited to be the product of the disruptivity signal $\mathrm{P}_{\mathrm{D}}(t)$ at each time step and the fractional time interval $t^{\prime}/n\Delta t_{\rm c}$ from (2.1). It is important to emphasize again that the time step must be less than or equal to the class time, i.e. $t^{\prime}/n\leq\Delta t_{\rm c}$, for this relation to hold. This can be achieved by the appropriate choice of $n$ for given $t^{\prime}$ and $\Delta t_{\rm c}$; in practice, the sampling time $\Delta t_{\rm s}=t^{\prime}/n$ of a disruption predictor is often much shorter than $\Delta t_{\rm c}$ anyway.

Because the survival function is always non-negative ($S(t)\geq 0$) and monotonically-decreasing ($\mathrm{d}S/\mathrm{d}t\leq 0$), the expected future lifetime, from (12), satisfies

[TABLE]

where $T$ is assumed finite. Before using (15) as the integrand of (16), note that $S(t)$ is minimal for $\mathrm{P}_{\mathrm{D}}(t)=1$; therefore, we can write

[TABLE]

Here, it is important to note that (i) the integral is over the domain $t^{\prime}\in[0,T]$ and (ii) $t^{\prime}/n\leq\Delta t_{\rm c}$ must still be satisfied. Thus, the upper limit is maximally $T=n\Delta t_{\rm c}$. If we wanted to evaluate the integral for $T>\Delta t_{\rm c}$, we would need to use a different form of (15) in accordance with (2.1), but that is not necessary here. The integral of (17) is evaluated to be

[TABLE]

For the lowest possible sampling rate, $n=1$ and then $\tau\geq\Delta t_{\rm c}/3$ for $t^{\prime}\leq\Delta t_{\rm c}$. However, for most cases of interest, $n\gg 1$ since $\Delta t_{\rm s}\ll\Delta t_{\rm c}$; then we arrive at the desired result $\tau\gtrsim\Delta t_{\rm c}$.

As a final remark, note that the survival function of (15) approaches the exponential $\exp(-t^{\prime}/\Delta t_{\rm c})$ in the limit $n\to\infty$ when $\mathrm{P}_{\mathrm{D}}(t)=1$. In this limit, the first integral of (16) can be carried out explicitly so that the expected future lifetime is exactly $\tau=\Delta t_{\rm c}$. It may seem counter-intuitive that the calculated future lifetime is not less than the class time, i.e. $\tau<\Delta t_{\rm c}$, when the plasma is always predicted to be in the disruptive state, i.e. when $\mathrm{P}_{\mathrm{D}}=1$. However, this is a consequence of (i) the uncertainty associated with a choice of class time $\Delta t_{\rm c}$ in binary classification of the data sets, (ii) a sampling time much shorter than $\Delta t_{\rm c}$ which provides frequent updates of the current plasma state, and (iii) the framework adopted in this paper which treats probabilities between time steps as independent. (Refer to sections 2 and 3 for further discussion.) Future work can explore the inclusion of past data in the prediction of future survival. For instance, perhaps a prediction of $\mathrm{P}_{\mathrm{D}}=1$ over the past $\Delta t_{\rm c}$ seconds should automatically trigger an alarm. In the end, $\tau$ may not even be the most appropriate metric, or quantity to monitor, in successive studies since the survival probability at future time $\tau$ can often be less 50%, i.e. $S(t+\tau|t)<0.5$, as seen in figures 2-6.

Appendix B Calculation of the maximum hazard

In this section, the (approximate) maximum value of the hazard function is calculated. From (4), the hazard function is computed

[TABLE]

Once again, the survival function is given by

[TABLE]

where it is assumed that the initial time is $t=0$ and $t/n\leq\Delta t_{\rm c}$. Taking the logarithm of (20) turns the product into a summation

[TABLE]

Here, we utilize the linear extrapolation from (10)

[TABLE]

where $\mathrm{P}_{\mathrm{D}0}=\mathrm{P}_{\mathrm{D}}(0)$ and $\mathrm{d}\mathrm{P}_{\mathrm{D}0}/\mathrm{d}t=\mathrm{d}\mathrm{P}_{\mathrm{D}}(0)/\mathrm{d}t$ are evaluated at the starting time, and the range of values is restricted to $\mathrm{P}_{\mathrm{D}}\in[0,1]$. Substituting (22) into (21) and taking the derivative gives

[TABLE]

This relation is maximal for $\mathrm{P}_{\mathrm{D}}(jt/n)\to 1$. Taking this limit and evaluating the summation gives

[TABLE]

Because $\mathrm{P}_{\mathrm{D}}\to 1$, it follows that the middle term of (24) is (maximally) 1, so that

[TABLE]

For the lowest sampling rate, $n=1$ and then $h(t)\leq 2/(\Delta t_{\rm c}-t)$, which can approach infinity for $t\to\Delta t_{\rm c}$. However, for more realistic cases with a sampling time $\Delta t_{\rm s}=t/n$ and $n\gg 1$, we find

[TABLE]