434 IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 13, NO. 3, MAY 2005
Finding the Source of Nonlinearity in a Process With Plant-Wide Oscillation
Nina F. Thornhill
Abstract—A plant-wide oscillation in a chemical process often
has an impact on product quality and running costs and there is,
thus, a motivation for automated diagnosis of the source of such a
disturbance. This brief describes a method of analyzing data from
routine operation to locate the root cause oscillation in a dynamic
system of interacting control loops and to distinguish it from prop-
agated secondary oscillations. The novel concept is the application
of a nonlinearity index that is strongest at the source. The index is
large for the nonsinusoidal oscillating time trends that are typical
of the output of a control loop with a limit cycle caused by nonlin-
earity. It is sensitive to limit cycles caused both by equipment and
by process nonlinearity. The performance of the index is studied
in detail and default settings for the parameters in the algorithm
are derived so that it can be applied in a large scale setting such
as a refinery or petrochemical plant. Issues arising from artifacts
in the nonlinearity test when applied to strongly cyclic data have
been addressed to provide a robust, reliable and practical method.
The technique is demonstrated with three industrial case studies.
Index Terms—Fault diagnosis, harmonics, limit cycles, nonlin-
earities, spectral analysis, surrogate data, time series.
I. INTRODUCTION
I
T IS important to diagnose and cure oscillations in a con-
trolled process because a system running steadily without
oscillation is more profitable and safer [1]. A feedback control
loop containing a nonlinearity such as a sticking valve often ex-
hibits self-generated and self-sustained limit cycle oscillation
[2]–[4], and many surveys have shown that these oscillations
are a significant industrial problem [5], [6]. The situation is
made worse when the oscillation propagates throughout a dy-
namic system such as a chemical plant where it can become
widespread due to physical coupling and recycles. Rapid de-
termination of the source of a system-wide oscillation allows
maintenance activity to focus on the root cause [7]. This article
presents a method aimed at that objective and presents three in-
dustrial case studies in which the method successfully found the
root cause.
The time trend of measurements from a limit cycle oscilla-
tion is a nonlinear time series, i.e., it cannot be described as the
output of a linear system driven by white noise. A nonlinearity
test from Kantz and Schreiber [8] has been adapted for the de-
tection of limit cycle oscillations and guidelines for its applica-
tion to process data have been devised. The underpinning idea
in root cause diagnosis is that the nonlinearity is greatest at the
Manuscript received February 18, 2004. Manuscript received in final form
June 14, 2004. Recommended by Associate Editor A. T. Vemuri. This work
was supported in part by the Royal Academy of Engineering under the Foresight
Award and in part by the NSERC-Matrikon-ASRA Industrial Research Chair in
Process Control at the University of Alberta.
The author is with the Imperial College/UCL Center for Process Systems
Engineering, Department of Electronic and Electrical Engineering, University
College London, London WC1E 7JE, U.K. (e-mail: n.thornhill@ee.ucl.ac.uk).
Digital Object Identifier 10.1109/TCST.2004.839570
source of the problem. By source is meant a measurement as-
sociated with the single-input–single-output controller that has
been caused to oscillate by a nonlinearity in the loop. Justifica-
tion of that assumption is given in the next section.
The possibility of a nonlinearity test was outlined in an earlier
conference publication [9] and demonstrated in an application
at Eastman Chemical Company [10]. The key advance in this
brief compared with [9] is the exploitation of the cyclic nature
of the measurements to optimize the method. In [10], the focus
was on the solution to a particular industrial application and a
set of parameters was selected for the algorithm but without a
discussion of their optimality. A contribution of this brief is to
give fundamental insights into the method and to extend it to
additional case studies involving the detection of valve, sensor
and process nonlinearities. It gives an in-depth explanation of
how the method works when applied to oscillating disturbances
and the criteria by which the parameters of the algorithm may be
selected. Default parameters are suggested to facilitate routine
application of the method to a large-scale plant.
II. D
IAGNOSIS OF NONLINEARITY
A. Nonlinear Time Series Analysis
The waveform in a limit cycle is periodic but nonsinusoidal
and therefore has harmonics. A distinctive characteristic of a
nonlinear time series is the presence of phase coupling which
creates coherence between the frequency bands occupied by the
harmonics such that the phases are nonrandom and form a reg-
ular pattern. Nonlinearity may, thus, be inferred from the pres-
ence of harmonics and phase coupling.
Methods for nonlinearity detection in the time series include
the techniques using surrogate data [8], [11] which have been
used in applications ranging from analysis of EEG recordings
of people with epilepsy [12] to the analysis of X-rays emitted
from a suspected astrophysical black hole [13]. Surrogate data
are times series having the same power spectrum as the time
series under test but with the phase coupling removed by ran-
domization of phases. A key property of the test time series is
compared to that of its surrogates and nonlinearity is diagnosed
if the property is significantly different in the test time series.
Another method of nonlinearity detection uses higher order
spectra because these are sensitive to certain types of phase
coupling. For instance, the bispectrum [14]–[16] responds to
quadratic phase coupling in a signal such as
below in which
the phase of the frequency component at
is ,but
there is no bispectral response if
is a random phase
The bispectrum and the related bicoherence have been used
to detect the presence of nonlinearity in process data [17], [18].
1063-6536/$20.00 © 2005 IEEE
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 13, NO. 3, MAY 2005 435
A potential disadvantage of the bispectrum for detection of non-
linear limit cycle oscillations is that limit cycles may have sym-
metrical waveforms (e.g., a square wave or triangular wave) and
the bispectrum of a symmetrical waveform is zero. Zang and
Howell [19] have investigated the types of limit cycles that are
amendable to bispectrum analysis.
The presence of harmonics in a time series has also been used
successfully for diagnosis of SISO control loop faults [7], [20],
[21]. Finding harmonics requires signal processing to isolate the
spectral frequencies of interest and inspection to confirm that the
frequencies are integer multiples of a fundamental. The inspec-
tion is often undertaken by visual examination of the spectra
and is therefore unsuitable for a large-scale implementation in-
volving several hundred or even a thousand or more plant mea-
surements. Moreover, it is possible that components at the har-
monic frequencies are not phase coupled in which case the har-
monic signature will be a misleading indicator of nonlinearity.
B. Propagation of a Nonlinear Limit Cycle
Repair of a faulty control loop requires that the engineer
knows which control loop should be maintained. In the case of
a plant-wide oscillation, it can be a very difficult problem to
know which loop to work on because the disturbance from a
control loop in a limit cycle typically propagates plant-wide to
cause numerous secondary oscillation in other control loops.
An automated means is therefore needed to determine which
among all the oscillating control loops is the root cause and
which are secondary oscillations. Successful studies have used
the presence of prominent harmonics to distinguish the source
of a limit cycle oscillation from the secondary oscillations in a
distillation column in a refinery [22] and in a pulp and paper
mill [23]. The reason why secondary oscillations have lower
nonlinearity is that as the signal propagates away from its
source it passes through physical processes which give linear
filtering and which generally add noise. (The filter may be
assumed linear if the system is oscillating around a fixed oper-
ating point). Such a filter destroys the phase coherence of the
time trends and often reduces the magnitudes of the harmonics.
Thus, nonlinearity reduces as the disturbance propagates away
from the source and the time trend with the highest nonlinearity
is the best candidate for the root cause. The nonlinearity statistic
to be discussed in Section III can be used for such root-cause
diagnosis.
C. Surrogate Data Analysis
A time series with phase coupling is more structured and
more predictable than a similar time series known as a surro-
gate having the same power spectrum but with random phases
[11]. The spread of values of some statistical property of a group
of surrogate data trends provides a reference distribution against
which the properties of the test time series can be evaluated.
The techniques of surrogate data analysis have been widely
applied for detection of nonlinearity in time series [12], [13]. In
the process area, Aldrich and Barkhuizen [24] detected nonlin-
earity in process data by comparing a singular spectrum analysis
of the test data with those from linear surrogate data. Barnard
et al. [25] showed that identification of systems is possible by
Fig. 1. Test data and typical surrogate. The time trends are mean centered and
scaled to unit standard deviation.
using surrogate methods to classify the data, as well as to vali-
date models derived from these data.
Issues have been identified with the use of surrogate data with
cyclic time series [26], [27]. The surrogate is derived by taking
the discrete Fourier transform (DFT) of the test data, randomiza-
tion of the arguments followed by an inverse DFT. Nonlinearity
testing based on strongly cyclic data can give rise to false de-
tection of nonlinearity because when the time trend is strongly
cyclic then artifacts in the DFT due to end-matching effects in-
fluence the surrogates. A demonstration of the consequences for
strongly cyclic data are demonstrated in this article although in
practical applications the effect was found to have a minimal
impact, as will be discussed later.
III. M
ETHOD
A. Overview
The basis of the test is a comparison of the predictability of
the time trend compared to that of its surrogates. Fig. 1 illus-
trates the concept. The top panel is an oscillatory time trend of a
steam flow measurement from a refinery. It has a clearly defined
pattern and a good prediction of where the trend will go after
reaching a given position, for example at one of ringed peaks,
can be achieved by finding similar peaks in the time trend and
observing where the trend went next on those occasions.
The lower panel shows a surrogate of the time trend. By con-
trast to the original time trend the surrogate lacks structure even
though it has the same power spectrum. The removal of phase
coherence has upset the regular pattern of peaks. For instance,
it is hard to anticipate where the trajectory will go next after
emerging from the region highlighted with a circle because there
are no other similar peaks.
Predictability of the time trend relative to the surrogate gives
the basis of a nonlinearity measure. Prediction errors for the sur-
rogates define a reference probability distribution under the null
hypothesis. A nonlinear time series is more predictable than its
surrogates and a prediction error for the test time series smaller
than the mean of the reference distribution by more than three
standard deviations suggests the time trend is nonlinear.
436 IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 13, NO. 3, MAY 2005
B. Construction of the Data Matrix
Nonlinear prediction of time series was described by Sug-
ihara and May [28] to distinguish determinism from random
noise, and the field of nonlinear time series analysis and predic-
tion has been reviewed by Schreiber [29]. Rhodes and Morari
[30] gave an early process application of nonlinear prediction
where the emphasis was on modeling of nonlinear systems when
noise corrupts a deterministic signal.
Nonlinear prediction uses a data matrix called an embedding
having
columns each of which is a copy of the original data
set delayed by one sampling interval. For instance, a data matrix
with
3is
Rows of the matrix Y represent time trajectories that are seg-
ments of the original trend. Since the original data formed a con-
tinuous time trend the trajectories in adjacent rows are similar.
They are called near-in-time neighbors. Also, if the time trend is
oscillatory then the trajectories in later rows of Y will be similar
to the earlier rows after one or more complete cycles of oscil-
lation. For instance, if the period of oscillation is 50 samples
per cycle then
will be small, where is the 51st
row vector of
and is the first. Those rows are called near
neighbors.
C. Calculation of Prediction Error
Predictions are generated from near neighbors. Near-in-time
neighbors are excluded so that the neighbors are only selected
from other cycles in the oscillation. When
nearest neighbors
have been identified then those near neighbors are used to make
an
step-ahead prediction. For instance, if row vector were
identified as a near neighbor of
and if were 3 then
would give a prediction of . A sequence of prediction er-
rors can, thus, be created by subtracting the average of the pre-
dictions of the
nearest neighbors from the observed value. The
overall prediction error is the rms value of the prediction error
sequence.
The analysis is noncausal and any element in the time series
may be predicted from both earlier and later values. Fig. 2 il-
lustrates the principle using a time series from the SE Asia re-
finery case study where the embedding dimension
is 16 and
the prediction is made 16 steps ahead. The upper panel shows
the 100th row of the data matrix Y which is a full cycle starting
at sample 100, marked with a heavy line. Rows of Y that are
nearest neighbors of that cycle begin at samples 67, 133, 166,
199, and 232 and are also shown as a heavy lines in the lower
panel. The average of the points marked
in the lower panel
are then used as a prediction for the value marked
.
D. Data Preprocessing
Detection of plant-wide oscillation is now a solved problem
and is starting to be offered by vendors [31], [32]. The periodic
nature of the detected oscillation may be exploited in order to
Fig. 2. Illustration of the nearest neighbor concept. The highlighted cycles in
the lower plot are the five nearest neighbors of the cycle in the upper plot. The
average of the points marked
o
gives a prediction for the point marked
x
.
give robust default settings for the parameters. A summary list
is presented here and the detailed reasoning behind the recom-
mendations will be presented in Section IV. With the data pre-
processing steps indicated here the default parameters can be
used for any oscillating time trend.
1) The period of the plant-wide oscillation is determined.
2) The number of samples per cycle
is adjusted to be
no more than 25. The time trends are subsampled if
necessary.
3) The number of cycles of oscillation in the data set should
be at least 12 for a reliable nonlinearity estimate.
4) The selected data are end-matched to find a subset of the
data containing an integer number of full cycles. The algo-
rithm and other issues associated with end-matching are
explained in detail in Section IV-F.
5) The end-matched data are mean centered and scaled to
unit standard deviation. The sequence
de-
notes end-matched and preprocessed data in the following
sections.
E. Surrogate Data
Surrogate data are derived from the preprocessed and end-
matched time trend. Surrogate data have the same power spec-
trum as the time trend under test. The magnitudes of the DFT
are the same in both cases but the arguments of the DFT of the
surrogate data set are randomized. Thus, if the DFT in frequency
channel
is
then the DFT of the surrogate is
where is a phase selected from a uniform random distribution
in the range
. The aliased frequency channels
above the Nyquist sampling frequency have the opposite phase
added. If the number of samples
is even and if the frequency
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 13, NO. 3, MAY 2005 437
channels are indexed as 1to the Nyquist frequency is in
channel
and the alias of the th frequency channel is
channel
. Then
and to
If is odd
and to ceil
where ceil is the rounded-up integer value of .
Finally, the surrogate data set is created from the inverse
Fourier transform of the phase randomized DFT.
F. Nonlinearity Test
The nonlinearity test requires the determination of mean
square prediction errors of surrogates. The statistical dis-
tribution of those errors gives a reference distribution. If the
test data prediction error lies on the lower tail of the reference
distribution then the test signal is more predictable and non-
linearity is diagnosed using the following statistic based on a
three-sigma test
where is the mean square error of the test data, is
the mean of the reference distribution and
its standard
deviation. If
then nonlinearity is inferred in the time
series. Larger values of
are interpreted as meaning the time
series has more nonlinearity, those with
are taken to be
linear.
It is possible for the test to give small negative values of
.
Negative values in the range
are not statistically
significant and arise from the stochastic nature of the test. Re-
sults giving
do not arise at all because the surrogate
sequences which have no phase coherence are always less pre-
dictable than a nonlinear time series with phase coherence.
G. Algorithm Summary
Step 1) Form the embedded matrix from a preprocessed and
end matched subset of the test data
Step 2) For each row of Y find the indexes
of nearest neighbor rows having the
smallest values of subject to a near-in-time
neighbor exclusion constraint
.
Step 3) Find the sum of squared prediction errors for the test
data
Step 4) Create surrogate prediction errors by ap-
plying steps 1 through 3 to
surrogate data sets.
TABLE I
S
UGGESTED DEFAULT
VALUES FOR
PARAMETERS
Step 5) Calculate the nonlinearity from
IV. DEFAULT PARAMETER
VALUES
A. Default Parameter Values
Empirical studies have been carried out to ascertain the sen-
sitivity of the nonlinearity index to the parameters of the al-
gorithm. Reliable results have been achieved using the default
values shown in Table I. The floor function in the third row
indicates that for noninteger values of
then is set to the
rounded-down integer value of
.
The next subsections explore each one of these recommen-
dations showing why they were selected. The time trends from
Fig. 3. were used for the evaluation (they are from the indus-
trial case study in Section V-B). Fig. 3 shows mean centered
data normalized to unit standard deviation while the spectra are
scaled to the same maximum peak height.
The time series of the first three measurements are nonlinear
because they are close to the root cause. Their spectra have har-
monics and the phase patterns are not random. The last two are
far from the root cause and are linear. The data have an os-
cillation period of 16.7 sampling intervals and the conditions
used were varied around default values of
12
16 8 and 50.
B. Number of Samples Per Cycle
It is practical to limit the number of samples per cycle
.
There is a tradeoff between the number of samples needed to
properly define the shape of a nonsinusoidal oscillation on the
one hand and the speed of the computation on the other. The
algorithm requires a distance measure to be ascertained between
every pair of rows in the embedded matrix and the time taken
for the computation increases as
, where is the
total number of samples in the time trend. Therefore the number
of samples per cycle
and the number of cycles cannot be
increased arbitrarily. Data sets of 200 samples
8.2
24 and 418 samples ( 16.7 and 25) gave successful
results in the industrial case studies reported in Section V.
It would be infeasible to operate with fewer than seven sam-
ples per cycle because harmonics would not be satisfactorily
captured. With
7, any third harmonic present is sampled at
438 IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 13, NO. 3, MAY 2005
Fig. 3. Time trends and spectra of the data used for detailed evaluation.
2.33 samples per cycle which just meets the Nyquist criterion of
two samples per cycle. The reason for focusing the recommen-
dation on the third harmonic is that it is the most prominent har-
monic in symmetrical oscillations having square or triangular
waveforms.
C. Embedding Dimension and Prediction Horizon
Fig. 4(a) shows the effect of changing
and . They were
kept equal to each other and both were varied together. The
threshold of nonlinearity
1 is also shown in the plot
(horizontal dashed line) as well as the
16 default for
the data set (vertical dashed line). Once
becomes larger than
half a cycle of the oscillation, in this case when
, the re-
sults for the nonlinearity index
become quite steady while for
small values of
the index falls toward the 1 threshold.
An aim of the work presented here is to give reliable default
values that are easy to determine. Determination of the period
of oscillation
is becoming a standard component of controller
performance tools [31], [32]. Therefore the recommendation to
set
floor is robust because it is in the steady region of
Fig. 4(a) and is easy to implement because
is already known.
Fig. 4. Effects of parameters of the algorithm on the calculated nonlinearity
index. The vertical dashed lines in (a)–(c) show the recommended default
values.
The poor performance with small values of arises because
of the phenomenon of false near neighbors [33], especially when
the time trends have high frequency features or noise. The upper
panel in Fig. 5 shows an example of what can happen when the
embedding dimension is small, in this case
2. The rows of
the Y matrix starting at sample 159 comprises just two samples,
159 and 160, shown as small square symbols. Near neighbors
are shown in the lower panel, these are two-sample segments
of the time trend whose values are similar to samples 159 and
160. However, some of them are false neighbors because they
are not from matching parts of the trend. The average of the
points marked
are used as a prediction for the value marked
, but some of them such sample 223 which is based on a false
neighbor are not accurate.
D. Number of Cycles and Near Neighbors
Fig. 4(b) shows a plot of the number of cycles of oscillation
presented for analysis versus the nonlinearity value.
