A STUDY ON ADVANCED CONTROL FOR AN
INDUSTRIAL SCALE DISTILLATION COLUMN:
MODEL DEVELOPMENT AND CONTROL SIMULATIONS

By

John Martin Wassick

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering
and Systems Science

1987

ABSTRACT

A STUDY ON ADVANCED CONTROL FOR AN
INDUSTRIAL SCALE DISTILLATION COLUMN:
MODEL DEVELOPMENT AND CONTROL SIMULATIONS

By

John Martin Wassick

This thesis describes the research that was performed to investigate
the modelling and control of an industrial scale distillation column
located at The Dow Chemical Company's Michigan Division. The overall
objective of this thesis was to determine if advanced control would be
beneficial to the column under study. This objective led to two major
themes in the research: 1) develop a dynamic, multivariable model of

the column, and 2) propose an alternate model based control scheme and

test it through simulation.

The model that is developed is a 2 input, 2 output matrix of
discrete time transfer functions. A novel identification procedure is
developed which greatly reduces the number of parameters necessary to
describe the behavior of multivariable systems with multiple time delays
and is based on a new ”delayed polynomial matrix" representation of
discrete systems. A computer algorithm for the delayed polynomial
matrix method is described in a way that makes it suitable for
interactive use. A simulated example demonstrates the ability of this

new' identification method to perform with up to 20% additive noise on

the data.

Least squares parameter estimation of the identified model is based
on real operating data. It was found that the concentration on the 57th
tray exhibited inverse response to changing reflux flow. This was
totally unexpected and has not been mentioned in the chemical
engineering literature. Discussion of the Linde column model shows that
other aspects of the model were very consistent with the experience of

the operating personnel.

The alternate control strategy selected is a feedfoward version of
Internal Model Control, IMC. Two peculiarities of the column required
extensions of existing theory on IMC, they are: l) multirate sampling of
the two product concentrations and 2) the inverse response already
mentioned. The multirate sampling problem is addressed by a unique
implementation of feedforward control. A reduced order controller

design technique is developed to handle the non-minimum phase behavior.

The multivariable controller out performs the conventional
controllers in load change simulations. Improved disturbance rejection

was achieved in both product streams.

To Maria, my loving wife, whose patience and sacrifice made this

accomplishment possible.

ii

ACKNOWLEDGMENTS

I am indebted to several individuals for help in preparation of this
thesis. I thank Lal Tummala, my thesis advisor, for his guidance and
his thoughtful comments on the manuscript. I am also thank the rest of
my guidance committee: Robert Barr, Robert Schlueter, Hassan Khalil,
Krishnamurthy Jayaraman, and David Yen for their timely suggestions
which saved me many hours of frustration. My heartfelt appreciation
goes to my wife, Maria, for the time she found in her hectic day to help

type the manuscript.

I must express my sincere gratitude to The Dow Chemical Company for
its generous financial support and for the use of its process equipment
and computer resources. In particular, I thank Ray Murphy and Fred
Swinehart for having the enough faith in me to muster the support for
such a long project. Finally, I thank my colleague Dave Camp for his

encouragement that ensured I finished this project.

iii

TABLE OF CONTENTS

LIST OF TABLES vii
LIST OF FIGURES viii
Chapter One, Introduction and Description of Thesis 1
1.1 Introduction 1

1.2 Research Objectives 2

1.3 Outline of the Thesis 4

1.4 Summary 5
Chapter Two, Description of the Process 6
2.1 Fundamentals of Distillation 6

2.2 Distillation Control Fundamentals 9

2.3 Description of the Linde Column 9

2.4 Current Control of the Linde Column 12

2.5 Summary of Chapter Two 12
Chapter Three, Linde Column Modelling 14
3.1 Literature Survey of Distillation Modelling 14

3.2 Current Model Identification Methods 16
3.2.1 Identification Literature 18

3.3 Delayed Polynomial Matrices 19

iv

3.4

3.5

3.6

3.7

Structural Identification
3.4.1 Phase 1
3.4.2 Phase 2

Computer Implementation

Singular Value Decomposition
Computer Algorithm

Noise Compensation

SVD Subroutines

Simulation Results

Summary of Computer Implementation

wwwwww
MMMU‘UIU‘
O‘U‘J-‘WNH

Column Model

Linde Column Operating Data

Least Squares Parameter Estimation
Transfer Function Identification
Results of Column Modelling
Discussion of the Linde Column Model

H‘

O‘C‘O‘O‘b‘fl
. . . . Q,

09030009001"
mwat-‘m

Summary of Chapter Three

Chapter Four, Advanced Control for the Linde Column

4.1
4.2

4.3

4.4

4.5

4.6

Current Control Strategy

Literature Survey on Distillation Control
Model Based Control, 8180

4.3.1 Minimum Phase Processes

4.3.2 Non-minimum Phase Processes

Model Based Control, MIMO

4.4.1 Feedforward Internal Model Control
4.4.2 IMC for the Linde Column

Results and Discussion

Summary of Chapter Four

Chapter Five, Summary and Recommendations

5.1

5.2

Summary of Linde Column Modelling

5.1.1 Results of Linde Column Modelling

5.1.2 Future Research on Distillation
Modelling

Summary of Advanced Control for the Linde

Column

5.2.1 Results of Advanced Control
Simulations

95

95

98
100
102
106
118
118
124
131

144

146

147
148

149

150

153

5.2.2 Future Control Research 154

5.3 Future Research Opportunities 156
APPENDIX 161
LIST OF REFERENCES 163

vi

LIST OF TABLES

Table 4-1 IMC Transfer Functions 135

vii

Figure

Figure

Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure

Figure
Figure

Figure

Figure

Figure

Figure

2-1

2-2

3-1

3-2

3-3

3-4

3-5

3-6

3-7

3-8

3-9
3-10
3-11

3-12

3-13

3-14

LIST OF FIGURES

Major components of the distillation
process

Linde Column instrumentation

2 input, 2 output system with time delays

Model order estimation with different
amounts of additive noise

Input delay test for y1 with 0% noise
Input delay test for y1 with 10% noise
Input delay test for y1 with 20% noise
Input delay test for y2 with 0% noise
Input delay test for y2 with 10% noise
Input delay test for y2 with 20% noise

Inputs and outputs of the Linde Column
Recorded reflux flow (data set 1)

Recorded concentration on the 57th tray
(data set 1)

Recorded bottoms concentration
(data set 1)

Recorded reboiler steam pressure
(data set 1)

Recorded reboiler steam pressure
(data set 2)

viii

11

34

56

58

59

60

61

62

63

66

67

67

68

68

70

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

3-15

3-16

3-17

3-18
3-19

3-20

3-21

3-22

3-23

3-24

3-25

3-26

3-27

3-28

3-29

3-30

3-31

Recorded bottoms concentration
(data set 2)

Recorded reflux flow (data set 2)

Recorded concentration on the 57th tray
(data set 2)

Recorded reflux flow (data set 3)

Recorded reboiler steam pressure
(data set 3)

Recorded concentration on the 57th tray
(data set 3)

Recorded bottoms concentration
(data set 3)

Results of the delayed polynomial
matrix test for 611(2)

Mean square residuals for G11 with a
time delay of 5

Auto-correlation of residuals of 4th
order model with a delay of 5 for 611

Response of estimated G11 to the re-
corded reflux flow of data set 1

A portion of the recorded 57th tray
concentration from data set 1

Order test for 612

Mean square residuals for 612 with a
time delay of zero

Auto-correlation of residuals of a 3rd
order model with a delay of l for G21(z)

Response of the estimated G21 to a

portion of the recorded reflux flow
of data set 1

A portion of the recorded bottoms
concentration from data set 1

ix

7O

71

71

72

72

73

73

81

82

82

84

84

85

85

87

88

88

Figure

Figure

Figure
Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

3-32

3-33

4-1
4-2

4-3

4-4

4-5

4-6

4-7

4-8

4-9

4-10

4-11

4-12

4-13

4-14

Response of estimated G22 to a portion
of the recorded steam pressure from
data set 2

A portion of the recorded bottoms
concentration from data set 2

Basic Internal Model Control structure
Filtered Internal Model Control structure

Response of 611(2) to a step load
disturbance

Response of 611(2) to a step load
disturbance, IMC control with Gc(z) -
1/G11(1) and no filtering

Response of 611(2) to a step load
disturbance, PI control

Response of G11(z) to a step load
disturbance, IMC control

Brosilow's feedforward Internal Model
Control structure

Revised feedforward Internal Model
Control structure

Simplified Decoupling with Revised
Feedforward IMC

Simplified Decoupling with Brosilow's
feedforward IMC

Critical Decoupling with Revised
Feedforward IMC

Critical Decoupling with Brosilow's
feedforward IMC

Response of the 57th tray concentration
to a step load disturbance

Response of the 57th tray concentration
to a step load disturbance in the
bottoms concentration

90

90

101

101

110

112

117

117

119

120

126

128

132

133

136

137

Figure

Figure

Figure

Figure

Figure

4-15

4-16

4-17

4-18

4-19

Response of bottoms concentration to a
step load disturbance on the 57th tray

Response of bottoms concentration to a
step load disturbance

Response of simply decoupled IMC
(Brosilow's version) to a step load
disturbance

Response of critically decoupled IMC
(Brosilow's version) to a step load
disturbance

Stable images of unstable zeros

xi

138

139

141

143

108

CHAPTER ONE

INTRODUCTION AND DESCRIPTION OF THE THESIS

1.1 Introduction

The research presented in this thesis has been motivated by the
technology gap that exists between "advanced control theory" and
industrial practice. The vast majority of literature on advanced
process control has concentrated on the derivation of various techniques
and the mathematical advantages and disadvantages of these techniques.
A survey of the literature will show that the few papers addressing the
application of advanced control achieve their results through lab scale
equipment or computer simulation. Very few papers investigate the
application of advanced control to full scale processes. In the last
few years several noted authors in the field have cited a need for

application oriented research [1]. This thesis addresses this problem.

In this thesis we will consider: lack of well defined process
models, non-uniform sampling, robustness of controllers, and the
economic incentive or disincentive for using advanced control
strategies. The specific application of this research is an industrial
scale distillation column at Dow Chemical Company's Michigan
Division. This type of process operation was selected because
distillation is one of the most important unit operations in the petro-
chemical industry. Also, distillation columns are one of the largest

energy consumers in a chemical plant and improved control has been shown

1

to be effective in reducing energy costs. For the column of this study,
reducing the energy consumption by one percent would result in a cost
savings of $30,000 in one year alone. Improved control could also be
used to increase the product quality of the column. Estimates predict
that a one percent reduction in the overheads product impurity would
produce a savings of $25,000 per year. Production supervision easily
agreed to support this research given these potential savings and the
admission by operating personnel that the current control showed
deficiencies. In addition, such research could develop solutions to
some of the practical problems involved in implementing advanced control
theory while solving the specific control problem presented by the
application. Also, a channel of communication between industry and
academia would be fostered by application oriented research. This

thesis describes such research.

1.2 Research Objectives

The overall objective of this research was to determine if advanced
control would be beneficial to a commercially successful industrial
scale distillation column. In other words, what potential does advanced
control have to improve control. Typically energy savings and greater
product quality of a distillation column are realized through improved
regulation of the .two product streams. Usually these streams are
dynamically coupled so that action taken to correct one will also affect
the other. This situation lends itself to a multivariable control

strategy which requires a multi-input multi-output model of the column.

From this perspective, the major goals of this research are:

1. Develop a column model (the major task of this
research).
2. Propose an advanced control scheme based on considerations of

the model and practical issues such as robustness of the control
and how well plant personnel would understand and therefore

accept a new control approach.

3. Test the proposed control scheme through simulation and compare

to conventional control.

The relatively straight forward approach outlined above was
complicated by problems that are not frequently discussed in the
literature. Due to the multiple stages of a distillation column, time
delay, often called deadtime, is an important parameter of the model.
The model constructed for the column is a matrix of linear discrete time
transfer functiOns. Model order estimation and time delay estimation
needed to be performed simultaneously. Also, the column could only be
studied as it produced a salable product, therefore meaningful data for
model building was not easily obtained. In addition, the measurement of
one of process variables had a fixed sampling rate while all others were
virtually continuously measured. This produced a multirate sampling
problem. Finally, part of the column exhibited unexpected non-minimum
phase behavior which led to a novel approach to the controller design.
The original contributions of this research contain solutions to these

problems.

1.3 Outline of the Thesis

The remainder of this thesis is divided into four chapters. In
Chapter Two, the basics operation and control of a distillation column
is introduced. A physical description of the Dow column and its control

system is also provided.

Chapter Three is concerned with the problem of identification and
parameter estimation of the model of the Dow column. A survey of the
literature regarding distillation column modelling is given. In
this chapter the problem of time delay plus model order estimation is
defined. Singular value decomposition is described and a method for
estimating time delay and model order is then developed which uses this
mathematical tool. Simulated data is used to show the ability of this
method. Real operating data of the column is presented next and the
method is applied to it. Other model order estimating techniques are
used to confirm the model order. Then parameter estimation is
performed. The resulting model is then discussed in the light of a

prior knowledge of the column.

Chapter Four deals with the design of an advanced control strategy
for the column using the model of Chapter Three. First the control
implications of the model are discussed. Next a survey of the
literature on distillation column control is presented. Then Internal
Model Control is introduced and a multirate feed forward version is
developed for the column. Also, a reduced order controller is presented
to handle the non-minimum phase dynamics in the column. Simulations of

the column under the standard proportional / integral control and IMC

are described and the results presented. Implications of the simulation
results are discussed while considering the practical matter of running

the column for production.

Chapter Five presents a summary of the results of the research
described in the previous chapters. Conclusions are drawn and

recommendations are made for future research.

1.4 Summary

The results of this thesis have been disseminated throughout The Dow
Chemical Company before the publication of this thesis. It has been
presented to a body of managers representing Dow's global operations and
to a conference of researchers and engineers working on advanced process
control from Dow's global manufacturing sites. Follow up research has
already begun at Dow Canada's Sarnia Ontario plant to further understand
the implications of the inverse response and its control in their Linde
type column. In addition, the technique developed for the design of IMC
for non-minimum phase processes has been adopted in the design of a
temperature controller for an extruder at Dow's Michigan Division. The
interest shown in this research and the additional research it has
spawned is a true measure of the success this work is as an applied

thesis.

CHAPTER TWO

DESCRIPTION OF THE PROCESS

The purpose of this chapter is to descibe the fundamentals of the
distillation process and how they relate to control. Also provided is a
description of the physical make-up of the Linde Column and its current

control system.

2.1 Fundamentals of Distillation

The purpose of a distillation column is to separate the chemical
components of a feed stream into more or less pure product streams. The
separation is based on the well known fact that pure liquids exhibit
different volatilities (tendency to vaporize). Thus, heat applied to a
mixture of substances will generate a vapor which is rich in the more
volatile substances and leave a liquid which is more rich in the less
volatile substances. If the vapor is condensed the remaining liquids

represent the purified products of the distillation process.

Typically a distillation column consists of a cylindrical vessel
containing a number of equally spaced trays inside, Figure 2-1. A weir
mounted on each tray maintains the liquid level on the tray. Liquid
overflowing the weir travels through a downcomer to the tray below. The
trays are equipped with a means to allow vapor to flow from below and
mix with the liquid on the tray. The bottom of the distillation column

is connected to a heating unit called a reboiler, which provides the

CONDENSOR

-a. coou~c
FLU”)

 

 

 

 

ACCUMULATOR
Z _ ; . OVERHEADS
f1::::3** ='§ . “ PRODUCT
____J
TRAY
|__ /
FEED
STREAM W;
L____w
____J
L____w
aorroms 4—"Efgge
PRODUCT

REBOILER

Figure 2-1 Major components of
the distillation process

energy to create a vapor flow up the column. A portion of the liquid at
the bottom not vaporized is removed from the column and is known as the
bottoms product. The top of the distillation column is connected to a
cooling device called a condensor, which condenses the vapors leaving
the column. A portion of the condensed vapors, called reflux, is
generally returned to the column to flow from top to bottom. What
remains of the condensed vapors is drawn off the column and is called
the overheads productor distillate. The feed stream of the column is

usually introduced near the middle of the column.

The operation of a distillation column is a series of heat and mass
transfer operations occurring on each tray. Consider the simplest
separation of a binary liquid mixture. Feed entering the column flows
down from tray to tray to the bottom. As it flows down it is contacted
by the vapor rising up. This contact initiates a heat and mass transfer
which results in the vaporization of the more volatile,1ighter,
component in the liquid and condensation of the less volatile, heavier,
component in the vapor. So as the feed travels down the column it
becomes more and more concentrated in the less volatile component of the
feed. The vapor becomes more and more concentrated in the more volatile
substance as it rises up the column. Above the feed tray the reflux
returned at the top of the column serves to mix with the rising vapor.
Each end of the column represents the extremes in concentration for both
components, the top for the lighter component, the bottom for the

heavier component.

2.2 Distillation Control Fundamentals

The variables used to control the concentrations of the product
streams are usually the reflux flow and the vapor flow (through the
amount of energy added to the reboiler). By far the most common control
strategy used in industry is to manipulate the heat added to the
reboiler to control the bottoms product concentration and adjust the
reflux flow to control the overheads product concentration. Simply put,
if the concentration of the lighter component in the bottom is too high
heat is added to boil it off, and if the concentration of the heavier
component is too high in the distillate more is returned as reflux to be
distilled again. This simple single input, single output strategy has
limitations since vapor used to control the bottom of the column must
reach the top of the column and has affects there. In the same way,
reflux used to control the top of the column must flow to the bottom and
have an affect there. Consideration of the column as a two input, two

output system can improve the control, which is part of this research.

2.3 Description of the Linde Column

The subject of this research was a distillation column located in
Dow Chemical Company's Michigan Division which is 12 feet in diameter
and 90 feet tall. The column performs a binary separation of styrene
and ethylbenzene. The distillation unit is equiped with 70 Linde type
trays giving rise to its name, "Linde Column”. Liquid is drawn off the

bottom of the column and pumped through a vertical reboiler using steam

10

as a heating medium. Vapor drawn off the top of the column is totally
condensed in six air cooled condensers. The feed flow is introduced

around the 45 tray.

A Dow designed process control computer manages the process
through direct digital control, therefore the process is fully equipped
with electronic instrumentation. Feed flow rate and reflux flow rate
are measured by orifice type flow meters. Steam pressure in the
reboiler is measured by an electronic pressure transmitter. Feed
temperature is measured by resistance temperature bulb. The
concentration of styrene on the 57th tray, which is used to control the
concentration in the overhead product, is measured by an online
refractive index analyzer. The concentration of ethylbenzene in the
bottom product is measured by samples taken in the reboiler loop and
analyzed by a gas chromatograph. The gas chromatograph sampling rate is
fixed at 8 minutes. All other measurements are continuous. Figure 2-2

illustrates the total instrumentation of the Linde column.

Operating data was logged on a PDPll/44 digital computer which
communicated with the Dow process control computer. The data was then
read to a flexible magnetic disk and transported to Michigan State
University where it was read into a PDP11/04 computer and transmitted to
the College of Engineering's Prime Computer. The bulk of the analysis
was performed on the Prime Computer. The data was also transmitted from
the PDP11/44 through a local area network to a VAX 11/785 for further

analysis.

11

CONDENSOR

 
  
 
 
 

AMBIENT
AIR

OVERHEADS
PRODUCT

FEED

 

0

{9
”EA"- I—— Q ©~
Mel—.1 @m

I [1: ® T°mp°'“‘"'°
@

ii;.T-%%§3;N§' Flowcontml

 

 

 

 

 

Bor'roms
PRODUCT 9

> 4 G
STEAM
L__J

REBOILER

 

 

Figure 2-2 Linde Column instrumentation

12

2.4 Current Control of the Linde Column

The current control strategy to maintain product quality is based
on digital Proportional-Integral controllers. The bottom composition is
regulated by changing steam flow to the reboiler through a cascade
control scheme. To maintain bottoms composition, the set point of a
slave steam pressure controller is adjusted. This slave controller in
turn regulates the flow of steam to the reboiler through a pneumatically

operated control valve.

The top product is also controlled by a cascading of control loops.
The overhead composition is maintained by regulating the concentration
on the 57th tray by indirectly manipulating reflux flow through changes
in overhead product draw off. Changes in over head product flow are
transformed into reflux flow changes through the overhead accumulator
which is level controlled by reflux flow. Although these controllers
are slightly more complex than the conventional control scheme already
described, the Linde Column controls still use steam to control the

bottom product and reflux to control the top product.

2.5 Summary of Chapter Two

This chapter has shown distillation to be a complex chemical
process. Even so, the generally applied control scheme of industry is
simple independent PI controllers applied to each product stream. In

the case of the Linde column, much more complex schemes can be applied

13

due to the full complement of electronic instrumentation coupled with
direct digital control. This situation offers good opportunities for

the application of more advanced control strategies.

CHAPTER THREE

LINDE COLUMN MODELLING

The purpose of this chapter is to develop a dynamic model of the
Linde Column. The model is a 2 input, 2 output matrix of discrete time
transfer functions. A survey of the literature on distillation column
modelling is given first. Then as a preliminary, the problem of
determining the structure and order of a multi-input, multi-output model
is addressed in the general sense. A new solution is presented. It is
then used on real data from the Linde Column. Parameters of the model
are determined by least squares estimation. Finally the identified

model is discussed in light of the operating experience of the column.

3.1 Literature Survey of Distillation Column Modelling

Distillation columns have been popular subjects of identification in
the literature. Industrial columns as well as pilot scale columns have
been considered. As in this research, several authors have derived
discrete time models of distillation columns using operating data.

However the research in this thesis differs significantly from each.

For example, Williams [3] modeled a six plate pilot scale column

using step and Pseudo Random Binary Signal (PRBS) inputs. He only

considered the dynamics between reflux flow and top product composition

14

15

and did not consider the steam flow to the reboiler as an input as was

done in this research.

Krishnamoorthy and Edgar [4] used pulse testing to model a pilot
scale column as a two input two output system. They
used least squares estimation and employed several methods to determine
model order. Their main concern was to properly prefilter the sampled
data before estimating the model order. Time delay was estimated by
shifting the data and then correlating the input output data and by
finding the time lag which minimized the resulting estimation error.
They concluded that prefiltering resulted in consistant estimation of
model order by the different techniques, but failed to report the models

they derived.

Foulard and Bornard [5] report on both steady state and dynamic
modelling of pilot scale and industrial distillation columns. Different
estimating techniques were used to develop low order dynamic models with
no time delay. Correct model order was determined by finding the order
that minimized the estimation error. Unlike this research, only single
input single output relationships were considered. Main conclusions
were that the type of input signal greatly affects the results of the
modelling procedure and that there exists a significant problem in

obtaining good experimental data from an industrial column.

A full scale column is modelled by Gauthier and Landau [6]. The
input applied to the column is two uncorrelated PRBS signals added to

the set points of the overhead and bottom temperature control loops,

16

rather than perturbing the inputs to the column as was done for this
research. A two input two output model was identified using a method
similar to that used in this research but time delay was not considered.
The instrumental variable method of estimating as well as Landau's
output error estimator were used to fit the model to the experimental

data.

Gustavsson [7] surveys other applications of discrete time modelling
to industrial scale distillation columns. The applications of his
survey that are the most similar to the present research are still
deficient in some area. Most did not evaluate any control strategy that
the derived model might have suggested. Some applications did not model
the column under study as a two input two output system. Other
applications considered tray temperatures rather than product
compositions as outputs. All authors did find modelling of an
industrial scale column a problem that is compounded by the lack of

control the experimenter has on the column as a whole.

3.2 Current Model Identification Methods

A prerequisite to the estimation of parameters of a dynamic model
of a process is the knowledge or at least an estimate of the model
structure and order. In some sense the identification of model
structure and order is a much more difficult task than the estimation of
parameters. In practice, the model structure and order are determined

partly through the existing knowledge of the process and partly through

17

a combination of statistical tests. The original tests of model order
were based on statistical evaluations of estimated models. Although
these tests are reliable, they require large amounts of computing
because all competing models must be estimated in order to apply the
tests. In recent years several techniques to estimate model order prior
to parameter estimation have been presented in the literature. A very
recent comparison of the most popular techniques is given by Freeman
[15] . The obvious advantage of these tests is that they significantly
reduce the computational burden of model development by eliminating the

need to estimate all competing models.

The aim here is to extend an existing multivariable technique to
include explicit consideration of an important structural parameter of
discrete linear models, namely time delay. The importance of time
delay in a model used for control is illustrated by the numerous
techniques that exist to compensate for it. In terms of modelling
chemical processes, time delay occurs very often because of actual
transport delays that exist in the process and because it is an
effective means to approximate higher order dynamics. Therefore any

effective modelling procedure must consider time delay .

Numerical robustness is another important attribute of an
identification procedure. This is because the computations of any
algorithm will be performed on a digital computer which has finite

accuracy limits. This is also addressed in this thesis.

18

3.2.1 Identification Literature

There are several statistical tests of the estimated model that have
seen wide spread use. Astrom [16] suggested using the statistical F-
test on models of increasing order. Tests based on the sum of residuals
squared or the sum of the reconstruction error are also popular, see
Gustavsson [17]. Akaike [18,19] has suggested two tests that combine
the variance of the residuals with the number of parameters into a
single measure. Tests based on the auto-correlations of the residuals
are also widely used. Although these tests can only be applied after a
model has been estimated, they are still useful in validating the model
structure which is selected by methods applied to the data prior to

parameter estimation.

Tests applied prior to parameter estimation are based on determining
the linear dependence between input and output data. Lee [20] was the
first to exploit this relationship by testing the singularity of the
product moment matrix. Woodside [21] then extended this idea by
developing three measures testing the singularity of the product moment
matrix. Wellstead [22] modified this idea by basing a test on an
instrumental product moment matrix and also included a procedure to

estimate time delay in a single input single output model.

The existing multivariable techniques for model identification,
suggested by Budin [24] and Guidorzi [25], do not explicitly estimate
the time delay of the model. Instead they estimate a model order which

is large enough to include the time delay. This results in a model

19

with many parameters which have values identically equal to zero. In
principle these methods rely on the parameter estimation scheme to
detect the magnitude of time delay by the number of parameters with
estimated values close to zero. In practice this methodology can break
down when applied to noisy data. A more reliable means to estimate time
delay would result if the zero valued parameters could be identified
prior to parameter estimation and then confirmed by other statistical

tests applied to the estimated model.

3.3 Delayed Polynomial Matrices

A method to identify multi-input, multi-output models with time
delay prior to parameter estimation is developed as follows. First a
theorem is proven which relates a system's state space realization with
input delays to a canonical input - output realization. Next the
results of the theorem are used to develop a method which estimates
model order and time delay in multivariable systems. Then singular
value decomposition is introduced as a means to detect near singularity
of matrices and it is then used as a measure to determine model order
and time delays in multivariable systems. Simulated data is then used
to demonstrate the ability of the method. Finally the method is applied

to real data from the Linde column.

20

Theorem 3.1
Let a completely observable discrete time system with input delays

be described by the following state space equations

x(k+1) - Ax(k) + Blu(k-d1) + Bzu(k-d2) +...+ Bsu(k-ds) (3.1)

y(k) - CX(k) (3.2)

where x e R“, u e Rr, and y e Rm.

Then there exists an equivalent input-output description of the

system with the following form

P(z)y(k) - Qlu(k-d1) + Q2u(k-d2) +...+ qu(k-ds) (3.3)

where P(z) and Qa(z), a - 1,2,...,s, are polynomial matrices whose

entries satisfy the following relations

deg{pi1(z)} > deglpij(z)}, j > i (3.4a)
deg{p11(z)} z deg{pij(z)), j < i (3.4b)
deg{p11(z)} > deg{pji(z)}, j i i (3.4c)

deg{p11(z)} > deglqa’ij(z)} (3.4d)

21

The degree of each of the polynomials is directly related to the

structure of A and C.

Proof:

First consider C as a matrix of row vectors

(3.5)

o-—<1 c:
Braror+4~i

Now the observability matrix may be constructed as

(n-l)
o - [cT, AICT, A: CT, ..., AT cT] (3.6)

And the columns of 0 searched for linearly independent columns in the

following order

(3.7)

It is well known that a equivalence transformation matrix T may be

constructed by arranging the independent columns in the following way

22

T T(”1'1) T(”m'1)

T - [01, A cl, ..., A 01,02, ..., A cm] (3.8)

Application of the transformation matrix T produces a new system of

equations

w(k+l) - A*w(k) + B:u(k-d1) + B:u(k-d2) +...+ B:u(k-ds)

(3.9)

y(k) - c*v(k) (3.10)

where w - Tx. The form of the matrix T imposes a special structure on

* *
the matrices A and C ,

A* - TAT'l - {A:j}, where 1,3 - 1,2,...,m and
0
* Q I
A11 - 6 u1_1 (3.11)
a ,a , ...,8
11,1 11,2 ii,v1 (V1 x vi)
0 0 0 0 ... ... ... ... 0 0
A? - 9 (3.12)
ij 0
aij,l’ aij,2’ ..., aij,v , 0 ... 0 (vi x u j)

13

23

,, _1 100 ...................... o
c - cr - 0 o 1 0 o ............. o (3.13)
60 00 ..... 0100 0
t t t
1 (u1+1)...(u1+ +um_1+l)

*
We see that C is composed of row vectors which have only one non-zero

entry (equal to unity). The position of this non-zero entry is shown
for the first, the second and the mth row. Meanwhile 3:, a - 1,2,...,s,

has the generic form

* * *T
b ........ b b
B: ;a,1l ;a,lr _ ;a,l (3.14)
* '* '*T
. ........ b b
“.111 a,“ a!

It is apparent from the structure of A* that the original system has

been decomposed into m interconnected subsystems. Because of the
complete observability of the system it follows that v1 + v2 ... + "m -
n. In fact, the integers Vi define the structure of A*and C*and are

invariant with respect to changes in state space coordinates, hence the

name "Kronecker invariants".

Consider now the vector w(k). For conciseness of notation define

i;k - u + . + y + k. Now w(k) may be expressed as

1 " i

24

 

 

V(k) - [ W1(k) 1 - [ w0.1(k) 1
wyik) l l w1;o(k) I
wu1+1(k) w1;1(k)
"u +u (k) l ’ w2;0(k) I
wv1+u2+1(k) "2;1(k)

I. "“00 .l I. wn<k) J

Using equations (3.9) and (3.10), the following equations may be written

for the ith subsystem

w1-1;1(k) - y1(k)

S
*1
"1-1;2(k) ' zy1(k) ' gglba,1-1;1 “(k-dc)

25

S
2 *T
”1-1;30‘) " z y1(k)’ 221%,1-13 “(k‘da)

S
*T

a-l
ui-l
w1;0(k) - Z Yi(k) ‘
u -2 s
i *T
2 E: ba,i-1;1 u(k-da) -
a—l
s
*T
- Z ba,i;-1 u(k-da) (3.15)
a-l

where yi(k) is the ith component of the vector y(k). Comparable

equations for each of the m subsystems can be written in a compact form

as

8
w(k) - V(z)y(k) - leac(z)u(k-da) (3.16)
a.

where

 

 

 

 

 

 

' 1 ......... o l
. z 0
V(z) - . ui_1 1 (3.17)
z z
o 'u _
_ 0 o ..... 0 z m 1 _
’ 0 0 . 0 0 0 ‘
*1
a,1 o ..... 0 o o
a _ ; - , (3.18)
b*T ...'b*T 0 0
a,v1-l a,l
0 0 0 o 0
*T
a,n-u +1 0 ° ° 0
. . m
1*1 *r
_ ba,n31 ..... a,n-v +1 3 [n x r(uM-l)]
I
21
0(2) - : (3.19)
"1'1 [r( - 1) x r]
z I "M
VM - maxi{vi)
Substitution of (3.16) into (3.9) produces the input-output equation
S
(21 - A*)V(z)y(k) - E: [ (zI - A*)W&G(z) + a: ]u(k-da) (3.20)

a-l

27

Of the n equations in (3.20) only m are significant, namely the V1
th, the (v1 +u2)th, ..., nth. These m significant equations can be

expressed as

S
P(z)y(k) - 2: Qa(Z)u(k-da) (3.21)
a-l
where
P (2) ... p (2)
P(z) - ;11 :1“ (3.22)
pm1(2) ... pmm(2)
q (Z) -.. q (Z)
Qa(z) - ga'll ga'lr (3.23)
qa’m1(z) "° qa,mr(z)

The polynomials of P(z) can be written directly from (3.20)

P110!) - Z - aii’y Z ' ... ' 811,22 " 811,1 (3.24)

pij(z) - - aij’yijz - ... - aij,22 - a”,1 (3.25)

Straight forward computations lead to the following forms of the

polynomials of the 00(2)

28

qa,ij(z) ' fia,i;0,jz + '°° + fla,i-1;1,j

(3.26)

Here the coefficients fla ij are obtained from the matrices

- MB - I I (3.27)

where the matrix M is given by

-a11’2....-a11’u1 .1 alm,2""' 1m,v1 1
’811,3 "11,4°':° .... a1m,3 ”‘1m,4":'
M - I .° 2 .'

-; 1" .... A 1'

111,3!1 1lm,v1m

-aml,2 aml,3:m .1 mm,2 -amm,v 1

-a a .... -a ...I'
ml,3 ml,4 . mm,3 mm,4 .

-$ 1 ' A 1
ml,v mm,u

- 0 m1 1 m J

 

 

Due to the structure of the matrix A* the following relations hold

for the polynomials of P(z) and the 00(2):

29

deg{pii(z)} > deg{pij(z)), j > i (3.29)
deglpii(2)} 2 deglP1j(Z)l. j < i (3 30)
deg{pii(z)} > deglpji(z)}, j fl 1 (3.31)
deg{pii(z)} > deg{qa,ij(z)} (3.32)

The validity of the theorem has thus been established.

The input-output description of (3.21) may be modified in a way that
is appropriate for the identification of the system from input and
output data. This modification is introduced by way of the following

lemma.

Lemma 3.1

The representation of (3.21) can be equivalently stated by a

composite input-output equation, the so called polynomial matrices

r°<z>y<k> - e°(z)u(k) (3.33)

where Pc(z) and Qc(z) are matrices of polynomials whose entries satisfy

the relations of (3.29) - (3.32).

Proof:

d
Multiplication of (3.21) by z 3 produces

30

ds ds-d1 . dS-d2
z P(z)y(k) - z 01(z)u(k) + z 02(z)u(k) +...

...+ Qs(z)u(k) (3.34)

By making the assignments

2°(z) - z SP(z) (3.35)

c ds.d1 ds.d2
Q (2) - 2 01(2) + 2 02(2) + ... + 08(2) (3.36)

equation (3.33) results.

Now since the degree of each element of Pc(z) is just the degree of

the corresponding element of P(z) increased by ds’ equations (3.29) -
(3.31) hold. The degree of any element of Qc(z) can be at most equal to
the degree of the corresponding element of 01(2) increased by ds - d1,

so equation (3.32) holds. The lemma has thus been proven.

The polynomial matrix description of (3.33) was first introduced by
Guidorzi [25]. The development found in [25] inspired its extension to

systems with input time delay found in this thesis.

31

REMARK 3.1

Often the identification of a model of a linear multivariable system
with multiple input delays is made from input and output data. It is
natural to consider an input-output description like that of (3.21) or
(3.33) coupled with the relations of (3.29) - (3.32). However these two
forms make inefficient use in the number of parameters. This can best

be seen by considering an arbitrary 2 input, 2 output system with 2

input delays.

[Sum p12<z>][::1<k)] _ [3141(2) q1,12(2>][u1(k-d1>] ,
21(2) p22<z) 20c) q1,21<z> q1,22<z> u2<k-d1>

[qznm q2,12<z)][u1<k-d2)]
q2,21(z) (12,2205) “2(k'd2)

Neglecting for the moment q1 12(2) and q2 12(z) the first equation

represented by (3.33) would be

u +d v +d -1 d d -1
1 2 1 2 2 2
(z -a11,ylz - ... -a11’lz -Oz -...- 0)y1(k) +
u +d u +d -1 d d -1
12 2 12 2 2 2
(z ‘312 ’Vlz - o o c -812’12 -oz - o o o - o>y2(k) -
u +d -1 u +d -d +1 V +d -d
l 2 l 2 l 1 2 1
(02 + ... +02 +fi1’y112 +...+52’11)u1(k)

A more compacted form for writing this equation would be

u +d -d d -d d -d -l
1 2 1 2 1 2 1
(z - ...-all’lz -Oz -...- 0)y1(k) +
v +d -d d -d d -d -1
12 2 l 2 1 2 1
(z - °"-al2,lz -Oz -...- 0)y2(k) -
u +d -d -d
1 2 l 1
(£1 V112 +...-+192 11)z u1(k)

where the d trailing coefficients of p11(z) and p12(z) and the d

1 1

leading coefficients of q11(z) have been taken care of by the delay term

-d
z 1. Of course this same reasoning can be applied to u2(k) and its

relation to y1(k) as well.

This more compact representation can be applied to each input/output
pair of a multivariable system with an arbitrary number of inputs and
outputs. Although the arguments made above only show the existance of

delay terms of magnitude equal to dl, in general a search can be made to

construct a delay term for each entry of 0(2) which is at least d1.

This produces the following delayed polynomial matrix description which

generally has a smaller number of parameters than (3.33), or (3.21)
* *
1’ (2)7(k) - Q (2)1100 (3.37)

where P*(z) has the same form as in (3.33) and 0*(2) ia a matrix of

dimension m x r of polynomials whose entries have the form

33

...-<11, 23b
qij(2) - bij a (3.38)

Equation (3.37) represents an input-output description which can
have a fewer number of parameters for multiple input, multiple output
systems with time delay than the classical transfer matrix description
or the polynomial matrix form of (3.33). This is especially true for

systems which have common modes in the outputs.

In the next section a procedure will be developed which will
estimate the order of the polynomial entries and delay terms of (3.37)
from input / output data. The following numerical example will

illustrate the usefulness of such a representation.
EXAMPLE 3.1

Consider the 2 input, 2 output system shown in Figure 3-1. A state

space equation that describes this system is

.8 O O 0 O 1 O O O
O .6 0 O O 0 1 O O
x(k+l) - 1 0 O 0 0 8(k) + 0 O u(k-1) + 0 O u(k-2)
0 l 2 O 0 0 0 O O
O 0 O 0 .3 0 O O 1
O l 1 O O

The observability matrix of this system is

34

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2
._.. -3
'z - 8 Z
1 1
' . z - 6 z
22 - .9

 

 

 

 

 

-2
u2'—’(z-.6) (z-.3)“’ z

 

 

 

 

Figure 3-1 2 input, 2 output system ..
with time delays

0 0 l 0 .8 2
l O .6 1 .36 .6
O - 1 0 0 2 0 0 ......
0 l 0 0 0 O
O 1 0 .3 0 .09

Looking first for the vectors of 0 associated with y1 produces the

following transformation matrix

0 1 1 0 0
l .6 0 0 0
T - .8 .36 0 0 0
0 0 0 l 1
0 1 2 0 .3

Here v1 - 3 and v2 - 2. Application of T via equations (3.11) and

(3.13) gives

0 1 0 00
* _1 o 0 1 00
a - TAT - 0 -.48 1.4 0 o
o o 0 01
o 1 2 0.3
o 1 00
* 1 .6 * 00
32- r32 8 .36 32 - 1'32 - 0 0
0 0 01
o 1 03

The matrix M is therefore given as

.48 -1.4 l 0 0

-l.4 l 0 0 0

M - l 0 0 0 0

0 -2.5 0 - 3 1

-2.5 0 0 1 0

which produces

-.6 0 0 0
l -.8 0 0
E: - an: - 0 1 E: - usz - 0 o
-2.5 -.5 0 0
0 - .5 0 1

The equivalent input-output description of (3.21) is therefore

23 - 1.422 + .482 + 0 0 y1(k)
.32 + 0

-2.522 + 02 - .6 22 - y2(k)

2
2 - .6 2 - .82 (k-l) 0 o (k-2)
[ 02- 2.5 -2.52 - .5 ] [ :;(k-1) ] + [ o 2 ] [ 3%(k-2) ]

The number of parameters of this representation is 22. The polynomial

matrix description of (3.33) is given by

24 - 1.423 + .4822 + 02 + 0 0 y1(k)
-2.523 + Oz2 - .62 + 0 23-.322+Oz + 0 y2(k)
023+2022+ 2 - .6 023 22 - .82 + 0 u1(k)

02 + Oz - 2.5 Oz -2.Sz + .5 u2(k)

The number of the parameters for this representation is 26.

37

As a contrast the delayed polynomial matrix representation of (3.37)

is
22 - 1.42 + .48 0 J [ y1(k) ]
-2.522 + 02 - .6 22 - .32 + 0 y2(k)
[ (z - .6)z'?1 (z - .8)z-]_’1 ] [ u1(k) ]
(Oz- 2.5)z (-2.52 + .5)z u2(k)

Here the number of the parameters is 16.

Example 3.1 demonstrates how representation of a system by delayed
polynomial matrices reduces the number of parameters necessary to
describe the system behavior. This can lead to more consistant
parameter estimation when noisy input-output data is used to fit the
parameter values. A method which identifies an delayed polynomial
matrix model of a system from input-output data prior to parameter

estimation is described next.

3.4 Structural Identification

In this section a method will be described which deduces the
structure of the delayed polynomial representation of a discrete system

from input and output observations. The term structure relates to the
* *
order of the polynomial entries of the matices P (z) and Q (2) and to

the magnitude of the delay terms of the entries of 0*(2) found in

38

equation (3.37). The method that follows can be applied prior to
parameter estimation and is suitable for interactive computer use. The

method is a two phase procedure:

PHASE 1 Estimate the polynomial matrix

representation, equation (3.33)

PHASE 2 Estimate the delay of each input
to each output to form the delayed
polynomial matrix model, equation

(3.37)

3.4.1 Phase 1

In order to identify the polynomial matrix model of a system the

Kronicker invariants, Vij’ must be estimated. The structural

identification will be performed by exploiting the linear dependence
relations imposed by (3.33) and the structural relations imposed by
(3.29) thru (3.32).

Consider first the i-th row of (3.33). It may be rewritten as

r V

1]

y.(k+u ) - X: a y (k+1-1) +
1 1 ‘12-]- 1_1 13,1 .1

39
”1

32);.

j-l _1 (u1+...+yi-1+1),juj(k+1‘1) (3-39)

where v - v ..
ii

Consider now the vectors of observations given as

ni(k) - [yi<k>. y,<k+1>. .... y,<k+N-1>1T (3.40)
T
nj(k) [uj(k), uj(k+1), ..., uj(k+N-l)] (3.41)

Because of (3.39) the following is true

r V

11
q (k+y ) - a q (k+1-1) +
i i §;1 121 ij,1 j
m "i
j-l 1-1/31_1;1’ij(k+1-1) (3.42)

Equation (3.42) suggests the following selection plan to determine

each u :

1.1

1. Select the observation vectors in the following

order

4O

r2100. 112(k). nr(k). mlck). «mm. 711(k+1). n2<k+1>.
nr(k+1), nl(k+l), ..., nm(k+l),n1(k+2), .....
(3.43)

2. Retain a vector if and only if it is linearly

independent from the previously selected vectors.

3. When a dependent vector ni(k+ui) is found it is
no longer necessary to test the vectors ni(k+yi+j), j > 0,

since they will also be linearly dependent. At this time,

note the selected sequence of vectors.

4. Continue until a dependent vector has been

found for each output yi. The Kronecker

indices, u have now been identified.

ij’

3.4.2 Phase 2

The structure of the delayed polynomial matrix representation can be
deduced «from input / output data by using the polynomial matrix form as
a starting point. To see this consider the ith row of equation (3.37).

It may be rewritten as

41

y*
r Vij

*
yi(k+ui) — §;H1-aij 1 yj (k+l- 1) +

*
7ij
j-112:b*ij 1 uj(k+l-l-rij) (3.44)

One can remove the time delay in the inputs by equivalently representing

the system in prediction form

*
ij

*

yi(k+yi+£mu1¥;:*aij 1 yj(k+l'1+£max,i) +

u<k+1-1+£max (3.45)

'7

,i

i“
3'1

Iii:
3-11-1

bij,1u 13)

where €max,i - max{£ij}, £13 - fij + 113. Furthermore, this system can

be equivalently represented by

y(k+v1) -
i: (a zyij;1 +a zsmax’1+Oz€maX}1+...+0)y (k)
J_1 ijvyij 1131 J
m v -1 v r -1
+ §:1(02 1 + +bij’71§ 1 ij+ ......
v -€ v -5 -1
..... +b i ij+02 1 13+ ..+0)u (k) (3.46)

J

42

'k

*
Where V - V ij + €max,i'

i i + g and v . - v

max,i ij
This last equation is the same as the ith row of the polynomial matrix

form of equation (3.33).

Equations (3.44) - (3.46) suggest that the delayed polynomial matrix
representation may be deduced from the normal polynomial matrix
representation by identifying the zero valued parameters of (3.46). The

following steps accomplish this identification.

STEP 1. Set 7 - 0, for i - l, 2, ..., r.
max,i

Set pmax i - 0, for i - l, 2, ..., r.

STEP 2. Select any output yi and start with the collection of vectors
which produced the linearly dependent vector qi(k+ui). Select any
output uj and remove groups of vectors from the collection in the

following order

"1(k)9 0"! "r(k)’ ~J(k+V1-1);
n1(k+1),..., nr(k+1),nj(k+ui-2);
..... (3.47)
STEP 3. Continue to remove vectors until it is found that qi(k+yi) is

no longer linearly dependent on the remaining vectors. If the removal

of n

3
-d

1'.

13

STEP 4

(k+u

. - l and
J

43

i-dj) eliminates the linear dependence of ni(k+ui) then set

- max(r

.}

r . ,r.
max,1 max,i 1]

. Repeat step 3 for all other inputs u1 by removing vectors from

the collection in the following order

STEP 5.

groups

order

- max,i);

"1(k+fmax,11 ), ., n (k+fmax,11 ),Ic1(k.+yi-rmax 1 );

rm

01(k+rmax’12 ). 2 );

Hon

., "r(k+fmax,12 ),n1(k,+u1-rmax

(3.48)

Having determined 7 1, now select any input u.j and remove

of vectors from the remaining collection using the following
~J(k).
xj(k+1),
nj(k+r ,1);
"1(k+rmax I1 ), , nr(k+r i1 ),xl(k+rmax 11);

44

"1(k+fmax 12 )’ ’ "r(k+'max 12 )’nl(k+'max 12);
...... (3.49)
STEP 6. Continue to remove vectors until it is found that ni(k+ui)

is no longer linearly dependent on the remainder. If the removal of ”j

(k+p ) eliminates the linear dependence of n (k+u ) then set 7. - v -
j i i 1j i

Tij - pij+ l and pmax,i - max(pmax,i’pj}'

STEP 7. Continue steps 5 and 6 for all other u1 by removing groups of

vectors in the following order

"1(k+'max,ipmax,i)’ ""

"r(k+’max,1pmax,i)"l(k+fmax,ipmax,i);

...... (3.50)

- p 7+ 1.

*
until all 7ij are identified. Set u - max 1

1 ”1 ' 'max,i

STEP 8. Go to step 2 and repeat for any previously unselected output.

45

For the above selection plan to work two conditions on the input
must be met. First the length of the observation vectors must be large

enough to accomodate (3.42) so

H

Also the vectors of input observations must be linearly independent.
This is the same thing as requiring that the input to the system be

persistently exciting.

3.5 Computer Implementation

The success of the selection plan discussed above rests on
accurately testing for the linear dependence among the observation
vectors. This problem, which is the same as determining the rank of a
matrix composed of the observation vectors, is well understood
mathematically but as a practical problem it remains a challenge. In
the present context the challenge lies in performing the mathematical
calculations on a finite precision machine using noisy input and output
observations of a dynamic system. As it shall be seen the errors
introduced by the computer are handled by a stable numerical algorithm,
singular value decomposition. While the noisy data is handled by proper

filtering techniques.

46

There exists several mathematically equivalent techniques to
determine the rank of a matrix. Unfortunately, different approaches

lead to computational methods which can give different results when

implemented on a digital computer. Several reasons exist for these
discrepancies. First the computations must be performed in bounded
arithmetic, bounded by finite precision and finite range. So

computations are executed with truncation errors and must lie within the
range of the computing machine. Furthermore, the mathematical problem
posed is an inexact representation of the underlying physical system.
In other words the linear relationship assumed among the observations is
only an approximation to the usually true nonlinear nature of the
process which produced them. For these reasons the best that we can
hope for is a numerically stable algorithm which quantifies the
closeness to singularity of a matrix composed of the set of considered

vectors .

3.5.1 Singular Value Decomposition

Singular value decomposition, SVD, is generally acknowledged as the
most reliable method for determining rank numerically [26]. As a
mathematical tool it has been around since the late nineteenth century.
Golub and Reinsch [27] developed an efficient numerical algorithm for
computing the SVD in 1970. SVD has seen some use in the control
literature since the middle 1970's. However it was overlooked by past
authors of techniques to estimate model order by testing the linear

dependence of the observations. In fact, the reason SVD was used in

47

this research is because of the poor performance that was experienced
when the competing methods of rank determination, namely determinant
calculation and eigen value calculation, were first tried. It is likely
that the application reported here will soon be one of the many uses of

SVD in the control field.

Singular value decomposition for real matrices can be defined in the
following way. Let A be a real matrix of dimensions n x m with rank r,

where n z m z r. Then there exists matrix, 0, matrix, V, and matrix, S,

such that
T
A - usv (3.51)
s1 0 0 0 0
S - 0 ' . 0 (3.52)
0 0 0 0' s
m
The numbers 51,...,sm are called the singular values of A and by
convention s z s z ... 2 s .
l 2 m

The well conditioned nature of the singular values is illustrated by

the following property. If A + F has singular values {1....,§m, then
[31- (1] S [IFII - {1. Thus, singular values are well conditioned with

respect to perturbations in a matrix of noisy input / ouput data.

Ideally, if r < m then, s z 32 2...: sr and s - ..... - s - 0.

l r+l m

80 the number of non-zero singular values gives the rank of the matrix.

48

However, 31 is clearly sensitive to scale so a better dimensionless
quantity is sr / s1 , which is sometimes called a condition number of

the matrix. Calculation of the condition number of a matrix of
observation vectors is useful when the values of the inputs and outputs

differ greatly in magnitude.

Even by using singular value decomposition one would be ill advised
to operate directly on the matrix constructed from the observation
vectors of equations (3.43), (3.47)-(3.50) because of its large

dimensions. A more efficient algorithm based on the principle that

rank(A) - rank(AIA) follows.

3.5.2 Computer Algorithm

Some notation must be explained before the computer algorithm to
determine the delayed polynomial matrix representation is presented.
Consider a single input, single output system and let a matrix of

observation vectors be designated as

0(#1=31.#2=32) - [0(k+31).n(k+31+1). .... n(k+p1).

n(k+82). n(k+62+1), ..., n(k+p2)] (3.53)

Now let

49

T
A(p1:81,p2:62) - O (p1:61,p2:62)0(p1:61,p2:62) (3.54)

Therefore a matrix of observation vectors from a system with m inputs

and r outputs would be designated as

0(#1=31. .... ”rzar’ ....#r m ar+m) - [n1(k+61). .... 01(k+#1).02(k+32).
., nr(k+pr), 21(k+ar+1), ..., nm(k+pr+m)] (3.55)
and
A(p1:61, ..., ”r+miar+m) -
0T( -a -a )o( -a -a ) (3 56)
”1' 1’ "" ”r+m' r+m “1' 1’ "°’ ”r+m' r+m '

The two phase method for the identification of delayed polynomial
matrices may now be implemented by the following steps on a digital

computer:

3.5.2.1 Phase 1

STEP 1. Construct the sequence of increasing dimension matrices

A(1:0,0:0,...,0:0),A(1:0,1:0,...,O:0),...,A(1:0,1:0,...,1:0)A(2:O,1:0,..

.,l:0),A(2:0,2:0,...,l:0),...,A(2:0,2:0,...,2:0)A(3:0,2:0,...,2),....

(3.57)

50

STEP 2. Test for the singularity of each matrix by calculating its
smallest singular value or its condition number. When a singular matrix
is found there will be a large drop in the evolution of the smallest

singular value or the condition numbers.

STEP 3. When a singular matrix is encountered, say A(p1:0,...,vi

:0,...,p :0), v has been identified. Let ”i - u - 1 and remain

r+m i i

constant. It also follows that Vij - ”j'
STEP 4. Continue to increase the other indices in the same manner until

all observability indices u are found, noting the structure of the A

1
matrix associated with each ”1'
3.5.2.2 Phase 2
STEP 1. Set 1 - 0, i - l, 2, ..., r.
max,i
Set pmax,i - 0, i - l, 2, ..., r.
STEP 2. Select any output yi and start with the singular matrix
A(v11:0,...,ui:0,...,vir:0,vi-l:0,...,v1-1:0) which identified its Vi .

Select any input u and construct the sequence of decreasing dimension

J

matrices

51

l i r r+l r+j r+m

i J» i l .1 Jr
A(uilzl,...,v1:1,. ,uir:1,vi:0,. .,vi-2:0,...,vi:0),
A(v11:2,...,vi:2,. ,Virz2,ui:0,. .,Vi-3:O,...,Vi:0),
A(vil:3,...,vi:3, ,virz3,ui:0,...,vi-4:0,...,ui:0),
..... (3.58)

STEP 3. Test for the non-singularity of each matrix by calculating its

smallest singular value or its condition number. When a non-singular

matrix is found, say A(v11:dj,...,u1:dj,...,uir:dj,ui:0,...,ui—l-dj:0,..
..,u1:0),dj:0,...,v1:0), there will be a large change in the evolution
of the smallest singular value or its condition number. Set 'ij - dj
- l and 'max,i - maxlrmax’i, 'ij"

STEP 4. Repeat steps 2 and 3 for all other inputs, ul, using the

sequence of decreasing dimension matrices

1 r r+l r+j r+m

i t J i t
A<Vilzrmax,i’""yir:'max,i’yi:o" .,vi-2:0,...,v1:0),
“(”11"max,1'°"'”1r"max,1'”1‘°"'°'”1'3‘°'°'°'”1‘°)'
A<Vil:'max,i""’Vir:fmax,i’yi:0’° ’Vi.'max,i:o’ ,u1:0),
A(Vil:'maxti’°'°’yir:fmaxT1’Vi:o’°"’yi"maxfizo""’yizo)’

52

STEP 5. Having determined 7 now select any input uj and

max,i’

construct the following sequence of decreasing dimension matrices

l r r+1 r+j
4 t t i
A(Vilzfmax,i"'"Vir:rmax,i’yi-ri1:o""’Vi-Tij:1""
.,ui-rir:0)
A<yil:fmax,i"°°’Vir:rmax,i’Vi-'ilzo"°°’Vi'rij:2’°"
...,Vi-rir:0)
A<Vilzrmax,i"'"yir:7max,i’Vi-rilzo"°°’Vi-fij:3’°°'
...,vi-rir:0)

A<Vilzfmax,i"°"Virzfmax,i’yi-ril:o’""yi-fij:fmax,i"'°
.,ui-rir:0)
A<yil:'max,11"°"yirzfmax,11’yi-ril:o"°"Vi-'ij:rmax,il
...,ui-rir:0)
........ (3.60)

STEP 6. Continue until a non-singular matrix is encountered, say

A(,...,ui-rijzpj,...). Then set 111 - 'ij-pij + l.

53

STEP 7. Continue steps 5 and 6 until all 7ij are identified. Set Vij

-p +1.

- V max,i

ij - fmax,i

STEP 8. Go to step 2 and repeat for any previously unselected output.

3.5.3 Noise Compensation

The above algorithm must be modified when the observation vectors are
constructed from noisy measurements. This is due to the fact that
because of the random noise, the matrices 0 and A are always near full
rank. This tends to reduce the significant changes in the evolution of

the condition numbers making detection of singularity difficult.

One way to handle noisy observations is to prefilter the data before
applying the identification algorithm. To do this, each data set for a
particular input or output could be passed through a digital filter
whose cut-off frequency is such that it passes only the power in the
true signal. A second order Butterworth type would be adequate. The
cut-off frequency of each of the filters could be determined by
inspecting the power spectrum of each signal. In general, the dynamics
of the process are much slower than the noise of its measurements so

that a cut-off frequency is easily determined.

If the covariance of the noise on each observation is known then the

algorithm can be performed using a compensated A matrix as described in

54

[25]. This method of noise compensation is much simpler than
prefiltering since it just involves subtracting the noise covariance
matrix from the A matrix constructed from the noisy observations. For
an uncorrelated noise structure the covariance matrix reduces to a
diagonal matrix whose diagonal elements are the variance of the noise on
each measurement. An estimate of the noise variances is easily
determined by holding all inputs to the process constant and monitoring

the deviations of the signals from their steady state values.

3.5.4 SVD Subroutines

There exists several reliable programs written to perform singular
value decomposition which may be used in the procedure just described.
The subroutine libraries EISPACK [28] and LINPACK [29] contain fortran
callable subroutines to perform singular value decomposition. The
Linear Algebraic Systems interpretor contains a function which performs
the decomposition. For the research of this thesis an interactive
computer subroutine was written utilizing the EISPACK routine SVD. This
subroutine steps the user through the algorithm just described and
reports the singular values at each step. Only 40 lines of fortran code

were necessary for the subroutine.

55

3.5.5 Simulation Results

As a means of demonstrating the ability of the proposed algorithm
under various conditions, consider again example 3.1. Simulated data
was obtained for this system by driving it with two pseudo random binary
sequences of 127 samples which were kept uncorrelated by delaying one by
half its length. A noise free data set of 500 samples was used to
determine the variance of both the outputs. Additive zero mean Gaussian
noise was added to the outputs. The ratio of the noise standard
deviation to the signal standard deviation was taken to be the noise to
signal ratio. Experimental results were obtained for noise to signal
ratios of 0%, 10% and 20%. The unfiltered data of length 450 samples
was used in all experiments in order to show the degrading effect of
noise and yet the ability of the method to perform in spite of it. The
test used for all cases was the evolution of the smallest singular value

because the magnitude of input and output were comparable.

Results of the identification of the two Kronecker indices for the
polynomial matrix model are illustrated in Figure 3-2. Clearly the
correct values are obtained for the noise free case. In the results for
the case of 10% additive noise it is apparent the significant change in
the evolution of the smallest singular value is very much decreased but

the correct values for ”l and ”2 are still obtained. Again for 20%

noise level, there is a decrease in the significant drOp of the
evolution of the smallest singular value but the value of the indices

can still be identified. Also note that in addition to looking for a

56

 

 

 

 

 

 

 

 

 

 

0% NOISE ORDEHTEST
15401
1E4” 3'
15-01 ‘-
IE-OZ '-
15.00 0
m .... .-
1E-05 ‘-
IE-OB '-
IE'OT ‘-
‘lE-08 0
18-09 : : c L
‘I 2 3 4 5 8
m
10% NOISE ORDERTEST
11
10'
8
8
7
W
smnuwubs '
5
‘
3
2
I
40
30
SHALLEST 2° 0
“VALUE

 

 

1° 0

 

 

 

FIGURE 3-2 MODEL ORDER ESTIMATION WITH
DIFFERENT AMOUNTS OF ADDITIVE
NOISE.

57

significant change in the evolution of the smallest singular value one

can look for the point at which the rate of change is small.

Figures 3-3 thru 3-8 report the value of the smallest singular
value as the observation vectors are removed from the matrix which
identified the Kronecker index for each output. For the noise free case
in Figures 3-3 and 3-6, the values of time delay are very evident due to
the enormous change in the smallest singular value. The same is true

for the determination of p11.

The tests for the 10% noise case are illustrated in Figures 3-4 and
3-7. The change in the singular values is very much reduced but the

true values of time delay and pij are obtained.

The 20% noise level case is shown in Figures 3-5 and 3-8. In this

case the delay of u1 to y2 is not detected due to the lack of a
significant rise in the singular value. The delay from u2 to y1 is also
questionable. However, correct values for pij are still obtained.

These results clearly demonstrate how additive noise prevents the

singularity of the A matrix.

58

1E+01

1E-O1

1 E-03

 

1 E-05

 

1 E-07

 

 

 

 

 

0 1 - 2 3
Magnitude of d (1 = 2)
I max,1

1 E+O1 -

1 E-01

1 503

1 8-05
1 E-07

 

Magnitude of 91 I

Figure 3-3 Input delay test for yl
with 0% noise

59

 

1O

 

 

 

 

 

, l
I

,0 1 2

Magnitude of 91,]

Figure 34 Input delay test folr y
with 10% noise

50

40

30

20

10

40

30

20

10

Figure 3-5 Input delay test for yl

llll_lllll

 

6O

  

 

IIIIJIIII

 

  

 

1 2
Magnitude of 01 J

with 20% noise

 

61

1E+O1 —

1E-O1

1 E-O3

 

1 E-OS

 

1 E-O7

 

 

1 E+01

1 E-O1

1 E-03

1 E-05

 

1 E-07

 

Magnitude of 92,]

Figure 3-6 Input delay test for y2
with 0% noise

62

10

 

 

 

   
    

Magnitude of d I (1

= ‘I )
max,2

 

 

 

Magnitude of 92, I

Figure 3-7 Input delay test for y2
with 10% noise

50

4O

30

20

10

50

4O

30

20

10

63

 

 

 

Magnitude of 92’]

Figure 3-8 Input delay test for y2

with 20% noise

64

3.5.6 Summary of Computer Implementation

We see that the delayed polynomial matrix method that has been
presented represents a unified approach to the structural identification
of multivariable systems even in the presence of relatively high noise.
The identified structure is directly linked to a minimal canonical state
space representation of the system. For systems with time delay in the
inputs, it reduces the number of parameters that need to be estimated
below that of previously reported methods of identification.
Furthermore, the approach is methodical and very little programming is
required to generate a interactive computer aided design tool using this
approach. In the next section we will see how this method compares to
conventional methods of model identification on real data from the Linde

Column.

3.6 Linde Column Model

This section describes the research that produced a dynamic model of
the Linde column. The data used for this work and the constraints it
imposes on the derived model are first described. Then least squares
estimation is explained. This is followed by a description of various
tests of model order. Results of model identification and parameter
estimation are then presented. Finally, the derived model is discussed

with respect to previous knowledge of the column.

65

3.6.1 Operating Data

In order to fully understand the model developed for the Linde
column, it is necessary to review the data used for estimation. For
this research it was necessary to use normal operating records of the
variables that make up the inputs and outputs shown in Figure 3-9. This
allowed for non-interrupted production. When data is collected under
these circumstances, the experimenter must be concerned with the quality
of information available from the data. The question is whether there
is enough transient behavior contained in the data for good model
estimation. This can be especially critical when the inputs to the
system in question are the result of feedback control as is the case
with the Linde Column. As a result, data from the Linde Column was
collected (at a sampling rate of one minute) periodically over the
course of several months. Each data set was reviewed for transient
behavior and only those that showed true excursions from steady state

above the signal noise were used for model estimation.

For instance, Data Set 1 contains a series of pulses in the reflux
flow, figure 3-10. These pulses are not the result induced by feedback
control, but probably occurred because of a sticky valve. The effects
of these pulses in the composition on the 57th tray, Figure 3-11, and in
the bottoms composition, Figure 3-12, are very pronounced. During this
same period the steam pressure on the reboiler remained constant. The

noisy recording of the steam pressure is shown in Figure 3-13.

66

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

7 G11 -
U1 .. . , ,_ . V1
reflu'x ,7 " '3
flow } +57th tray
G ‘ concentration
12 ..
r
\
4-
U2 V2
reboiler steam
pressure ;{ bottoms.
, G22 ,, . concentration

 

 

 

 

 

Figure 3-9 Inputs and outputs of the Linde Column

m

LBS/HR
160
158
«a

120

67

 

lllL

 

 

[11+ 11!!

 

 

 

 

 

 

 

 

    
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

II I I I I ' ' '
mil .:.. «L .I. .3. ..L. ..L. ....
FIGURE 3-10 RECORDED REFLUX FLOW TI"E ("I"-)
(DATA SET 1)
30
% STYRENE -
23
26 I 1
_r ,1 w , ,
E J I
24
22
scene? ' ' I | r u u .
° 3" 400 600 as me 1200 14..
FIGURE 3-11 RECORDED CONCENTRATION ON TIME (MIN.)

THE 57th TRAY (DATA SET 1)

68

 

8.12

% ETHYL-
BENZENE

 

 

 

 

 

 

 

 

 

 

 

“381;; ' I I I I I I I I I

am 400 680 m was 1280 1480
TIME (MIN.)

FIGURE 3-12 RECORDED BOTTOMS CONCENTRATION
(DATA SET 1)

 

PSIG

 

 

 

“1.1...“ ILIIII I

 

 

 

 

 

25.4

 

“I I" || rll

 

 

 

 

 

 

 

 

mg“: l I I I I I I

II a 18 4i 8 608 808 1800 1200 1480
TIME (MIN.)

FIGURE 3-13 RECORDED REBOILER STEAM

PRESSURE (DATA SET 1)

69

In Data Set 2 the changes in steam pressure occurred due to manual
control of the pressure controller set point as depicted in Figure 3-14.
The response of the column to these changes was the classical step
responses in the bottoms concentration shown in Figure 3-15. The reflux
flow during this time period shows a slight drift and several small
pulses, Figure 3-16. The response of the concentration on the 57th tray

' is shown in Figure 3-l7.

Finally in Data Set 3 a series of pulses occurred in the reflux
flow, figure 3-18, while the steam pressure of the reboiler was
constant, Figure 3-19. The response of the concentration on the 57th
tray, Figure 3-20, is similar to that of Figure 3-11. The response at

the bottom of the column is shown in Figure 3-21.

Although the inputs in these data sets was sufficient to produce
responses in the outputs with large deviations from steady state, the
characteristics of the data still put restrictions on the form of the
model that could be estimated. A close inspection of the data will show
that when the column was most excited from steady state only one input
was the cause. Because of this the polynomial matrix representation of
equation (1) could not be estimated directly since it requires both
inputs to be sufficiently exciting at the same time. The logical choice

in this case is a matrix of transfer functions of the form

[VIM ] - [611(2) 612(2) J [111(2) J (3.61)
y2(2) 621(2) 622(2) u2(2)

70

 

26.8
PSIG

 

 

 

Illl

 

 

 

 

 

 

24.5

llllIllll

 

 

 

 

 

 

 

 

5”. I I I I I I

0 200 400 000 000 1000 1200
TIME (MIN.)
FIGURE 3-14 RECORDED REBOILER STEAM
PRESSURE (DATA SET 2)

 

% ETHYL-

BENZENE

0.12 l
0.11

 

 

 

 

 

0.07 “hiflmrj
gigs I I I I r I I

e an 400 can see was 1288
TIME (MIN.)

[Ill lllLIIJIIIIIILLll
I
I

 

 

Illl

 

 

 

 

 

 

 

FIGURE 3-15 RECORDED BOTTOMS
CONCENTRATION (DATA SET 2)

130

LBS/HR
125

120

I05

27

% STYRENE

24-

71

 

llLl

 

 

 

 

 

 

 

 

 

 

Ill]

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I —I I
a see 400 can ' can 1000 1208
TIME (MIN.)
FIGURE 3-16 RECORDED REFLUX FLOW
(DATA SET 2)
J
1
CI
I I r I
8 ass «on can sea 1000 120!
TIME (HIN.)

FIGURE 3-17 RECORDED CONCENTRATION

ON THE 57th TRAY
(DATA SET 2)

72

 

140

LBS/HR
13S

 

Illl IJLL

 

  
  
 
 

 

130

Illl

F'—
=1:—

 

 

 

 

LllJ

' I IN MI! '1 I * ~ I"
. I I ..

I1IITIII llll ITII lllI
303 I

too see 300 m 500
TIME (MIN.)

 

115

.1
d
‘1
q
q
q
d

 

 

 

 

 

 

FIGURE 3-18 RECORDED REFLUX FLOW
(DATA SET 3)

 

 

 

 

   

mm In". ., 1

 

 

 

 

 

 

 

 

 

240‘ It
_ I
£413 I I I I I I I 1 I I I I I I I I I I I I
I0. 200 3“ 400 5“

TIME (MIN.)
FIGURE 3-19 RECORDED REBOILER
STEAM PRESSURE
(DATA SET 3)

73

a?
% STYRENE

 

llll

 

 

 

Ll_ll

 

 

 

 

 

 

 

 

a
lllJ
H

—

 

 

 

 

III
A
III
fEIIIEE
‘____
~lIE:;,

 

 

 

 

 

 

 

3alzafiIII III—I IIII IIII IIII
L 1“ a!” 300 4“ 5“
TIME (MIN.)
FIGURE 3-20 RECORDED CONCENTRATION
ON THE 57th TRAY
(DATA SET 3)

 

0.105

% ETHvLé
BENZENE
0.1”

[III

 

 

=E;l‘.::=
b:
J

. A n‘
...,, I In HI
I“ .4

0 IIIII IIIFIIIIIIIIIIII
{£300 1” 200 ace 4“ 5”

TIME (MIN.)

 

 

 

-r“l

 

__j:=

 

 

 

 

 

 

 

 

 

 

 

 

 

FIGURE 3-21 RECORDED BOTTOMS
CONCENTRATION (DATA SET 3)

74

Here u1(k) is the deviation from steady state of the reflux flow into
the column and u2(k) is the deviation from steady state of the steam
pressure on the reboiler. Likewise, y1(k) is the deviation from
steady state of the concentration on the 57th tray and y2(k) is the

deviation from steady state of the concentration in the bottoms product.

In general each individual transfer function has the form

n
z" X b 2.0
a
a-l

-a
1 + 2:1 aaz

G(z) - (3.62)

 

Another restrictive characteristic of the data is the longer
sampling rate of the bottoms concentration. The 8 minute sample and
hold function is apparent in the recordings of the bottoms

concentration. This led to G21(z) and 622(2) having sample rates of 8

minutes while 611(2) and G12(z) had sample rates of 1 minute.

A final characteristic of the data that restricted the model was
the type of inputs. The type of input excitation developed by the
reboiler steam pressure was that of a step change. With such an input

transfer functions G12(z) and 622(2) must be modified as

75

-r-1
boz

G(z) - (3.63)

n
1 + 21 aazm‘

where bo - b1 + b1 +...+bn. To see this, consider (3.62) rewritten as a

difference equation

11

n
y(z) - -Z aay(k-a) + Z bu(k-a-r) (3.64)
1 a-1 a

a-
Now in the case of a step input u(k) - c for k > 0. So at any time k,

n n

y(z) - -Z aay(k-a) + Z bac (3.65)
l a-l

a-

Obviously in (3.65)the parameters b1,...bn appear as a sum and

separation of them is impossible for a constant input.
3.6.2 Least Squares Parameter Estimation

Having established the form of the model the parameter estimation

problem may now be stated. The generic transfer function of (3.62) may

be rewritten as

 

76

y(k) - -a1y(k-l) - ... - any(k-n) + b1(k-1-r) + ...

. + bnu(k-n-r) + 6(k) (3.66)

where 6(k) accounts for model error and measurement noise. This can be

rewritten more concisely as

y(k) - ¢T<k)o + e<k> (3.67)

b .. anT and 0 - [y(k-l). y(k-z).

where 9 - [~a1, ..., -an, ..., 1,

., y(k-n),u(k-1-r), u(k-n-r)]T.

Now given the model of (3.67) determine the parameters a1 and bi that

minimize the least squares criterion

N
J (o) - Z 62(k;0) (3.68)
n
k-n
over the data {y(k), u(k), k - l, 2, ..., N}.

A

The least squares estimate of the parameter vector, 0, may be found
by processing the data in batch form or by proceeding through the data
sequentially. The least squares solution for processing the data as a

whole has the form

A

a - [0T01'10TY (3.69)

77

where

e - [¢<n>. ¢(n+1>. .... ¢(N)]T

Y - [y(n). y(n+1>. .... y(N>]T (3.70)

For large data sets it may be simpler to process the data a sample
at a time thus avoiding having to deal with a large O matrix. The
recursive least squares algorithm performs such a task. The least
squares estimate of 0 may be computed iteratively from the following

equations

9(k> - o<k-1) + x<k>ty<k> - ¢T(k)0(k-1)] (3.71)
where [(k) is the time varying gain matrix computed as

P(k-1)¢(k)
K(k) - T (3.72)
l + ¢ P(k-l)¢(k)

 

Matrix P(k) is computed recursively from

P(k-1>¢(k>¢T<k>r<k-1>
P(k) - P(k-l) -

 

T (3.73)
1 + ¢ (k)P(k-1)¢(k)

78

The recursive equations (3.71) through (3.73) require starting

values for P(O) and 0(0) . For large data sets the values can simply be

P(O) - pI, where u is large, and 0(0) - [0,...,0]T.

3.6.3 Transfer Function Identification

The model order n and the time delay 7 of each of the four transfer
functions in (3.61) were determined by a combination of methods. Where
appropriate, the delayed polynomial matrix method presented in the
previous section was. applied for the single input single output case.

Next, models of increasing order and varying time delay were estimated
and the mean square residuals, Jn(0), of each was computed. The value

of the mean square residualsshould decrease as the order of the models
is increased toward n + r. Ideally, when the correct time delay is
chosen and the order increased above the true order, the magnitude of
the mean square residuals should remain relatively constant. Likewise

when the value of the time delay is changed from the true value there

A

should be a corresponding increase in Jn(0).

A statistical method was used to see if the mean square residuals

changed significantly for an increase in model order, say from n - n1 to

n - n2. This so called F-test assumes that the residuals are Guassian

and that Jn(0) and Jn(0) - Jn(0) are statistically independent and are
2 l 2

Chi- square distributed with (N - 2n2) and (2n2 - 2n1) degrees of

79

freedom respectively. To test the null hypothesis that n1 is the true

model order, the test statistic

 

Jn1 - an N - 2n2
t - (3.74)
Jn 2(n2 - n1)
2
was computed. The statistic t has an F{(2n2- 2n1),(N - 2n2)

distribution under the null hypothesis. After a risk level was defined
*
a corresponding t was taken from a table of the F distribution. The

null hypothesis was rejected for t > t*.

As a final check the sample auto-correlation function of the

residuals

n+N
Re(r) - _1_ X: e(k)e(k+r) (3.75)
N k-n+l
r - O, 1, 2,

was computed for the selected model. Assuming the errors in (3.67) are
white, and the correct model is chosen, the residuals should also be
white. When this is the case the values of the lags in the sample auto-

correlation should all be nearly zero except for the zero lag (r - 0).

80

3.6.4 Results of Column Modelling

It was found that nearly all the model order tests for transfer

function G11(z) agreed. The results of the delayed polynomial matrix

model order test described in the previous section are shown in Figure
3—22. Part a of Figure 3-22 shows the Kronecker index to most likely be
9. Part b indicates that the time delay to be 5. These estimates are
supported by the magnitude of of the mean square residuals shown in
Figure 3-23. It is clear that for a time delay of 5 there is little

reduction in the model error for order larger than 4 but a large
increase in error for order less than 4. Likewise the value of Jn(0)

for order 4 and varying delays were:
delay 4 17.2
delay 5 14.0
delay 6 16.0

indicating that a time delay of 5 is the most appropriate.

*
Using the F-test with a 10% risk factor the values of t are:

n2 - III 12*
1 2.30
2 1.94
3 1.77
4 1.67

For n1 - 6 with a delay of 5 the following values of t resulted:

81

G(1 1) ORDER TEST

SING. VALUE

 

 

 

0.0 : : 2 E : % E : 4
2 3 4 5 8 7 8 9 10 11
cm
G(1 1) TIME DELAY TEST

SIHG. VALUE

 

 

 

0 i 2 3 4 s 6 7
most."

FIGURE 3-22 RESULTS OF THE DELAYED POLYNOMIAL
MATRIX TEST FOR 611(2).

82

G(11) M. S. RESIDUALS

IIEANSGUARE 3° ..
RESIDUALS

 

 

 

 

10 L

FIGURE 3-23 MEAN SQUARE RESIDUALS
FOR G11 WITH A TIME

DELAY OF 5.

 

LL11

 

[Ill

 

l_l

 

0.01

 

 

 

 

 

 

.. AAA ....
'I/vv U \/\

‘230‘ I I I r I I I I I I I I I Ti I I I

 

 

 

 

 

FIGURE 3-24 AUTO CORRELATION OF RESIDUALS
OF 4TH ORDER MODEL WITH A
DELAY OF 5 FOR G11.

83

n2 t

3 105.91
4 11.39
5 1.20

The F-test indicates that for the given risk level a model order of 5 is
correct.

However the residuals of a fourth order model with time delay of 5
appear to be white as demonstrated by the sample auto correlation

function shown in Figure 3-24.

This model structure was accepted and parameter estimation was
performed using recursive least squares on Data Set 1. The estimated

model is

l 2 3

(.0332' - .022' + .00242' - .0512’4)z'5

6(2)-
11 1 - .832'1 + .392”2 - .97z'3

 

+ .482-4

A portion of the output of this model driven by the recorded reflux flow
is shown in Figure 3-25 and can be compared to the true column response

in Figure 3-26.

For transfer function 621(2) Data Set 1 was used. The model order

tests were found to be in agreement but some were more definite than
others. Applying the delayed polynomial matrix method of the previous

section to determine the Kronecker index produced Figure 3-27. The

% STYRENE

% STYRENE

84

 

 

 

 

 

 

 

 

 

 

 

 

._' [W «I
- H
.. u

TIME (MIN.)

FIGURE 3-25 RESPONSE OF ESTIMATED G11

TO THE RECORDED REFLUX
FLOW OF DATA SET 1.

 

 

 

 

 

 

 

.6 I'll IIVI ITII erflIl,TUllI

500 800 700 800 300 1000 1100
TIME (MIN.)

FIGURE 3-26 A PORTION OF THE RECORDED
57TH TRAY CONCENTRATION
FROM DATA SET 1.

 

 

 

 

 

 

 

SIIALLEST
SING. VALUE

0.00055 «-

0.00050 ‘-

0.00045 --

0.00040 ‘-

0.00035 '-

0.00030 --

 

0.00025

85

G(21) ORDER TEST

 

 

FIGURE 3-27 ORDER TEST FOR 612.

MEANSOUARE

0.00030

0.00030

0.00034 '

0.00032 ‘-
0.00030 --

0.00023 ‘-

0.00026 -

 

0.00024

6(21) M. s. RESIDUALS

 

FIGURE 3-28 MEAN SQUARE RESIDUAS'

FOR G12 WITH A TIME
DELAY OF ZERO.

86

evolution of singular values indicates that a third or fourth order
model could be appropriate. The mean square residuals, Figure 3-28,
clearly show the model order to be 4. The estimated 4th order model

with no time delay was

1 2 3 4

.0000332' + .000162' + .000222' + .000152-

G (2) -
21 1 - .582'1 - .362'2 + .IIz'3 - .0132'4

 

The relatively low magnitude of a4 and b1 in the model suggests that

they might be set identically to zero resulting in a third order model

with time delay of 1. This alternative model was estimated. For these

competing models t - 0.279 and t* - 2.39 for a 10% risk factor. This

suggests that a model with order 3 and a time delay of 1 is acceptable.

A model order of 3 with time delay 1 was accepted and least squares

estimation lead to the following model

1 2

(.000162' + .000222' + .00152‘3)z'1

21 1 - .592’1 - .352.“2 + 0.942’3

 

The auto correlation of the residuals of this model are shown in Figure
3-29 and appear to be white. A portion of the output of this model
driven by the recorded reflux flow is shown in Figure 3-30 and can be

compared to the true column response in Figure 3-31.

0.2

éfika

87

 

 

 

 

 

 

 

 

 

\//\ V VA\/
I I I I I I I I fl I I I I I I I
0 10 15
FIGURE 3-29 AUTO CORRELATION OF RESIDUALS.

OF A 3rd ORDER MODEL WITH A
DELAY OFil FOR 621(2).

 

% ETHYL-
BENZENE

x 10-2

% ETHYL-
BENZENE

0.4

FICURE 3-30

0.02
0.01

IV“

88

 

 

 

 

 

 

 

 

 

 

 

Gz

 

 

 

 

. /\ A‘
‘IH AIM
-VJ\ /

I \M

- L/

 

1200

TIME (MIN.)

RESPONSE’OF THE ESTIMATED

TO A PORTION OF THE

RECORDED REFLUX FLOW OR
DATA SET 1.

 

 

l_lll

 

[III

A

 

IIIJ

 

 

 

 

IIII

 

IIII

 

FIGURE 3-31

TWII

 

 

 

 

IIII

1000

III]

1100

 

TIII

TIME (MIN.)

A' FORTION OF THE RECORDED
BOTTOMS CONCENTRATION FROM
DATA SET 1.

89

In determining transfer function 022(2) the delayed polynomial

matrix method of the previous section was not used because the shape of
the response of the bottoms concentration to a step change in steam
pressure strongly suggests that at most a second order denominator with
perhaps at most one period of time delay would adequately fit this

response. The mean square residuals for these competing models were:

lst order, 0 delay .000076
2nd order, 0 delay .000072
lst order, 1 delay .000078
2nd order, 1 delay .000074

It is clear that for a delay of l the model error is increased. For

comparing the second order model with no delay to the first order model

with no delay the F-test statistic t - 2.61 which is lower than the t* -

2.83 for this comparison. Therefore G22(z) was taken to be lst order

with no delay. The resulting model derived by least squares estimation

is

-.0026z'1
G (2) -
22 1 - .942'1

 

A portion of this model driven by the recorded steam pressure is
shown in Figure 3-32 and can be compared to the true column response in

Figure 3-33.

90

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mat
% Ethyl- ‘
Benzene : /
. Him
aJm“ V' . \\
-001- .\\
-a.ea \
A \\\\~\\
“0.03 rIIr Ile IIII IIIIIIII IIII IIII
«no 500 see no see see moo nee
Time (Min.)
Figure 3-32 Response of estimated C22
to a portion of the recorded
steam flow from Data set 2.
one
: (I
_, V
% Ethyl- _
Benzene 0.09

 

L
' WWW '1.

.. WM

 

 

 

 

 

 

 

 

 

 

 

 

POP IIIFIIII IIIFfIIrIIIITIII rIII
400 500 600 700 800 900 1000 11"

Time (Min.)
Figure 3-33 A portion of the recorded
bottoms concentration from
Data set 2.

91

The construction of transfer function 612(2) proved to be the most

difficult due to the data available for model estimation. Data Set 2
was used because of the step changes in steam pressure. However, the
drift in the reflux flow could not be discounted because of its strong

influence, as evidenced by 611(2). This influence was removed before
model estimation by subtracting the output of G11(z) driven by the

recorded reflux flow from the record of the 57th tray concentration.
The model order tests did not give clear results but it was found that
the mean square residuals were minimized for a model order of 3 with
time delay of 5. Recursive least squares estimation resulted in the

following model

(.112’1)z'5

G (2) -
12 1 - .aaz'l + .392'2 - .872'3

 

3.6.5 Discussion of the Linde Column Model

Several observations may now be stated in regard to the derived
model of the Linde Column. The most surprising is the non-minimum phase
behavior of the concentration on the 57th tray to a positive change in
reflux flow as seen in the inverse response in Data Sets 1 and 3. This

behavior is present in 611(2) in the form of the zeroes outside the unit

circle. This is the first time such behavior has been reported in the

literature. The control implications of this discovery have been

92

outlined in the introduction and will be discussed in detail in the next

chapter.

Data Set 3 was also used to estimate 611(2). Model order tests were

fairly consistent with those applied to Data Set 1. Parameter
estimation lead to a model which showed a large shift in some of the

poles and zeros of transfer function 611(2). The poles and zeros

obtained from both data sets are

POLES ZEROS

Data Set 1 0.59 1.37
0.93 -0.38 + j.98
0.34 + j0.87

Data Set 3 0.12 1.01
0.94 -0.55 + j0.81
0.01 + j0.78

The reason for these inconsistencies is the fact that Data Set 3 does
not have an obvious steady state value for the recording of the 57th
tray concentration. Only controlled experiments on the column could
resolve this issue. However, it will be seen in the next chapter that a
model based controller can be developed which explicitly handles process

/ model mismatch.

93

The first order model obtained for 622(2) is consistent with the

perception the production operators have concerning the dynamics in this
part of the column. They had learned that if the bottoms concentration
was too high a large increase in steam pressure would bring it back

without oscillations. The two hour time constant obtained for 62

2(2) is also consistent with the ability of the operators to maintain

bottoms concentration specification using manual control on the reboiler
steam pressure controller while taking manual samples of the bottom
product every 4 hours. Only a system with such a large time constant

would make such a manual control scheme possible.

The existence of the coupling terms G12(z) and 621(2) confirms the

suspicion of process engineers associated with the Linde Column. These
engineers had experienced the effects due to such coupling during start-
ups and abnormal operations of the column. The level of coupling

represented by 612(2) and 621(2) is also consistent with the ability of

two independent PI controllers to maintain control of the column. This

will be demonstrated in the next chapter.

3.7 Summary of Chapter Three

This chapter has dealt with the development of a dynamic model of

the Linde Column. The research presented in this chapter led to two

significant contributions to the literature. A new identification

94

procedure which includes time delay explicitly was developed through
adaptation of current theory on determining the structure and order of a
multi-input, multi-output model. In addition to this theoretical
contribution, this Chapter also describes a form of non-minimum phase
behavior in the column which has not been mentioned before in the

literature.

The new identification procedure is significant because it greatly
reduces the number of parameters necessary to describe multivariable
systems with multiple time delays. The reduction in parameters is
realized through a new representation of multi-input, multi-output
systems called the "delayed polynomial matrix” representation. This new
representation explicitly accounts for multiple delays. The procedure
first estimates the order of a polynomial matrix model of the system
which serves as an upper bound on the number of parameters necessary to
describe the system behavior. It then derives the delayed polynomial
matrix model in a methodical and straightforward manner. An integral
part of the procedure is the use of singular value decomposition as a
means of testing for singular matrices. The procedure is well suited

for interactive computer use.

CHAPTER FOUR

ADVANCED CONTROL FOR THE LINDE COLUMN

The purpose of this chapter is to investigate the implications that
the model of the Linde Column has on control strategy. Of course
practical considerations will weigh heavily in the investigation since
the Linde Column is a production scale process. No attempt is made to
find the best control strategy out of the set of all possible control
strategies. Instead, the work described here is focused on a particular
alternative model based control strategy suggested by the constraints

imposed by the model and practical concerns.

4.1 Current Control Strategy

As described in Chapter Two, the current control strategy of the
Linde Column is two independent PI controllers, one controlling bottoms
concentration through manipulation of reboiler steam pressure and the
other controlling the concentration on the 57th tray through
manipulation of the reflux flow. This present -scheme has been
successful. The column operates in a relatively stable manner and
produces a product which meets specifications. However, several dynamic

characteristics brought to light in the previous chapter have

95

96

detrimental effects on such a control strategy and explain why

production engineers have never been truly satisfied with it.

Probably the least interesting characteristics, but
certainly one that is seen often in industrial processes, is the
apparent sticky valve manipulating the reflux flow. Obviously a sticky
valve would have adverse affects on any control strategy. Currently the
reflux flow valve does not have a positioner. A positioner should be

installed even if nothing else is changed in the control strategy.

A characteristic of more theoretical interest is the interaction

between the top and bottom of the column described by 612(2) and
621(2). This interaction reduces the achievable performance for the

independent controllers of the current control scheme. The controller
response to process disturbances must be made sluggish enough to avoid
having the control action in one loop cause an appreciable disturbance
in the other loop. A multivariable control strategy could address this

problem.

The non-minimum phase behavior at the top of the column noted in

Chapter Three via G11(z) is troublesome for the current PI controller

used there. A PI controller is not well suited for either time delay or
inverse response. A control strategy which takes explicit account of
the time delay and inverse response could improve the control of the

concentration at the top of the column.

97

In spite of the disadvantages just outlined, there exists strong
practical reasons to use the current control strategy. First of all PI
control is well understood by operating personnel. This is important
because a control strategy only has a chance to succeed if it is
understood and therefore accepted by the people who must live with it on
a daily basis. Another advantage PI control has is its smooth
operation. It rarely causes abrupt changes in the process input unlike
many high performance control algorithms. Probably the most important
reason to use PI control is its robustness. This simple control
strategy can usually be implemented with little knowledge of the process
model. The parameters of the algorithm can be determined experimentally
online. Then if the process dynamics change enough to cause a
noticeable affect in the controller performance, operating personnel
know enough to ”detune' the parameter values. Certainly any alternative
to the present strategy would have to have these favorable

characteristics.

Alternative strategies that would possess the advantages of PI
control would be those which maintain the PI controller as an element of
the overall scheme. The essence of these approaches is to compensate
the PI controller. These schemes remain familiar to the operating
personnel while they address problems such as time delay. The Smith
Predictor is an example. While these strategies have been reported to
‘work well for the control of distillation columns, a novel approach will

be taken in this thesis for reasons that will become clear.

98

4.2 Literature Survey on Distillation Control

It is useful to review past research on distillation column control
before describing the control strategy developed for the Linde Column.

This will put the contributions of this thesis in perspective.

Simultaneous control of both product concentrations of a
distillation column has received considerable attention in recent years
because of the importance of the process. Many different multivariable
strategies have been tried either through simulation or application to
an operating process. However all these investigations were concerned

with the dynamics of pilot scale columns.

For instance, Wood and Berry [8] implemented two different
multivariable schemes on an eight tray 9 inch diameter distillation
column. They confirmed experimentally what earlier researchers had
found through simulations, ratio control (Rijnsdorp [9]) and decoupling
compensators (Luyben [10]) were an improvement over the conventional
control described in Chapter Two. They conclude that application to

industrial columns would show even more improvement.

Dahlqvist [11] investigated the use of self-tuning regulators for
control of top and bottoms product compositions. He was more concerned
with the details necessary to implement the scheme than the improvement

such a scheme might have over conventional control.

99

Ogunnake and Ray [12] used the model from Wood and Berry to simulate
on a digital computer the control obtained from their multi-time delay
compensator. The multi-time delay compensator was used with and without
decoupling Compensation. In both cases the multivariable approach was
an improvement over conventional control and they concluded that even
better results could be obtained with industrial columns where time

delays would even be greater.

Ogunnaike et al [13] then applied Ogunnaike and Ray's multi-delay
compensator to a 19 plate 12 inch diameter experimental distillation
column. They tested disturbance rejection as well as set point
following. They found that the multi-delay compensator out performed

conventional control for most tests.

Weischedel and McAvoy [14] studied the application of three
different decoupling schemes to three different columns of 13, 17 and 19
trays. Digital simulation, was used to test the various control
strategies. They concluded that for moderate to high purity product
separations decoupling schemes failed to fully decouple the overhead and
bottom controllers but the multivariable approaches were a definite

improvement over conventional control.

Garcia and Morari [30] applied Internal Model Control to a
simulation of the pilot scale column first introduced by Wood and Berry.
They too found that a multivariable controller out performed
conventional control in simultaneously controlling overhead and bottoms

product compositions. Although robust control is one of the attributes

100

of Internal Model Control Garcia And Morari failed to test this in their

simulations.

The research presented in this thesis extends the work of past
research on multivariable control of distillation columns by applying
one of the multivariable strategies to a model of a full scale column.
Several new and useful modifications were made to Internal Model Control
to tailor it for the Linde column. The control studies of the Linde
column presented in this thesis provide insight into the transition of

advanced control to full scale columns.

4.3 Model Based Control, 5150

Internal Model Control (IMC) is one of the few modern, model based
control strategies that retains the advantages of PI control and is
capable of overcoming its deficiencies. For this reason it was selected
as the alternate control strategy for the control studies of this
chapter. The basic structure of IMC that was used in this research is
shown in Figure 4-1. The structure of IMC and its relation to other
control schemes was developed by Garcia and Morari [2]. Their work was
based on Brosilow's [31] earlier work on inferential control. It will
be shown that this simple structure addresses performance and robustness

in a very straight forward manner.

101

 

 

controller process

E 7*

 

 

 

9.1

Internal model

 

 

Figure 4-1 Basic Internal Model control structure

 

 

 

+ ,
s—DQ—o P(z)

fllter controller process

 

 

 

 

internal model

Figure 4—2 Filtered Internal Model control stucture

102

4.3.1 Minimum Phase Processes

The similarity of IMC to conventional feedback control is evident in
Figure 4-1. This feature makes it understandable to production
personnel. In fact the IMC structure is analogous to the way an

experienced operator might control a process. In IMC the control action

A

, m, is applied to a model of the process, G(z), as well as the process
itself, G(z). A prediction of the process output is made by the model
and that is compared to the measured process output. In the same way an
operator uses his mental model of the process based on his experience.
By comparing instrument readings with his expected process response he
can identify faulty instruments or large process disturbances and
respond accordingly. IMC is actually more intuitively appealing than

conventional feedback control.

The biggest reason IMC overcomes the deficiencies of PI control is

the special form of its feedback signal 3(2)
3(2) - (1 + (G(z) - G(z))cc<z>51 d<z> (4.1)
Note that for a perfect model,

5(2) - G(z). (4.2)

103
the feedback signal is the disturbance d(z). This gives automatic time
delay compensation. Also note that differences in the model and process

are contained in 3(2). As we shall see this information can be used in

a straight foward manner to address robustness.

The other desirable features of IMC can be illustrated by referring

to the closed loop equations

Gc(2)
(3(2) - d(2)) (4.3)

 

m(2) - .
1 + (Gc(2)(G(Z) - 6(2))

G (z)G(z)
° (s<z> - d<z)) + d<z> (4.4)

 

Y(Z) - K
1 + GC(Z)(G(Z) ' G(z))

It is easy to show that for a controller Gc(z) which satisfies

Gc(l) - 1 / 0(1) (4.5)

there is automatic integral action leading to zero steady state offset.

Since for a constant set point s(z) - s*

“ -1
6(1) 6(1) (9* - d) + d - 3*

 

y A A
” 1+cafhmn-ca»

104

for constant disturbance d. Notice that steady state error is zero even

for an imperfect model.

From equation (4.4) it is evident that for a stable minimum phase

process perfect control is possible. Let

N(z) z-

G(z) ———D(z)

Where N(z) and D(z) are polynomials with roots inside the unit circle.

Then if Gc(z) - 0-1(2) and 8(2) - G(z)

D(z) N(z) z-l
N(z) D(z)
y(Z) ' D(z) .
1 + N—(z)—(G(z) - G(z))

 

(5(2) - d(2)) + d(2)

 

1 + (1 - z'l)d(z)

- s(z)z-
Notice that the single period of delay is unavoidable since it is
inherent in all physical discrete time systems. The case of non-minimum

phase systems will be discussed in the next section

Robustness, the ability to control in the event of parameter changes
in the process, is addressed by the insertion of a filter in the IMC
structure, Figure 4-2. The addition of a filter changes equations (4.3)

and (4.4) to

105

 

F(z)Gc(z)
m(z) - A (5(2) ' d(Z)) (4.6)
1 + F(2)Gc(2)(c(z) - G(z))
F(z)Gc(z)G(z) ‘
y(z) - (8(2) - d(z)) + d(z) (4.7)

 

1 + F(2)Gc(2)(G(Z) - 0(2))

The filter P(z) can be chosen to ensure stable input and output
sequences for a given process-model mismatch. So as the model deviates
farther from the true nature of the process, the filter can be modified

to slow the controller down to maintain stability.

The filter can also be used to control the closed loop response of

the system for a well matched internal model, G(z) - G(z), since
“2’ - P(z)
8(2)

Now if the filter is first order its time constant can be used by
operators as a single online tuning parameter to shape the system
response or to adjust to changing process dynamics. This is an
improvement over PI control since only one parameter needs to be

adjusted rather than two.

106

Other desirable features, such as, control effort limiting and
bumpless transfer between manual and automatic, have been described by
Parrish and Brosilow [34]. These features together with the ones
detailed here make IMC a very appealing control scheme from an
industrial point of view. IMC was used in the control studies of this
research since it retains the advantages of PI control (zero steady
state error, robustness, easy tuning, and conceptually simple) and also
offers improved control, especially for processes with time delay or

inverse response.

4.3.2 Non-minimum Phase Processes

We have seen in the previous section that for minimum phase
processes (stable with no time delay or zeroes outside the unit circle)
IMC can produce perfect control based on a controller designed with
equations (4.2) and (4.5). When dealing with a non-minimum phase
process equation (4.5) cannot be used directly since G(z) is not
invertible. In this section several techniques are described to replace

equation (4.5)

4.3.2.1 Factorization

One approach to handle a non-invertible process model is to apply

the following factorization

107

G(z) - G_(z)G+(z) (4.8)

where G+(z) contains all the zeroes outside the unit circle and all the
time delays. G-(z) is then used in equation (4.5) in place of G(z).

The factorization (4.8) is not unique. Garcia and Morari [2] offer the

following form

2 - V1 1 - u

~r-l
G+(Z) - Z A (4.9)

i-l z - u l - v

 

where r is the time delay of the process and "i are the p zeroes of G(z)

A

and ”i are the images inside the unit circle, i.e.

A

”1

Vi, [vi] 5 l

A

vi - l / v1, luil > 1

Figure 4-19 illustrates this procedure .

The use of the exact model inverse or its factored form of (4.8) for

the design of the controller can lead to excessive control action for

108

Typical unstable
zero with magnitude
greater than 1
O
4 Hell}

0 - Image of the
unstable zero

  

 

Figure 4-19 Stable images of unstable zeros

109

minimum or non-minimum phase processes. This is demonstrated for non-

 

minimum phase processes by 011(2) of the Linde Column model,
3 2 -5
(0.0332 - 0.01992 - 0.002382 - 0.0507) 2
G (2) -
11 z“ - 0.826723 + 0.388222 - 0.96642 + 0.4814

Using (4.8) and (4.9) to factor G11(z) produces

2'5(2 - 1.3698)(2 - 0.3831iJO.9874)

(2 - 0.7300)(2 - 0.3417ijo.8808) x

 

G+(Z) -

(l + 0.7300)(l + 0.6835 + 0.8927)
(1 - 1.3698)(1 + 0.7662 + 1.1218)

 

(z - 0.7300)(2 - 0.3417i30.8808)

4 3 2 x
2 - 0.82672 + 0.38832 - 0.96652 + 0.4143

G_(2) -

 

(1 - 1.3698)(1 + 0.7662 + 1.1218)
(1 - 0.7300)(1 + 0.6835 + 0.8972)

 

The resulting controller is

24 - 0.826723 + 0.388322 - 0.96652 + 0.4143

Gc(z) - 3 2
(z - 0.04652 + 0.39372 - 0.6517)(-0.0506)

 

Figure 4-3a shows the results of a simulation in which G11(z) was

considered as a single input, single output system controlled by an IMC

with the filter time constant set equal to zero and a control block

-1

2.5

2.0

1.0

0.5

$05

110

 

llLL

 

 

L111

 

LL11

 

 

1 LLJ

 

 

11L_L

 

 

IIII IIII rrIT fiIII

0 1 00 200 300 . 400

FIGURE 4-3a RESPONSE OF 611(2) TO A STEP

LOAD DISTURBANCE, IMC CONTROL
WITH NO FILTERING.

 

 

 

 

 

 

1114 _Lll1
”‘21

 

 

[J 11
M

 

 

JLJLJLIJI
M

 

 

1J 1.1

 

I I I I I I I I I I I I I I I I

0 1% 200 300 400
FIGURE 4-3b RESPONSE OF 611(2) TO A STEP

LOAD DISTURBANCE, IMC CONTROL
WITH FILTER TIME CONSTANT = 0.85

 

 

 

 

 

lll

given as Gc(z). In this simulation a step load disturbance was

introduced at time sample 50. We see that disturbance is compensated
for quickly but that there is a large initial spike and oscillations
around the set point. Figure 4-3b shows much smoother operation with a
large filter time constant of 0.85. As a comparison Figure 4-4 shows

the response of the system with the control block set to a constant,

A

i.e. Gc(2) - l / G(l), and the filter time constant set equal to

zero. Notice the smooth operation without the use of the tuning filter.

4.3.2.2 Predictive Controller

Figure 4-4 suggests that a way to avoid a hyperactive controller is
to use an approximate inverse of the process other than the one
prescribed by (4.9). One possibility is the method presented in [2] in
which an approximate inverse is found by solving the following

predictive control problem

Consider the impulse response representation of the process

0
0(2) _ z-TZ giz-i-l-l
i-l

At each discrete time k find the solution to

2.5

2.0

1.5

1.0

0.5

112

 

 

J111 Alli
=3

 

1111
PM

 

L111
fl

 

L111

 

 

1111

 

 

FIGURE 4-4

 

 

 

 

100' 200 300 400

RESPONSE OF 011(2) TO A STEP

LOAD DISTURBANCE, IMC CONTROL
WITH (Gc(z) = l/Gll(1) AND NO

FILTERING.

113

P
min 2: 7:[yd(k+r+n) - y(k+r+n|k)]2 + B:m(k+n-1)2
n-l

over the set {m(i): i-k, k+1, ..., k+M-l}, subject to the constraints

y(k+r|k) - ym(k+r) + d(k+r|k)

- g1m(k-l) + g2m(k-2) + ... + gNm(k-N) + d(k+r|k)
m(k+M-1) - m(k+M) - - m(k+P-1)
fig - 0, n > M

Here yd(k+r+n) is the future set point, P is the prediction horizon, 7:

are the time varying weights on the output error, 3: are the time

varying weights on the input, M is the input suppression parameter which
specifies the number of intervals into the future during which m(k) is
allowed to vary (m(k) remaining constant afterwards), r is the system

time delay, ym(k+r) is the output of the internal model, y(k+r|k) is the

predicted output, and d(k+r|k) is the predicted disturbance.

This method allows the control engineer considerable flexibility in
designing the closed loop response. However, the solution depends on

several tuning parameters, P, M, 7 fin, and the solution is

n!

accomplished by a least squares approximating prdblem. In addition, if

114

on-line tuning is desired a matrix equation must be solved each time a
tuning parameter is Changed. To avoid the complexity of the predictive
control problem, a simpler procedure was developed and used in this

research.

4.3.2.3 Reduced Order Controller

The method that follows is based on formulating the controller from
an inverse of a reduced order model of the process. This approach is
motivated from three observations. (1) The abrupt behavior demonstrated
in Figure 4-3a is a manifestation of high frequency components in the
control system. Therefore a reduced order model approximating the low
frequency dynamics of the process can reduce this high frequency
behavior. (2) The real benefit of the IMC structure lies not in the
controller but in the internal model. This means that for invertible
process models (minimum phase) little performance is sacrificed when
lower order models are used for controller design. In fact satisfactory
performance is usually obtained even with a proportional controller
(i.e. a constant gain equal to the inverse of the model steady state
gain) as shown in Figure 4-4. For non-invertible process models (non-
minimum phase) the loss of performance can be even less. (3) A reduced
order controller would decrease the complexity of the control program.
Although it is true that the calculations of the program are usually
transparent to the production operators, they are not charged with
maintaining the program or updating it for changes in the process. The

production engineers who would be responsible for program maintenance

115

tend to relinquish that responsibility after a rather short time. So it

is important to keep the program as simple as possible so that it is

understandable for the new engineer.

The concept of reduced order models for controller design is not
new. The problem of model order reduction has been addressed in the
literature, see Jamshidi [33] for an overview and a good list of
references. There are two basic approaches for reducing transfer
functions. One deals with manipulations of the transfer function using
methods such as continued fraction expansion or Pade' approximation.
The other is based the minimization of the difference in the frequency

or time responses of the reduced and full models.

The approach taken for G11(2) in this research was of the second
type. A reduced order model, Gr(z), with the same time delay as 611(z)

was sought which minimized the difference between its time response and

that of 011(2). Since there was a level of uncertainty associated with
the coefficients of G11(2) due to the noisy data used for estimation, it

was natural to use as reduced order models those that were obtained in
Chapter Three. This is in contrast to the method proposed by Sinha and

colleagues [35] and [36] in which a model is found which approximates

the step response of the full model.

Models of four different orders exist for a reduced order model of

011(2): 3rd, 2nd, lst, and zero or constant. The criteria used to

116

select the model for controller design was the simulated time response
of the closed loop system to the same step load Change that was used to
generate Figure 4-3. It was found that a first order model resulted in
a controller which produced a quick return to set point without
overshoot. The 2nd and 3rd order models produced overshoot. The first
order model had the added advantage that it did not have a zero outside
the unit circle so that it could be inverted directly. The
factorization of (4.8) had to be invoked for both high order models to

produce stable Controllers. For these reasons the following first order

controller was selected

1 - 0.92332‘1
60(2) ' ' 0.3532 (4 1°)

 

For this controller the numerator was adjusted slightly so that its

steady state gain was the reciprocal of 611(1).

A filter time constant of 0.5 was determined experimentally to
produce an acceptable level of controller action. Figure 4-6 shows the
response of the system to the same disturbance as Figure 4-3 and 4-4.
Notice how the reduced order controller is able to nearly match that of
the full order controller, Figure 4-3. Notice also the improvement the
selected controller shows over the constant controller, Figure 4-4, and

a well tuned PI controller, Figure 4-5.

1.5

1.0

1.5

‘ 1.0

117

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

- ll

3 l

j \

~ by

3

q l I I l l

0 100 400
FIGURE 4-5 RESPONSE OF 611(2) TO A STEP LOAD

DISTURBANCE, PI CONTROL

1 l

.l

I

d l l l I I —r

0 100 200 400

 

 

 

 

FIGURE 4-6 RESPONSE OF 011(2) TO A STEP LOAD

DISTURBANCE , IMC CONTROL

118

4.4 Model Based Control, MIMO

For a two input, two output system like the Linde Column it is

natural to consider the multivariable form of the IMC structure. In

A

this case the internal model G(z) becomes a matrix of transfer functions
of the form of equation (3-12). Just like the case of a single input,
single output system, the controller can be designed as the inverse of
the process model and a filter used in series with the controller to
address robustness and to shape closed loop response. Garcia and Morari
[30] give a detailed design procedure for the multivariable case.
However, that structure will not be considered here since it requires
uniform sampling of the inputs and outputs. Instead two feedforward

controllers each with its own sampling rate will be investigated.

4.4.1 Feedforward Internal Model Control

There are two different structures of feedforward compensated IMC
that are applicable to a decoupling problem. The scheme offered by
Brosilow [37] is shown in Figure 4-7. An alternative is a modified
version of scheme suggested by Garcia and Morari [30] shown in Figure 4-
8. As we shall see the design equation for the feedforward compensator

fo(z) is the same for both forms, however the dynamic response is not

the same.

119

feedforward d
compensator

 

 

 

disturb-

ance

process
+ m

controller process

internal model

 

 

 

Q1

 

 

Figure 4-7 Brosilow's Feedforward
Internal Model control structure

120

feedforward d
compensator l

'9" (le
L"
+ .' + m +
s O O 6(2) O
+ 4-

controller process

di ingernal
stur ance 2
model Gd ( )

 

 

disturb-
ance
process

———>y

 

 

-
a + '

lntemal model

 

 

 

Figure 4-8 Revised Feedforward
Internal Model control structure

121

Consider first the structure of feedforward IMC shown in Figure 4-8.
In this scheme the effect of measurable disturbances Gd(z)d are
compensated for by fo(z) so that it may cancel their effect on the
process. The disturbance effects are also accounted for by an internal
model Gd(2). The modification developed in this research was to add the
output of fo(z) to the output of Gc(z) rather than to the input of
Gc(z) as presented in [30]. This modification was made so that the
calculation of the Gc(z) could be made independently of fo(z) and thus

allow for simplified decoupling of two control loops (simplified
decoupling will be explained in the next section). This structure is a
unique form of feedfoward IMC, and to this author's knowledge has not
been presented in the literature. Like other forms of feedfoward
control it is able to completely cancel the effects of the disturbance

in some situations.

Since this form of IMC is new it is necessary to derive the form of

the feedfoward compensator fo(z). Straight forward algebra leads to

the following closed loop transfer functions

Gc(2)

 

m(z) - [5(2) ‘ d(z)(Gd(z) ' Gd(z))]

1+gexmn-5e»

G (2)
+ ff A d(z) (4.11)

1 + Gc(2)(G(2) - G(2))

122

GC(Z)G(Z) A
y(z) _ K, [5(2) ’ d(Z)(Gd(Z) ‘ Gd(z))]

1 + GC(Z)(G(Z) - 6(2))
G(2)fo(2)
+ A d(2) + d(z)Gd(z) (4.12)
1 + GC(Z)(G(Z) - 6(2))

Now assuming that G(z) - 0(2) and Gd(z) - Gd(z)

y(Z) - Gc(2)G(2)S(2) + (Gd(2) + G(Z)fo(Z))d(2)

For the disturbance effects to be completely canceled requires
fo(2) - 'Gd(z) / G(Z)

Note that G(z) is not guaranteed to be invertible since it may be non-

minimum phase. In this case fo(z) is designed as

cff(2) - -cd(2) / c (2) (4.13)

As we see, the question of model invertibility is also relevant to
the design of feedfoward compensators for the IMC controller. For non-
minimum phase systems an approximate inverse must be found for G(z) in
order to apply the design equation (4.13). One can use either the

method of [2] or the one described in this thesis. The obvious choice

123

is to use the same approximation employed in the design of the feedback
controller. However, equation (4.13) allows for the cancellation of

time delays. Therefore for non-minimum phase processes fo(z) may be

reduced to

0ff(2) - -c;d(2)c:_(2)’1 2"* (4.132)

where 7* - minlrd , r} and r is the time delay of G(z) and 1d is the
time delay of Gd(z). For such a fo(z) and with Gc(z) - l / G (2) the

resulting closed loop response is given by
-f* -f*
y(z) - G+(z)z 3(2) + d(z)Gd(z)(1 - G+(z)z ) (4.14)

The design equation for the feedforward compensator in Brosilow's

scheme can be easily derived starting from the its closed loop equation,

Gc(z)G(z)
[8(2) - d(2)(fo(2)G(2)+Gd(2)l

 

y(Z) - .
1 + Gc(2)(G(2) - 6(2))

+ d(2)(fo(z)G(z) + Gd(2)) (4.15)

Now if G(z) - G(z) then

y(Z) - Gc(2)G(2)S(2) + d(2)Gc(2)G(2)(fo(2)G(2) - Gd(2))

124

+ d(2)(fo(Z)G(2) - Gd(2))

So as with the previous form of feedfoward IMC the disturbance affects

are completely canceled if fo(z) - ~Gd(2) / G(z).

As before, fo(z) - -Gd(2)G_(z)-lz-f* for non-invertible G(z).This

leads to a different closed loop response
- * - *
y(z) - 0+(2)2 ’ 5(2) - d(z)Gd(z)(1 - 0+(2)2 ' )2 (4.16)

It will be shown that these different closed loop responses are

noticeable in simulations of the Linde Column.

4.4.2 IMC for the Linde Column

In applying these schemes to the Linde Column, the 57th tray
concentration was controlled by manipulating the reflux flow while
changes in reboiler steam pressure were considered as measurable
disturbances. The bottoms concentration was controlled by manipulating
the reboiler steam pressure while changes in the reflux flow were
considered as measurable disturbances. This is similar to the
distillation control strategy described by Shunta and Luyben [32] in
which they used conventional feedback controllers with feedfoward
compensation. The results described here were achieved by simulating

the bottoms controller with a 8 minute sampling interval to coincide

125

with the sample rate of the bottoms concentration measurement while the

other controller had a 1 minute sampling period.

The modification presented in this thesis to feedforward IMC scheme
was made to allow for simplified decoupling of the Linde Column
controllers. Simplified decoupling was first introduced by Buckley [38]
for continuous systems with conventional feedback controllers. A block
diagram of this new IMC structure to support simplified decoupling is
shown in Figure 4-9. In this diagram the Linde Column is depicted by
the transfer functions G

G12, G21, and 622. The top controller

11’

contains the control block GCT whose output is summed with that of the

feedforward block GFT to produce the total top control action mT (reflux

flow). This output signal is also sent to the internal model GT' A
filter, FT’ is placed ahead of the top control block for tuning

purposes. The feedforward signal sent to GFT is the output of the

bottom control block GCB’ This signal is also sent to the internal

A

feedforward model GdT' The output of the feedforward internal model is

summed with the output of the standard internal model, G and the

T’
result is subtracted from the measured top controlled variable
(concentration on the 57th tray) to produce the feedback signal for the
top controller. The feedback signal is compared to the top set point

ST. The difference in these two signals is the input to the top filter.

126

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4-9 Simplified Decoupling with
Revised Feedforward IMC

127

In a like manner, the bottom control action, mB (reboiler steam

pressure), is calculated in the bottom controller using the bottom

control block G the bottom feedforward compensator G and the

CB’ FB’

bottom tuning filter FB' The input to the bottom feedforward

compensator is the output of the top control block. The bottom feedback

A

signal is calculated from the bottom internal model GB, the bottom

A

disturbance model GdB’ and the measured bottom controlled variable

(bottoms concentration). The bottom feedback signal is compared to the

bottom set point 83.
This scheme allows the feedback control calculations for both

controllers to be made independently and then their results processed

and added to each others's to give the overall control signals.

Brosilow's version allows for the same simplified decoupling. A block
diagram is shown in Figure 4-10. The difference with Brosilow's version
is the lack of an internal disturbance model and the signal applied to
the standard internal model. As before, the top controller contains the

control block GCT whose output is summed with that of the feedforward
block GFT to produce the total top control action mT (reflux flow).

However, only the output of the control block GCT is sent to the

internal model GT’ A filter, F , is placed ahead of the top control

block for tuning purposes. The feedforward signal sent to GFT is the

128

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4-10 Simplified Decoupling with
Brosilow's Feedforward IMC

129

output of the bottom control block GCB' The output of the internal

A

model, GT’ is subtracted from the measured top controlled variable

(concentration on the 57th tray) to produce the feedback signal for the
top controller. The feedback signal is compared to the top set point

ST. The difference in these two signals is the input to the top filter.

The bottom controller operates in the same fashion which again allows
independent feedback calculations of both controllers. These feedback

control signals are again fed forward to the other controller for

compensation.

The multirate sampling nature of the two controllers gives rise to a
question of how to use, if at all, the values of the feedback

calculation of the overheads controller, mCT(kT), which occur in between

the sampling instants of the bottoms controller, 8T, 16T, ect. In the
simulations of the Linde Column it can be observed that much smoother

control is achieved if the feedforward portion of the bottoms controller
was based on the trailing average, “CT(kT) rather than only on

mCT(8T), mCT(16T)’ ect. The trailing average of mCT(kT) can be

expressed as

1 k
X
mET(kT.N) - _N— m

1-k-N91

(1T) k
N

130

The effect of the trailing average (when N - 8) is to better approximate
the assumed 8 minute zero order hold in the reflux flow that is inherent

in the pulse transfer function representation of G21(z).

Simplified decoupling as detailed above fails to fully decouple the

two control loops for two reasons. (1) GFT(2) cannot be designed as a
perfect feedforward compensator because G11(z) is non-invertible. (2)

Only the feedback portion of the total control signal from either
controller is fed forward to the other controller. Nothing can be done
to eliminate the first problem, however the multirate nature of the two

controllers allows for nearly eliminating the second problem.

The basic idea behind the following strategy, which will be called
critical decoupling, is to take advantage of the fact that the output of
the bottoms controller remains constant for seven intervals of the
overheads controller because the former has an 8 minute sampling period
while the latter has a 1 minute sampling period. In addition, this
strategy takes advantage of the trailing average used for the

feedforward signal of the bottom controller.

Because the bottom controller is held constant for seven periods of
the overhead controller, its total control signal may be used in the
feedforward portion of the overheads controller for those seven periods.
In fact its total control signal may be fed forward at all times if the
bottoms controller calculations are executed before the calculations of

the overheads controller.

131

Likewise, nearly all of the overheads controller output may be used
in the trailing average calculation of the feedforward portion of the
bottom controller. The trailing average can be based on the previous 7
overheads controller outputs which are obviously known and the current
feedback portion of the overheads controller. The trailing average of

the overheads controller may now be expressed as

1

-1
mT(iT) + -E- mCT(kT)

1 k
X

(sz8) - —
"I 8 1k

This more exact form of decoupling is illustrated for both types of

feedforward IMC in Figures 4-11 and 4-12.

4.5 Results and Discussion

The proposed control strategies were tested through simulations of
the Linde Column using the model developed in the last chapter. The
simulations tested 4 different strategies for rejection of step load
disturbances in both the overheads composition and the bottoms

composition. The 4 strategies tested were:

Independent PI ................... Strategy 1
Independent IMC .................. Strategy 2
Simple Decoupled IMC ............. Strategy 3

Critically Decoupled IMC ......... Strategy 4

132

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4-11 Critical Decoupling with
Revised Feedforward IMC

133

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4-12 Critical Decoupling with
Brosilow's Feedforward IMC

134

Table 4-1 lists the transfer functions for the various components of the

IMC strategies. The two PI controllers used are

-2.1 + 2.02'1

-1
z

 

Top: GC(z) - 1 -

-378.7 + 355.02'1

1

 

Bottom: GC(z) - _
1 - 2

It is interesting to note that because G22(2) is first order with no

time delay the PI controller used for bottoms composition control is
equivalent to an IMC controller without feedforward compensation. See

the Appendix for the proof.

It was found that the Linde Column remained stable under PI control
for the disturbances considered. In fact the level of interaction
between the two control loops was consistant with the experience of the
operating personnel. The response of the 2 independent PI controllers
is illustrated in Figures 4-l3a and 4-15a for a step load disturbance on
the 57th tray, and in Figures 4-l4a and 4-l6a for a step load
disturbance in the bottoms concentration. The interaction between both
control loops is quite evident for both types of disturbances. An
example of the interaction is seen by the deviation from set point for
the bottoms concentration when the load disturbance occurs at the 57th

tray in Figure 4-15a.

0

C,T

F,T

07> 67> C)

D,T

0

0,0

F,B

8,0

135

Table 4-1 IMC Transfer Functions

2 - 0.923

 

- 0.03532

3.0872 - 2.850

3 2
Z - 0.4772 + 0.3882 - 0.866

3 2
0.0332 - 0.0202 +0.0002382 - 0.0507

4 3 2
2 - 0.8272 10.3882 - 0.9672 + 0.481

0.1 0902

3 2
Z - 0.4772 + 0.3882 - 0.866

2 - 0.937

 

-0.002642

0.06022 3+ 0.02612 2. 0.02212 - 0.518

3 2
Z - 0.5882 - 0.3512 + 0.938

-0.00264

 

2 - 0.937

0.0001592 2+ 0.0002182 + 0.000146

3 2
2 - 0.5882 - 0.3512 + 0.0938

 

mquEPFm—c c<c._ awhm < o...
zc—...<~_._.zuuzcu ><~_._. :hnm 5:. ...c mmzcmmu: .m—Ic maze—...

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

136

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.8 Sn 3. = .3 ss S. .3 .8 3. :

LIPIPIPIIPIEI npnlbI—D-DIDII’e? .... pun. bb-b nub. b... I... 1's...
“up. No... lly/ H...-
m... ,4 m...

H... 8: . . H... 2:

u .0...— ” .05

H . .emau . H . mnaz_m

D. .-

m ... . m .6
H... . ..-
n... . I n..-

.8 2s 3. = .3 e3 e3 .8 Se 3. ._
...p ...!F-.pF-.#F.0..OI. .... in. sun. .b-. .... --—lu.--.

l’I/ m... l/ W...
H?- //= lll mun.-
... A...

l 11w... .9: I .lw... E

m .n H ....
In . ... W...
a: 2 He.-

 

137

PDPD

D...

but.

zo~h<mhzuuzou machhom mzh z— muz<mm=hmna o<cn amen < Op
><zh gusm m1» zo zo~h<zhzmuzou mc mmzoammz

IPIPDD

8-
supp-LL».

 

I

a...

 

IIII

 

 

_
<
_

IIII

 

 

 

 

 

 

 

 

 

DPDD

D...

RD.»

0....

 

 

IIII

 

IIII

 

 

III.

 

IIII

 

 

 

 

 

 

 

 

 

...
e
.

u:

.3.

.tamu
...

Us"
.2

«fine umawhm

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

E. Ian .0: E ‘0 E. :
P-Db p-b- mph» Dublb mph-.5...

— I

T

I

0

I

n

I — I

I

j

I

'0 0n 8.. i ‘0 s. :
b-bh pbbb bib. DP-h Dhbb-nbrb

I

"AA, A

. I

0

T

1

r

. e

0

I

1

I

 

 

 

 

 

 

 

 

0.01

0.01

uzm
.u.o
mndzmm
.0

Ha
.m

138

.><~:. :hnm E: zc uuz<ma=pm8
a5: dwhm < E zo~h<mhzmuzcu mzotbm mo umzoamum 21¢ maze...

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

on. can 8.. ! Dan 0.. : 8. Con .09 i SO 0.. :
...? pp». p... -.L- ..pp p... 00.01 uppp .PL. up». .... .... ..b- g...
T I e
l 'o. 1? I .
< . < . 2 e
8.. .uw. - . 2.. 2.:
H u a . dd
. :mu - mix—m
T . .... .c I . . . . 2.. .u
T 0
. a... . - a e...
8e a: 8' 3n 8. 8. : Se 3. 8' e8 e8 8. =
hb-P ppbblbbPIRn-bP b-blhlh-bb ‘e... DDLP Dub. but bRbblb-n- nPrr g...
I '6 111111. 3..
I :0. ‘0.
n 2: . . H E
H .a f .m
0... . O...
I . . a... A I a...

 

 

 

 

 

 

 

139

muz<m~=fim3 23 3.5 < o...
zoCSfizmuzou mzotbm ...—c mmzodmmm 31¢ mass...

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

S. ... .... .8 .... . .2 : ... ... .... ... ... ... :
.... .... .... .... .... .... ‘e? ...m P... .m.. .PP. .m.. P... ”e,
11 . . .... ... . .18..
H O f .H ‘0. U.
.1... s: . a...”
. .u.o - 2382..
H ..Cmu . L H
1.... .u . . ...a... .u
. .... .. . 2...
... ... ... ... ... ... : ... ... .... ... ... ... :
IPLIIPIPI.... .... P... .... ..L. g..l r..L P.P. .... ..P. .P.. .... ‘..I
1.111! ; - .... ll...) . ....
T I
K.“ . . .... . ....
g n 8: .. H E
. .... .n .12.. .m
f I
- ...... - ....

140

Control Strategy 2 demonstrated that the major hurdle to improved
control is the coupling between the two control loops. Figures 4-13b
through 4-l6b show the response of control strategy 2. These
simulations indicate that a slight improvement is made in the
performance of the controller. A comparison of Figure 4-l4a with Figure
4-l4b shows that the oscillations are reduced and set point is achieved
sooner. A comparison of Figure 4-13a and Figure 4-l3b leads to the same
conclusion. However, no noticeable improvement is seen in the
performance of the bottom control loop. Compare Figure 4-lSa to Figure
4-15b and 4-l6a to Figure 4-16b. This is to be expected since the PI
controller on the bottoms composition is equivalent to an IMC

controller.

The addition of simplified decoupling in Control Strategy 3 greatly
improved the simultaneous control of both product concentrations. The
response of Control Strategy 3 is shown in Figures 4-13c through 4-l6c
for the feedforward IMC introduced in this' thesis. Figure 4-17
illustrates the results using Brosilow's version of feedforward IMC.
The improvement in performance due to the feedforward compensation is
very evident. This is especially true for the controller which is only
responding to the interaction between the two controllers. For example,
when the bottom controller reacts to the disturbance produced by the
overhead controller correcting for a step load change in the 57th tray
concentration. This is to be expected since the majority of the
overheads control signal is due to the feedback control block. It is
precisely this signal that is used by the feedforward compensator in the

bottom controller. Figure 4-15c and 4-l7c both show dramatic reduction

 

141

.zoaeaaezmuzou ><se seem Au.sv

 

 

 

 

 

 

 

 

 

 

 

 

 

m..-I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

can can no. nun can no. :
.... .... .... a... P...

< r

f

j

r

T

r

I

.06 com on. non can on. .—

.... .r.. ..L. .... .... I

1

T

I

1

T

m

I

r

n

r.

T

.I

r

...

0

I

 

 

“3
9

as

an

zoaeemezmezou mzcetom 18.2V
HE muzemmstmso a<oS anew
< oe Azoemmm> m.3¢323ommv 0:2 omsazoumo >sazam Lo mmzoammm A2-s assume

L

.

.|.

.P..

....

P...

OI—

..L.

 

TIIT

 

III

 

jIIr

 

 

 

 

 

 

 

I

 

 

r...-

....

.P...

 

IIII

 

 

IIIT

 

IIII

 

 

IIII

 

 

 

 

 

 

 

 

IIII

 

 

2... u.

a...
me.-
.zc2e<mezmuzou
... mzoeeom
.6
a...
n..:
e.-
u.-
0..
.zo2eemezmuzou
><st
seam
e.~
A.
u.~

142

in the deviation from set point for the bottoms concentration compared
to classical PI control, Figure 4-15a and independent IMC control,
Figure 4-15b. The improvement in the overheads controller is not as
dramatic. Compare Figure 4-l4b with Figure 4-l4c and Figure 4-l7b. The

lack of improvement is due to the non-minimum phase nature of G11(z).

When the controller which is correcting for the step load
disturbance is considered, the response of Control Strategy 3 is nearly
identical to that of Control Strategy 1 and Control Strategy 2. For
example, compare Figure 4-l6c to Figure 4-l6a and Figure 4-16b for the
bottoms controller. This is to be expected. Under these circumstances
the controller in question does not receive much feedforward
compensation because the controller at the opposite end of the column is
almost totally operating as a feedforward controller. Control Strategy

3 only feeds forward the feed back portion of the other controller.

Control Strategy 4 which makes use of the most information about the
column shows the best performance of all the strategies considered. The
response of Control Strategy 4 is shown in Figure 4-l3d through Figure
4-l6d for the revised feedforward IMC and in Figure 4-18 for Brosilow's
version. A marked improvement is seen for the case of either controller
responding to a step load change at its end of the column. For example,
compare the response of the overheads controller to a step load change
in the 57th tray concentration, Figure 4-13d, to that of Control
Strategy 3, Figure 4-l3c. For this case the addition of the feedforward

portion of the bottoms controller to the feedforward signal of the over

Era-“1

143

+...

.zo2t<mezmezou ><mt seam 16.63

....

...-

SM
hep.

....

 

WIT!

 

II‘II

 

 

 

 

 

 

 

 

 

 

 

....

....

....

....

 

 

IIII

I
In
9

 

 

IIII

 

IIII

 

IIII

 

IIII

 

 

 

 

 

IIII

 

 

 

 

zo_e<aezmuzou mzcteom “2.22

Hz. uuz<mm=hmfio o<on dmhm < op
Aonmmm> m.3onfimommv oz. ownd=oomo >44<unh~mu mo mmzoammm wane mmzwmm

As

....

P...

..t

.a
...L

 

 

 

 

 

 

 

 

 

 

 

 

...l.l

 

..Pr

 

IIII

 

IIII

 

IIII

 

II II

 

 

IIII

 

 

 

 

 

 

 

 

IIII

 

 

M...

D..-

9..
.

zomh<mhzmuzou
mZCFFOm

Au

onh<mhzmuzou

><mh
:unm

A.

144

heads controller has allowed a much quicker return to set point of the
57th tray concentration. This same improvement is seen in the bottoms
controller (compare Figure 4-l6d to Figure 4-16c). This control

strategy has nearly decoupled the two control loops.

4.6 Summary of Chapter Four

This chapter has dealt with the design of a multivariable, model
based control system for the Linde Column. The design was influenced by
a combination of practical considerations and the dynamic model of the
column. Internal Model Control was selected as the alternate control
strategy due to the practical benefits of the current PI control
strategy. Adapting IMC for the specific characteristics of the Linde
Column led to three significant contributions to the literature on the
theory of Internal Model Control. These contributions are: l) a
technique of reduced order controller design for non-minimum phase
systems, 2) a new feedforward IMC structure, 3) a strategy for

multirate, multivariable use of IMC.

The contribution of reduced order controller design is significant
because it simplifies the implementation of IMC for non-minimum phase
systems. This allows IMC to retain its simplicity for all types of
systems. It was found that in the case of the Linde Column there was
very little difference in the performance of the simple reduced order

controller and the more complicated designs that exist in the

145

literature. This new approach is very appealing from a practical point

of view.

The contribution of the new feedforward IMC structure was
significant because it offers an alternative to the existing feedforward
structures to allow simplified decoupling of control loops. It was
found that this new form of IMC offered improved control of the Linde

Column over the existing form.

The contribution of multirate IMC control was significant because it
addresses a common problem found in the process industries, namely non-
uniform sampling of process variables. This allows IMC to be used for
an even greater variety of process control problems. This extension is

an example of the contributions that can be made from application

oriented research.

Finally the results of this chapter demonstrate that improved
control for the Linde Column is possible through an alternative strategy
that is appealing to the Operating personnel. It is this final
contribution that has spawned further work within the Dow Chemical

Company.

Chapter Five

SUMMARY AND RECOMMENDATIONS

This thesis has described research regarding the modelling and
control of an industrial scale distillation column, named the "Linde
Column", located at the Dow Chemical Company's Michigan Division.
Justificaton for this research can be found in two areas. One is the
possible economic impact, improved control would have on the operation
of the column. The second is the contribution this research makes to
bridging the gap between advanced control theory and its practice in the
process control industry. The overall objective of this thesis was to
determine if advanced control would be beneficial to the Linde column.
This objective led to two major themes in the research: 1) develop a
dynamic, multivariable model of the column, and 2) propose an alternate
control scheme and test it through simulation. The remainder of this
chapter will
separately summarize the results of both phases of research. Each
summary will be followed by suggested future research on particular
topics that were addressed. Recommendations on future opportunities for

research on a broader scale will then follow.

146

,F H —.5

147

5.1 Summary of Linde Column Modelling

The research that led to a dynamic model of the Linde column was
described in Chapter Three. The model developed is a 2 input, 2 output
matrix of discrete time transfer functions. The preliminary problem of
estimating the order and time delay of a multi-input, multi-output model
was
first addressed. A novel identification procedure was developed to
solve this problem. Parameter estimation of the identified model used
real Operating data. Non-minimum phase behavior was discovered in a
part of the column which is not mentioned before in the literature. The
new identification procedure and the discovery of non-minimum phase

behavior are the two major contributions of Chapter Three.

The novel identification procedure presented in Chapter Three is a
significant contribution to the literature on model identification of
multivariable systems. This new approach greatly reduces the number of
parameters necessary to describe the behavior of multivariable systems
with multiple time delays and is based on a new "delayed polynomial
matrix" representation of discrete systems. The delayed polynomial
matrix representation explicitly accounts for multiple delays in a
multivariable system. This is in contrast to the conventional
polynomial matrix model which represents time delay by parameters with

values equal to zero.

The delayed polynomial matrix method was developed by starting with

a generic state space description of a discrete time multivariable

148

system with multiple input delays. Theorem 3.1 derived the relationship
of the state space description to an input/output description with
multiple input delays. Lemma 3.1 then demonstrated the relationship of
this input/output description to the normal polynomial matrix
representation. The delayed polynomial matrix representation was
proposed by Remark 3.1 which showed how it reduces the number of
parameters in an input/output model of a discrete system. A two phase
identification procedure was then outlined. This procedure is based on
testing the linear dependence between input and output observations. It
first identifies the normal polynomial matrix model from input/output

data and then derives the delayed polynomial matrix model.

A computer algorithm for the delayed polynomial matrix method was
then described in a way that made it suitable for interactive use.
Singular value decomposition was chosen as the numerical tool for
testing linear dependence. This is the first time singular value
decomposition has been used in an identification method. Two techniques
were presented to filter noisy data. Finally, a simulated example
demonstrated the ability of this new identification method to perform

with up to 20% additive noise on the data.
5.1.1 Results of Linde Column Modelling
Chapter Three then led to developing a model of the Linde column

from real operating records of the process. The problem of obtaining

good data for model building was first addressed. A review of the data

149

showed that it contained sufficient information for modelling. However,
it also put constraints on the form of the model. The delayed
polynomial matrix method of model identification was used where
appropriate. It was found to be in close agreement with other
identification techniques. Least squares estimation was then used to
fit the parameters of the identified model. Discussion of the Linde
column model showed that certain aspects of the model were very

consistent with the experience of the operating personnel.

This research has produced an important discovery concerning the
dynamics of distillation columns in addition to the original
contributions to the theory of model identification and controller
design. It was found that the concentration on the 57th tray exhibited
inverse response to changing reflux flow. This was totally unexpected

and has not been mentioned in the chemical engineering literature.

5.1.2 Future Research on Distillation Modelling

The modelling effort of this thesis provided much insight on the
dynamics of the Linde column. However, some questions remain
unanswered. The differences in the model parameter values obtained from
different data sets should be resolved. This can only be done by using
data collected during controlled experiments. The present model of the
column can be used to determine the amplitude and frequency content of
the experimental inputs. The model developed in this thesis is an

excellent starting point for future research.

m
l

150

The delayed polynomial matrix method of identification requires
further investigation. Alternative noise compensation techniques should
be researched. Incorporating the instrumental variable approach to
handle correlated noise would enhance the application of the delayed
polynOmial matrix method. Extending the method to identify a noise
model would be another way to compensate for colored noise. Since noisy
process data is normal in the process industry, these suggestions would

go far toward improving the practical application of this method.

5.2 Summary of Advanced Control for the Linde Column

An alternative control scheme for the Linde column was considered in
Chapter Four. The alternate control strategy was selected after a
careful examination of the advantages of the current control. Tailoring
the new control approach to the model of the Linde column led to several

original contributions.

The alternate control strategy selected is a feedfoward version of
Internal Model Control, IMC. The essence of IMC is the concept that
control is applied to a model of the process as well as the process.
Assuming the model and the process are well matched, differences in
output of the two represent a measure of the process disturbance. Such

a scheme offers several appealing features:

1. Time delay compensation is inherent since thecontroller is only

acting upon processdisturbances.

151

2. IMC is understandable to operating personnel because they can
identify with the way that the internal model in the IMC

strategy predicts the behavior of the process.

3. Robustness is also addressed in the structure of IMC by the
addition of a filter ahead of the control calculation. The
filter time constant can function as a single tuning constant,
thus the operator can respond to unexpected process shifts

without resorting to manual control.

4. Easy controller design is another appealing feature of IMC.
For well behaved processes, the controller is just the inverse

of the process model.

Because IMC addresses these practical issues in a straight forward
manner, it was selected over controllers based on more abstract

approaches, such as state space or optimal control theory.

Two peculiarities of the column required extensions of existing
theory on Internal Model Control. First, one of the two controlled
variables, the bottoms concentration, could only be sampled every 8
minutes due to the instrument used. The other controlled variable, the
concentration on the 57th tray, which could be sampled as often as
wanted, was sampled once a minute to adequately cover its dynamic
response. Therefore the multivariable version of IMC could not be
directly applied to control the column, because it assumes uniform

sampling for all outputs. A multirate feedfoward strategy was developed

 

 

 

152

to address this problem. Secondly, the concentration of the 57th tray
exhibited unexpected inverse response with respect to the reflux flow.
As a result, the model of this input / output relationship has at least
one zero outside the unit circle. A controller designed as the inverse
of the model would be unstable. To address this problem, a straight
forward approach was developed using a reduced order model of the
process. This simpler approach maintains the overall simplicity of IMC

for non-minimum phase processes.

Reduced order controller design for non-minimum phase systems is an
important contribution of Chapter Four. Compared to the literature,
this approach simplifies the practical implementation of IMC for non-
minimum phase systems. One solution cited in the literature is to
eliminate the unstable zero before inversion [2]. This was tried
through simulation. The result was good regulation but the input
applied to the process was highly oscillatory and therefore
unacceptable. Another solution cited [2] solves a predictive control
problem. This solution relies on the formulation of a least squares
approximation problem and its solution and involves several tuning
parameters. The approach taken in this thesis was to approximate the
process by a simple minimum phase system which was then inverted for the
controller design. Simulation showed little difference in performance

between the reduced order controller and the more complex designs.

Another contribution of Chapter Four was a modification to
feedforward IMC to allow for simplified decoupling control. This

modification is significant because it offers an alternative feedforward

153

structure. The closed loop response was mathematically derived to Show
how it differed from existing feedforward IMC. Simulations showed that
it produces improved control for some types of disturbances on the Linde

column.

Decoupling control of both product streams on the Linde column was
accomplished by applying the new feedforward structure in a multirate
fashion. This is a significant contribution because it offers a way to
apply multivariable IMC to non-uniform sampling problems. This class of
problems is common in the process industry. The multirate extension of
IMC is an example of the contributions made by application oriented

research.

5.2.1 Results of Advanced Control Simulations

Once the multivariable control strategy was designed, simulations
were used to compare it to the conventional control using independent PI
controllers at the bottom and the at 57th tray. The simulations tested

the response of the controller to step load changes.

The multivariable controller out performed the conventional
controllers in the load change simulations. Improved disturbance
rejection was achieved in both product streams by the multivariable
controller. This is not surprising since the multivariable controller
made use of the knowledge of the dynamic coupling and non-minimum phase

behavior in the distillation column. This is not to say that PI control

154

is bad or that it could not be improved. In fact, the model of the

column has interesting implications for the current PI control.

5.2.2 Future Control Research

At the very least the present controllers could be tuned using
the model developed in this research. The present controllers were
tuned manually with the main concern of getting the process to run
acceptably but not necessarily optimally. So there is room for
improvement here. Of course, techniques exist to compensate PI
controllers for deadtime or coupling. These could be used to improve

control.

The inverse response on the 57th tray suggests an alternate
controlled variable. It was originally reasoned that the control system
has an advantage if it controls at the 57th tray since any disturbance
detected there must propagate up another 15 trays to affect the output
of the column. Therefore, there would be extra time to respond to
disturbances. However, the performance of a PI controller is adversely
affected by inverse response, so such a controller applied at the 57th
tray may not reject disturbances as well as a controller applied at the
top of the column where minimum phase behavior can be expected. Control
of reflux flow based on distillate concentration using the concentration
on the 57th tray as a feedforward term might perform better. Control
improvements are available in several different forms because of the

knowledge obtained in this research.

155

The performance demonstrated by decouplng control in the simulations
of the Linde Column suggest more research be conducted on the control of
moderate to high purity distillation columns by using an interior tray
concentration as a control point. Weischedel and McAvoy [l4] argued,
and then demonstrated with simulations, that decoupling of moderate to
high purity columns is not feasible with linear decouplers. The
simulations of the Linde Column show that decoupling is feasible if the
concentration on an interior tray is controlled. Decoupling is possible
for this situation because the separation at an interior tray is low
purity. The ability to control a high purity column with this approach
does come at a cost; that cost is the non-minimum phase behavior that is
present at the interior tray. However, simulations of the Linde Column
show that a model based control approach can effectively handle this
type of behavior. Future research on a three point control strategy
which maintains the concentration in the bottoms product, the overheads
product, and an interior tray as an intermediate should prove fruitful

for high purity columns.

The control concepts presented in Chapter Four offer new
opportunities for further research on IMC control and for the
application of IMC to other process operations. Other types of
predictive control should be considered to handle non-minimum phase
systems. For example, the controller assuming no change in its output,
could predict the position of the process beyond the initial inverse
response. It could then calculate a change in its output necessary to
drive the process to set point over the prediction horizon. The

prediction could be based on a first order model with sufficient time

156

dealy or on a transfer function equivalent to the internal model. The
estimate of the process disturbance available in the IMC structure could
also be further exploited. Methods to test the whiteness of the
disturbance (and therefore not respond) or to predict the disturbance
(if it is correlated) would improve the performance of IMC. The
trailing average used in the multirate version of IMC could be further
investigated. Alternative averages such as a weighted average should be
tried. This thesis, like all successful research, raises new questions

even as it obtains its objectives.

5.3 Future Research Opportunities

This thesis quite clearly demonstrates that efficient control of
industrial processes is a difficult and challenging problem. Global
competition is enforcing tighter requirements on yields and production
costs. The ability of computer automation to improve yields and reduce
production costs, using standard control techniques, has made it a
standard entity in every competitive process operation. This thesis
illustrated how a closer look at the process can reveal dynamic
characteristics which have major implications on process operation. It
also demonstrated that a model based control approach could be modified
so that it more effectively addresses industrial control problems. The
task of the research community is now to demonstrate ways that greater
improvements can be realized through more sophisticated use of the
computer/automation base that is part of every serious process

operation.

It. '.-_I '.r.- .s "’ E
.'

157

Many opportunities exist for research and development of computer
based methods which compliment or extend the standard control techniques
now in use throughout industry. Some of the promising areas for

industrial application will now be mentioned.

Artificial intelligence and expert systems offer solutions to
problems that are not easily addressed by standard control programming.
In many situations of industrial control, the objective of the
controller is not well defined, or the control now implemented by the
human operator does not translate well into precise mathematical rules.
An expert system based on heuristic rules can provide control strategy
for these situations. Monitoring of standard control algorithms and
logic is another use of artificial intelligence. The area of process
diagnostics has seen a lot of development activity recently.
Supervision of feedback controllers by an expert system appears to be a
fruitful area for research, especially for adaptive controllers.
Research on the use of artificial intelligence to assist human operators
during situations of multiple alarms also has a high potential for
payoff. In the future, it will be the marriage of artificial
intelligence and mathematically based monitoring and control that will

prove to be the most effective use of both approaches.

A drawback of expert system programs and standard control programs
is that they are application specific. The programs must be developed
individually for each production plant. Duplication of similar programs
is a partial solution to this problem, but the real answer lies in

computer automated program generation. The ideal program generator

158

would use graphic representation of the process as input and produce a
useable control program (which can be further customized) as output.
Research in this area would be a real benefit to reducing the cost and

time needed to automate production processes.

Greater knowledge of process dynamics and process disturbances is
required if automation is to move beyond the standard control
techniques. This usually requires a mathematical model of the process
for the vast majority of advanced process control techniques. Excellent
research opportunities exist in the area of online process modelling and
simulation. In the future, process automation systems will have
dedicated computer resources for predicting process behavior for control
just as in Internal Model Control, but on a much broader scale. The
routine application of advanced control is very much dependent on an
integrated approach to model building, online process simulation, and

updating of process models.

Supporting technologies of advanced process control, namely
computers and instrumentation, are as important to its future as its
methodologies and algorithms. The complexity of the processing required
for advanced process control coupled with the falling prices and
increasing power of micro-computers is leading to distributed processing
and distributed control. Complete plant wide automation from sales to
shipment will require reliable and cost effective computer systems which
include: 1) transparent high speed communications, 2) shared computer
memory, 3) fault tolerance and 4) application programs which can be

designed and tailored by the user. Sensing instruments which have

159

always been the subject of aggressive research and development will
continue to offer fruitful opportunities for research. Instruments
which can offer self calibration and self diagnostics, multiple
measurements, and digital communications will no doubt be the standard
fare in the future. A comprehensive hardware base is necessary for the

routine application of advanced process control.

The discussion above has only outlined a few areas of research that
are desirable from an industrial point of view. An important area,
which was the subject of this research, is the study of current advanced
control techniques on industrial applications. Investigations of this
kind can be the basis for productive joint ventures between industry and
academia. It is industry which can offer challenging control problems
and it is academia which can provide fresh approaches where the standard

ones have performed unsatifactorily.

A paradox of computer based control is that one of its greatest
assets, flexibility, can also be one of its greatest detriments. This
stems from the ease at which ad hoc modifications can be made to
standard PID controllers in an attempt to address its short comings.
This often leads to overly complex control programming which only works
marginally at best. The greatest contribution made by this thesis to
its industrial sponsor, The Dow Chemical Company, is not the potential
improved control of the Linde column, rather it is the investigative
approach it used to solve the problem. Already several significant
contributions to the profitability of the company have been made by

solving control problems by first gaining a deeper understanding of the

160

process through experimentation and model building; rather than adding
more 'fixes' to the existing control scheme. These early successes

indicate a positive future for the industrial application of advanced

control.

APPENDIX

APPENDIX

The IMC control structure can be converted to conventional feedback

form by defining C(z), the feedback controller as

Gc(2)
G(z) -

 

l - Gc(z)G(z)

where Gc(z) is the control block in the IMC structure and G(z) is the

internal model.

Now a first order system

 

 

 

. 621
0(2) - G(Z) - 1
1-22
produces
1.221 1 a _1
b ‘b—‘Tz
1-221 621 1-21
1-
b 1 _ ail

161

162

which is a discrete form of a proportional / integral controller.

LIST 0]? REFERENCES

10.

11.

LIST OF REFERENCES

McAvoy, T.J., "On the Difference Between Distillation Column
Control Problems", presented at the Second Engineering Foundation
Conference on Process Control 1981.

Garcia, C.E. and Morari M., "Internal Model Control. 1. A Unifying
Review and Some New Results”, Ind. Eng. Chem. Proc. Des. Dev., 21,
308-323, (1982).

Williams, B.J., “Statistical and Theoretical Modelling of a Binary
Distillation Column", Conference on Industrial Applications of
Dynamic Modelling, Durham, IEE Conference Publication No. 57,
(1969)

Krishnamoorthy, V. and Edgar, T.F., ”Identification and Estimation
of Discrete Models for a Distillation Column", University of Texas
Research Report No. WP28- 3:50

Foulard, C. and Bornard, 6., “Identification and System Parameter
Estimation of Distillation Processes - A Case Study", Proc. 3rd
IFAC Symp. on Identification and System Parameter Estimation, The
Hague, (1973)

Gauthier, A. and Landau, I.D., ”On the Recursive Identification of
Multi-Input, Multi-Output Systems", Automatica, 14, 609-614, (1978)

Gustavsson, I., "Survey of Applications of Identification in
Chemical and Physical Processes”, Automatica, 11, 3-24, (1975)

Wood, R.K. and Berry, M{W., “Terminal Composition Control of a
Binary Distillation Column", Chem. Eng. Sci., 28, 1707-1717, (1973)

Rijnsdorf, J.E., "Interaction in Two-Variable Control Systems for
Distillation Columns - I", Automatica, 3, 15-28, (1965)

Luyben, W.L., ”Distillation Decoupling", A.I.Ch.E. Journal, 16,
198-203, (1970) ‘

Dahlqvist, S.A., ”Control of a Distillation Column Using Self-

Tuning Regulators", Canadian Journal of Chem. Eng., 59, 118-127,
(1981)

163

 

12.

13.

14.

15.

l6.

l7.

l8.

19.

20.

21.

22.

23.

24.

25.

26.

27.

164

Ogunnaike, B.A. and Ray, W.H., "Multivariable Controller Design for
Linear Systems Having Multiple Time Delays", A.I.Ch.E. Journal, 25,
1043-1057, (1979)

Ogunnaike, B.A., Lemaire, J.P., Morari, M., Ray, W.H., ”Advanced
Multivariable Control of a Pilot-Plant Distillation Column",
A.I.Ch.E. Journal, 29, 632-640, (1983)

Wieschedel, H., McAvoy, T.J., "Feasibility of Decoupling in
Conventionally Controlled Distillation Columns", Ind. Eng. Chem.
Fundam., 19, 379-384, (1980)

Freeman, T.G., "Selecting the Best Linear Transfer Function Model",
Automatica, 21, 361-370, (1985)

Astrom, K.J., "System Identification - A Survey", Automatica, 7,
123-162, (1971)

Gustavsson, I., "Comparison of Different Methods for Identification
of Industrial Processes", Automatica, 8, 127-142, (1972)

Akaike, H., "Statistical Predictor Identification", Ann. Inst.
Statist. Math., 22, 203-217, (1970)

Akaike, H., "Use of an Information Theoretic Quantity for
Statistical Model Identification”, Proc. 5th Hawaii Int. Conf. Sys.
Sci., 249-250, (1972)

Lee, R.C.K., ”Optimal Estimation, Identification and Control",
M.I.T. Press, 1964

Woodside, C.M., "Estimation of the Order of Linear Systems",
Automatica, 7, 727-733, (1971)

Wellstead, P.E. and Rojas, R.A., "Instrumental Product Moment
Model-Order Testing; Extensions and Application", Int. J. Control,
35, 1013-1027, (1982)

Buckley, P.S., Luyben, W.L. and Shunta J.P., "Design of
Distillation Column Control Systems”, Instr. Soc. America, (1985)

Budin, M.A., "Minimal Realization of Discrete Linear Systems from
Input-Output Observations", I.E.E.E. Trans. Auto. Control, AC-l6,
395-401, (1971)

Guidorzi, R., "Canonical Structures in the Identification of
Multivariable Systems", Automatica, 11, 361-374, (1975)

Klema, V.C., and Laub, A.J., ”The Singular Value Decomposition: Its
Computation and Some Applications", I.E.E.E. Trans. Auto. Control,
AC-25, 164-176, (1980)

Golub, G.H. and Reinsch, C., ”Singular Value Decomposition and
Least Squares Solutions", Numer. Math., 14, 403-420, (1970)

 

28.

29.

30.

31.

32.

33.

34.

35.

36.

37.

38.

165

Smith, B.T., et al., "Matrix Eigensystem Routines - EISPACK Guide",
Lecture Notes in Computer Science, 6, 2nd ed., 1976

Dongarra, J.J., Moler, C:B., Bunch, J.R. and Stewart, G.W.,
"LINPACK, Users' Guide", 1979

Gacia, C.E. and Morari, M., "Internal Model Control. 2. Design
Procedure for Multivariable Systems", Ind. Eng. Chem. Process Des.
Dev., 24, 472-484, (1985)

Brosilow, C.B., "The Structure and Design of Smith Predictors from
the Viewpoint of Inferential Control”, Joint Auto. Con. Conf.
Proc., Denver, (1979)

Shunta, J.P. and Luyben, W.L., "Sampled-Data Noninteracting Control
for Distillation Columns”, Chem. Eng. Sci., 27, 1325-1335, (1972)

Jamshidi, M., "Large-Scale Systems; Modeling and Control", North-
Holland, 1985

Parrish, J.R. and Brosilow, C.B., "Inferential Control
Applications". Automatica, 21, 527-538, (1985)

Sinha, N.K. and Pille, W., "A New Method for Reduction of Dynamic
Systems", Int. J. Control, 14, 111-118, (1971)

Sinha, N.K. and Bereznai, G.T., ”Optimum Approximation of High-
Order Systems by Low-Order Models”, Int. J. Control, 14, 951-959,
(1971)

Brosilow, C.B., Inferential Control Seminar Presented at Dow
Chemeical Company, 1985

Buckley, P.S., Chemical Engineering Seminar Presented at Ohio
University, 1967

 

"I7'11[llILlI'llllITllllllIlIT