Full text of DOI:10.1002/aic.690461010

PROCESS SYSTEMS ENGINEERING
Two-Step Procedure for Data-Based Modeling for
Inferential Control Applications
Raja Amirthalingam, Su Whan Sung, and Jay H. Lee
School of Chemical Engineering, Purdue University, West Lafayette, IN 47907
A two-step procedure for building an inferential control model, which uses both his-
torical operation data and plant test data, is proposed. Moti®ation for using the two
types of data is gi®en, and a systematic way to combine them in the model-identification
step is proposed. Some potential problems associated with the procedure in practice and
their solutions are discussed. The efficacy of the procedure is demonstrated in a case
study in®ol®ing a multicomponent distillation column simulated in HYSYS.
Introduction
Ž
ventional feedback controller, and the indirect or model-
For many chemical processes, control is greatly limited by
process delays and measurement difficulties. For example,
measurement of a stream composition often requires an on-
line gas chromatograph, which introduces a large delay to the
system and suffers from a high failure ratio. Most product
property variables are similar in that they can only be ana-
lyzed in laboratories, requiring substantial processing time
and limiting the frequency of analysis. Moreover, there are
cases wherein key variables are measured in the downstream
but can be affected only in the upstream. Some control of the
variables of these types can be achieved by keeping the mea-
.
based approach in which one builds a stochastic system
model for the overall process and disturbances and designs a
state estimator and a state-feedback controller based on it.
These two approaches are depicted in Figure 2. Thus far, the
direct approach has seen some application in the industry,
but the indirect approach has received little attention from
the practitioners. This can be attributed to the fact that
building a model for the inferential controller design requires
substantially more information and effort than that for the
conventional control. Most importantly, such a model must
quantify how the measured outputs are related to the unmea-
sured outputs. In the context of fundamental modeling, this
Ž
surable process conditions such as temperatures, pressures,
.
flow rates stable and implementing feedforward controllers,
Ž .
but the variability can still be quite substantial given the of-
ten-significant feed variations and other unmeasured distur-
bances. For these cases, inferential control, a control strategy
based on the inferring of the composition or property vari-
means that the source s of disturbance must be identified
and their temporal behavior characterized as a stochastic
process. In the empirical modeling’s context, this requirement
calls for a multiple-input, multiple-output MIMO stochastic
system identification, something that is inherently difficult and
requires a high-level of expertise.
(
)
Ž . Ž
Ž . .
able s referred to as the ‘‘primary variable s ’’ based on the
Ž
measurement of other variables called ‘‘secondary
.
variables’’ , can be an attractive, viable alternative. A pictorial
Despite the practical barriers, there is a significant incen-
tive to further develop and tailor the indirect approach so
that it can be practiced. The recent acceptance by the indus-
representation of the candidate process for inferential con-
trol is given in Figure 1.
Ž
.
The idea of inferential control has existed for a long time,
but its industrial application is still at a nascent stage. Broadly
speaking, there are two approaches to inferential control; the
direct approach in which one builds a converter between the
try of model-predictive control MPC as the standard ad-
vanced control method is one reason. The indirect approach
fits naturally into the scheme of MPC, as a multistep predic-
tor of the primary variable that uses all the available mea-
surements in an optimal manner can be easily developed given
a stochastic system model. At a more general level, it pro-
vides a framework to predict and control the system behavior
optimally, without the need to categorize the variables explic-
itly into the categories of primary and secondary variables.
Multirate measurements and variable time delays that are
common in practice can be handled straightforwardly and
rigorously within this setting.
Ž
secondary measurements and the primary variables usually
.
through data regression of sorts and couples it with a con-
Correspondence concerning this article should be addressed to J. H. Lee.
Present addresses for: R. Amirthalingham, ProSys, Inc., Baton Rouge, LA; S. W.
Sung, Department of Chemical Engineering, Pohang University of Science and Tech-
nology, Pohang 790-784, Korea; J. H. Lee, School of Chemical Engineering, Georgia
Institute of Technology, Atlanta, GA 30332.
1974
October Vol. 46, No. 10
AIChE Journal

Figure 1. Candidate proces s for inferential control.
Since we pointed to the modeling difficulty as the main
the effect of unmeasured disturbances and system noises on
the future behavior of primary variables must be predicted
from secondary measurements, we need to capture accu-
rately into the model the autocorrelation and cross-correlation
information. This is unlike in the conventional feedback-con-
trol case, where the quick and frequent feedback of errors
for all output channels lessens the demand on the accuracy
of model prediction.
bottleneck in taking the indirect approach, it is natural to
seek a systematic framework for a building a requisite model
using plant data alone. The importance of such a framework
in order for a widespread use of the approach is supported
by the fact that almost all model-predictive controllers imple-
mented today are designed based on data from plant tests.
However, both data-collection and system identification
methodologies for building an inferential control model need
to be substantially different from those currently used. Since
Development of a data-based modeling framework for in-
ferential control is the main topic of this article. The ensuing
Figure 2. Two approaches to inferential control s ys tem des ign.
October Vol. 46, No. 10
AIChE Journal
1975

development stems from the issue of data requirement for
system identification. In order to capture the requisite corre-
lation information in the model, we need a large data set that
is rich in effects of unmeasured disturbances and noises. Data
from plant tests, the basis for model-building in current MPC
applications, do not meet this criterion, as data records from
plant tests are usually quite short and external conditions
the main ideas and compare the two. For a more detailed
Ž
discussion and references, see Amirthalingam and Lee 1997,
.
1999 .
Direct Approach. In the direct approach, an estimator for
the primary variable is constructed and coupled with a regu-
lar feedback controller, as shown in Figure 2a. Typically, a
static linear estimator is employed, which is either designed
using the revelant steady-state gain matrices or identified di-
Ž
.
such as disturbance variables are made artificially ‘‘quiet’’
Ž
.
during the test period in an effort to ensure good signal-to-
noise ratio. Hence, they seldom contain a sufficient amount
of disturbance effects needed to build a good correlation
model. For this, historical operation data are better suited, as
they come in a more abundant quantity and usually contain
disturbance effects that are more representative of real oper-
ation. On the other hand, one cannot build the entire control
model with historical operation data alone. Since the manip-
ulated inputs are likely to show very different correlation be-
havior once the inferential control loop is closed, we need a
good causal model between the manipulated inputs and the
outputs. Historical operation data rarely contain a sufficient
amount of uncorrelated input changes necessary to build a
good causal model. This information must be generated
through a plant test.
rectly through data regression. Least-squares LS or partial
Ž
.
least-squares PLS regression techniques have been most
commonly used for this purpose. Recently, the use of artifi-
Ž
.
cial neural networks ANNs for the estimator development
has become popular in industrial applications in which the
nonlinearity is a major issue. To account for the dynamics,
time-series-based regression has also been suggested.
The main attraction for the direct approach is the concep-
tual simplicity. On the other hand, this approach does not
render itself naturally to the implementation of advanced
control methods like MPC. To implement a predictive con-
troller using the direct approach, one has to build a separate
prediction model that uses the output of the estimator as the
main feedback variable. A potential problem with this is that,
when the loop is closed and the inputs are manipulated, the
Ž
.
These considerations motivated us to consider a two-step
approach, in which a stochastic system model is built first
using historical operation data, after which the model is used
as a prefilter in deterministic model identification conducted
with plant test data. We point out a number of issues that
arise within the two-step approach and propose some practi-
cal solutions. Finally, we apply the two-step approach to a
multicomponent distillation column in order to demonstrate
the feasibility of the approach and to bring to the surface
some additional issues that are specific to this type of pro-
cess.
estimator designed in prior to the loop closure may no longer
provide good estimates. This is because the closing of the
loop can significantly change the correlation behavior of the
Ž
.
input and output signals such as their power spectra .
(
)
Indirect Model-Based Approach. The root of this ap-
proach can be traced all the way to the general problem of
designing an output feedback controller for a stochastic sys-
tem. Theories for optimal prediction and control of linear
stochastic systems are well established, mostly through the
Ž
pioneering work of Kalman Kalman, 1960; Kalman and Bucy,
1961 . For articles that specifically address the topic of
.
The rest of the article is organized as follows. In the fol-
lowing section, two different inferential control approaches
are introduced, with a higher emphasis given to the model-
based approach. We also highlight some relevant develop-
ments in system identification and model-predictive control.
In the third section, we present the two-step identification
procedure. We first discuss a method for identifying a
stochastic system model using historical operation data and
the issues therein; we then present a method to obtain the
overall system model for the inferential controller design, uti-
lizing the developed stochastic system model and the data
from the plant test. Again, some potential issues are identi-
fied and addressed. We briefly discuss how the identified
model can be transformed into a prediction model that can
utilize measurements of variable delays and multiple or
nonuniform sampling rates. We also show how the state-space
model can be implemented as a predictor in the conventional
step-response-model-based MPC framework. In the fourth
section, we present the result of the case study. We conclude
the article in the fifth section with some final comments and
suggestions.
model-based inferential-control-system design, see the recent
Ž
.
review by Doyle 1998 . In the context of model-predictive
Ž
control, several researchers Lee and Morari, 1992; Lee et
.
al., 1992; Ricker, 1991; Muske and Rawlings, 1993 suggested
combining the Kalman filter with MPC formulated in state
space. Similar approaches have been reported in various ap-
Ž
plications’ context Banerjee et al., 1997; Cambell and Rawl-
.
ings, 1998; Yu et al., 1994 . There also have been numerous
articles reporting the success of the model-based approach in
Ž
batch processes Quintero-Marmol and Luyben, 1992; Lee
and Datta, 1994; Russell et al., 1998; Tatiraju and Souroush,
.
1997; Tatiraju et al., 1999 . These works employed the ex-
tended Kalman filters and other nonlinear observers to cope
with the strong nonlinear behavior of most batch processes.
Despite the academicians’ enthusiasm, the indirect ap-
proach has so far received little attention from the practition-
ers. The reason for this is that few articles shed light on the
requisite traits of the model for inferential control and even
fewer hint at any realistic way to obtain such a model when
the necessary fundamental knowledge is lacking. An excep-
Ž
tion is the recent work by Amirthalingam and Lee Amirtha-
.
lingam and Lee, 1997, 1999 , who discuss the data-based con-
struction of an inferential model for a pulp digester.
Preliminaries
Inferential control
Model-predicti©e control for inferential control problems
A convenient system representation for inferential con-
Inferential control techniques can be broadly classified into
direct and indirect approaches. In the following we present
1976
October Vol. 46, No. 10
AIChE Journal

Ž
.
troller design within the MPC framework is given by
x kq1 s Ax k q Bu k qw k
for the available data Ljung, 1987 . This method requires
some canonical parameterization of state-space matrices,
which require significant prior knowledge, especially in the
multiple output case. It also results in a nonconvex optimiza-
tion. An alternative is the method of subspace identification,
Ž
.
Ž .
Ž .
Ž .
y
k
Cc
®
®
k
cŽ .
cŽ .
s
x k q
Ž .
,
1
Ž .
Ž .
which consists of 1 the construction of data for two consec-
utive Kalman state vectors through an oblique projection in-
volving appropriate inputroutput matrices, and 2 the fitting
of model matrices through linear least squares Larimore,
Cs
^`_
y
sŽ .
k
k
sŽ .
^ ` _
^ ` _
Ž .
C
y k
Ž .
␯ k
Ž .
Ž
where u represents the manipulated input vector; yc and ys
Ž . Ž .
1990; Van Overschee and De Moor, 1994; Verhaegen, 1991;
.
the primary and secondary output vectors, and w k , ®c k ,
Ljung and McKelvey, 1996 . This method only involves nonit-
Ž .
and ®s k are white noise sequences. Measured disturbances
erative linear algebra operations and does not require any
prior knowledge. The subspace identification is a major de-
velopment in that it enables an easy, automatic construction
of a multivariable stochastic system model needed for the in-
ferential control system design. The prediction error method
and the subspace identification method can be used in a
complementary sense, for example, by using the result of sub-
space identification to obtain an appropriate parameteriza-
tion and starting values for the prediction error minimiza-
tion, which further improves the model fit.
can be included in ys without loss of generality.The preced-
ing model can be thought of as a superposition of two mod-
els, the deterministic model
xd kq1 s Axd k q Bu k
Ž .
Ž
.
Ž .
yd k sCxd k ,
2
Ž .
Ž .
Ž .
expressing how the manipulated inputs affect the outputs, and
the stochastic model
xs kq1 s Axs k qw k
Ž .
Two-Step Identification Procedure
Ž
.
Ž .
Why the two-step procedure?
ys k sCxs k q® k ,
3
Ž .
Ž .
Ž .
Ž .
Though the basic identification methods are already in
place, practicing them to obtain an inferential control model
is still not very straightforward. Most importantly, the issue
of data requirement deserves our attention. Traditionally,
data used for system identification are collected through a
plant test, designed to generate information about all the rel-
evant parts of system dynamics with an adequate signal-to-
noise ratio. This practice is justifiable when the sole aim of
the identification is to obtain a good causal model between
the manipulated inputs and the measured outputs. However,
modeling for inferential control has another dimension in that
Ž
which summarizes the statistics of the residual that is, the
.
total effect of all the unknown inputs and model error . Here
w and ® do not correspond to any physical input; rather, they
are just mathematical instruments used to express the auto-
correlation and cross-correlation of the residual. The super-
scripts d and s stand for ‘‘deterministic’’ and ‘‘stochastic’’
portions of the variables.
Ž
With such a model, an optimal state estimator the Kalman
.
filter and the optimal prediction equation for y can be easily
developed. Treatment of multiple sample rates, measure-
ments delays, etc., is straightforward within the Kalman fil-
tering. The future control input trajectory can be optimized
at each sample time based on the prediction equation, which
is easily derived from the Kalman filter equation, and the
typical quadratic criterion used in MPC. Of course, any rele-
vant constraints on the inputs and outputs can be entered
into the optimization. The details of the approach can be
Ž
we must also capture the correlations of the residuals due to
.
unmeasured disturbances, system noise, etc. . For this, plant
test data are seldom adequate, as they are typically short in
length and rarely contain sufficient disturbance information,
which is necessary to extract the correlation information that
is representative of normal operation.
In many cases, such information is better drawn from his-
torical operation data, which are available in more abundant
quantity and contain bountiful disturbance effects. The
premise here is that useful correlation information for pre-
dicting the behavior of future operations are contained in the
data collected from previous operations. This is reasonable in
view of the fact that autocorrelations and cross-correlations
seen in the outputs are mostly due to the process dynamics
and interactions. Barring any significant change in the pro-
Ž
found in several articles Lee et al., 1992; Muske and Rawl-
.
ings, 1993 .
The point we stress here is that the state-space model of
Eq. 1 represents a general and convenient vehicle for imple-
menting an inferential control strategy within the already
proven framework of model-predictive control. The main
challenge in practice then lies in obtaining a model of such
form, especially when the requisite knowledge for building a
fundamental process model is lacking.
Ž
cess dynamics due to large changes in the operating condi-
.
tion, etc. , a correlation model obtained through statistical
averaging of previous operation data should be relevant for
future operations.
Rele©ant identification methods
There are several ways to obtain a model of the form in
Eq. 1 using inputroutput data, but two methods are most rel-
evant in the present context. The first is the prediction error
method, where the model is first put in a predictor form and
parameters are estimated by minimizing the prediction error
This reasoning serves as the motivation for the two-step
approach, wherein we identify the stochastic part of the sys-
Ž
.
tem model given by Eq. 3 using historical data and the de-
Ž
.
terministic part given by Eq. 2 using plant test data. It still
remains to decide the best way to combine the two types of
AIChE Journal
October Vol. 46, No. 10
1977

identification. One can envision two different paths with re-
spect to this. In the first path, one identifies the deterministic
part of the model first using the plant test data. Then, using
the model, one subtracts the effect of the deterministic in-
puts from the historical operation data and fits a stochastic
model to the residual. In the second path, one would first
identify the stochastic part of the model using the normal
operation data and then follow this by identification of the
deterministic part using the test data. In this article, we pur-
sue the second path.
be correlated with the residuals. The standard N4SID method
gives a biased result in such a case.
In addition to these, there is the issue of data selection;
not all historical data are appropriate for the stated purpose.
Then, how do we select appropriate data from a huge
database? We discuss these issues next.
Choosing Slices of Data from Historical Data Base. In gen-
eral, historical data refer to huge piles of data stored in the
computer for days, months, and years. Most importantly, slices
of data corresponding to periods of normal operation should
be selected. Recall that we are trying to draw from these data
a correlation structure that can be used to predict the behav-
ior of future operations. Any abnormal situations that are
either one-time events or rare events should not be reflected
in the statistical model. In addition, periods corresponding to
unusually large upsets should be avoided, since they tend to
make correlations nonlinear. These large or abnormal upsets
may also be beyond the realm of automatic control, as they
often require special actions by the operator. The following
list can be used as a guideline in choosing the historical data.
Identification of stochastic system model using historical
operation data
As previously mentioned, the first of the two-step proce-
dure is to capture the autocorrelation and cross-correlation
information contained in the historical operation data into a
stochastic system model of Eq. 3. The method of subspace
identification is convenient for this, as it enables us to cap-
ture this information without a difficult numerical optimiza-
tion. It is also user-friendly in that it requires little input from
the user. For example, the N4SID method, one of the most
ⅷ
Select a period in which the system’s production rate did
not deviate significantly from the nominal value.
Ž
popular subspace identification algorithms Van Overschee
ⅷ
Select a period in which the reliability of the critical
.
and De Moor, 1994 , requires the user to give only an upper
Ž
measurements were high such as after a shutdown mainte-
bound of the system order and it automatically yields the sys-
tem model of the form
.
nance .
ⅷ
Select a period that includes no serious plant upset.
Select a period where the disturbances were likely to be
ⅷ
x kq1 s Ax k q Bu k q K⑀ k
Ž .
Ž
.
Ž .
Ž .
representative of most other times.
ⅷ
Select a period where input manipulations were not un-
y k sCx k q⑀ k .
4
Ž .
Ž .
Ž .
Ž .
usually high.
Nonstationary Disturbances. Chemical processes com-
monly experience disturbances that make the process vari-
From the preceding, we retain only
xs kq1 s Axs k q K⑀ k
Ž .
Ž
.
ables drift or jump . This violates the basic requirement of
the N4SID method that the stochastic part of the system be
Ž
.
Ž .
Ž
.
stationary Van Overschee and De Moor, 1994 . When the
data used do not meet this requirement, these algorithms can
behave badly, even creating artificial, unstable modes.
Assuming that the changes are not large, this nonstation-
ys k sCxs k q⑀ k ,
5
Ž .
Ž .
Ž .
Ž .
since the deterministic part of the model is hardly adequate
given the usual lack of systematic input excitation in histori-
cal operation data. The subspace identification method re-
quires a relatively large amount of data, a major drawback
for regular system identification, but not so for our purpose,
since historical operation data are usually available in bounti-
ful quantity and at no cost.
Even though the methodology is well established in the lit-
erature, we found that applying it to chemical processes gives
rise to some additional issues not covered in the literature,
and thus requires some tailoring and care. These issues in-
clude:
Ž .
ary behavior is well approximated by including integrator s
Ž
to the usual time-invariant stochastic system model Morari
.
and Stephanopoulos, 1980 . Within this context, the presence
of integrators is what makes the system nonstationary. A
simple way to resolve this conflict is to ‘‘difference’’ the in-
Ž
putroutput data to express the data in terms of their incre-
.
mental change from the previous samples which effectively
removes the integrator from the system. In other words, dif-
ferencing of the data before applying the subspace identifica-
tion algorithm makes the assumption of stationarity more re-
alistic and hence makes the algorithm better behaved. In fact,
this feature of data-differencing can be found in the impulse-
response identification modules of some commercial MPC
ⅷ
The N4SID method requires that the stochastic part of
the system be stationary. On the other hand, most chemical
processes show nonstationary characteristics, most notably the
process variables drifting or experiencing jumps.
Ž
.
packages Cutler and Yocum, 1991 and is referred to as ‘‘de-
trending.’’ Along with the differencing, some low-pass filter-
ing may be necessary, as the differencing tends to overaccen-
tuate the portions of the data pertaining to the system’s
high-frequency dynamics.
ⅷ
Most inferential control problems involve multirate mea-
surements, with the primary variables sampled at a substan-
tially slower rate than the secondary variables. Measurements
can also be made with substantial delays, which may vary with
time.
The obtained model will be in the form of
ⅷ
Most historical data contain input adjustments, made by
x kq1 s Ax k q B⌬u k q K⑀ k
Ž .
Ž
.
Ž .
Ž .
operators or existing supervisory loops in response to some
upsets or major disturbances. This means input data would
⌬ y k sCx k q⑀ k .
6
Ž .
Ž .
Ž .
Ž .
1978
October Vol. 46, No. 10
AIChE Journal

Note that the inputs and outputs of the system model are
incremental changes in the process inputs and outputs at each
sample time instead of the usual deviation values. We will
discuss later how we can use this model to form a predictor
for MPC.
Let us support the preceding assertion by presenting a sim-
ple example. An output data set consisting of 1000 samples
was generated from the following stochastic state-space sys-
tem:
use ‘‘lifting’’ which is a technique to transform a multirate
Ž
sampled-data system into a single-rate system Friedland,
.
1960; Kranc, 1957; Lie et al., 1999 .
Handling Closed-Loop Historical Data. Historical data of-
ten include changes in the manipulated inputs that are made
by the operator in response to some observed output behav-
ior. If such manipulations are recorded and the records are
available in the historical data, it should be possible to re-
move their effect from the data. Because these moves are
based on output behavior, however, they are correlated with
the output data. With such correlated data, the conventional
subspace identification method gives a biased result.
If input adjustments are minor and infrequent, they can be
ignored without any substantial problem. If input adjust-
ments are significant, there are extensions of the subspace
identification method that can handle closed-loop data. For
x kq1 s0.8 x k q0.4␩ k
Ž .
Ž
.
Ž .
y k s0.2 x k q␩ k ,
7
Ž .
Ž .
Ž .
Ž .
Ž .
where ␩ k was obtained by integrating a white-noise se-
quence of variance 0.01. The application of the N4SID algo-
rithm with y data resulted in the following stochastic system
model:
Ž
.
example, Ljung and McKelvey 1996 propose a method that
Ž
uses multistep predictions of the output based on a high-
.
order ARX model instead of actual output data in defining
the state.
1.0000 y0.0136
0.3918
y7.0867
x kq1 s
x k q
⑀ k
Ž
.
Ž .
Ž .
In order to further demonstrate this point, let us consider
the simple system
0.0070 1.1836
w
x
Ž .
y k s 1.5144 y0.0037 x k q⑀ k .
8
Ž .
Ž .
Ž .
x kq1 s0.8 x k q0.7u k q0.4⑀ k
Ž .
Ž
.
Ž .
Ž .
Ž
.
Note the presence of unstable poles 1.005 and 1.1836 , which
resulted from the integrating nature of the external input ␩.
We also applied the algorithm with ⌬ y data, that is, after
differencing the data. The resulting model was
y k s0.2 x k q⑀ k ,
10
Ž .
Ž .
Ž .
Ž .
Ž .
where ⑀ k is a white-noise sequence of variance 0.01. We
Ž
assumed that normal operation data consisting of 1000 sam-
x kq1 s0.7720 x k q0.4640⑀ k
Ž .
.
Ž
.
Ž .
ples are collected under the following intermittent closed-
loop control strategy:
⌬ y k s0.2 x k q⑀ k ,
9
Ž .
Ž .
Ž .
Ž .
u k sy0.5= y k
for ks10, 20,30, etc.
Ž .
Ž .
Ž .
where ⑀ k is a white-noise sequence. We can see the clear
improvement in the model fit.
u k su ky1 otherwise.
11
Ž .
Ž
.
Ž .
Measurement Delays and Multiple Sample Rates. It is best
that known delays be taken out before the identification al-
gorithm is applied. Delays increase the system order, which
makes system identification more difficult. This is especially
true if the length of the delay is significant in relation to the
sample time. The delay taken out can always be added back
after the model is identified.
Application of the regular N4SID algorithm yielded the fol-
lowing stochastic system model:
x kq1 s0.7751 x k y1.2033u k y0.925⑀ k
Ž .
Ž
.
Ž .
Ž .
y k s0.2 x k q⑀ k .
12
Ž .
Ž .
Ž .
Ž .
To eliminate delays due to the sampling and analysis, data
Ž
can be marked with the time of sampling rather than the
Clearly the identified model is extremely biased, having a gain
of wrong sign. With the algorithm modified to handle
closed-loop data, we obtained
.
time when the analysis is completed . This eliminates mea-
surement delays from the data regardless of whether the de-
lays are constant or time-varying. Any known process delay
Ž
.
can also be taken out by shifting the output input data by
that amount. If the length of the delay is not exactly known
or time-varying, the data can still be shifted by a lower bound.
In addition, data for the primary variables can be available
at much slower rates than those for the secondary variables.
The model, on the other hand, should be developed at the
faster of the two, in order to realize the benefits of the sec-
ondary measurements. One may be able to increase the sam-
pling rate just for the purpose of system identification in some
cases, but this may not always be possible. If the sampling
rate for the primary variable cannot be increased, interpola-
tion can be used to fill in the missing points. From our expe-
rience, simple linear interpolation works adequately given
that the measurements are not very noisy. An alternative is to
x kq1 s0.6635x k q0.7198u k q0.3333⑀ k
Ž .
Ž
.
Ž .
Ž .
y k s0.2 x k q⑀ k ,
13
Ž .
Ž .
Ž .
Ž .
which is a clear improvement.
Closed-loop identification is still at an evolving stage and
many issues have not been resolved to a satisfactory level at
Ž
this point. Several recent articles Ljung and Forssell, 1999;
.
Forssell and Ljung, 1999 address the issue of closed identifi-
cation and its importance. Undoubtedly, better closed-loop
identification algorithms for stochastic systems will come
along in the future, and these algorithms can be readily used
in the overall model-building strategy.
AIChE Journal
October Vol. 46, No. 10
1979

Identification of the deterministic part using plant test data
Ž
.
This is an output error OE structure. In prediction error
Ž .
minimization, Gf q can be estimated by minimizing the out-
The next step after identifying the correlation structure is
to generate a data set that is rich in input excitation by con-
ducting a plant test. The usual rules for plant testing apply
here as well. However, there is an issue that arises because
we difference our data before model fitting. We will also dis-
cuss what role, if any, the stochastic part of the system model
Ž .
Ž . Ž .
put error ⌬ yf k yGf q ⌬u k . The data from a plant test
are to be used for this purpose.
Ž .
The preceding requires the inverse of H q to be stable.
Recall that the stochastic system model is obtained through
the N4SID method, which in effect identifies the model in
Ž
.
the form of a Kalman filter the innovation form . It can be
Ž .
Ž
.
which is already fixed should play in obtaining the deter-
Ž
shown that H q generated from a Kalman filter or any sta-
ministic part.
.
ble observer of a stochastic system according to Eq. 16 is
always stably invertible. If H q happens to contain unstable
zeros due to some identification error, one can always design
the Kalman filter based on the model to readjust H q , which
Test Signal Design. Since we have obtained the stochastic
Ž .
Ž
.
model in terms of the differenced variable ⌬ y , it is conve-
nient also to identify the deterministic model in terms of the
Ž .
Ž
.
differenced variables ⌬ y and ⌬u . This should also be done
in view of the fact that process means do drift even during a
plant test, which often lasts for several days. While the use of
gives the desired result.
The OE identification can be conducted with the regular
polynomial transfer function structure or the finite impulse
response FIR structure. The use of a MIMO structure is
not required since correlation information necessary for the
inferential control design has already been captured. Once
Ž
.
random white-noise signals or pseudorandom binary-se-
Ž
.
Ž
.
quence PRBS signals is recommended in the literature,
these signals, once differenced, contain zero power at the zero
frequency. The implication is poor accuracy for the low-
frequency part of the model.
Ž .
Ž .
Gf q is identified, G q can be extracted out, using the rela-
tionship
Ž .
To combat this problem, we can generate ⌬u k as a PRBS
sequence. Integrating the sequence then gives the desired re-
sult. The generated sequence must be scaled to fit within the
allowed magnitude bounds for the test signal. Different real-
izations of PRBSs give different maximum deviations from
the mean, once integrated. Here, one should choose the real-
ization that yields the lowest maximum deviation in order to
maximize the power of ⌬u.
Another common problem encountered in identifying
strongly interactive systems such as a high-purity distillation
column is the problem of ill-conditioning. Directional input
design, which requires information on the directionality of
G q s H q G q .
18
Ž .
Ž . f Ž .
Ž .
We can then perform a minimal realization on the overall
model of Eq. 15 to obtain a state-space model of the form
x kq1 s Ax k q B⌬u k q K⑀ k
Ž .
Ž
.
Ž .
Ž .
⌬ y k sCx k q⑀ k .
19
Ž .
Ž .
Ž .
Ž .
Why Output Error Identification? One might ask why we
advocate OE identification here. After all, it seems more
convenient to apply the subspace identification just one more
time. One could also consider fitting an ARMAX structure
afterward to further improve the result. A key dilemma here
is what to do with the stochastic part of the system model. It
is tempting to simply ‘‘throw away’’ the stochastic part and
retain only the deterministic part. However, there are some
potential risks associated with this. To demonstrate the risk,
we have applied N4SID, ARMAX, and OE identification
methods to the following system:
Ž
the system, has been suggested as a remedy Andersen et al.,
1989; Jacobsen and Skogestad, 1994; Li and Lee, 1995; Coo-
.
ley and Lee, 2000 .
Identification of the Deterministic Part Using Output Error
Identification. After applying the N4SID method to histori-
cal operation data, we are left with the stochastic system
model in the form of
xs kq1 s Axs k q K⑀ k
Ž .
Ž
.
Ž .
⌬ ys k sCxs k q⑀ k .
14
Ž .
Ž .
Ž .
Ž .
1
2
0.5qy
1y0.4qy1 q0.04qy
The overall system model we are seeking can be reexpressed
in an inputroutput operator form as
y k s
Ž .
u k q
Ž .
⑀ k ,
Ž .
2
2
1y0.98qy1
1y0.98qy1
Ž
.
Ž
.
20
Ž
.
⌬ y k sG q ⌬u k q H q ⑀ k ,
15
Ž .
Ž .
Ž . Ž .
Ž . Ž .
Ž .
Ž .
Ž .
where G q is the transfer matrix between u k and y k .
Ž .
where ⑀ k is a white-noise sequence of variance 0.1 and the
test signal is a uniformly distributed random signal between
"0.5. The results we obtained with a data set consisting of
500 samples were as follows:
Ž .
The noise model H q is already set by the stochastic system
model of Eq. 14 according to
y
H q sC qIy A 1 Kq I.
16
ARMAX
Ž .
Ž
.
Ž
.
1Ž .
Premultiplying the Eq. 15 with Hy q gives
2
0.3535
1y1.0024qy
1q0.8104qy1 q0.59qy
1y1.0024qy1
y k s
u k q
Ž .
⑀ k .
Ž .
Ž .
1
1
1
Hy q ⌬ y k s Hy q G q ⌬u k q⑀ k . 17
Ž . Ž .
Ž . Ž . Ž .
Ž .
Ž
.
21
Ž
.
^
`
_ ^ `
Gf q
_
⌬ y
k
f Ž .
Ž .
1980
October Vol. 46, No. 10
AIChE Journal

N4SID
analysis delay can vary from one sample to the next. In this
case, the number of delayed outputs augmented should cor-
respond to the maximum delay. This way, one always has the
delayed measurement appearing as a part of the state. One
1
0.7071
1y1.0033qy
1y1.2294qy
y k s
Ž .
u k q
⑀ k . 22
Ž .
Ž .
Ž
.
.
1y1.0033qy1
1
Ž .
simply has to choose ⌶ k at each time to pick out the appro-
priate output. Even with the augmentation, the form of the
equation will still remain the same as Eq. 25.
The preceding is in the standard state-space form and is
ready for the design of the Kalman filter. Also, one can de-
OE
3.877
y k s
Ž .
u k q⑀ k .
23
Ž .
Ž .
Ž
1y0.9842qy
1
Ž .
velop a multistep prediction equation for y k based on the
Kalman filter equation that serves as the basis for MPC cal-
culation. For details, the readers are referred to Lee et al.
Notice that both the ARMAX identification and the N4SID
method created an unstable pole in the deterministic model.
This is not a problem of any practical relevance if the
stochastic system model is used together with the determinis-
tic model in the predictive controller design. In that case, the
unstable pole gets stabilized by the same pole in the noise
model in forming the optimal predictor. However, if the
stochastic part of the system model is abandoned and re-
Ž
.
Ž
.
1992 , Amirthalingam and Lee 1997, 1999 .
Implementation as a Step-Response-Based Predictor. Most
commercial MPC packages are based on the step-response
representation of system dynamics. Hence, it is of interest to
find a simple way to implement the model that results after
the identification into these packages. For this, we show that
Eq. 19 can be converted directly into a step-response predic-
tor in the form of
Ž
placed with something else such as the one obtained in the
.
previous step , the optimal predictive controller becomes in-
Y kq1 s MY k qSu⌬u k q K⑀ k
26
27
Ž
.
Ž .
Ž .
Ž .
Ž
Ž
.
.
ternally unstable. The same problem cannot occur in the OE
identification. Although this is a simplistic example designed
for the sole purpose of demonstrating the point, we experi-
enced similar problems in identifying higher-order systems,
including the one we discuss later in the case study.
y k s NY k q⑀ k ,
Ž .
Ž .
Ž .
Ž .
where Y k is the state vector of the step-response model
with the meaning of
In summary, fixing the noise model ahead of identification
Ž
means that we are a priori fixing the observer or the state
.
T
estimator that will be used in the real application. It is then
⌬
Y k s yT k , ... , ... , yT kqny1
for ⌬u kqi s0
Ž .
Ž .
Ž .
Ž
.
sensible to minimize the prediction error under that very ob-
server. And this objective leads to OE identification.
and ⑀ kqi s0, iG0, 28
Ž .
Ž
.
Turning the model into a multistep predictor for MPC
M represents the usual shift operator for the step-response
Ž .
Ž .
Y k , and
Ž
.
model Lee et al., 1994 , N is the matrix extracting y k from
Implementation in State Space. The two-step approach
yields a combined model of the form given in Eq. 19. Since
Ž
.
our goal is to predict the behavior of y rather than ⌬ y , it is
CBu
convenient to use the following augmented form of Eq. 19:
CK q I
Ý1 CAiK q I
Ý1is0CAiBu
Ž
.
is0
.
.
.
.
.
.
x kq1
y kq1
Ž
x k
Ž
.
.
Ž .
A
CA
0
I
Su s
,
Ks
,
29
s
Ž
.
y k
Ž .
.
.
.
.
.
.
B
CB
0
I
k
Ck
ny1
Ý
CAiBu
Ýny1CAiK q I
q
⌬u k q
⑀ kq1 q
⑀ k
Ž .
Ž .
Ž
.
Ž
.
is0
is0
where n should be chosen large enough so that CAnBu f0
and CAnK f0. Note that this model can be implemented di-
x k
Ž .
ym k sw0 H k x
Ž . Ž .
,
24
Ž
.
y k
Ž .
Ž .
Ž .
Ž .
rectly as a state estimator, since ⑀ k s y k y NY k .
The prediction equation can also be developed easily based
on the preceding model and takes the form of
Ž .
where ym k represents the subset of y measured at the kth
Ž .
Ž .
sample time and H k is the matrix that extracts ym k from
Ž .
y k . For convenience, let us denote the preceding as
Y kq1Nk s M Y k q S ⑀ w y k y NY k xq S U⌬U k .
Ž
.
Ž .
Ž .
Ž .
Ž .
p
z kq1 s⌽z k q⌫⌬ u k qw k
Ž .
Ž
.
Ž .
Ž .
30
Ž
.
ym k s⌶ k z k .
25
Ž .
Ž .
Ž . Ž .
Ž
.
Basically, this equation adds to Y kq1 the effect of future
w
Ž
.
Ž
.
x
Measurement delays not included in Eq. 19 can be incorpo-
rated into the preceding model using the standard procedure,
input moves ⌬u kq1 , ⌬u kq2 , etc. . The future values of
⑀ are assumed to be zero since it is a zero-mean variable.
The symbol Mp is used instead of M for the shift operator,
since the prediction horizon p may be different from n.
Ž
that is, by augmenting the state with delayed outputs see
˚
.
Astrom and Wittenmark, 1997, pp. 38᎐41 . In some cases,
¨
AIChE Journal
October Vol. 46, No. 10
1981

Table 1. The Nominal Operating Condition of the
Depropanizer
measurements were assumed to be available every 30 min,
with a 30-min delay.
There are several reasons why the multicomponent distilla-
tion column makes for a good case study. First, unlike binary
columns, a direct one-to-one relationship between tray tem-
peratures and tray compositions is missing. In the presence
of feed disturbances, the well-known industrial practice of
regulating selected tray temperatures, no matter how well
placed the thermocouples are, may not deliver a satisfactory
result. This point is illustrated in Figure 3, which shows the
time profiles of the endpoint compositions of the lean key
components when the seventh and twenty-fourth trays are
controlled. When a step disturbance occurs in the feed com-
position, despite the fact that the temperatures are success-
fully returned to their respective setpoints, the endpoint com-
positions show significant offsets. There are some case-
specific solutions in the literature, such as those involving dif-
ferentials and double differentials of certain tray tempera-
Feed Compositions
Distillate Compositions
Ethane
Propene
Propane
i-Butane
cis2-Butene
0.0300 Ethane
0.4000 Propene
0.1500 Propane
0.1500 i-Butane
0.0517
0.6863
0.2562
0.0050
0.2700 cis2-Butene
0.0007
Bottoms Compositions
Flow Rates
Ethane
0.0000 Feed
3,000 kmolrh
1,740 kmolrh
1,260 kmolrh
6,105 kmolrh
Propene
Propane
i-Butane
cis2-Butene
0.0048 Distillate
0.0033 Bottoms
0.3501 Reflux
0.6417 Boilup
7,622 kmolrh
7
Ž
.
Reboiler dutyᎏ8.516=10 kJrh
Temperatures
90ЊC
Others
Feed
Column pressure
2,600 kPa
Distillate
Bottoms
59.27ЊC Feed enthalpy
125.1ЊC
y1.259=108 kJrh
Ž
.
tures Yu and Luyben, 1984 , but these are merely point solu-
tions that do not generalize.
Cas e Study
Description of the case study
Data generation
In this section, we present the result of applying the two-
step approach to a multicomponent distillation column, a de-
Historical Operation Data. In generating the historical op-
eration data, we assumed the maximum feed flow-rate distur-
bance to be "150 kmolrh. To simulate the feed-composition
variations, we added an external stream of propylene with a
maximum flow rate of 10 kmolrh to the feed stream. This
external stream was added to the feed stream before the feed
flow-rate control valve. Feed-enthalpy variations introduced
were "1=105 kJrh, which changed the vapor fractions of
the feed stream from 0 to 0.0046 the nominal value of 0.0016 .
These values for the disturbances were chosen after studying
the effect of the individual disturbances on the controlled
variables. The bounds were chosen such that all disturbances
contribute about the same to the output variations. All three
disturbances were simulated as integrated white-noise se-
quences.
Ž
.
propanizer, studied by Yiu et al. 1990 . The column has 29
trays and is designed to fractionate a feed stream consisting
of 3% ethane, 40% propene, 15% propane, 15% i-butane,
Ž
.
and 27% cis-2-butene on molar basis into a distillate stream
containing the lighter components, ethane, propene, and
propane, and a bottom stream containing the heavier compo-
nents, i-butane and cis-2-butene. We simulated this column
using HYSYS, a dynamic simulation software. The input and
output stream compositions and column conditions are sum-
marized in Table 1.
Ž
.
Given the goal of reducing the heavier components in the
distillate stream and lighter components in the bottom stream,
we defined our controlled variables to be the concentration
Ž
.
of the light-key component propane in the bottom and the
In order to simulate the operation more realistically, a pro-
Ž
.
concentration of the heavy-key component i-butane in the
distillate. We assumed that three types of unmeasured distur-
bances can occur at the feed, namely, feed-composition varia-
tions, feed-enthalpy variations, and feed flow-rate variations.
The inferential control was aimed at reducing the variability
in the controlled variables in the presence of these distur-
bances. The manipulated variables were the reflux flow rate
and the reboiler duty. The secondary measurements for in-
ferential control were the temperature measurements at the
first, fifth, tenth, fifteenth, twentieth, twenty-fifth, and
twenty-ninth trays.
w Ž .
Ž .x
portional controller u k sy60,000 y k was assumed to
Ž
have been in place prior to the inferential controller imple-
.
mentation , which controls the composition of i-butane in the
distillate by manipulating the reflux flow rate every 30 min.
Inputroutput data were collected for a duration of 3,000 min,
with a sample time of 3 min. The data matrix used in the
closed-loop subspace identification algorithm was a 1000=11
matrix containing 7 temperature measurements, 2 composi-
Ž
tion measurements, and 2 inputs the reflux flow rate and the
.
reboiler duty .
Plant Test Data Generation. In generating the plant test
data, we faced an issue that is specific to ill-conditioned sys-
tems like the high-purity distillation column. The process gain
for such a system depends strongly on the input direction. In
HYSYS and MATLAB were interfaced using the digital
Ž
.
data exchange DDE capability of the two software pro-
grams. Hence, all the mathematical calculations except for
the simulation of the column were carried out in MATLAB.
At every sample time, the values of the output variables of
the depropanizer column were transferred to MATLAB from
HYSYS and the calculated input moves and the stored dis-
turbance sequences were transferred back to HYSYS. The
basic sample interval was chosen as 3 min. This means all the
temperature measurements were read in and new control
moves were implemented every 3 min. The end composition
w
x Ž
the distillation column, the direction of 1 y1 that is, in-
creasing one of the manipulated variables and decreasing the
.
other induces a large response in the compositions, while
w
x Ž
the direction of 1 1 that is, a simultaneous increase or de-
crease of both manipulated variables induces a much smaller
.
response. If both inputs are perturbed independently with an
equal amount of energy, the output response will be predom-
1982
October Vol. 46, No. 10
AIChE Journal

the presence of the proportional control loop, we pro-
grammed the modified N4SID algorithm suggested by Ljung
Ž
.
and McKelvey 1996 in MATLAB and used it to obtain a
stochastic system model. In the identification, the maximum
order we gave was 20, and the model order of 10 was chosen
after the singular values of the predicted output block-Hankel
Ž
.
matrix were inspected Van Overschee and De Moor, 1994 .
A model in the form of Eq. 6 was then obtained, and the B
matrix was discarded, since the level of input excitation in
these data was judged too low. The model was then scaled
back to the original units using the inputroutput scaling fac-
tors.
Determinstic model identification
In the deterministic model identification, our strategy in-
volved filtering the differenced output data with the stochas-
tic model, followed by an OE identification for each output.
In the present case study, we needed to identify nine MISO
models using the plant test data. First, the input and output
data were differenced. Then, the inverse of the noise model
was computed and the differenced output data were filtered
by it. The filtered output vectors and input vectors were scaled
Figure 3. Behavior of compos itions under the control of
tray temperatures for a s tep change in feed
compos ition of I0.05 in propylene and + 0.05
in cis obutene.
Ž .
similarly as before. We then identified Gf q using the scaled
Ž .
Ž .
data for ⌬u k and ⌬ yf k . For this, we used the following
Ž
.
MISO OE structure for each output:
Steady-state feed composition: ethane 0.03 , propylene
Ž
.
Ž
.
Ž
.
Ž
.
0.40 , propane 0.15 , isobutane 0.15 , cisobutene 0.27 .
⌬ yf
k
iŽ .
4
4
b1 qy1 qb2 qy2 qb3 qy3 qb4 qy
inantly along the high gain direction, causing a poor fit of the
dynamics along the low gain direction.
i, 1
i, 1
i, 1
i, 1
s
q
⌬u k
1Ž .
1qa1 qy1 qa2 qy2
i
i
For the depropanizer column, the condition number was
about 25. Although this is not a severely ill-conditioned sys-
tem, by experience, we realized that the process’s gain direc-
tionality must be taken into account in the test signal design
in order to identify a good causal model needed for two-point
composition control. Since we knew that the large gain direc-
b1 qy1 qb2 qy2 qb3 qy3 qb4 qy
i, 2
i, 2
i, 2
i, 2
⌬u2
k
Ž .
1qa1 qy1 qa2 qy2
i
i
qe k , is1,... , 9. 31
Ž .
iŽ .
w
x
w
x
tion is 1 y1 and the small gain direction 1 1 , we first gen-
w
x
erated an integrated white-noise sequence in the 1 y1 di-
rection, with the constraints on the reflux rate and the re-
boiler duty set as "45 kmolrh and "1=105 kJrh. After gen-
erating another integrated white-noise input sequence in the
This structure was chosen by comparing the mean square of
the output error for several different structures. After scaling
back Gf q to correspond to the original variables, we deter-
Ž .
Ž .
mined the deterministic transfer matrix G q by premultiply-
w
x
1 1 direction, with the constraints relaxed to "225 kmolrh
Ž .
Ž . Ž
.
Ž .
ing Gf q with H q as in Eq. 15 . The transfer matrix G q
was then realized as a state-space system and put together
with the stochastic system to form the overall system model.
and "2.5=106 kJrh, we superimposed the two and used that
as the test input. The disturbances introduced during the
plant test were similar to those introduced during the histori-
cal data generation. The plant test was conducted for 600
min, and a data set consisting of 200 data points was ob-
tained.
Estimation results
Using the identified model for the development of a multi-
rate Kalman filter was already discussed in the third section.
After augmenting the state with the outputs and 10 addi-
Stochastic model identification
Ž
tional delayed outputs to account for the 30-min delays in
.
The historical operation data matrix included 1000 samples
of the 7 temperatures, the 2 compositions, the reflux flow
rate, and the reboiler duty. After differencing the data to
the primary measurements , we designed a multirate Kalman
filter and tested its estimation performance through simula-
tion. Since the primary measurements were assumed to be
available every 30 min and the temperature measurements
every 3 min, a periodic Kalman filter with the period of 10
sample steps was used. We compared the estimation perfor-
mance of the multirate Kalman filter against that of the sin-
gle-rate Kalman filter which was only the temperature mea-
surements and that which uses only the delayed composition
Ž .
Ž .
create a data matrix containing ⌬ y k and ⌬u k , we scaled
the data by dividing the elements of each differenced input
and output vector by its standard deviation. The scaling is
critical since the choice of state space and model fitting in
the N4SID algorithm is based on the total prediction error
for the given data. Because these data had been generated in
AIChE Journal
October Vol. 46, No. 10
1983

measurements. For the latter, we identified a separate state-
Ž
.
space model of sample time of 30 min , which had the com-
position variables as the only outputs. The estimation was
performed only once every 30 min in this case.
First, the performances of the three estimators were simu-
lated when deterministic step disturbances occurred in the
feed enthalpy, feed composition, and feed flow rates. Estima-
tion results were also generated by varying all the three dis-
Ž
turbances in a stochastic fashion that is, as independent in-
.
tegrated white-noise sequences . The results are shown in
Figures 4᎐7 and Tables 2᎐5.
The simulation results show that reasonable predictions of
the composition variables can be made using the temperature
measurements alone, but the predictions can be improved
further by using the delayed infrequent composition mea-
surements in addition to the temperature measurements. Us-
ing the delayed composition measurements alone, on the
other hand, gives a poor estimation performance.
Closed-loop control results
The next step was to design a model-predictive controller
and simulate its closed-loop performance. At each sample
time, the controller used the result from the Kalman filter to
build a prediction of future outputs and calculated a se-
quence of optimal input moves based on the quadratic crite-
Figure 5. Evaluation of the es timators ’ performance for
(
a s tep change in the compos ition external
)
propylene s tream flow rate + 10 kmol/h .
Ž
rion Lee et al., 1992; Lee et al., 1994; Muske and Rawlings,
Ž
.
where ⌳ y was chosen as diag 1, 1, 0, 0, 0, 0, 0, 0, 0 , given the
equal importance of the compositions and the absence of
control requirements for the temperatures. Since the system’s
gain exhibited a strong dependence on the input direction,
we tuned the input weighting matrix with the structure
.
1993 of
p
min
yT kqiNk ⌳ y kqiNk
Ž
.
Ž
.
Ý
y
Ž
.
Ž
.
⌬ u k , . . . ,⌬ u k q m y1
is1
m y1
␭
0
q
⌬uT kq j ⌳ ⌬u kq j , 32
Ž . Ž . Ž .
Ý
u1
u
1
1
1
y1
1
1
1
y1
⌳u s
.
33
Ž
.
js 0
0
␭
u2
Figure 4. Evaluation of the es timators ’performances for
5
(
)
a s tep change 1=10 kJ /h in the feed en-
thalpy with a change in vapor fraction from
0.016 to 0.019.
Figure 6. Evaluation of the es timators ’performances for
(
a s tep change in the feed flow rate + 150
)
kmol/h .
1984
October Vol. 46, No. 10
AIChE Journal

Table 5. Evaluation of the Estimators’ Performance for
Stochastic Changes in the Feed Composition, Feed
Temperature, and Feed Flow Rate
Mean Square Error in Estimation
y1
y2
With delayed comp. meas.
With temp. meas. only
With delayed comp. meas.
and temp. meas.
1.4606ey008
2.4351ey009
4.8160ey010
3.6739ey007
7.5129ey010
7.2304ey011
Ž
.
see Table 6 . This is logical, as the choice of input weight
needed to achieve a certain speed of response directly de-
pends on the size of the gain. The values we used for the
tuning parameters are summarized in Table 6.
Ž
We simulated the same three scenarios that is, with the
.
three different measurement combinations as before. In the
cases where the temperature measurements were used either
by themselves or together with the delayed composition mea-
surements, the input calculations were performed and the
optimal input moves were implemented every 3 min. When
only the composition measurements were used, this was done
once every 30 min. As before, we first simulated the closed-
loop responses for various deterministic changes in the feed
and then introduced stochastic disturbances. The simulation
results for a step change in the feed temperature, feed com-
position, and feed flow rate, and for stochastic changes in all
three disturbances, are shown in Figures 8᎐11 and Tables
7᎐10 .
Figure 7. Evaluation of the es timators ’ performance for
s tochas tic changes in the feed compos ition,
feed temperature, and feed flow rate.
Table 2. Evaluation of the Estimators’ Performance for
5
(
)
a Step Change 1=10 kJrh in the Feed Enthalpy with
a Change in Vapor Fraction from 0.016 to 0.019
Mean Square Error in Estimation
Although not shown here, we made several painstaking at-
tempts to build an inferential control model using the plant
test data alone. In general, the combined deterministicrsto-
chastic model, which we obtained from applying the subspace
identification method to the plant test data, showed poor es-
timation and control performance. Even when we gave a very
high signal-to-noise ratio by eliminating most of the distur-
bances during the plant test, the results did not improve. In
this case, even though the deterministic part of the obtained
model was fairly accurate, the captured correlation between
the temperatures and the compositions, critical for inferen-
tial control, was poor since such data contained little distur-
bance information. On the other hand, when we increased
the level of disturbances during the plant test, identification
of the deterministic model became a problem. The estima-
tion performance improved, but the control performance be-
came extremely poor, often showing instability. Only by com-
bining the historical operation data and the plant test data in
the proposed manner could we derive an inferential control
model that delivered satisfactory estimation and control re-
sults.
y1
y2
With delayed comp. meas.
With temp. meas. only
With delayed comp. meas.
and temp. meas.
7.7741ey008
1.8137ey009
2.6687ey010
1.3020ey006
3.4625ey010
2.0517ey011
Table 3. Evaluation of the Estimators’ Performance for
a Step Change in the Composition
Mean Square Error in Estimation
y1
y2
With delayed comp. meas.
With temp. meas. only
With delayed comp. meas.
and temp. meas.
1.6038ey007
9.5158ey011
2.6637ey006
1.0132e-010
9.9340ey011
2.0562ey011
Note: External propylene stream flow rate q10 kmolrh.
Table 4. Evaluation of the Estimators’ Performance for
a Step Change in the Feed Flow Rate
Mean Square Error in Estimation
y1
y2
Conclus ion
With delayed comp. meas.
With temp. meas. only
With delayed comp. meas.
and temp. meas.
1.7402ey007
9.8302ey008
6.9290ey009
4.3230ey006
9.6829ey009
1.0766ey009
In this article, we examined the feasibility of using plant
data to build a stochastic system model that can be used to
Note: q150 kmolrh.
Table 6. MPC Tuning Parameters
Control horizon: ms3
Prediction horizon: ps10
This way, we could penalize input movements along the low-
gain direction less than those along the high-gain direction
Ž
.
u
2
Ž
.
Input weight:
␭
, ␭ s 0.02, 0.1
u
1
AIChE Journal
October Vol. 46, No. 10
1985

Figure 8. Evaluation of the controller performance for a
Figure 10. Evaluation of the controller performance for
5
(
)
(
s tep change 1=10 kJ /h in the feed en-
a s tep change in the feed flow rate + 150
)
thalpy.
kmol h .
/
design an inferential controller. The focus was given to the
requirements for the data, specifically what information the
data are required to possess and which data types are the
most appropriate. Our conclusion is that, for many systems,
both historical operation data and plant test data are needed
to build an adequate model for inferential control. This is in
contrast to the conventional feedback controller design for
which plant test data alone are usually adequate. We pro-
posed a two-step procedure for building a model using both
Figure 11. Evaluation of the controller performance for
s tochas tic changes in the feed compos ition,
feed temperature, and feed flow rate.
Table 7. Evaluation of the Controller Performance for a Step
5
(
)
Change 1=10 kJrh in the Feed Enthalpy
Mean Square Error
y1
y2
Open loop
1.8717ey006
9.3407ey007
1.5466ey008
1.5481ey008
3.9123ey008
1.8915ey008
5.5139ey010
5.5108ey010
With delayed comp. meas.
With temp. meas. only
With del. comp. meas.
and temp. meas.
Figure 9. Evaluation of the controller performance for a
(
s tep change in the compos ition external
)
propylene s tream flow rate + 10 kmol/h .
1986
October Vol. 46, No. 10
AIChE Journal

Andersen, H. W., M. Kummel, and S. B. Jorgensen, ‘‘Dynamics and
Identification of a Binary Distillation Column,’’ Chem. Eng. Sci.,
Table 8. Evaluation of the Controller Performance for a Step
Change in the Composition
Ž
.
44, 2571 1989 .
Mean Square Error
˚
Astrom, K. J., and B. Wittenmark, Computer-Controlled Systems: The-
¨
y1
y2
Ž
.
ory and Design, Prentice Hall, Upper Saddle River, NJ 1997 .
Banerjee, A., Y. Arkun, B. Ogunnaike, and R. Pearson, ‘‘Estimation
of Nonlinear Systems Using Linear Multiple Models,’’ AIChE J.,
Open loop
3.7981ey006
1.8714ey006
1.3615ey007
4.7393ey008
8.5512ey008
4.4615ey008
5.1916ey009
1.9616ey009
with delayed comp. meas.
With temp. meas. only
With del. comp. meas.
and temp. meas.
Ž
.
43, 1204 1997 .
Cambell, H. C., and J. B. Rawlings, ‘‘Predictive Control of Sheet and
Film Forming Processes,’’ AIChE J., 44, 1162 1998 .
Cooley, B., and J. H. Lee, ‘‘Control-Relevant Experiment Design,’’
AIChE J. 2000 .
Cutler, C. R., and F. H. Yocum, ‘‘Experience with the DMC Inverse
for Identification,’’ CPC-IV, AIChE Symp. Ser., P-67, 297 1991 .
Doyle, F. J., III, ‘‘Nonlinear Inferential Control for Process Applica-
tions,’’ J. Process Control, 8, 339 1998 .
Forssell, U., and L. Ljung, ‘‘Closed-Loop Identification Revisited,’’
Automatica, 35, 1215 1999 .
Friedland, B., ‘‘Sampled-Data Control Systems Containing Periodi-
cally Varying Members,’’ Proc. IFAC Cong., 1960 .
Ž
.
Ž
.
Note: External propylene stream flow rate q10 kmolrh.
Ž
.
Table 9. Evaluation of the Controller Performance for a Step
Change in the Feed Flow Rate
Ž
.
Ž
.
Mean Square Error
y1
y2
Ž
.
Jacobson, E. W., and S. Skogestad, ‘‘Inconsistencies in Dynamic
Models for Ill-Conditioned Plants: Application to Low-Order
Models of Distillation Columns,’’ Ind. Eng. Chem. Res., 33, 631
Open loop
4.4710ey006
2.8145ey006
6.1400ey008
4.1067ey008
6.5296ey007
3.3486ey007
2.3039ey008
2.0032ey008
With delayed comp. meas.
With temp. meas. only
With del. comp. meas.
and temp. meas.
Ž
.
1994 .
Kalman, R. E., ‘‘A New Approach to Linear Filtering and Prediction
Ž
.
Problems,’’ Trans. ASME, Ser. D: J. Basic Eng., 82, 35 1960 .
Kalman, R. E., and R. S. Bucy, ‘‘New Results in Linear Filtering and
Prediction Theory,’’ Trans. ASME, Ser. D: J. Basic Eng., 83, 95
Note: q150 kmolrh.
Ž
.
1961 .
Kranc, G. M., ‘‘Input-Output Analysis of Multi-Rate Feedback Sys-
Table 10. Evaluation of the Controller Performance
for Stochastic Changes in the Feed Composition, Feed
Temperature, and Feed Flow Rate
Ž
.
tems,’’ IRE Trans. Automat. Contr., AC-3, 21 1957 .
Larimore, W. E., ‘‘Canonical Variate Analysis in Identification, Fil-
tering, and Adaptive Control,’’ Proc. IEEE Conf. on Decision and
Ž
.
Control, Hawaii, 596 1990 .
Mean Square Error
Lee, J. H., and A. K. Datta, ‘‘Nonlinear Inferential Control of Pulp
y1
y2
Ž
.
Digester,’’ AIChE J., 40, 50 1994 .
Lee, J. H., M. Morari, and C. E. Garcıa, ‘‘State Space Interpretation
´
Open loop
8.5288ey007
1.1139ey008
7.0732ey009
2.5188ey008
1.2428ey009
8.0135ey010
Ž
.
of Model Predictive Controll,’’ Automatica, 30, 707 1994 .
Lee, J. H., M. S. Gelormino, and M. Morari, ‘‘Model Predictive Con-
trol of Multi-Rate Sampled-Data Systems: A State-Space Ap-
With delayed comp. meas.
With temp. meas. only
With del. comp. meas.
and temp. meas.
Ž
.
proach,’’ Int. J. Control, 55, 153 1992 .
9.4616ey010
1.3526ey010
Lee, J. H., and M. Morari, ‘‘Robust Inferential Control of Multi-rate
Ž
.
Sample-data Systems,’’ Chem. Eng. Sci., 47, 865 1992 .
Li, D., S. L. Shah, T. Chen, and R. Patwardhan, ‘‘System Identifica-
tion and Long-Range Predictive Control of Multi-Rate Systems,’’
types of data. We also brought out some practical issues that
may arise during an application of the procedure and pro-
posed solutions. Finally, we showed the efficacy of the method
in a case study involving a mulitcomponent distillation col-
umn simulated in HYSYS.
As industries continue to invest heavily in new sensors and
information-management systems, the volume of plant data
being collected is rising rapidly. While data may come cheap,
unlocking values hidden in them takes a great many re-
sources. The proposed method can serve as a tool for extract-
ing value from collecting data.
Ž
.
Proc. Amer. Control Conf., San Diego, California, p. 336 1999 .
Li, W., and J. H. Lee, ‘‘Control Relevant Identification of Ill-Condi-
tioned Systems: Estimation of Gain Directionality,’’ Comp. Chem.
Ž .
Ž
.
Eng., 20 8 , 1023 1995 .
Ljung, L., System Identification: Theory for the User., Prentice Hall,
Ž
.
Englewood Cliffs, NJ 1987 .
Ljung, L., and T. McKelvey, ‘‘Subspace Identification from Closed
Ž
.
Loop Data,’’ Signal Process., 52, 209 1996 .
Ljung, L., and U. Forssell, ‘‘An Alternative Motivation for the Indi-
rect Approach to Closed-Loop Identification,’’ IEEE Trans. Au-
Ž
.
tomat. Contr., AC-44, 2206 1999 .
Morari, M., and G. Stephanopoulos, ‘‘Minimizing Unobservability in
Ž
.
Inferential Control Schemes,’’ Int. J. Control, 31, 367 1980 .
Muske, K. R., and J. B. Rawlings, ‘‘Model Predictive Control with
Ž
.
Linear Models,’’ AIChE J., 39, 262 1993 .
Acknowledgment
We gratefully acknowledge the financial support from the National
Science Foundation’s Young Investigator Program through the Grant
CTS 9357827.
Quintero-Marmol, E., and W. L. Luyben, ‘‘Inferential Model-Based
Control of Multicomponent Batch Distillation,’’ Chem. Eng. Sci.,
Ž
.
47, 887 1992 .
Ricker, N. L., ‘‘Model Predictive Control: State of the Art,’’ Proc. Int.
Conf. on Chemical Process ControlᎏCPC IV, Y. Arkun and W.
Ž
.
Ray, eds., South Padre Island, TX, p. 271 1991 .
Literature Cited
Amirthalingam, R., and J. H. Lee, ‘‘Subspace Identification Based
Russell, S. A., D. G. Robertson, J. H. Lee, and B. A. Ogunnaike,
‘‘Quality Monitoring and Control of Nylon 6,6 Autoclave,’’ Chem.
Inferential Control of a Continuous Pulp Digester,’’Comput. Chem.
Ž
.
Eng. Sci., 53, 3685 1998 .
Ž
.
Eng., 21, S1143 1997 .
Tatiraju, S., M. Souroush, and R. Mutharasan, ‘‘Multi-Rate Nonlin-
ear State and Parameter Estimation in a Bioreactor,’’ Biotechnol.
Amirthalingam, R., and J. H. Lee, ‘‘Subspace Identification Based
Inferential Control Applied to a Continuous Pulp Digester,’’ J.
Ž .
Ž
.
Bioeng., 63 1 , 22 1999 .
Ž
.
Process Control, 9, 397 1999 .
Tatiraju, S., and M. Souroush, ‘‘Nonlinear State Estimation in a
AIChE Journal
October Vol. 46, No. 10
1987

Ž
.
Polymerization Reactor,’’ Ind. Eng. Chem. Res., 36, 2679 1997 .
Van Overschee, P., and B. De Moor, ‘‘N4SID: Subspace Algorithms
for the Identification of Combined Deterministic Stochastic Sys-
Yu, C. C., and W. L. Luyben, ‘‘Use of Multiple Temperatures for the
Control of Multicomponent Distillation Columns,’’ Ind. Eng. Chem.
Ž
.
Process Des. De®., 23, 590 1984 .
Ž
.
tems,’’ Automatica, 30, 75 1994 .
Yu, Z., W. Li, J. H. Lee, and M. Morari, ‘‘State Estimation Based
Model Predictive Control Applied to Shell Control Process: A Case
Study,’’ Chem. Eng. Sci., 18, 15 1994 .
Verhaegen, M., ‘‘A Novel Non-iterative MIMO State Space Identifi-
cation Technique,’’ Preprints 9th IFACrIFORS Symp. on Identifica-
Ž
.
Ž
.
tion and System Parameter Estimation, Budapest, Hungary 1991 .
Yiu, Y., G. A. Carling, and R. K. Wood, ‘‘Dynamic and Control Be-
haviour of Multicomponent Distillation Columns,’’ Simulation, 201
Ž
.
1990 .
Manuscript recei®ed July 8, 1999, and re®ision recei®ed Apr. 21, 2000.
1988
October Vol. 46, No. 10
AIChE Journal