Forecasting Industrial Aging Processes with Machine Learning Methods

2020·Arxiv

Abstract

Abstract

Accurately predicting industrial aging processes makes it possible to schedule maintenance events further in advance, ensuring a cost-efficient and reliable operation of the plant. So far, these degradation processes were usually described by mechanistic or simple empirical prediction models. In this paper, we evaluate a wider range of data-driven models, comparing some traditional stateless models (linear and kernel ridge regression, feed-forward neural networks) to more complex recurrent neural networks (echo state networks and LSTMs). We first examine how much historical data is needed to train each of the models on a synthetic dataset with known dynamics. Next, the models are tested on real-world data from a large scale chemical plant. Our results show that recurrent models produce near perfect predictions when trained on larger datasets, and maintain a good performance even when trained on smaller datasets with domain shifts, while the simpler models only performed comparably on the smaller datasets.

Keywords: Machine Learning, Time Series Prediction, Predictive Maintenance, Catalyst Degradation.

1. Introduction

Aging of critical assets is an omnipresent phenomenon in any production environment, causing significant maintenance expenditures or leading to production losses. The understanding and anticipation of the underlying degradation processes is therefore of great importance for a reliable and economic plant operation, both in discrete manufacturing and in the process industry.

With a focus on the chemical industry, notorious aging phenomena include the deactivation of heterogeneous catalysts [1] due to coking [2], sintering [3], or poisoning [4]; plugging of process equipment, such as heat exchangers or pipes, on process side due to coke layer formation [5] or polymerization [6]; fouling of heat exchangers on water side due to microbial or crystalline deposits [7]; erosion of installed equipment, such as injection nozzles or pipes, in fluidized bed reactors [8, 9]; corrosion inside pipes or vessels [10]; and more.

Despite the large variety of affected asset types in these examples, and the completely different physical or chemical degradation processes that underlie them, all of these phenomena share some essential characteristics:

1. The considered critical asset has one or more key performance indicators (KPIs), which quantify the progress of degradation.

2. On a time scale much longer than the typical production time scales (i.e., batch time for discontinuous processes; typical time between set point changes for continuous processes), the KPIs drift more or less monotonically to ever higher or lower values, indicating the occurrence of an irreversible degradation phenomenon. (On shorter time scales, the KPIs may exhibit fluctuations that are not driven by the degradation process itself, but rather by varying process conditions or background variables such as, e.g., the ambient temperature.)

3. The KPIs return towards their baseline after maintenance events, such as cleaning of a fouled heat exchanger, replacement of an inactive catalyst, etc.

4. The degradation is no ‘bolt from the blue’ – such as, e.g., the bursting of a flawed pipe –, but is rather driven by creeping, inevitable wear and tear of process equipment.

Any aging phenomenon with these general properties is addressed by the present work.

Property (4) suggests that the evolution of a degradation KPI is to a large extent determined by the process conditions, and not by uncontrolled, external factors. This sets the central goal of the present work: To forecast the evolution of the degradation KPI over a certain time horizon, given the planned process conditions in this time frame. If instead external factors such as humidity, ambient temperature, human interventions, etc. are the dominating driving forces of an aging process, the presented approach is not applicable.

For virtually any important aging phenomenon in chemical engineering, the respective scientific community has developed a detailed understanding of their microscopic and macroscopic driving forces. This understanding has commonly been condensed into sophisticated mathematical models. Examples of such mechanistic degradation models deal with coking of steamcracker furnaces [11, 12, 13], sintering [14, 15] or coking [16] of heterogeneous catalysts, or crystallization fouling of heat exchangers [17].

While these models give valuable insights into the dynamics of experimentally non-accessible quantities, and can help to verify or falsify hypotheses about the degradation mechanism in general, they are usually not (or only with sig-nificant modeling effort) transferable to the specific environment in a real-world apparatus: Broadly speaking, they often describe ‘clean’ observations of the degradation process in a lab environment, and do not reflect the ‘dirty’ reality in production, where additional effects come into play that are hard or impossible to model mechanistically. To mention only one example, sintering dynamics of supported metal catalysts are hard to model quantitatively even in the ‘clean’ system of Wulff-shaped particles on a flat surface [15] – while in real heterogeneous catalysts, surface morphology and particle shape may deviate strongly from this assumption. Another disadvantage of mechanistic models is that their numerical solution can be computationally expensive or even intractable.

While the latter issue of computational complexity can be mitigated by surrogate modeling methods [18, 19], the former issue of real-world complexity is the main reason why mechanistic models of degradation dynamics are rarely used in a production environment.

This weakness is addressed by hybrid process models (also referred to as gray-box models), which combine mechanistic with data driven modeling approaches to bridge the gap between idealized mechanistic models and the real world [20, 21, 22, 23, 24]. However, to the best of our knowledge, these models have not yet been used to forecast industrial aging processes in chemical plants.

When dealing with the complexity of modeling real-world aging processes, statistical approaches have proven successful in a variety of applications. For example, data-driven methods for fault detection of chemical plants [25], such as multivariate anomaly detection with fisher discriminant [26] or principal component analysis [27], are routinely applied nowadays to monitor process equipment. However, we emphasize that most of these applications focus on the detection or monitoring of degradation, not on the prediction of its progression.

Publications of data-driven models that predict degradation dynamics predominantly address specific aging phenomena, such as batch-to-batch fouling of heat exchangers as a function of the polymer type produced in the respective batch [6], or fouling of the crude preheat train in petroleum refineries [28, 29]. To the best of our knowledge, all published work in this field is either based on classical statistical regression methods, such as ordinary or partial least squares [6], or on small-scale machine learning (ML) methods, such as small feed-forward neural networks (FFNN) trained with limited datasets [28, 29]. So far, advanced ML algorithms, such as recurrent neural networks (RNN), trained with years or even decades of historical plant data, have not been studied in depth in the context of predicting degradation of chemical process equipment [30]. It is the aim of this work to investigate the prospects of advanced ML methods for this problem, compare them to classical regression methods and understand potential limitations.

The rest of this paper is structured as follows: First, we formalize the general

IAP problem setting (Section 2). Then we describe our two datasets (Section 3), as well as quickly introduce the five ML models that we evaluated for this task (Section 4). Finally, we present the prediction results of the different models on both datasets (Section 5) and conclude the paper with a discussion (Section 6).

2. Problem Deﬁnition

The general industrial aging process (IAP) forecasting problem is illustrated in Fig. 1: The aim is to model the evolution of one or several degradation KPIs y(t) as a function of the planned process conditions x(t) over the time frame of an entire degradation cycle, i.e., until the next required maintenance event. Formally, the problem statement reads:

Given. A time window [0] between two maintenance events, referred to as the i-th degradation cycle, and the planned process conditions within this time window.

Determine. An estimate ˆof the evolution of one or several degradation KPIs ) in the i-th cycle.

Objective. Minimal discrepancy between the estimated and true KPI values, as measured by mean squared error (MSE), cf. Eq. (2).

Figure 1: Illustration of the industrial aging process (IAP) forecasting problem. The degradation KPI (e.g. pressure drop ∆p in a fixed bed) increases over time (e.g. due to coking), influenced by the (manually controlled) process conditions (e.g. reaction temperature T and flow rate F). The KPI recovers after a maintenance event, which segments the time axis into different degradation cycles. The IAP forecasting task is to predict the evolution of the KPI, i.e., the target (dependent) variable ), in the current cycle i, given the upcoming process conditions, i.e., the input (independent) variables

Some peculiarities inherent to IAPs make this problem more challenging than it first appears. First, degradation phenomena can exhibit pronounced memory effects, which means that a certain input pattern x(t) may affect the output ) only at much later times . In addition, these memory effects can also occur across multiple time scales, which makes these processes notoriously hard to model. As an example, consider a heat exchanger suffering from coking of the inner tube walls. The observed heat transfer coefficient serves as KPI ), and process conditions ) comprise mass flow, chemical composition and temperature of the processed fluid. The time horizon is one cycle between two cleaning procedures (e.g. burn-off). If at an early time in the cycle an unfavorable combination of low mass flow, high content of coke precursors, and high temperature occurs, first coke patches can form at the wall, which are not yet big enough to impact heat transfer significantly. However, they serve as a nuclei for further coke formation later in the cycle, so that ) drops faster at compared to a cycle where the process conditions were not unfavorable around , but with very similar process conditions throughout the rest of the cycle.

An additional complication arises from the fact that in real application cases, the distinction between degradation KPI y, process conditions x, and uncontrolled influencing factors is not always clear-cut. Consider, for example, the case of a heterogeneous catalyst subject to deactivation, where the loss of catalytic activity leads to a decreased conversion rate. In this case, the conversion rate could serve as a target degradation KPI y, while process conditions, such as the temperature, which are manually controlled by the plant operators, would be considered input variables x for the model. However, the plant operators might try to keep the conversion rate at a certain set point, which can be achieved by raising the temperature to counteract the effects of the catalyst degradation. This introduces a feedback loop between the conversion rate and the temperature, which means the temperature can not be considered an independent variable anymore, as its actual value (partially) depends on the target. Therefore, care has to be taken, since including such a dependent variable as an input x in a model could lead one to report overly optimistic prediction errors that would not hold up when the model is later used in reality.

3. Datasets

To gain insights into and evaluate different ML models for the IAP forecasting problem, we consider two datasets: one synthetic, which we generated ourselves using a mechanistic model, and one containing real-world data from a large plant at BASF. Both datasets are described in more detail below.

The reason for working with synthetic data is that this allows us control two important aspects of the problem: data quantity and data quality. Data quantity is measured, e.g., by the number of catalyst lifecycles in the dataset, which can be chosen as large as we want for synthetic data, to test even the most data-hungry ML methods. Data quality refers to the level of noise in the dataset, or, in other words, the degree to which the degradation KPI y(t) is uniquely determined by the provided process conditions x(t) in the dataset. In a synthetic dataset based on a deterministic degradation model, we know that there is a functional mapping between x and y, i.e., there exists no fundamental reason that could prevent a ML model from learning this relation with vanishing prediction errors. In contrast, with real data, a bad prediction error can either be a problem of the method, and/or of the dataset, which might not contain sufficient information on the input side x to accurately predict the output quantity y.

Please note that the synthetic dataset refers to a different chemical process than the real-world dataset. As a consequence, models trained on the synthetic data have to be re-trained when applying them to the real-world dataset.

3.1. Synthetic dataset

For the synthetic dataset, we modeled the wide-spread phenomenon of slow, but steady loss of catalytic activity in a continuously operated fixed-bed reactor. Ultimately, the catalyst deactivation leads to unacceptable conversion or selectivity rates in the process, necessitating a catalyst regeneration or replacement, which marks the end of one cycle.

The chemical process in the reactor under consideration is the gas-phase oxidation of an olefine. To generate the time series for all variables, we used a mechanistic process model with the following ingredients (further details can be found in Section A in the supplement):

• Mass balance equations for all five relevant chemical species (olefinic reactant, oxygen, oxidized product, CO, water) in the reactor, which is, for simplicity, modeled as an isothermal plug flow reactor, assuming ideal gas law. The reaction network consists of the main reaction (olefine + Oproduct) and one side reaction (combustion of olefine to CO).

• A highly non-linear deactivation law for the catalyst activity, which depends on reaction temperature, flow rate, and inflowing oxygen, as well as the activity itself.

• Kinetic laws for the reaction rates.

• A stochastic process determining the process conditions (temperature, flow rate, etc.).

Based on the current process conditions and hidden states of the system, we used the mechanistic model to generate a multivariate time series [x(t), y(t)] comprising 50 years of historical data, separated into 2153 degradation cycles with a total of 435917 time points. The final dataset consists of six input features x(t) at each time point t, namely the five process parameters (mass flow rate, reactor pressure, temperature, and mass fractions of the two reactants olefine and O) and the time since the last maintenance event, while the two degradation KPIs (conversion and selectivity) constitute the target variables y(t) (Table 1).

To give an impression of the simulated time series, one month of data is shown in Fig. 2. The duration of deactivation cycles is around 8-10 days. The catalyst activity A(t) is a hidden state and therefore not part of the dataset, but is only shown to illustrate the dynamics of the problem: System output y(t) (selectivity and conversion) is not only affected by the current process parameters x(t), but also the current catalyst activity A(t), which is non-linearly decreasing over each cycle.

Table 1: Variables in the synthetic dataset used for machine learning. The variable type indicates if a quantity is an input (x) or output (y) in the IAP forecasting problem.

Figure 2: One month of the synthetic dataset, showing loss of catalytic activity in a fixed-bed reactor. At each time point t, the vector of process conditions x(t) includes the reactor temperature T, mass flow rate F, reactor pressure p, and mass fractions of the reactants at the reactor inlet. Degradation KPIs y(t) are conversion and selectivity of the process.

3.2. Real-world dataset

The second dataset contains process data for the production of an organic substance in a continuous world-scale production plant at BASF. The process is a gas phase oxidation in a multi-tubular fixed-bed reactor.

The catalyst particles in the reactor suffer from coking, i.e., surface depo-

sition of elementary carbon in the form of graphite. This leads to reduced catalytic activity and increased fluid resistance. The latter is the more severe consequence and leads to an increasing pressure drop over the reactor, as measured by the difference ∆p of gas pressure before and after the reactor. The physical reason for the increasing pressure drop is the reduction of void fraction in the catalyst bed due to the build-up of coke residues, as described by Ergun’s equation [31] and its extensions [32].

When ∆p exceeds a pre-defined threshold, the so-called end-of-run (EOR) criterion is reached. Then, the coke layer is burned off in a dedicated regeneration procedure, by inserting air and additional nitrogen into the reactor at elevated temperatures for a variable number of hours. Operational reasons can lead to a delayed burn-off with ∆p exceeding the EOR threshold, or, vice versa, a premature burn-off when ∆p has not yet reached the EOR threshold. Some exemplary cycles for ∆p are shown in Fig. 3.

Figure 3: One month of historic data of the real-world dataset, showing the pressure loss ∆p over the reactor, which is the degradation KPI y(t) in this IAP forecasting problem. When ∆p reaches a value in the order of the EOR threshold of 70 mbar, the cokes deposit is burned off, which marks the end of a cycle.

Since coke is not removed perfectly by this burn-off procedure, coke residues accumulate from regeneration to regeneration, making the pressure drop issue ever more severe. Therefore, the entire catalyst bed must be replaced every 6-24 months.

Suspected influencing factors for the coking rate are:

1. mass flow F through the reactor (“feed load”) 2. ratio of organic reactant to oxygen in the feed 3. intensity of previous regeneration procedures

4. length of the previous degradation cycle

The dataset includes seven years of process data from the four most relevant sensors, extracted from the plant information management system (PIMS) of the plant. Given the time scale of 4 to 7 days between two burn-off procedures, this corresponds to 375 degradation cycles belonging to three different catalyst batches. The sampling rate is 1/hour for all variables with a linear interpolation to that time grid. After removing some outlier cycles (shorter than 50 hours), the final size of the dataset is 327 cycles for a total of 36058 time points.

The task is to predict, at an intermediate moment during a degradation cycle, the coking-induced pressure drop ∆p over the entire remaining duration of the cycle. Of particular interest is a prediction of the time point at which the EOR threshold ∆= 70 mbar is reached. The input variables x(t) of the model consist of the process parameters, as well as several engineered features (Table 2), either derived from the process parameters themselves or from the degradation KPI ∆p in the previous cycles, yielding a total of 10 input features.

Table 2: Variables in the real-world dataset used for machine learning. The variable type indicates if a quantity is an input (x) or output (y) in the IAP forecasting problem. The first four input variables are the process parameters (measured with sensors in the plant), all others are additional engineered features. The total inflowing reactant F R is a mixture of fresh reactant from a tank reservoir (F FRESH), and recycled reactant from downstream separation units. The variable F FRESH is not direclty used as an input feature, but is necessary for the definition of F FRESH / F R.

4. Machine Learning Methods

We now frame the IAP forecasting problem described in Section 2 as a machine learning problem, by defining a concrete function f that returns ˆ), an estimate of the KPIs at a time point t in the i-th degradation cycle, based on the process conditions at this time point as well as possibly up to k hours before t:

The task is to predict ) for the complete cycle (i.e., up to ), typically starting from about 24 hours after the last maintenance event that concluded the previous cycle.For simplicity, in many cases we only write x and y, omitting the reference to the current cycle i and time points t in questions, while x might include the process conditions for multiple time points from a fixed time window in the past (i.e. up to ).

In this paper we examine five different machine learning models, which can be divided into two main subgroups: stateless and stateful models (Fig. 4).

Figure 4: A comparison of stateless and stateful models for time series forecasting. (a) shows a stateless model, which bases the predictions on the information contained in a fixed time window in the past, while (b) displays a stateful model, where information about the past is maintained and propagated using a hidden state.

Stateless models directly predict the target given the inputs from a fixed time window in the past, independent of the predictions for previous time points. The stateless models used in this paper are linear ridge regression (LRR), kernel ridge regression (KRR), and feed-forward neural networks (FFNN), i.e., one linear and two non-linear prediction models. Stateful models, on the other hand, maintain an internal hidden state of the system that encodes information about the past and which is utilized in addition to the current process conditions when making a prediction. The stateful models we explore in this paper are two recurrent neural network (RNN) models [33]: echo state networks (ESN) and long short term memory (LSTM) networks. The five ML models are briefly introduced in the following subsections, while a more detailed description can be found in Section B in the supplement.

As the models utilize the inputs differently when making their predictions, the kind of machine learning method that is chosen for the forecasting task determines the exact form of the function f from Eq. (1). Yet, while the functional form of f is predetermined, its exact parameters need to be adapted to fit the dataset at hand in order to yield accurate predictions. For this, the available data is first split into training and test sets, where each of the two sets contains the entire multivariate time series from several mutually exclusive degradation cycles from the original dataset, i.e., multiple input-output pairs consisting of the planned conditions x and degradation KPIs y of the given process. Then, using the data in the training set, the ML algorithm learns the optimal parameters of f by minimizing the expected error between the predicted KPIs ˆ) and the true KPIs ) [34, 35, 36]. After the ML model has been trained, i.e., when f predicts ) as accurately as possible on the training set, the model is evaluated on the test set, i.e., previously unseen data, to give an indication of its performance when later used in reality. If the performance on the training set is much better than on the test set, the model does not generalize well to new data and is said to have “overfit” on the training data.

In addition to the regular parameters of f, many ML models also require setting some hyperparameters, that, for example, determine the degree of regularization (i.e., how much influence possible outliers in the training set can have on the model parameters). The values for these hyperparameters are determined by training models with different settings on a subset of the training data and evaluating them on the remaining data from the training set (also called the validation set) and then choosing those hyperparameters with the best performance on this validation split. The final model with the selected hyperparameters is then trained on the complete training set and evaluated on the test set as described above. Specific details about the hyperparameters of the different models and their chosen values for each dataset can be found in Section C in the supplement.

4.1. Linear ridge regression (LRR)

LRR is based on the standard linear regression model, but with an added regularization term that prevents the weights from taking on extreme values due to outliers in the training set. The target variables y are predicted as a linear combination of the input variables x, i.e.,

where is a weight matrix, i.e., the model parameters of f that are learned from the training data. The simple model architecture, globally optimal solution, and regularization of LRR all contribute to reducing overfitting of the model. Additionally, training and evaluating the model is computationally cheap, making it a viable model for large amounts of data as well.

The regularization parameter is the only hyperparameter of LRR that needs optimization, which we select by searching for the best performing hyperparameter on a grid from 0.001 to 10 using a 10-fold cross-validation [37].

For the synthetic dataset, the input vector at time point t consists of the process parameters for the past 24 hours, giving the model a time window into the past, i.e. ) = [1); 24)]. For the real-world dataset, the input at time point t only consists of the process conditions at that time point x(t), as a larger time window would reduce the number of data points in the training set (if k hours from the past are taken, the inputs for each cycle have to start k hours later, leading to the loss of k samples per cycle) while increasing the number of input features, therefore making overfitting more likely.

4.2. Kernel ridge regression (KRR)

KRR is a non-linear regression model that can be derived from LRR using the so called ‘kernel trick’ [38, 35, 34, 39, 40]. Instead of using the regular input features x, the features are mapped to a high (and possibly infinite) dimensional space using a feature map , corresponding to some kernel function k such that ) = ). By computing the non-linear similarity k between a new data point x and the training examples for j = 1, . . . , N, the targets y can be predicted as

where are the learned model parameters.

The non-linear KRR model can adapt to more complex data compared to LRR while the globally optimal solution can still be obtained analytically, which has made KRR one of the most commonly used non-linear regression algorithms in the field of machine learning [41, 42]. However, the performance of the model is also more sensitive to the choice of hyperparameters, so a careful selection and optimization of the hyperparameters is necessary. Additionally, the fact that solving the optimization problem scales cubically in runtime and quadratically in memory with the number of training examples N makes it difficult to apply KRR to problems with large training sets. The scaling issue is exactly one of the challenges we encountered when applying KRR to the synthetic dataset. The linear algebra library we used only supported operations with matrices of a maximum size of 50 000, so when applying KRR to the synthetic dataset we only used 10 percent of the training set, consisting of around 40 000 points. The KRR model has two hyperparameters that need to be optimized, the regularization parameter and the kernel width . The optimal hyperparameters were selected with a 10-fold cross-validation from an exponential grid, with values ranging from 10to 1 for and from 0.1 to 1000 for .

Analogously to LRR, the KRR model for the synthetic dataset was trained using inputs from a 24-hour time window into the past, while the model for the real-world dataset only used inputs at the current time point.

4.3. Feed-forward neural networks (FFNN)

FFNNs were the first and most straightforward neural network architecture to be conceived, yet, due to their flexibility, they are still successfully applied to many different types of machine learning problems [43, 44, 45]. Analogously to LRR, FFNNs learn a direct mapping f between some input features x and some output values y. However, unlike a linear model, FFNNs can approximate also highly non-linear dependencies between inputs and outputs. This is achieved by transforming the input using a succession of “layers”, where each layer is usually composed of a linear transformation with a weight matrix W followed by a non-linear operation :

The values of the weights are optimized to minimize the squared loss between the predicted value and the target using error backpropagation [46]. Since the loss function is highly non-convex, the optimization procedure only finds a local minimum, in contrast to the globally optimal solution found by LRR and KRR, though this is usually irrelevant to the performance of FFNNs [47] (see supplementary Section B.1 for more details).

The FFNN has a series of hyperparameters that need to be optimized, ranging from hyperparameters about the model architecture itself (e.g. size and number of layers), to hyperparameters relating to the training procedure, where we used stochastic gradient descent (SGD) with Nesterov momentum [48].

Due to the larger number of hyperparameters to optimize, we selected the best hyperparameters by performing a grid search on a small number of values for each hyperparameter. The optimal hyperparameters were determined based on the performance on a validation set consisting of a random selection of 15% of the cycles in the training set. The number of training epochs was chosen using early stopping, with training being stopped if the validation error had not improved in the last six epochs in the case of the synthetic dataset, while a stopping criterion of 30 epochs was used for the real-world dataset due to the overall smaller size of the training set.

As with the other stateless models, the FFNN was trained using inputs from a 24-hour time window into the past for the synthetic dataset, and inputs only at the current time point for the real-world dataset.

4.4. Echo state networks (ESN)

ESNs are an alternative RNN architecture that can alleviate some of the training related problems of RNNs (see Section B.2 in the supplement) by not using error backpropagation for training at all [49]. Instead, ESNs use very large randomly initialized weight matrices combined with a recurrent mapping of the past inputs; collectively called the “reservoir”. This way, ESNs can keep track of the hidden state (with ) of the system by updating h(t) as a combination of the previous hidden state 1) and the current input vector x(t). The final prediction of the output is then computed using LRR on the inputs and hidden state, i.e.,

In general, echo state networks are a very powerful type of RNN, whose performance on dynamical system forecasting is often on par with or even better than that of other, more popular and complex RNN models (LSTM, GRU, etc.) [49, 50]. Since the only learned parameters are the weights of the linear model used for the final prediction, ESNs can also be trained on smaller datasets without risking too much overfitting.

For the ESN model, which had the largest number of hyperparameters (see supplementary Sections B.2 and C), we used a random search to select candidate hyperparameter sets, which were then evaluated using a 10-fold cross-validation on the training set.

As the ESN is a stateful model, capable of encoding the past in its hidden state, the input vector at any time point t only consists of the process conditions at the current time point, i.e. x(t).

4.5. LSTM networks

Another very popular architecture for dealing with the vanishing gradients problem of RNNs is the long short term memory (LSTM) architecture, which was developed specifically for this purpose [51]. LSTMs are trained using error backpropagation as usual, but avoid the problem of vanishing gradients by using an additional state vector called the “cell state”, alongside the usual hidden state, which is updated slowly via special layers called “gates”, and thus is capable of preserving long-term information. While the updates of the hidden state h(t) of an LSTM network are much more complex compared to ESNs, the final prediction is again only a linear transformation of the network’s internal hidden state:

However, in this case, the parameter values of are optimized together with the other parameters of the LSTM network, instead of using a separate LRR model.

Due to the multiple layers needed to model the gates that regulate the cell state, the LSTM typically requires larger amounts of training data to avoid overfitting. Despite its complexity, however, the stability of the gradients of the LSTM make it very well suited for time series problems with long ranging dependencies [52, 53, 54, 55, 56].

The LSTM hyperparameters are very similar to that of the FFNN and were again selected with a grid search on a small number of values. Analogously to the FFNN, the different hyperparameter sets were evaluated on a validation set consisting of a random selection of 15% of the cycles in the training set, and the number of training epochs was determined using early stopping, i.e., training stopped when the validation error did not improve for six epochs for the synthetic dataset and 30 epochs for the real-world dataset.

Since the LSTM is also capable of encoding the past into a hidden state, the input vector at any time point t only consists of the process conditions at the current time point, i.e. x(t).

5. Results

In this section, we report our evaluation of the five different ML models introduced in Section 4 using the synthetic and real-world datasets described in Section 3. To measure the prediction errors of the ML models, we use the mean squared error (MSE), which, due to the subdivision of our datasets into cycles, we define slightly differently than usual: Let the dataset D be composed of N cycles, and let ) denote the KPIs at time point within the i-th cycle, where is the length of the i-th cycle. Then, given the corresponding model predictions ˆ), the MSE of a model for the entire dataset is calculated as

Since the synthetic and real-world datasets are very different, they were used to examine different aspects of the models. The synthetic dataset was used to examine how the models perform in a nearly ideal scenario, where data is freely available and the noise is very low or even non-existent. On the other hand, the real-world dataset was used to test the robustness of the models, since it contains only a limited amount of training samples and a relatively high noise level.

5.1. Synthetic dataset

The synthetic dataset was randomly split into a training set (1938 cycles consisting of all together 391876 time points) and a test set (215 cycles / 44041 time points). For all models, the inputs were all scaled to have zero mean and unit variance, while the KPIs were left unscaled. Only results for conversion as a degradation KPI are discussed; results for selectivity follow the same trends.

Fig. 5 shows the normalized mean squared errors (MSEs) for each of the models on the training and test splits when trained on differently sized subsets of the full training set. For most of the models, the error converges relatively early, meaning that even with a fraction of the complete dataset, the models manage to learn an accurate approximation of the dynamics of the synthetic dataset, as far as the respective model complexity permits. This also indicates that the existing errors in the models are largely due to the limitations on the flexibility of the models themselves, and not due to the training set not being

Figure 5: Training and test set MSEs for the five different models, where the models were trained on differently sized subsets of the training set from the synthetic dataset. The performance for KRR cuts off early as it is too computationally expensive to train KRR with the roughly 2000 degradation cycles, corresponding to 400000 time points.

large enough. This is clearly evident with LRR, which essentially achieves its maximum performance using 5% of the total dataset size. Since LRR is a linear model, it can only learn the linear relations between the inputs and outputs. While this high bias prevents the model from learning most of the non-linear dynamics regardless of the training set size, this also means that the model has low variance, i.e., it tends not to overfit on the training data [57]. For the FFNN, the error slowly declines as the number of samples increases, though at an ever slower rate, with the error using the full training dataset being significantly lower that LRR.

As for ESN and LSTM, both methods seem to somewhat overfit for the smaller training set sizes, as indicated by the differences between training and test errors, however, even in the small data regime the test errors are much lower compared to the three stateless models. The errors of both stateful models converge at around 50% of the full dataset, after which there is virtually no overfitting and no significant improvement of the performance for larger dataset sizes.

The general lack of overfitting can be partially attributed to the lack of noise in the synthetic dataset, since overfitting usually involves the model fitting the noise instead of the actual signal/patterns. Additionally, the low amount of overfitting in the small data regime also points to the fact that such time-dependent stateful models are well-suited for this problem setting.

Across all dataset sizes, the LSTM model is clearly the best performing, with its error when using the full dataset being 5 times smaller than the error of the ESN model. This trend is also confirmed by other performance metrics (Table 3). Additionally, the scale invariant measures (normalized MSE and ) show that all of the models have relatively small errors compared to the variance of the dataset.

Given the great performance of the ESN and especially the LSTM model, these experiment clearly demonstrate that even with smaller amounts of highquality data, entire degradation cycles can in principle be predicted with very high accuracy.

Table 3: Regression performance metrics for the different models on the synthetic dataset. ∗Note that the KRR model was trained using only 10% of the full training data.

Fig. 6 shows plots of the true and predicted conversion rates of the different models for some randomly selected cycles from the training and test sets. These show that all the models are capable of accurately predicting the instantaneous effects of the input parameters on the output, since this relation is largely linear and not time dependent. However, where the models differ the most is in the non-linear long term degradation, where the stateless models only predict a roughly linear trend, with FFNN coming slightly closer to the actual degradation trend due to its non-linearity, while the ESN model predicts the degradation better but fails to capture the rapid decline near the end of each cycle. The LSTM model, on the other hand, manages to capture the short and long term effects almost perfectly, with only small errors at the very ends of the cycles where there is smaller amounts of data, due to the varying length of the cycles. This trend is also noticeable on the average true and predicted KPIs in Figure 7, where the LRR and FFNN models have larger differences especially in the last part of the curves, where the memory effects would be more pronounced, while the mean of the ESN and LSTM predictions have almost perfect alignment with the mean KPI.

5.2. Real-world dataset

The real-world dataset is much smaller than the synthetic, consisting of a total of 327 cycles (after outlier removal, see Section 3) with all together 36058 time points. As the real-world dataset stretches over 3 time periods with different catalyst charges in the reactor, we test the performance in a realistic manner by selecting the third catalyst charge as the test set, which makes it possible to see to what extent the models are able to extrapolate across the different conditions caused by the catalyst exchange. This resulted in a training set that is more than 10 times smaller than the training set of the synthetic dataset, consisting of 256 cycles (28503 time points), while the test set consists of 71 cycles (7555 time points). Analogously as with the synthetic dataset, the

Figure 6: Comparison of the true and predicted conversion rates across four ML models for some random training and test cycles from the synthetic dataset.

Figure 7: Comparison of the means of the true and predicted conversion rates across four ML models for the training and test cycles of the synthetic dataset.

Figure 8: Training and test set MSEs for the five different models evaluated on the real-world dataset.

input variables were all scaled to have zero mean and unit variance, while the KPIs were left unscaled.

Fig. 8 shows the normalized mean squared errors for each of the models on the training and test sets. Due to the larger amounts noise and smaller dataset size, as well as the subtle differences between the training and test distributions due to the different catalyst beds, the results here are different compared to the synthetic dataset: The more complex models show more overfitting, indicated by the test errors being significantly larger than the corresponding training errors, especially for KRR, which also has the largest test error of all models. On the other hand, LRR shows almost no overfit and its performance on the test set is much closer to that of the other models. Once again, ESNs and LSTMs outperform the stateless models, but this time the margin is slimmer and both models show a very similar performance. The values for multiple other performance metrics (Table 4) confirm this conclusion, with the minor differences being found in the mean absolute errors of the model, where the difference in performance between the recurrent models and the rest is more clearly pronounced, and where now the ESN has a slightly better, but still very comparable performance to the LSTM model.

Table 4: Regression performance metrics for the different models on the real-world dataset.

Fig. 9 shows the plots of the true and predicted pressure differences of the different models for some randomly selected cycles from the training and test sets of the real-world dataset. In this case, the outputs are much noisier and none of the models captures the dynamics perfectly. Once again all of the models capture the instantaneous dynamics fairly accurately, though not as well as in the synthetic dataset, and all of the models struggle with the non-linear degradation trend. The ESN and LSTM models capture the dynamics in the training set fairly accurately, but due to the differences in dynamics between the training and test sets, the accuracy of both stateful models is somewhat lower when predicting the degradation trend at the end of the cycles for the selected test cycles. The average predicted and true KPIs (Fig. 10) reveal that both LRR and LSTM seem to overestimate the trend at the end of the training cycles on average, while they underestimate it in the test cycles, which is likely a consequence of the covariate shift between the training and test sets of the real-world dataset. All other models also underestimate the exponential trend at the end of the test cycles on average, with the ESN predictions following the trend of the true KPIs most closely. This can point to better a generalization of the ESN in the case of longer sequences.

Overall, these results confirm that more complex stateful models still perform well in a small data regime, even in the presence of noise and covariate shifts.

Figure 9: Comparison of the true and predicted pressure differences across four ML models for some random training and test cycles from the real-world dataset.

Figure 10: Comparison of the means of the true and predicted pressure differences across four ML models for the training and test cycles of the synthetic dataset.

6. Discussion

Formulating accurate mathematical models of industrial aging processes (IAP) is essential for predicting when critical assets need to be replaced or restored. In world-scale chemical plants such predictions can be of great economic value, as they increase plant reliability and efficiency. While mechanistic models are useful for elucidating the influencing factors of degradation processes under laboratory conditions, it is notoriously difficult to adapt them to the spe-cific circumstances of individual plants. Data-driven machine learning methods, on the other hand, are able to learn a model and make predictions based on the historical data from a specific plant and are therefore capable of adapting effortlessly to a multitude of conditions, provided enough data is available. While simpler, especially linear prediction models have previously been studied in the context of predictive maintenance [30], a detailed examination of more recent and complex ML models, such as recurrent neural networks, was missing so far.

In this paper, we address the task of predicting a KPI, which indicates the slow degradation of critical equipment, over the time frame of an entire degradation cycle, based solely on the initial process conditions and how the process will be operated in this period. To this end, we have compared a total of five different prediction models: three stateless models, namely linear ridge regression (LRR), non-linear kernel ridge regression (KRR) and feed-forward neural networks (FFNN), and two recurrent neural network (RNN) based stateful models, echo state networks (ESN) and LSTMs. To assess the importance of the amount of available historical data on the models’ predictions, we have first tested them on a synthetic dataset, which contained essentially unlimited, noise-free data points. In a second step, we examined how well these results translate to real-world data from a large-scale chemical plant at BASF.

While the stateless models (LRR, KRR, and FFNN) accurately captured instantaneous changes in the KPIs resulting from changing process conditions, they were mostly unable to pick up on the underlying trend caused by the slower degradation effects. ESN and LSTMs, on the other hand, were able to additionally correctly predict long-term changes and obtained higher accuracy without requiring larger amounts of training data compared to the stateless models. With more parameters to tune, the stateful models can overfit on spe-cific patterns observed in the training data, however, the errors associated with overfitting are offset by the increased expressive power and the presence of temporal context, which enable these models to more accurately predict temporally dependent dynamics such as the long term degradation trend. In general, all models, especially those based on RNNs, yielded very promising predictions, which are accurate enough to improve scheduling decisions for maintenance events in production plants. The choice of the optimal model in a particular case depends on the amount of available data. For larger datasets, we found that LSTMs can yield almost perfect forecasts over long horizons. In the case of a smaller, noisy dataset, the stateful models still performed better than the simpler stateless models, disproving the common notion that complex models are always ‘data hungry’. Nevertheless, if only a small number of cycles are available for training (e.g. less than 100) or the data is very noisy, it can be advantageous to resort to a simpler regression model. In these cases, more extensive feature engineering [58] or inclusion of domain knowledge [20, 21, 22, 23, 24] could help to improve the prediction. Overall, ESNs prove to be a reasonable compromise, as they automatically expand the feature space and keep track of the internal hidden state using the randomly parametrized “reservoir” and therefore only require training to fit the set of output weights.

We believe that our work can serve as a basis for further research aimed at forecasting industrial aging processes with machine learning models. While the models examined in this manuscript exhibit good performance in the synthetic and real-world scenarios, there is still a wide range of research directions that can improve the performance.

While machine learning models are very good at interpolating, i.e., predicting values for data points from the same distribution as they were trained on, extrapolating beyond the patterns observed in the historical data is much harder and clearly a limitation of most machine learning models. Unfortunately, this is often what is required when dealing with real-world data. For example, on our real-world dataset, the machine learning models struggled to transfer their model of the training data to the test data, which came from a different catalyst charge and thus from a different distribution than the training data. These types of continuous changes and improvements to the setup of a production plant make it very difficult to compile a large enough dataset containing the relevant information for a model to predict future events. These time dependent changes in the data also require extra care as not to overfit on the past, e.g., when selecting hyperparameters for the models on the training set. Some of these effects might be tackled by explicitly addressing such non-stationarity between training and test set [59, 60]. It might furthermore be interesting to examine the effects of training on a larger collection of more diverse historical data and using transfer learning to be able to adapt the models to new datasets with smaller amounts of training data [61, 62]. This way, more expressive models such as LSTMs could be trained to learn some general patterns on a larger dataset and then fine-tuned on the more recent time points (which should be closer to the current degradation dynamics) to yield more accurate predictions for the future.

Especially when there is reason to believe that the predictions of a machine learning model may not always be perfectly on point, for example, when the model was trained on a very small dataset and predictions need to be made for inputs from a sightly different distribution, a predictive maintenance system could be further improved by providing confidence intervals, which quantify the uncertainty of the model for individual predictions [63]. This could, for example, be achieved by incorporating probabilistic approaches [64, 65, 66, 67]. Predicting intervals instead of point estimates may further facilitate planing by making it possible to assess the best and worst case scenarios.

While accurate predictions of IAPs will improve the production process by allowing for longer planing horizons, ensuring an economic and reliable operation of the plant, the ultimate goal is of course to gain a better understanding of and subsequently minimize the degradation effects themselves. While mechanistic and linear models are fairly straightforward to interpret, neural network models have long been shunned for their nontransparent predictions. However, this is changing thanks to novel interpretation techniques such as layer-wise relevance propagation (LRP) [68, 69, 70, 71, 72, 73], which make it possible to visualize the contributions of individual input dimensions to the final prediction. With such a method, the forecasts of RNNs such as LSTMs could be made more transparent, therefore shedding light on the influencing factors and production conditions contributing to the aging process under investigation, which could furthermore be used to help improve the underlying process engineering [74].

7. Acknowledgments

KRM acknowledges partial financial support by the German Ministry for Education and Research (BMBF) under Grants 01IS14013A-E, 01GQ1115 and 01GQ0850; Deutsche Forschungsgemeinschaft (DFG) under Grant Math+, EXC 2046/1, Project ID 390685689 and by the Technology Promotion (IITP) grant funded by the Korea government (No. 2017-0-00451, No. 2017-0-01779). MB acknowledges financial support by BASF under Project ID 10044628.

References

[1] P. Forzatti, L. Lietti, Catalyst deactivation, Catalysis today 52 (2-3) (1999) 165–181.

[2] J. Barbier, Deactivation of reforming catalysts by coking-a review, Applied Catalysis 23 (2) (1986) 225–243.

[3] P. Harris, Growth and structure of supported metal catalyst particles, In- ternational materials reviews 40 (3) (1995) 97–115.

[4] P. H. Nielsen, Poisoning of ammonia synthesis catalysts, in: Ammonia, Springer, 1995, pp. 191–198.

[5] H. Cai, A. Krzywicki, M. C. Oballa, Coke formation in steam crackers for ethylene production, Chemical Engineering and Processing: Process Intensification 41 (3) (2002) 199–214.

[6] O. Wu, A. E. Bouaswaig, S. M. Schneider, F. M. Leira, L. Imsland, M. Roth, Data-driven degradation model for batch processes: a case study on heat exchanger fouling, in: Computer Aided Chemical Engineering, Vol. 43, Elsevier, 2018, pp. 139–144.

[7] H. M¨uller-Steinhagen, Heat exchanger fouling: Mitigation and cleaning techniques, IChemE, 2000.

[8] B. Wang, Erosion-corrosion of thermal sprayed coatings in fbc boilers, Wear 199 (1) (1996) 24–32.

[9] J. Werther, Fluidized-bed reactors, Ullmann’s encyclopedia of industrial chemistry.

[13] M. Berreni, M. Wang, Modelling and dynamic optimization of thermal cracking of propane for ethylene manufacturing, Computers & Chemical Engineering 35 (12) (2011) 2876–2885.

[14] E. Ruckenstein, B. Pulvermacher, Growth kinetics and the size distribu- tions of supported metal crystallites, Journal of Catalysis 29 (2) (1973) 224–245.

[15] L. Li, P. N. Plessow, M. Rieger, S. Sauer, R. S. Sanchez-Carrera, A. Schae- fer, F. Abild-Pedersen, Modeling the migration of platinum nanoparticles on surfaces using a kinetic monte carlo approach, The Journal of Physical Chemistry C 121 (8) (2017) 4261–4269.

[16] G. Froment, Modeling of catalyst deactivation, Applied Catalysis A: Gen- eral 212 (1-2) (2001) 117–128.

[17] F. Brahim, W. Augustin, M. Bohnet, Numerical simulation of the fouling process, International Journal of Thermal Sciences 42 (3) (2003) 323–334.

[18] Z. Wang, M. Ierapetritou, Constrained optimization of black-box stochastic systems using a novel feasibility enhanced kriging-based method, Computers & Chemical Engineering 118 (2018) 210–223.

[19] S. H. Kim, F. Boukouvala, Surrogate-based optimization for mixed-integer nonlinear problems, Computers & Chemical Engineering (2020) 106847.

[20] M. Von Stosch, R. Oliveira, J. Peres, S. F. de Azevedo, Hybrid semi- parametric modeling in process systems engineering: Past, present and future, Computers & Chemical Engineering 60 (2014) 86–101.

[21] M. J. Willis, M. von Stosch, Simultaneous parameter identification and discrimination of the nonparametric structure of hybrid semi-parametric models, Computers & Chemical Engineering 104 (2017) 366–376.

[22] N. Asprion, R. B¨ottcher, R. Pack, M.-E. Stavrou, J. H¨oller, J. Schwientek, M. Bortz, Gray-box modeling for the optimization of chemical processes, Chemie Ingenieur Technik 91 (3) (2019) 305–313.

[23] S. Zendehboudi, N. Rezaei, A. Lohi, Applications of hybrid models in chem- ical, petroleum, and energy systems: A systematic review, Applied energy 228 (2018) 2539–2566.

[24] J. Glassey, M. von Stosch, Hybrid Modeling in Process Industries, CRC Press, 2018.

[25] E. Russell, L. Chiang, R. Braatz, Data-driven Methods for Fault Detec- tion and Diagnosis in Chemical Processes, Advances in Industrial Control, Springer London, 2012.

[26] L. H. Chiang, E. L. Russell, R. D. Braatz, Fault diagnosis in chemical pro- cesses using fisher discriminant analysis, discriminant partial least squares, and principal component analysis, Chemometrics and intelligent laboratory systems 50 (2) (2000) 243–252.

[27] J. V. Kresta, J. F. Macgregor, T. E. Marlin, Multivariate statistical moni- toring of process operating performance, The Canadian Journal of Chemical Engineering 69 (1) (1991) 35–47.

[28] V. Radhakrishnan, M. Ramasamy, H. Zabiri, V. Do Thanh, N. Tahir, H. Mukhtar, M. Hamdi, N. Ramli, Heat exchanger fouling model and preventive maintenance scheduling tool, Applied Thermal Engineering 27 (17-18) (2007) 2791–2802.

[29] J. Aminian, S. Shahhosseini, Evaluation of ann modeling for prediction of crude oil fouling behavior, Applied thermal engineering 28 (7) (2008) 668–674.

[30] J. H. Lee, J. Shin, M. J. Realff, Machine learning: Overview of the re- cent progresses and implications for the process systems engineering field, Computers & Chemical Engineering 114 (2018) 111–121.

[31] S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (1952) 89–94.

[32] D. Nemec, J. Levec, Flow through packed bed reactors: 1. single-phase flow, Chemical Engineering Science 60 (24) (2005) 6947 – 6957.

[33] D. P. Mandic, J. Chambers, Recurrent neural networks for prediction: learning algorithms, architectures and stability, John Wiley & Sons, Inc., 2001.

[34] V. N. Vapnik, The nature of statistical learning theory, Springer, 1995.

[35] K.-R. M¨uller, S. Mika, G. R¨atsch, K. Tsuda, B. Sch¨olkopf, An introduction to kernel-based learning algorithms, IEEE transactions on neural networks 12 (2).

[36] K. Hansen, G. Montavon, F. Biegler, S. Fazli, M. Rupp, M. Scheffler, O. A. Von Lilienfeld, A. Tkatchenko, K.-R. M¨uller, Assessment and validation of machine learning methods for predicting molecular atomization energies, Journal of Chemical Theory and Computation 9 (8) (2013) 3404–3419.

[37] M. Stone, Cross-validatory choice and assessment of statistical predictions, Journal of the royal statistical society. Series B (Methodological) (1974) 111–147.

[38] B. Sch¨olkopf, A. Smola, K.-R. M¨uller, Nonlinear component analysis as a kernel eigenvalue problem, Neural computation 10 (5) (1998) 1299–1319.

[39] B. Sch¨olkopf, A. J. Smola, Learning with kernels: support vector machines, regularization, optimization, and beyond, MIT press, 2002.

[40] K.-R. M¨uller, A. J. Smola, G. R¨atsch, B. Sch¨olkopf, J. Kohlmorgen, V. Vap- nik, Predicting time series with support vector machines, in: International Conference on Artificial Neural Networks, Springer, 1997, pp. 999–1004.

[41] C. E. Rasmussen, C. K. Williams, Regression, in: Gaussian processes for machine learning, Vol. 2, MIT press Cambridge, MA, 2006.

[42] A. J. Smola, B. Sch¨olkopf, Bayesian kernel methods, in: Advanced lectures on machine learning, Springer, 2003, pp. 65–117.

[43] C. M. Bishop, et al., Neural networks for pattern recognition, Oxford Uni- versity Press, 1995.

[44] I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016.

[45] J. Schmidhuber, Deep learning in neural networks: An overview, Neural networks 61 (2015) 85–117.

[46] R. Hecht-Nielsen, Theory of the backpropagation neural network, in: Neu- ral networks for perception, Elsevier, 1992, pp. 65–93.

[47] A. Choromanska, M. Henaff, M. Mathieu, G. B. Arous, Y. LeCun, The loss surfaces of multilayer networks, in: Artificial Intelligence and Statistics, 2015, pp. 192–204.

[48] Y. Bengio, N. Boulanger-Lewandowski, R. Pascanu, Advances in optimizing recurrent networks, in: 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, IEEE, 2013, pp. 8624–8628.

[49] H. Jaeger, H. Haas, Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication, Science 304 (5667) (2004) 78–80.

[50] F. M. Bianchi, E. Maiorino, M. C. Kampffmeyer, A. Rizzi, R. Jenssen, An overview and comparative analysis of recurrent neural networks for short term load forecasting, arXiv preprint arXiv:1705.04378 (2017) 78–80.

[51] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural computa- tion 9 (8) (1997) 1735–1780.

[52] A. Graves, N. Jaitly, Towards end-to-end speech recognition with recurrent neural networks, in: International Conference on Machine Learning, 2014, pp. 1764–1772.

[53] I. Sutskever, O. Vinyals, Q. V. Le, Sequence to sequence learning with neural networks, in: Advances in neural information processing systems, 2014, pp. 3104–3112.

[54] J. Donahue, L. Anne Hendricks, S. Guadarrama, M. Rohrbach, S. Venu- gopalan, K. Saenko, T. Darrell, Long-term recurrent convolutional networks for visual recognition and description, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 2625– 2634.

[55] A. Karpathy, The unreasonable effectiveness of recurrent neural networks, Andrej Karpathy blog 21.

[56] A. van den Oord, N. Kalchbrenner, K. Kavukcuoglu, Pixel recurrent neural networks, in: International Conference on Machine Learning, 2016, pp. 1747–1756.

[57] J. Friedman, T. Hastie, R. Tibshirani, The elements of statistical learning, Vol. 1, Springer series in statistics New York, NY, USA:, 2001.

[58] F. Horn, R. Pack, M. Rieger, The autofeat python library for automatic feature engineering and selection, in: P. Cellier, K. Driessens (Eds.), Machine Learning and Knowledge Discovery in Databases. ECML PKDD 2019, Springer International Publishing, Cham, 2020, pp. 111–120.

[59] M. Sugiyama, M. Krauledat, K.-R. M¨uller, Covariate shift adaptation by importance weighted cross validation, Journal of Machine Learning Research 8 (May) (2007) 985–1005.

[63] K. Efthymiou, N. Papakostas, D. Mourtzis, G. Chryssolouris, On a predic- tive maintenance platform for production systems, Procedia CIRP 3 (2012) 221–226.

[64] D. M. Blei, A. Kucukelbir, J. D. McAuliffe, Variational inference: A review for statisticians, Journal of the American statistical Association 112 (518) (2017) 859–877.

[65] S. Kolassa, Evaluating predictive count data distributions in retail sales forecasting, International Journal of Forecasting 32 (3) (2016) 788–803.

[66] K. Kasiviswanathan, K. Sudheer, Quantification of the predictive uncer- tainty of artificial neural network based river flow forecast models, Stochastic environmental research and risk assessment 27 (1) (2013) 137–146.

[67] M. Feindt, U. Kerzel, The neurobayes neural network package, Nuclear In- struments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 559 (1) (2006) 190–194.

[68] S. Bach, A. Binder, G. Montavon, F. Klauschen, K.-R. M¨uller, W. Samek, On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation, PLOS ONE 10 (7) (2015) e0130140.

[69] G. Montavon, W. Samek, K.-R. M¨uller, Methods for interpreting and un- derstanding deep neural networks, Digital Signal Processing 73 (2018) 1–15.

[70] G. Montavon, S. Lapuschkin, A. Binder, W. Samek, K.-R. M¨uller, Ex- plaining nonlinear classification decisions with deep taylor decomposition, Pattern Recognition 65 (2017) 211–222.

[71] L. Arras, G. Montavon, K.-R. M¨uller, W. Samek, Explaining recurrent neural network predictions in sentiment analysis, in: Proceedings of the 8th Workshop on Computational Approaches to Subjectivity, Sentiment and Social Media Analysis, 2017, pp. 159–168.

[72] P.-J. Kindermans, K. T. Sch¨utt, M. Alber, K.-R. M¨uller, D. Erhan, B. Kim, S. D¨ahne, Learning how to explain neural networks: Patternnet and patternattribution, in: International Conference on Learning Representations, 2018.

[73] S. Lapuschkin, S. W¨aldchen, A. Binder, G. Montavon, W. Samek, K.-R. M¨uller, Unmasking clever hans predictors and assessing what machines really learn, Nature communications 10 (1) (2019) 1096.

[74] Q. Zhao, T. Hastie, Causal interpretations of black-box models, Journal of Business & Economic Statistics (just-accepted) (2019) 1–19.

[75] N. R. Draper, H. Smith, Applied regression analysis, Vol. 326, John Wiley & Sons, 2014.

[76] C. M. Bishop, Pattern recognition and machine learning, Springer, 2006.

[77] Y. A. LeCun, L. Bottou, G. B. Orr, K.-R. M¨uller, Efficient backprop, in: Neural networks: Tricks of the trade, Springer, 2012, pp. 9–48.

[78] K. Hornik, M. Stinchcombe, H. White, Multilayer feedforward networks are universal approximators, Neural networks 2 (5) (1989) 359–366.

[79] K. Doya, Bifurcations in the learning of recurrent neural networks, in: Cir- cuits and Systems, 1992. ISCAS’92. Proceedings., 1992 IEEE International Symposium on, Vol. 6, IEEE, 1992, pp. 2777–2780.

[80] R. Pascanu, T. Mikolov, Y. Bengio, On the difficulty of training recurrent neural networks, in: International Conference on Machine Learning, 2013, pp. 1310–1318.

[81] H. T. Siegelmann, E. D. Sontag, On the computational power of neural nets, Journal of computer and system sciences 50 (1) (1995) 132–150.

[82] M. Lukoˇseviˇcius, A practical guide to applying echo state networks, in: Neural networks: Tricks of the trade, Springer, 2012, pp. 659–686.

Appendices

A. Mechanistic model to generate the synthetic dataset

The mechanistic process model generating the dataset has the following ingredients:

• Mass balance equations for all five relevant chemical species (olefinic reactant, oxygen, oxidized product, CO, water) in the reactor, which is for simplicity modeled as an isothermal plug flow reactor, assuming ideal gas law. The reaction network consists of the main reaction (olefine + Oproduct) and one side reaction (combustion of olefine to CO). With this, the mass flow of each species i = 1...5 at the reactor outlet is determined by the volumetric reaction rates , stoichiometric matrix , molar weights , reactor volume V , and residence time as follows:

Here, F = is the total mass flow through the reactor, which is conserved, and are the molar concentrations at the reactor inlet. With known reactor temperature T and pressure mass flow rates can be converted into mass fractions and molar concentrations :

• A highly non-linear deactivation law for catalyst activity A(t), which depends on the reaction temperature T (in Kelvin), flow rate F, inflowing oxygen , activity A(t) itself, and kinetic parameters and :

Catalyst activity is expressed in relative units (0% to 100%). As the deactivation can not be observed directly, A(t) is a hidden state variable of the system.

• Kinetic laws for the reaction rates , with kinetic parameters and :

• The relation for selectivity and conversion of the process:

Parameter values for are provided in Table A.1.

Table A.1: Parameter values for the generation of the synthetic dataset.

The cumulative feed of olefine for each cycle is calculated as:

where is the beginning of the current deactivation cycle, i.e., the last catalyst regeneration event with .

With the mechanistic model, the time series of the process dynamics are generated in the following way:

1. Beginning of new deactivation cycle: Set catalyst activity to A(t) = 100%.

2. Pick discrete random values for the process parameters , according to the procedure described below, and calculate the residence time , Eq. (A.2).

3. Set the duration of the next time window of constant process parameters [in hours] by a random integer number between 1 and 24.

4. For each hourly interval in this time window

(b) Calculate the mass flow at reactor outlet, Eq. (A.1), convert into molar concentrations, Eq. (A.4), and infer selectivity S and conversion C from it, Eq. (A.6).

(c) If C < 75%: End of cycle criterion met. Catalyst is regenerated. Begin new cycle (step 1).

(d) Reduce catalyst activity A(t) according to the deactivation dynamics, Eq. (A.5).

5. Go to step 2.

Generation of random process parameters. Each of the process parameters = , and is generated independently through the following random process:

1. At initial time t = 0, pick = 0) uniformly from the range of potential values [] for that parameter. This range is [3500, 4000] kg/h for F; p is in [1.25, 1.45] bar; T is in [500, 510] C; and in [47.5%, 52.5%].

2. At the time of the k-th change of process parameters, pick ) from a Gaussian distribution that is centered at the previous value ) and has a standard deviation of (10.

3. After generating the random values in this way, round them them to six equidistant values in []. This way, e.g., the potential values of mass flows F are (3500 + 100) kg/h, with n = 0, 1, ..., 5.

Finally, set = 1 . The remaining mass fractions at the reactor inlet are zero.

Stateless models are machine learning models that base their forecast only on the inputs within a fixed time window in the past, i.e., exactly as stated in Eq. (1).

Stateless models include most typical machine learning regression models, ranging from linear regression models to many types of neural networks [75, 76]. The stateless regression models explored in this paper are linear ridge regression (LRR), kernel ridge regression (KRR), and feed-forward neural networks (FFNN).

Linear ridge regression (LRR). The LRR model assumes that the process parameters x and KPIs y at a time t are linearly related. For this purpose, a weight matrix is used to model how each of the process parameters affects each KPI, which can be used to predict the KPIs in y(t) up to some noise ) that can not be explained by the model:

In order to reduce the influence of outliers on the model and to prevent over-fitting, we use a linear ridge regression model, which imposes a L2-regularization on the weight matrix of the standard linear regression model. The amount of regularization is controlled using the regularization parameter . Let denote the matrix containing the inputs of all cycles concatenated along the time dimension, while denotes the corresponding output matrix constructed analogously to X. The optimal weight matrix W can then be found by solving the following optimization problem:

The advantage of ridge regression is that the optimization problem can be solved analytically to obtain the globally optimal weight matrix W as

Despite its relative simplicity, LRR is widely used in many application scenarios and can often be used to approximate real-world processes at fairly high accuracies, especially if additional (non-linear) (hand-)engineered features are available [58]. Furthermore, considering the limited amount of training data that is usually available for real-world IAP problems, reliably estimating the parameters of more complex non-linear prediction models such as deep neural networks needs to be done with great care [77], while linear models provide a more robust solution as they provide a globally optimal solution and are less likely to overfit given their linear nature.

Kernel ridge regression. The linear ridge regression optimization problem in Eq. (B.1) is rewritten by applying the feature map on the inputs to get

leading to the analogous analytic solution

Since the feature space is not known explicitly, W can not be calculated directly. However, using some algebraic transformations and the kernel trick [38, 35, 39], the following solution can be derived:

The kernel function used in this paper is the radial basis function (RBF) kernel, also known as the Gaussian kernel:

The kernel width parameter , along with the regularization parameter results in a total of two hyperparameters to optimize for the KRR model.

Feed-forward neural networks (FFNN). Let denote the feature vector obtained after the transformation at the l-th layer, and let be the linear transformation matrix, the bias and the non-linear function applied at layer l. Then a generic l-layer feed-forward neural network can be defined by the following sequence of operations:

An illustration of this typical FFNN architecture is also shown in Fig. B.1.

Figure B.1: An overview of a basic FFNN architecture with multiple hidden layers with non-linearities.

The 1st layer is usually called the input layer and the l-th layer the output layer, while all the layers in-between are known as hidden layers. By changing the values of the weight matrices of the different layers, the network can be adjusted to approximate the target function. In fact, multi-layer FFNN are proven to be universal function approximators, meaning that the weights can be adjusted to fit any continuous function [78]. The most common method for training neural networks to fit a given function is to use error backpropagation to iteratively update the values of the parameters of the network, usually consisting of the weight matrices and the bias vectors [77]. Error backpropagation is a way to efficiently perform gradient descent on the FFNN parameters by minimizing the error/loss of the network’s outputs with respect to some target labels. One of the most commonly used loss functions for regression problems, and also the one we use in this paper, is the squared error (SE), which for a given FFNN, denoted by , where denotes the network’s parameters, an input time series x, and corresponding target labels y is given by:

By taking the gradient of the loss function with respect to the FFNN’s parameters, using backpropagation, one can efficiently perform gradient descent to update the parameters and minimize the loss of the network [43]. One apparent weakness of gradient descent it that it often converges at local minima, in contrast to LRR and KRR, where the globally optimal solution to the optimization problem can be obtained analytically. However, in neural networks the losses in these local minima are often similar to the global optimum [47], so this properties does not significantly impact the performance of a properly trained neural network. Additionally, due to a FFNN’s large number of parameters () and high flexibility, if not properly trained it may overfit, especially when using smaller training sets.

B.2. Stateful models

In contrast to stateless models, stateful models only explicitly use the input x(t), not the past inputs 1)), to forecast the output y(t) for some time point t. Instead, they maintain a hidden state h(t) of the system that is continuously updated with each new time step and thus contains information about the entire past of the time series. The output can then be predicted utilizing both the current input conditions, as well as the hidden state of the model: ˆy(t) = f(x(t); h(t)).

The two stateful models that we are considering for this paper both belong to the class of recurrent neural networks (RNNs). RNNs are a powerful method for modeling time series, however they can be difficult to train since their depth increases with the length of the time series. If training is not performed carefully, this can lead to bifurcations of the gradient during the error backpropagation training procedure, which can result in a very slow convergence (“vanishing gradients problem”), if the optimization converges at all [79, 80].

The most commonly used stateful models for the modeling of sequential data are recurrent neural networks (RNNs) [33]. While RNNs are some of the most powerful neural networks, capable of approximating any function or algorithm [81], they are also more involved to train [79, 80]. Consequently, in this paper we chose to model IAPs using two different RNN architectures that are designed precisely to deal with the problems arising while training regular RNNs: echo state networks (ESN) and long short term memory (LSTM) networks.

Figure B.2: An overview of the ESN architecture.

Echo state networks (ESN). The basic architecture of an ESN is displayed in Fig. B.2. The ESN is composed of two randomly initialized weight matrices, the input weight matrix and the (usually sparse) recurrent weight matrix , where m is the dimensionality of the hidden state vector and is usually much larger than the dimensionality of the input features . We will refer to these two matrices collectively as the reservoir. Let be the hidden state vector at time t and ˜h(t) its update, while (0, 1] is called the leaking rate determining trade-off between adding new and keeping old information in the hidden state. Additionally, let [] denote vertical vector/matrix concatenation. Now given an input vector x(t), the update for the hidden state is given by

Note that tanh is applied element-wise here. It may be counter-intuitive that the randomly generated matrices and can produce useful features without being trained, but because m >> d the transformation in Eq. (B.2) can be interpreted as a random feature map that produces a non-linear, high dimensional feature expansion with memory of the previous inputs. In this sense, one can draw a parallel to kernel methods, which also map the input features to a non-linear and high dimensional features space, where the features are more easily linearly separable [82].

Finally, the predictions are made based on the time-dependent features generated by Eq. (B.2) using a final output matrix as

where is trained to minimize the MSE of ˆy(t), in our case using LRR as described in Section B.1. Using LRR provides a globally optimal analytical solution for , making training the ESN simple and avoiding the problems that arise when training RNNs with error backpropagation. Due to the high dimensionality and non-linearity of the hidden state features, a simple model like LRR can still fit complex and non-linear dependencies between the input and the output features. Additionally, since LRR is very robust against overfitting, ESNs are also very well suited for problems with smaller amounts of data.

One downside to ESN is that their performance is strongly dependent on the choice of hyperparameters, which mostly regulate the initialization of the reservoir. Since the reservoir weights are not trained, this initialization is crucial in order to ensure that the ESN has desirable properties. The main hyperparameters of the ESN include the dimensionality of the hidden state m, the sparsity of the matrix , the distribution from which the non-zero elements of the reservoir are sampled, the spectral radii of each of the reservoir matrices, the scaling of the input data and finally the leakage rate . These hyperparameters, as well as general recommendations for their choice are summarized in Table B.1 (for more details see [82]).

One of the most important hyperparameters is the spectral radius of the matrix, because keeping this close to 1 helps maintain the echo state property of the network, which is essential for the ESN to be well-behaved. In a nutshell, the echo state property means that the hidden state h(t) should be uniquely defined by the fading history of the input x, i.e., for a long enough

Table B.1: Summary of the main ESN hyperparameters alongside general recommendations for choosing their values.

Connectivity of A small fixed number of non-zero elements (e.g. 10) per row on average, irrespective of the hidden state size.

Spectral radius of Larger for highly non-linear and small for linear tasks

Leakage rate Adjust according to the dynamics of the time signal; if signal is changing slowly set close to zero, for fast dynamics set close to 1

input x(t), the hidden state h(t) should not depend on the initialization of the reservoir anymore.

Figure B.3: An overview of the LSTM architecture.

LSTM networks. The cell state is a unique and core component of the LSTM and runs through the entire recurrent chain. It can only be updated slowly at each time step using only linear updates, making it capable of preserving long term dependencies in the data and maintaining a stable gradient over long sequences. The inclusion of new or removal of old information to the cell state is carefully regulated by special neural network layers called gates. Fig. B.3 shows the full architecture of one recurrent LSTM cell. The first gate of the LSTM is the forget gate, which is trained to regulate which information is to be removed or ‘forgotten’ from the cell state. The forget gate produces a vector f(t) that is calculated by applying the sigmoid function, denoted by , element-wise to a linear transformation of the current input x(t) and the previous hidden state 1), resulting in the equation

The non-linearity ensures that the values of f(t) are between 0 and 1, with a value of 0 of ) meaning forget everything at position i in the cell state, and a value of 1 meaning preserve the information at this position completely. Analogously to the forget gate, the input gate creates a vector i(t), which regulates which values of the cell state will be updated. Additionally, another layer is used to create a vector ˜c(t) of candidate update values that could be added to the cell state based on i(t). This results in the following two equations:

where the sigmoid non-linearity is once again used to obtain values between 0 and 1 for i(t), while the candidate update is generated using a tanh non-linearity to ensure that the values in the cell state can be updated both by increasing and decreasing them. In the next step, the cell state is updated based on the vectors produced by the forget and input gates using element-wise multiplication, denoted by . This update is given by

Then, the next short-term hidden state needs to be generated, which is based on the now updated cell state. A sigmoid non-linearity is once again used to regulate to which degree the values of the cell state will be included in the hidden state, and a tanh non-linearity is used to compress the cell state in the interval (-1, 1). The new hidden state is thus calculated as

Finally, the output vector is calculated from the hidden state. This part is not always required, as many time series problems require only one output for a whole sequence. In our case however, we need an output for each time point, which is then calculated using a simple linear layer as

The interplay of the different gates contributes to the stability of the gradients, making the LSTM very well suited for time series problems with long-term dependencies. This is why LSTMs have become one of the most widely used machine learning models for sequential data, such as speech-to-text recognition [52], machine translation [53], video classification [54], text [55] and image generation [56].

C. Model hyperparameters

All of the models used in this paper have certain hyperparameters that need to be optimized and the correct choice of these hyperparameters can in some cases strongly influence the performance of the model. Of the models used in this paper, the performance of the KRR and ESN models was strongly influenced by the choice of hyperparameters, while the performance of the other models remained stable over a wider range of hyperparameter values. For the ESN model, which had the largest number of hyperparameters, we used random search to select sets of hyperparameters for validation, while for the other models the hyperparameters were optimized using a grid search. In the following we report the hyperparameters that resulted in the lowest validation errors for each of our models.

LRR only has one hyperparameter, the regularization parameter , which was selected as 0.01 for the synthetic dataset and 0.5 for the real-world dataset (i.e., the noisier real-world dataset required a stronger protection against outliers compared to the relatively clean synthetic dataset). KRR has two hyperparameters, the regularization parameter and the kernel width . The resulting optimal hyperparameters from our cross-validation procedure were = 0.00046 and = 121.99 for the synthetic dataset and = 0.00183 and = 1.75 for the real-world dataset. Since the other models have multiple hyperparameters that need to be optimized, the optimal values for each dataset are presented in Table C.1 for the FFNN, Table C.2 for the ESN model, and Table C.3 for the LSTM model hyperparameters.

Table C.1: The optimized hyperparameters for the FFNN model for both datasets.

Table C.2: The optimized hyperparameters for the ESN model for both datasets.

Hidden state size 343 58 Connectivity of 32 90 Distribution non-zero values Standard normal Standard normal Spectral radius of 1.18 0.77 Spectral radius of 0.004 0.0016 Bias scale for 0.68 0.08 Leakage rate 0.05 0.012

Table C.3: The optimized hyperparameters for the LSTM model for both datasets.

designed for accessibility and to further open science