DEVELOPING COMPUTER SOFTWARE FOR
CONTROLLER PERFORMANCE
MONITORING
By
QING LI
Bachelor of Science
Zhejiang University
Hangzhou, China
1985
Doctor of Philosophy
Oklahoma State Cniversity
Stillwater, Oklahoma
2002
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
in partial fulfillment of
the requirements for
the Degree of
MASTER OF SCIENCE
December, 2003
DEVELOPING COMPUTER SOFTWARE FOR
CONTROLLER PERFORMANCE
MONITORING
Thesis Approved:
II
ACKNOWLEDGMENTS
I would like to express my sincere appreciation to my major advisor Dr. John P.
Chandler, who inspired my interest in computer science, and has been a persistent
supporter ofmy research work in software development for controller performance
monitoring. My appreciation also goes to Dr. R. Russell Rhinehart, who introduced me
to the exciting area of controller performance monitoring, and provided me with
intelligent and constructive guidance in my previous study as a Ph.D. student in School of
Chemical Engineering at OSu. I thank Dr. Nohpill Park and Dr. Douglas Heisterkamp
for not only teaching me courses but also being on my Graduate Advisory Committee and
giving me valuable advice.
My special thanks go to Dr. Thomas W. Gedra in the School of Electrical and
Computer Engineering at OSU, who not only provided me with the financial support by
hiring me as a research assistant or a teaching assistant throughout the whole period of
my study in Computer Science, hut also provided me with the chance to write VI (Virtual
Instruments) software for the OSU Power Engineering Laboratory, and the chance to
solve real problems in load flow and optimization of combined electric networks and
natural gas networks. Thanks to Dr. Gedra, my years as a computer science graduate
student have been productive and enjoyable.
111
I also want to thank my friend, soontobe Dr. Seungwon An, who is Dr. Gedra's
Ph.D. student, for his help when I worked as a research assistant and teaching assistant in
the School of Electrical and Computer Engineering.
I thank my wife, Wenxia Peng, for her love, encouragement, and precious support
in my study and life. I also wish to express my gratitude to my parents, Zewen Li and
Guohua Liu, who have encouraged and supported me since my childhood. Great thanks
go to my mother for her love and numerous assistances throughout my whole life.
Finally, my two precious children, Grace Jiayun Li and Benjamin Jiayi Li, deserve my
special thanks because they brighten my life.
IV
Chapter
TABLE OF CONTENTS
Page
1. INTRODUCTION 1
1.1 Literature Review 2
1.2 Proposed Work 4
2. CONTROLLER PERFORMANCE MONITOR 5
2.1 Brief Description 5
2.2 Run Length (RL) 7
3. STATISTICAL TESTS ON RUN LENGTH (RL) DISTRIBUTIONS 8
3.1 Differences in RL Distributions 8
3.2 Statistical GoodnessOfFit Tests 9
3.2.1 ChiSquare Test. 10
3.2.2 KolmogorovSmimov Test (KS Test) 11
3.2.3 Comparison of ChiSquare Tests and KS Test on RL Distribution 14
4. ALGORITHM ANALYSIS .20
4.1 Old Algorithm for Controller Perfonnance Monitor. .20
4.2 New Algorithm for Controller Performance Monitor 23
4.2.1 Recurrent Equations for ChiSquare Test Statistic Calculation 25
4.2.2 Recurrent Algorithm for KS Test Statistic Calculation 30
4.3 Comparison ofNew Algorithms and Old Algorithm 37
4.3.1 New Algorithm Using ChiSquare Tests 37
v
Chapter Page
4.3.2 New Algorithm Using KS Tests 38
4.3.3 Runtime Comparisons 39
5. STATISTICAL MODELS FOR RL DISTRIBUTIONS .41
5.1 RL Distribution Model for White Noise Random Processes .42
5.2 RL Distribution Models for Correlated Processes .43
5.3 RL Distribution Probability Functions for Moving Average (MA) Processes ..47
5.3.1 RL Distribution Probability Functions for 1stOrder MA Processes .47
5.3.2 RL Distribution Probability Functions for 2ndOrder MA Processes 57
5.3.1 RL Distribution Probability Functions for HighOrder MA Processes 62
6. DEMONSTRATION OF CONTROLLER PERFORMANCE MONITOR 66
5.1 Description of Experimental Unit 66
5.2 Experimental Procedure 68
5.3 Results and Discussion 68
7. CONCLUSIONS AND RECOMMENDATIONS 74
BIBLIOGRAPHY 75
APPENDICES 77
A.I Code for Old Algorithm of Controller Performance Monitor 77
A.2 Code for New Algorithm Using ChiSquare Test 79
A.3 Code for ~ew Algorithm Using KS Test.. 83
A.4 Code for Calculating Theoretical RL Distribution Probability Functions for
1stOrder MA Processes 87
A.4 Code for Calculating Theoretical RL Distribution Probability Functions for
VI
Chapter Page
2ndOrder MA Processes , 89
vii
LIST Of TABLES
Table Page
2.1 Runtime Comparison ofNew Algorithms and Old Algorithm 39
5.1 Theoretical RL Distributions of 151Order MA Processes 54
6.1 Reference RL Distribution Bin Groups, Water Flow Control 69
Vlli
LIST OF FIGURES
Figure Page
2.1 Flow Chart of Controller Performance Monitor 7
2.2 Zero Crossing and Run Length 8
3.1 RL Distributions with and without Oscillations 10
3.2 Discrete Cumulative RL Histograms with and without Oscillations 15
3.3 ZeroMean Random Data, ChiSquare Test Statistic, KS Test Statistic,
and Their Critical Va1ues 17
3.4 NonZeroMean Random Data, ChiSquare Test Statistic, KS Test Statistic,
and Their Critical Values 18
3.5 ZeroMean Random Data Plus Oscillations, ChiSquare Test Statistic, KS Test
Statistic, and Their Critical Values 19
5.1 Comparison of Theoretical RL Distribution and Sample RL Distribution of
a 151Order MA Process with bdbo= 0.5 56
5.2 Comparison of Theoretical RL Distribution and Sample RL Distribution of
a 151Order MA Process with b]/bo= 1 57
5.3 Comparison of Theoretical RL Distribution and Sample RL Distribution of
a 151Order MA Process with b]/bo = 2 58
5.4 Comparison of Theoretical RL Distribution and Sample RL Distribution of
a 2ndOrder MA Process with bl/bo = 1 and b1/bO= 1 62
6.1 TwoPhase Flow Experimental Unit with One Water Flow Control Loop and
IX
~~ p~
Two Air Flow Control Loops 68
6.2 Sample Reference RL Distribution for Water Flow Control Loop 70
6.3 Experimental Data, ChiSquare Test Statistic, and Flagging When a Sample
RL Distribution is Used as Reference Distribution 72
6.4 Experimental Data, ChiSquare Test Statistic, and Flagging When a Theoretical
RL Distribution is Used as Reference Distribution 73
6.5 Experimental Data, ChiSquare Test Statistic, and Flagging When a Grace
Peliod is Used Along with a Theoretical Reference RL Distribution 74
x
CHAPTER 1
INTRODUCTION
1.1 Background
Chemical plants are controlled by a large number (102 or 103
) of process
controllers to maintain safe and economic operation. What a controller does is to keep a
process variable (called the controlled variable, or CV), such as temperature or flow rate,
at a desired value (called the setpoint, or SP). Due to measurement errors, external
disturbances or setpoint changes, there will exist differences (called errors) between the
setpoint and the measured controlled variable. As long as the error is not zero, the
controller will take some control action, that is, adjust a variable (called the manipulated
variable, or MV), to try to make the error zero. In reality, the errors will rarely be exactly
zero because there always ex ist external disturbances and random measurement noises,
which make the measured controlled variable deviate from the setpoint.
A good controller takes reasonable control actions to keep the errors small. A
poor controller may take either too sluggish control actions, which reduce the errors too
slowly, or tooaggressive control actions, which cause oscillations around zero value in
the errors due to overreaction. Both sluggish and aggressive controllers cause larger
errors than a good controller.
Human evaluation of the control loop performance not only is tedious and timeconsuming,
but also requires expertise and experience. Some poor controller
perfomIance is not easy to detect visually. For example, when the measurements are very
noisy, such as flow rate measurements, the oscillations due to a tooaggressive controller
may be hidden in the noise and not visually obvious to operators. In practice, poor
controller performance often has existed for a quite long time before being detected.
An automatic controller performance monitor is desirable in process industry
because (1) a poorly performing controller will increase operating costs, and (2) manual
detection ofpoor controllers is timeconsuming, tedious, and often timedelayed.
The objective ofthis work is to develop a software tool for automatic detection of
poorly performing controllers used in the process industry.
1.2 Literature Review
Traditionally, a controller's performance is evaluated by step tests [1]. After
introducing a step change in setpoint or disturbance, the errors between the setpoint (SP)
and the controlled variable (CY) are recorded and used to measure its performance. The
common measures are: mean or integral of absolute error (rAE), squared error (ISE), OT
timeabsolute error (ITAE). Other measures include overshoot, settling time, decay ratio,
etc. [2]. The traditional measures are simple to compute, but they require plant tests,
which disturb the normal plant operation and thus are undesirable.
The process industry needs the monitoring techniques that do not require plant
tests, i.e., use routine plant operation data only. The minimum variance control (MVC)
based performance monitoring techniques, which were first introduced in 1989 by Harris
[3], use routine plant data only without requiring plant tests, and thus became popular in
industry. However, most of these MYCbased techniques require knowledge of the
process delay, which may not be easy to obtain without disturbing the plant, and is
2
subject to change during plant operation. Reviews of the MVCbased techniques can be
found in [46].
There are other purely databased monitoring techniques. An automated online
goodnessofcontrol performance monitor was proposed by Rhinehart et aI. [79]. The
method uses a computationally simple, robust statistic, called the rstatistic, which is
defined as the ratio of the expected variance of the deviation of the controlled variable
from the setpoint to the expected variance based on the deviation between two
consecutive process measurements.
Many methods to detect oscillations in poor control loops have been proposed,
such as methods using the integral of absolute error between zero crossings (lAEZC) [10,
11], using the autocovariance function (ACF) [12], and using zero crossings of ACF [13].
The reference model approach is used in a relative controller performance
monitor proposed by [14]. The monitor compares the perfonnance of a control loop
under monitoring to that of a reference model  to measure the performance of the actual
control loop relative to the good performance represented by the reference model. It
provides a measure (relative performance index  RPI) of control loop performance
relative to a reference model of acceptable control. The reference model simulates the
controlled variable output of a userdefined, acceptably tuned control system.
Recently, a new technique, on which this work is based, uses a distributioncomparison
approach to monitor the controller performance [15]. This technique detects
poor controller performance based on statistical differences in runlength (RL)
distributions. The RL is defined as the time period (number of sampling periods)
between two consecutive zerocrossings (sign changes) of the controller actuating error
3
signals (setpoint minus controlled variable). The histogram ofRL values under different
control performance, such as sluggish control, aggressive control (oscillations) and good
control, are significantly different, and therefore are used to detect poor controller
performance. The monitor does not require either process knowledge or process model,
and it uses only routine plant operation data.
1.3 Proposed Work
In this work, we develop a software tool for automatic detection of poor controller
performance, using the basic idea of the RLbased controller performance monitoring
teclmique proposed in [15J, but making several significant improvements upon it.
The work can be divided into the following three major parts:
• Design and implement an efficient new algorithm for the monitoring software
that maximizes the computation speed with minimal memory;
• Build statistical models, and derive theoretical probability functions for the
reference RL distributions for correlated data;
• Compare the two goodnessoffit test methods for RL distribution
comparisons: the chisquare test and the Kolmogorov  Smimov test (KS
test).
4
CHAPTER 2
CONTROLLER PERFORMANCE MONITOR
2.1 Brief Description
The controller performance monitor proposed in this work is an improvement
upon the original controller performance monitor in [15]. Figure 2.1 shows the
simplified flow chart of the original monitor.
The monitor performs the following steps:
(1) At each sampling time, measure the process controlled variable (CY) and
setpoint (SP), and calculate the error (SP  CV).
(2) Calcul.ate the RL distribution (to be explained later) from a window ofthe past
N error values.
(3) Apply the statistical goodnessoffit test (the chisquare test, or the KS test) to
detennine whether the obtained RL distribution is significantly different from a reference
RL distribution, which has previously been built from data collected during a good
control period as judged by human experts.
(4) Detect poor performance based on the above statistical test results and a grace
period. (One choice for the grace period is the window length plus the settling time of a
welltuned controller after a step change.)
(5) At the next sampling time, repeat the whole procedure.
5
Reference RL
distribution
Start
Sample CY, SP
Calculate run length (RL)
distribution, using past N
actuating errors,
e = (SP  CY).
Significant
difference by
statistical test?
No
Reset violation
counter
Yes
Increment violation
counter
> Grace
period?
No
Repeat for
next data
Figure 2.1 Flow Chart of Control Performance Monitor
6
2.2 Run Length (RL)
Run length (RL), in this work, is defined as number of samples between two
consecutive zero crossings. A zero crossing is defined as each time the error signal
changes its sign from + to , or  to +. For a window of a data sample sequence, the
samples before the first zero crossing and the samples after the last zero crossing may
also be treated as a run and thus have a run length.
Figure 2.2 illustrates the zero crossing and run length concepts. A vertical bar
represents a zero crossing, and number of samples between two consecuti ve zero
crossings is the run length.
ample
Error signal
Each bar represents a zero crossing
Run Length, RL
/ 
0 0 I 0 0
00 0 ~o o 0
0 0 0 0 k> 0
o c o 000 0
000
p 0 0 0
0 00 S
0 0 0
o
Figure 2.2 Zero Crossing and Run Length
7
CHAPTER 3
STATISTICAL TESTS ON RUN LENGTH (RL) DISTRIBUTIONS
3.1 Differences in RL Distributions
The RL distributions are different when the controller perfonnance is different.
Figure 3.1 shows three RL distribution histograms. The straight line represents the
theoretical RL distribution (to be discussed later) of a random signal with nonnal
distribution with zero mean and unit variance. The second distribution (represented by *)
is obtained from a random signal with normal distribution with zero mean and unit
variance. This RL distribution is similar to the RL distributions when the controller is
performing well The third RL distribution (represented by o) is calculated from a
combined signal, which is the sum of the random signal (the same as the second
distribution) and a sinusoidal signal with an amplitude equal to the standard deviation of
the random signal and a period equal to 10 samplings. This RL distribution is similar to
that obtained when there are oscillations in the controlled variable due to a tooaggressive,
poor controller. The error bars represent the 95% confidence intervals below
and above the mean values the histograms obtained from 1000 sets of samples with each
set of sample having 1000 run lengths.
8
Run Length Histogram
10
3
r,.,.,,;::::J::====c::::==::r.:==::::;l
 Theoretical
+ Random Sampies
<) Random+Oscillalions
"J.:
l
2 3 4 5 6 7 8 9
Run Length, samples
Figure 3.1 RL Distributions with and without Oscillations. The RL distribution with
oscillations is obtained using the data generated from e = en + crsin(27tk/lO), where eo is a
random variable with a nonnal distribution N(O, cr2
), and the sinusoidal osci Ilations have
an amplitude equal to eo's standard deviation cr and an oscillation period of 10 samplings.
The straight line represents theoretical RL distribution of the errors eo without oscillations
or offsets, and is shown for comparison. The error bars are 95% confidence intervals
around the mean values of the histograms obtained from 1000 sets of samples with each
sample set having 1000 run lengths.
We perform a statistical goodnessoffit test to see if a sample RL distribution is
statistically different from a reference distribution.
3.2 Statistical GoodnessOrFit Tests
Statistical goodnessoffit tests are commonly used to decide whether a data
sample comes from a population with a specific distribution [16, 17]. Goodnessoffit
tests are a form of hypothesis testing, where the null and alternative hypotheses are flo:
9
Sample data come from the stated distribution, and HI: Sample data do not come from
the stated distribution. Goodnessoffit tests indicate whether or not it is reasonable to
assume that a random sample comes from a specific distribution.
There are several ways to test goodnessoffit, such as the chisquare test, the
KolrnogorovSmirnov test (KS test), and the AndersonDarling test [16]. The KS test is
usually used for a continuous distribution, but there is a corresponding KS test for a
discrete distribution function [18]. The AndersonDarling test is only applicable to a
limited number of continuous distributions.
3.2.1 ChiSquare Test
The original controller performance monitor in [15] uses the statistical chisquare
goodnessoffit test on the differences in two RL distributions to determine whether the
RL distribution under current operation is significantly different from a reference RL
distribution that represents good control performance.
For n independent observations of runs, let OJ be the observed number of runs
with a run length of i, or RL =i, and Ej be the expected number (Ej =n*pj, where Pi is the
probability of get a run with length i). Assuming OJ is a Poisson distrihuted random
variable, then, the following statistic approximately follows the chisquare distribution
with (k  1  p) degrees of freedom:
2 f.(O _E.)2
X =~ , I
;=1 E,.
where OJ is the observed number of run lengths in bin i of a RL histogram, and 0 1 is
(3.1)
assumed to be Poissondistributed; Ej is the expected number of run lengths in bin i; k is
the number of bins of the RL histogram; p is the number distribution parameters
10
estimated from the sample, which is zero because our expected distribution is
predetennined. This chisquare test is also called Pearson's goodnessof fittest.
The advantage of the chisquare goodnessoffit test is that it can be applied to
any univariate distribution for which you can calculate its distribution function [19]. This
property that the test statistic is independent of the underlying distribution is important
because the RL distribution in practice may not follow any specific distribution function.
The major disadvantages of using the chisquare test on the RL distributions are:
(1) some bins of the RL histogram may need to be grouped because the bins
corresponding to large RL values may have a very small number of expected
observations, and in order for the test statistic to have a valid chisquare disttibution, it is
recommended that each group have at least 3 to 5 expected observations [16, 17]; (2)
there are no wellagreed rules for grouping the bins; (3) the value of the test statistic will
be affected by how we group the bins, and thus the test statistic is not unique because bingrouping
is not unique. We use a rule ofthumb to group all bins into about 5 classes with
each class having approximately equal probability.
3.2.2 KolmogorovSmirnov Test (KS Test)
The KS test statistic S [18] is defined as
j
S = max I(O;  E;)
1:s; j :s; m ;=\
(3.2)
The S statistic is the maximum difference between the observed cumulative distribution
function and the expected one. Thus, the S statistic can also be written as
11
where
S = max IFo (j)  FE (j)1
l5,j5,m
j
Fo(j) = LO;
;=1
i
FEU) = LE;
;=1
(3.3)
(3.4)
(3.5)
FoU) is the observed, cumulative histogram for RL 5, j, FEU) is the expected, cumulative
j
histogram for RL 5, j, and Fo (j)  FE (j) =L(0;  Ej ) is the differences between the
;=1
observed and expected cumulative histograms for RL 5,j. Therefore, S is the maximum
differences between the observed and expected cumulative histograms fOT all RLs. The
maximum positive difference, S+, and the maximum negative difference, S, are used fOT
onesided tests.
Although the most common use of the KS test is to test continuous distributions,
there is a corresponding discrete KS test for discrete distributions [18,20]. The discrete
KS test is one ofthe goodnessoffit methods based on the empirical distribution
function (EDF) for discrete data. Note that the RL distribution is a discrete distribution.
An attractive feature of the KS test is that it is an exact test (the chisquare
goodnessoffit test depends on an adequate sample size for the approximations to be
valid) [21 J. Also, the distribution of the KS test statistic itself does not depend on the
underlying cumulative distribution function being tested, so, the discrete KS test can be
performed on any RL distributions, which in general do not follow any specific known
distribution.
12
The null hypothesis Ho is that the sample RL distribution follows the reference
distribution (representing good controller perfonnance), and the altemative hypothesis HI
is that the sample RL distribution does not follow the reference distribution. The null
hypothesis Hocan be expressed as
Prob(an observed RL = i) =Pi, i = 1,2, ... , m
where Pi is the probability of RL = i in the underlying population distribution, and m is
the upper limit ofnm lengths considered for the statistical test. In the controller
performance monjtor, Pi is completely specified by the reference RL distribution
representing good controller performance.
The reference RL distribution can be established from the actual historical data
collected during a good control period, or from some theoretical RL distribution, which
will be discussed in a later chapter.
If the calculated statistic S is greater than a critical val.ue with certain significance
level, the null hypothesis wi 11 be rejected. The tables for critical values of the statistic S
can be found in many references [18, 20].
Figure 3.2 shows the cumulative histograms for the three RL distributions (whose
histograms were already shown in Figure 3.1), including the cumulative theoretical RL
histogram of a random signal, the cumulative sample RL histogram of 1000 RLs for a
random signal, and the cumulative sample RL histogram of 1000 RLs for a combined
random and sinusoidal signal. Also shown are the 95% confidence intervals for the two
sample histograms. The figure shows the differences in the cumulative RL histograms
due to existence of oscillations. It can been seen that for small run lengths (RL ::; 5) the
13
theoretical values of the cumulative histogram are outside of the 95% confidence
intervals with oscillations, but inside those without oscillations,
Cumulative Histogram of 1000 Run Lengths
2 8 9 '10
 .0. ,'J'
~
I
ill
I
ID
r
...L 1

I ' 1
1 • Theoretical
Random Samples
i 0 Random+Oscillations
1
p 1I
300o
1100
1000
900
>>
(j
c 800 0)
~
0
~
lL
0) 700
.~
(C
::)
E 600
::)
u
500
400
3 4 567
Run Length, samples
Figure 3.2 Discrete Cumulative RL Histograms with and without Oscillations.
These cumulative histograms are calculated from the histograms shown in Figure 3. ] .
The error bars are the 95% confidence intervals.
3.2.3 A Comparison of ChiSquare Tests and KS Tests on RL Distributions
A random signal with zero mean and unit variance has the following theoretical
run length probability function (see Chapter 5 for details):
Pr(RL = r) = (~J r = 1,2, ... (3.6)
14
Among all runs, Y2 have a run length of 1, ~ have a run length of2, and so on. If we take
64 random samples of runs, we expect to observe the following numbers of runs with
respective run lengths:
Run Length
1
2
3
4
5
6
7
8
9
10
Expected Observations, Ej
32
16
8
4
2
1
0.5
0.25
0.125
0.0625
To make Equation (3.1) a valid chisquare test statistic, we follow the recommendation
that each bin have 3 or more expected observations [16, 17]. Therefore, bins with RL = 5
or more must be grouped together as one class before performing chisquare tests.
Figure 3.3 shows zeromean random data (with a n0n11al distribution with mean =
oand variance = 1), and the values of the chisquare test statistic and KS test statistic,
along with the critical values with a significance level of 0.01 and a data window size of
50. We can see that most of the time, both statistics are below its critical values,
indicating both tests show no significant differences between the sample RL distributions
and the theoretical reference RL distributions.
15
5
Random Samples, ChiSquare Test Statistic. KS Test Statistic
,r,.,,
5 \
50 100 150 200 250 300 350 400 450 500
A5 uen 7i'i r"i ';
~ 10 '! II (f) ':'1 'I .Q...) " 1)1 11 .1, r, co " : ,,.:'1 1 f r If •\; I ~ (\ r\l'l~
5
r',l
:' '; I.
:::> I I' '. ,
(', 'i'~ '. }l', I ! \'\' ,', I I I; 0 r : i\/., i j :'/'n.\ It,i '~, .. ,
'I, J
(f) I . d'
! T '" !~ ,10\ r ,fl,
, ; 'll'",(1'l i
\/.......'.f· \\ (fll) " I:.'
... 1(~,,,. J f 'V\~"I
..c .......1:'' ~ . ~'...'.'_., ." '.' ,
u 0 
50 100 150 200 250 300 350 400 450 500
10
u
~
~ 5
(f)
u:
~
I:',
, '
,.: i
,
200 250 300 350
Sample Time. k
400 450 500
Figure 3.3 ZeroMean Random Data (with a nonna! distribution with mean = 0 and
valiance = 1), ChiSquare Test Statistic, KS Test Statistic, and Their Critical Values
(with a Significance Level of 0.01 and a Data Window Size of 50).
Figure 3.4 shows nonzeromean random data (with a norma! distribution with
mean = 1 and variance = 1), and the values of the chisquare test statistic and KS test
statistic along with the critical values with a significance level of 0.0 1 and a data window
size of 50. We can see that many val lies of the chi square statistic are above its critical
values, indicating significant differences in RL distribution and poor control
performance. However, values of the KS statistic are still below its critical value,
indicating no significant differences, which is wrong.
16
2 ~."__'__L__'__' L__'__l_~
50 100 150 200 250 300 350 400 450 500
u 150
~
2 100
150 200 350 400 450 500
,. I,
~J I II
I I ..,
'. ,
"I
300
," ...... / ~t .I
250
, ,
.: I
..l7 ..'
i
".l',
,..,,
. 0
·'i;
__".J'! j!1 _
50
(f)
~
ro
::l
0(
f)
I
..c
U
u
Ul
5
~j'1:.1.' •
, . I;"
.....,. ,.'
100 150 200 250 300 350 400 450 500
o'_Jc__'__'__'__'__'__'_.L..._J
50
Sample Time. k
Figure 3.4 NonZeroMean Random Data (with a nonnaI distribution with mean = 1
and variance = 1), ChiSquare Test Statistic, KS Test Statistic, and their Critical
Values (with a Significance Level of 0.0 I and a Data Window Size of 50).
Figure 3.5 shows zeromean random data (with a nonnaI distribution with mean :..c
1 and variance = 1) plus sinusoidal osci Ilations (period = 10 samples and amplitude = 2),
and the values of the chisquare test statistic and KS test statistic along with the critical
values with a significance level of 0.0 1 and a data window size of 50. We can see that
values of the chisquare statistic are above its critical values, indicating significant
differences in RL distribution and poor control performance. The values of the KS
17
statistic are also above its critical value most of the time, indicating significant
differences in RL distribution, i.e., poor control performance.
Random Samples. ChiSquare Test Statistic, KS Test Statistic
5r,rr,~_.____._____.__r____,
fl ',!lli'l,!~rl\~11;1.1~. 'l!I'~II~\,1 rll!I'~jj), 'I ~ ;i;1 :1 Il'\ I! II!\ !I!;! 1\ rill :I' i ': !J~ :i I; :I~: Ii ;.'1, i \~ I ,Ii;i 1\ \'~ I, 1\/ jq! 1\ r\ , J 1\ 1\ ! I~
~ 0 llJ u~ ,I~r'\.' itI iJ \1\ I'\~ \; Ij 1(1,lti1'
\1 IIUIj! 'i II I'Ii :~ i\1 \' Ii \1 \/1)' '(' ~i\ fu II 'dIllr~~ 'I!II:I':i l\iJ'
I Ii " \1 I f '.I;'!' :1 W';'1 ;' ;{ Ii \1 ! 1 I! f' ,: I ,I II ' I, '. I '" I II ~ I';".' ~ I i,
I. , ': {' I I I ii, ' r', I '! I I 'I I ,
5 r I I I I L I !
50 100 150 200 250 300 350 400 450 500
350 400 450 500
I'
,'...J
150 200 250 300
'j
" '
" , ,'f j, I' II ;" . 'l': \
.,l,:/, \.'..... " ,r'd ;,,", \
("/1
___.__________ _ __~""."';c.;.:'I oL_'__' L__". '__'__'_' _
50 100
~
~ 100
£o
uen
4 __'_''_'_ .1. _ _'_''_'__', _
50 100 150 200 250 300 350 400 450 500
Sample Time, k
Figure 3.5 ZeroMean Random Data (with a nonnal distribution with mean = 0 and
variance = 1) plus Oscillations (period = 10 samples and amplitude = 2), and Values of
tbe CbiSquare Test Statistic and KS Test Statistic Along with the Critical Values
with a Significance LeveJ of 0.01 and a Data Window Size of 50.
This comparison shows that the chisquare tests are more sensitive for detecting
the common poor control perfonnance: oscillations and longtime offsets.
The major reason why the chisquare test shows more testing power than the KS
test is that the chisquare test statistic we lise here is based on the histogram, but the KS
IS
test statistic is based on the cumulative histogram. The nature of cumulative sum for
each bin averages out the variations, and thus makes the KS test statistic much less
sensitive to changes in RL distribution than the chisquare test statistic, which is based on
histogram not cumulative histogram. Although many references such as [16] state that
the KS test is in general more powerful than the corresponding chisquare test, it is said
under the context that the same cumulative histograms are being tested by both statistics
since the KS test statistic is defined on cumulative histograms only. Therefore it does
not apply here because the chisquare statistic we lIse is based on histograms and the KS
test we use is based on cumulative histograms. Therefore, our finding that the chisquare
test is more powerful than the KS test on RL distributions is reasonable.
19
CHAPTER 4
ALGORITHM ANALYSIS
The controller perfonnance monitor perfonns two major computation steps: RL
distribution computation and goodnessoffit statistical tests on RL distributions. The
input to the controller performance monitor is one data sample at each sampling time, and
the output is the statistical test results.
The original algorithm in [15], shown in Figure 2.1, has a lot of room for
improvement in terms of computation time and memory space. For example, in step (2),
the RL distribution is computed from the past N errors, and the whole computation is
repeated each time a new sample is available. A lot of computation time could be saved
if we exploit the fact that the N data used to build the new RL distribution is just a shifted
version ofN old data used to build the previous RL distribution, with only one new data
point added and the oldest data point discarded. That is, if we could use some recursive
way to compute the new RL distribution from the previous RL distribution, we can avoid
repeating the whole procedure, which is timeconsuming. In addition, memory use could
also be optimized. A new and efficient algorithm in terms of execution time and memory
space wi II be designed and implemented in this work.
4.1 Old Algorithm for Controller Performance Monitor
Before proposing the new algorithm to compute the RL distribution, the old
algorithm is shown first. The new algorithm will be compared with the old one to show
the improvements in computation speedup and memory reduction.
20
The old algorithm for the controller performance mor.ltor is described by the
following highlevel pseudocode.
for each data sample, i = 1 to n
take a window of the past N data
for each data in the window of size N I/compute RL histogram
if zero crossing occurs
increment the histogram of current run length
reset the cu rrent RL to 1
else
increment the current RL by 1
perform statistical goodnessoffit test
output test result
The computational complexity of the above old algorithm is O(n*(N+m)), where
n is the input data sample size, N is the data window size to build RL histogram, and m is
the upper limit of run lengths considered for the goodnessoffit test (all run lengths with
RL 2': m will be treated as RL = m).
A more detailed pseudocode for the old algOlithm is shown below (the complete
code in Matlab is shown in Appendix A.I):
Inputs:
Output:
x, a sequence of real numbers of size n, one data per sampling,
N, an integer, data window size, number of data samples used to
build one RL distribution,
m, an integer, upper limit of ron length considered, (m <= N).
S, a real number, the goodnessoffit test result.
21
CPM_Old(x, N, m) /lold algorithm for control performance monitor (CPM)
{
Line
1
2
real w [N};
int y[m];
lIa window of data of size N
lIthe RL histogram
34
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25 }
for k = 0 to N  1 //fill in the first window of size N,
w[k] = x[k];
for i = N to n  1 IIi is the sample time,
{ Ilstart computing RL histogram after read N data.
for k = 0 to N  2 Ilupdate the window of data
Ildiscard the oldest data by shifting left w by one place
w[k] = w[k + 1];
w[N  1] =xli]; Iladd the new data to the right end
int RL =1; Ilinitial value of RL =1
for k =0 to N  1 Ilinitialize histogram to 0 for each window
y[k] =0;
for j =1 to N  1 Ilfrom window of past data of size N
{ /Ifind all runs and RLs within the window of size N
if w[j] crosses zero
y[RL}++; !/increment histogram of current RL
RL = 1; Ilreset the RL to initial value 1.
else
RL++; lIincrement RL for the current run
}
y[RL]++; lIincrement histogram for the last, incomplete RL.
s =goodnessoffit test(y); Ilcompare RL histograms
output test result s;
}
The number of executable statements of the old algorithm, is, approximately,
Told = N + (n N)*(N + 1 + N + 2N + 01)
=N+(nN)*(4N+m+ 1)
= 4nN  4N2 + nm  Nm + n
22
Here, we just count the body ofloops (ignoring the loop head and statements outside
loops), and assume that the goodnessoffit test has a linear time in the number of bins
(with an upper limit of m).
The memory space required (ignoring the constant size, temporary memory), is
one array ofreal values of size N for storing the data window, and one anay of integers
of size m for storing the RL histogram, i.e.)
Mold =N*Sreal + m*Sinl>
where N is the data window size, m the upper limit of run lengths considered, Silll the
memory size for storing an integer datum, and Sreal the memory size required to store a
real number datum (float or double).
4.2 New Algorithm for Controller Performance Monitor
The major improvements of the proposed new algorithm are (1) remove the ilmer
forloop (between lines 13 and 20 in CPM_Old), (2) perfOInl tbe statistical tests only
when a new run is completed, rather than at each sampling, amI (3) perfonn statistical
tests recurrently, i.e., calculating the next test statistic from the previous RL distribution
and result, rather than doing it all over again for each new RL distribution. The
improvements in computation time and memory space are discllssed later.
The new algorithm uses an array (R[] of integers of size N) to store the sequence
ofRLs within the current data window of size N. Since the number of runs within a data
window of size N is less than or equal to N) the array R of size N is large enough to store
all RLs in the data window. Two pointers (or index variables) are used to point to the
begil1Jling RL (the oldest RL) and the last RL (the newest RL) within the data window.
".
As the data window is shifting, new RLs are added to the array R and old RLs are
removed from the array R. The pointers to the oldest and newest RLs are circulating
within the alTay R.
The highlevel pseudocode of the new algorithm is shown below:
for each data sample, i =1 to n
//update histogram & RL sequence after RLoldest passed the window
if the oldest, incomplete run has passed out of the data window
decrement histogram for RLoldest
runDropped = true
update the RLoldest pointer to the next oldest RL in R
if the newest data sample crosses zero
increment the histogram of current RL
ru nAdded = true
save current RL into R as RLnewest
update the RLnewest pointer to the next empty space in R
reset the current RL to 1
else
increment the current RL by 1
perform the statistical test
The last line "perform the statistical test" represents the code sections (to be
shown later) for performing the chisquare test or the KS test. Perfonning the statistical
test is necessary only in these three cases: (1) when one nm has dropped out of the data
window, and no nms have been added during window shifting to include the new data
sample; (2) when one run has been added to the window, and no runs have been dropped;
(3) when one run has been dropped and one run has been added at the same time. Jn
24
other cases, the test statistic remains at its previous value. Note that the incomplete runs
at the both ends of the data window are not used to build the RL histograms, and thus do
not affect the test statistic.
4.2.1/ Recurrent Equations for ChiSquare Test Statistic Calculation
For the chisquare tests, to avoid recalculating the chisquare test statistic for each
nm added or dropped, we need to change the fonnula shown in Equation (3.1). The
purpose is to calculate the statistic recunently. From Equation (3.1), we have
x? = f 0 2 20£ +£2 ~ ttl I
i=1 E,.
Since
and
K K 2: 0 ,.= IEi =nl{
i=1 i=1
(4.1)
(4.2)
(4.3)
where DR is sample size of runs, and Pi is the probability of a run falling in bin i, we have
(4.4)
If at time k, the test statistic is X.k
2
, and at time k + l, a mn within bin j is dropped
out of the data window, which is the first case of an RL distribution change, then we have
25
K 0 2 0 2
Xk2 =(1L~.1)+, I1Rk
lI Rk ;=1 Pi nRkPi
i~j ,
(4.5)
(4.6)
Since nR(k+l) = nRk  1 when a run is dropped, we have
(4.7)
test statistic for the first case of RL distribution changes, when one run is dropped out of
i' e." .".,...
(4.8)
K 2 (0 _1)2
2 ~ 0, j. 2
(lI Rk 1)%"+1 =(L,..)+ (nRIr. 1)
;=1 Pi Pi
iot i
Therefore, we have the following recurrent formula for updating the chisquare
and
Subtracting Equation (4.6) from Equation (4.8), we have
the data window,
,2 _ n Rk 2  20j + I 2n Rk  1
X.~+I  Xk + + ''=
(nRIr.  1) (n Rk  1)PinRIr.  I
(4.9)
For the second case, where a run is added to the data window, using a similar
method, we can obtain the following recurrent fonnula for updating the chisquare test
statistic,
26
or
for the third case, where exact one run is dropped in bin i and one run is added in
bin j, then, from Equation (4.6), we have
K 0 2 0 2 0 2
nRk X; = (L_s) +' +' ~ (n Rk ) 2
s~l. Ps Pi PI
s;<,
s~j
K 0 2 (0.1)2 (0.+1)2
nR(k+I)X;+1 =(I<)+' +'  (n R(k+I)2
s~1 Ps Pi Pi
5""
s$Oj
Since the total nm number does not change, i.e., IlR(k+l) ,.., nRk, subtracting the first
equation above from the second equation, we have
Therefore, the pseudocode for recurrent calculation of the chisquare test statistic is
if runDropped =true and nmAdded = false lIease 1
calculate new x2 from old one using Equation (4.9);
else if runDropped = false and nunAdded =true Ilcase 2
calculate new x2 from old one using Equation (4.10);
else if run Dropped = true and runAdded =true Ilcase 2
calculate new x2 from old one using Equation (4.11);
27
(4.11 )
?\fate that it takes a constant time to calculate new X2 values using Equations (4.9)
 (4.10), so if the chisquare test is used, the new algorithm for the controller
performance monitor has a linear runtime complexity, i.e., O(n).
A more detailed pseudocode of the new algorithm using the chisquare test is
shown below (see Appendix A.2 for the complete code of the new algorithm using the
chisquare test).
CPM_Chi2_New(x, N, m) IINew algorithm for controHer performance
monitor
else
}
RL++; Ilincrement the current (newest) RL
for k = 0 to N  1
R[k] =0; Ilinitialize the RL sequence to o.
for k =0 to m  1
y[k] = 0; lIinitialize RL histogram to 0
Ilbuild RL distribution for the 1st data window
Iistores sequence of RL within a window of size N
Ilmaximum number of runs within the window <= N.
int y[m]; lithe RL histogram, m is upper limit of RL considered
int R [N];
int RL = 1;
int IOldestRL = 0; lIthe index of oldest RL within the data window
IIR[loldestRLl is the leftmost RL in data window
int InewestRL= 0; lIthe index of the newest RL within the data window
IIR[lnewestRLl is the rightmost RL within window
R[lnewestRLl = 1; lithe initial value for any run length = 1
if x[k] crosses zero Ila zero crossing occurs.
y[RL]++; Ilincrement the histogram of current RL
R[lnewesIRL] = RL; Ilrecord the current RL
IneweslRL++; Ilincrement index for next (newest) RL
RL = 1; IIreset the initial value of new RL to 1
for k = 1 to N  1
{
Line {
1
2
34
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
Ilcomputing RL distribution after read initial N data, general case
int RLpartOld = R[loldestRL]; IIRLpartold is partial, oldest RL within window
28
22
23
24
25
26
27
28
29
30
31
booi runDropped =false;
bool runAdded = false;
for i =N to n  1 Iii is the sample time,
{
RLpartOld = RLpartOld  1; IIshift out one data from partial RL
if RLpartOkl =0; Iloidest partial RL =0, dropped from window
y[R[loldestRL]] ; Ildecrement oldest RL histogram
runDropped =true;
loldestRL =(loldestRL + 1) mod N; Ilnext oldest RL index
RLpartOld =R[loldestRLJ;
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46 }
47 }
Ilupdate newest RL after adding one newest data to window
if x[i] crosses zero lIa zero crossing occurs.
y[RL]++; Ilincrement the histogram of current RL
runAdded = true;
R[lnewestRd = RL; Ilrecord the current RL
InewestRL =(InewestRL + 1) mod N; Ilupdate the index for
Ilnext (newest) RL, wrap up to empty front space.
RL = 1; Ilreset the initial value of new RL to 1
else
RL++; lJincrement the current (newest) RL
if runDropped = true and runAdded = false Ilcase 1
calculate new x2 from old one using Equation (4.9);
else if runDropped = false and runAdded = true Ilcase 2
calculate new X2 from old one using Equation (4.10);
else if runDropped =true and runAdded =true Ilcase 2
calculate new l from old one using Equation (4.11);
... .".'
The number of executable statements (again, just counting the body of loops and
ignoring the loop head and statements outside the loops, assuming the worst case run
time for ifstatements) in the new algorithm is, approximately,
Tnew ::::::N + m + 4N + (n  N)* 16
= 16nllN+m
where n is the input data size, N the data window size, and m the upper limit ofRL
considered.
29
1.
The required memory space of the new algorithm (again, ignoring the temporary,
constant memory) is one array of integers of size N for storing the RL sequence within
the data window, and one arrays of integers of size m, for storing the RL histogram, i.e.,
Mnew = (N + m)*Sinb
where Sint is the memory size for one integer.
/
4.2.2 Recurrent Algorithm for KS Test Statistic Calculation
The run time for the KS test statistic calculation is more complex because the KS
test uses the cumulative histogram and therefore dropping or adding one run to the RL
histogram may change the cumulative histogram for many bins.
The pseudocode for the first case, where one run has dropped out of the data
window due to window shifting, is shown below:
if runDropped =true and runAdded =false lido statistical test
for k = RLctroppect to m Ilupdate cumulative histogram
cumHist[k] = cumHist[k]  1
if abs(cumHist[k]  cumRef[k]) > maxD1
IImaxD1 is max difference in cumHist in affected bins
maxD1 =abs(cumHist[k]  cumRef[k])
if maxD1 > S
liS is the max difference in cumHist, i.e., KS test statistic
S =maxD1 Ilupdate KS test result
else if maxD1 < Sand RLmaxD 2:: RLctropped
Ilprevious maxD bin frequency is affected by dropping a run
Ilunaffected bin should also be searched for maxD, i.e., S.
S = maxD1
for k = 1 to RLdropped  1
30
if abs(curnHjst[k] cumRef[k]) > S
S = abs(cumHist[k]  cum Ref[k])
runDropped =false
The pseudocode for the second case, where one run has been added to the data
window due to window shifting, is shown below:
if runAdded =true and runDropped =false lIdo statistical test
for k =RLadded to m Ilupdate cumulative histogram
cumHist[k] =cumHist[k] + 1
if abs(cumHist[k]  cumRef[k]) > maxD1
IImaxD1 is max difference in cumHist in affected bins
maxD1 = abs(cumHist[k]  cum Ref[k])
if maxD1 > S
lIS is the max difference in cumHist, i.e., KS test statistic
S = maxD1 Ilupdate KS test result
else if maxD1 < Sand RLmaxD 2: RLadded
Ilprevious maxD bin frequency is affected by adding a run
Ilunaffected bin should also be searched for maxD, i.e., S.
S = maxD1
for k = 1 to RLadded  1
if abs(cumHist[k]  cumRef[k]) > S
S = abs(cumHist[k]  cumRef[k])
runAdded = false
The pseudocode for the third case, where one run has dropped out of, and one run
has been added to, the data window due to window shifting, is shown below:
31

.. .... ~,
'" ".' ~.
if runDropped =true and runAdded =true lIdo statistical test
if RLdropped =RLadded Iisame run length is dropped and added
return Iinothing is changed
if RLdropped < RLadded
for k = RLdropped to RLadded  1 Ilupdate cumHistogram
cumHist[k] = cumHist[k]  1
if abs(cumHist[k]  cumRef[k]) > maxD1
IImaxD1 is maxDiff in cumHist in affected bins
maxD1 = abs(cumHist[k]  cumRef[k])
else if RLadded < RLdropped
for k = RLadded to RLdropped  1 Ilupdate cumHistogram
cumHist[k] = cumHist[k] + 1
if abs(cumHist[k]  cumRef[k]) > maxD1
IImaxD1 is maxDiff in cumHist in affected bins
maxD1 = abs(cumHist[k]  cumRef[k])
if maxD1 > 8
118 is the max difference in cumHist, i.e., K8 test statistic
8 = maxD1 Ilupdate KS test result
else if maxD1 < 8 and min(RLdrOPped, RLadded) :s: RLmClxO
and RLmaxD < max(RLadded, RLadded)
Ilprevious maxD bin is affected by run dropping & adding
Ilunaffected bin should also be searched for maxD, i.e., 8.
8 =maxD1
for k = 1 to min(RLdroPPed, RLadded)  1
if abs(cumHist[k]  cumRef[k]) > 8
8 =abs(cumHist[k]  cumRef[k])
for k =max(RLdropped, RLadded) to m
if abs(cumHist[k]  cumRef[k]) > 8
8 = abs(cumHist[k]  cumRef[k])
runDropped = false
runAdded =false
32
The runtime of the new algorithm depends on the code sections for performing
the statistical KS test, which have a worst case run time of Oem), and a best case nm
time of 0(1).
The best case happens when neither any run has dropped out of, nor any run has
been added to, the data window during window shift to include a new data sample, or
when exactly one run has been dropped and at the same time another run with the same
run length has been added to the data window. In the best case, the RL histogram has not
changed, so the test statistic is the same as the previous value, and there is no need to
perform the statistical test, and hence performing the statistical test has a constant time in
the best case.
The worst case happens when all m bins of the cumulative histogram need to be
visited to find the bin with the maximum difference between the current histogram and
reference histogram, i.e., the KS test statistic S. More specifically, i l happens when the
bin with the maximum difference between the previous RL histogram and the reference
histogram has a different histogram due to one run dropping out of, or one run being
added to, the data window or both, AND the maximum difference in the bins affected by
run dropping or adding is less than the previous statistic, i.e., the previous maximum
difference in all bins. In the worst case, all bins of the histogram need to he visited to
find the maximum difference because the maximum difference may be in one of the
unaffected bins, and hence the worst case run time for performing the statistical test is
linear, i.e., Oem).
33
There are cases between the best and the worst. Many times, not aJ I bins need to
be visited, and only part of the bins whose cumulative frequencies have changed need to
be visited. For example, if the frequency of the bin that has the maximum difference in
the previous cumulative histogram, called RLmaxD, is not affected by the run dropping or
adding due to the new data sample, then only the bins affected by run dropping or adding
(i.e., the bins with RL ~RLdropped or RL ~Ladded) need to be visited to find the
maximum difference, and there is no need to visit the remaining bins (i.e., the bins with
RL < RLdropped or RL < RLadded), which are not affected by run changes. Another
example, if the maximum difference in the affected bins, called maxD1, is larger than the
previous maximum difference in all bins, called maxD, then this maxD, is the maximum
difference in all bins no matter whether or not the bin with the previous maxD is affected
by run changes, and therefore there is no need to check the remaining bins, either, in this
case. The only case that needs to visit all bins, i.e., the worst case, is when the bin with
previous maxD is affected by run changes and the maxD1 value in the affected bins is less
than the previous maxD. All other cases should take less time than a linear sequential
traversal of all m bins of the cumulative histogram.
Although in terms ofbigOh notation, the run time for the statistical KS test is
still O(m), we should note that this is the worst case run time, and the actual, average case
nm time should fall between the constant time 0(1) and the linear time Oem).
The run time of the overall new algorithm using the KS test is O(n*m) (in the
worst case), which is an improvement over the old algorithm with O(n*(N+m)) because
N» m, in general.
34
A more detailed pseudocode of the new algorithm using the KS test is shown
below (See Appendix A.3 for the complete code of the new algorithm using the KS test).
CPM_KS_New(x, N, m) IINew algorithm for performance monitor
Line {
1
2
int R [N];
int y[m];
Ilstores sequence of RL within a window of size N
Ilmaximum number of runs within the window <= N.
lIthe RL histogram, m is Lipper limit of RL considered
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
for k = 0 to N  1
R[k] = 0; //initialize the RL sequence to O.
for k = 0 to m  1
y[k] =0; Ilinitialize RL histogram to a
int RL =1;
int loldeslRl = 0; lIthe index of oldest RL within the data window
IIR[loldeslRLl is the leftmost RL in data window
int InewestRl= 0; lithe index of the newest RL within the data window
IIR[lnewestRd is the rightmost RL within window
R[lnewes!Rd =1; lithe initial value for any run length =1
for k = 1 to N  1 Ilbuild RL distribution for the 1st data window
{
if x[k] crosses zero Iia zero crossing occurs.
yI[RL]++; Ilincrement the histogram of current RL
R[lnewestRl] =RL; I/record the current RL
IneweslRl++; Ilincrement index for next (newest) RL
RL = 1; Ilreset the initial value of new RL to 1
else
RL++; llincrement the cLlrrent (newest) RL
}
Ilcomputing RL distribution after read initial N data, general case
int RLpartOld = R[loldestRd; IIRLpartOld is partial, oldest RL within window
bool runDropped = false;
bool runAdded = false;
for i =N to n  1 IIi is the sample time,
{
RLpartOld = RLpartOld  1; Iishift out one data from partial RL
if RLpartOld = 0 Iioidest partial RL = 0, dropped from window
y[R[loldestRLl] ; Ildecrement oldest RL histogram
runDropped = true;
loldestRl = (IoldestRl + 1) mod N; Ilnext oldest RL index
35
_.
The number of executable statements (again, just counting the body of loops and
ignoring the loop head and statements outside the loops, assuming the worstcase runtime
31
32
33
34
35
36
37
38
39
40
41
42
43
44 }
45 }
RLpartOld = R{loldestRd;
Ilupdate newest RL after adding one newest data to window
if x[i] crosses zero lIa zero crossing occurs.
y[RL]++; Ilincrement the histogram of current RL
runAdded =true;
R[lneweslRd = RL; Ilrecord the current RL
InewestRL = (InewestRL + 1) mod N; Ilupdate the index for
Ilnext (newest) RL, wrap up to empty front space.
RL =: 1; Ilreset the initial value of new RL to 1
else
RL++; Ilincrement the current (newest) RL
if run Dropped ="true or runAdded =true
S = K_S_test(y); IIrefer previous 3case pseudocode
output test result S;
O(m) for the KS test) in the new algorithm is, approximately,
Tnew :::::: N + m + 4N + (n  N)*(12 + m)
= nm  Nm + 12n  8N .J_ m
where n is the input data size, N the data window size, and m the upper limit ofRL of the
KS test.
The required memory space of the new algorithm (again, ignoring the temporary,
constant memory) is one array of integers of size N for storing the RL sequence within
the data window, and two arrays of integers of size ill, one for storing the RL histogram,
and the other for the cumulative RL histogram, i.e.,
Mnew = (N + 2m)*Smt,
where Sint is the memory size for one integer.
36
.~...,.'.,,
..
"'
."" t.
4.3 Comparison of New Algorithms and Old Algorithm
4.3.1 New AJgorithm Using ChiSquare Tests
In terms of the number of executable statements, compared with the count for the
old algorithm, 4nN  4N2 + urn  Nm + n, the speedup ratio ofthe new algorithm using
the chisquare test to the old one, Rsc, is
R = Told = 4nN  4N2 + nm  Nm r n
sc T I6n lIN + m new
Since the number of input data, n » N, the window length, and N » m, the upper limit
of the RL for the KS test, the speedup is, approximately,
Compared with the required memory space of the old algorithm, Mold = N*Sreal +
m*Sinl (again, ignoring the temporary, constant memory), the memory ratio of the new
algorithm to the old one, in terms of memory req uirement, RMc, is
R = MTiew = eN \ m) * Sinl
Me M N*S +m*S old r",,1 inl
Since N » m and Sreal ~ Sint, the memory ratio is, approximately,
Since the memory required to store an integer number is usually less than or equal
to that for a real number, i.e., Sint/Sreal :::; 1, therefore, the memory ratio of the new
algorithm to the old one, is also less than or equal to 1.0, i.e., RMC :::; 1. That is, the new
algorithm requires no more memory than the old one. For a specific machine, if the same
37
storage is used for an integer and a real number (float or double), then the new algorithm
requires (almost) the same memory as the old one; if a real number requires twice as
much memory as an integer, then the new algorithm requires about half of the memory
required by the old algorithm, thus a close to 50% reduction in memory requirements.
Further reduction in memory requirement for the new algorithm is possible if the
machine has a type of integer with one byte. Because a run length is rarely larger than
laO in practice and is always positive, and one byte is enough to store a run length
between a and 255, we can save memory by using an array of Ibyte integers of size m
instead of the normal 4byte integers to store the RL sequence. However, note that the
reduction in memory is not on the order of magnitude, and the new algorithm just
requires no more memory than the old one.
4.3.2 New Algorithm Using KS Tests
The speedup ratio of the new algorithm, using the KS test, to the old one, RSK, is
4nN4N2 =+nmNm+n
nm  Nm + I2n  8N + m
Since the number of input data, n» N, the window length, and N» m, the upper limit
of the RL for the KS test, the speedup is, approximately,
4N RSK
::::;
m
Compared with the required memory space of the old algorithm, Mold = N*Sreal +
m*Sinl (again, ignoring the temporary, constant memory), the memory ratio of the new
algorithm to the old one, in tenns of memory requirement, RMK, is
38
R _ MileII' = (N+2m)*Silll
,14K M N*S *S
oM renl + m im
Since N»m and Sreal ~ Sint, the memory ratio is, approximately,
R
_ Sml
,14K 
Sren/
4.3.3 Runtime Comparisons
Table 2.1 shows the runtime comparisons of the old algorithm using chisquare
test and the two new algorithms (one uses the chisquare test and the other uses the KS
test) of the controller perfonnance monitor. The average runtimes for different input data
sizes and window sizes and the speedup ralios are listed. The input data are random
samples with a zero mean and unit variance. The runtimes are the average of 10 runs for
each condition. A PC with 95 megabytes of RAM, a 650 MHz CPU and the Windows 98
operating system is lIsed. The code is written in the Matlab language.
T able 2 1 RuntIme companson 0 f tbe New AII~on.tI1m and the Old AI120n'01m
Input data Data Old New algorithm
SIze window size algorithm
Using chis( uare tests Using KS tesLs
n ::J Average Average Speedup Average Speedup
Runtime, Runtime, ratio Runtime, ratio
sec sec sec
1,000 50 3.686 0.241 15.3 0.379 9.7
1,000 100 5.817 0.231 25.2 0.368 15.8
1,000 200 9.299 0.208 44.7 0.335 27.8
1,000 400 13.133 0.170 77.3 0.269 48.8
10,000 50 29.550 1.774 I 16.7 2.917 10.1
10,000 100 47.955 1.769 I 27.1 2.977 16.1
10,000 200 84.459 1.752 48.2 2.993 28.2
10,000 400 155.13 1.724 90.3 2.983 52.1
39
We can see that when the input data size n is fixed, as the window size N
increases, the nmtime of the old algoritlun increases in a linear rate, but the new
algorithm using the chisquare test does not change, indicating that its runtime is
independent of the window size N. Therefore, as window size n increases, the speedup
ratio is increase proportionally, which means that the larger the window size is the higher
the speedup ratio of the new algorithm using the chisquare test is. The speedup ratio of
the new algorithm using the chisquare test is independent of the input data size n. For
the nominal case with a window size of 100, the speedup ratio is above 20.
The runtime of the new algorithm using the KS test is also an order of magnitude
less than that of the old algorithm, but the improvement is not as large as the new
algorithm using the chisquare test. The runtime of neither of the two new algorithms is
dependent on the window size, while the runtime of the old algorithm is linear in the
window size n, which indicates an order of magnitude improvement. Similarly, the
higher the window size is, the higher speedup ratio ofthe new algorithm using the KS
test is. For the nominal case with a window size 0 f 100, the speedup ratio of the new
algorithm using the KS test is about 15, while the speedup ratio of the new algorithm
using the chisquare test is above 20.
We can see that the new algoritlun using the chisquare test is computationally
more efficient than the new algoritlun using the KS test, the larger the window size is,
the higher speedup is achieved by using the new algorithms, and both new algorithms are
much faster than the original algorithm.
40
;
::.'
.:5,
.'~.. i~
"
CHAPTERS
STATISTICAL MODELS FOR RL DISTRIBUTIONS
The controller perfonnance monitor uses a statistical goodnessoffit test to see
whether a1sample RL distribution is significantly different from a specific reference RL
distribution. The controller perfonnance monitor in (15] uses the actual data samples
collected during a representative good control period, judged by plant engineers, to build
the reference distribution. This reference distribution is totally databased and empirical,
and is only applicable for the specific controlled variable for which the data samples were
taken.
It would be useful to build theoretical statistical models for the reference RL
distributions, and investigate the relations between the RL distribution models and the
controller performance. Knowing the RL distribution models and their relations to the
controller performance, we may directly specify a reference RL distribution from a
theoretical model, rather than using an empirical RL distribution. Using databased,
empirical reference RL distributions requires finding a representative good control period
and collecting enough data samples for each controlled variable, which is not only time
consuming, but also requires expertise and experience in controller performance
judgment.
A model for the RL distribution of a random process with an independent,
identical distribution (iid) is established in [15]. However, no models of RL distributions
for correlated data have been found. There are many statistical models that deal with
runs, such as, waiting for the first run occurrence in Bernoulli trials (independent and
41
,......
~)
..,
.'..7.., ..'
identical binary trials), which results in a geometric distribution (see [22] for a review),
but no models are found dealing with data samples that are correlated.
This work derives theoretical statistical models for RL distributions for common,
practical, correlated processes, such as moving average (MA) processes, autoregressive
(AR) processes, and autoregressive and moving average (ARMA) processes. The
\
relations of the models to the controller performance will be investigated. First the
theoretical models for white noise are reviewed, and then the tlleoretical models for
correlated processes are derived.
5.1 RL Distribution Model for Wbite Noise Random Processes
A white noise random process is a sequence of random variables that are
independent and identically distributed. Gaussian noise is a white noise random process
with independent and identical Gaussian (or normal) probability density functions (pdf).
In the ideal, perfect control condition, where the controlled variable is always exactly at
the setpoint, the error signal (controlled variable measurements minus the setpoint)
contains only the measurement noise, which is commonly assumed to be Gaussian noise.
The theoretical statistical model for the RL distribution of a white noise has been
derived in [15J, and the derived RL probability mass function (or simply called the
probability function or distribution function) of a zero mean white noise is,
r...., .. )
~
....
dO ;r
"
r = 1,2, ... (5, I)
where R is a random variable representing a run length, r is any possible value of the run
length R, and Pr(R = r) is the probabi lity of having a run with run length equal to r.
42
Under ideal, perfect control condition, where the true value afthe controlled
variable is always exactly at the setpoint, the error signal (controlled variable
measurements minus the setpoint, i.e., e = CV· SP) is equal to the measurement noise,
which is commonly assumed to be a Gaussian white noise. So, after assuming that the
measurement noise is Gaussian noise, Equation (5.1) represents the theoretical model of
the RL distribution, i.e., the RL probability function, under ideal, perfect control
condition.
5.2 RL Distribution Model for Correlated Processes
The controlled variable measurements in practice are often correlated, so the
models that assume independent, identically distributed random variables are
inappropriate to represent the reference RL distribution, although this assumption is close
to being true in some cases.
The theoretical statistical models of RL distIibutions for a general, possibly autocorrelated,
process are derived below.
A run with a length of I has a sign pattern of either  +  or +  +, where ""
represents a negative value of the process output, and "+" represents a positive one.
Similarly, a run with a length of 2 has a sign pattern of either  + +  or +   r. In
general, a run with a length r has a sign pattern of eilher  (+y  or + (l +, where the
superscript r on (+Y means that there are r consecutive plus signs (+). A run with a sign
pattern of  (+Y  is called a positive side run, and + ()' + a negative run.
A zero crossing has a sign pattern of either  + or + . In general, a nm starts
with a zero crossing. A positive side run starts with a zero crossing with a sign
43
pattern +, and a negative side nm starts with a zero crossing with a sign pattern +.
If we use y(k) to represent the measured error signals at time k, we can represent
the event of having a positive side run starting at time k as
(y(k  I) < 0 It y(k) > 0 It y(k + 1) < 0)
where It represents a logical and operation. Therefore, the probability of a positive side
run with length 1 is
Since the probability functions of nm lengths are the probabilities for each nm
sign pattern of  +, therefore, the probability function [or a positive side run of length 1,
length among all run lengths, and a positive side run starts with a zero crossing with a
Pr(y(k  1) < 0 It y(k) > 0 It y(k + 1) < 0)
And the probability of a negative side run with length I is
Pr(y(k  1) > 0 It y(k) < 0 It y(k + 1) > 0)
(5.2)
(5.3) ~
r:r
)3
,)
.f.) )...
)
among all positive side runs, is
Pr(Rp = 1) = Pr(y(k  1) < a It y(k) > 0
.... ,...)
\
I
It y(k + 1) < 0 I y(k  I) < 0 It y(k) > 0) (5.4)
.o...J
which is a conditional probability, and the condition is the part after the sign I, i.e., y(k
I) < a () y(k) > 0, which represents a zero crossing with a sign pattem of  +.
Since the conditional probability has the following property,
Pr(AIB) = Pr(A () B)/Pr(B)
therefore, Pr(Rp = 1) in Equation (5.4) becomes
Pr(R
p
= I) = Pr(y(k 1) < Ony(k) > Ony(k + 1) < 0)
Pr(y(k 1) < Ony(k) > 0)
(5.5)
(5.6)
Similarly, the probability function for a negative side run of length I, among all
negative side runs, is
44
Pr(RN = 1) = Pr(y(k  1) > 0 n y(k) < 0
n y(k + 1) > 0 Iy(k  1) > 0 n y(k) < 0)
Pr(R
N
= 1) = Pr(y(k 1) > Ony(k) < Ony(k + 1) > 0)
Pr(y(k 1) > 0 n y(k) < 0)
(5.7)
(5.8)
Since a run must be either on the positive side or on the negative side, not both,
,
then the positive side runs and the negative ones are disjoint sets. And since the positive
side nll1S and the negative side runs always appear aitemately, one type following the
other type, then ignoring end effects the total number of positive side nms is equal to that
of negative side runs. Therefore, the total probability ofa run of length I among all run
lengths, i.e., the probability function, Pre R = 1), is:
Pr(R =1) = Pr(Rp = 1) +Pr(R N =1)
2
Pr(R = 1) = Pr(y(k 1) < Ony(k) > Ony(k+ 1) < 0)
2Pr(y(k 1) < 0 n y(k) > 0)
+ Pr(y(k 1) > 0 n y(k) < 0 n y(k + I) > 0)
2Pr(y(k 1) > 0 n y(k) < 0)
Similarly, the probability functions for run lengths or2 or more are
Pr(R = 2) = Pr(y(k 1) < Ony(k) > Ony(k+ 1) > On y(k + 2) < 0)
2Pr(y(k 1) < On y(k) > 0)
(5.9)
(5.10)
+ Pr(y(k 1) > Ony(k) < On y(k + 1) < Ony(k + 2» 0) (5.1l)
2Pr(y(kl) > Ony(k) < 0)
45
Pr(R = r) =
Pr(y(k 1) < On y(k) > 0" y(k + 1) > 0,,··· " y(k + r l) > 0 n y(k + r) < 0)
2Pr(y(k 1) < 0" y(k) > 0)
+ _PI'·C'y('k__l')_>_O_":y'(ek:..)<_O_"=:y'C_k_+:1)'<_0_,,_·._."=:...Y,(k_+_r_=1)_<_O_"=:...y'(k_+,r)_>'0)
2Pr(y(k 1) > O"y(k) < 0)
(5.12)
We can simplify the above functions. Without loss of generality, we can start
time counting at time 0, or, k  1 = 0 or k = I, which means that a zero crossing occurs
In addition, we are dealing with random variables with probability density
positive side runs are equal to those associated the negative ones, and therefore total
probability is just twice that associated with the positive side runs. Therefore, the RL
functions (pdf) that are symmetric to zero, so the probabilities associated with the
~
...
)'
or
)3)"
.. )
)'
.),
.~ ,
?
,
PR(r =1) =Pr(R =1) = Pr(y(O) < O"y(l) > 0" y(2) < 0) (5.13)
Pr(y(O) < 0" y(l) > 0)
probability functions become
between time k = 0 and time k = 1, and a new runs starts at time k = I.
PR(r =2) = Pr(R =2) = Pr(y(O) < O"y(l) > 0"y(2) > 0" )1(3) < 0) (5.14)
Pr(y(O) < 0" y(l) > 0)
PR(r) =Pr(R =r) =
Pr(y(O) < 0 n y(l) > 0" y(2) > 0"···,, y(r) > 0" yer + 1) < 0)
Pr(y(O) < 0 n y(l) > 0)
(5.15)
To derive the RL probability functions, we need to know r.ow the output y(k) at
different times k are related. A moving average process model is Llsed in the next section
to derive the RL distribution of correlated processes.
46
5.3 RL Distribution Probability Functions for Moving Average (MA) Processes
Moving average processes are important for controller perfoTIl1ancc monitoring
because the controlled variable measurements under minimum variance control (MVC),
on which many performance assessment techniques are based, is exactly an MA process,
and its order is equal to the plant delay.
An nthorder moving average process can be represented as follows:
where k is the discrete time, y the controlled variable measurements, n the order of the
MA process, bo, b l , .. _, b" the constant model parameters, and w(k) a zeromean white
y(k) = bow(k) + b)w(k  1) + ... + b"w(k  n)
noise random process.
5.3.1 RL Distribution Probability Functions for 151Order MA Processes
A firstorder MA process can be represented as,
y(k) = bow(k) + b1w(k  1)
The output sequence is
yeO) = bow(O) + b]w(I)
y(1) = bow(l) + b1w(O)
y(2) = bow(2) + bIw(l)
Note that the y(2) and yO) are correlated because both are dependent on w(2).
(5.16)
(5.17)
~..
)'
T)
J
)
I)
4
)...
),.' ,.)
\
J,.
~
To obtain the RL distribution for this correlated MA process, we need to calculate
the probability functions for different RLs.
47
For a firstorder MA process shown in Equation (5.17), we have
Pr(y(O) < 0 (1 y(l) > 0)
= Pr(bow(O) I b1w(1) < 0 (1 bow(l) I b]w(O) > 0)
= JJJfW(I)W(O)W(I) (w( 1), w(O), w(1))dw( l)dw(O)dw(l) (5.18)
bow(O)+b,w(  1)<0
bow(I)+b,w(O»O
where fW(_l)w(O)W(l)(w(l), w(O), w(l) is the joint pdfofthe three random variables w(I),
w(O) and w(1), and the integration region is over the 3D region that satisfies the two
variables w(1), w(O) and w(l):
Therefore, the joint pdf is the product of each individual pdf of each of the three random
inequalities, bow(O) I b,w(l) > 0 and bow(1) I b1w(0) < O. Since w(k) are white noise,
all w(k) are independent and have the identical pdf, fw(w)., i.e.,
~
•)
r
)
J
)
.'.) )
to
i
fW(k,(w(k)) = fw(w(k)) for all integer k (5.19)
fW(I)W(O)W(l)(w(I), w(O), w(1)) = fW(_I)(w(I)) fw(olw(O))fw(I)(w(l))
= fw(w(l))fw(w(O))fw(w(l)) (5.20)
Then Equation (5.18) becomes
Pr(y(O) < 0 (1 y(l) > 0) =
HJfw (w( I))dw(l)fw (w(O»dw(O)fw (w(1»)dw(l)
bow(O)rb,w( 1)<0
bow(l)+b,w(O»O
(5.21)
To satisfy the two inequalities, we have the integra] limits
b
W(O) <' weI)
bo
(5.22)
b
weI) > _I w(O)
bo
(5.23)
Then Equation (5.2]) becomes
48
Pr(y(O) < 0 n y(l) > 0) =
b,
r:fw (11'( 1 »dw(1) L~"(I)fw (w(O»dw(O) rb:W(O) .(,., (w(l»dw(1) (5.24)
bo
Similarly, we can obtain
Pr(y(O) > 0 n y(l) < 0) =
Ii r:fw (11'( l»dw(1) J~:w(I)~v (w(O»dw(O) L~'W(O) ~v (w(l»dw(l) (5.25)
"0
Because the pdf ofy(k) is symmetric around zero, we have
Pr(y(O) < 0 n y(l) > 0) = Pr(y(O) > 0 n y(l) < 0) =
b, r:fw (w( l»dw(1) r~~W(I) fw (w(O»dw(O) J~ w(O) ~v (w(l»dw(l) (5.26)
bo
Pr(y(O) < 0 n y(l) > 0 n y(2) < 0)
= Pr(bow(O) + b1w(I) < 0 n bow(l) + b]w(O) > 0
HJfw(J)w(0)W(,)1Y(2) (11'(1), 11'(0), w(1), w(2»dw(1)dw(O)dw(1)dw(2)
how(O)+b,,,,( I )<0
bowl 1)+bjw(O).>O
bo",(2)+b,,,,(I)<0
= HffW(J) (w(l»dw(1)fw(0) (w(0»dw(O)j~V(,)dw(l)fw(2)dw(2)
bow(O),·b,,,,( I)<(j
bowl J )·>b,w( 0).>0
"ow(2)+b,",(J)<O
~•
)r)!)
'..) )
to
)0
~
I
"
or, Pr(y(O) < 0 n y(l) > 0 n y(2) < 0)
Ii
= r:fw (11'( 1»dw(1) L~'W(I) fw (w(O»dw(O)
b,
J:; fw (w(1»dw(1) J"0",(1) fw (w(2»dw(2)
11'(0) '"
bo
49
(5.27)

Similarly, we can derive the following equations,
Pr(y(O) < 0 (l yO) > 0 (l y(2) > 0 (l y(3) < 0)
. ~
= r:fw(w(l)dw(l) L~ W(')fw(w(O))dw(O)
b
r~w(o)fw(w(l»dw(l) r~W(') fw (w(2»dw(2) L~lW(2)fw (w(3»dw(3) (5.28)
"0 "0
Pr(y(O) < 0 (l y(J.) > 0 (l y(2) > 0 (l ... n y(r) > 0 (l y(r + 1) < 0)
b
= r~fw (w(l»dw(l) f~'H'(l) fw (w(O»dw(O)
I: fw(w(1»dw(l)~: II .!,¥(w(2»dw(2)J... w(O) j 'w(l)
bo bQ
b
I:" fw (w(r»d'vv(r) ~".(r) fw (w(r + l»dw(r +1)
0",(rl) l~
bo
(5.29)
~
•)r)J
)
'..).
}
lt
The probability functions for all run lengths shown in Equations (5.13), (5.14),
and (5.15) can be represented by the integrals shown in Equations (5.24)  (5.29) for
computation.
To simplify the probabjlity function expressions, we define the following
jntegrals,
b,
ilO(x) = L~x fw (w)dw
iI2 (X) = r~/II(W)fw(w)dw
bo
50
(5.30)
(5.31)
(5.32)
Putting Equations (5.30) and (5.31) into Equation (5.26),
Pr(y(O) < a11 y(l) > 0)
= r:!w (w(l»dw(l) r:IV(_I/IO (w(O»!w (w(O»dw(O)
bo
From Equation (5.27), we have
Pr(y(O) < 011 y(l) > 0 11 y(2) < 0)
b
= r:fw (w( l»dw(1) r~IIV(I) fw (w(O»dw(O)
r:W(O) i JO (w(l»hv (w(l»dw(l)
bo
b
=r~°j;j! (w( l»dw(1) L~CIV(l)i, J (w(O»fw (w(O»dw(O)
From Equation (5.28), we have
Pr(y(O) < 0 11 y(l) > 0 11 y(2) > 0 11 y(3) < 0)
h
=r:fw (w(l)dw(1) t~'IV(J)fw (w(O»dw(O)
J::w(O) fw (w(l»dw(l) r:IV(,) i lO (w(2»fw (w(2»dw(2)
bo bo
b
= r:fw(w(l»dw(l) J~~'W(I)fw(w(O»dw(O)
f:; il J (w(l»fw (w(l»dw(l)
 w(O)
b"
51
(5.33)
(5.34)
(5.35)
~
•»r,
I)
,..
~
,".
..
!
"
Or,
Pr(y(O) < 0 n y(l) > 0 n y(2) > 0 n y(3) < 0)
61 = r:Iw (w(1»dw(1) L~ w(I) i I2 (w(On/w (w(O»dw(O)
Pr(y(O) < 0 n y(l) > 0 n y(2) > 0 n ... '1 y(r) > 0 n y(r + 1) < 0)
I;
= r~jw (w(1»dw(1) L~' 11'(1) ilrlw (w(O»dw(O)
(5.36)
(5.37)
Therefore, the RL probability functions for a 151order MA process, shown in
Equations (5.13)  (5.15), become
b, J+oolw (w(l»dw(1) JboW(I) ill (w(O»j;y (w(O»)dw(O)
PR (r = 1) = ·00 +00 <:/J
too ill ( w(1» 1w ( w(1»)dw(1)
b r:/w (w(l»dw(1) L~' wII) i l2 (w(O»/w (w(O»dw(O)
P (r = 2) = ~
R r:iIl(w(l»lw (w( Inc/w(1)
h,
J+oo/w (w(I»dw(1) Jb~W(I)i lr (w(O»fw (w(O»c/w(O)
P (r) = 00 <:/J
R J+<:/J
_00 ill ( w(1» 1w( w(1))dw(1)
·)... I
;t
(5.38) I•
,..
•I
(5.39)
(5.40)
To compute the values of the RL probability functions, PR(r), we need to assume
the pdf of the white noise w(k).
To derive the theoretical RL distributions, we assume the white noise is a
Gaussian white noise (with zero mean and unit variance). A Gaussian distribution is also
called a normal distribution, and its probability density function (pdf) is,
52
(x Ji)' I ,
IN (X;,u,cr) =,J"2; c cO'
2Jfcr
oo<x<oo (5.41)
where fl is its mean, and cr2 its variance. A nonnal distribution is represented as N(fl, cr\
N(O, (J2) is a zeromean nonnal distribution with fl = 0, and its pdf is,
oo<X<oo (5.42)
N(O, 1) is the standard normal distribution, whose mean is zero and whose
variance is one (fl =°and (J2 = 1). The pdf ofN(O, 1) is,
)
Since all w(k) are independent, and have the same pdfN(O, 1), shown in Equation
(5.43), then we have
_Xl lIN
(x;O,l) = ~ e 2
...;2Jf
oo<x<oo (5.43)
.,,..,
I
,
)
...
j,.
1 :w(i)2
Iw (w(i» = IN (w(i);O,l) = ~ e 2
...;2Jf
i = 0, 1, ... (5.44)
After specifying fw(w(k» using Equation (5.44), we can calculate the theoretical RL
probability functions for any 1storder MA process using Equations (5.38)  (5.40).
Numerical methods are used to calculate the integrations in these equations (see
Appendix A for details).
Table 5.1 lists the computed values of the theoretical RL probability functions for
some lSIorder MA processes with coefficient ratio bt/bo. See Appendix A.4 for the
complete code for calculating the RL probability functions for 151Order MA processes.
53
Table 5.1 Theoretical RL Distributions of 1sl0rder MA Processes
b,/bo Pr(l) Pl~2) Pr(3) Pr(4) PreS) Pr(6) Pr(7) Pr(8) Pr(9) Pr( 10) I
1 .6250 .2750 .0792 I .0173 .0030 .0005 I .000] .0000 .0000 .0000
0.8 .6225 .2735 .0808 .0187 .0037 .0006 .0001 .0000 .0000 .0000
0.6 .6127 I .2686 .0864 .0240 .0062 .0016 I .0004 .0001 .0000 .0000
0.5 .6038 .2650 .0908 .0283 .0085 .0025 .0007 .0002 .0001 .0000
0.4 .5916 .2610 .0962 .0336 .0116 .0040 .0014 .0005 .0002 .0001
0.2 .5549 .2534 .1095 .0470 .0201 .0086 .0037 .0016 .0007 .0003  r
0.1 .5297 .2509 .1171 .0546 .0254 .0119 .0055 .0026 .0012 .0006
0 ' .5000 .2500 .1250 .0625 .0312 .0156 I .0078 .0039 .0020 .0010
0.1 .4663 .2511 .1329 .0704 .0373 .0198 .0105 .0055 .0029 .0016
0.2 .4298 .2544 .1403 .0780 .0433 .0241 .0134 .0074 .0041 .0023
0.3 .3921 .2600 .1468 .0850 .0491 .0283 .0164 ,.0094 .0054 .0031
0.4 .3556 .2674 .1519 .0913 .0542 .0323 .0192 .0114 .0068 .0041
0.5 .3225 .2757 .1553 ' .0965 .0585 .0357 .0218 .0133 .0081 .0049
0.6 .2949 .2839 .1573 .1007 .0619 .0385 .0239 .0148 .0092 .0057 
0.7 .2740 .2910 . I581 .1037 .0643 .0405 .0254 .0160 .0]00 .0063
0.8 .2599 .2961 .1584 .1057 .0658 .0419 .0265 .0168 .0106 0067 I
0.9 .2523 .2991 .1584 .1068 : .0666 .0426 .0270 .0172 .0109 .0070
1 .2500 .3000 .1583 .1071 .0668 .0428 .0272 0.173 .0110 .0070
1.2 .2567 .2974 .1584 .1062 .0661 .0422 .0267 .0170 .0107 .0068
1.4 .2715 .2918 .1582 .1041 .0646 .0408 .0256 .0161 .0101 .0064
1.5 .2802 .2888 .1580 .1028 .0636 .0399 .0250 .0156 .0098 .0061
 1.6 .2890 .2858 .1576 .1015 .0626 .0391 .0243 .0151 .0094 .0059
1.8 .3064 .2803 .1566 .0989 .0605 .0373 .0230 .0]42 .0087 .0054
2 .3225 .2757 .1553 .0965 .0585 .0357 .0218 .0133 .0081 .0049
3 .3797 .2623 .1486 .0872 .0509 .0297 .0173 .0101 .0059 .0034
4 .4109 .2569 .1437 .0816 .0463 .0262 .0149 .0084 .0048 .0027
5 .4298 .2544 .1403 .0780 .0433 .0241 .0134 .0074 .0041 .0023
6 .4422 .2530 .1379 .0755 .0413 .0226 .0124 .0068 .0037 .0020
7 .4509 .2522 .1362 .0737 .0399 .0216 .0117 .0063 .0034 ! .0019 ,
8 .4574 .2517 .1348 .0723 .0388 .0208 .0112 .0060 .0032 .0017
9 .4624 .2513 .1337 .0713 .0380 .0202 .0108 .0057 .0031 .0016
10 .4663 .2511 .1329 .0704 .0373 .0198 .0105 .0055 .0029 .0016
15 .4779 .2505 .1303 .0678 .0353 .0184 .0095 .0050 .0026 .0013
"'
20 .4836 .2503 .1290 .0665 .0343 .0177 .0091 .0047 f.0. 024 .0012
30 .4892 .250] .1276 .0652 .0333 .0170 .0087 .0044 .0023 .0012
I
40 .4919 .2501 .1270 .0645 .0327 .0166 .0084 .0043 .0022 .0011
I
50 .4936 .2500 .1266 .0641 .0324 .0164 .0083 .0042 .0021 .0011
70 .4954 .2500 .1261 .0636 .0321 .0162 .0082 .0041 .002] .0010 
" 100 .4968 .2500 .1258 .0633 .0318 .0160 .0081 .0041 .0020 .0010
500 .4994 .2500 .1252 .0627 .0314 .0157 .0079 .0039 .0020 .0010 ..
1000 ' .4997 .2500 .1251 .0626 .0313.0157 .0078 .0039 .0020 .0010
S4
..)... ,i
,I
,
)
...
Computer simulations were perfonned to verify the theoretical values of the
derived RL probability functions. Figures 5.1 to 5.3 show the comparisons between the
theoretical RL distributions and sampIe distributions of some MA processes. We can
find that the sample RL distributions match the theoretical distribution well.
Run Length Histogram
r~,:~.,~ ' I Theoret~alll
. MA Sample=!  
300
250
>uc~
200
IT
i!?
LL
150
100
50
.}...
.~l
o _'__L~__,
1 234
,.~.'.
/ ,
'1;1 (.) r ,';,' 
5 6 7 8 9 10
Run Length, samples
Figure 5.1 Comparison of Theoretical RL Distribution and Sample RL Distribution
of a 1slOrder MA process with btfbo = 0.5.
55
Run Length Histogram
400 ,r,....,..,.;:=::=:c==.==:l::::=:::;1
350
Theoretical
MA Samples
300
>.
u
c
~ 200
0
~
u..
150
)
,'.'. I
'<'. ' .. 1 '
:;J
~~ J
(1)
.. ;
2 3 4 5 6 7 8 9 10
Run Length. samples
50
01L_'_'1
100
Figure 5.2 Comparison of TlIeoretical RL Distribution and Sample RL Distribution
of a 1510rder MA process with b1/bo = 1.
56
Run Length Histogram
::: f, .~,,
l~
300' i
250
>.
()
c
~ 200
CT
~
LL
150
~)
)
8 9 10
"" _., ,,
_1 __ .i_: 2_:_'J'
4 5 6 7
Run Length, samples
3
5°L o ~
2
100
Figure 5.3 Comparison of Theoretical RL Distribution and Sample RL Distribution
of a 1510rder MA process with b]/bo= 2.
5.3.2 RL Distribution Probability Functions for 2ndOrder MA Processes
The outputs of a 2ndorder moving average process shown in Equation (5.16) are:
yeO) = bow(O) + b lw(l) + b2w(2)
y(2) = bow(2) + b]w(1) + b2w(O)
y(k) = bow(k) + b,w(k  1) + b2w(k  2)
Then, compared with the 1storder case, we only need to change the integral limits and
the number of integral levels. Using the same procedure as before, we can derive the
following equations,
57
Pr(y(O) < 0 n y(l) > 0)
= IffI};fI(2)W(I»W(0)W(I)dw(2)dw(l)dw(O)dw(l)
bQw(O)+b,w( 1)+b1 w(2)<0
bo"'(I)+6,,,,(0)+b,,,,(1»0
= r:};¥ (w(2))dw(2) r:111' (w(1) )dw(1)
Pr(y(O) > 0 n y(1) < 0)
= r:'1w (w(2»dw(2) r:1w (w(l»dw(1)
Because the pdf of y(k) is symmetric around zero, we have
Pr(y(O) < 0 n y(l ) > 0) = Pr(y(O) > 0 n y( 1) < 0)
= r:111' (w(  2))dw( 2) r~1) f,y (w(  J»dw( 1)
 bl b,
r:\V(I)+ 6, w(2) 111' (w(O»dw(O) L~ w(O)+ .\;'",(1) 1w (w(l»dw(l)
b. bo
Now define the following integrals:
bl b,
i 20 (xl' x 2
) = t~'I+b;Xl 111' (w)dw
58
(5.45)
(5.46)
(5.47)
(5.48)
(5.49)
(5.50)
h.
Then Equation (5.47) becomes
Pr(y(O) < 0 n y(l) > 0) = Pr(y(O) > 0 n y(1) < 0)
= r:fw(w(2»dw(2) r:~v(w(l))dw(l)
r~W(_I)+::b)11'(2) i 20 (W(O), W( l»)fw (w(O))dw(O)
bo 00
= r:fw (W(2»dw(2) r:i21 (w( 1), w(2»fw (w(l»dw(1)
Similarly we can obtain
Pr(y(O) < 0 n y(l) > 0 n y(2) < 0)
= HIHfw(2)W(I)W(0)!l'(!)W(2)dw(2)dw(1)dw(0)dw(l)dw(2)
bow(O)+b,'>( 1)+b2w( 2)<0
bow(1 )+b, ",(0)+b2w( )>0
oow(2)+~w(l)+b) 11'(0)<0
= r:./;v (w( 2»dw(2) r:fw (w( l»)dw(l)
b, b,
Jb;W(')'b;W(2) Ill' (w(O))dw(O) [:, b
1
fw (w(l)dw(l)
00 w(O)+·w(I)
Va bo
b v
t~1 w(l)+~\V(O) fw (w(2»)dw(2)
'' r</~ f w(w(2»dw(2) r:./;v (w(l»c/w(I)
(5.51)
(5.52) •..
b, V2 L~W(')+boW(2) Iw (w(O»dw(O) f~b, 11'(0)+ b
1
11'(1) 120 (w(l), w(O»lw (w(l ))dw(l)
bo b"
=; r:fw(w(2»dw(2) r:fw(w(l»dw(l)
59
bl b,
t~ ",(1)+,;;",(2) i
21
(W(O), w(I »J;v (w(O»dw(O)
Pr(y(O) < 0 (l y(1) > 0 (l y(2) > 0 (l y(3) < 0)
= r:"fw(w(2»dw(2) r:fw(w(l»dw(J)
h b
L~'IV('}+~IV(2)in (w(O), w(l»J;jI (w(O»dw(O)
Pr(y(0) < 0 (l y(l) > 0 (l y(2) > 0 (l .. , (l y(r) > 0 (l y( r + 1) < 0)
= r:fw (w(2»dw(2) r:J;y (w(l»dw(1)
b, b,
L~ W(Il+ bo ·W
(2) i
2r
(w(O), w(l»fw (w(O»dw(O)
(5.53)
(5.54)
(5.55)
Then, from Equations (5.13)  (5.15) and Equations (5.52)  (5.55), the RL probability
functions for a 2ndorder MA process become
hi h1 r~fw (w( 2»dw(2) r:fw (w(l»dw(1) r~W(I)+~",(2) in (w(O), w(l»fw (w(O»dw(O)
r:fw (w(2»dw(2) r~i21 (w( 1), w(2»fw (w( l)dw( 1)
(5.56)
r: b, b, fw (w(2»dw(2) r:fw (w(l»dw(I) L::W(')''b;w(2} in (w(O), w(l»fw (w(O»dw(O) r:fw (w( 2»dw(2) r:i 2J (w( 1), w(2»fw (w(l»dw(1)
(5.57)
60
~ b, r:Ill' (W(2)dw(2) r<Y~Iw (w(l»dw(1) L~W(')+b;;0Il'(2)i2r (w(O), w(1»lw (w(O»dw(O)
r:Iw (w(2»dw(2) r:i21 (w(1), w(2»/11' (w(l»)dw(1)
(5.58)
Figure 5.4 shows the theoretical RL distributions and sample distributions of a
2ndorder MA process. We can see that the sample RL distribution matches the
theoretical distribution well.
Run Length Histogram
~. 0.
3 4 5 6 7 8 9 10
L_' ~__j_____l._ _L.. ' . __ ,•.•
400
350
300 
250t
>.
0c
~ 200
g
~
l.L
150
100
50
0
1 2
Run Length, samples
Figure 5.4 Comparison of TheoreticaJ RL Distribution and Sample RL Distribution
of a 2ndOrder MA process with b)/bo = 1 and b2/bo = 1.
61
5.3.3 RL Distribution Probability Functions for HighOrder MA Processes
The outputs of an nthorder moving average process shown in Equation (5.16)
are:
y(l) = bow(1) + b1w(0) + b2w(1) + ... + bnw(l  n)
y(k) = bow(k) + b1w(k 1) + b2w(k  2) + ... + bnw(k n)
Then, compared with the 1s'order case, we only need to change the integral limits and
the number of integral levels. Using the same procedure as before, we can derive the
following equations:
Pr(y(O) < 0 n y(l) > 0)
Jff .. Jfw(nJw(2Pllll'(oJw(I)dw(n)dw(1 n)···dw(O)dw(l)
00"'(0)+0,,,( I )+021"/( 2)+·· ·+o"w( PI )<0
bow( 1)+ I>, w( 0) +b1w( I )+.. ·+b" w(J PI»0
= r~fw (w( n»dw(11) [~fw (w(l n»dw(l  n) r:···r:fw(w( l»dw(1)
(5.59)
Pr(y(O) > 0 n y(l) < 0) = Pr(y(O) < 0 n y(1) > 0)
=r:fw (w(n»dw(n) r:fw (w(1 n»)dw(l n) r:··· [~fw (w(l»dw(1)
(5.60)
61
Pr(y(O) < 0 (l y(l) > 0 (l y(2) < 0)
= JH· JfWlII)WIIn)WII)1l'12)dw(n)dw(1  n)· ··dw(1)dw(2)
bow(O)+b,,,,( I )+b,w(2)+· ··+bnw( 11)<0
bOlO' I )+b,"'(O)+b,w(I )+.. ·+b" w(I1I»0
bowl 2)+b, \I'(I)+I.>,w(O)+.· ·+b"w(2n)<0
=r:fw (W( n»dw(/1) r:fw (w(l n»)dw(w(l n» r:···r:lw (W(l»dw(1)
Pr(y(O) < 0 (l y(l) > 0 (\ y(2) > 0 (\ y(3) < 0)
(5.61)
=r:./;y (w(n»dw(n) r~l)J;y (w(l n»)dw(w(1 n» r:···r:fw (w(l»dw(1)
(5.61)
Pr(y(O) < 0 (l y(l) > 0 (l y(2) > 0 (l ... (\ y(r) > 0 (\ Y(T + 1) < 0)
= r:lw (w(n»dw(n) r:Iw(w(1 n»)dw(w{1 n» r:~··· r~'lw (w( l»dw(1)
b, b, I.>"
J~"'(I)+~W(2)+'+;;;;W(II) Iw (w(O)dw(O) r: I.> b Iw (w(l»dw(l)
<Y) .!wlO)+.lw(I)+···+'!. ...(III)
bo 00 b"
63
b, b, b
1~",(r)b~\I(rI)..t."\I'(rn+l) fw (w(r + 1»dw(r + 1)
Now define the following integrals:
Then Equations (5.59)  (5.62) become
Pr(y(O) < 0 n y(l) > 0) = Pr(y(O) > 0 n yO) < 0)
(5.62)
(5.63)
(5.64)
(5.65)
(5,66)
= r:fw (w(n»dw(n) r:fw (w(1n»dw(1  n) r:··r:.~¥ (w(·2»dw(2)
r:in' (w( 1), w(2),··" w(n»fw (w(l»dw(I)
Pr(y(O) < 0 n y(1) > 0 n y(2) < 0)
(5.67)
= L:fw (w(n»dw(n) r:fw (w(l n»dw(w(l n» r:···r~J)fw (w(1»dw(I)
b ~ b
L~~' w(I)+ &;w(2)++ ,,:w(n) i
nl
(w(I), w(2)," " w(n»fw (w(O»dw(O) (5.68)
Pr(y(O) < 0 n y(l) > 0 n y(2) > 0 n y(3) < 0)
= r:};1/ (w( n»dw( n) r~f w(w( I  n) )dw( w(l  n» r:...r:~fw (w(1» dw( 1 )
64
 b. n, b
L~~' 1I(ll/'OW(2l+ +b~W(II)i
ll2
(w(I), w(2),·", w(n))fw (w(O»dw(O) (5.69)
Pr(y(O) < 0 (j y( I) > 0 (j y(2) > 0 (j ... (j y(r) > 0 (j y(r + 1) < 0)
= L~jw (w(n)}dw(n) r:fw (w(1 n»)dw(w(l n» r~ ..·r:./;v (w(I»dw(1)
t:&1 b1 b lV(I)+~lV(2l~ . + bn" w( II) i"r (w{1), w( 2), .. ·, w(n»fw (w{O»dw(O) (5.70)
Using Equations (5.13)  (5.15) and Equations (5.67)· (5.70), we can calculate the RL
probability functions for aU run lengths for any order MA processes.
65
CHAPTER 6
Demonstration of Controller Performance Monitor
The controller performance monitor is demonstrated using the experimental data
from reference [15].
6.1 Description of Experimental Unit and Procedure
The experimental unit is illustrated in Figure 6.1. Two air streams (large and
small) and one water stream are mixed and then enter a vertical tube (2.5 cm nominal
diameter,S m tall), where the water phase and the air phase interact and flow upward.
The three streams are controlled independently by three independent PI controllers with
the CamileTG 2000 hardware and software (version 4.0) for automatic data acquisition
and control (DAC). Orifice differential pressure transducers transmit 420 rnA signals to
the DAC where software calibrations reconstmcl flow rates, and send 420 rnA signals to
the i/p devices which provide instrument air for the flow control valve actuators.
66
Airwater twophase flow
column
:~~
FT
3
:~~
FT
1
Water flow
Large air flow
SmaJl air flow
Figure 6.1 TwoPhase Flow Experimental Unit with ODe Water Flow Control Loop
and Two Air Flow Control Loops.
These goodnessofcontrol experiments are conducted on the water flow control
loop. The process exhibits firstorder plus deadtime dynamics, modest valve stiction and
noticeable nonlinear characteristics. When valve stick happens, the valve does not move
until the change in input signal to the valve is larger than a significant threshold value.
67
6.2 Experimental Procedure
Tune the PI controller at the nominal operating point. After tuning, observe the
closedloop settling time from a setpoint change. Collect the CV and SP data during a
representative period of good control. From the data, build the reference RL distribution,
and determine the sampling window length, N, and the grace period (if a grace period is
desired). Specify the level of significance or confidence coefficient for the chisquare
test. Start to run the perfonnance monitor. Observe whether the performance monitor
can detect poor control performance.
6.3 Results and Discussion
Figure 6.2 shows the sample reference RL distribution built from a representative
period of good control with 1000 data samples and 295 RL. Table 6.1 shows how the
bins of the RL distribution are grouped into classes for the chisquare test, and the
theoretical RL distribution probabilities of a 151order MA process, and a 2ndorder MA
process. We can see that the theoretical RL probabi liti es of selected MA processes are
close to the sample RL probability distribution. We can use the MA process model
parameters, b ,Ibn for Islorder processes, or b)lboand b2/bo for 2ndorder processes, to
characterize the reference RL distribution. The advantage is that the MA process model
parameters indicate the autocorrelations among the process data. In an ideal, perfect
control case, the process output should not exhibit any autocorrelation.
68
Reference RL Distribution, 1000 CV Samples (295 RL Observations)
120 I I i 
,.1 _c:::L.._ _ _
10 15 20 25 30
RL, samplings
5
(/)
c
.Qm~
Q)
(/)
.0 o
....J
0::: o
'
Q)
.0
E
::J
Z
Figure 6.2 Sample Reference RL Distribution for Water Flow Control Loop. Use
1000 data samples with 295 RL observations. (Sampling period Ts = O. J second)
Table 6.1 Reference RL Distribution Bin Groups, Water Flow Control
I
Class # I 1 2 J I 4 )
 I
RL, samplings 1 2 3 4,5,6 7,8, .. ,
Sample Probability 0.386 0.180 0.146 0.156 0.132
r
Theoretical Probability of 1slorder 0.3846 0.2613 0.1479 0.1657 0.0405
MA with b,/bo = 0.32
Theoretical Probabi lity of 2ndorder 0.3810 0.] 818 0.1586 0.2029 0.0757
MA with b1/bo = 0.5, bl/bo = 0.3
Theoretical Probability of 2ndorder 0.3901 0.1794 0.1564 0.2003 0.0738
MA with bl/bo = 0.48, b1/bo = 0.3
69
The data window size is N = 129 samplings. Select the chisquare test
significance level (onesided) of ex = 0.003. Since there are k == 5 classes of RL values
and no distribution parameters are estimated from the data samples, the number of degree
of freedom of the chisquare statistic is 5  1 = 4. The critical value i u.kJ is 16.17.
Figure 6.3 shows the experimental data, CV, SP, MV, the chisquare statistic, and
the flagging when the sample reference RL distribution is used. From Figure 6.3, it can
be seen that when the plant operates near the nominal operating point with a water flow
rate of 35 kgfhr, where the controller is tuned, the monitor does not flag even though
there are setpoint changes. As the water flow rate setpoint decreases and moves away
from the nominal operating point, control performance deteriorates. The monitor starts
flagging when the flow rate setpoint is reduced to 15 kg/hr from 20 kg/hr, when poor
control performance with excessive oscillations occurs. The monitor flags most of the
time when the water flow rate is less than 15 kg/hr and the controller performs poorly.
After the flow rate moves back to the region near the nominal operating point and the
controller performance becomes good again, the monitor stops nagging.
70
...................
Experimental Data (CV, SP, MV), Chi2 Statistic, and Flagging
40 (.,
(f)Go.._ 20 f""" '.' c:· :..<~"'"'C ;~
~. ~~~~'
o1' 1.._·;;_.·~..'_'1',.,'':;:."'.:...•.:..''_:__......_' .J
a 5000 10000 15000
5000
lime, sample
10000
10000
15000
IJ
15000
Figure 6.3 Experimental Data, CbiSquare Test Statistic, and Flagging Wben a
Sample RL Distribution is Used as Reference Distribution.
When the theoretical RL distribution of a 2ndorder MA process with b I/bo= 0.5
and b2fbo= 0.3 is used as the reference RL distribution, the chisquare test statistic and
the flagging are shown in Figure 6.4.
71
Experimental Data (CV, SP, MV), Chi 2 Statistic, and Flagging
i), ::r'~"~,'. '. ',._.,,_. :~1
> ~~.:>w ;,........u...~
U . 7"';..;. ~'........\ ...~. '~~... 7,:'" :..... ,..;...
O ~
'"''.:c.:'.'__ L _
o 5000 10000 15000
15000
10000 15000
5000 10000
5000
1r..r~' ..:,'•..,;...c
~ 0: l' """"" ,.;\•••;>j.<;{.~.'!'i.';':;'.::;".)';/"")"')"'r.".,_i:_,:J_""'· ".'A"""~j
o
ro 200 , ,r _
U I 'I 'C
u , ,
N 100 . r '. i N~ 0 .',_ ••" ••J"J , ..',v" \';';'''.!.j! .:.•"'.!'.·;i~·.'Y;"L •...,......:o~:,J
o 0
1
I ..
g' 0.5
0:::
; i
o ~
o 5000 10000 15000
tirre, sample
Figure 6.4 Experimental Data, ChiSquare Statistic, and Flagging When a
Theoretical RL Distribution (of a 2ndOrder MA Process with Parameters b)/bo = 0.48
and b2/bo = 0.3) is Used as Reference Distribution.
From Figure 6.4, we can see that using a theoretical reference RL distribution
introduces more flagging than a sample reference RL distribution. It is because the
theoretical RL distribution does not exactly match the real distribution, and therefore
introduces errors and shorttenn flagging. To reduce this effect, we can use a grace
period, and flag only when the chisquare statistic exceeds its critical value continuously
for a period longer than the grace period. The flagging with a grace period of a window
72
length (129 samples) is shown in Figure 6.5. We can see that using a grace period
prevents shortterm flagging.
Experimental Data (CV. SP, MV), Chi2 Statistic, and Flagging
15000
15000
15000
.:•..•,::•••v •. ;",._•.•::1

10000
10000
10000
:.'.
,, ..
5000
5000
5000
s..:.>;,~
.. '.
o ::0..,' .~.h......,;':.,_r. F'J....}.)\;~}\\,j)
o
1
40 ~.~
0
(/) 20
>
0
0
0
N
.c o
o
<0 200..,, _
u
N
B
100 I r=
o
1l r 
'I I ! :. ,Ii I I I 1" ! Ii!
. 'I I .
: , I o i· 'J'" I . I ~ ._;,l.J~ .J _"_
5000 '0000 15000
time, sample
Figure 6.5 Experimental Data, ChiSquare Statistic, aDd Flagging When a Grace
Period is Used Along with a Theoretical Reference RL Distribution (of a 2ndOrder
MA Process with Parameters bJ/bo = 0.48 and b2/bo=0.3). The Grace Period is Equal to
the Window Length (129 samples).
The run time of the new algorithm using the chisquare test for this input data size
of 15,000 and a window size of 129 is 4.45 seconds, while the old algorithm takes 118.20
seconds. The speedup ratio of the new algorithm to the old one is 26.56.
73
CHAPTER 7
CONCLUSIONS AND RECOMMENDATIONS
Conclusions
(1) An efficient new computer algorithm has been created for monitoring (process)
controller performance based on RL distributions. The algorithm analysis shows
that the new algorithm reduces the runtime complexity significantly, from O(n*(N
+ m)) to O(n), compared with the original algorithm, and the memory required is
slightly less (where n is the input data size, N the data window size, and ill the
maximum run length considered).
(2) Statistical models ofRL distributions for correlated data were built, and the
theoretical RL probability functions for any order MA processes were derived.
Simulations show that the derived theoretical distributions match those built from
simulated data.
(3) The statistical chisquare goodnessof fit test method provides more testing power
and uses less computation time than the KS test.
Recommendations
The relationships between the model parameters ofMA processes and the
controller perfonnance are worth further investigation so that we can directly specify the
reference RL distribution by detennining the MA process parameters accord ing to the
desired controller perfonnance.
74
BIBLIOGRAPHY
]. Ogunnaike, B.A., Ray, W.H. Process Dynamics, Modeling, and Control. Oxford
University Press: New York, 1994.
2. Marlin, T.E. Process Control  Designing Processes and Control Systems for
Dynamic Performance. McGrawHill: New York, 1995.
3. Harris, TJ. Assessment of control loop performance. Canadian Journal of
Chemical Engineering 1989; 67:856861.
4. Harris, TJ., Seppala, c.T., Desborough, L.D. A review of perfOlmance
monitoring and assessment techniques for univariate and multivariate control
systems. Journal ofProcess Control 1999; 9(1): 117.
5. Qin, SJ. Control perfOlmance monitoring  a review and assessment. Computers
& Chemical Engineering 1999; 23(2):] 73186.
6. Huang, B., Shah, S.L. Performance Assessment ofControl Loops: Theory and
Applications. Springer Verlag: New York, ] 999.
7. Rhinehart, R.R. A watch dog for controller performance monitoring, in
Proceedings ofthe 1995 American Control Conference: Seattle, WA, 1995; 22392240.
8. Venkataramanan, G., Shukla, V., Saini, R., Rhinehart, R.R. An automated online
monitor of control system performance, in Proceedings ofthe American Control
Conference: Albuquerque, NM, 1997; 13551359.
9. Narayanaswamy, G. Refinement and experimental verification ofa statistica!hased
goodnessofcontrol monitor. M.S. Thesis. Department of Chemical
Engineering. Texas Tech University: Lubbock, Texas, 1998.
10. Hagglund, T. A controlloop performance monitor. Control Engineering Practice
1995; 3(1):1543]551.
11. Thornhill, N.F., Hagglund, T. Detection and diagnosis of osci l1ation in conlrol
loops. Control Engineering Practice 1997; 5(10):13431354.
12. Miao, T., Seborg, D.E. Automatic detection of excessively oscillatory feedback
control loops, in Proceedings of 1999 IEEE International Conference on Conlro!
Applications: Hawaii, 1999; 359364.
13. Thornhill, N.F., Huang, B., Zhang, H. Detection of multiple oscillations in control
loops. Journal ofProcess Control 2003; 13:91100.
75
14. Li, Q., Whiteley, J.R., Rhinehart, R.R. A relative perfol1nance monitor for process
controLlers. International Journal ofAdaptive Control and Signal Processing
2003; 17(9):685708.
15. Li, Q., Whiteley, l.R., Rhinehart, R.R. An automated perfonnance monitor for
process controllers. Control Engineering Practice 2003 :in press, accepted in June
2003.
16. D'Agostino, R.B., Stephens, M.A. GoodnessOlFit Techniques. Marcel Dekker:
New York, 1986.
17. Kotz, S., lOMson, N.L. Encyclopedia ofStatistical Sciences. Jolm Wiley & Sons,
1988.
18. Stephens, M.A. Tests based on EDF statistics, in GoodnessOlFit Techniques,
D'Agostino, R.B., et aI., Editors. Marcel Dekker: New York, 1986; 97193.
19. Moore, D.S. Tests of ChiSquared Type, in GoodnessOlFit Techniques,
D'Agostino, R.B., et aI., Editors. Marcel Dekker: New York, 1986; 6395.
20. Pettitt, A.N., Stephens, M.A. The KolmogorovSmimov goodnessoffit statistic
with discrete and grouped data Technometrics 1977; 19:205210.
21. NIST/SEMATECH. Engineering Statistics. eHandbook ofStatistical Methods.
NIST/SEMATECH: b.!.!.R:l1www.itl.nisl.gov/di v898/handbook/, date of latest
access: June 27, 2003.
22. Balakrishnan, N., Koutras, M.V. Runs and SCClns with Applications. John Wiley &
Sons, Inc., 2002.
76
APPENDICES
A.I Code for tbe Old Algorithm of the Controller Performance Monitor
% CPM_Chi2_0Id.m, Old algorithm of Controller Performance Monjtor
% by Qing Li, Last modified: 11/23/2003
% inputs: x is a vector containing all data (actuating errors)
% N is the window size, # of data within data window
% maxRL is the maximum RL considered
% output: y is a vector containing a sequence of chisquare test
% stacistic resuit at each sample time
funct~on y = CPM Chi2 Old(x, N, maxRL)
tic;
w zeros (N, 1);
w x(::N);
n length (x) ;
y zeros (n, 1);
for k = (N + 1) : n
w = [w(2:end); x(k)];
hRL = zeros(l, maxRL);
sPrev = sign(w(l); %sign of previous data
if sPrev == 0
sPrev = 1;
end
runLen = 1;
for i = 2:N
s = sign (w(i» ;
if s = 0 & s = sPrev
if runLen <= maxRL
hRL(runLen) = hRL(runLen) + 1;
else
hRL(maxRL) = hRL(maxRL) + 1;
end
sPrev = s;
runLen I;
else
runLen runLen + 1;
end
end
% goodnessoTftttest(hRL);
y(k) = chi2Test(hRL);
end
t = toe;
global executionTimeChi20ld;
executionTimeChi201d = executionTimeChi20ld + t;
function y = chi2Test(hRL)
% reference,distributlon
binWidth = [1, 1, 2, 3];
binNum = length(binWidth) + 1;
refHist = zeros(l, binNum);
77
icum = 1;
for i = l:binNum  1
for j = l:binWidth(i)
refHist(i) = refHist(i) + 2"(icum);
icum = icum + 1;
end
end
refHist(end) = 1  sum(refHist(l:end  1»);
~ current histogra~
curHist = zeros(l, binNum);
icum = 1;
for i = l:binNum  1
for j = l:binWidth(i)
curHist(i) = curHist(i) + hRL(icum);
icum = icum + 1;
end
end
totalRLs = sum(hRL);
curHist(end) = totalRLs  sum(curHist(l:end  1);
chi2Statistic = 0;
for i = l:binNum
if refHist(i) = 0
chi2Statistic chi2Statistic ...
+ (curHist(i)  refHist(i) *totalRLs)"2
/(refHist(i)*totalRLs) ;
else
disp('Warning in criisquare test. ,);
end
end
chi2Critical = 9.488; ~critical value
hRL (1: 7) ;
[refHist*sum(curHist); curHist];
y = sign(chi2Statistic  chi2Critical);
78
A.2 Code for the New Algorithm Using ChiSquare Test
% CPM Chi2 ~e~_{n  
% Ne''I algorithm for Concroller Perfolmance IVlonitor using C11i2 test
% by Oin9 Li, Last modified: 11/~3!2003
% inputs: x is a vector cOltaining all data (actuating erlors)
% N is the window size, ff of data within data window
% maxRL is the max~mum RL considered
~ output: y is a vector containing a sequence of chi square test
% statistic at each sample time
function y = CPM_Chi2_New(x, N, maxRL)
tic;
RLS = zeros(N, 1); %RL sequence in window
hRL = zeros(l, maxRL);
runLen = 1;
InewestRL = 1; ~index of newest run in window
RLS(InewestRL) = 1; ~initial value for new ~un
% build RL distribution for 1st data window
sPrev = sign(x(I); %sign of previous data
if sPrev == 0
sPrev = 1;
end
for k = 2:N
s = sign (x (k) )
if s = 0 & S = sPrev '~runLe'1 is true RI~
if runLen <= maxRL
hRL(runLen) = hRL(runLen) + 1;
else
hRL(maxRL) = hRL(maxRL) + 1; ~but count as maxR~
end
RLS(InewestRL) = runLen;
InewestRL = InewestRL + 1;
sPrev = s;
runLen 1;
else
runLen runLen + 1;
end
end
% chi2TestlhRL)
% reference distribution
refDist = zeros(l, maxRL);
refDist = O,5,A(I:maxRL);
refDist(end) = 1  sum(refDist(l:end  1));
binWidth [I, 1, 2, 3];
binNum = length (binWidth) + 1;
~refHist is binGrouped h~st, wh:le refDist is not grouped.
refHist = zeros(l, binNum);
classBin = binNum*ones(maxRL, 1); '"indey. is PL, content lS class bin ;;
icum = 1;
for i = l:binNum  1
for j = ~:binWidth(i)
refHist(i) = refHist(i) + refDist(icuml
classBin(icum) = i;
79
icum icum + 1;
end
end
refHist(end) = 1  sum(refHist(l:end  I});
~; cun:ent histogram
curHist = zeros(l, binNum);
icum = 1;
for i = l:binNum  1
for j = l:binWidth(i)
curHist(i} = curHist(i} + hRL(icum);
icum = icum + 1;
end
end
totalRLs = sum(hRL);
curHist(end) = totalRLs  sum(curHist(l:end  1));
chi2 = 0;
for i = l:binNum
if refHist(i) = 0
chi2 = chi2 + (curHist(i)  refHist(i}*totalRLs)A 2 ...
I (refHist (i) *totalRLs) ;
else
disp('Warning in chi·square test. ,);
end
end
~build RL distributions tor subsequent windows
IoldestRL = 1; ~index of oldest run in window
RLToDrop = RLS(IoldestRL);
RLpartOld = RLTo~rop;
runDropped = 0;
runAdded = 0;
n = length (x) ;
y = zeros (n, 1);
for k = (N + 1) : n
RLpartOld = RLpartOld  1; ~shj ft au 1 data
if RLpartOld == 0 ~oldest RL dropped out
hRL(min(maxRL, RLToDrop) = r.RL(min(maxRL, RLToDrop))  1;
curHist(classSin(min(maxRL, RLToDrop))) =
curHist(classBin(min(maxRL, RLToDrop))}  1;
runDropped = 1;
totalRLs = totalRLs  1;
IoldestRL = IoldestRL + 1;
if. IoldestRL > N
IoldestRL = IoldestRL  N; %nex~ oldest
end
RLDropped = RLToDrop;
RLToDrop = RLS(IoldestRL);
RLpartOld = RLToDrop;
end
~update newest RL after adding one new dat3
s = sign (IX. (k) ) ;
if s = 0 & s = sPrev
if runLen <= maxRL
hRL(runLen) = hRL(runLen) + 1;
80
eise
hRL(maxRL) = hRL(maxRL) + 1;
end
curHist{classBin(min(maxRL, runLen»)
curHist(classBin(min{maxRL, runLen»} + 1;
totalRLs = totalRLs + 1;
RLS(InewestRL) = runLen;
RLAdded = runLen;
InewestRL = InewestRL + 1;
if InewestRL > N
InewestRL = InewestRL  N;
end
sPrev = s;
runLen = 1;
runAdded = 1;
else
runLen runLen + 1;
end
% chisquare goodnessoffittest;
if runDropped == 1 & runAdded == 0 ;,case 1
nRk = totalRLs + 1;
j = classBin{minlmaxRL, RLDropped);
OJ = curHist(j) + 1;
pj = refHist(j};
chi2 = chi2*nRk!(nRk  1) + 12*Oj + 1)!(nRk  l)!pj
+ 12*nRk  I)! (nRk  1);
runDropped = 0;
elseif runDropped == 0 & runAdded == 1 %case 2
nRk = totalRLs  1;
j = classBin(min(maxRL, RLAdded);
OJ = curHistlj)  1;
pj = refHistlj);
chi2 = chi2*nRk! (nRk + 1) + {2*Oj + l)!(nRk + l)!pj
+ (2*nRk  1)!(nRk + 1)
runAdded = 0;
elseif runDropped == 1 & runAdded == 1 ~ case 3
if classBin(min(maxRL, RLDropped» =
classBin(minlmaxRL, RLAdded»
nRk = totalRLs;
i = classBin(minlmaxRL, RLDropped»;
Oi = curHistli) + 1;
ppi = refHistli);
j = classBinlmin(maxRL, RLAdded);
OJ = curHist(j)  1;
pj = reEHist (j) ;
chi2 = chi2 + (2*Oi + l)!(nRk*ppi) + 12*Oj + l)!(nRk*pj);
end
runDropped = 0;
rUnAdded = 0;
end
y(k) = ehi2;
end
t = toe;
81
giobal executionTjmeChi2Newi
executionTimeChi2New = executionTimeChi2New + ti
82
A.3 Code for the New Algorithm Using KS Test
9., CPM KS lITe"!. m
% New algorithm Eor Controller Performance Mo . or u ing KS test
% by Qing Li, Last m0dified: 11/23/2003
• inputs: x is a vector containing all data (actuating errors)
% N is the window size, # of data within data window
~ maxRL is the maximum RL considered
% output: y is a vector containing a sequence of chisquare t0st
% 3tatistic at each sample time
function y = CPM KS_New(x, N, maxRL)
tic;
% reference distribution
refHist = zeros(l, maxRL);
refHist = 0.5. A (1:maxRL);
refHist(end) = 1  sum(refHist(l:end  1));
cumRef = cumsum(refHist);
RLS = zeros(N, 1); ~RL sequence in window
hRL = zeros(l, maxRL);
runLen = 1;
InewestRL = 1; %index of newest run i.n ~indow
RLS(InewestRL) = 1; %initial value fo~ new run
% build RL distribution for 1st data window
sPrev = sign(x(I)); %sign of previous data
if sPrev == 0
sPrev = 1;
end
for k = 2:N
s = s i9n (x (k) ) ;
if s = 0 & s = sPrev
if runLen <= maxRL
hRL(runLen) = hRL(runLen) + 1;
else
hRL(maxRL) = hRL(maxRL) + 1;
end
RLS(InewestRL) = runLen;
InewestRL = InewestRL + 1;
sPrev = s;
runLen 1;
else
runLen runLen + 1;
end
end
cumHist cumsum(hRL) ;
[maxCumD, RLmaxCumD] = max (abs (cumHist  cumRef));
ibuild RL distributions for subsequent WIndows
IoldestRL = 1; ~indez of oldest run i~ window
RLToDrop = RLS(IoldestRL);
RLpartOld = RLToDrop;
runDropped =' 0;
runAdded = 0;
83
n = length (x) ;
y = zeros (n, 1);
for k = (N + l):n
RLpartOld = RLpartOld  1; ~shift out 1 data
if RLpartOld == 0 %oldest RL dropped out
hRL(min(maxRL, RLS(IoldestRL))) =
hRL (min (maxRL, RLS (IoldestRL) ))  1;
runDropped = 1;
IoldestRL = IoldestRL + 1;
if IoldestRL > N
IoldestRL = IoldestRL  N; ~next oldest
end
RLDropped = RLToDrop;
RLToDrop = RLS(IoldestRL);
RLpartOld = RLToDrop;
end
%update newest RL after adding one new data
s = sign(x(k));
if s = 0 & s : sPrev
if runLen <= maxRL
hRL(runLen) = hRL(runLen) + 1;
else
hRL(maxRL) = hRL(maxRL) + 1;
end
RLS(InewestRL) = runLen;
RLAdded = runLen;
InewestRL = InewestRL + 1;
~f InewestRL > N
InewestRL = InewestRL  N;
end
sPrev = s;
runLen = 1;
runAdded = 1;
else
runLen runLen + 1;
end
% goodnessoffittest;
%case 1
if runDropped == 1 & runAdded 0
maxD1 = 0;
for i = RLDropped:maxRL
cumHist (i) : cumHist (i)  1;
if abs(cumHist(i)  cumRef(i)) > maxDl
t maxDl is the max difference in
% affected bins of cLinHisl:
maxDl = abs (cumHist (i)  cumRef (i));
RLmaxDl = i;
end
end
~ f maxDl > maxCumD %part mazDl is ne..,! ma.%CumD
maxCumD = maxD1 ;
RLmaxCumD = RLmaxDl;
elseif.maxDl < maxCumD & RLmaxCumD >= RLDropped
~pr~1J maxCumD bi.n [lequert':/ lS affected b"j
~droppi.ng a run, S0 unfcf:ct",d bins also
~n!';ed to bE: sea rch",d f0," 1Ti<l;~CI.i "D
84
maxCumD = maxDl;
RLmaxCumD = RLmaxDl;
for j = l:RLDropped  1
i~ abs(cumHist(j)  cumRef (j}) > maxCumD
maxCumD = abs(cumHist{j}  cumRef(j);
RLmaxCumD = j;
end
end
end
runDropped 0;
$;case 2
1
o·,
elseif runDropped == 0 & runAdded
maxDl = 0;
for i = RLAdded:maxRL
cumHist(i} = cumHist(i} + 1;
:; f abs (cumHist (i)  cumRef (i)) > maxD1
maxDl = abs(cumHist(i}  cumRef(i»);
RLmaxD1 = i;
end
end
if maxD1 > maxCumD ~"PC:l:'':: meUO: ~s \lev) maxCumD
maxCumD = maxDl;
RLmaxCumD = RLmaxDl;
else~f maxD1 < maxCum~ & RLmaxCumD >= RLAdded
maxCumD = maxD1;
RLmaxCumD = RLmaxD1;
for j = l:RLAdded  1
if abs(cumHist(j}  cumRef(j» > maxCumD
maxCumD = abs(cumHist(j}  cumRef(j);
RLmaxCumD = j;
end
e~d
end
runAdded
(; case 3
elseif runDropped 1 & runAdded 1 ...
& RLDropped RLAdded
maxD1 = 0;
if RLDropped < RLAdded
for i = min (maxRL, RLDropped) :min{maxR~, RLAdded)  1
cumHist(i) = cumHist{i}  1;
if abs {cumHist (i)  cumRef (i)} > maxD1
maxD1 =' abs(curnHist{i}  cumRef(i)};
RLmaxD1 = i;
end
end
elseif RLDropped > RLAdded
for i = min(maxRL, RLAdded} :min(maxRL, RLDropped)  1
cumHist(i} = cumHist(i) + 1;
5 f abs (cumHi st (i)  cumRef (i)} > maxD1
maxD1 =' abs{cumHist(i}  cumRef(i);
RLmaxD1 =' i;
end
end
end
8S
if maxDl > maxCumD
maxCumD = maxD 1 ;
RLmaxCumD = RLmaxDl;
elseif maxDl < maxCumD .
& min( [RLDropped RLAdded]) <= RLmaxCumD ...
& RLmaxCumD < max( [RLDropped RLAddedJ)
;;,prev maxCumD Din is affected by run dropping/adding
%50 J~~ffected bins also needs co be searched
maxCumD = maxDl;
for i = l:min( [RLDropped RLAdded])  1
if abs (cumHist (i)  cumRef (i)) > maxDl
maxDl = abs(cumHist(i)  cumRef(i));
RLmaxDl = i;
end
end
for i
if
max( [RLDropped RLAdded]) :maxRL
abs(cumHist(i)  cumRef(i)) > maxDl
maxDl = abs(cumHist(i)  cumRef(i));
RLmaxDl = i;
end
end
end
runDropped = 0;
runAdded = 0;
end
%maxCUlnD, RLmaxCumD
y(k) = maxCumDi
er.d
t = toci
global executionTimeKSNew;
executionTimeKSNew = executionTimeKSNew + t;
86
A.4 Code for Calculating the Theoretical RL Distribution Probability Functions for
1stOrder MA Processes
% pfRLlstMA.m calculate theoretica RL probability functions
% for ~3'order MA Processes
~; by Qing Li, Last modified: 09/17/2003
% input: b is vecto:':" of the iY1A j:'rocess lfIode 1 parameters
% output: y is vector of probabili.ty functLons for RL =;., 2, .,., ll\a}.RL
function y = pfRLlstMA(b)
bO = b (1) ;
bl = b(2);
xlimit = 5;
dx = le4;
rmax = 10;
pf = zeros (rmax, 1);
pfl pf;
pf2 = pf;
x = (xlimit:dx:xlimit)';
xlen = length(x);
'i;ib is the index of bl/bOx,
ib = bl/bO*(xlimit + dx*(O:xlenl) ,);
ib = round ( (ib + xlimit)/dx) + 1;
%[ib(1:4} I ;ib(xlen3 :xlen)' 1
ib max (1, ib); ~;set to 1 all ib<J
ib = min (xlen, ib); %set to xlen all ib::xlen
%caJ~ulate LOa, Otl'1 integ<<11 arri.J'j
ipdfa '= pdf (:{) ;
pdfa = exp(x. A 2/2)/sqrt(2*pi);
';siOa '= i.O (x); %x is integra ljrnit ( .lnr:, x)
iOa = (1 + ed (x/sqrt (2))) /2;
~);iO,/, iv is iOa, ia in the order of ~b
iOv = iOa (ib) ;
%i'.} 's ia values in (bl/bO*/., infl order
'hvhi'e x i.2 in (inf, in£) orde,
*calculate ia, integral arral
ia = zeros (xlen, rmax);
% pdfi is pdf ~or integral
pdfi = pdfa.*iOv*dx;
cum cumsum(pdfi) ;
for k l:rmax ~k is PL, ika
'1:tmp = pdfa."iv(:, kj.*dz;
'i:tlllP = CUl1lsUll\{pdfi (:, X»);
ia ( :, k) = cum (end)  cum;
iv '= ia(ib, k);
pdfi = pdfa.*iv*dx;
cum = cumsum(pdfi);
87
pf(k) = 2*sum(pdfa.*cum(ib)*dx);
pf1 (k) 2*cum(end) ;
pf2 (k) = 2*ia(floor(xlen/2) + 1, k);
end
(pC pfl, pf2J
y = pf/pf1(1)
sm = sum(y)
plot (y) ;
88
A.5 Code for Calculating the Theoretical RL Distribution Probability Functions for
2lldOrder MA Processes
~ pfR':..,2ndMA.m calculate theoreticaJ RL probabi ~ity functions
~ for 2ndordeL MA processes
% by Qing Li, Last modified: 09/22/2003
% input: b is vector of tile r~A process '11odel parameters
% output: y is vector of probability functions (or RL ~ I, 2, _., maxRL
function y = pfRL2ndMA(b)
~6b [lll]
bO b (I) ;
b1 b (2) ;
b2 b (3) ;
xlimit = 4;
dx = IeI;
rmax = 10;
pf = zeros (rmax, I};
pf1 pf;
pf2 = pf;
x = {xlimit:dx:xlimit} ';
xlen length(x) ;
pdfa exp(x. A 2j2}/sqrt(2*pi} ;
% 20(x) single variable
i20x = (1 + erf(x/sqrt{2)) )/2;
%i20x1Y.2
for i = l:xlen
for j l:xlen
xi b1/bO.*x(i}  b2/bO.*x(j);
Xl min(max{xi. xlimit}, xlimit);
xi round((xi  (xlimit»)/dx) + 1;
i20x1x2(i, j} = i20x(xi);
end
end
i2xlx2i = zeros(xlen, xlen, rmax);
for i = l:xlen
for j l:xlen
xi bl/bO.*x(i}  b2/bO.*x(j};
xi min(max(xi, xlimit}, xlimit};
xi round((xi  (xlimit})/dx) + 1;
tmp = i20xlx2(:, i).*pdfa.*dx;
tmp = cumsum (tmp) ; .
i2xlx2i(i,'j, 1) = tmp(end}  tmp{xl};
end
end
for k = 2:rmax
si21xlx::' zeros (xlen. xlE:~1);
Ear i = 1: xlen
89
for j = l:xlen
tmp = i2xlx2i(:, i, k  l).*pdfa.*dx;
tmp = cumsum(tmp);
xi bl/bO.*x(i)  b2/bO.*x(j);
xi = min(max(xi, xlim.it), xlimit);
xi = round((xi  (xlimit)/dx) + 1;
i2xlx2i(i, j, k) = tmp(end)  tmp(xi);
end
end
end
't;denumerator of PRl
iin = zeros(size(x»); 'integral inside 2 'nf integrals
for i = l:xlen
iin(i) = sum(i2xlx2i(:, i, 1).*pdfa.*dx);
end
den = sum(iin.*pdfa.*dx);
~ numerator integral inside two inf integrals
iNumeratorln = zeros (xlen, xlen, rmax);
for k = l:rmax
for i = 1:xlen
for j = 1:xlen
tmp = i2xlx2i ( :, j, k). *pdfa _*dx;
tmp = cumsum(tmp);
xi b1/bO.*x(j) b2/bO.*x(i);
xi = min(max(xi, xlimit), x:.Limit);
xi = round ( (xi  (xlimit)/dx) + 1;
iNumeratorln(i, j, k) = tmp(xi); ttmp(end)  tmp(xii;
end
end
er.d
~nurnerator of PRi
numR = zeros (rmax, 1);
iinf1 = zeros(xlen, rmax); %integral after 1st inf integral
toc k = l:rmax
tor i = l:xlen
iinfl(i, k) = sum(iNumeratorln(i, ., k)'.*pdfa.*dx);
end
numR (k) = sum (i inf1 ( :, k). *pdfa. *dx) ;
end
~~Rl, probability functions
PRi numR/den
y = PRi;
5m2 = sum(y)
plot(y)
90
VITA (i).
QingLi
Candidate for the Degree of
Master of Science
Thesis: DEVELOPING COMPUTER SOFTWARE FOR CONTROLLER
PERFOMANCE MONITORING
Major Field: Computer Science
Biographical:
Education: Graduated from Zhuzhou NO.2 High School, ZhuzhoLl, Hunan
Province, China in July, 1981. Received a Bachelor of Science degree and
a Master of Science degree both in Chemical Engineeritlg from Zhejiang
University, HangzhoLl, Zhcjiang Province, China, in July, 1985 and July,
1988 respectively. Received a Doctor of Philosophy degree in Cllemical
Engineering from Oklahoma State University, Stillwater, Oklahoma in
May, 2002. Completed the requirements for the Master of Science degree
at Oklahoma State University in December, 2003.
Professional Experience: Employed as a process engineer by Sinopec Shanghai
Petrochemical Company, Shanghai, China from June, 1988 to July, 1996.
Employed as a graduate research assistant in Advanced Process Control,
School of Chemical Engjneering, Oklahoma State University, from
August, 1996 to May, 2002. Employed as a research assistant in School of
Electrical and Computer Engineering, Oklahoma State University, from
August, 2002 to December, 2002, and from May, 2003 to August, 2003.
Employed as a teaching assistant in School of Electrical and Computer
Engineering from January, 2003 to present.