Ubc 1999-463419

CONTROL LOOP P E R F O R M A N C E ASSESSMENT A N D OSCILLATION
DETECTION
By
Lahoucine Ettaleb
B. Sc. Ecole Hassania des T. P. Casablanca, 1987
M . Sc. A. Ecole Polytechnique Montreal, 1994
A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF
THE REQUIREMENTS FOR T H E DEGREE OF
DOCTOR OF PHILOSOPHY
in
THE FACULTY OF GRADUATE STUDIES
DEPARTMENT OF ELECTRICAL A N DCOMPUTER ENGINEERING
We accept this thesis as conforming

to the required standard
T H E UNIVERSITY OF BRITISH C O L U M B I A
July 1999
© Lahoucine Ettaleb, 1999
In presenting this thesis in partial fulfilment of the requirements for an advanced
degree at the University of British Columbia, I agree that the Library shall make it
freely available for reference and study. I further agree that permission for extensive
copying of this thesis for scholarly purposes may be granted by the head of my
department or by his or her representatives. It is understood that copying or
publication of this thesis for financial gain shall not be allowed without my written
permission.
Department of ^lg-^VW <~»X gu^«k Cs^ôvW— wxWi^
The University of British Columbia

Vancouver, Canada
DE-6 (2788)
Abstract
In this thesis an analysis of process control loop performance based on the Harris's performance
index is given. Good and poor performances will be defined and it will be shown that poor
performance may be obtained when some loops are cycling. Then a procedure for localizing
the oscillating loops in a cascade multiloop system will be introduced. The procedure uses
data collected under normal operation and assumes knowledge of the different frequencies of
the periodic component of the output signal. An oscillation index is introduced to characterize
the oscillating loops.
The estimation of time delay is a major step in assessing control loop performance. To
overcome some difficulties of the available methods for its estimation, a new off line extended
Least-Squares technique for parameter identification and time delay estimation will be intro-
duced for processes where acting disturbances are white noise (i.e. C(q^ ) = 1).
1
The Harris performance index for single-input single-output (SISO) systems is a very use-
ful tool, therefore its extension to multivariable systems becomes necessary. The multi-input
multi-output (MIMO) performance indices available in the literature involve use of the MIMO
time delay matrix, also called the interactor matrix, to perform the assessment. The method
proposed in this thesis does not require knowledge of the interactor matrix. The method uses
knowledge of the delays between different input/output pairs of the process and defines an
absolute lower bound on the achievable output variance for each output. This bound is then
used to define a performance index associated with each output. The sum of these bounds is
used to characterize the overall control loop performance. The results are independent of the
matrix time delay representation, or interactor. This method is evaluated by application to the
loops controlling a lime kiln in a Kraft pulp mill.
ii
Table of Contents
Abstract ii
List of Tables vi
List of Figures vii
Acknowledgment ix
1 Introduction 1
1.1 Control Loop Performance Monitoring 1
1.2 Outline of this thesis 8
2 Control Loop Performance Monitoring: SISO case 11
2.1 Minimum Variance Control (MVC) 11
2.2 Evaluation of performance 15
2.2.1 Evaluation of I using time series

v 15
2.2.2 Evaluation of I using Laguerre network

p 18
2.2.3 Time delay estimation 19
3 A n Extended off-line Least-Squares Method and Time Delay Estimation 22
3.1 Introduction 22
3.2 Least squares parameter estimation 24
3.2.1 Least-squares method 24
3.2.2 Generalized least-squares 26
3.2.3 Extended least-squares 29
3.3 An ARMA model for a closed-loop system 30
iii
3.4 A n off-line extended least-squares method 34
3.5 Off-line time delay estimation 39
3.6 Case of colored noise 45
3.7 Conclusion 46
4 Realistic Performance Assessment in the SISO Case 47
4.1 Introduction 47
4.2 Problem formulation 47
4.3 Satisfactory and unsatisfactory performance 50
4.4 Control chart 54
4.5 Tuning procedure 55
4.6 Simulation Example . 57
4.7 Conclusions 60
5 Monitoring Oscillations in a Multiloop System 61

i
5.1 Introduction 61
5.2 Background 63
5.3 Monitoring the oscillations 65
5.4 Accuracy of the Describing Function Approach 72
5.5 Some Limitations 73
5.6 Simulation example 73
5.7 Conclusion 83
6 Control Loop Performance Monitoring in the M I M O case 84
6.1 Assessment strategies based on MIMO minimum variance control (MVC) . . . . 85

6.2 Ultimate output variance bound 88
6.3 MIMO Performance Indices 91
6.4 Estimation of the lower bounds 95
6.5 Simulation examples 98
iv
6.5.1 First example 98
6.5.2 Second example 101
6.6 Conclusions 104
7 Performance Assessment of Control Loops on a Lime K i l n 106
7.1 Process description 106
7.2 Lime kiln process 108
7.3 Loop performance assessment for the lime kiln 112
7.4 Conclusion 115
8 Summary 116
8.1 Contributions of this thesis 116
8.2 Recommendation for future work 117
Bibliography 118
v
List of Tables
3.1 Results 43
6.2 C i = 0.9 <TI = 1 CT = 1

2 95
6.3 Reference values for example 2 101
vi
L i s t of Figures
2.1 Process control 12
3.3 Covariance of the residual 33
3.4 White noise and its estimate 34
3.5 Estimation value of the coefficient a 41
3.6 Estimation value of the coefficient b 42
3.7 Estimation value of the time delay d 42
3.8 Loss function V 43
3.9 Estimation value of the coefficient a by ordinary least-squares 44
3.10 Estimation value of the coefficient b by ordinary least-squares 44
3.11 Estimation value of the delay d by ordinary least-squares 45
4.13 I and steady state error versus P-gain

p 49
4.14 A single nonlinear loop 51
4.15 Saturation characteristic 51
4.16 Nyquist Diagram 52
4.17 Partition of Ip curve 52
4.18 Performance analysis 54
4.19 Process approximation 55
4.22 Comparison of true and approximate model 58
vii
4.23 Estimated curves versus Controller parameters 59
4.24 Estimated curves versus Controller parameters 60
5.4 The output y(t) of the outer loop, example(5.6.1)-test 1 75
5.5 The output yi(t) of the inner loop, example(5.6.1)-test 1 76
5.6 The output out\(t), example(5.6.1)-test 1 76
5.7 The output outfit), example(5.6.1)-test 1 77
5.8 The output y{t) of the outer loop, example(5.6.1)-test 2 78
5.9 The output y\(t) of the inner loop, example(5.6.l)-test 2 78
5.10 The output out\(t), example(5.6.1)-test 2 . 79
5.11 The output out2(t), example(5.6.1)-test 2 79
5.12 The output y(t) of the outer loop, example(5.6.2) 80
5.13 The output yi(t) of the inner loop, example(5.6.2) 81
5.14 The output outi(t), example(5.6.2) 81
5.15 The output out (t), example(5.6.2)

2 82
5.16 The output outfit) versus iri2(t), example(5.6.2) 82
6.17 On line performance index estimation kl=0.08 and k2=0.05 104
7.18 Lime cycle 107
7.19 Lime Kiln 108
7.20 First bump test 110
7.21 Second bump test Ill
7.22 Validation of the model 112
7.23 On-line performance indices 114
viii
Acknowledgment
I am grateful to a number of individuals who greatly assisted me throughout the course of this
thesis. Without their support this work would not have been be achieved. In particular I wish
to thank:
Both Dr. Michael S. Davies and Dr. Guy A. Dumont my thesis supervisors for giving me
the opportunity to work under their guidance. The research reported in this thesis reflects true
interactions, inspiration and encouragement from both of them. I am very grateful for their
direction and their financial support.
Dr. Ezra Kwok for his valuable suggestions and help with my research work.
Rita Penco, the magic Librarian for providing any required documentation right on time.
Georgina White, Ken Wong, Brenda Dutka , Lisa Brandly, Peter Taylor and Tim Paterson
for their outstanding assistance. Thanks are also due to Brian D. McMillan for computer
support and to my student colleagues and friends at the Pulp and Paper Centre.
Ministere des T.P. of Morocco, Natural Sciences and Engineering Research of Canada and
Networks of Centres of Excellence of Canada for providing the financial support.
Dr. Abdelhaq Elhakimi, Dr. Mohamed Y. Tachafine, my colleagues at Ecole Hassania des
T.P. for their moral support.
John Ball from Canfor pulp and paper mills and Bruce Allison from Paprican for providing
lime kiln data used in chapter 6.
My wife Fatima Elmourabit, my children: Sarah, Fatima Azzahra and Hamza for their
unconditional love and unselfish support.
My large family, brothers, sisters and friends back in Morocco.
I pay special tribute to Sandra Davies for her kindness and care. She was always there when
myself or my family needed some one to help. I am indebted to remain grateful to her.
ix
Finally, it is to my mother Hajja Zahra Irramed that I dedicated this thesis for her eternal
love and encouragement.
x
Chapter 1
Introduction
1.1 Control Loop Performance Monitoring
The purpose of on-line surveillance of control loops is to obtain information on the performance
of plant components while the plant is operating, without disturbing its normal operation.
This information is used for a number of purposes such as 1) to provide information on the
performance deterioration with time and its correlation with operating conditions so as to help
the diagnosis of a problem; 2) to allow improved control of normal operating conditions so that
the plant can run closer to its limits. The analysis of this information guides decisions on what
repair work or other modifications will be needed at the next shutdown and whether continued
operation is safe under the normal conditions.
On-line surveillance is a general concept, and research in the area is widely treated in the
literature because the problems addressed are as old as the process industry itself Mah [53].
When process measurements and data acquisition were regularly carried out by plant operators
and engineers, the accuracy and consistency of the data were checked on the basis of their
knowledge and experience of the process. Now these functions are performed automatically by
computers; hence, a far greater volume of process data may be checked and processed without
requiring of human vigilance and intervention.
The surveillance in a wide sense includes many activities such as performance monitoring,
fault detection and state or parameter estimation. A common characteristic among these tech-
niques is the use of collected data. But the objective of each method is different. Performance
monitoring attempts to state whether satisfactory system performance standards are met. If
not, then fault detection procedure is to be initiated to locate the faulty component. If there
1
Chapter 1. Introduction 2
is none, then an engineering modification must be undertaken.
Degradations in control performance can be caused by an inaccurate sensor reading due
to bias, or an actuator which is not responding properly to control signals due to a non-linear
phenomenon such as backlash or stiction, or any leaks in pipelines. These malfunctions are
refered to as failures or faults. In some cases these faults can be modeled as an additive
disturbance, and therefore they can be detected using Kalman filter-based state estimators.
Observing changes in the mean and variance of the main process variable can also indicate that
a failure has occured because faults introduce changes in the properties of system components.
Taking this fact into account, some recent methods of fault detection are based on that principle
and many detection algorithms using recursive parameter estimation are now available in the
literature as detailed in [8]. A review of the most significant methods for fault detection is given
in [46, 7].
Fault detection and parameter estimation are two wide areas of research well treated in the
literature. The main focus in this thesis will be on control loop performance monitoring based
on the use of the minimum variance controller as benchmark.
Control loop performance has recently received great attention from researchers and man-
ufacturers who realize that it is a key to manufacturing high quality product. In the past,
statistical such as Shewart control charts were used. Recently, focus has been placed on the de-
velopment of multivariate statistical control charts[55]. Principal Component Analysis (PCA)
and Partial Least Squares are among the best known methods for processes with a large number
of measurements. The idea behind these methods is to compress the highly correlated informa-
tion contained in the data into a lower dimensional space which allows use of univariate control
charts. These charts can also be generated using methods from chemometrics[69].
As mentioned before, the ultimate goal of performance monitoring is to improve product
quality by reducing process variation, Smith [60].
The most commonly used measure of performance is the variance or standard deviation
of key process variables. Standard deviation is used for monitoring because reduction of this
quantity implies control and manufacture of a uniform product in the face of randomly varying
raw materials and chemicals[9]. It was shown in Astrom et al. [3] that by using the concept of
minimum variance control (MVC), it is possible to control the process variability, so that the
remaining variations are due only to a random noise which, when dead time is present, may
be correlated. In some cases it is not possible to implement MVC, and time series analysis
tests exist which identify whether minimum variance has been approached. Pryor [57] discusses
some aspects of the use of autocorrelation and power spectra to analyze paper machine dynamic
data. The analysis of control performance monitoring from dynamic data was also addressed
by Staton [63]. The approach of Harris [35] for control loop performance assessment is based on
the use of routine operating data and little prior knowledge of the process. Although the control
objective may not be to minimize process output variance cr , Harris proposes the use of the
2
minimum achievable output variance as a benchmark to assess the performance of a regulator.
To determine the minimum variance <T^„ , the process time delay of Td sampling intervals must
be known and a model of the closed loop output data must be estimated. An A R M A model of
the form y(t) = H(q~ )e(t) is used. The first Td terms of the impulse response of the filter H
1
(hi,..,hT -i)
d are feedback-invariant and are used to estimate c r ^ as 2
a 2
mv = (l + h + ... + h _ )al
2
1
2
Td 1 (1.1)
This last equation requires knowledge of the process time delay Td- The estimation of this entity
is a well known problem in the literature and, many methods have been developed, among them
the techniques described in Zheng et al.[71] or in Elnaggar and Dumont [22]. The method of the
last reference is used in many applications and proved to give good results[52, 30, 56]. But, it
has also some limitations[52]. In the third chapter, for processes where the acting disturbances
are white noise i.e. the polynomial C(q~ ) = 1 a new method will be proposed.
l
Other performance measures have also been used in the literature. Notable is that used by
Astrom [2] and Astrom et al.[6] for assessing achievable performance using PID control.
Harris's technique has been used successfully to improve quality and uniformity of products
in many applications, for example [56, 52]. Unfortunately, this method has some limitations.
One of them is that the Harris's performance index is only valid for linear minimum phase
systems. Tyler and Morari[66] proposed an extension to a nonminimum phase systems by
taken into account unstable system poles and zeros.
Another limitation is that when comparing the output variance to the minimum output
variance cr^v o m
Y the extreme situations are well defined: good variance performance when
the performance index (ratio of the actual output variance to the minimum achievable output
variance noted I ) is equal to one, and poor performance when the index is large. In practice,
p
the index takes on values between these extremes and it is necessary to determine whether the
actual controller performance is acceptable, acceptable with a small tuning correction, or so
unacceptable that redesign of the controller is needed.
Another issue related to the minimum variance benchmark limitations is whether the per-
formance index itself is appropriate. In practice, the ability of the controller to reduce the effect
of noise on the output variance is only one of several desirable attributes. Therefore, having
I ~ 1 for a certain controller may be too narrow as assessment of performance. In similar fash-
p
ion, I = 1 is obtained when implementing a minimum variance control strategy (MVC) for a
p
minimum phase system. As discussed by Isermann [45], minimum variance strategy may result
in slow output tracking of setpoint changes if the process zeros are too slow. This strategy may
therefore, not always be useful. To overcome some limitations of the minimum variance bench-
mark, Kozub and Garcia[49] propose a performance index which allows an exponential decay
of the closed-loop response. Horch and Isaksson[39] introduce a modified performance index
which does not require more knowledge than the original one. Huang and Shah[41] propose a
"user-defined benchmark" which requires knowledge of the disturbance model. They assume
that the disturbance model is either an integrator or has to be estimated from data.
However, a good control strategy is often obtained when compromising between output
variance as a performance criterion and other measures such as steady-state rejection of distur-
bances. These points will be discussed in more detail in chapter 4 where it is also shown that
poor performance may result when loops are cycling. Oscillations have an impact on product
uniformity, and can cause increase energy consumption and waste of raw material. In a recent
paper Clarke [16] showed that by eliminating oscillations in control loops, 6% saving in the
steam consumption of a batch temperature controlled process was possible. Bialkowski [9], in a
survey, indicated that 30% of all control loops in Canadian paper mills were oscillating because
of sensing and valve nonlinearity problems. Most of the valve nonlinearities encountered in the
pulp and paper industry are caused by backlash or stiction. Backlash is due to the difference in
motion response between an increasing and a decreasing output usually caused by mechanical
gearing. Backlash does not typically result in cycling, except in integrating processes, but it
does degrade control performance by limiting the loop's ability to attenuate load disturbances
and respond to set point changes. Stiction, on the other hand, affects the valve control element
and is more severe than backlash, almost always resulting in loop cycling[68]. Stiction causes
the actuator force to continue to increase without a change in the valve position until the force
overcomes the friction force. Then, the valve will move suddenly to a new position. Stiction
causes reduced valve resolution and results in a limit cycle due to the controller's inability to
position the valve to maintain a setpoint[68]. It is important to mention that although these
valve problems cannot be eliminated by controller tuning, they can however be minimized. But
only repairing or replacing the valve can eliminate them. Oscillations may also be due to bad
control tuning.
Dealing with detection of oscillations due to non-linear problems is a subject of many articles
in the literature. In Horch and Isaksson[38], extended Kalman filter was used to generate
two sets of residuals: One describing a valve with an acceptable level of friction and another
describing a valve with unacceptable amount of stiction. A generalized likelihood ratio is then
used to decide which of the sets of residuals is more likely to correspond to the measured data.
In Taha et al.[65] a friction index is proposed which measures the deviation of the valve's input-
output characteristic in normal operating conditions from its static characteristic established by
the manufacturer. The index measures the degradation of the valve due to its use. In chapter
5, the oscillating loops are first located and then the cause of the oscillations is found. The
procedure uses prior knowledge of different frequencies present in the periodic component of the
output signal. A n oscillation index will be introduced to characterize the oscillating loops. The
prior knowledge of oscillation in the loops can be obtained by techniques such as the Discrete
Fast Fourier Transform (DFFT) or Singular Value Decomposition (SVD)[47].
The Harris method has been extended to Multi-Input Single-Output (MISO) in Desborough
et al.[18] and very recently to Multi-Input Multi-Output system (MIMO) by Huang et al.[42]
and Harris et al.[36]. The last two references use minimum variance control as a benchmark. In
fact with no explicit time delay structure, the key steps in the derivation of SISO and MIMO
minimum variance control strategies are:
• A predictor for the time shifted variable £(q)y(t), where £(q) stands for the matrix time
delay representation, is developed.
• The noise term ^(q)N(q~ ) is factored into future and past noise terms as
1
Z^Niq- )1
= F(q) + Riq' )
1
(1.2)
where F(q) = F qd
d
+ ... + F q and ^ ( g ) is a proper transfer function.
t
-1
• Under M V control, and in steady state, the tracking error is shown to be a finite moving
average polynomial of the form
y(t)-y*(t) = cHq)F(qHt), (1.3)
where {y*{t)} is a given reference trajectory and £(q) is the time delay representation
considered. Note that the expression for this tracking error (equation 6.183) does not
depend on controller parameters.
For SISO systems £~ {q)

l
= q~ , therefore equation (1.3) reduces to equation (1.1) which
Td
expresses the minimum achievable output variance for SISO systems.

Following from SISO assessment, one might be tempted to use the predicted performance
of MIMO minimum variance controller as an absolute lower bound against which current per-
formance can be assessed. A bound is indeed obtained when taking the variance of both sides
of equation (1.3). However, this lower bound will depend on the type of matrix time delay
representation used and so is not unique.
The MIMO time delay representation was first proposed by Elliott et al. [21], in 1982, who
showed that the notion of delay in the MIMO case is related to the system interactor matrix
£(g) introduced by Wolovich et al.[70], in 1976, which satisfies the following two properties for
a given system with the transfer function T:
1. l i m _ £ ( q ) T ( q ' ~ ) = K for any non singular matrix

g 00
1
2- £(<?) = H{q)dia,g{q , ...,q ) where H(q) is a unimodular matrix.

dl dn
The interactor matrix was introduced initially by Wolovich with no apparent justification other
than the resulting computational procedure for its calculation. It was Elliot and Wolovich [21]
who showed the usefulness of this entity and related it to time delay in the MIMO case.
The interactor proposed by Bittanti et al.[ll], in 1994, satisfies
S(q)Z(q- ) = i,
1 T
(1.4)
and also satisfies property (1) above. This interactor is called the spectral unitary interactor
matrix.
As can be deduced for the same system one can calculate at least two inter actors.
The procedure for control loop performance assessment proposed by Huang[40] uses the
unitary interactor. In Harris et al. [36], the procedure is based on the interactor
- 7(9)^(9)
where £(q) is the lower triangular form as defined by Wolovich, and 7(g) is obtained by solving
the following spectral factorization equation
b W f Q i W r 1
) ] = m-TQiM*- )- ]
1 1
(1.5)
where Q\ is a diagonal matrix. In fact, it is easy to show that (1.5) is a more general form of
equation (1.4). Writing (1.5) as
(1.6)
£(9)Qi£VT = Qi
we let Q i = I then f ( g ) f ( g )
_1 T
= I. If Qi # J and |(g) satisfies (1.4) i.e. | ( g ) | ( g )
_1 T
= J,
then by writing Q i = QI QI the interactor matrix
£'(q) = Ql£{q)Qi h
(1.7)
is also a solution to (1.6). On the other hand if £'(g) satisfies (1.6) then the interactor £(g)
defined in (1.7) satisfies (1.4). Therefore (1.6) serves as a link between (1.4) and (1.6), showing
the equivalence of the inter actors.
Knowledge of the interactor matrix is a drawback in the approaches of Huang and Har-
ris. Though Huang et al.[43] showed that only a few Markov parameters are needed for its
determination. If the assessment is meant to be on-line the estimation of the first few Markov
parameters is not straightforward for the following reasons:
• Most processes that require monitoring are under feedback control.
• For most processes their linear model depends on the operating point which keeps chang-
ing therefore, the interactor matrix also does.
For these and other reasons many authors avoid use of the interactor matrix for adaptive
control (see for example Dugard et al.[19], Shah et al.[58], ...etc). However, the approach which
is proposed in chapter 6 is based on defining an absolute lower bound on the achievable value for
each output variance, thus avoiding the need to identify the interactor matrix. It also removes
the ambiguity in the lowest achievable performance.
1.2 Outline of this thesis
As mentioned before, this thesis is concerned with control loop performance monitoring based on
minimum variance control as a benchmark. The research accomplished will be presented in six
chapters. In the second chapter, the Harris approach for control loop performance monitoring
is described along with different techniques for its implementation for single-input single-output
systems. The Harris method requires knowledge of process time delay, therefore, in the third
chapter, an off-line extended least-squares method to identify the parameters of a linear system
under closed-loop conditions for which the acting disturbances are white noise (i.e. C{q~ ) = 1)
l
is proposed. An off-line time delay estimation method, as an application of the proposed
method, will also be presented . These results have been reported in [24].
For some applications the Harris performance index does not lead to a clear cut assessment.
In the fourth chapter, an analysis of process control loop performance based on Harris's per-
formance index is given. Good and poor performance will be defined and it will be shown that
poor performance may be obtained when some loops are cycling. A controller tuning procedure
is given that is based only partly on output variance. These results have been reported in [27].
The important result from this third chapter is that poor performance results mainly from the
fact that some loops are cycling. The fifth chapter, describes a procedure for localizing the os-
cillating loops in a cascade multiloop system. The procedure uses data collected under normal
operation and assumes knowledge of the different frequencies of the periodic component of the
output signal. An oscillation index is introduced to characterize the oscillating loops. These
results have been reported in [26]. The detection and removal of oscillations caused by non
linear effects is necessary if the linear methods of chapter 3 are to be used. For saturating gains
the system can be brought into the linear region and then identified by controller adjustment.
For more complex uncertainties it may not be possible to revert to linear operation.
In the sixth chapter, a unified method for assessing single and multivariable control loop
performance which does not require knowledge of the interactor matrix is presented. The
method uses knowledge of the delays between different input-output pairs of the process and
defines an absolute lower bound on the achievable output variance for each output. This bound
is then used to define a performance index associated with each output. The sum of these bounds
is used to characterize the overall control loop performance. The results are independent of
the matrix time delay representation, or interactor. The application of the proposed method
to an industrial process is presented in the seventh chapter, the proposed method is evaluated
by application to the loops controlling a lime kiln in a Kraft pulp mill. In the last chapter
conclusions and recommendations for future work will be given.
The chronological use of the results of this thesis in a procedure to to improve control loop
performance for a given system is as follows:
1. Oscillation detection: Investigate oscillating loops as a first step because oscillations
cause poor performance and may lead to an incorrect calculation of the performance
index. In the case of cascade loops, the results from chapter 5 are applied.
2. Performance indices calculation: Two cases are to be distinguished.
• SISO systems for which methods from chapter 2 and chapter 3 are applied
• MIMO systems for which results from chapter 6 are applied.
3. Interpretation of performance indices: For SISO systems, chapter 4 describes a
realistic interpretation of the performance index and a tuning method that are based on
the use of the performance index I po which can be estimated by the relay method.
Chapter 2
Control Loop Performance Monitoring: SISO case
The objective of this chapter is to describe the Harris single-input/single-output (SISO) ap-
proach for control loop performance monitoring and discuss the available methods for its on-line
implementation. As mentioned before this technique relays on the use of minimum variance
control as benchmark therefore the organization of this chapter is as follows. In the next sec-
tion a brief discussion of minimum variance control will be given followed in section 2 by the
definition of the performance index and how to determine it from operating data. The last
section deals with time delay estimation.
2.1 Minimum Variance Control ( M V C )
The concept of minimum variance control was initially developed by Astrom et al.[3] to reduce
process variability due to random disturbances. The derivation of M V C is outlined below.
Let us consider the basic process control shown in figure (2.1). The linear system with delay
is represented by the discrete time transfer function
G
^ = ^&)^ < 2 8
>
where A; is a positive integer, B and A\ are polynomials in q~ , the backward shift operator
l
(i.e q y ( i ) = y{t — 1)). The order UA of A\ is considered to be the same as or higher than n
_1
B
that of B. The disturbance model represented by
G
^ > = -£fk
_1
<->2 9
The operator A is defined as A = 1 - q~~ . The polynomials A^q )

l -1
and C(q~ ) are of order
1
p and q respectively, and are stable.
11
Chapter 2. Control Loop Performance Monitoring: SISO case 12
u(t) + -f-
y(t)
Gp(q') o —
Figure 2.1: Process control.
The expression for the output is
m
- l ^ ^ w +
< 2
- 1 0 )
By long division ^^-IJA c a n

be expressed in the following form (leading to a Diophantine
equation)
ntn-i\ . . n(n~ \ l
(2.11)
where a unique solution is obtained if the polynomial F(q~ ) is monic and of order k — 1 1
(F(q~ ) = 1 + /i<Z +... + fk-iq~ ),

1 -1 k+1
and G(q~ ) is polynomial of order max[q — 1, p +rf— 1].
1
Now, by using equations (2.10) and (2.11), the following expression gives the output at time
t+k
y ( t + k ) =
(^^< ) t +
%V®) (
+Fe t + k
) (- )
2 12
= y(t + k) + Fe(t + k). (2.13)
At time t + k, the best estimate of y(t + k) is y(t + k) because the term Fe{t + k) represents
the information yet to come, as e is zero mean and e and y are uncorrelated.
If the control law
BA AFu(t)
2 + GAiy(t) = 0 (2.14)
is used, then, from (2.12),
y(t + k) = Fe{t + k)
i.e. y(t + k) is a linear combination of the future values of the white disturbance sequence. It is
not difficult to show that E[y (t + k)] is minimized in that case. The output variance is given
2
by
a'mv = E[y\t + k)} = (1 + ft + ... + fU)<- (2.15)
The control law defined by (2.14) represents the minimum variance control law.
In closed loop form (Figure 2.2), the transfer function between the white noise e and the
output y is given by
y(t) (2.16)
i | r t(Bg«)(r')
= H{q~Mt). (2.17)
By using the same manipulation as before one can easily derive the following expression of the
output at time t
rw x / MG - BA FG A 2 C \ , ,
= H(q-i)e(t) - (219)
which shows that the first k Markov coefficients of H(q~ ) are feedback invariant. Then, it
1
follows that
ol > o .
2
mv (2.20)
Under (MVC) the variance of the output is minimized and it is given by
2 _ 2
which expresses the minimum output variance that can be achieved (limit of performance).
The limit is independent of controller parameters. It depends only on the process noise model
and the present system time delay. Therefore, the evaluation of a^ av gives insight into what
performance may be expected from the actual system structure. If the performance is not
satisfactory, then by referring to equation (2.15), change to the system structure can be made
by either reducing the time delay by re-locating the sensors or reducing the disturbance variance.
Equation (2.20) shows that minimum variance control is the best possible control in the
sense that no regulator can have a lower variance. The connection between M V C and other
strategies such as the PID controller and Dahlin's algorithm is discussed in Harris [37].
From equation (2.14) we can deduce that MVC control includes the inverse process model
in the controller. It follows that for a non minimum phase process such a controller will result in
an unstable control system. The problem can, however, be overcome by the technique discussed
in detail by Astrom et al.[4]. It should also be noted that it is often undesirable to use M V C due
to the excessive control action required. This matter will be discussed extensively in chapter 4.
2.2 Evaluation of performance
Even though the MVC cannot always be implemented, it still remains a convenient bound
on achievable performance against which other controllers can be assessed. An index can be
obtained by defining a performance index I based on the theoretical limit as

p
1 < Ip < oo.
Then the following comments can be made:
• I ~ 1 : the actual performance is close to the best achievable performance.

p
• Ip « oo : the actual performance is significantly worse than optimum.
To evaluate I one needs to identify the closed loop transfer function (equation 2.16) and
p
express it as a Diophantine equation to obtain the F(q~ ) required to determine the minimum
1
achievable output variance. I can then be computed [52]. An alternative to closed loop transfer
p
identification is to determine the transfer function between white noise and the output using
Laguerre network[52] or using time series analysis [18]. These two methods are summarized in
the following subsections.
2.2.1 Evaluation of I using time series

p
The procedure for evaluating the performance index I using time series analysis is detailed
p
in Desborough et al. [18]. This procesure is described below. It consists of fitting y(t) to an
autoregressive moving average (ARMA) time series of the form
y(t)-py = H{q- )e(t)-

l
where p is the average of y(t). From the first section it was established (equation (2.18) that
y
y(t)-py = H{q- )e(t)

l
(2.21)
(2.22)
= iPi{q- )e{t) + i}>2{q- )e(t - k)

l l
(2.23)
(2.24)
Equation (2.24) is also equivalent to

oo
y(t) - ^ = e'(t) + oii(y(t - k - i + 1) - p ).
y
i=l
In practice, the infinite series is truncated after m terms

m
y(t) - p y = e'(t) + oti(y(t - k - i + 1) - p)
y
i=l
= m
and, using matrix notation
Y = Xa + e
where
y{n)
y(n - 1)
Y = , a =
y(k + m) C*m
y(n — k) y(n — k — 1) y(n — k — m + 1)

y(n — k — 1) y(n — k — 2) y(n — k — m)
X =
y{m) y(m - 1) y(i)
An estimate for is given by the residual mean square error

By using the mean square error of y(t), an estimate of I is given by

p
s2
-(n — k — m + 1). (2.25)
YY
T
+ f4
If fiy is not known, it has to be estimated. For that, the closed loop model behaviour is modified
to include a constant term an. The parameter vector a is augmented by cio a n (

i the matrix
X is also augmented by a leading column of ones. These transformations give the following
equations:
Y = Xa + e (2.26)
where
y(n) OtQ
y(n - 1)
Y = , a= a.2
y(k + m)
1 y(n — k) y(n — k- 1 ) ... y(n- -k-m + l)
1 y{n — k —1) y(n - k - 2 ) ... y{n — k — m)

X=
1 y(m) y(m - 1) y(i)
The method described above gives an estimate off-line of Ip. The on line estimate can be
obtained by using a recursive least squares algorithm to estimate the vector a in equation
(2.26) and the following recursive equations to estimate the minimum variance performance as
well as the moving mean square error of the output
Smv(t) = \S v(t-l)
m + (l-\)e\t)
S {t) =
y XS {t - 1) + (1 - A)y(t)
y
2
where A is a forgetting factor (0.98 < A < 1) and e(t) is the prediction error at time t.
An alternative method for evaluating I , to the one above is to use a Laguerre Network
p
Lynch and Dumont [52] as is summarized below.
2.2.2 Evaluation of I using Laguerre network

p
When using a state space representation of the relation between the white noise and the output
x(t +1) = Ax(t) + be(t) (2.27)
y(t) = Cx(t) + e(t), (2.28)
the minimum achievable output variance is expressed in terms of the Markov parameters of the
process Cb, CAb, CA b, 2
... as
k-1
-2 _ 2 l + îpA^b) 2
(2.29)
i=l
The same formula is still valid when using a Laguerre network instead of state space representa-
tion [52]. In this case, the system of equation (2.27)-(2.28) denotes the discrete Laguerre state
space representation of Gunarsson [33], where the elements of matrix A and those of b and C
are given by
a,ij = a i = 3
A = = 0 i < j
«u = (-«r (l-« ) i>3 j_1 2
bi = (-of-Vl-a 2
i = ltoN
C = [gi ... g \- N
The parameter a is the time scale, N is the number of Laguerre filters used (The optimal value
of N is derived by Fu et al.[29]) and gw represents the N th

Laguerre gain. The Nth
Laguerre
filter xj^{q) is given by
q-a \ q—aJ
In the above analysis it was assumed that the time delay was known, but in most cases it
is either unknown or time varying. The next section deals with the time delay's definition and
estimation.
2.2.3 Time delay estimation
The estimation of time delay is a familiar problem in control system design, in acoustics,
and in radar remote control. Several good techniques have been developed for extracting
this quantity from measurements. Most of those methods are based on analyzing the cross
correlation function of u and y as the two signals of interest
Ru (k) = E[u(t - k)y(t)}

V (2.30)
where E[.] denotes the expectation over time t to determine the time delay between u and y
[48].
In some cases, when the signals are highly coloured, the cross-correlation function does
not provide the correct time delay. This problem can be overcome by whitening the signals
before performing the cross correlation Carter [14]. In Elnaggar et al. [22] the concept of
variable regression estimation (VRE) to estimate the time delay is used. It is shown that
this estimation can be added to any standard parameter estimation with minimal additional
computation, which means faster convergence of the identification procedure. The model used
to describe the process time delay is a first order given by

_T
a = —e T
e-T sd bq- ~
k 1
< b= 1 - e "
l+TS 1 + aq- 1
then
y(t) = ay(t - 1) + bu(t - k - 1)
the prediction error e is defined as
e(t,k) = y{t)-y{t, k).

where y(t, k) is the prediction of y given the estimate k. The argument k that minimizes
E[e(k) ] provides an estimate of time delay.

2
For b > 0 it is shown in [22] that this is equivalent to maximizing E\ given by
Ex = Ru {k) - aRuy(k - 1).

V
The choice of a = 1 gives unbiased estimation of delay. E\ becomes
E = E{{y{t +
1 l)-y(t))u{t-k)}.
The on-line delay estimation procedure uses the following algorithm:
• Ei(t, ki) = \Ei(t - 1, hi) + [y{t) - y(t - l)]u(t - ki) Vki e [k , min k ] max
• choose the delay that maximizes E\
Ei(t, ki) = max Ei(t, ki) Vfcj € [k

mini kmax\•
In practice this procedure gives good results in open loop [52, 30, 56]. Unfortunately this
technique has also some limitations [52]. In the next chapter a new method for time delay
estimation will be proposed.
Some methods for time delay estimation are based on Pade approximations which is often
convenient as the approximation of the time delay is given in term of a rational transfer function
instead of e~ TdS
. The approximation can be obtained from the following Pade approximation
of order n
- s
Td „ 2 - T s + (-T s) /2 + ( - r a ) / 3 ! + .. +
d d
2
(-T s) /n\
d
3
d
n
6
~ 2 + T s 4- {T sf/2 + (r a)3/3! + .. + {T s) /n\
d d d ' d
n
The reverse problem of the Pade approximation is to obtain the equivalent time delay of a
complex system described by

\ TA + «i« + a s + .. 2
G(s) = K —
2
—5 .
The equivalent time delay was derived by Matsubara [54] and is given by
D = b\ — a\
n m
i=l t=l
where pi and Zi are respectively poles and zeros of the transfer function G(s). This powerful
result was used by Isermann [44] to show that the following standard models
k ke~ TdS
G l ( s ) =
(TT7^ G 2 ( s )=
lTT7^ or
Gs(s) k
ES=i(i + ^ ) n
can be used for describing a wide range of different low-pass models.
These approximations or model reduction methods are useful for:
• improving delay estimation.
• Model reduction (reduction of parameter to be estimated).
In the following chapter, an alternative method for time delay estimation will be presented, as
an application of an off-line extended least-squares method.

Chapter 3
A n E x t e n d e d off-line Least-Squares M e t h o d a n d T i m e Delay E s t i m a t i o n
Most control loop performance assessment methods rely on some prior knowledge of the process
under investigation. This knowledge may be the process time delay or the complete process
model. Therefore, a major difficulty with these methods of assessment resides in the determina-
tion of the process model or the time delay. Because most performance assessment methods are
expected to work with normal operating data, it is not clear that it is possible to identify the
time delay in the absence of an external perturbation, e.g. setpoint change to ensure closed-loop
identifiability
In this chapter, it will be developed a technique that guarantees closed-loop identifiability
of a process without using an external signal, as long as the process disturbance is assumed to
be white. The case of colored disturbances although briefly discussed will not be studied here,
this is reserved for future work.
The organization of this chapter is as follows. In the following section, an introduction
to closed loop parameter identification is given. In section 2, some existing variants of the
least squares method are presented and an A R M A model is discussed for a regulatory closed
loop system in section 3. In section 4, the proposed extended off-line least-squares method is
presented. In section 5, the proposed method is used to determine the parameters and time
delay of a system. Section 5 draws some conclusions.
3.1 Introduction
In order to predict the future response of a dynamic process or to manipulate its outputs to fol-
low desired trajectories, a mathematical model which describes the dynamics of the underlying
22
Chapter 3. An Extended off-line Least-Squares Method and Time Delay Estimation 23
process is necessary. In some cases this model can be obtained by analyzing the internal mech-
anism or physics of the system. In other cases the internal mechanism is not well understood
or is highly complex, and an identification experiment has to be carried out.
The model may be linear or non-linear though, in practice, linear models are more common
since they are simple to analyze for identification, prediction or control purposes. Much work
has been done in linear identification and control, and the theory is well known and documented.
The estimation of the system parameters or equivalent of the process model can be performed
in open-loop by many methods. A detailed description of these techniques is given, for example,
in Ljung [51], Soderstrom and Stoica [62]. Under certain assumptions these methods can also
be used to estimate the parameters of a process operating in closed-loop. If the process is open
loop unstable, or the controller is adaptive, then it is preferable that the process model be found
by closed-loop identification.
Furthermore, it is now known that if the purpose of modeling is to design a controller, then
it is advantageous to perform the identification in closed-loop[50].
For these reasons, closed-loop identification has become an important research area, and
many methods dedicated to closed-loop identification are now available [28].
Most methods of closed-loop identification use a sufficiently exciting external signal added to
the control signal or to the reference signal. In Van den Hof et al.[67], the sensitivity function
is estimated first, then the process transfer function is determined. A recursive method is
proposed by Landau and Karimi[50] using the sensitivity function and an external signal, but
requires perfect knowledge of the controller.
The use of an external signal to achieve identifiability under feedback conditions was pro-
posed by Gustavsson et al.[34]. An alternative to using an external signal is to switch between
different regulators as proposed by Soderstrom et al. [62].
The use of an excitation signal is the best known way to guarantee the estimation of the
process parameters. However, for performance assessment, the use of such signals is rarely
allowed. In that case, it will be shown in this chapter that in the case of a regulatory process
where the acting disturbance can be described by white noise, it is possible to estimate the
process parameter without an external excitation signal.
In the following an overview of least squares methods is given.
3.2 Least squares parameter estimation
In the field of parameter estimation the least squares technique has reached a significant level of
popularity and perfection. This technique was first introduced by Karl Friedrich Gauss in 1795
who used it for astronomical computations. Some variants of the least-squares method will be
reviewed based on a survey article by Strejc[64]. More details and references are to be found
in that article. The objective of this section is to show the difference between the proposed
extended least-squares, which will be presented in the following section, and the existing least-
squares approaches.
3.2.1 Least-squares method
Consider the problem of estimating the polynomials A(q~ ) and B(q~ ) of the system described
1 l
by the following equation
A{q- )y(t)l
= B{q- )u{t)+e{t),
l
(3.31)
where
A(q~ ) l
= l + aiq- l
+ ...+a q- ,
nA
nA
(3.32)
Biq' )
1
= hq' 1
+ ... + b q-
nB
nB
(3.33)
y(t), u(t) and e(t) are respectively the process output, the process input and the process dis-
turbance which is assumed to be a sequence of independent and identically distributed random
variables with zero mean and variance a . Equation (3.31) is also equivalent to
2
y(t) = <p{t)0 + e(t) (3.34)

where
ip{t) = [-y(t - 1)... - y(t - n ) u(t - l)...u(t - n )}

A B (3.35)
6 = [ai ... anA bi ... b ) .

nB
T
(3.36)
Given the measurements y(l), y(t) and the estimates £(1), e(t) of e(l), e(t) the
least-squares estimate 9 of the vector 9 is determined by minimizing the following loss
function
N
V(9) = We\t) =ee T
(3.37)
where e(t) is the equation error defined by
e(t) = y(t) - v(t)9. (3.38)
The estimate 9 is given by
0 = {^NYÎYN (3.39)
if 4>%<J>N is invertible, which expresses the excitation condition (Astrom and Wittenmark[3]),
where
' y(i) ' V(l) '
YN = and (J>N = (3.40)
y(N) <p(N)
Some statistical properties of 9 are (Astrom and Wittenmark[3]):
E{9) = f5and (3.41)
cov(9) = <%{f <t> )-\ N N

(3.42)
Recursive equations can be derived for the case when the observations are obtained
sequentially as follows. Let 9(N) denote the least-squares estimate based on 7Y measure-
ments, and assume that (<frJf(f>N) is a

non-singular matrix. The least-squares estimate 9
then satisfies the recursive equations [64]:
§(N + 1) = 9(N) + K [y(N N + l)-<p {N T

+ l)9(N)} (3.43)
K N = P <p{N + N + <p (N + l)P <p(N + 1)]"

T
N
1
(3.44)
PN+1 = [I-K <p {N N

T
+ l)]P . N (3.45)
This method is based on the assumption that e(t) is white noise.
3.2.2 Generalized least-squares
The procedure of the generalized least squares considers the case when the additive noise
is a colored noise i.e. the mathematical model of the process is described by
A(q-i)y(t) = B( -i)„(t)g + J ^ L (3.46)
- Biq-ûty + wit). (3.47)
where A(q~ ) and B(q~ ) are as before and

l l
Dig- )
1
= l + d - A...
iq
l
+ d q- , nD
nD
(3.48)
and w[t) = £>(g-i) is the colored noise. w(t) can also be expressed as:
w(t) = -diw{t - 1) - ... - d w{t - n ) + e(t). nD D (3.49)
Multiplying both sides of equation (3.46) by D(q~ ) yields 1
A(q- )y {t)
l
f = B{q- )u (t) l
f + e{t), (3.50)
where
y (t)
f = Diq-^yit) and (3.51)
u (t)
f = D{q- )u(t).
l
(3.52)
If the autocorrelations of w(t) are known, then, define the covariance matrix W of
the colored noise w(t) as:
E[w(t)w(t)] E{w(t)w(t - n )} D
E[w(t)w(t - 1)] E[w(t - l)w(t - n )}

D
W =R W = (3.53)
E[w(t)w(t — n )\ D ... E[w(t — n )w{t — n ))

D D
The inverse of this matrix can be expressed as
W~ = V V,
l T
where V being triangular matrix. The estimate 9 of the vector 9 is determined by-
minimizing the following loss function
V{8) = (Ve) (Ve).

T
(3.54)
In this case it can be shown [64] that
e= {{V(i> yv<t> )-\v<i> YVY .

N N N N
• If the polynomial D(q~ ) is known then 9 can be determined using the least-squares
1
method based on filtered data yj and uj of equations (3.51)-(3.52) above.
• On the other hand when prior knowledge of the noise is lacking, the polynomial
D(q~ ) has to be estimated first.

1
Clarke [15] suggested the following off line estimate D given by
D = -(E E)- E HT 1 T
(3.55)
where
e(N) = y'N)-(p(N)§ and (3.56)
S = [e(n + 1) e(n
D D + 2) ... e(N)} (3.57)
and
e(nc) ... e(l)

E (3.58)
e(N-l) ... e{n -N) D
The filtered sequences
no
yf(t) = y(t) + Y, iy( - ) d t i (3.59)
i=i
Uf(t) = u(t) + E diu(t — i) (3.60)

i=l
are used for calculation of the estimate
9
= [Myf> f)<l>N{yf,Uf)]
u
(f) (yf,u )
N f Y (y ,u )
N f f (3.61)
which can be considered as the first approximation of 6 in equation (3.56). This procedure
is repeated to improve the estimation of polynomials A, B and D. It is noted that the
calculation of 6 gives A and B while D is calculated separately. Therefore, the estimation
of D depends on that of A and B and vice-versa.
A recursive solution exists and can be found in [64].

3.2.3 Extended least-squares
In the extended least-squares method the mathematical model considered is given by
A(q-')y(t) = B(q-i)u(t) +^ê(t). (3.62)
For simplicity we consider D(q~ ) = 1, therefore the model considered is

1
Aiq-^yit) = B{q- )u(t) + C{q- )e{t),

l l
(3-63)
where A, B and C are polynomials in the backward shift operator q~ . The estimation of x
9 can be done as before except the changes in the arrays 0and <p which are now expressed
as follows:
<p(t) = [—y(t — 1) ... — y(t — TIA) u(t — 1) ... u(t — n ) e(t — 1) ... e(t — n )] and
B c
9 = [ai ... a nA bi ... b ng c ...

t c ].
nc
T
The calculation of the error estimate e(t) by means of past estimates of parameters 9
which is the same, by means of the following relation
e(t) = y(t) - <p(t)9. (3.64)
The recursive algorithm for this method is similar to that of a recursive least squares
except the calculation of the prediction error. It is given in the following.
§(N + 1) = 9(N) + K [y(n+l)-(p {N N

T
+ l)9{N)} (3.65)
K N = P ^ ( N + l)[l + (^ (A^+l)Piv¥'(^ + l ) ] "

N
T 1
(3.66)
P^ + 1 = [I-K <p (NN

T
+ l)]P N (3.67)
e(N) = y(N) - (p(N)9(N) (3.68)
The difference between these three methods resides in how the prediction error, which
is an estimate of the noise, is taken into account. In the following a modified version of
the extended least squares method for a regulatory closed loop system is proposed.
3.3 A n A R M A model for a closed-loop system
In many cases, a regulatory closed-loop system can be described by an A R M A model.
Let us consider the plant model described by the following A R M A X model.
Aiq-^yit) = Biq-ûW + Ciq-êit), (3.69)
where
Aiq' )
1
= l + a q- + ... + a q- *,
1
1
nA
n
(3.70)
B{q~ ) l
= b^- 1
+ ... + b q~ nB
nB
and (3.71)
Ciq' )
1
= l + c q- + ...-rC q- .
1
1
nc
nc
(3.72)
y(t), u(t) and e(t) are respectively the process output, the process input and the process
disturbance which is assumed to be a sequence of independent and identically distributed
random variables with zero mean and variance cr . e
The process is assumed closed-loop stable under the feedback control law given by
«(*) = -f^foW-r), (- )
3 73
where r is a constant reference.
Replacing u(t) in equation (3.69) by the expression given by equation (3.73) yields
B(q-')R(q-i)
y(t) =
A(q-i)S(q-i) + B(q-i)R(q-i)
Ciq-^Siq- )- 1
A(q-l)S(q-l) + B(q-l)R(q-l)
with
B(q^)R{q- ) 1
A(q-l)S(q-l) + B(q-l)R(q-l)
This equation can be rewritten as
C(q-i)S(q-i)
e(t) (3.74)
A( -i)S(q-i)
q + B(q-i)R( -i) q
(3.75)
where y(t) = y(t) — y(t). This equation shows that the closed-loop system described
by equations (3.69) and (3.73) can be represented by an autoregressive moving-average
model ( A R M A ) defined by equation (3.75).
Therefore it is then proposed to estimate the past noise sequence e(t) and use it as
available information to estimate the process parameters A and B. The estimate of the
past noise is obtained once the parameters of the A R M A model of equation (3.75) are
determined. In the ideal case
Ciq-^Siq' ) 1
and (3.76)
Aiq-^Siq-^ + Biq-^Riq- ). 1
(3.77)
A R M A models are widely used in the field of time series analysis because they are a
general representation of a linear, stationary, stochastic process. Their parameters can
be determined in a number of ways described in Box and Jenkins[13], Ljung[51] or in
Astrom and Soderstrom[5]. In this last reference, the maximum likehood method is used
to estimate the polynomials ct{q~ ) and P(q~ ) which are defined as follows:
l l
a(q~ ) l
= 1 + c^q- 1
+ ... + a q- na
na
and
Piq- )
1
= l+p q-
1
1
+ ...+P q- i>.
n(i
n
The cost function to be minimized is given by
(3.78)
where the residual e(t) is a function of the observations y(l), y(2), y(t) defined by
a(? )e(*) = kQ~ )y(t)

_1 l
(3-79)
where
&{q~ )
l
= l + + ... + a q~ n&
n&
and
•3{q- ) =1
l + jhq- 1
+ ... + Pn q- ».
fi
n
d ( 9 ) and /?(g ) are assumed to be stable. Astrom and Soderstrom [5] showed that a,
_1 _1
/3 which minimize V(a,(3) satisfy
<% ) _1
aiq' ) 1
(3.80)
and
al = E{e {t))>E{e\t))
2
= ol (3.81)
where E(.) denotes the expectation function, given that
ria > n a and > rip. (3.82)
The equality o — o is obtained if n„ = n

2 2
a and rip = rip. It is also shown that the
residual e(t) is a good estimate of e(t) satisfying
E(e(t)e(t -k)) = 0 for k > 1. (3.83)
The estimation of e(t) can be very accurate. In the following, this accuracy is illustrated
by an example.
Example 3.3.1 ]52] The process considered is described by
_ 0.2q-'u{t) (l-0.5q-' + 0.6q- )e(t), 2
V { t ) _
(1 - 0.89") 1 +
(1 - 0.8(ri) 1
'
where the process disturbance e(t) is a sequence of independent and identically distributed
random variables with zero mean and variance o~ = 1. The control law is given by
e
^ = ^ 6 - 0 ^ ( 3 g 5 )
w
0.2 - 0.2?- 1 v y K
" K 1
The estimated A R M A model for this process is given by ^ p r j where
a'q- )
1
= 1 - 1.3955g + 0.3835g~ + 0.4585?- - 0.7784?- + 0.4394g~ - 0.0997g~
-1 2 3 4 5 6
p(q-i) = l - 1.1387?- - 0.3336?- + 0.8466g" + 0.2867g- .

1 2 3 4
The residual e(t) which is an estimate of e(t) is obtained by
e(t) = - !)• (- )
3 86
The variance obtained for e(t) is of = 1.0186 which is very close to cr = 1. 2

Its
•y tj* 5- <b © o o e> © o ©
Figure 3.3: Covariance of the residual.
autocorrelation is shown in figure (3.3). It can be seen that the condition of equation
(3.83) is satisfied. Figure (3.4) shows the near perfect match between the residual (dashed
line) and the white noise (solid line). Clearly the estimate e(t) oie(t) can be very accurate,
which confirms what is already known from the references cited above. For these reasons,
in the following, an off-line extended least-squares method to estimate the parameters
of a process under closed-loop is proposed. The proposed technique is a batch method,
and is based upon the estimate e(t) of e{t) as a first step, followed by estimates of the
process parameters.
2 .O I 1 1 1 1 1 1 r
- 2 .S 1 1 1 1 1 1 1 1
—
1 1 1
eoo eoo eio eio 020 e2S

Tim*
eao eas 84o e-as eoo
Figure 3.4: White noise and its estimate.
3.4 An off-line extended least-squares method
As mentioned above, it is proposed to estimate first the past noise which is obtained
by estimating the residual e(t), and then to use it as prior knowledge to determine the
process parameters. This method is restricted to processes that can be described by the
following equations:
y® = f f e l T ^ ) + (3-87)
<t) = -|^l(y(t)-r), (3.88)
where r is a constant reference (regulatory control), e(t) is a white noise and d is the
process time delay.
In the following it is first considered the case of d = 0 and r = 0. B y defining a
new variable yd(t) = y{t) — e(t) the following set of equations is deduced from the above
equations
Vd{t) = fjjrry«(*) (3-89)
u{t) = -j^j(yd(t) + e(t)). (3.90)
These last two equations represent a closed loop system where the reference signal is
the white noise —e(t) and the output signal is yd(t). Because —e(t) is a white noise we
are guaranteed to satisfy the closed-loop identifiability condition for this loop. There
we can approximate yd(t) as y(t) — e(t), where e(i) is an approximation of e(t) obtained
above. Therefore, it appears that there is no need for an external signal to be applied to
determine the transfer function which relates the input u to the output y . In this case d
both processes represented by (3.87)-(3.88) and (3.89)-(3.90) have the same closed loop
transfer function. The estimation of the polynomials A(q~ ) and B(q~ ) is done as will l 1
be described in the following. Equation (3.89) is also equivalent to
y (t) - ip(t)d
d (3.91)
where
ip(t) = [-y (t-l)...-y (t-n )u(t-l)...u{t-n )}

d d A B (3.92)
9 = [ai ... a nA h ... b ) .

nB
T
(3.93)
Given the measurements y ( l ) , y(t) and the estimates e(l), e(t) of e(l), e(t) the
variable yd(t) is reconstructed from the measurements of y(t) and the estimates e(t) of
e(t) as:
y (t)=y(t)-e(t),
d (3.94)
then the extended least-squares estimate 9 of the vector 9 is determined by minimizing
the following loss function

i N
(3.95)
L
t=i
where e(t) is the equation error defined by
e(t) = y (t) - <p(t)9.

d (3.96)
The estimate 9 is given by
9= (^V)- ^ 1
(3.97)
if (j) (j) is invertible. This is an excitation condition see Astrom and Wittenmark[3], where
T
Vd(l)
and (j) = (3.98)
Vd(N) ^(iV)
Some statistical properties of 9 are given by the following result.
Proposition 3.4.1 Consider the estimate 9 defined by equation (3.97). Assume that
e(t) is white noise and e is an estimate of e which satisfies equation (3.83) with zero
mean. Then the following properties hold:
1. 9 is the unique minimum point of V(9).

2. E(§) = 9.
3. cov{6) = o-l_ {4?4>)-\

e)
Proof: To prove this result, first write V{9) as
v{e) = hy -<t e] [Y -ct>e],

d )
T
d (3.99)
then setting the gradient of the loss function V{9) to zero yields
dV
(3.100)
O = - ^ = - F / 0 + 0 (0 0), T T
therefore 9 is obtained as in equation (3.97) which is also equivalent to

(3.101)
§ = (<j> ft- <l> (<f>9 +
T 1 T
e-e)
where
" e(l) ' ' e(l) '
e= and e = (3.102)
e(AQ e(N)
Therefore
0 = 0 + f> 0)-V (e-e), T

T (3.103)
which implies that E(6) = 9, because E(e — e) = 0. Equation (3.103) shows that if e is a
good estimate of e then (f/> <^) 0 (e — e) —• 0 and 9^9.

T _1 T
Rearranging Equation (3.103)
yields
9-9 = ((j) <f))- (t> (e - e)

T l T
(3.104)
Then taking the expectation of both sides of the last equation yields
E[(6 - 6)(6 - 6) ] = ol_ {<t> <t>r .

T
e
T l
(3.105)
The estimate 6 depends on the nature of the matrix (<j) <f)). T

This matrix should
be positive definite [4, 28] which is not always the case specially for process that is
under linear feedback of sufficiently low order. For such processes there exist some linear
dependencies among the columns of the matrix (b which means that the parameters can
not be determined uniquely [4]. The problem with lack of identifiability due to feedback
disappears if a linear feedback of sufficiently high order is used. Then the columns of the
matrix (b are no longer linearly dependent [4]. This, however requires the knowledge of the
true plant order, a practical impossibility. Another possibility is to consider an external
input signal that is persistently exciting. In Forssell and Ljung[28] it is shown that this
method leads to (<b (b) being a positive definite matrix. Furthermore, according to the
T
same reference[28]:"the general condition for (4> 4>) to be positive definite is that there
T
should not be a linear, time-invariant, and noise-free relationship between u and y. With
an external reference signal this is automatically satisfied ". Therefore, for the process
described by equation (3.89)-(3.90) the resulting matrix {(b (b) is positive definite. In
T
deed, it can be viewed as a closed loop system where the reference signal is a white noise
and also because in the expression of u(t) = — g(g-i) A e(t)) there is e(t).
The following section uses the method introduced here to estimate the parameters of
a process with time delay.

3.5 Off-line time delay estimation
The process considered is described by the following transfer function where the notations
and assumptions are those of the previous sections.
V(t) = | | ^ 9 - w ( t ) + e(t),
d
(3.106)
d is the process time delay for which a maximum and a minimum value are known i.e.
d g [d in dmax]. A(q~ ), B(q~ ) and yd(t) are as defined in the last section. The control
m
1 1
law is given in equation (3.73). Equation (3.106) is equivalent to
y (t)
d = tp(t, d)6
where
<p(t, d) = [-yd(t - 1) ... - yd(t - n ) u(t - d) ...u(t-d

A - n )}
B and
9 = [ai...a bi...b,
T
nA
Assuming the same conditions as before except for the presence of d, the loss function is
V(e,d) = lJ2e (t), 2

(3.107)
z
i=i
where e(t) is the equation errors defined by
e(t) = y -<p(t,d)0.
d (3.108)
The estimates 9 and d are obtained from the following minimization problem
[9,d\= min V(9,d). (3.109)
This problem is solved in two stages. First solve equation (3.109) for a given d g
[d i d ax] to obtain an estimate 9 defined by

m n m d
9 = {<t> \d) j {d)y f{d)Yd

d
1
( )
l
(3.110)
if (b (d)(b{d) is invertible, where Yd and (b(d) are defined by

T
VdiX) <p{l,d)
Y =
d and (b(d) (3.111)
y (N)
d <p{N,d)
then the estimates 9 and d are those which minimize V(9,d) for all d e [d i d \.
m n max
Proposition 3.5.1 Consider the estimated process time delay d and the estimated 6
defined by equation (3.110). Assume that e(t) is white noise and s(t) is an estimate of
e(t) which satisfies equation (3.83) with zero mean. Then [d, 9] is the unique minimum
point of V(9,d).
Proof: For a fixed d, it is shown in the previous section that 9d defined in equation
(3.110) is the unique minimum point of V(6d). d can only take values between d i and m n
dmax, therefore d which minimizes V(9 ) exists and is included in [d

d min d \. •
max
In the following, the results of this chapter are illustrated by the following example.
Example 3.5.1 [52] The process considered is described by
y(t) (3.112)
where the process disturbance e(rj) is white noise of zero mean and variance a = 1. The
e
control law is given by the following proportional controller:
u(t) = l-y(t). (3.113)
The first step consists in determining the closed loop A R M A model or the filter H(q ) :
which will be used to estimate the noise sequences. The estimated A R M A model for this
process using A R M A X command of M A T L A B is given by H(q *) = ^ p r j where
a{q- )
1
= 1 + 0.0843?- - 0.24939
1 -2
and
fitq- )
1
= 1 + 0.1132?- - 0.2359g - 0.0504g- + 0.0208?-
1 -2 3 4
- 0 . 1 4 6 9 g - 0.1466g - 0.1179g .
-5 -6 -7
Figure 3.5: Estimation value of the coefficient a.
Using 500 different noise sequences, the estimates value of the parameters a and b
are shown respectively in Figure (3.5) and Figure (3.6). The estimated values of the
time delay are plotted in Figure (3.7). Figure (3.8) shows the values of the loss function
V(6,d). The average estimated model is given by:
(0.1842)?-
y(t) = (3.114)
It can be seen from comparing this model to the one described by equation (3.112)
that the time delay is exactly identified and the process model reasonably approximated.
0.2
0.1051-
0.165' 1 1 1 1 1 1 1
1 1 1
O SO IOO 1 SO 200 2SO 300 35a 400 -J50 SCO
N u m b e r o1 & l m u l a 1 i o r m
Figure 3.6: Estimation value of the coefficient b.
e
SB -
5JS -
fc*5.4
-
£
a 5
9
.a -».B
-
•w
E
SU*
4.4 -
42
O 50 IOO ISO 200 2SO 3 OO 350 400 4SO BOO

N u m b * r o1 & l m u l s ' lo n&
Figure 3.7: Estimation value of the time delay d.

o.ozs —
0.0£4 -
0.023 -
O SO TOO I SO 2CO ZSO 300 3SO -*00 -ISO BOO

N u m b e r o l » l m u l a l l o rm
Figure 3.8: Loss function V .
Applying ordinary least-squares to the example above, and by considering 500 differ-
ent noise sequences, the estimates value of the parameters a and b are shown respectively
in Figure (3.9) and Figure (3.10). T h estimated value of the parameter a is completly
wrong. The estimated values of the time delay are plotted in Figure (3.11). The average
estimated model is given by:
V(t) = ( f l ^ ^ ^ + e^), (3.115)
It can be seen from comparing this model to the one described by equation (3.112) that
the ordinary least-squares method has difficulty to approximate the process model. The
following table summarizes the obtained results from using the proposed method and
ordinary least-squares.
Real values a=-0.8 b=0.2 d=5

Parameters a 4 b oi d
Extended least-squares -0.775 3.85 IO" 5
0.1842 1.99 10~ 5
5 0
Ordinary least-squares -0.24 5.68 10~ &
0.1842 3.2 IO" 4
5 10.51
Table 3.1: Results.
Chapter 3. A n Extended off-line Least-Squares M e t h o d and Time Delay Estimation 44
- 0 . 2 2 \-
—0.2T7 ' " " ' 1 1 1 1 « I I

O SO IOO 1 SO ZOO 2SO 3QO 3SO 400 450 500
Number oi ilmuta'iorm
Figure 3.9: Estimation value of the coefficient a by ordinary least-squares.
0.2B
0.24 \-
0.14' 1 1 1 1 " 1 — 1 1 1 1
O SO IOO 1 SO ZOO ZSO 300 350 400 A SO SCO
N u m b s r o l « i m u l s l l o rm
Figure 3.10: Estimation value of the coefficient b by ordinary least-squares.

350 -*00 -*SO SCO
Figure 3.11: Estimation value of the delay d by ordinary least-squares.
3.6 Case of colored noise
The case where the disturbance is colored is more complex to solve, as knowledge of the
disturbance dynamics is required to compute yd = y — w. Consider the following process
= T7?^9~V*) + (- )
3 116
(3.117)
where r = 0 is a constant reference (regulatory control), w(t) is a colored noise i.e.
w(t) = H i (q~ )e(t)

no se
l
and d is the process time delay. Equation (3.116) and (3.117) are
also equivalent to
y (t)
d = y(t)-w(t) = ^^q- u(t) d
(3.118)
u{t) = -j^(y ,(t)-w(t)).

d (3.119)
These last two equations are similar to equation (3.89)-(3.90) except that w(t) is colored
noise. To apply the previous analysis, it is needed to estimate e(t) as described above
and determine the filter Hnm^q" ) then reconstruct the two signals:
1
w(t) = H (q- )e(t)

noise
l
and (3.120)
y (t)
d = y(t)-w(t). (3.121)
However, if no external signal is considered, the determination of the filter H (q~ )

noise
l
can only be obtained if the process considered operates for a while in an open loop man-
ner. Therefore, the general case of non white disturbance is not solved by the proposed
method.
3.7 Conclusion
The method presented here considered the case of a regulatory process where the poly-
nomial C(q~ ) = 1. For that special case it was shown that the process is identifiable
1
without an external output signal. The general case though briefly discussed was not
solved.
The following chapter addresses some problems related to the use of minimum variance
control as benchmark in assessing the performance of controllers.

Chapter 4
Realistic Performance Assessment in the SISO Case
4.1 Introduction
Harris's approach [35] for control loop performance assessment is simple and effective.
The method requires routine operating data and little prior knowledge of the process.
When using this method to assess control loop performance, problems may arise, some
of which are addressed in this chapter.
The organization of this chapter is as follows. In section 2, two existing performance
index definitions are used to illustrate some of the problems related to the use of this
method. In section 3, satisfactory versus unsatisfactory performance is investigated and
the transition between them is established as the point at which the process starts to
oscillate. In section 4, a tuning procedure to bring the output variance to an acceptable
level is described and shown to be a trade-off between other criteria of interest.
4.2 Problem formulation
The comparison of the output variance to the minimum output variance a^ v can be
carried out by using one of two performance index definitions encountered in the literature
[18, 35].
Here the performance index is defined as:
(y-r)
1 < I < oo r is the setpoint.
p (4.122)
mv
47
Chapter 4. Realistic Performance Assessment in the SISO Case 48
When the actual performance is close to the best achievable performance, I fa 1. When
p
the actual performance is significantly worse than optimum, I p >> 1. To bound the
performance measure between zero and one, a performance index can also be defined as
V
In the definition used here, only the extreme situations are well defined: good variance
performance when the performance index is equal to one and poor performance when the
index is large. In practice, the index takes on values between these extremes, and it is
necessary to determine whether the actual controller performance is acceptable, accept-
able with a small tuning correction, or so unacceptable that redesign of the controller is
needed.
As discussed in the introduction, the minimum variance control strategy has some
limitations; therefore it is not always useful to implement. However, a good control
strategy is often obtained when compromising between output variance as a performance
criterion and other measures such as steady-state rejection of disturbances. To illustrate
this idea, consider the following example.
Example 4.2.1 The process is described by Figure 4.12 where the variance, cr , of the 2
white noise disturbance is equal to 0.36.
Figure 4.13 shows I and the steady state error S

p se ( defined as E[^f^} in % ) versus
the proportional controller gain K. It can be seen from that figure that for small gains
the output variance is relatively low but the steady state error is very high. For K > 1.5,
the output variance becomes high and the steady state error is reduced significantly. A
compromise between these two criteria may be needed if a proportional controller is the
only available design. Thus, although I has advantages as a performance index, a better
p
overall design is obtained when combining it with other performance criteria.

Ref
) k Controller
w
l-0.67q J
Figure 4.13: I and steady state error versus P-gain.

p
Recognizing that minimum variance performance is not always desirable, the next sec-
tion discusses how good operating performance can be distinguished from unsatisfactory
operating performance.
4.3 Satisfactory and unsatisfactory performance
From the performance index definition, poor performance is obtained when the perfor-
mance index takes large values. This extreme situation corresponds to a nearly unstable
loop. In reality, poor performance will be obtained long before this extreme situation
occurs and, for a process with no integration, only a small range of values of I actuallyP
represents realistic performance alternatives. In practice, the control signal inside the
loop is bounded between u i m n and umax and will saturate if values outside this range are
called for. Now, consider the single loop system, without integration shown in Figure
4.14. The nonlinear component (N. L.) is a saturated amplifier with a linear gain equal
to 1 represented by its equivalent gain ^ y . The remainder of the process is modeled by
a linear component represented by KP(jui). The saturation characteristic is shown in
Figure(4.15). The equivalent gain is defined[17] as
Ni(e) = ——— = — / f(ecos(u)t))(cos(cot)+jsm(ut))d(ujt).

e ire Jo
The stability of such a system can be investigated in the Nyquist plane, Figure 4.16, by
using the Describing Function method. When no intersection between KP(ju) (curve
1) and - ^ y occurs, the loop is stable and the performance can be assessed by use of the
output variance criteria. Increasing the gain K, the curve KP(ju) moves to the left to
intersect - ^ y (curve 2) for the first time, causing the loop to oscillate at frequency and
amplitude (cJn, e ) determined by solving P(JUJ) = Increasing the gain K causes
the loop to oscillate at a constant frequency, but with increasing amplitude (curve 3).
The conclusions to draw here are:

e(t)
lit) vlt)
N(E) p> KP(jffi)
Figure 4.14: A single nonlinear loop.
Figure 4.15: Saturation characteristic.

Chapter 4. Realistic Performance Assessment in the SISO Case
Figure 4.17: Partition of Ip curve.

1. in this case there are two situations which correspond to either an oscillating or a
non-oscillating loop,
2. the performance may or may not be considered acceptable up to the point I o, at p
which the loop starts to oscillate,
3. performance is definitely poor from Ipo up to IpUmit which is obtained when replacing
a saturated amplifier by a relayor k —• oo.
This analysis is not restricted to a saturated amplifier but can be extended to other types
of nonlinearity. These conclusions are summarized by Figure 4.17.
This approach provides a link between oscillation and performance. It recognizes
that whenever oscillation is present in a process, the output variance represents poor
performance. It is also noted that the beginning of unsatisfactory performance is obtained
at the latest for I = I o which is finite << co and corresponds to an oscillatory system.
p p
The loops causing these oscillations should be localized in order to remedy the problem.
In the following chapter the problem of localizing the oscillating loops will be addressed.
IpO plays a major role in some tuning methods because it defines a limit between
acceptable performance and unacceptable performance. The Ziegler-Nichols method, for
example, uses (indirectly) Ipo to determine the controller tuning parameters as functions
of the gain k and the period P .

u u The ultimate gain k and period P are obtained
u u
by setting integral and derivative control action to zero, and increasing the proportional
gain k until the measured variable oscillates at period P for a certain gain k (Shinskey
c u u
[59]).
Figure 4.18 from Fu and Dumont [30], shows an oscillating loop. Tuning was per-
formed at t = 27 minute following which the oscillation disappeared and the performance
index dropped to a lower value.

Mill Trial Result:

A Dryer Pressure output
Controller Tuning was performed

based o n the evaluation performance
Figure 4.18: Performance analysis.
4.4 Control chart
The curve of figure (4.17) can be used for supervisory purposes. Given IpQ, determined
during the tuning procedure and the estimate cr^„ of the minimum achievable output
variance, it is proposed to rescale the calculated performance index as given in the fol-
lowing.
I (%)
P ^-^100. (4.123)
I pO
This new performance index is expressed in percentage (% ). This establishes a new scale
for I , which can be interpreted as how far the actual controller performance is from that
p
of a minimum variance controller and how close it is to the worst performance. The
alarm should be given as soon as the performance reaches the range 85 ~ 100.
The performance index within the non-oscillatory region (figure 4.17) is satisfactory
only if it results from a useful trade-off between all criteria considered. In the next
section, a tuning procedure is presented.

\
4.5 Tuning procedure
Before compromising between the criteria it is necessary to estimate the useful portion
of I curve, corresponding to I G [1 Ipo], in terms of the controller parameters.

p p Note
that I o can be estimated by the relay method.

p
Although I itself may be estimated without knowledge of the process (except any
p
delay), the variation of I with controller parameters can only be predicted once process
p
dynamics have been estimated.
To estimate the I curve, a design method is proposed in which:

p
1. The process is approximated, Figure 4.19, by one of the following models

o-T s d
M ^ s ) = f ^ , M (s) = ~ 2 or M ( s ) = 9 6 T d S
3
1 + rs' s + as + b 2
s(s + as + b)' 2
•TdS
ge
l+TS
•T s
ge d
s + as +b
2
ge- TdS
s(s âs+b)
Noise
Ref
C ontroller Process
Figure 4.19: Process approximation.
This choice covers the most common process behaviours encountered in process
industries [44]. The model chosen is that minimizing the error between its output
and that of the process, i.e., minimizes the cost function:

N
1 il2
J{Mi)
The models can be obtained using the method of chapter 3. They can also be
obtained by estimating the time delay and dynamics separately:
• the time delay Td is estimated by one of the two techniques described in
Elnaggar and Dumont [22] or in Zheng et al.[71].
• g, T, a and b are estimated for each model using a least squares method
assuming sufficient excitation.
2. a noise model is then estimated by one of the two following techniques: Time series
analysis[35] or Laguerre network[52].
3. the controller is adjusted in order to obtain the I curve and corresponding curves
p
for other performance criteria of interest.

e(t)
l-OV
l-O.CTcf 1
In considering these curves one can:
• establish the compromise controller tuning parameters,
• be guided by these curves to define "satisfactory" performance, and
• decide the direction in which the controller should be adjusted.

The controller which results from this procedure is, therefore, a trade-off between all
criteria considered.
The following section illustrates this method by means of an example.
4.6 Simulation Example
The process considered is shown in Figure 4 . 2 1 . A proportional controller is taken in the
first test, and a proportional plus integral (P.I) in the second test. The estimation of the
I curve for the interval [1 7 ] is carried out following the approach described in section
p p o
4. The following model provides the closest representation:
0.232g- 4
, , 1 - 0.5339- 1
, .
y(t) = 1 - 0 . 6 2 3 4 ? " u(t)
T H — e(t). 7
y\ > 1 ncoo^ -i w -r 1 _ 0.62349- 1 w
Ref 0.33(1-0.6^^-*
Controller
A - (1-0.6q-)(\-0.2ci)
l l
Figure 4 . 2 1 : Process control.
F i r s t test: Cbntroller=P-gain i.e. proportional only. The two criteria, output variance
and steady state error S er are considered (S er is defined as E[ ' ^] r ef

in %). Figure 4 . 2 2
show the simulation results. The solid curves represent the true process, and dashed
curves are obtained using the method described above. It can be seen that the estima-
tion of I p curve and steady states error on the basis of the simplified process model is
very close to true values, especially at low gains. From these curves, the best choice of
tuning parameter can be determined as a compromise between steady state response and
disturbance rejection.
Figure 4.22: Comparison of true and approximate model.
Second test: Controller^ P i.e. Proportional plus Integral control. Interest is in
the process step response, and in output variance. I , tr (tr is defined as the time needed
p
for output y(t) to reach 90% of y(oo)) versus controller gains for the actual process, and
its approximation is shown in Figure 4.23 and Figure 4.24.
From these curves it can be seen that:
• The estimated curves are close to the true ones and in this case lead to satisfactory
controller design.
• Combining the curves Figure 4.23-4.24 one can establish the best compromise con-
troller tuning parameters. These curves also indicate the direction in which the
controller parameters should be changed to improve a certain performance crite-
rion.
• Similar curves may be obtained for any controller candidate with a small number
of adjustable parameters.
Figure 4.23: Estimated curves versus Controller parameters .
Remark: This technique uses the minimum variance o mv = (1 + ft + ... + fx -i)4

d to
normalize curves of the output variance. Note that although incorrect delay estimation
may lead to one or more terms, such as f^-i^l being added or missing in the estimation
of o ,
mv the effect is to move the I curve up or down. This estimation error, thus, need not
p
invalidate the method since the choice of the tuning parameters is always a compromise
among all the criteria considered, and the actual value of Ipo is usually less important
than its variability.

O O .3 1 1 -S 2 O O .3 1 1 .3 2
Figure 4.24: Estimated curves versus Controller parameters .
4.7 Conclusions
In this chapter it was linked oscillations to poor performance and shown that by broaden-
ing Harris's performance index as a measure of output variance to consider other criteria
a better overall assessment of control loop performance can be obtained. A n approach
to controller tuning based on a combination of performance criteria and output variance
was also proposed.
Because oscillations amputate the control loop performance, oscillating loops should
be localized and controllers should be tuned. In the following chapter the problem of
localizing the oscillating loops is investigated.

Chapter 5
Monitoring Oscillations in a Multiloop System
5.1 Introduction
Poor performance of process control in the pulp and paper industry may be attributed
to three causes:
1. The first one is due to nonlinear behavior of control valve and sensors [9]. For
example at the Mead Fine Paper mill in Chillicothe, O H , about one loop in five
cycles at steady state[23] because of nonlinearity in control valves.
2. The second cause is simply use of the control strategy or poor process design.
3. The third cause results from the fact that most mills have a great number of control
loops and a very limited staff of instrument technicians and maintain to check the
adequate functioning of each loop.
A cycling loop affects the performance of others as it spreads the oscillations to its sur-
rounding loops; therefore, it is critical for overall loop performance to detect and eliminate
oscillations. To illustrate this fact, consider a loop where the process is described by
(1 - 0.74 - )y(t) = 0.26q- u(t)

9
1
A
+ (1 - O.bq' + O.Qq~ )e(t)
1 2
and the control law is given by
0.5 - 0.37?- 1
u(t) = 7(r(t) - y(t)),

w
0.26-0.13?- - 0 . 1 3 g -
1 4 V y ) y y > h
61
Chapter 5. Monitoring Oscillations in a Multiloop System 62
where r = 5 is the setpoint or reference signal, y(t) is the process output variable and e(t)
is a white noise sequence (^(t) =

0-81) • Assume that this loop is oscillating at frequency
OJ = 0.314 rad/sec as result of its interaction with an oscillating loop. This situation
0
is simulated by adding a sine wave disturbance signal equal to 2.5sm(0.6287r£) at the
output of the loop. The output signal y(t) can be expressed as:
y(t) = y + % ) e ( £ ) + A sin{u t
_ 1
0 0 + <f> )
0 (5.124)
where ft(<? ) is the transfer function which relates e{t) to y(t). A n estimated A = 3.207
-1
a
is obtained using equation (5.127)-(5.130). The minimum achievable output variance is
°~mv =
1-145. The actual output variance is o ^ 2
= 7.84. The performance index is
Ip{v{t)) =
6.877. Based on this value of I , the performance of the actual controller
p
may be judged poor. However, this loop is not generating the oscillation. Then, by
subtracting the contribution of the frequency u> to the output variance, which may be
0
A2
estimated as -f = 5.1424, the estimated output variance without the frequency u would 0
be
A2
_ 2 O _— r, n
°~(y(t)-oscillation) — y(t)
a
^~ --
Z D
Therefore, a more realistic performance index for this loop would be I i = 2.3, which is rea
significantly less than the one calculated before.
This example shows how the variability of a non-oscillating loop is increased when
interacting with other oscillating loops. It also shows that not taking into account oscil-
lations may yield a wrong performance index calculation. Therefore, it is of interest to
localize oscillating loops and, using an appropriate controller tuning method, to reduce
the amplitudes or eliminate the oscillations completely. Knowledge of the presence of os-
cillations in loops may be obtained by techniques such as Discrete Fast Fourier Transform
( D F F T ) or Singular Value Decomposition (SVD)[47].

The issue of localizing the oscillating loops in a cascade multiloop system is addressed
in this chapter. The organization of this chapter is as follows. In section 2, background on
the concept of equivalent gain for nonlinear systems is given. In section 3, an oscillation
index to help localize the oscillating loops, in a cascade multi-loop system, is define. In
section 4, the accuracy of using the describing function method is discussed, followed by
a simulation example, and then a conclusion.
Note that the analysis in this chapter is performed in the continuous time domain.
5.2 Background
A n industrial process consists of a combination of linear and nonlinear components.
When it is possible to separate the process into a nonlinear gain and linear dynamics,
then such process can be represented by Figure (5.1) where e(t) is a white or colored
noise. The behavior of such systems has been extensively discussed in the literature
[17, 10].
This investigation discusses, for any loop with (or without) nonlinear components,
conditions under which oscillations are present and, if so, the form of their characteristics.
One method commonly used to address these questions, for systems of the form of Figure
5.1, is the describing function (DF)[17], introduced during the late 1940s. The D F is an
approximation procedure that recognizes, for the simple feedback loop of Figure 5.1, that
when limit cycles occur the input to the nonlinearity N is almost sinusoidal due to the
low pass filtering action of a typical control system plant, G(s). The nonlinearity can
therefore be approximated by an equivalent gain or describing function, N(e), defined as
the complex ratio of the fundamental of the output of the nonlinearity to the amplitude
of a sinusoidal input at that frequency [10]. Assuming that the nonlinearity input is
e(t) = ecos(cui), its output is thus approximated by u(t) ~ uiCOs(tvt) + U2sin(ut). If
u = f(e) represents the nonlinearity input/output relation, the equivalent gain is given
by
Ni(e) = = — / f (e cos(ut))(cos(ut) + jsm(ut))d(ut).

e ns Jo
For the autonomous case (r(t) = 0), the prediction of the limit cycle rests upon finding
a solution to the equation
l + N(e)G(s)\ „ J = 0, (5.125)
which is equivalent to a unit loop gain for the fundamental component of the signals e(t)
and y(t) denoted respectively
Ef(t) = efsm(ut + <fi )

£f
Vfit) = yfsm(wt + <l)y ) f
where (f> — <f) = —180°, resulting in

yf ef
= -Vf(t)-
5.3 Monitoring the oscillations
The reverse problem, that of limit cycle prediction, involves recognizing which loop in a
multiloop system is causing oscillation at a particular frequency. For that, an oscillation
index is defined as follows:
Definition 5.3.1 For the linear feedback loop of Figure 5.2 where an oscillation of fre-
quency LU is detected, define the oscillation index, at frequency LU, by
/ 0 S C H = |1-||, (5.126)
where ^ is the gain of the considered loop at frequency LU. This gain is determined by
extracting the magnitude and the phase of the oscillation LU from both signals y(t) and
e(t) using the following set of equations:
y = y/yi+y*, (5.127)
y\ = j,J y{t) cos (ut)dt, (5.128)

2 [T
2/2 = j; J Q y(t) sin (Lut)dt, and (5.129)
<f> y = arctan(^). (5.130)

e is obtained from the following set of equations:
e = \Jel + sl (5.131)
2 fT
= — f (t) cos (ut)dt, (5.132)
T Jo
£ l £
2 r T
e = — / e(t)sm(ut)dt, and (5.133)

T Jo
2
<f> £ = arctan(—). (5.134)

£2
For a cascade multiloop system, where loops are defined by sets of a closed paths according
to Mason's loop rule, the following result is proposed:

e(tj
r(t)
rL[5)t<\ f r M
IJ\5 }
P r o p o s i t i o n 5.3.1 For oscillation at frequency to at the output of loop U in a multiloop
system, the loops associated with the oscillation u are those with /^ (cu) ~ 0. c
Proof: To prove this result, first consider the single loop of Figure 5.2. This loop is
oscillating at frequency u if ju> is a solution to the following characteristic equation
(5.135)
1 + C( )G( )U = 0.
5 S
This last equation is in fact a set of two equations which are
s=ju) — 1 gain equation, (5.136)
lC(s)G(s)\ ,s=JU = - 1 8 0 ° phase equation. (5.137)
The frequency UJ was detected at output of that loop; therefore, it remains to extract
the magnitude and the phase of that frequency from both signals e(t) and y(t) to find
out whether equations (5.136)-(5.137) are met. In fact this is what is done by equation
(5.126) and equations (5.127) to (5.134).

The equation (5.135) is a necessary and sufficient condition for any loop to generate
oscillations at frequency UJ; therefore, for the loop k taking from the considered multiloop
system, to generate oscillations at frequency UJ equation (5.135) should hold which implies
that Iii (UJ) ~ 0. •

osc
The extension of these results to a nonlinear multiloop system is proposed by the
following result.
Proposition 5.3.2 For oscillation at frequency UJ at the output of a nonlinear cascade
multiloop system, the loops associated with the oscillation UJ are those with /o (cu) — 0.l
sc
Proof: To prove this result, first consider the oscillating nonlinear feedback loop of
Figure 5.1, where the input r(t) = 0, and then consider the case of r(t) = RQ+RI sin(a;it).
The extension of these results to a multiloop system will be discussed in the following.
The autonomous case, r(t)=0
From the above, any limit cycle (UJ,O) satisfies Equation 5.125. To obtain the equiva-
lent gain for the nonlinear component, the D F approach of [10] is followed. Write the
approximation at the output of the nonlinear component as
U(t) Ui cos (cut) + « 2 sin (ujt) (5.138)
u sin (ut + (b ) u (5.139)
where
(5.140)
(5.141)
1 Jo
arctan (—) (5.142)
u 2
then, u is given by
(5.143)
u= y/ul + v%.
For the nonlinearity input a similar approximation gives
e(t) ~ £i cos(ut) + £ sm(cot) 2 (5.144)
= esm(ut +(j) ) £ (5.145)
where
£ l = —t £(t) cos(ut)dt (5.146)

T Jo
e 2 = | ; f e(t)sm{tot)dt
T
(5.147)
1 Jo
(j> = e arctg{—) (5.148)
£2
and £ is given by
£ = yjel + el (5.149)
The estimated equivalent gain N for (cu, a) is given by
\N\= ^^1^, Z0 -0e-

u (5.150)
For the linear component, because it satisfies the superposition theorem, the gain at co
can be obtained by applying the same reasoning to the fundamental component:
y(t) = yi cos(ut) + y sm(ut) + terms 2 (5.151)
= ysm(cot + (j) ) + terms, y (5.152)
where y\ and y are obtained by the same technique as above, y is given by

2
y= \fyJ+?2- (5.153)
The estimated linear gain at u> is given by
\G(u)\ = ^ A L^-K. (5-154)
It should be noted that the integrations considered above are affected by the presence of
noise, n(t), introducing the following integral terms
/ n(t)cos(wt)dt, (5.155)
Jo
[ n(t) sm(wt)dt (5.156)
Jo
in the expressions above. To minimize the contribution of the band limited noise terms,
the integral in equation (5.155)-(5.156) should be taken over a number of periods, kT,
instead of one period.
Replacing G and ./V by their estimated value in equation (5.125) yields, for this simple
case:
- = 1 (5.157)
s
I (u)
osc = | 1 - | | = 0. (5.158)
Equation (5.125) is a set of two equations: the magnitude equation which is used for
equation (5.157) and the phase equation:
</> -<f> = - 1 8 0 ° .
y e (5.159)
Equation (5.159) must be satisfied for a uniformly maintained oscillation OJ, both the
total gain must be equal to one, and there must be a total of 180° of dynamic phase lag
in the loop.
R e s u l t 1: If the oscillation to is generated inside the loop of Figure 5.1, then
Io8c{u) = 0.
The above analysis considered the case of an autonomous system. To generalize this
result, consider the cases of an external input, and of an external input, combined with
auto-oscillation.
Forced system
The input into the system of Figure 5.1 is now assumed:
r(t) = R + Ri sm(wt).
0 (5.160)
The most obvious possibility is that all signals in the feedback loop can be approximated
by
e(t) ~ eo + eis'm(ut + (j) ), £ (5.161)
and y(t) by
y(t) ~ G(0)N e 0 0 + e i | G f » | W ( e i ) sm(wt + </>').
Since
e{t) = r(t)-y(t),
it follows that
e(t) ~ £ + £isin(ut
0 + <j> ) = r(t) - y(t)
e (5.162)
~ R^ + R^m^-G^Noeo-exlGijLJ^Nie^smicot+ (/)') (5.163)
~ Ro - e G(0)No + yjR\ + [e^G^N{e^

0 sm(ut + 0,) (5.164)
and by equating coefficients results
R0 = e (l + G(0)N )
o o and (5.165)
R\ = (l-[|G0u;)|Ar( )] )4 £l
2
(5.166)
Since R\ ^ 0 implies [\G(JUJ)\N(EI)] 2

^ 1, the following result is established.
Result 2: If the oscillation UJ is not generated inside the loop then
The estimated total gain \G(JUJ)N(EI)\ is obtained from data by the procedure given
above.
Now, consider the case of an external input r(t) at frequency UJ into the loop of Figure
5.1 combined with internal oscillation at frequency UJ . Then the input into the system C
of Figure 5.1 is
r(t) = R + 0 R sm(ujt).
1
(5.167)
The most obvious possibility is that all signals in the feedback loop can be approximated
by
e(t) ~ €Q + e\sin(ujt + <j>i) + E sin(uj t + (b ) c c c
and y(t) by
y(t) = G{0)N e +\G{juj )\N{e )e sm{ujt

0 o 1 1 1 + f> )
( 1
+\G(juj )\N(e )e sm(ujt

c c c + <f> ). c
Rearranging this equation and equating coefficients as before results in:
Ro = e (l o + G(0)N ) o (5.168)
R\ = ( l - O G W I ^ e O H e ? (5.169)
0 = (l-[\G(juj )\N(e )] )el

c c
2
(5.170)
Repeating the same analysis leads to the following conclusion:
I c(u)
OS ^ 0 and
Iosc(u ) c = 0.
Thus the evaluation of the oscillation index for each detected oscillation frequency at the
output of a loop will distinguish between oscillation generated inside the loop and that
coming from outside into the loop.
Now return to the multiloop system, the condition for limit cycles at frequency UJ and
amplitude a to occur in a loop (i) is that the fundamental component of the signals y (t) l
and e (t)
l
Ef(t) = efsm(ut + <f) )ef
y){t) = y}sm(ujt + <j) ) yf
where <pyf — (b£f verify
E){t) = -y}(t). (5.171)
The reasoning in the above analysis completes the proof of the proposition (5.3.2). •
The approximation considered here can be applied still to a signal containing more
than one frequency [17]. Note that when a loop is externally excited, auto-oscillation
may be hidden. This behavior does not cause trouble once the remaining detected fre-
quencies have been assigned to loops. B y tuning oscillating loops in sequence, the hidden
oscillations will appear and the tuning procedure may be repeated until all loops are
tuned.
In the following section the effectiveness of the describing function method will be
discussed.
5.4 Accuracy of the Describing Function Approach
The number of terms that should be included in the oscillatory signal component, xap(t)
can be determined using the normalized Mean Square (MS) error given by
2 E{\x(t) - x (t)\ }
ap
2
6
- E{\x(tW} • ( 5
- 1 7 2 )
Taking e to have a value less than 0.05 implies that the approximation accounts for 95%
of normalized MS variation of x(t).
5.5 Some Limitations
The proposed method is limited to cascade loops. Note that a given loop may oscillate
at frequency LO either if ju> is solution to the equation (5.135) or as result of a periodic
disturbances at frequency LO acting on that loop. If periodic disturbances may be present
these two cases can not be distinguished directly with the present method.
5.6 Simulation example
disturbances
i»,(t) mi)
N.LI NI.2
sfl
To illustrated the proposed method of oscillating loop detection it is considered in
the following two examples.
Example 5.6.1 Consider the cascade control system of Figure 5.3 where the distur-
bances are white noise (aj isturbances = 0.4),
1
G(s) = and (5.173)
s + 2s + 4s
3
2
(5.174)
The nonlinearity N . L . 1 is described by
5 if ini(t) > 5
outi(t) k ini(t) if - 5 < ini(t) < 5
p
-5 if ini(t) < - 5
while the nonlinearity N . L . 2 is described by
5 ifin (t)>5
2
out (t)
2 k in (t)
2 if - 5 < in (t) < 5
2
-5 ifin (t)<-5
2
The process considered here, is formed from two loops. Therefore, if an oscillation is de-
tected at the output y(t), and assuming that it is not resulting from an external periodic
disturbance, which is not handled in this study, then this could not happen unless at
least one of the two loops is oscillating. In such a case, and provided the scheme of the
two loops, the oscillation index of the outer loop at the detected oscillation frequency
will be equal to one and that of the inner loop at the same frequency will be equal to
one only if it is causing the oscillation. Therefore, only the oscillation index for the inner
loop at the detected frequency will be calculated in the following two test.
Test 1: First, set linear gain k = 4 and k = 10. This corresponds to an oscillation
p
frequency at to = 0.589 rad/sec and a period of T = 10.667 sec. The oscillation is caused
by the outer loop. The output of the outer loop is shown in Figure 5.4. The output
of the inner loop is shown in Figure 5.5. The output of the nonlinearities N . L . 1 and
N . L . 2 are shown in figure(5.6) and figure(5.7) respectively. From these figures it can
be seen that the oscillation drives the system into saturation regions for N . L . 1 and N . L . 2.
The estimated total gain from the input i n to the output y at oscillation to = 0.589 rad/sec,
2 2
is
total gain (0.589) = 0.67. (5.175)
Figure 5.4: The output y(t) of the outer loop, example(5.6.1)-test 1.
Therefore the oscillation index for the inner loop is
7 (0.589) = |1 - 0.67] = 0.33 f 0.

OSC (5.176)
It can be concluded that the inner loop is not generating the oscillation at frequency
to = 0.589 rad/sec. It is generated by the outer loop.
Test 2: Now consider k = 10 and k = 0.5. This corresponds to an oscillation frequency

p
at LO = 2.47 rad/sec period of T = 2.53 sec. The oscillation is caused by the inner
loop. The output of the outer loop is shown in Figure 5.8. The output of the inner loop
is shown in Figure 5.9. The output of the nonlinearities N.L. 1 and N . L . 2 are shown
in figure(5.10) and figure(5.11) respectively. From these figures it can be seen that the
Figure 5.5: The output yi(t) of the inner loop, example(5.6.1)-test 1.
" nnnnn ~
• i i —— i —
.n n" "n r ~1
JUL
tOO 150 200 2SO 300 350 <00
Tlrm(Mc)
Figure 5.6: The output outi(t), example(5.6.1)-test 1.

—r —r-
- - - -
,
- - - - — r
-1
-
- _ - - _ _ - _ - - - - J - - -
61 1 1 1 i_1 I I
IOO 1 5 0 2 0 0 2 S O 3 0 0 3 5 0 I O O
Tlrrm(S9c)
Figure 5.7: The output out (t), example(5.6.1)-test 1. 2
oscillation drives the system into saturation regions for N . L . 2.
The estimated total gain from the input m to the output y 2 2 at frequency u =
2A7rad/sec, is The oscillation index for the inner loop is
I (2.47) = |1 - 1.0097| = 0.0003 ~ 0,

osc (5.177)
therefore, the oscillation is caused by the inner loop.
E x a m p l e 5.6.2 The process considered is defined in the example (5.6.1) where the non-
linearity N . L . I is as defined above where k p = 1 and the nonlinearity for N . L . 2 is
described by the following equation.
out (t) = sign(in (t))[|in (t)| + 3].

2 2 2 (5.178)
The discontinuity offset of 3 at zero models coulomb friction and the linear gain of one
models viscous friction.

Figure 5.9: The output yi(t) of the inner loop, example(5.6.1)-test 2.

0.2 -
0 ' " ' ' ' 1

100 1 SO 200 250 300 350 -400
Figure 5.10: The output outi(t), example(5.6.1)-test 2.

The simulation results of this new process showed an oscillation frequency at co =
2.25 rad/sec and a period of T = 2.78 sec. The output of the outer loop is shown in
Figure 5.12. The output of the inner loop is shown in Figure 5.13. The output of N . L . 1
and N . L , 2 are shown respectively in figure (5.14) and figure (5.15). Figure (5.16) shows
in2(t) versus out2(t) which indicates that the inner loop is driven to the non linear region
of its non linear component. The estimated total gain from the input i n to the output 2
3 -
2.5 -
_, I , . . , , 1
IOO ISO 200 250 300 3SO AOO
Tlm0(MC)
Figure 5.12: The output y(t) of the outer loop, example(5.6.2).
y2 at frequency LO = 2A7rad/sec, is
total gain = 1.039. (5.179)
The oscillation index for the inner loop is
7 (2.25) = |1 - 1.039| = 0.0039 ~ 0,

osc (5.180)
therefore, the oscillation is caused by the inner loop.

Figure 5.13: The output y\(t) of the inner loop, example(5.6.2).
Figure 5.14: The output outi(t), example(5.6.2).

Figure 5.16: The output out (t) versus in-2(t), example(5.6.2).

2
What can be learned from this study, is that oscillations are transmitted between
interacting loops, cause excess variability, and may yield wrong performance index cal-
culation.
5.7 Conclusion
In this chapter, it was shown that oscillations transmitted between interacting loops,
cause excess variability, and may yield wrong performance index calculation. Therefore,
for cascade loops a unified oscillation index that can be used for testing linear and
nonlinear loops to help localize any loops oscillating at a given frequency is introduced.
The procedure does not require knowledge of the transfer functions; however, it requires
knowledge of the frequencies involved.
Knowledge of the limit cycle conditions can be used to tune the controller (for example
using Ziegler and Nichols method [59]) and so eliminate oscillation and guarantee better
performance.
Up to now only SISO systems were considered; however, many processes are M I M O
systems. The assessment of process control of such multivariable systems will be pre-
sented in the next chapter.

Chapter 6
Control Loop Performance Monitoring in the M I M O case
This chapter is concerned with the extension of Harris's control loop performance assess-
ment technique for SISO systems to Multi-Input Multi-Output (MIMO) systems. This
extension was considered in Harris et al.[36] and in Huang et al.[40]. The approaches
developed in these two references are quite similar as both require the knowledge of the
interactor matrix to perform the assessment. This knowledge may be considered the main
drawback for these two techniques as the interactor is not unique and requires significant
knowledge of the system for its determination.
The approach presented in this chapter is based on defining an absolute lower bound
on the achievable value for each output variance, thus avoiding the need to identify the
interactor matrix.
The organization of this chapter is as follows. In the first section, strategies of control
loop performance assessment based on multi-input multi-output (MIMO) minimum vari-
ance control (MVC) are discussed. This is followed by the definition of ultimate output
variance bounds. These bounds are used in section 3 to define a performance index asso-
ciated with each process output. The accuracy of the proposed method is compared to
the existing ones and is investigated by means of an example. In section 4, a procedure
for estimating the performance index is given, followed by a simulation example, and,
finally, some conclusions are drawn.
84
Chapter 6. Control Loop Performance Monitoring in the MIMO case 85
6.1 Assessment strategies based on M I M O minimum variance control ( M V C )
In this section M I M O minimum variance control strategies are discussed together with
methods of assessing performance based on M I M O M V C . Consider a M I M O linear
discrete-time system described by the following transfer function
y(t) = T{q- )u(t)

l
+ Niq-êit) (6.181)
where y(t), u(t) are m x l output and r x 1 input vectors, respectively. e(t) denotes the
m x 1 innovations white noise sequence. T(q~ ) and N(q~ ) are rational matrices in the
1 1
unit delay operator q" . 1
In recent years, many algorithms have been described for the adaptive control of linear
stochastic systems. For more discussion of these algorithms see, for instance, Goodwin et
al.[31] (1984). Early work on (MVC) by Astr6m[l] focused on Single-Input/Single-Output
(SISO) systems. Borisson [12] (1979) considered the Multi-Input/Multi-Output (MIMO)
M V C where a common delay was associated with all system outputs. In Goodwin et
al.[32], a different delay was associated with each output. A more complicated M I M O
time delay representation was proposed by Elliott et al. [21] (1982), who showed that
the notion of delay in the M I M O case is related to the system interactor matrix £(g)
introduced by Wolovich et al.[70] (1976) which satisfies the following two properties:
1. l i m ^ o o £(q)T(q~ )
l
= K for any non singular matrix.
2- = #(9)diag(<7 , ...,q ) where H(q) is a unimodular matrix.

dl
dn
where H(q) is a unimodular matrix of the form
1 0 0
h21 1 0
H(z) 0
— h m—l 1
m
and hij(q) is either divisible by q or is zero. B y considering the Markov parameters
representation and by explaining the meaning of property (2), Shah et al. [58] gave an
interesting interpretation of the interactor matrix as follows:
•t-i
TV ) 1
= Y.Giq
i=0
= foq* + Ziq*- + ... + U-iq 1 1
where d is the order of f (q)
lim T(q)t(q) = £ G - i + ^G -
0 d d 2 + ••• + td-iG 0 = K.
K represents, then, the linear combination of the rows of the first d Markov parameters.
That means that the interactor matrix satisfies the following algebraic equations:
CoGo = 0
+ = 0
ioGd-i + £,iGd-2 + ... + £d-iG 0

=
Solving these equations gives another way of evaluating the interactor matrix. The first
method to evaluate £(q) is the one given originally by Wolovich and Falb [70] and used
in [31].
W i t h the above interactor, which is in lower triangular form, the first output variable
is decoupled from the other outputs.
The spectral unitary interactor matrix introduced by Bittanti et al.[ll] (1994)
is an interactor matrix preserving the spectral properties of the underlying system. It
satisfies
t(q)Z{q- )
1 T
= i, (6-182)
and also satisfies property (1) above.
Regardless of the interactor matrix adopted, it can be shown that the minimum
achievable output variance, for a non minimum phase system, is a finite moving average
polynomial of the form
y(t)-y*(t)=C\q)F(q)e(t), (6.183)
where {y*(t)} is a given reference trajectory and £(q) is the time delay representation
considered. F(q) and R(q~ ) are obtained from the following equation
1
Z(q)N(q- )e(t) = F(q)e(t) + R(q- )e(t),

1
1
(6.184)
where F(q) = F q d
d
+ ... + F q and i ? ( g ) is a proper transfer function. The terms
x
_1
F(q)e(t) and R(q~ )e(t) are respectively, future and past noise sequences.
1
Note that the expression for this tracking error (equation 6.183) does not depend on
controller parameters.
As for SISO assessment, one might be tempted to use the predicted performance of
a M I M O minimum variance controller as an absolute lower bound against which current
performance can be assessed. A bound is obtained when taking the variance of both sides
of equation (6.183). However, this lower bound will depend on the type of interactor used.
In the next section a new procedure for performance assessment without knowledge
of the interactor matrix is proposed.

6.2 Ultimate output variance bound
In the last section it was shown that the minimum achievable output variance depends
on the type of time delay matrix used. This non-uniqueness, together with the fact that
the structure of an "optimal" interactor matrix is not yet known, leads to the following
technique in which an absolute bound is independently established for each output.
Proposition 6.2.1 Consider a minimum phase system. The minimum achievable out-
put variance for output y.i(t) is given by:
K)L = EWCFjEf, (6.185)

k=l
where F l
is a vector containing the coefficients of the polynomial F (q) obtained from l
the following Diophantine equation
q d L n N ^ q - 1
) = F \ q ) + R ^ q ' 1
) , (6.186)
where d l
m i n ) N \ F l
and R l
are respectively, the minimum delay in row i, the i th
row of
N , the i th
row of F and that of R.
Proof: First consider the output ordering y(t) = {yi(t), y2(t), ...,y (t)) m and consider
the lower triangular interactor matrix as a M I M O time delay representation. Dugard
et al.[19] (1984) showed that an absolute bound or an unconditional minimum output
variance is achieved for the first component of the process output, i.e. for yi(t). For
the other components, the optimization is achieved subject to constraints. To determine
the lower bound achieved by output (yi(t) — y*(t)), return to the major steps in the
derivation of the M V C strategy.
• First note that the first row of the interactor matrix is
6(9) = ..,o),
where e ^ is the minimum delay in the first row of

in T ( q ~ 1
) .
• Then, the first row of equation (6.184) is
^ L n j V 1
^ - 1
) = F ^ q ) + R ^ q ' 1
) ,
where F 1
( q ) = q F \ + ... + q d l
m i n F i , N (q),
x
and R 1
^ 1
) are respectively the first
min
row of F ( q ) , the first row of N(q~ ) l

and that of R ( q ~ 1
) .
• Under the M V control and in steady state the tracking error is a finite moving
average polynomial given by equation (6.183). Therefore,
yi{t)-y{{t) = q - d l
^ F \ q ) e ( t ) ,
A n d the minimum achievable output variance (<r^) „ for output yi(t) m is
d . 1
E{(vi(t) - yKt)) } = E 2
F}Q(F}) - T
(6-187)
Now, for an arbitrary output yi(t), consider the following output ordering
yif) = (yt(t), y2(t),...,yi(t),.,y (t)), m
then replace the subscript 1 by the subscript i in the above equations to complete the
proof. •
It is shown in Dugard et al. (1984) [19] that the unconditional minimum output
variance ( C y . ) m v is achieved for all outputs if and only if the system is minimum phase
and the interactor matrix is diagonal. The test for a diagonal interactor is established
by the following result.
Proposition 6.2.2 Define d m i n = mini<j< dkj, i.e. the minimum delay on row i of the
r
matrix T ( q " 1
) . Define the matrix D ( q ~ 1
) as
D ^ 1
) = { k i j q - d
- }
D(q~ ) is the first non-zero Markov parameter of T(q ).

1 1
Define D associated with
D(q~ ) as D = {kijSij} where

1
6^ = lifd ij = dl nin (6.188)
= 0 if not. (6.189)
T(q~ ) has a diagonal interactor matrix

1
£(q) = diag(q ™™, ...,q ™™)

d d
(6.190)
if and only if the associated matrix D has full rank.
Proof: To prove this result let us express the transfer function as:
T(q~ ) = {hjq-^h^q- )},

l 1
(6.191)
where kij and dij are constant and delay between the input j and the output i.
, / -IA _ i +M - +- + M ~ 1 m
n i j [ q
>~ l + a^-1 + ... + M ^
is a transfer function without delay between the input j and the output i. £{q) =
diag(q ,q )
ai am
is an interactor for T(q~ ) if and only if
1
limefaOTfar ) 1
= Um {kijq'^^hijiq- )} 1
=D (6.192)
D is nonsingular only if a* = min, d^. Otherwise if cn, < min^ d^ then the i th
row of D
is a zero row therefore D is singular. If instead OJ; > min^ d^ then at least one infinite
zero is present in D; therefore, it is nonsingular. •
The next section deals with performance indices.

6.3 M I M O Performance Indices
Proposition 6.2.1 defines a unique absolute bound associated with the output yiit). If
the process is minimum phase then this bound is achievable by using the lower triangle
interactor and considering the output ordering (yi(t), ...,y (t)). n Furthermore, if the pro-
cess is minimum phase with a diagonal interactor matrix, all outputs can achieve their
lower bounds simultaneously. Therefore an absolute measure of controller performance
may be built based on these bounds by defining the performance index associated with
output yiit) as:

2
W ) ) = ( 2 ^ • (6.193)
y yi{t))mv
a
Where CT .^ is the actual output variance of the output yi(t), and (o- (t))
2
yi mv is the min-
imum achievable output variance bound defined by equation (6.185). Note that this
unique bound is achievable by at least one output. For an estimation of the overall per-
formance, it is proposed the following reference bound which is unique and can only be
achieved if and only if the interactor is diagonal.
m
(° y(t))ref = Y,K(t))rnv.
2
(6.194)
i=l
The M I M O performance index is defined as
W ) ) = J^l2 ? M(
• (6-195)
It is clear that the procedure of assessment will be mainly based on equation (6.193).
In many cases, not all outputs have the same importance, and it is thus useful to know
what would be the achievable limit for an individual output.
The benefits of the proposed performance indices over existing ones will be shown by
means of the following example.

Example 6.3.1 [19] (Dugard et al. (1984)) The process considered is described by
,-1
yi(t) 0 Ui(t) 1 + Cl q~ l
0 ei(t)
+
q-1
- 2q -2 .-3
u (t)
2
0 1 e (t)
2
where £[ei(*) ei(*)] = of and E[e {t) e {t)}

T
2
T
2 = o\.
In the following, three different minimum variance control strategies, MVCi, MVC 2
and MVC$, obtained for three different interactor matrices, are considered. MVC\ and
MVC 2 are obtained respectively when considering the lower triangular interactor matrix
with the output order (yi(t), y2(t)) and (y2(t), yi(t)). MVC 3 is obtained when a unitary
interactor matrix is considered.
MVC\'. For the output ordering y(t) = (yi(t) y2(t)), the interactor is
q 0
-q + 2q
2
3
q3
Based on this interactor matrix, the output variances are
E{bn(t)-ym?} = oi,
E{[ (t)-y* (t)} }
y2 2
2
= (l + ^l)al + o 2
2
and E{[y(t) - y*(t)] [y(t) - y*(t)]} = (2 + 5cj)a + o\.

T
2
MVC :
2 If instead y(t) = (y (t) Vi(t)), then the interactor matrix is
2
q 0
6(?) =
+ 2q )
2
-(q 3
q 3
Based on this interactor matrix the output variances are
E{[yi(t)-yKt)} } 2
= (i + cl)o , 2
Emt) - y* (t)?} = 4 2
and E{[y(t) - y*(t)f[y(t) - y*(t)}} = (1 + c\)a\ + o\.

MVCz' The unitary interactor is determined as
[ 0.35? + 0.35? - 0.35q

2 3
0.35q + 0.35g
3
-0.35q - 0.35<? 3
-0.35g + 0.35g + 0.35g
2 3
The M V C based on this interactor matrix (equation 6.183) leads to
E{[yi(t)-ym} } = (l + l4)ol 2
EiM^-ymnÂ+o-i
and E{[y(t) - y*(t)] [y(t) - y*(t)}} = (1 + \c\)ol T
+ cx .
2
As can be seen, this control strategy achieves better overall performance than M V C l and
M V C 2 . However, it does not provide an absolute minimum output variance for any of
the individual outputs. Further, the unconditional minimum variance control is achieved
for (t)
Vl by M V C l and for y {t) by M V C 2 .
2
This example shows that the results of the M V C strategies differ depending on the
type of interactor matrix considered. Now, consider them as three different controllers
and compare them using the performance index as introduced in this chapter and also
that proposed in [40, 36] (Huang (1997); Harris et al. (1996)). If these controllers ( M V C l ,
M V C 2 and M V C 3 ) were implemented on the process considered, the output variances
would be those determined above. Therefore, according to the proposed approach these
output variances should be compared to the following performance reference bounds:
(°~yi )ref =
°\ f ° output yi and
r
{vyzYref =
a
2 f o r
output y. 2
The overall performance will be compared to (o~ (t))f f = o~\ + o~\. y e
The control performance of M V C l is first analyzed. Using the previous performance

index definition, the following results are obtained:
\°~y\)ref
h\W)) = 7 = 2 f O T
J/2,
\ y2)ref
a
°2
(2 + 5c )a + o\ (2 + 5c?)<r + a\
WW = ,
2 2 2
N2
W r e /
- a„2 , J
°2 l+
for an overall performance assessment.
When applying the approach described in [36, 40], the reference output variances are
those obtained when M V C 3 is implemented which are:
(°yi)ref = (! + g i > i c
for output y i ,
= ^ i°"i+ 2
c a
for output y . . 2
The overall performance is then compared to (cr (t))f f y e

=
(1 + f i ) ° " i + °2- Based on that
c
the performance indices would be
cr 2
5
^P(?/I(*)) = 7—T2~ = 1
+ QCI > 1 for
output yu
\ yi)ref
a
°
r^ w\ a
v2 (1 + 5c )a + a\ 2 2
t
IpW)) = , \ 2 = — i 2 o , o — for output 2/2,

r ( aw a
v (2 + 5c )a + <r 2 2 2
2
The same analysis is carried out for the two remaining controllers (MVC2 and M V C 3 ) .
Results are shown in Tables 6.2, where the notation used is:
• U . I.: the performance indices based on the use of unitary interactor.
• U . B . : the performance indices as defined in equation 6.193 and 6.195.
Note that in Table (6.2) the values of I are expected to be greater or equal to one, but
p
in column 2 and 4 of this table, I (yi(t)) p and I (y2(t)) show values less than one. This
p
Ip i (yi(t))
P
ipMt)) i (y(t))
P
MVd 0.66 1 5.5 6 2.7 3.5

MVC 2 1.2 1.8 0.9 1 1.07 1.4
MVC 3 1 1.5 1 1.1 1 1.3
Approach U . I. U. B. U . I. U. B. U . I. U. B.
Table 6.2: d = 0.9 o = \ a = 1 x 2
occurs because the M V C reference used is based on a unitary interactor which is not
optimal when applied to individual outputs. The overall performance of MVC3 is better
than that of MVC\ and MVC 2 and is very close to the proposed performance references
(compare column 6 with column 7). This example shows that assessment based on the
bounds defined in Section 3 is accurate. In the next section the steps for estimating the
control loop performance based on previous bounds are given.
6.4 Estimation of the lower bounds
Proposition 6.4.1 The bounds defined by equation (6.185) are feedback invariant.
They can be estimated from routine closed-loop operating data.
Proof: Express T(q~ ) as T(q~ ) = [2 ...T ...T ] and i V f a ) = [Wi.../y ..JV ] where
1 x
1 i m
-1
i m
Ti and Ni are the i th

row respectively of T and N. Define y(t) = (yi(t), ...,y (t)) m and
y(t) = (yi(t),y (t),

2 ...,y (i)).
m It is easy to see that y(t) = My(t) where the (1, i)
component of M is equal to one (i.e M\ = 1), M ^ i = 1 and ti
0 0 0 .. 1 0
0 1 0 .. 0 0
M
0 0 0 .. 0 ..
then let the feedback control law be u(t) = —Ky(t). Substituting this equation into
equation 6.181 yields
y(t) = ((/ + TKy'NYq-êit). (6.196)
Multiplying equation (6.181) by M yields
y(t) = My(t) = MTiq-ûit) + MN^q-êit). (6.197)
The new output ordering is y(t) = (yi(t), y2(t), ...,y (t)). m The new transfer function is
MT{q~ ) l
= [T ...T ...T ]
i 1 m and the new noise term is MN^ ) 1
= [N ...N ...N ].
i 1 m The
lower triangular interactor matrix £ (q) is determined for the new transfer function MT.
m
Multiplying equation (6.197) by £ (q) yields m
U{q)My{t) = U(.q)MT(q- )u(t) l

+ U(q)MN{q- )e(t)- l
(6-198)
The noise term is expressed as
U(Q)MN(a- ) 1
= F(q) + R(q- ). 1
(6.199)
Substituting equation (6.196) and equation (6.199) into equation (6.198) yields
U{q)My{t) = [R- UMTK{I + T K ) " i V ] ( g - ) ( t ) + F(q)e(t).

1 1
e (6.200)
Because the matrix H(q~ ) = (R-£ MTK(I x

m + TK)- N)(q~ )1 1
is proper, the first term
on the right hand side of equation (6.200), which depends on the controller, contains
only present and past values of e(t). The second term, which does not depend on the
controller, contains only future values of e(t). Therefore, the two terms are independent.
The first row of equation (6.200) is given by
qdl .(t)
iny = [R- UMTK(I + TKy'NUq-êit) + F (q)e(t),
1
(6.201)
where [R- £ MTK(Im + T K ) ' ^ is the first row of the matrix H = R - £ MTK(I
1
m +
TK)~ N
l
and F = F (q) is defined in equation 6.186. Equation 6.201 is equivalent to
1 l
yi{t) = (F* + ... + Fk . <THe(t) + (a di . q~<^

+1 + ....) (t). e (6.202)
mm mm
From the right hand side of equation (6.202) the first term, which is the i th
process output
under minimum variance control, is not affected by any feedback control. Therefore,
modelling the process output y(t) by a multivariate M A process
y(t) = $ a(t) + $ q- a(t-l)

1 2
1
+... = # ( g _ 1
K (6.203)
where a(t) is zero mean white noise, the i th

process output is
Vi(t) = ((ft)* + - + . )*q- ^)a(t)

d
+ . + 1 ) <--
g
1
+ ...)a(t). (6.204)
mm min
Thus, the first terms on the right hand side of equation (6.204) can be estimated as
(<j>iy = E[y (t)al }Q- ,

i j
1
where Q = E[a a[], then the lower bound can be determined as

t
«( u = (tiro-mr+...+(4.
^
t)
v
' mtn
YQMS .
mm
rr (6.205)
Based on this, the procedure of control loop performance assessment consists of the
following steps.
Algorithm for M I M O conrol loop assessment:
1. Determine the delays between process inputs and outputs.
2. Determine the process output M A or A R model to obtain the estimated sequence
a(t). (equation (6.203))
3. Determine the estimated lower output variance bounds of different outputs by use
of equation (6.205).
4. Determine the performance indices associated with different outputs (equations
(6.193) and (6.195)).

To know whether those bounds are all achievable it is necessary to determine the matrix
D and determine whether the interactor is diagonal. A l l steps in this procedure can be
performed on line with a small amount of computation. In the following, two simulation
examples are considered.
6.5 Simulation examples
6.5.1 First example
This first example is from Harris et al.[36]. The process considered is a catalytic packed-
bed reactor having the following transfer function
5.158<r -6.643<r +1.8<r

5 4 3
0.5384g-
l+0.3412 - -1.04g- 2 1
l+0.5203g- -1.262q-2 1
G{q- ) =
9
X
-3 0.1143<7- +0.1592<?-
6 5
-1.012?- - 3.398?- - 0.7498g"

5 4
1-0.6143?- 1
The disturbance process N(t) is a multivariable Integrated-Moving Average process of
order (1,1):
-l
1 + 0.126?- 1
-0.1249?- 1
1 - q'1
0
N(t) = e(t),
0.212?- 1
1 + 0.126?- 1
0 1 - q- 1
where e(t) is a Gaussian noise sequence with covariance matrix given by:
0.134 0.043
Q=
0.043 0.2
Though this process is nonminimum phase, the assessment described above can still
be used for comparison purposes [36]. The objective of this example is to compare the
bounds obtained when the unitary interactor matrix is used to those obtained when using
equation (6.185).
The unitary interactor associated with this process is given by:
1 0.355 - 0.355 q~ l
q 3
0
-0.4166 0.8521 + 0.1479 q' 1

0.4167 q 4
q 4
0.852071 q + 0.147928 q 3 4
-0.355 q + 0.355 q 3 4
-0.35497 q + 0.35507 q 3 4
0.1479 q + 0.8521 q 3 4
then F(q), i?(<? ) defined in equation 6.184 and ( £ F ) ( < 7 ) defined in equation 6.183
_1 -1 -1
are given by
1.126g + 1.126<7 + 1.093q + 0.148<? -0.125q - 0.125q - 0.157q + 0.355q

2
3
4
2 3 4
F(q) =
0.211 q + 0.211 q + 0.225 q + 0.355 q
2 3 4
0.609 q + 0.609 q + 0.62 q + 0.852 q 2 3 4
1.1265 + ^ -0.1249 0.1249

' (-1+9)
RiQ- ) =
1
0.6096
0.2113+ 0.6096 + (-1+9)
and
1 + 0.241 q~ + 1.126 g
3
- 2
+ 1.126g" 1
0.197q~ - 0.124q~ - 0.125?-
3 2
1
0.58 q~ + 0.211 q~ + 0.211 q-

3 2 1
1 + 0.475 q~ + 0.609 q~ + 0.6099-3 2
1
Therefore, the output covariance matrix is
0.551 0.127
E{[y(t)-y*(t)}[y(t)- *(t)] } y
T
0.232 0.420
Based on this strategy the output variances to take as references against which the control
loop performance will be compared are:
KUf = E{( (t) yi - yl(t)) } 2

(6.206)
= 0.551, (6.207)
« W = E{(y (t)2 - y* (t)f)

2
(6.208)
0.420 (6.209)
and
KW
2
= E{[y{t)-y*{t)) [y{t)-y\t))}
T (6.210)
= 0.551 + 0.420 (6.211)
= 0.971. (6.212)
When applying the approach introduced above to determine the output variance refer-
ences, the minimum delay in each row is three; therefore, equation 6.186 for each output
is
q'N^q) = 0.211 0.211

q + 1.126? + 1.126? -0.125? - 0.125? +
2 2
3
1-g- 1
1-g - 1
q N (q) = 0.211 0.6096

0.212? + 0.212? ? + 0.6096? + 0.6096? +
3 2 3 2
2
1-q-1
l-q- 1
Therefore the lower output variance bound for yi(t) is
K)ref = E{( (t)

Vl - y{(t)) } 2 (6.213)
= 0.4558, (6.214)
and, for y (t) 2
(<4W = E{(y (t)2 -y* (t)) }

2
2 (6.215)
= 0.3827. (6.216)
The overall output variance bound would be
KW2
= E{(y(t)-y*(t))*} (6.217)
= 0.8385. (6.218)
The results of this example are summarized in Table 6.3. It can be seen that the references
given in equation (6.207)-(6.209) are bigger than those established in equation (6.214)-
(6.216). That is because the M V C based on the unitary interactor does not achieve an
ultimate minimum variance for the individual output. Note, also that these results are
very close to those obtained using the approach introduced here where the knowledge of
the interactor matrix was not used.
Table 6.3: Reference values for example 2.

references for yi(t) references for y (t) references for y(t) 2
0.551 0.4558 0.420 0.3827 0.971 0.8385

Unitary Int. U . bound Unitary Int. U . bound Unitary Int. U . bound
6.5.2 Second example
For this example the process considered is described by
-0.2(?-
i -0.12( -
—yj.J-t.il 0.06(7- 3
1 2
+0.12g- +0.2g- 2
?
3
l-0.86(j- -|-0.12q- +0.2g-
1 2 3 ui(t)
l-0.86g- +0.12 - +0.2g-

1
g
2 3
l-0.86(j- -f0.12q- +0.2(j-
1 2 3
J u {t)
2
1-0.369- 1
l-0.86q- +0.12g- +0.2(j-

1 2 3 u ei(t)
+ 1-OAq- 1
0 l-0.86g- +0.12 - +0.2 -

1
g
2
g
3
where
ei(t) 0.42 0
E ei(t) e (t)
2 )\=Q =
0 0.38
By choosing the output ordering (yi(t), y (t)), the interactor matrix is

2
q 0
-2.5g - 0.5? 3 2
-
<t
3
The M V C based on this interactor leads to the output tracking error
yi(t)-yl(t) 1 0
, (6.219)
y2(t)- my -1.250- + 1.02g"
1 2
1 + 0.469- + 0.27g- 1 2
C2(*)
and the output variances are
E{[yi(t)-yl(t)} } 2
= 0.42 and
E{\jfy(t) - y»(t)?} = 1-59-
To apply the proposed control loop performance assessment procedure, it is needed to
know the minimum delay for each output and a M A or A R model for process outputs.
In this case the minimum delay for both outputs is one, and equation (6.219) is the
M A model needed. Therefore, according to equation (6.185), the minimum achievable
output variances for outputs y\ and y are 2
d . 1
mtn
E{(m(t) - yl(t)) } = 2
J2 jQ( j) F F T
= 0-4229 and
j=i
d? .
E{(y (t)
2 - y* (t)) }
2
2
= £" F Q(F ) 2 2 T
= 0.3802.
j=i
Therefore, the performance index for y\ is I (yi(t)) p = 1 and that of y is I (y (t))2 p 2 = 4.19.
If instead y(t) = (y (t)

2 yi(t)) then, the interactor matrix is
o
=
-0.40? + 0.080g 3 2
The M V C based on this interactor leads to the output tracking error
1 + O.bq' + O.Slq- 1 2
-0.18Q- + 0.07q-
1 2
vi(t)-yi(t)
(6.220)
y2(t)-y* (t)
2 0 1 e (*)
2
and the output variances are
E{[yi(t) - y*M } 2
= 0.5838 and
E{[y2(t)-y* (t)?} 2 = 0.3802.
In this case the minimum delay for both outputs is still one and equation (6.220) is the
M A model to be used in the same way as before to obtain that the minimum achievable
output variances for outputs yi and y are 2
d . 1
min
E{( (t) yi - yl(t)) } 2

= £ F}Q(F}) T
= 0.4229 and
d . 2
E{(y2(t) - y* (t)) } 2
2
= 2:^(^ = 0.3802.
Therefore, the performance index for y\ is I (yi(t)) p = 1.39 and that of y is I (y (t)) 2 p 2 = 1.
The process described above is now simulated with the following controller.
" k 0.04a~ 2
-0,0344q~ 3
+0.0048Q~ 4
4-0.0080a~ 5
fc - 0 . 0 6 O ~ 2
+ 0 . 0 5 1 6 Q ~ 3
- 0 0 0 7 2 q - 4
- 0 . 0 1 2 q ~ 5
, 1
(0.02+0.019<l-2)(l_q-l) 1
(0.02+0.019q-2 ) ( 1 - q - 1 )
C(q 1
) =
. 0.5 — O . O a q - 1
-0.284q~ +0.148q~
2
3
+O.Q8q~ 4
. —0 . 2 + 0 . 0 5 2 q ~ 1
+0.0792q~ 2
— 0 . 0 5 4 4 q - 3
- 0 . 0 2 4 q ~ 4
2
(0.02+0.019q-2)(l-q-l) 2
(0.02+0.019q-^ ) ( l - q - l ) -
The design of this controller was done in two stages. The first stage dealt with decoupling;
then a proportional plus integral controller was designed for each loop.
The simulation results are obtained when applying the procedure of assessment given
above. Delays are considered known. The A R filter is determined using "arx" command
of Matlab. Once the parameters of the filter are known, its recurrent form is used to
obtain the residuals on line then updating the minimum output variances and output
variances are updated considering a forgetting factor in the range of [0.9 0.99]. These
simulation results are shown in Figure 6.17.
It can be seen from combining all the results obtained above with different control
strategies applied to the considered process that at least one output can achieve a per-
formance index of one (first two controllers). The performance index for all outputs can
be close to one for a good controller such as the above last one. The proposed procedure
is able to track on-line the control loop performance.
-i i 1 i r -
~I 1 ! 1
_l I 1_ I I I
200 300 4oo eoo eoo
Figure 6.17: On line performance index estimation kl=0.08 and k2=0.05.
6.6 Conclusions
It was shown that the proposed multivariable procedure for control loop performance as-
sessment is an improvement over existing techniques. No prior knowledge of the interactor
matrix is needed, and only knowledge of the minimum delay for each output is needed.
The method is then a natural extension of SISO control loop performance assessment
to M I M O control loop performance assessment. This procedure can be used to identify
the loops showing a large output variance and so in need of special attention, from those
having an output variance close to the M V C benchmark. A significant advantage of the
bounds derived here is that they are independent of the interactor matrix.
In the following chapter the proposed method is evaluated by application to the loops
controlling a lime kiln in a Kraft pulp mill.

Chapter 7
Performance Assessment of Control Loops on a Lime Kiln
In this chapter, the method proposed in the last chapter is evaluated by application to the
loops controlling a lime kiln in a Kraft pulp mill. The data was obtained in collaboration
with engineering staff at Canfor, Prince George and P A P R I C A N , Vancouver Laboratory.
The organization of this chapter is as follows. In the first section the lime cycle is
described. This is followed by the development of a dynamic model of the lime kiln. The
results of the assessment procedure are then given in section 4, followed by a conclusion
in the last section.
7.1 Process description
The lime cycle is an important part of the kraft recovery and recausticizing system in a
pulp mill. The purpose of the cycle is to regenerate the alkaline white liquor from the
spent green liquor recovered from the pulping process, see Figure (7.18) from Dumont
[20]. The white liquor, which is a solution of sodium hydroxide (NaOH) and sodium
sulfate (Na2S), is used to break down the lignin in the wood. The green liquor is a
solution of sodium carbonate (Na2COz), sodium hydroxide (NaOH), and sodium sulfite
(Na S).
2
The first step in the lime cycle is to add lime (CaO) to the green liquor where it reacts
with water in solution to produce calcium hydroxide. This reaction is called slaking and
106
Chapter 7. Performance Assessment of Control Loops on a Lime Kiln 107
SLMUU WHITE
WHITE LIQUOR
C A U 5 T I C I Z CR
L14 U 4 ft CL A RIME It LIQUOR
fUEL LIHE
HUD WASHER
HUD
Figure 7.18: Lime cycle.
can be written as:
CaO + H 0 2 — • Ca(OH) 2 + heat.
The calcium hydroxide then reacts with the sodium carbonate in the green liquor to
produce sodium hydroxide and calcium carbonate. This reaction is called causticizing
and can be written as
Na C0
2 3 + Ca(OH) 2 — • 2NaOH + CaC0 . 3
The white liquor (NaOH) is separated from the calcium carbonate (lime mud) by passing
the solution through a clarifier, where the calcium carbonate settles out as an insoluble
mud. To complete the cycle, the calcium carbonate or lime mud is washed and filtered,
then dried and heated in large rotary kilns. The heat causes the calcium carbonate to
break down, giving off carbon dioxide gas and calcium oxide or lime. This reaction can
be written as
CaC0 3 + heat —• CaO + C0 2 T.

The lime produced is then used in slaking. From the above, it is evident that the lime
kiln is a crucial element of the lime cycle.
7.2 Lime kiln process
FUME S ID FAN
LIME MUD
HO T END
Figure 7.19: Lime Kiln.
The lime kiln, depicted in Figure (7.19), is a long rotating cylinder, sloping typically
at four centimetres per meter, and rotating at one revolution per minute. The lower end
of the kiln is referred to as the hot end or the feed end. The fuel and air for combustion
enter at this end. The lime mud enters into the kiln at the feed end or the cold end. It
is then dried and converted into hot lime as it flows down the kiln toward the hot end,
which is maintained at a temperature of about 1000°C. Detailed descriptions of the lime
kiln can be found in [61].
Clearly, lime kilns are multivariable distributed systems. The objective of a lime
kiln control system is to produce lime of good reactivity while minimizing the energy
consumption. Underburned or overburned lime from an ineffectively controlled kiln can
be very costly in terms of reduced pulping efficiency and increased energy costs. A review
of the principles and current practice of lime kiln control is given in Dumont [20].
The manipulated variables may include mud flow rate, fuel rate, rotational speed of
the kiln, primary air flow (FD Fan), and draft (ID Fan). However, setting of the mud
flow rate depends on the availability of stored mud and the desired lime production rate
and so may not be adjustable for control purposes. Typically, the controlled variables
include the feed end and the back end temperatures. The quality of the lime produced
depends upon the temperature profile in the longitudinal direction of the kiln. Another
control concern is the oxygen content of the fume gasses, which indicates how efficiently
fuel has been burned.
For the K i l n under consideration, the controlled variables (outputs) are: feed end
temperature (FET), hot end temperature (HET), and the excess of oxygen (XS 02)-The
rotational speed of the kiln is adjusted according to the load. The manipulated variables
(inputs) are the fuel rate and air flow (ID fan). The lime kiln is considered here as a
multivariable system with two inputs and three outputs. The models are assumed to be
first-order with time delay determined from the two bump tests shown i n figure (7.20)
and figure (7.21). The resulting transfer function model is shown below:
0.52e- 67s
-0.002e- 12s
HET 87s+l lOs+1

FUEL
= 0.062e- 0.228e- - (7.221)
13s 5 7s
FET 10.52s+l 23s+l

FAN
-0.0088e-°- 67s
0.018e-° 1 7 s
Xs0 2
2.5s+l 1.75s+l
where all time constants and delays are expressed in minutes. Dynamic operating re-
sponse data were collected and compared to the response of the model determined above.
The model delays show that the excess of oxygen can be adjusted with less than a minute
of delay. However, a delay of at least 12 minutes from input to any change in the hot end
temperature output and 5.7 minutes for the feed end temperature. In Figure (7.22), the
Q.
< 10 20 30 40
Tim e ( m n )
Figure 7.20: First bump test.

d
LL 150 200
1000 —1
E
Q.
980
z
< 960
100
Tim e (mn)
Figure 7.21: Second bump test.

true response curves from mill data in solid line, while those obtained from the model are
in dashed line. It can be seen that the model provides a reasonably good approximation
of the kiln's response curves.
In the next section the control performance of this process is investigated.
1 100
III \ / ~*
500 1000 1500 2000 2500 3000 3500 4000
1
380 - i i i i
O>360
3-340
£±320
300 i i i i I 1
500 1000 1500 2000 2500 3000 3500 4000

Tim e ( m n )
Figure 7.22: Validation of the model.
7.3 Loop performance assessment for the lime kiln
The control loops of this process are assessed using the approach described in the last
chapter. It consists of the following steps:
1. Determine the delays between process inputs and outputs (from the model above).
2. Determine the multivariate M A time series model for the process outputs
(y - y)(t) = $ ( t ) + ^qâit - 1) + ... =

lfl H(q- )a ,
1
t
where a(t) is a zero mean white noise disturbance, q is the unit delay operator and
y is the mean of y. The i th

process output is
(lU ~ Vi)(t) = ((<!>{)* + ... + (tit . Tq- ^)a(t) + ((cb .

d
di + 1 ) q - « - - » + ...)a(fi).222)
m i n min '
3. Determine the estimated lower output variance bounds of different outputs by use
of the following equation
("J-*)™ = WG[($)T + +
3 1 W l
m mm
. -i)*QW i. -,Y) ,
mm
d
T
(7.223)
where Q = E[a(t)a(tY] the covariance matrix of the filtered data and ((b))* can be
estimated by
((b)Y = E[(y -y ) al }Q-\i i t j
4. Determine the performance indices associated with different outputs using
W ) ) = jH r , (7-224)
where 1 < I with higher value indicating poor control.

p
These steps can easily be performed on-line. Neither the nature of the control strategy
nor the process model is needed, but the delays between different input output pairs of
the process are required. In this case these delays, expressed in minutes, are obtained
from the model above. The delays required by the assessment procedure are the smallest
ones from each row of the above model. The on-line output variance is computed using
the following equation
K)N = A(tT )jv_i

2
+ (i - x)( (N)
Vl - m)?,
where A e [0.98 1] is the forgetting factor. The resulting on-line performance index for
each output is shown in figure (7.23). This preliminary study shows that the performance
41 1 1 1 1 r
I i i i i i
1500 2000 2500 3000 3500 4000
30 I 1 1 1 1 1—
I 1 I I I I I
1 5 0 0 2 0 0 0 2 5 0 0 3 0 0 0 3 5 0 0 4 0 0 0
Tim e (m n)
Figure 7.23: On-line performance indices.
index of the hot end temperature is in the range of [1.8 3.2] , that of the feed end
temperature is within [4 8] and that of excess of 0 2 is within [10 20].
The conclusion to draw at this point is that the hot end and the feed end temper-
atures are relatively well controlled while the control of the excess of oxygen may not
be. Note however that there will be a contribution to the performance index arising
from the frequent setpoint change. Based on this study it was decided to design a new
controller based on MPC (model predictive controller) to replace the one in place which
is a combination of PI and Smith predictor.

7.4 Conclusion
The application of control loop performance monitoring presented here shows the impor-
tant role that the monitoring tools can play in maintaining high process performance.
By comparing the minimum variance level that each output can achieve with the actual
variance, a decision can be made to accept, retune, or replace the control algorithm. The
assessment was performed using routine operating data and little prior knowledge of the
process. The approach proposed in [25] follows the spirit of control loop performance
assessment initiated by Harris [35] for single-input single-output systems.

Chapter 8
Summary
8.1 Contributions of this thesis
The purpose of this thesis is to develop tools and techniques for monitoring and diagnosing
the performance of control loops based primarily on operating data. The contributions
of this thesis to this topic are:
• In chapter 3, the development of a Least-Squares method for parameter identifi-
cation and time delay estimation for processes where the acting disturbances are
white noise, is presented in an attempt to overcome the difficulties of the existing
method for time delay estimation.
• In chapter 4, Harris's performance index was broadened as a measure of output vari-
ance to consider other criteria, so obtaining a better overall assessment of control
loop performance. It was shown that oscillations result in poor performance and
proposed an approach to controller tuning based on a combination of performance
criteria and output variance.
• In chapter 5, a unified oscillation index is proposed that can be used for linear
and nonlinear cascade loops to help localize the loops causing oscillation at a given
frequency. The procedure does not require knowledge of the transfer functions
involved. Analysis of responses is used to determine the set of frequencies at which
further invistigation is needed.
116
Chapter 8. Summary 117
• In chapter 6, a multivariable procedure for control loop performance assessment is
proposed. This method leads to good results. No prior knowledge of the interactor
matrix is needed and only knowledge of the minimum delay for each output is
required.
8.2 Recommendation for future work
The issues addressed in this thesis can be viewed as background for several interesting
research directions. In chapter 3, the closed-loop method developed needs further study
specially when the acting disturbances are described by colored noise. In chapter 4,
by considering the presence of a nonlinear element in a given loop, it was possible to
narrow the useful range for the performance index and to show that the performance
index under normal operation is finite. The need to develop a framework for nonlinear
systems is still an open problem. In Chapter 5, an oscillation index was introduced for
SISO systems (cascade loops only); therefore, an extension of that study to multivariable
systems is needed. In chapter 6, the M I M O case was addressed, but there remains a need
to investigate the weight assigned to each performance index for practical use.
The topic of control loop performance monitoring is relatively new and is an active
area of research. Some other areas for future work would be to study the effect of the
sampling rate on the estimation of the minimum achievable output variance, and to
develop a confidence interval for performance index in both cases SISO and M I M O . This
is needed for on-line performance assessment to compensate for varying time delay or
parameter uncertainties and therefore, for robustness.
Another area of special interest to paper industry is the machine cross directional
control and it would interesting to extend the M I M O performance index defined in this
thesis to cover such processes as well.

Bibliography
Astrom K . J . "Computer Control of a Paper Machine an Application of Linear

Stochastic Control Theory". IBM J. Res. Dev., 11:389-405, 1967.
Astrom K . J . "Assessment of Achievable Performance of Simple Feedback Loops".

Int. J. Adapt. Cont, 5:3-19, 1991.
Astrom K . J . and B . Wittenmark. "Computer Control System, Theory and Design".

Prentice-Hall, Inc, 1986.
Astrom K . J . and B . Wittenmark. "Adaptive Control". Addison-Wesley, 1995.
Astrom K . J . and T. Soderstrom. "Uniqueness of the Maximum Likelihood Esti-

mates of the Parameters of an A R M A Model". IEEE Trans. Automat. Contr., vol.
AC-19, No. 6:769-773, 1974.
Astrom K . J., C. C. Hang, P. Persson and W . K . Ho. "Towards Intelligent P I D

Control". Automatica, 28:1-9, 1992.
Basseville M . "Detection changes in signals and systems - a survey ". Automatica,

24(3):309-326, 1988.
Basseville M . and I. Nikiforov. "Detection of abrupt changes: Theory and Applica-

tions". Prentice-Hall, Inc, 1993.
Bialkowski W . L . "Process Variability, Control, Standards and Competitive Pro-

gram". Proc. ISA/92 Canada, Toronto, Canada, pages 395-418, April 1992.
Billings S. A., J.O. Gray and D.H. Owens . "Nonlinear system design". I E E Control
Engineering Series 25, 1984.
Bittanti S., P. Colaneri and M . F. Mongiovi. "The Spectral Interactor matrix for
the Singular Ricati Equation". Proc. 33rd CDC, pages 2165-2169, 1994.
Borisson U . "Self-Tuning Regulators for a Class of Multivariable Systems". Auto-

matica, 15:209-215, 1979.
Box G . E . P. and G . W . Jenkins. "Time Series Analysis forcasting and Control".

Holden-Day, 1976.
Carter G . C , A . H . Nuttall, and P. G . Cable. "The Smoothed Coherence Trans-

form" . Proc. of the IEEE, Oct. 1973.
118
Bibliography 119
Clarke D. W . "Generalized-least-squares estimation of the parameters of a dynamic

model". In Identification in Automatic Control Systems. IFAC Symposium, Prague,
1967.
Clarke D. W . "Sensor, Actuator, and Loop Validation". IEEE Control Systems

Magazine, pages 39-45, Aug 1995.
Cook P. A . "Nonlinear dynamical systems". Prentice-Hall, Inc, 1986.
Desborough L . and T. Harris. "Performance Assessment Measures for Univariate

Feedforward/Feedback Control". Can. J. Chem. Eng., 71:605-616, 1993.
Dugard L . , G . C . Goodwin and X . Xianya. "The Role of Interactor Matrix in Mul-

tivariable Stochastic Adaptive Control". Automatica, 20:701-709, 1984.
Dumont G. A . " A Review of the Principles and Current Practice of Lime K i l n

Control". Course notes.
Elliott H . and W . A . Wolovich. " A Parameter Adaptive Control Structure for Linear
Multivariable Systems". IEEE Trans. Aut. Cont., AC-27:340-352, 1982.
Elnaggar A . , G . A . Dumont and A . L . Elshafei. "Adaptive Control with Direct

Delay Estimation". Control Systems '92: Dream vs. Reality, Whistler, B.C., pages
13-16, Sept. 1992.
Ender D. B . "Implementation of Deadband Reset Scheduling for the Elimination

of Stick-Slip Cycling in Control Valves". In Tappi Process Control, Electrical, &
Information Conference, pages 83-88, 1997.
Ettaleb L . , G. A . Dumont and , M . S. Davies. " A n Extended off-line Least-Squares

Method for Parameter Identification and Time Delay Estimation". Procedings of
the 37th IEEE Conference on Decision and Control, Dec. 16-18, 1998.
Ettaleb L., G. A . Dumont and M . S. Davies. "Performance Assessment of Single and

Multivariable Feedback Controllers: a Unified Approach". submitted to Automatica.
Ettaleb L . , M . S. Davies, G. A . Dumont and E . E . Kwok. "Monitoring Oscillations

in a Multiloop Systems ". Procedings of the 1996 IEEE International Conference
on Control Applications, pages 859-863, 1996.
Ettaleb L., M . S. Davies, G. A . Dumont and E . E . Kwok. "Realistic Performance As-

sessment in the SISO Case". Procedings of Mediterranean Conference on Electronics
and Automatic Control, pages 114-117, Sep. 17-19 1998.
Forssell U . and L . Ljung. "Closed loop identification revisited". Automatica,

35:1215-1241, 1999.
Bibliography 120
Fu Y . and G. A . Dumont. "Optimum Laguerre Time scale and its On line Estima-
tion". Int. J. adpt. Cont. Sig. Proc, 5:313-333, 1991.
Fu Y . and G. A . Dumont. "On-line Evaluation of Control Loop Performance".

Procedings of 4th IEEE cca, pages 655-656, 1995.
Goodwin G. C. and K . S. Sin. "Adaptive Filtering Prediction and Control". Prentice-

Hall, Inc, 1984.
Goodwin G. C , P. J. Ramadge and P. E. Caines. "Discrete time multivariable

adaptive control". IEEE Trans. Aut. Cont, vol. AC-25:449-456, 1980.
Gunnarsson S., and B . Wahlberg. "Some Asymptotic results in Recursive Identifi-

cation Using Laguerre Models.". Int. J. adpt. Cont. Sig. Proc, 5:313-333, 1991.
Gustavsson I., L . Ljung and T. Soderstrom. "Survey Paper: Identification of Pro-

cesses in Closed Loop -Identifiability and Accuracy Aspects". Automatica, 13:59-75,
1977.
Harris T. J. "Assessment of Control Loop Performance". Can. J. Chem. Eng.,

67:856-861, 1989.
Harris T. J., F . Boudreau and J. F. MacGregor. "Performance Assessment of Mul-

tivariable Feedback Controllers". Automatica, 32(11):1505-1518, 1996.
Harris T. J., J.F Macgregor and J. D. Wright. " A n Overview of Discrete Stochastic
Controllers: Generalized PID Algorithms with Dead-Time Compensation". Can. J.
Chem. Eng., 60:425-432, 1981.
Horch A . and A . J. Isaksson. " A methode for detecting of stiction in control valves".
In IFAC Workshop on On-line Fault Detection and Supervision in the Chemical
Process Industries, Lyon France, 1998.
Horch A . and A . J. Isaksson. "A modified index for control performance assessment".
In American Control Conference, Philadelphia, PA, pages 3430-3434, 1998.
Huang B . "Multivariate statistical methods for control loop performance assess-

ment". Ph.D thesis, Dep. of Chemical Eng., University of Alberta, 1997.
Huang B . and S. L . Shah. "Practical Issues in Multivariable Feedback Control Per-

formance Assessment". In IFAC ADCHEM, Banff, Canada, pages 429-434, 1997.
Huang B . , S. L . Shah and E . E. Kwok. "Good, Bad or Optimal? Performance

Assessment of M I M O processes". Automatica, 33(11):1175-1183, 1997.
Huang B . , S. L . Shah and H . Fujii. "Identification of the Time delay /Interactor

matrix for M I M O systems using Closed Loop data". IFAC World Congress, 1996.
Bibliography 121
[44] Isermann R. "Results on the Simplification of Dynamic Process Models". Int. J.

Cont, 19(1):149-159, 1974.
Isermann R. "Digital Control System". Springer-Verlag, 1981.
Isermann R. "Fault diagnosis of machines via parameter estimation and knowledge

processing - tutorial paper". Automatica, 29(4):815-835, 1993.
Kanjilal P. P. and S. Palit. " On Multiple Pattern Extraction Using Singular Value
Decomposition". IEEE Trans, on Automatic Cont, 43:1536-1540, 1995.
Knapp H . C. and G. C. Carter. "The Generalized Correlation Method for Estimation

of Time Delay". IEEE Trans. Acoust. Speech Sig. Process, 27:320-327, 1976.
Kozuband D. J. and C . E . Garcia. "Monitoring and diagnosis of automated con-

trollers in the chemical process industries ". In AIChE Meeting. St. Louise, MO,
1993.
Landau I. D. and A . Karimi. " A n Output Error Recursive Algorithm for Unbaised
Identification in Closed Loop". Automatica, 33(5):933-938, 1997.
Ljung L . . "System Identification". Prentice-Hall, 1998.
Lynch C. and G. A . Dumont. "Control Loop Performance Monitoring". IEEE

Trans. Cont. Sys. Tech., 4,2:185-192, 1992.
Mah R. S. H . . "Design and Analyse of Process Performance Monitoring Systems".

In Chemical process control, Proceedings of the Engineering Fondation Conference,
January 18-23, 1981.
Matsubara M . "On the Equivalent Dead Time". IEEE Trans, on Automatic Cont.,
pages 464-466, oct. 1965.
Nomikos P. and J. F . MacGregor. "Multivariate S P C charts for monitoring batch

processes". Technometrics, 37(l):41-59, 1995.
Perrier M . and A . Roche. "Towards Mill-Wide Evaluation of Control Loop Perfor-

mance". Control Systems, Whistler, B.C., pages 51-60, Sept. 1992.
Pryor C. "Autocovariance and Power Spectrum Analysis- Derive New Information

from Process Data". Control Eng. 2, pages 103-106, Oct. 1982.
Shah S., C. Mohtadi and D. W . Clarke. "Multivariable Adaptive Control without a

prior knowledge of the Delay Matrix". Syst. Cont. Let, pages 295-306, 1987.
Shinskey F. G . "Process Control System". McGraw-Hill, New York, 1988.

Bibliography 122
[60] Smith K . W . " A Methodology for Applying Time Series Analysis Techniques". In
Proc. ISA/92 Canada, Toronto, Canada, pages 419-435, 1992.
Smook G. A . "Handbook for Pulp & Paper Technologists". T A P P I , C P P A , 1986.
Soderstrom T. and P. Stoica. "System Identification". Prentice-Hall, 1989.
Staton B . "Characteristics and Control of Process Disturbances". Hydrocarbon

Process, 63:51-60, 1984.
Strejc V . "Least Squares Parameter Estimation". Automatica, 16:535-550, 1980.
Taha O., G . A . Dumont and M . S. Davies. "Detection and diagnosis of oscillation

in cotrol loops". In 35th IEEE Conference on Decision and Control Kobe, Japan,
pages 2432-2437, 1996.
Tyler M . and M . Morari. "Performance assessment for unstable and nonminimum

phase systems". In IFAC Workshop on On-line Fault Detection and Supervision in
the Chemical Process Industries, Newcastle upon Tyne, England, 1995.
Van den Hof P. M . and R. J. P. Schrama. "An Indirect Method for Transfer Function
Estimation from Closed Loop Data". Automatica, 29(6):1523-1527, 1993.
Weldon A. D. "Control valves- The key component in managing process variability".

Tappi journal, 80(8):71-74, 1997.
Wise B . M . and N . B . Gallagher. "The process chemometrics approach to process

monitoring and fault detection". Journal, of Process Control, 6(6):329-348, 1996.
Wolovich W . A . and P. L . Falb. "Invariants and Canonical Forms Under Dynamic

Compensation". SIAMJ Cont. and Opt., 14:996-1008, Nov. 1976.
Zheng W . X . and C. B . Feng. "Identification of Stochastic Time Lag Systems in the

Presence of Colored Noise". Automatica, 26:769-779, 1990.

Ubc 1999-463419

Uploaded by

Copyright:

Available Formats

Ubc 1999-463419

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Ubc 1999-463419

Uploaded by

Copyright:

Available Formats

CONTROL LOOP P E R F O R M A N C E ASSESSMENT A N D OSCILLATION

A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF

THE REQUIREMENTS FOR T H E DEGREE OF

DEPARTMENT OF ELECTRICAL A N DCOMPUTER ENGINEERING

We accept this thesis as conforming

Department of ^lg-^VW <~»X gu^«k Cs^^ovW— wxWi^

The University of British Columbia

List of Figures vii

1.1 Control Loop Performance Monitoring 1

1.2 Outline of this thesis 8

2 Control Loop Performance Monitoring: SISO case 11

2.1 Minimum Variance Control (MVC) 11

2.2 Evaluation of performance 15

2.2.1 Evaluation of I using time series

2.2.2 Evaluation of I using Laguerre network

2.2.3 Time delay estimation 19

3 A n Extended off-line Least-Squares Method and Time Delay Estimation 22

3.2 Least squares parameter estimation 24

3.2.1 Least-squares method 24

3.2.2 Generalized least-squares 26

3.2.3 Extended least-squares 29

3.3 An ARMA model for a closed-loop system 30

4 Realistic Performance Assessment in the SISO Case 47

5 Monitoring Oscillations in a Multiloop System 61

6 Control Loop Performance Monitoring in the M I M O case 84

6.1 Assessment strategies based on MIMO minimum variance control (MVC) . . . . 85

6.5.2 Second example 101

6.6 Conclusions 104

7 Performance Assessment of Control Loops on a Lime K i l n 106

7.1 Process description 106

7.2 Lime kiln process 108

7.3 Loop performance assessment for the lime kiln 112

7.4 Conclusion 115

8.1 Contributions of this thesis 116

8.2 Recommendation for future work 117

6.2 C i = 0.9 <TI = 1 CT = 1

2.1 Process control 12

2.2 Process control 13

3.3 Covariance of the residual 33

3.4 White noise and its estimate 34

3.5 Estimation value of the coefficient a 41

3.6 Estimation value of the coefficient b 42

3.7 Estimation value of the time delay d 42

3.8 Loss function V 43

3.9 Estimation value of the coefficient a by ordinary least-squares 44

3.10 Estimation value of the coefficient b by ordinary least-squares 44

3.11 Estimation value of the delay d by ordinary least-squares 45

4.12 Process control 49

4.13 I and steady state error versus P-gain

4.14 A single nonlinear loop 51

4.15 Saturation characteristic 51

4.16 Nyquist Diagram 52

4.17 Partition of Ip curve 52

4.18 Performance analysis 54

4.19 Process approximation 55

4.20 Process control 56

4.21 Process control 57

4.22 Comparison of true and approximate model 58

4.24 Estimated curves versus Controller parameters 60