Transfer Function Modeling of Processes with Dynamic Inputs

Time series structures, which are common occurrences with data in many industrial processes, complicate a quality practitioner's efforts to accurately position control chart limits. ARIMA modeling and a variety of control charting methods have been recommended for monitoring process data with a time series structure. Estimates of ARIMA model parameters may not be reliable, however, if assignable causes of variation are present in the process data used to fit the time series model. Control limits may also be misplaced if the process inputs are dynamic and exhibiting a time series structure. Our purpose in this paper is to explore the ability of a transfer function model to identify assignable causes of variation and to model dynamic relationships between process inputs and outputs. A transfer function model is developed to monitor biochemical oxygen demand output from a wastewater treatment process, a process with dynamic inputs. This model is used to identify periods of disturbance to the wastewater process and to capture the relationship between the variable nature of the input to the process and the resulting output. Simulation results are included in this study to measure the sensitivity of transfer function models and to highlight conditions where transfer function modeling is critical.


Introduction
li N an extensive survey, Alwan and Roberts (1995) found that more than 85% of industrial process control applications resulted in charts with possibly misplaced control limits. In many instances, the misplaced control limits result from autocorrelation of the process observations, which violates a basic assumption often associated with the Shewhart chart (Woodall (2000)). Autocorrelation of process observations has been reported in many industries, including cast steel (Alwan (1992)), blast furnace operations (Notohardjono and Ermer (1986)), wastewater treatment plants (Berthouex, Hunter, and Pallesen Dr. West is an Assistant Professor in the Department of Decision Sciences. He is a Member of ASQ. His email address is westd@mail.ecu.edu Dr. Dellana is an Associate Professor in the Department of Decision Sciences. He is a Member of ASQ. Dr. Jarrett is a Professor in the Department of Management Science and Statistics. (1978)), chemical processes industries (Montgomery and Mastrangelo (1991)), semiconductor manufacturing (Kim and May (1994)), injection molding (Smith (1993)), and basic rolling operations (Xia, Rao, Shan, and Shu (1994)).
Several models have been proposed to monitor processes with autocorrelated observations. Alwan and Roberts (1988) suggest using an autoregressive integrated moving average (ARIMA) residuals chart, which they referred to as a special cause chart. For subsample control applications, Alwan and Radson (1992) describe a fixed limit control chart, where the original observations are plotted with control limit distances determined by the variance of the subsample mean series. Montgomery and Mastrangelo (1991) use an adaptive exponentially weighted moving average (EWMA) centerline approach, where the control limits are adaptive in nature and determined by a smoothed estimate of process variability. Lu and Reynolds (2001) investigate the steady state average run length of cumulative sum (CUSUM), EWMA, and Shewhart control charts for autocorrelated data Vol. 34, No.3, July 2002 316 DAVID WEST, SCOTT DELLANA AND JEFFREY JARRETT modeled as a first order autoregressive process plus an additional random error term.

Transfer Function Modeling of Process Data
If a process quality characteristic, Zt, has a time series structure, an ARIMA model of the following general form can represent the undisturbed or natural process variation: We propose a more general transfer function: an ARIMA model that accounts for both outliers in process output and dynamic effects from process input. In the following section, we briefly describe the relevant theory of time series analysis used in this paper. We then analyze the transfer function model terms to identify disturbances in a wastewater treatment process. We follow in later sections with supporting empirical evidence on the sensitivity of these methods. The paper concludes with a discussion of the implications for quality practitioners who may be monitoring processes which produce data with time series structures and which have dynamic inputs. In this paper, we refer to autocorrelated processes as either autoregressive (AR) or as moving average (MA) (as defined by Box, Jenkins, and Reinsel (1994)). (1)

¢(B)a(B)Zt = B(B)at.
In Equation (1) Finally, at is a white noise series with distribution N(O,(]"~). This model is described by Chen and Liu (1993b). If the series Zt is contaminated by periods of external disturbances to the process, the ARIMA model may be incorrectly specified, the variability of the residuals overestimated, and the resulting control limits incorrectly placed.
The following transfer function model of Box and Tiao (1973) describes the observed quality characteristic, Yt, as a function of three sources of variability: The ARIMA and intervention models are appropriate for autocorrelated processes whose input streams are closely controlled. However, there are quality applications, which we refer to as "dynamic input processes," where this is not a valid assumption. The treatment of wastewater is one example of a dynamic process that must accommodate highly fluctuating input conditions. In the health care sector, the modeling of emergency room service must also deal with highly variable inputs. The dynamic nature of the input creates an additional source of variability in the system, namely the time series structure of the process input. For these applications, modeling the dynamic relationship between process inputs and outputs can be used to obtain im-A problem with all of these control models is that the estimate of the process variance is sensitive to outliers. If assignable causes are present in the data used to fit the model, the model may be incorrectly identified and the estimators of model parameters may be biased, resulting in loose or invalid control limits (Boyles (2000)). To justify the use of these methods, researchers have made the assumption that a period of "clean data" exists to estimate control limits. Therefore, methods are needed to assure that parameter estimates are free of contamination from assignable causes of variation. Intervention analysis, with an iterative identification of outliers, has been proposed for this purpose. The reader interested in more detail should see Alwan (2000, pp. 301-307), Atienza, Tang, and Ang (1998), and Box, Jenkins, and Reinsel (1994, pp. 473-474). Atienza, Tang, and Ang (1998) recommend the use of a control procedure based on an intervention test statistic, -X, and show that their procedure is more sensitive than ARIMA residual charts for process applications with high levels of positive autocorrelation. They limit their investigation of intervention analysis, however, to the detection of a single level disturbance in a process with high levels of first order autocorrelation. Wright, Booth, and Hu (2001) propose a joint estimation method capable of detecting outliers in an autocorrelated process where the data available is limited to as few as 9 to 25 process observations. Since intervention analysis is crucial to model identification and estimation, we investigate varying levels of autocorrelation, autoregressive and moving average processes, different types of disturbances, and multiple process disturbances.

317
where Wo is a constant. A step disturbance to the process is modeled as a level-shift outlier (a form of innovational outlier or 10) in the form, Chang, Tiao, and Chen (1988) and Chen and Liu (1993b) discuss both types of disturbance. Chang, Tiao, and Chen (1988) extended the concepts of Box and Tiao (1975) to an iterative method for detecting the location and nature of outliers at unknown points in the time series. They define a procedure for detecting innovational outliers and additive outliers, and for jointly estimating time series parameters. Their work also demonstrates that the The first term, v(B)Xt-b' is the dynamic input term and represents an impulse function, v(B), applied to the input Xt-b with a lag of b time periods. If a dynamic relationship between the input and output time series exists, lagged values of process inputs can be modeled, resulting in considerable reduction of unexplained variance. The second term, (w(B)j5(B))It, is the intervention term and identifies periods of time when assignable causes are present in the process. Here, It is an indicator variable with a value of zero when the process is undisturbed and a value of one when a disturbance is present in the process. See, for example, Box, Jenkins, and Reinsel (1994, p. 392) for the development of the transfer function term, and Box, Jenkins, and Reinsel (1994, p. 462) for details of the intervention term. The rational coefficient term of It is a ratio of polynomials that defines the nature of the disturbance as detailed in Box, Jenkins, and Reinsel (1994, p. 464). The third term, (B(B)/¢(B))at, is the basic ARIMA model of the undisturbed process from Equation (1). We refer to Equation (2) as the "transfer function" model throughout this paper.
Different types of disturbances can be modeled by the proper design of the intervention term. The two most common disturbances for quality applications are a point disturbance, with an impact observed for only a single time period, and a step disturbance, with an impact persisting undiminished through several subsequent observations. The point disturbance is modeled as an additive outlier (AO). An AO impacts the observed process at one observation. The AO is modeled in the form

Application of Transfer Function Modeling to Wastewater Treatment
Wastewater contains a variety of substances that must be removed, such as human wastes, food scraps, oils, soaps, chemicals, microorganisms, phosphorous compounds, nitrogen compounds, suspended solids, and organic matter. A critical step in the typical wastewater treatment process, shown in Figure 1, is the use of microorganisms to decompose the organic matter. The amount of oxygen consumed by the microorganisms, called the biochemical oxygen demand (BOD), is used as a process control measure. The challenges of controlling the BOD levels of a wastewater treatment plant result from the complexity of the process and the variances in the composition and flow rate of the input stream as discussed by Wen and Vassiliadis (1998). The value of transfer function modeling is illustrated in the following example, where we analyze 527 daily measurements of BOD from an urban wastewater treatment plant as reported by Poch, Bejar, and Cortes (1993). The dynamic nature of the input and output BOD is evident from the first 100 daily measurements plotted in Figures 2 and 3.
In the following subsections, we sequentially analyze the three sources of variability in the transfer function model of Equation (2). We first evaluate an ARIMA residual chart to serve as a benchmark. The second subsection details an analysis that includes intervention terms in the model, and the final subsection includes the dynamic input term. In all cases, model parameters are estimated from the presence of outliers may cause serious bias in the estimation of ARIMA model parameters. The iterative method of Chang, Tiao, and Chen (1988) is effective for large, isolated outliers, but a masking effect may occur when multiple outliers cluster together. To overcome these problems, Chen and Liu (1993a) proposed an iterative outlier detection and adjustment methodology designed to handle multiple outliers of various types in a time series. Box, Jenkins, and Reinsel (1994, p. 473) report the significance of the intervention term in the estimation of (Ya for two chemical process applications. For a chemical process temperature reading, three disturbances were identified by the intervention terms, resulting in a 26% reduction in &a. For chemical process viscosity readings, a single disturbance was identified by the intervention term and modeled with a resulting reduction in &a of 6 percent.  I  I  I  I  I , first 100 days of operations and the remaining 427 days serve as an independent holdout sample. The ARIMA analysis and estimation was performed using SAS software; the intervention test statistic, .x, was calculated with an Excel spreadsheet using a method similar to the approach of Atienza, Tang, and Ang (1998).    (1)). 2. From the model estimated in step 1, compute the residuals, et, and Ua' 3. Estimate A1T and A2T to determine the presence of AO and 10 outliers, respectively, for t = 1 to n observations. The first statistic, A1T' is used to test the hypothesis that H o: Wo = 0, H 1: Wo =I=-0 in Equation (3). The second statistic, A2T' is used to test H o : Wo = 0, H 2 : Wo =I=-0 in Equation (4). These estimators are defined as i Jenkins, and Reinsel (1994, pp. 183-223) and Montgomery and Johnson (1976, pp. 193-200) discuss such an analysis. The autocorrelation function for the first 100 observations cuts off abruptly at a lag of one, while the partial autocorrelation function follows a sinusoidal decay pattern. This information suggests that the process is stationary, and that an MA(I) model is appropriate for the BOD output of the wastewater treatment plant. Estimation of MA(I) model parameters yields the following result:

Dumped in
If the fitted model is an adequate representation of the sample data, then the model residuals shown in Figure 4 should approximate a white noise process, an important property for the success of the residuals control chart. The obvious disturbance beginning in day 60 can affect model identification and parameter estimates. The magnitude of this impact will be clear in the next subsection, where we identify and model these events. For the moment, we will ignore the disturbance and estimate (Jain Equation (5). The results, shown in Table 1, are that u a = 28.1 and that the 3-sigma limits are ±84.3. The residuals control chart identifies days 60 and 61 of Figure 4 as potential process disturbances. Applying these control limits to holdout observations for days 101-527 identifies a single disturbance at day 465.

Addition of Intervention Analysis to the Wastewater Process Model
The use of the iterative method of Chang, Tiao, and Chen (1988) can improve on the ARIMA residuals results by identifying and explicitly modeling wastewater process disturbances. This method is based on the likelihood ratio criteria, A, defined by Fox (1972), and uses the following algorithm: The seven outliers detected by intervention analysis did not result in the identification of a different form of model; in this case the model remains an MA(1) model. However, the value of~is reduced For the wastewater treatment process, the result from step 1 of the iterative method is the model estimated in Equation (5). Applying steps 2-6 with C = 3 yields the following ARIMA intervention model: In Equations (8) and (9), the 7fi terms are the coefficients of the following rational polynomial with terms defined in Equation (1)

Sensitivity of Intervention Analysis
The quality practitioner may be concerned with the power of intervention analysis to detect and model process disturbances and in understanding those process conditions where it is most important to use transfer models. We explore these issues next using simulated data.
from 0.582 (estimated for the basic ARIMA residuals chart) to 0.337. Most importantly, the estimate (Ta, which is the basis for placing control limits, is reduced. The estimate (Ta = 28.1 of the basic ARIMA residuals is reduced significantly to aa = 6.24 with the inclusion of intervention terms (see Table 1). This results in tighter placement of control limits for a residuals chart and more sensitivity to identify additional potential disturbances to the wastewater process. Potential disturbances are now identified for days 31, 60-64, 90, and 95 and for days 465-467 of the holdout sample.

Impact of Process Disturbances on Variance Estimation
We next include a term for the dynamic regression of the input BOD in the control model. Autocorrelation and partial autocorrelation functions calculated for the input BOD suggest that this time series is also an MA(I) process. We fit a transfer function model of the form of Equation (2) yielding the following result:

Addition of Dynamic Input Analysis to the Wastewater Process Model
where Xt is the BOD input to the process at time t and It represents the same seven intervention terms defined in Equation (10). By modeling the variability of the BOD input to the wastewater process, the estimate (Ta is reduced to 5.84 (see Table 1). The placement of the control limits is now at ±17.5, the tightest control limits of all three models investigated. Including the dynamic input term in the model provides increased sensitivity, identifying a potential disturbance at day 114, which use of the model in Equation (10) did not detect.
To understand the impact of outliers on the estimation of the noise variance in ARIMA models, we estimate aa for simulated time series ranging in length from 50 to 200 observations for 8 1 and ¢1 values of 0.3, 0.5, and 0.7. The impact of a single point 7]AO = max 1).2,t I 7]10 = max l).l,tl t = 1, ... ,no If either or both of these values exceeds a threshold, C, which typically varies from 2.5 to 4, it is determined that an outlier exists. If the maximum value is 7]10, then the outlier is an 10. Conversely, if the maximum is 7]AO, then the outlier is an AO. The residuals are adjusted by a method that is dependent on the nature of the outlier. For an 10 outlier, the model residual for that single period is set to zero. For an AO outlier the model residuals are adjusted using After adjusting the model residuals, a new value of (Ta is computed.
5. Continue to identify additional outliers by repeating steps 3 and 4 with modified residuals and re-estimating aa until no further outlier candidates exist.
6. Treat the disturbance times for all outliers identified above as known and simultaneously estimate the ARIMA model parameters and interventions. Repeat steps 1 to 6 with this model and stop when no additional outliers are identified.  Table 2 for AR(1) models and Table 3 for MA(1) models are averages of 1,000 repetitions. Since the simulated data at is constructed with aa = 1, the 3sigma control limits should be at ±3. To the extent that the model estimate of aa exceeds 1, the control limits are misplaced due to the failure to identify and model the process disturbance. The results suggest that for time series of length greater than 200 observations, a single disturbance of 4a a or less will have a modest impact of 4% or less on the placement of control limits. For shorter time series of approximately 50 observations, disturbances of 2aa or greater will significantly bias the estimators of control limits. These limits will be overestimated from Vol. 34, No.3, July 2002 4% to 16%. These conclusions hold for both AR(1) and IVIA(1) time series structures and for all levels of e 1 and 1;1 investigated. The effect of two disturbances in the time series is summarized in Table 4 for an AR( 1) time series with 1;1 = 0.5. In this case, the shorter time series control chart limits will be overestimated by 9% to 30%, and, for the longer series of 200 observations, overestimated by 2% to 9%.
Power of Iterative Intervention Analysis Chang, Tiao, and Chen (1988) report the power of their iterative method to detect large outliers varying from 3a a to 5a a using simulated AR(1) (1;1 = 0.6) and rVIA(1) (e 1 = 0.6) time series consisting of 150 observations. They found that the iterative method achieves a power of 0.975 for identification of an 10, and a power of 0.99 for identification of an AO in both AR(l) and MA(l) time series with SeTa disturbances. For the smaller 3eT a disturbance, the power is reduced to 0.55 for an 10 and 0.69 for an AO for both time series structures.
We extend the work of Chang, Tiao, and Chen (1988) by investigating the power of the iterative detection method for outliers of leTa' 2eT a, and 3eT a in both AR(l) and MA(l) time series to characterize more typical process-related outliers. The experi-  Step, Point 0.0, 1.0, 2.0, 3.0 (multiples of u a ) sented by either ¢l for the autoregressive model or B 1 for the moving average model), the nature of the assignable cause of variation, and the magnitude of the disturbance. The parameters of the time series structure include a low level of 0.3, a moderate level of 0.5, and a high level of 0.7. Four different levels of disturbances are studied: O.Ou a , 1.0u a , 2.Ou a , and 3.0u a . We investigate both point disturbances and step disturbances. The disturbance is generated at random locations between observation 2 and observation 150, which is the end of the simulated data.
Reported results are averages of 1000 repetitions for each condition. The two aspects of the power of intervention models in process control applications are measured in this study by the rate of true positives and the rate of false positives (or false alarms).  The experimental results of the intervention analysis are reported in Table 6 for the AR(1) models and in Table 7 for the MA(l) models. These tables summarize the proportion of times that a disturbance was detected by the intervention analysis for each experimental condition. Our operational definition for detecting a disturbance is based on using a threshold value of C = 3. We also report operating characteristic curves for other values of the threshold C in Figure 5.
The results of the AR(l) time series of Table 6 indicate that a point disturbance is easier to identify than a step disturbance. This is likely caused by the fact that the persistent nature of the step disturbance Journal of Quality Technology results in significantly overestimating aa and thereby understates the values of >. calculated in Equations (6) and (7). The results also suggest that it is easier to detect a disturbance at higher values of autocorrelation. The 3a a point disturbance is detected 50.5% of the time at (Pt = 0.3, 61.8% at <PI = 0.5, and 70.2% at (Pt = 0.7. The corresponding results for the step disturbance are 20.1% for <PI = 0.3, 26.5% for <PI = 0.5, and 30.9% for <PI = 0.7. Detection of disturbances smaller than 3a a is unreliable. From 1% to 23% of disturbances of size 1aa and 2aa are detected. From the results of Table 7, it is evident that outlier detection is more difficult with the MA(l) than The need to analyze quality measurements for a time series structure is well documented in the quality literature. Perhaps less well understood are the benefits of checking process inputs for time series structure. If process inputs vary over time, they produce an additional source of variability that should be modeled. When analyzing a wastewater treatment plant, for example, we found that the input BOD had an MA(l) structure. By incorporating this information into the transfer function model for wastewater treatment, we were able to reduce our estimate of aa by 6.4%.
A disadvantage of the application of time series methods to process control is the loss of simplicity ity to identify process disturbances resulting from assignable sources of variation prior to the estimation of model parameters, and its ability to explicitly model relationships between dynamic process inputs and output quality levels. This allows the quality practitioner to more accurately estimate the process variability and minimize the problem of misplaced residual chart control limits. The use of a transfer function model is most beneficial for relatively short time series with 50 to 200 observations. For these short time series, control limits can be overestimated by as much as 15% from a single reasonably small point disturbance and overestimated by 30% for two such disturbances. Small changes in control limits can have large effects on the statistical performance of the chart. For time series with more than 200 observations, transfer function modeling is most important under conditions of large (greater than 4a a ) or multiple disturbances.
Our simulation studies indicate that point disturbances, confined to a single observation, are more easily detected than are step disturbances, with an effect persistent over many observations. It also appears to be easier to detect a given disturbance in an AR(l) time series than in an MA(l) time series, and detection rates increase for AR(l) models at higher values of the autoregressive parameter cPl. The sensitivity of the transfer function model to detect outliers is dependent on the choice of the threshold, C, used in the iterative procedure. Operating characteristic curves constructed from the simulated results suggest that a threshold value of C = 2.5 provides a good tradeoff, resulting in high sensitivity while maintaining a false positive rate of less than 1%. For a threshold of 2.5, approximately 79% of 3a a disturbances are detected in AR(l) series and 66% in MA(l) series.

Implications for Quality Practice
Transfer function modeling is shown to be useful in monitoring process quality because of its abil-The false positive rate can be inferred from Tables 6 and 7 by the proportion of disturbances detected for the experimental conditions with a zero disturbance size. Since there is no disturbance for this condition, the distinction between the step and point disturbance is irrelevant. By averaging across all three parameter levels, we estimate that the intervention model detected a "false disturbance" about 0.37% of the time. A similar analysis for the MA(l) results yields a false positive rate of 0.38%.
The selection of an appropriate threshold, C, to use in identifying potential outliers is dependent on the economics of the process. To guide this decision we plot operating characteristic curves in Figure 5 that define the tradeoff between true positive and false positive rates for the detection of a single disturbance in AR(l) and MA(l) time series with threshold values of C = [2,2.5,3,3.5,4]. Higher values of C correspond to the lower left of the figure and lower values to the upper right. For both AR(l) and MA(l) time series, a significant increase in sensitivity occurs by lowering the threshold from C = 4 to C = 2.5. Below C = 2.5, the gains in sensitivity are at the expense of a substantially higher number of false positives. At C = 2.5, approximately 79% of point disturbances are identified in AR(l) time series and 66% in MA(l) series. The corresponding false positive rate is slightly less than 1%.
which is characteristic of Shewhart charts. We would like to emphasize that the computational burden of transfer function modeling is not as formidable as it may at first appear. The determination of ARIMA models with intervention terms and dynamic relationships is only performed periodically, when the identification of a new model is appropriate. Commercial software products are available that have automated routines for this analysis. Weller (1994) and Kusters (1995) identify PC-Expert from Scientific Computing Associates (www.scausa.com) and Autobox 3.0 (www.autobox.com) as two of the leading commercial products with capabilities for intervention analysis and transfer function modeling. The daily monitoring of the process can be accomplished by transposing the transfer function model into a spreadsheet format. For example, if Equation (1l) were the basis of wastewater process monitoring, the operator would be required to enter only an input and output BOD rate. With this information, an appropriate monitoring strategy can be employed. For example, a residual can be determined for the residual chart of Alwan and Roberts (1988) or an intervention statistic calculated for the model of Atienza, Tang, and Ang (1998).