University of Southern California Dissertations and Theses

Performance monitoring and disturbance adaptation for model predictive control

(USC Thesis Other)

Performance monitoring and disturbance adaptation for model predictive control

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content PERFORMANCE MONITORING AND DISTURBANCE
ADAPTATION FOR MODEL PREDICTIVE CONTROL
by
Zhijie Sun
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(CHEMICAL ENGINEERING)
December 2012
Copyright 2012 Zhijie Sun
ii
Dedication
To my family.
iii
Acknowledgments
I would like to thank everyone who encouraged me, supported me
and inspired me during my doctoral work. In fact, this work would not be
finished without these people influencing me in every aspects of my life.
Foremost, I would like to offer my sincerest gratitude to my advi-
sor, Prof. S. Joe Qin. Throughout my PhD research, he guided me to new
research areas, encouraged me to develop new ideas and to think innova-
tively, helped me improve my analytical and technical skills, and assisted
me with documenting the results. All of these are the essential elements in
conducting research. In addition, I learned much from him on professional-
ism and communications, which I believe are more important in developing
my career in the future. I would also thank Mork Family Graduate Fellow-
ship, Praxair Inc., and the Center for Interactive Smart Oilfield Technologies
(CiSoft) for their financial support during my graduate studies.
I am grateful for having an exceptional dissertation committee and
would like to thank other members of the committee, Prof. Michael Safonov,
and Prof. Pin Wang for their guidance and insights. I would also acknowl-
edge Prof. James Rawlings and Prof. Thomas Edgar for their valuable com-
ments and suggestions during the semiannual Texas-Wisconsin-California
iv
Control Consortium (TWCCC) meetings.
The experience with the entire Qin group was enjoyable and mem-
orable. I feel very fortunate to have the former members who offered me
useful advices and helpful suggestions, including Carlos, Quan, Yamin, Bo,
and visiting scholars Prof. Li and Prof. Yang. I would also like to thank
the current members of the group, Tao, Yu, Jingran, Hu, Yingying, Johnny,
Yining, Prof. Sun and Prof. Liu, for their friendliness and consistent help.
Especially I appreciate the time spending with Yu Zhao on our joint work.
He had a deep understanding about Kalman filter theory and taught me
everything he knew about system identification.
Deep gratitude also goes to Dr. Ashish Singhal, Dr. Larry Megan
and engineers of Praxair. They provided suggestions from the view of in-
dustry and offered me the platform to validate my research work. I was also
lucky to have many great instructors of the courses I took, particularly Prof.
Thomas Marlin and Dr. Stephen Stoyan. Advised by Prof. Marlin, I worked
on the class project that was related to my research, and finally the results
leads to a conference paper. Dr. Stoyan was a great instructor who offered
me the basic research tools in the field of optimization; I also acquired a lot
of practical teaching skills from him when I was a teaching assistant of him
for three semesters.
I have been so fortunate to have my lasting friendships over my jour-
ney to PhD. Sometimes our discussions were on academic-related topics
v
than on a personal basis. I would like to thank Shilin, Yue and Peng for
the inspiring discussions on mathematical fundamentals, statistical models
and optimizations. Their viewpoints from other disciplines refreshed my
understanding on my ongoing research work. Thanks also go to my friends
like Ming, Wei, Jianshu, Jiaying and old friends in China for their constant
support.
Finally, I am indebted to thank my family for their continuous en-
couragement, support and advice.
vi
Table of Contents
Dedication ii
Acknowledgments iii
List of Tables ix
List of Figures x
Abstract xii
Chapter 1. Introduction 1
1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Chapter 2. Adaptive EWMA-based Run-to-Run Control 14
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2 EWMA-Based Controllers . . . . . . . . . . . . . . . . . . . . . 16
2.2.1 Process description . . . . . . . . . . . . . . . . . . . . . . 16
2.2.2 EWMA controllers . . . . . . . . . . . . . . . . . . . . . . 16
2.3 Adaptive Run-to-Run Control Design . . . . . . . . . . . . . . . 18
2.3.1 Conventional adaptive control design . . . . . . . . . . . 18
2.3.2 Proposed adaptive run-to-run control scheme . . . . . . 20
2.4 Convergence of the Adaptive Run-to-Run control . . . . . . . . 23
2.5 Simulation Examples . . . . . . . . . . . . . . . . . . . . . . . . 25
2.5.1 Linear process with IMA(1,1) disturbance . . . . . . . . 26
2.5.2 Nonlinear chemical mechanical polishing (CMP) process 31
2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
Chapter 3. MPC Disturbance Modeling and Adaptation for an In-
dustrial Process 35
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
vii
3.2 Offset-free Model Predictive Control and Disturbance Model . 37
3.2.1 Internal model and disturbance . . . . . . . . . . . . . . 37
3.2.2 Dynamic matrix control . . . . . . . . . . . . . . . . . . . 38
3.3 Disturbance Model Estimation and Adaptation . . . . . . . . . 42
3.3.1 IMA(1,1) model . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.2 General disturbance model . . . . . . . . . . . . . . . . . 43
3.3.3 Disturbance Adaptation for MPC . . . . . . . . . . . . . 46
3.4 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.4.1 Simulated distillation column . . . . . . . . . . . . . . . 48
3.4.2 An industrial air separation process . . . . . . . . . . . . 54
3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
Chapter 4. Improving Industrial MPC Performance with Data-
Driven Disturbance Modeling 60
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2 Dynamic Matrix Control . . . . . . . . . . . . . . . . . . . . . . 62
4.2.1 Plant model in step response form . . . . . . . . . . . . . 62
4.2.2 DMC predictions and offset-free control . . . . . . . . . 64
4.2.3 Plant model in state-space form . . . . . . . . . . . . . . 65
4.3 Subspace Identification of a Disturbance Model . . . . . . . . . 66
4.3.1 Identification of the observer Markov parameters . . . . 66
4.3.2 System Markov parameters of the disturbance model . . 68
4.4 DMC with Moving Average Disturbance Model . . . . . . . . . 70
4.5 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
Chapter 5. Performance Monitoring of Model-Based Controllers via
Model Residual Assessment 79
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.2 Feedback Invariance of Disturbance Innovations . . . . . . . . 84
5.2.1 Motivating SISO processes . . . . . . . . . . . . . . . . . 84
5.2.2 General MIMO processes . . . . . . . . . . . . . . . . . . 85
5.3 Estimating Disturbance Innovations from Routine Closed-
Loop Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.4 Estimating Disturbance Innovations from Mismatched Model
Residuals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
viii
5.5 Data-Driven Model Quality Assessment in MPC . . . . . . . . 93
5.5.1 MPC objective and key performance index . . . . . . . . 93
5.5.2 Model residual assessment . . . . . . . . . . . . . . . . . 95
5.6 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.6.1 Wood-Berry distillation column . . . . . . . . . . . . . . 96
5.6.2 Tennessee Eastman challenge problem . . . . . . . . . . 100
5.6.3 An air separation plant of Praxair . . . . . . . . . . . . . 104
5.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
Chapter 6. Control Performance Monitoring of LP-MPC Cascade
Systems 108
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
6.2 Revisit of MIMO Interactor and Minimum Variance Benchmark110
6.3 Introduction of Block Lower Triangular Interactor Matrix . . . 113
6.4 LP-based MPC Control Performance Monitoring . . . . . . . . 118
6.4.1 LP problem and active CV constraints . . . . . . . . . . . 118
6.4.2 Conditional MVC law and performance benchmark
under block lower triangular interactor . . . . . . . . . . 120
6.5 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . 123
6.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
Chapter 7. Conclusions 128
Bibliography 131
ix
List of Tables
Table 2.1 Control performance of the controllers regulating CMP pro-
cess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Table 3.1 Output Variance of Different MPCs . . . . . . . . . . . . . . . 50
Table 3.2 Effects of Disturbance Model Adaptation . . . . . . . . . . . 52
Table 3.3 Model Quality Index of Results of the Industrial Air Sepa-
ration Process . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Table 4.1 Control Performance of Different DMC Schemes . . . . . . . 77
Table 5.1 Model Residual Assessment with Perfect Plant Model but
Mismatched Disturbance Models . . . . . . . . . . . . . . . . 99
Table 5.2 Model Residual Assessment with Mismatched Plant and
Disturbance Models . . . . . . . . . . . . . . . . . . . . . . . . 100
Table 5.3 MQI of the TEP Example . . . . . . . . . . . . . . . . . . . . . 103
Table 5.4 Model Quality Index of Results of the Industrial Air Sepa-
ration Process . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
x
List of Figures
Figure 2.1
k
by PEM with b = 2 and past data window size of 50,
100 and 200. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
Figure 2.2
k
by PLR withb = 2 and past data window size of 50, 100
and 200. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
Figure 2.3
k
by PEM withb = 3 and window size of 50, 100 and 200. 28
Figure 2.4
k
by PLR withb = 3 and window size of 50, 100 and 200. . 28
Figure 2.5 Time-varying disturbance parameter in the simulation. . 29
Figure 2.6
k
with window size of 100 under time-varying disturbance. 30
Figure 2.7
k
with window size of 200 under time-varying disturbance. 30
Figure 2.8
k
for CMP process with window size of 200. . . . . . . . . 32
Figure 3.1 Schematic diagram of MPC disturbance adaptation. . . . . 47
Figure 3.2 Simplified diagram of an industrial air separation process . 54
Figure 3.3 Control results of the air separation process by MATLAB
simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
Figure 4.1 Identified system Markov parameters of the disturbance
model in case I. . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Figure 4.2 Identified system Markov parameters of the disturbance
model in case II. . . . . . . . . . . . . . . . . . . . . . . . . . 74
xi
Figure 4.3 Identified system Markov parameters of the disturbance
model in case III a). . . . . . . . . . . . . . . . . . . . . . . . 76
Figure 4.4 Identified system Markov parameters of the plant and dis-
turbance model in case III b). . . . . . . . . . . . . . . . . . . 77
Figure 5.1 Schematic diagram of a closed-loop controlled process. . . 81
Figure 5.2 Relationship of measurements, predictions, simulation er-
rors and their orthogonal projections. . . . . . . . . . . . . . 93
Figure 5.3 Autocorrelation of the first element ofe

(k), lag 1. . . . . 98
Figure 5.4 Schematic diagram of the Tennessee Eastman process. . . . 100
Figure 5.5 Step responses of plant model of MPC1 controlling TEP . . . 102
Figure 5.6 Step responses of plant model of MPC2 controlling TEP . . . 102
Figure 5.7 Control results of TEP by MPC1 and MPC2. . . . . . . . . . 104
Figure 5.8 Simplified diagram of an industrial air separation process . 105
Figure 6.1 Control results of Example 6.2 iny
1
-y
2
plane. . . . . . . . . 126
xii
Abstract
Model predictive control (MPC) is a widely used advanced process
control technique in the process industry. According to the internal model
principle, the internal model of MPC has to include both exact plant and
disturbance models to be optimal. However, in practice, the MPC usually
assumes a step-like disturbance or a fixed disturbance model. As a result,
the MPC will be suboptimal when disturbance changes slowly. Moreover, it
lacks a tool assessing the optimality of control performance in terms of the
MPC model.
In this dissertation, a new MPC disturbance adaptation method is
presented. Starting from a single-input-single-output (SISO) semiconductor
manufacturing process, we replaced the conventional run-to-run controller
by an adaptive EWMA controller. It is shown that the plant model mis-
match can be compensated by adapting the disturbance model. Analysis
has been done to show that the adaptive controller is stable and converges
to the optimal controller.
The proposed method is then extended to multi-input-multi-output
(MIMO) systems. For the ease of practical applications, the integrated mov-
ing average (IMA) model with order (1,1) is recommended. The equiva-
xiii
lence between the IMA(1,1) parameter and the prediction error filter con-
stant in commercial MPC has been established. Implementation of distur-
bance adaptation is explained.
Another disturbance modeling tool is presented. It focuses on the
closed-loop identification of a nonparametric disturbance model. The
method incorporates the plant model information during the conversion
from observer Markov parameters to system Markov parameters.
A new control performance assessment method evaluating MPC
model quality is then presented. Feedback invariant principle is intro-
duced, based on which a method estimating disturbance innovations is
given. A model quality index is developed as the performance benchmark,
which compares prediction errors with disturbance innovations. It is
shown that the model quality index related to the MPC performance index.
Most industrial processes are optimized by a linear programming
(LP) problem on top of the MPC. A new control performance monitoring
method for cascaded LP-MPC system is developed. The block-lower-
triangular interactor matrix is introduced to form a new method that is
able to determine the performance benchmark based on controlled variable
(CV) priorities coming from the LP results.
1
Chapter 1
Introduction
1.1 Overview
Being one of the most successful advanced process control schemes,
model predictive control (MPC) is capable of handling multivariable con-
strained control problems. At each control interval, MPC calculates the op-
timal future control signals based on the predictions of future outputs by
the process model. The first input of control sequence is fed into the plant.
At next control interval, the optimization is performed again to get the new
control signals. In terms of industrial applications, it has been reported by
Qin and Badgwell [68] that there had been more than 4,500 industrial linear
MPCs in different areas from refining to aerospace a decade ago. Recent sur-
vey on process control practice in Japan reveals that linear MPC is the most
popular among the advanced process control methods that are standard-
ized or applied [39]. Other research in emerging technologies also adopt
MPC as their control scheme. For example, in fuel cell systems, MPC has
been successfully applied to avoid the lack of oxygen [82].
An MPC model is one of the key factors affecting the control perfor-
mance. The average life of an MPC is between 3 and 5 years in the process
2
industry because of the slow degradation of the performance [36]. This is
partly due to the slowly time-varying factors of the process or process re-
vamp. The MPC model includes two parts: the process and disturbance
model. The process model describes the relationship between manipulated
variables (MV) or inputs and controlled variables (CV) or outputs, and the
disturbance model is the shaping filter of the unmeasured disturbances. The
slow time-varying changes of the process may arise in both the process and
disturbance models. Under such circumstances, there are two categories of
MPC schemes to achieve a good control performance.
One type is the robust MPC [8, 91, 72]. It takes every model in the
uncertainty set into account when calculating the optimal control sequence.
Bemporad and Mosca [6] uses reference trajectory in predictive control. It
has been shown by [4] that offset-free tracking can be achieved when distur-
bance converges to zero. Wang and Rawlings [85] proposed an offset-free
robust MPC based on predictions of model uncertainties. One of the lim-
itation of robust MPC is that it requires disturbances to be bounded such
that uncertainty set is bounded. An interesting adaptive strategy proposed
by DeHaan et al. [12] adapts parametric uncertainty sets in constrained
nonlinear robust MPC. This method gives more flexibility in uncertainty
boundaries. In real plants, however, some disturbances are nonstationary
and hence unbounded. Therefore, the conventional robust MPC is unable
3
to account all possible disturbances into uncertainty set.
The other type of MPC augments plant model by disturbance model,
and output feedback is used. The differences between measured outputs
and predicted outputs are considered to be the disturbance, which may con-
sists of unmeasured disturbances and mismatch in process and disturbance
model. Since in practice a priori information about disturbances as well as
true process model is often unavailable, MPC such as dynamic matrix con-
trol (DMC) [11, 25] usually adopts a step-like disturbance model, which is
considered to be the best model [22, 23]. That is, after obtaining the mea-
sured CV , MPC treats the disturbance to be the same over the prediction
horizon. In the calculation, all outputs predicted by MPC process model
will be offset by the amount of disturbance. Therefore, disturbances serve
as feedback for controller, and MPC is corrected by measurements at each
control interval.
MPC could take advantage of the augmentation of the step distur-
bance model. The integral action of MPC can be generated by including
output step disturbance for stable systems [57] or state step disturbance
model for unstable systems in the model [56]. Pannocchia and Rawlings
[63] proved that for square systems without integrating modes offset-free
tracking can be achieved when using a step disturbance model. They also
argue that for non-square systems the number of output step disturbance
4
models must be the same as number of measured variables to guarantee
removing of offset when no constraint is active. Recent work by Maeder et
al. [51] provides design guidelines for offset-free tracking using step dis-
turbance when the number of tracked variables is less than the number of
measured variables.
However, the step disturbance model does not always give MPC best
predictions. According to internal model principle [19], to track reference
signals without offset disturbance dynamics need to be included in the con-
troller. Moreover, since the control performance of MPC depends on the
internal model quality, it calls for a general disturbance model instead of
step-like disturbance. Morari and Stephanopoulos [54] shows that state fil-
ter with a nonstationary disturbance model could result in superior control
performance with the existence of unmeasured variables. It is suggested
by Garcia and Morari [21] that a diagonal structure of first order filter be
used to filter updates. They prove that with appropriate filter parameter
chosen the controller would be capable to stabilize systems with virtually
any plant/model mismatch. Later in the survey by Prett et al. [25] such
an update filter is converted to the form of a first-order disturbance model
which is then related to the robustness of the controller. More complex dis-
turbance modeling is further suggested by Lee et al. [47]. By modeling
the disturbance using a state-space form, an optimal state estimator can be
5
obtained and gives best predictions. In addition, computational disadvan-
tages of step disturbance model has been addressed in the paper. However,
this work mostly focuses on the first order filter since it emphasizes the con-
nection between the disturbance model and the traditional MPC techniques
such as DMC. In Muske and Rawlings’ [57] formulation, the noise model
for measurement noise is added to output disturbance model. Muske and
Badgwell [55] studied a general form of disturbance model in offset-free
control. In their work, input, output and state disturbances are all included
in disturbance model. They derived requirements for disturbance models,
say, all states of the augmented system should be detectable. The same cri-
terion is developed independently by Pannocchia and Rawlings [63]. In
addition to the linear model, hidden Markov model is introduced recently
to model disturbances with time-varying characteristics such as intermit-
tent drifts and abrupt changes [86]. Therefore, the disturbance model needs
to be further refined for systems with complicated disturbances.
Choice of disturbance structure is an important issue in developing
disturbance models for MPC. An arbitrarily chosen disturbance structure
may contribute little improvement to the MPC performance. One exam-
ple is that improper selection of output step disturbance, the most com-
mon configuration in industrial MPCs, may even make controller not be
able to stabilize systems with unstable poles [56]. Shinskey [77] has pointed
6
out that MPC can outperform PID controllers on setpoint changes but not
on disturbances entering input of a process upstream with dominant lag.
Input disturbance models are introduced as one solution to this problem
[57, 55]. Instead of outputs, inputs are assumed to have additive distur-
bances. Their models are augmented to system in a similar way as output
disturbance. Some comparisons of different disturbance models are shown
in [63]. It is further suggested by Tenny and Rawlings [80] that an input
disturbance model performs better on systems with perturbed states. The
mechanism of effects of input and output disturbance model on model mis-
match is explained graphically. Which model to use can be determined in-
tuitively from the graph if a priori knowledge of process is available. Recent
work by Pannocchia [61] studies disturbance modeling for ill-conditioned
systems. Contrast to conclusions in [22], robust performance is achieved
by implementation of input disturbance model rather than output distur-
bance model, if system exhibits strong interactions. Rajamani et al. [70]
derived that input and output disturbance model are equivalent in terms of
closed-loop performance and differ by a similarity transform if covariances
are estimated appropriately. This result can be used to simplify disturbance
models. Hence, proper identification of the disturbance model structure
needs to be further investigated.
Adaptive strategies receive less attention as compared with robust
7
strategies in MPC area for the handling of model uncertainty. In the pio-
neer work by
˚
Astr¨ om and Wittenmark [3], self tuning regulator has been ap-
plied to minimum variance control. However, little progress has been made
specifically for MPC. It had been reported by Mayne et al. [52] that by 2000
there is no stabilizing adaptive MPC developed for constrained systems.
Recently, some MPC with adaptation is published. Karra et al. [40] pro-
posed a two-tier modeling scheme to develop adaptive MPC. They use out-
put error model to characterize plant and ARMA model to characterize un-
measured disturbance. Dougherty and Cooper [16] developed an adaptive
scheme based on multiple linear models. In their approach, future predic-
tions are determined by adaptive interpolation of different model outputs.
It is shown such strategy is effective for nonlinear systems. Fukushima et
al. [20] carry out an adaptive MPC algorithm for constrained system, where
process model is adapted in a moving horizon fashion. Unfortunately, their
implementation of MPC is extended from robust methods to handle uncer-
tainties.
Traditional model-based adaptive control methods try to model all
aspects of a system including process and disturbance model. However,
as aforementioned, industrial processes are usually changing slowly. Be-
sides, they are also so complicated that a large number of parameters are
needed to build the process models. In addition, Safonov and Tsao [75]
8
state that the data may conflict with the adaptive control law when existing
performance constraints. Therefore, methods adapting both process and
disturbance models could lead to unnecessary flexibility and problematic
controllers, since the disturbance usually varies faster than the process. Fur-
thermore, changes in the process can be counted towards disturbances, as
revealed in [22] that a first-order disturbance model improves robustness.
These issues calls for the research in adapting disturbance models to filter
system states optimally.
One example of MPC augmenting disturbances is run-to-run con-
trol in semiconductor manufacturing. The manufacturing process is repre-
sented by a linear single-input-single-out (SISO) model in which the input
is the recipe for some batch and the output is the quality variable. Run-to-
run control provides the recipe update from run to run since in situ mea-
surements are not available. The exponentially-weighted-moving-average
(EWMA) filter is augmented as disturbance model [74, 9]. According to Box
and Jenkins [7], when disturbance behaves as integrated moving average of
first order (IMA(1; 1)), the EWMA controller provides minimum variance
control, with EWMA weight equivalent to IMA disturbance parameter. The
same conclusion can be drawn from internal model principle [19]. Practi-
cally, the EWMA controller gives quite satisfactory results even when dis-
turbances are not necessarily IMA models. Due to its simplicity, run-to-run
9
control can be a starting point for studying MPC disturbance adaptation
techniques.
In practice, control performance of MPC is continuously monitored.
Control performance monitoring approaches compare some statistic from
input/output data with a benchmark, and then evaluate controller perfor-
mance by an index. A common type of monitoring methods is based on
minimum variance benchmark [27] whose multivariate version is devel-
oped in [28]. However, the MPC objective may have input reference penalty
and input move suppression term. Minimum variance metrics may not be
achieved. Huang and Shah [34] compare MPC with linear quadratic reg-
ulator (LQR) and propose the LQ benchmark. The disadvantage of this
method is similar to minimum variance benchmark, namely, LQ benchmark
cannot be attained when there are strict constraints. Some model based
approach such as design-case method [65] utilize model information and
evaluate MPC objective function value. Ko and Edgar [42] calculate bench-
mark using minimum variance prediction over future horizon. Recently,
Harrison and Qin [30] propose minimum variance map considering effects
of constraints. However, for MPC with disturbance model, the key factor
is quality of disturbance model since control performance is directly related
to model accuracy. The method discriminating Kalman gain and process
model mismatch is developed by Harrison and Qin [29], but it lacks the
10
performance evaluation in terms of MPC model. Furthermore, most of the
existing performance monitoring methods are lumped indices consisting
of all factors impacting the control performance, e.g., constraint setup and
MPC tuning. However, few methods have been found to monitor the model
quality for the purpose of MPC revamp or disturbance model adaptation.
Ji et al. [36] make a valuable attempted on this issue. Their approach is ca-
pable of quantifying the model errors through frequency analysis. Despite
this, the method is intrusive that special test signals are required. It is also
an empirical method which needs the choice of excitation frequencies by
the designer.
The steady-state values of MPC setpoints are usually optimized for
better economics. Linear programming (LP) is a typical steady-state opti-
mizer. MPC setpoints are determined by LP solutions. This forms a cascade
structure of LP-MPC. Ying and Joseph [89] investigated stability and im-
plementation issue of LP-MPC. It is stated that for a constrained two-stage
MPC, end condition at control horizon guarantees global stability. Rao and
Rawlings [71] provide a strategy handling active constraints when calcu-
lating steady-state target. Nikandrov and Swartz [58] perform sensitivity
analysis on LP stage constraints whose results help one explain the interac-
tions between LP and MPC. It is reported [43, 68] that an important under-
lying problem of LP-MPC scheme is that steady-state optimum may bounce
11
around. The reason is that LP solutions are always at vertex of feasible re-
gion. These unnecessary fluctuations may lead to poor control performance
as a result of the interactions between two layers. In addition, the LP so-
lutions indicate the priority of CVs by checking the active set. Therefore,
it is desirable to evaluate control performance of the LP-MPC system as a
whole.
1.2 Outline
To address and attack the problems reviewed in that last section, the
dissertation is organized as follows.
In Chapter 2, we start with the SISO EWMA-based run-to-run con-
trol with disturbance adaptation. Two estimation methods for disturbance
model identification are studied. A convergence analysis is performed to
show that the adaptive controller is stable and converges to the optimal
control law even with an incorrect process gain. The proposed adaptive
control scheme is tested on two simulation examples including an indus-
trial chemical-mechanical polishing process.
Chapter 3 extends the method to MIMO MPC. The plant and distur-
bance model as MPC internal model is introduced in the state-space form.
A practical IMA(1,1) disturbance model is presented, which is later shown
to be equivalent to MPC prediction error filtering. The steps of implemen-
tation of disturbance adaption is provided. One simulation example and
12
one industrial case are shown to illustrate the advantages of adapting dis-
turbance.
Chapter 4 presents a non-parametric closed-loop identification
method. By HOARX modeling of closed-loop data, the observer Markov
parameters can be obtained first. Then, a new approach is proposed to
convert the observer Markov parameters to the system Markov parameters
by utilizing the MPC model. It is shown that due to the integrator in MPC
the conversion can be performed using the process model only, instead
of the full augmented model. Integration of the proposed method with
DMC is demonstrated. Simulation of a distillation column is performed to
compare the proposed method with conventional DMC.
Chapter 5 gives a control performance monitoring approach in terms
of model quality. The feedback invariant principle is discussed first. The al-
gorithm for calculation of disturbance innovations is provided, which can
be based on either routine closed-loop data or mismatched model residuals.
The model quality assessment can be performed by comparison of the dis-
turbance innovations and the model residuals. Results from simulations of
a distillation column and Tennessee Eastman process as well as data of an
air separation process of Praxair are presented.
Control performance monitoring of LP-MPC systems is developed
in Chapter 6. CVs are prioritized based on LP solutions. The block lower-
triangular interactor matrix is introduced to describe the delay structure of
13
a MIMO process with CV priorities. The minimum variance for each CV
group with the same priority can be calculated from closed-loop data. The
control performance is evaluated by taking the ratio of CV variance and
minimum variance within each group. Simulation examples are presented
for comparison.
Chapter 7 gives some concluding remarks for the dissertation.
14
Chapter 2
Adaptive EWMA-based Run-to-Run Control
2.1 Introduction
Recent research in semiconductor manufacturing process has shown
a great interest in run-to-run control. It provides the recipe update from
run to run when in situ measurements are not available. The exponentially-
weighted-moving-average (EWMA) statistic is introduced to the run-to-run
control scheme for estimating process disturbances [74, 13]. The pioneering
work of Box and Jenkins [7] shows that when disturbance behaves as inte-
grated moving average (IMA) with order (1,1), the EWMA controller pro-
vides minimum variance control with the EWMA filter parameter relating
to the IMA disturbance parameter. The stability region of EWMA controller
was studied in [74].
However, most of existing EWMA controllers assume a time-
invariant process disturbance and a fixed weight factor. This is generally
not the case. For instance, after many runs, some parameters such process
gain and disturbance dynamics may change significantly as compared to
their values at the initial stage. The model parameters may vary slowly
with surrounding conditions, leading to time-varying disturbance models.
15
Recent research by [81] indicates the weakness of fixed weight factor,
which only focuses on the initial stage. Moreover, a fixed weight factor
may destabilize the process, as is pointed out in [14] that under-estimation
of process gain may results in a unstable system. Therefore, an adaptive
EWMA controller is needed to make the system more robust and, further-
more, have a better performance. Some examples of self-tuning adaptive
EWMA controllers were studied in [13]. Using a recursive least squares
(RLS) method, online estimation of the process gain is implemented to
achieve adaptive tuning. Unfortunately, such self-tuning controllers are
focused on estimation of the process gain. Patel and Jenkins [64] presented
an adaptive scheme for EWMA-based control. Their method successfully
tracks the optimal controller gain, but it is based on heuristics and the
convergence speed can be slow. In this chapter, we present an adaptive
EWMA controller whose weight factor is tuned adaptively according to
disturbance changes. Adapting the disturbance enables the controller to
vary the optimal weight factor accordingly to compensate for the incorrect
estimate of the process gain, and to reduce the performance degradation by
complex disturbance dynamics.
16
2.2 EWMA-Based Controllers
2.2.1 Process description
We assume that the process under run-to-run control is represented
by the following linear model.
y
k
=u
k1
+
k
+!
k
(2.1)
wherey
k
andu
k
are the controlled and manipulated variables, respectively,
and is the process gain. Specifically, in run-to-run control,y
k
is the mea-
sured response describing the process property at the end of kth run.
k
here is the deterministic process bias, and!
k
is assumed to be white noise.
2.2.2 EWMA controllers
In the EWMA controller based on (2.1), the deterministic bias
k
and
!
k
are considered together by defining

k
=
k
+!
k
:
The parameter are estimated from run-to-run by the EWMA expression:
a
k
, ^
k
=(y
k
bu
k1
) + (1)a
k1
(2.2)
where 0 < < 1 is the EWMA discount or weight factor, and b is the
estimated process gain used by the controller. Note that the estimateb may
not be close to the actual process gain due to noise or nonlinearity. A
larger value of implies the current bias estimatea
k
relies more on the past
estimates.
17
The feedback control law based on (2.2), known as EWMA control,
is
u
k
=
Ta
k
b
(2.3)
where T is the targeted response value. In practice, the weight factor is
tuned manually. The objective of the EWMA-based controller is to minimize
the mean-squared output tracking error
lim
k!1
E

e
2
k

(2.4)
by defininge
k
=y
k
T .
Rewriting (2.2), we obtain
a
k
=a
k1
+(y
k
bu
k1
a
k1
) (2.5)
Here, e
k
= y
k
bu
k
a
k1
is the one-step-ahead prediction error of the
EWMA controller. Therefore,
k
can be expressed as

k
=a
k
+e
k
=

1q
1
e
k
+e
k
=
1 (1)q
1
1q
1
e
k
(2.6)
Equation (2.6) implies that the EWMA-based run-to-run controller assumes
a disturbance described by the integrating moving average (IMA) model
with order (1,1) if the prediction error is white noise.
18
2.3 Adaptive Run-to-Run Control Design
2.3.1 Conventional adaptive control design
Consider the process (2.1) with IMA(1,1) disturbance (2.6). Without
loss of generality, we assume the setpoint isT = 0. One may first differenti-
ate both sides to eliminate the integrator in the disturbance:
y
k
=u
k1
+

1 (1)q
1

e
k
(2.7)
where is the real model parameter of the disturbance and y
k
=y
k
y
k1
.
The estimate of, denoted as
^
, can be either obtained by prediction error
method (PEM) assuming a moving average (MA) model with order (1,1), or
the pseudo linear regression (PLR).
Least-squares cannot be applied to (2.7) directly to estimate . The
pseudo linear regression, however, converts (2.7) to a linear model by letting
e
k1
=y
k1
:
y
k
=u
k1
(1)y
k1
+e
k
y
k1
=u
k1
(1)y
k2
+e
k1
.
.
.
y
kp+1
=u
kp
(1)y
kp
+e
kp+1
As a consequence,
^
b
k
and
^

k
, the estimates of and at time k, can be
obtained by

^
b
k
1
^

k

=

Z
T
p;k
Z
p;k

1
Z
T
p;k
y
p;k
(2.8)
19
whereZ
p;k
is defined as
Z
p;k
=
2
6
6
6
4
u
k1
y
k1
u
k2
y
k2
.
.
.
.
.
.
u
kp
y
kp
3
7
7
7
5
andy
p;k
= [y
k1
y
k2
y
kp
]
T
The identification method gives the estimates of both process gain
and disturbance parameter . Note that the estimates may be biased due
to lack of excitation. In the conventional adaptive control, we substitute the
estimates into the EWMA control law (2.2) and (2.3), i.e., lettingb =
^
b
k
and
=
^

k
:
a
k
=a
k1
+
^

k
(y
k

^
b
k
u
k
a
k1
)
u
k+1
=
a
k
^
b
k
However, estimating both and is redundant, since the EWMA-
based control is actually an integrating controller [64]:
u
k
=u
k1

b
y
k
(2.9)
It can be concluded from (2.9) that the controller is optimal as long as the
following relationship is satisfied.

b

=

(2.10)
where

and b

are the optimal parameters for the EWMA controller. In
addition, sinceu
k
can be expressed by past inputsu
i
(i = 1;:::;k 1) and
20
outputs y
i
(i = 1;:::;k 1), there may exist the problem of insufficient
excitation when implementing closed-loop identification. Therefore, it is
unnecessary to estimate both parameters. The controller (2.9) is optimal as
long as (2.10) holds. In consequence, we can setb

to some pre-determined
b and estimate

such that (2.10) holds.
2.3.2 Proposed adaptive run-to-run control scheme
In the proposed adaptive run-to-run control method, we choose to
fix the process gain estimate and to adapt the disturbance parameter. It
is equivalent to fixing the disturbance parameter and adapting the process
gain. Based on this idea, the adaptive EWMA-based run-to-run control will
be presented.
2.3.2.1 Determining the model gain
First, it is required to choose a combination of the process gain es-
timate b
0
and the initial guess
0
of the disturbance. From Kalman filter
theory, the IMA(1,1) process is convertible if and only if1< 1< 1. It
can be inferred that EWMA control is stable if and only if 0<< 2 . Note
that from (2.10) we have

=b=. Therefore, the fixedb should satisfy
0<
b

< 2: (2.11)
Without loss of generality, we assume > 0. Suppose the upper
and lower bounds of process gain are known as

and, respectively, i.e.,
21
0 < < <

. Since there is no more knowledge about and , it is
recommended to pick b = . Otherwise the optimal EWMA weight

is
more likely to destabilize the process.
2.3.2.2 Disturbance model adaptation
As aforementioned, adapting disturbance model only is sufficient to
guarantee the optimality of the adaptive EWMA control. To achieve distur-
bance adaptation, define the process disturbance as
d
k
=y
k
bu
k1
(2.12)
As the IMA(1,1) model is shown to be the assumption of EWMA based
control, we model the disturbance as
d
k
=
1 (1)q
1
1q
1
e
k
(2.13)
which is equivalent to
d
k
=ce
k1
+e
k
(2.14)
wherec = 1 is defined for simplicity.
Equation (2.14) suggests that d
k
is a moving average process. There
are two ways to obtain the estimate ofc at timek, denoted by ^ c
k
, similar to
the conventional adaptive control approaches:
Use PEM [48] to calculate ^ c
k
. The method involves nonlinear opti-
mization that minimizes the prediction errors. It converges quickly,
but more computation time for each iteration is expected.
22
Use the PLR. The prediction errore(k) equalsy(k) due to the deadbeat
characteristic of run-to-run control. Hence,
d
k
=cy
k1
+e
k
d
k1
=cy
k2
+e
k1
.
.
.
d
kp+1
=cy
kp
+e
kp+1
where p is the window size of past data. The ^ c
k
can be obtained by
least-squares:
^ c
k
=

y
T
p;k
y
p;k

1
y
T
p;k
d
p;k
(2.15)
wherey
p;k
= [y
k1
y
k2
y
kp
]
T
and d
p
(kp + 1) =
[d
k
d
k1
d
kp+1
]
T
.
Remark 2.1. Although the disturbance described by (2.13) is assumed to be
driven by the white noise, the prediction error e
k
may not be white when
suboptimal control is applied. However, e
k
is zero mean due to the inte-
grating action of the EWMA control. Therefore, the least-squares by (2.15)
is still unbiased.
With the estimate ^ c
k
, the EWMA weight becomes
k
= 1 ^ c
k
, and
the EWMA control law at timek can be expressed as
a
k
= (1
k
)a
k1
+
k
(y
k
bu
k
) (2.16a)
u
k+1
=
a
k
b
(2.16b)
23
2.4 Convergence of the Adaptive Run-to-Run control
Equation (2.7) suggests that the differentiated process follows a
MA(1) model. The convergence analysis of adaptive run-to-run control
law can then be performed based on the stochastic approximation theory
[44]. As aforementioned, conventional adaptive control updates both
process gain and disturbance parameter, while the proposed method only
adapts the controller gain. Therefore, the following lemma focusing on one
parameter is presented to show the optimal control parameter given by
(2.10) is an equilibrium point.
Lemma 2.1. Consider the system (2.1) with the adaptive EWMA-based controller
(2.16). Assume the model gain b has the same sign as the process gain . Then,

=
b

is an equilibrium point of the EWMA weight sequence
k
.
Proof. Define the unknown parameter as
=

b

Let the estimate of at timek be
^

k
=

b

k
Since the estimator is designed to minimize the variance of one-step-ahead
prediction errors
min
1
2
Ef(y
k
^ y
k
)
2
g
24
Following the derivations in [49], the recursions can be expressed by
^

k
=
^

k1
+

k
k
1
r
t1

k
"
k
(2.17a)
r
k
=r
k1
+

k
k

2
k
r
k1

(2.17b)
where
k
=k is the weight sequence and"
k
=y
k
^ y
k
. If estimation is achieved
by PEM,
k
is the gradient of ^ y
k
() with respect to; if estimation is achieved
by methods of least-squares type,
k
is simplyy
k
.
The ODE corresponding to (2.17) is
d()
d
=
()
r()
E

f
k
"
k
g (2.18a)
dr()
d
=()

E

f
2
k
gr()

(2.18b)
In order to obtain an equilibrium point for
^
,E

f (k)"(k)g is needed first.
DenoteC(q) = 1 (1)q
1
. Substituting the control law (2.9) into
the differentiated system (2.7) gives
y
k+1
=y
k
+y
k
(1)e
k
+e
k+1
= (1)y
k
+ ( +)y
k
+C(q)e
k+1
= [1C(q)]y
k+1
+ ( +)y
k
+C(q)e
k+1
Therefore,
C(q)y
k+1
= ( +)y
k
+C(q)e
k+1
which is equivalent to
y
k+1
=C
1
y
k
( +) +e
k+1
25
Since "
k
= y
k
^ y
k
and ^ y
k
= 0 under run-to-run control, the ODE (2.18a)
becomes
d()
d
=
()
r()
E

f (k)C
1
y
k
( +)g
showing that d()=d = 0 when + = 0. Therefore,

=
b

is an
equilibrium point of
k
.
To examine the stability of this equilibrium point, we need to check
the Lyapunov function
V () =
1
2
( +)
2
(2.19)
When PLR is used, the derivative of Lyapunov function (2.19) is always neg-
ative sinceC(q) is strictly positive real [44] as long as 0 << 2. When pa-
rameter is estimated by PEM method, the negative property ofdV (())=d
is shown in [49]. Combination of these two cases leads to the following
lemma.
Lemma 2.2. Consider the process (2.1) controlled by the adaptive run-to-run con-
troller (2.16). Assume the disturbance parameter satisfies 0 <
b

< 2. Then
the estimate of controller weight
k
converges to the optimal weight

, and the
adaptive control law converges to the minimum variance control.
2.5 Simulation Examples
In this section, two simulation examples will be presented to illus-
trate the proposed method. Both PEM and PLR are applied to identify the
26
controller gain. The first example shows how the run-to-run controller
adapts when process consists of only a linear gain. The second exam-
ple gives the case when process is nonlinear and studies how adaptive
controller handles nonlinearities.
2.5.1 Linear process with IMA(1,1) disturbance
Suppose the process observes the model (2.1). The process gain is
= 2 and the disturbance model parameter is = 0:5. A total of 1000 data
samples are collected.
First, the speed of convergence is studied, since the more data used
for identification result in less variance of EWMA weights while the less
data for identification imply a better tracking of true values of disturbance
parameters or optimal weights. The adaptation scheme is not applied until
the number of runs reaches the window size. The initial guess is
0
=

=
0:5.
Fig. 2.1 and 2.2 show the converging speed of PEM and PLR with
different past data windows, respectively, when the process gain is known
correctly, i.e.,b = 2. As can be seen, with the past window size of 50, either
method is unable to obtain a good estimate for most of the time. While the
window size is increased to 200, the estimates are quite close to the true
value with less variations.
However, when the process gain is known incorrectly asb = 3, the
27
0 200 400 600 800 1000
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

window=50
window=100
window=200
Figure 2.1:
k
by PEM withb = 2 and past data window size of 50, 100 and
200.
0 200 400 600 800 1000
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

window=50
window=100
window=200
Figure 2.2:
k
by PLR withb = 2 and past data window size of 50, 100 and
200.
28
0 200 400 600 800 1000
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2

window=50
window=100
window=200
Figure 2.3:
k
by PEM withb = 3 and window size of 50, 100 and 200.
0 200 400 600 800 1000
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2

window=50
window=100
window=200
Figure 2.4:
k
by PLR withb = 3 and window size of 50, 100 and 200.
29
0 200 400 600 800 1000
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figure 2.5: Time-varying disturbance parameter in the simulation.
optimal weight

=
b

= 0:75. Fig. 2.3 and 2.4 show the weights
k
based
on the disturbance adaptation method. It can be verified that weights
k
by
both adaptive schemes are able to converge to

without prior knowledge
of the true process gain.
Then, we will compare estimation with different window size under
slow time-varying disturbances. Suppose the disturbance parameter is
sinusoidal over time, as shown in Fig. 2.5. It varies so slow that within a
window of 200 the changes of its value do not exceed 0.3. The results of
EWMA weights by adaptive run-to-run control are illustrated by Fig. 2.6
and 2.7. It is demonstrated that the adaptive weights are capable to track
the slow time-varying parameters. Using a longer past data horizon usually
leads to a weaker tracking ability or a larger phase difference between esti-
30
0 200 400 600 800 1000
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

φ
λ
k
by PEM
λ
k
by PLR
Figure 2.6:
k
with window size of 100 under time-varying disturbance.
0 200 400 600 800 1000
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

φ
λ
k
by PEM
λ
k
by PLR
Figure 2.7:
k
with window size of 200 under time-varying disturbance.
31
mated and true values in this case, while a shorter past data horizon results
in a stronger tracking ability but larger errors in estimates.
The choice of parameter identification method depends on the dis-
turbance characteristics. The underlying assumption of PLR is the ARMAX
model (or MA(1,1) model in the adaptive run-to-run control); the PEM is
based on the stochastic gradient method [49] to find the estimate that min-
imizes the prediction error. Therefore, if the disturbance can be described
by the IMA(1,1) very well, then PLR can serve as a good candidate, as re-
vealed by the first case with time-invariant disturbances; otherwise, PEM is
recommended, for example, in the time-varying disturbance case.
Therefore, the proposed method has been validated that it is able to
adaptively adjust the controller gain to its optimum only by closed-loop
data and to adapt the controller to slow time-varying disturbances.
2.5.2 Nonlinear chemical mechanical polishing (CMP) process
A chemical mechanical polishing process modeled by SEMATECH
[15, 46] is simulated. The process has four manipulated variables and two
controlled variables. Among these, the manipulated variables are plate
speed (u
1
), back pressure (u
2
), polishing downforce (u
3
), and the profile of
the conditioning system (u
4
); the controlled variables are removing rate (y
1
)
and with-in wafer non-uniformity (y
2
). The target values are 2000 and 100
fory
1
andy
2
, respectively. A total of 1000 runs are simulated.
32
0 500 1000
−0.2
0
0.2
0.4
0.6
(a)
0 500 1000
−0.2
0
0.2
0.4
0.6
(c)
0 500 1000
−0.2
0
0.2
0.4
0.6
(b)
0 500 1000
−0.2
0
0.2
0.4
0.6
(d)
Figure 2.8:
k
for CMP process with window size of 200. (a): PEM method
withb = 300; (b): PLR method withb = 300; (c): PEM method withb = 600;
(d): PLR method withb = 600.
The process is simplified to test the proposed adaptive SISO run-to-
run control scheme. We assume only the removing ratey
1
is under control
and plate speed u
1
is adjustable; other three inputs are fixed at u
2
=2,
u
3
= 2, andu
4
=0:1. Moreover, the trends are removed so that the process
is controllable by EWMA-based controller. Therefore, the simplified process
is described by
y
1
= 159:3u
1
38:2u
2
+178:9u
3
+24:9u
4
67:2u
1
u
2
46:2u
2
1
19:2u
2
2
28:9u
2
3
+e
1
(2.20)
wheree
1
N(0; 60
2
).
According to the results of the previous example, the past window
33
Table 2.1: Control performance of the controllers regulating CMP process
Output variance
EWMA withb = 600 and fixed = 0:1 3951.5
EWMA withb = 600 and fixed = 0:05 3913.0
Adaptive PEM-based EWMA withb = 300 4008.0
Adaptive PEM-based EWMA withb = 600 3961.9
Adaptive PLR-based EWMA withb = 300 4235.6
Adaptive PLR-based EWMA withb = 600 3972.3
size is selected as 200. Intuitively, the process is quadratic with respect to
the inputs, and the additional disturbance is Gaussian. The process can be
considered linear when it is close to its operating point. Therefore, theoret-
ically, the controller gain needs to be zero after the steady-state operating
point is found.
Fig. 2.8 gives the adaptive EMWA weights by PEM and PLR methods
with different process gain assumed. The initial guess is
^

0
= 0:5. The
results show that the PEM-based estimation performs better than the PLR-
based estimation, which even yields some negative
^

k
values leading to
destabilizing controller. The possible reason, as argued before, is that the
model residual y
k
bu
k
does not follow IMA(1,1) model. It can also be
observed that assuming a larger (linear) process gainb improves the results
by PLR since the controller gains at initial stage are more conservative. But
note that a too largeb may cause instability.
Table 2.1 summarizes the control performance of the controllers. The
34
PEM-based adaptive run-to-run controller leads to a control performance
similar to the nearly-optimal EWMA controller. For PLR-based method,
its performance is not as satisfactory as PEM-based method. This is partly
due to the fact that the process disturbance including nonlinearities and
gain mismatch does not fit IMA(1,1) model very well. Therefore, for non-
linear time-invariant processes, adaptive run-to-run control may be able to
achieve similar performance as optimal controllers.
2.6 Summary
In this chapter, we present an EWMA-based run-to-run controllers
adapting to process disturbances only. The EWMA weight is related to the
disturbance parameter that is estimated online. Two estimation approaches
are suggested: PLR and PEM. Convergence analysis of the proposed adap-
tive scheme is performed. It is proved that the estimated controller gain
converges to the optimal controller gain, which is an equilibrium point. This
point is shown to be stable under almost all circumstances. In the simula-
tion, choice of identification method and window size is discussed. The re-
sults demonstrate that the adaptive controller outperforms the fixed-weight
controllers under the case of slow time-varying disturbances and nonlinear
processes for most of the time.
35
Chapter 3
MPC Disturbance Modeling and Adaptation for an
Industrial Process
3.1 Introduction
Model predictive control (MPC) is one of the advanced process con-
trol techniques regulating industrial MIMO processes. It calculates optimal
values for manipulated variables based on future behaviors predicted by
a model it maintains. When there are slow time-varying unmodeled dis-
turbances or model mismatches, predictions of controlled variables made
by the model may not be accurate. The bias between the current measured
and estimated outputs is usually added to future predicted values to reduce
prediction errors. It is known that this strategy could minimize prediction
errors only when disturbance is a sequence of integrating white noise. An
alternative approach would rather consider these disturbances as a result
of model mismatch, and then applies robust MPC [5]. A limitation of the
method is that unmeasured disturbances often bear dynamics; the uncer-
tainty range could be very large. Recently, some adaptive MPC schemes
have been proposed. Fukushima et al. [20] combines robust MPC with on-
line parameter identification. Dougherty and Cooper [16] construct a con-
36
troller composed of multiple linear MPCs adapting to changing operating
points of nonlinear systems. However, most adaptive MPC methods focus
on entire MPC model including both process and disturbance models. Ex-
cessive adaptation of the entire model may not improve the model quality
at all times. In contrast, with limited amount of data, the uncertainties of
estimated model increases with the complexity of model, leading to a poor
model.
In this chapter, we extended the disturbance adaptation technique
for SISO model-based controllers to MIMO MPCs, where process models
are usually more complicated than disturbance dynamics. Therefore, in
our approach, only the disturbance model is adapted without the need for
external perturbations. The disturbance model is updated by re-estimating
the Kalman gain to give better predictions of future disturbances. MPC
then takes advantage of these predictions to calculate optimal manipulated
variables. Implementation of the proposed MPC strategy in commercial
MPC packages such as DMC is discussed. In this case the objective is the
ease of implementation rather than trying to be the most general form of
formulation. Performance of the proposed MPC strategy will be studied by
simulations based on the Wood-Berry distillation model.
37
3.2 Oset-free Model Predictive Control and Disturbance
Model
3.2.1 Internal model and disturbance
We assume that the MPC uses the following discrete linear time-
invariant model in the state space form:
x(k + 1) =Ax(k) +Bu(k) (3.1a)
y(k) =Cx(k) (3.1b)
wherex(k)2R
n
is the state of the system,u(k)2R
nu
is the input or manip-
ulated variable, andy(k)2R
ny
is the output or controlled variable.
If there is no plant-model mismatch or other process disturbances,
the MPC based on process model (3.1) is able to achieve optimal control
in open loop. However, this is generally not the case. In order to remove
offset caused by disturbance, the process model is usually augmented by an
additional disturbance model with integrators [63]:

x(k + 1)
d(k + 1)

=

A B
d
0 I

x(k)
d(k)

+

B
0

u(k) +w(k) (3.2a)
y(k) =

C C
d

x(k)
d(k)

+v(k) (3.2b)
whereB
d
2R
nn
d
,C
d
2R
nyn
d
, andx
d
2R
n
d
denotes the disturbance state.
Process and measurement noise vectors w(k) and v(k) observe Gaussian
distribution with zero mean. The augmented system is detectable if and
38
only if the following condition holds [62]
rank

IA B
d
C C
d

=n +n
d
(3.3)
The model (3.2) is called the internal model of MPC.
As long as (3.3) is satisfied, the corresponding Kalman filter for (3.2)
can be written as

^ x(k + 1jk)
^
d(k + 1jk)

=

A B
d
0 I

^ x(kjk 1)
^
d(kjk 1)

+

B
0

u(k)
+

L
x
L
d

h
y(k)C^ x(kjk 1)C
d
^
d(kjk 1)
i
(3.4)
in which filter gain is partitioned asL
x
2R
n
andL
d
2R
n
d
for plant and dis-
turbance states, respectively. Theoretically, the filter gain [L
x
; L
d
]
T
should
be steady-state Kalman gain calculated from discrete algebraic Riccati equa-
tion.
Difference between the measured output and predicted output by
process model is usually referred as internal model error:
^
d(k) =y(k)C^ x(kjk 1) (3.5)
It is caused by process model mismatch and/or process and measurement
noise.
3.2.2 Dynamic matrix control
DMC uses the input-output model rather than the state-space form:
y(k) =
N
X
i=1
S
i
u(ki) +y
0
+d(k)
39
wherey
0
is the initial output, u(j) =u(j)u(j1),d(t) is the disturbance,
andS
i
is theith step response matrix
S
i
=
2
6
6
6
4
s
1;1;i
s
1;2;i
s
1;nu;i
s
2;1;i
s
2;2;i
s
2;nu;i
.
.
.
.
.
.
.
.
.
.
.
.
s
ny;1;i
s
ny;2;i
s
ny;nu;i
3
7
7
7
5
(3.6)
where s
j;k;i
denotes the ith step response coefficient from the kth input to
thejth output. The step response coefficient matrixS
i
can be derived from
the state-space model (3.1):
S
i
=
i1
X
l=0
CA
l
B; i = 1;:::;N
As many other industrial MPCs do, DMC uses the step-like distur-
bance model, which is equivalent to
B
d
= 0; C
d
=I (3.7)
in (3.2) and (3.4) in state-space form. Such disturbance model assumes that
the disturbance is a constant for all the future outputs. The filter gain in
DMC is set to be a constant value
L
x
= 0; L
d
=I (3.8)
rather than the Kalman gain computed from Riccati equation.
However, some industrial controllers, e.g., DMCplus, have the flexi-
bility for users to apply part of the prediction error
e(k) =y(k)C^ x(kjk 1)C
d
^
d(kjk 1); (3.9)
40
to the future predictions. The equivalent Kalman gain in state-space form is
L
x
= 0; L
d
=K
d
: (3.10)
Thus, the state estimator of the augmented system (3.2) can be rewritten as

^ x(k + 1jk)
^
d(k + 1jk)

=

A 0
0 I

^ x(kjk 1)
^
d(kjk 1)

+

B
0

u(k)
+

0
K
d

h
y(k)C^ x(kjk 1)
^
d(kjk 1)
i
(3.11)
The one-step-ahead prediction of the output can be obtained as
^ y(k + 1jk) =C^ x(k + 1jk) +
^
d(k + 1jk) (3.12)
It is inferred from (3.11) that the plant state is independent of one-step-
ahead prediction errors. The process model mismatch is lumped into dis-
turbance states, and feedback action is performed by disturbances.
Combining (3.9), (3.11) and (3.12), we can obtain the input-output
model for theith output:
y
i
(k) =G
i
(q)u(k) +
1 (1K
d;i
)q
1
1q
1
e
i
(k) +y
0;i
(3.13)
where G
i
(q) =
P
N
i=1
S
i
q
i
is the transfer function model, and K
d;i
is the
ith diagonal element of K
d
. Note that K
d
2 R
nyny
is diagonal due to the
assumption of output disturbance model for each CV . The second term of
41
(3.13) suggests that the prediction error filtering in industrial MPC is equiv-
alent to the assumption of IMA(1,1) disturbance model. Especially, the tra-
ditional step-like disturbance model, where K
d
= I, is suitable for inte-
grated white noise.
Although the IMA(1,1) model is good enough under many cases, it
is still desirable to have a more general disturbance model as follows,

x(k + 1)
d(k + 1)

=

A 0
0 A
d

x(k)
d(k)

+

B
0

u(k) +w(k) (3.14a)
y(k) =

C C
d

x(k)
d(k)

+v(k) (3.14b)
Compared with (3.2), (3.14) does not require disturbance share the same
modes with plant model. The corresponding state estimator is

^ x(k + 1jk)
^
d(k + 1jk)

=

A 0
0 A
d

^ x(kjk 1)
^
d(kjk 1)

+

B
0

u(k)
+

L
x
L
d

h
y(k)C^ x(kjk 1)
^
d(kjk 1)
i
(3.15)
This type of disturbance is the same as Box-Jenkins model in the input-
output form.
Thep-steps-ahead DMC predictions can be computed by
^ y(k +pjk) =y(k) +
p
X
i=1
S
i
u(k +pi) +
^
d(k + 1jk)C
d
d(k) (3.16)
where
^
d(k + 1jk) depends on the disturbance model used. Based
on (3.16), DMC determines the optimal future input move sequence
42
u(k);:::; u(k +M), whereM is the control horizon, that tracks reference
signal while satisfying all the constraints.
3.3 Disturbance Model Estimation and Adaptation
3.3.1 IMA(1,1) model
We will first explore a practical method estimating the filter gainK
d
for IMA(1,1) model, which can be implemented in many industrial MPCs
through prediction error filtering. Similar to (3.5), the internal model error
^
d(k) ofith output is
^
d
i
(k) =y
i
(k)G
i
(q)u(k): (3.17)
When the IMA(1,1) disturbance model (3.13) is used, (3.17) becomes
^
d
i
(k) =
1
i
q
1
1q
1
e(k): (3.18)
where
i
= 1K
d;i
. Here, innovation or prediction errore(k) is regarded as
white noise. However, when plant or disturbance model mismatch exists,
this is generally not the case. We will discuss this issue later.
Differentiating (3.18), we obtain

^
d
i
(k) = (1
i
q
1
)e(k) (3.19)
whose right hand side is a moving average (MA) (1) process. Its lag-one
autocorrelation is

i;1
=

i
1 +
2
i
(3.20)
43
Therefore, an estimate of
i
can be
^

i
=
1
q
1 4
2
i;1
2
i;1
(3.21)
and the Kalman gain to use is then
K
d;i
= 1
^

i
The estimation is unbiased ife(t) is white [7]. It is implied in (3.21)
that
i;1
should be between0:5 and +0:5; if it falls out of this region, it can
be concluded that IMA(1,1) is unsuitable for the disturbance. Furthermore,
a positive indicates a negative
^

i
. This results in an under-damped filter,
which may cause oscillation problems in practice. Thus, it is undesirable to
apply IMA(1,1) model when lag-one autocorrelation is positive.
3.3.2 General disturbance model
When it is insufficient to use IMA(1,1) to describe disturbance char-
acteristics, a general disturbance model is required. The advantages of hav-
ing a general model for disturbance can be summarized as
the restriction of diagonal structure of output disturbance is removed;
model order can be determined by subspace methods;
identification algorithm is more efficient than I/O model based
method such as prediction error method.
44
To identify the disturbance model in state-space form, the system
equation of disturbance states are extracted from (3.14):
x
d
(k + 1) =A
d
x
d
(k) +L
d
e(k)
d(k) =C
d
x
d
(k) +e(k)
By stacking thed(k), we can define
d
f
(k) =
2
6
6
6
4
d(k)
d(k + 1)
.
.
.
d(k +f 1)
3
7
7
7
5
;d
p
(k) =
2
6
6
6
4
d(k 1)
d(k 2)
.
.
.
d(kp)
3
7
7
7
5
The future disturbances can be written as
d
f
(k) =

H
fp
y
p
(k) +

H
f
y
f
(k) +e
f
(k)
To enforce causality,

H
f
must have the lower triangular form:

H
f
=
2
6
6
6
6
6
4
0 0 0 0

h
1
0 0 0

h
2

h
1
0 0
.
.
.
.
.
.
.
.
.
.
.
.

h
f1

h
f2

h
1
0
3
7
7
7
7
7
5
where

h
i
is obtained from high-order ARX (HOARX) model
d
k
=
p
X
i=1

h
i
d
ki
+e
k
Eliminating the effects by future disturbances
~ y
f
(k) =y
f
(k)

H
f
y
f
(k)
=

H
fp
y
p
(k) +e
f
(k) (3.22)
45
theH
fp
matrix can be estimated using lease squares non-parametrically. The
non-parametric model (3.22) can be used to do predictions of future dis-
turbance data based on past window of disturbances. However, one can
go further to gain the system matricesA
d
;L
d
;C
d
so that predictions can be
made in the state-space framework.
The state vector can be represented by past disturbance data:
^ x
k
=

L
p
d
p
(k)
We have
~
d
f
(k) =

f

L
p
d
p
(k) +e
f
(k) (3.23)
It can be seen from (3.22) and (3.23) thatH
fp
=

f

L
p
. To estimate

f
and

L
p
,
it is recommended to perform canonical variate analysis (CVA) by Larimore.
After obtaining
R
fp
=
1
N
N
X
k=1
~ y
f
(k)z
T
p
(k)
R
ff
=
1
N
N
X
k=1
~ y
f
(k)~ y
T
f
(k)
R
pp
=
1
N
N
X
k=1
z
p
(k)z
T
p
(k)
singular value decomposition can be done on
R

1
2
ff
R
fp
R

1
2
pp
=U
n
D
n
V
T
n
+
n
46
Thus, the CVA estimates give the CVA estimates are

CVA
f
=R
1
2
ff
U
n
D
1
2
n

L
CVA
p
=D
1
2
n
V
T
n
R

1
2
pp
:
There are several ways to get system matrices A
d
;L
d
;C
d
for distur-
bance model. Two of the most commonly used methods are
1. Reconstruct state sequencef^ x(k)g using past disturbance data d
p
(k)
along with

L
p
. Then getA
d
;L
d
;C
d
by least squares.
2. ObtainC
d
andA
d
from

f
, and then estimateL
d
by linear regression.
With the system matrices identified, we can use the general disturbance
model (3.11) for MPC control.
3.3.3 Disturbance Adaptation for MPC
When there is no model mismatch and process and disturbance are
time invariant, MPC is able to achieve an optimal estimation of the states
and hence to provide best control actions. However, under the existence of
plant-model mismatch and/or slow time-varying process and disturbance,
one may have to adapt the model to have a better control performance.
Adaptive modeling usually includes both plant and disturbance
model adaptation. However, adapting disturbance model only is recom-
mended for the following reasons:
47
o
y(k)
e(k)
u(k)
+
+
MPC
_
+
r(k)
Figure 3.1: Schematic diagram of MPC disturbance adaptation.
1. most of process models of MIMO processes are very complicated;
adaptation of such process models in either input-output or state-
space formulation requires a larger number of parameters, and thus
considerable amount of closed-loop data are needed for identification;
2. some industrial MPCs like DMC is based on non-parametric model,
e.g., FIR model, which increases uncertainties in the estimated models;
and
3. disturbance model can partially compensate mismatch in plant model.
The steps of disturbance adaptation can be generalized in three steps,
as shown in Fig. 3.1. First, the frequency of adaptation needs to be deter-
48
mined. Adapting too frequently is unnecessary and may cause instability
due to the estimation error. While adaptation can be performed at a fixed
frequency, it is suggested adaptation be determined by the control perfor-
mance. That is, when control performance is satisfactory, the MPC model
is left unchanged; when control performance is degraded significantly, the
adaptive mechanism becomes active. Second, the disturbance model is es-
timated using closed-loop data. It depends on the disturbance model struc-
ture of MPC. Third, the identified disturbance model is converted to the
form that can be applied to the MPC. For example, the IMA(1,1) disturbance
model is translated to the prediction error filtering.
3.4 Case Studies
3.4.1 Simulated distillation column
The disturbance adaptation method is demonstrated by a distillation
column. Due to the interactions and time delay, distillation columns are not
easy to control. Wood and Berry [87] model an methanol/water column. Its
transfer function model is represented by
G

(s) =
2
6
6
4
12:8e
s
16:7s + 1
18:9e
3s
21:0s + 1
6:6e
7s
10:9s + 1
19:4e
3s
14:4s + 1
3
7
7
5
:
With a sampling time of 1 minute, it is discretized as
G

(q) =
2
6
4
q
2
0:7440
1 0:9419q
1
q
4
0:8789
1 0:9535q
1
q
8
0:5786
1 0:9123q
1
q
4
1:3020
1 0:9329q
1
3
7
5
49

y
1;k
y
2;k

=G

(q)

u
1;k
u
2;k

where inputu
1
is reflux rate in lb/min, inputu
2
is steam flow rate in lb/min,
outputy
1
is composition of top product of column in mol %, and outputy
2
is composition of bottom product of column in mol %. The setpoints are
chosen as y
1
= 96:25% methanol and y
2
= 1% methanol.
We simulated several different MPCs for this process in MATLAB.
The DMC formulation is used. The prediction and control horizon are se-
lected to be 200 and 8, respectively. The output errors are penalized in the
objective function by the weight matrixQ = diag([1; 196]). No other terms
are included in the objective so that minimum variance control (MVC) is
achieved. The reason why MVC law is applied is that the (weighted) out-
put variance can serve as the control performance metrics. A smaller num-
ber of control performance implies better control performance. We perform
Monte Carlo simulation for 30 different noise source sequences to better
demonstrate the results.
3.4.1.1 Case I: No plant model mismatch, IMA(1,1) disturbance
In this case, it is assumed that the process model is known exactly.
The following model is assumed for stochastic disturbances, which is
adopted from [34]. Slight changes are made for the disturbances to make
50
them diagonal, i.e., they are added to outputs individually:
d
k
=
2
6
6
4
(1 0:5q
1
)
1q
1
0:14(1 0:7q
1
)
1q
1
3
7
7
5
e
k
(3.24)
wheree(k)2R
2
is independent white noise with unit variance.
Table 3.1: Output Variance of Different MPCs
DMC DMC with IMA(1,1) mdl. DMC with general mdl.
Case I 8.1933 6.1080 6.3491
Case II 13.8279 13.4074 13.6102
Case III 11.0291 7.7190 6.6566
We compare the DMC to the MPCs with disturbance adaptation.
The control performance in terms of output variance is summarized in
Table 3.1. We can see that the step-like disturbance model used by DMC
yields the largest variance, meaning the worst control performance. The
performance is improved by having a diagonal IMA(1,1) disturbance
model, whose Kalman gain is identified as
L
d
=

K
d;1
0
0 K
d;2

=

0:5375 0
0 0:2736

(3.25)
for one of the noise sources. This is quite close to the true Kalman gain of
(3.24)K
d;1
= 0:5;K
d;2
= 0:3, and thus reduction of CV variance is expected.
When general disturbance model is used, a similar control performance is
observed. The output variance of results by general disturbance model is
even slightly larger than the output variance of results by IMA(1,1) model.
51
It is partly due to the fact that the disturbance in this case fits the diago-
nal structure assumed in IMA(1,1) model. The general disturbance model,
which is in state-space form, also accounts for off-diagonal elements and is
overfit in this case.
3.4.1.2 Case II: No plant model mismatch, higher-order disturbance
Some dynamics are added to the disturbance (3.24) so that it follows
a higher order model than IMA(1,1):
d
k
=
2
6
6
4
(1 0:5q
1
)
1q
1
1 0:6q
1
1 0:7q
1
0:14(1 0:7q
1
)
1q
1
1 0:2q
1
1 0:8q
1
3
7
7
5
e
k
(3.26)
The control performance is compared in Table 3.1. It can be seen
that the diagonal IMA(1,1) structure is able to handle higher-order distur-
bance dynamics well. One of the Kalman gain from identified models is
K
d;1
= 0:2864;K
d;2
= 0:8506. This type of reduced-order disturbance mod-
eling combines the additional poles and zeros to the only zero, which is
Kalman gain in IMA(1,1) model. Moreover, we observe the similar results
as in Case I that the performance of controller with general disturbance
model is poorer than the performance of controller with diagonal model.
This is the result of higher-order diagonal disturbance added; the general
model in state-space model attempts to model complex dynamics on the di-
agonal, and off-diagonal elements are also modeled in higher order by the
52
identification. Therefore, the parametric method using IMA(1,1) model out-
performs the non-parametric method using a general model, although the
model order of the parametric method is insufficient.
3.4.1.3 Case III: Plant model mismatch, IMA(1,1) disturbance
In this case, how the proposed approaches perform under plant
model mismatch is explored. This is a more realistic case since mismatch
between the MPC model and the process is inevitable in the plant. It is
assumed that the MPC model is mistaken by
G(q) =
2
6
4
1:2q
2
0:7665
1 0:9419q
1
q
4
0:9000
1 0:9535q
1
0:8q
8
0:6055
1 0:9123q
3
0:9q
4
1:3472
1 0:9329q
1
3
7
5
The disturbance to be added is (3.24) which is the same as Case I.
Table 3.2: Effects of Disturbance Model Adaptation
DMC DMC with IMA(1,1) mdl. DMC with general mdl.
Iter. 1 11.3653 8.1029 7.2630
Iter. 2 11.3653 8.0904 7.8331
Iter. 3 11.3653 8.0898 7.2965
Iter. 4 11.3653 8.0898 7.3389
Iter. 5 11.3653 8.0898 7.3244
As indicated by (3.17), with the existence of plant model mismatch,
the internal model error is
^
d(k) = Gu(k)+d(k). Therefore, the disturbance
includes the plant model mismatch term in addition to noises. Adaptation
of disturbance model is needed; in this LTI system case, disturbance model
53
has to be updated constantly until the control performance reaches a fixed
point. The convergence results of disturbance adaptation are summarized
in Table 3.2. Since DMC always assumes step disturbance, only the last
two columns are focused. It is found that the control performance of meth-
ods based on either IMA(1,1) model or general model converge to a fixed
value. Since IMA(1,1) model is parametric, it converges more quickly than
the general disturbance model.
The results of Case III are also listed in Table 3.1. In each of the
Monte-Carlo simulation, five iterations of adaptations are executed. Only
the runs in which disturbance model converges in terms of performance
metrics are considered. It is observed that DMC using general disturbance
model performs better than DMC with IMA(1,1) model. The reason is
that under plant model mismatch the internal model error
^
d(k) consists of
G(q)u(k), suggesting that off-diagonal elements should be included. And
the state-space form is capable of modeling a full matrix of disturbance
dynamics. Therefore, general disturbance model is more appropriate in
this case to compensate the effects by model mismatch.
Now we summarize all three cases. The MPC with IMA(1,1) distur-
bance has a diagonal structure, which is easy to be implemented in com-
mercial MPCs. It utilizes a parametric disturbance model. This enables
faster convergence and less computation. Although under the case of model
54
mismatch it is inferior to MPC with general disturbance model, significant
improvements are observed as compared to DMC with step disturbance
model, and the mismatch is partly compensated by the IMA(1,1) model.
Therefore, disturbance adaptation using IMA(1,1) is recommended in prac-
tice if the model mismatch is not significant.
3.4.2 An industrial air separation process
An air separation process of Praxair is studied to evaluate the pro-
posed disturbance modeling technique against the traditional DMC with
step-like disturbance model.
Figure 3.2: Simplified diagram of an industrial air separation process [84].
55
3.4.2.1 Process Description
The industrial case in this section is a large air separation plant of
Praxair Inc. Air separation is a process that separates air into its compo-
nents, mostly nitrogen and oxygen. In most large commercial air separa-
tion plants that yield high purity products, the distillation and cryogenic
processes are used. A schematic diagram is shown in Fig. 3.2. First, atmo-
spheric air is filtered to remove dirts and dusts, compressed and purified
to eliminate residual dirt. The air stream is then fed into the lower col-
umn (high pressure column) to produce high purity nitrogen and oxygen-
enriched liquid air, which is passed into the crude argon condenser before
subcooling in the nitrogen superheaters. The stream then enters the up-
per column (low pressure column) for final separation. A stream of nitro-
gen vapor withdrawn from the top of the lower column is liquefied in the
main condenser. The liquid nitrogen from the condenser is divided into
two streams: one returned to the top of the lower column and the other to
be subcooled in the nitrogen superheaters.
3.4.2.2 MATLAB simulation
The process is simulated in the MATLAB with Simulink first. A sub-
system consisting of 7 manipulated variables and 8 controlled variables are
included in the simulation. The purpose of disturbance adaptation is to im-
prove MPC predictions, i.e., the prediction error of MPC is expected to be
56
reduced by adapting disturbance models. However, the MPC KPI is not
a direct measurement of model quality. Other factors such as constraints
may also affect the KPI. For example, for a constrained minimum variance
control, an MPC with good model may make MV hit constraint frequently,
causing increase in output variance; but an MPC with smoothing may make
MV hit constraint less, although it is not optimal. Therefore, the advantage
of MATLAB simulation is that we can evaluate MPC with no constraints,
which is impractical in real plants, and hence we can simply use LQ metric
to evaluate control performance.
The simulated plant and disturbance are identified from the past rou-
tine operation data. We also identify the transfer function model from the
DMC model so that MATLAB MPC toolbox can be used. The simulated
plant model and MPC model are not the same, thus indicating the existence
of plant model mismatch. The control results are shown in Fig. 3.3. The LQ
objective of DMC with step disturbance is 1.6787, while the LQ objective of
DMC with IMA(1,1) disturbance is 1.5402. We observe 8.25% improvement
in terms of control peroformance. This verifies that the IMA(1,1) model is
able to characterize disturbances and to partially compensate DMC model
errors.
57
0 500 1000 1500 2000
−1
0
1
2
3
4
0 500 1000 1500 2000
−1
0
1
2
3
4
Figure 3.3: Control results by MATLAB simulation. The CVs are scaled. The
upper one is by DMC with step disturbance model and the lower one is by
DMC with IMA(1,1) model.
3.4.2.3 Plant data study
The proposed disturbance adaptation method is applied to the real
plant of Praxair. As aforementioned, factors other than model quality may
contribute to the MPC KPI. Therefore, an alternative performance metric
has to be adopted to assess improvements of the proposed method. We
choose the model quality index (MQI) [79] as the control performance met-
ric.
In the plant, a commercial MPC software is applied to control a sub-
58
process contains 22 manipulated variables and 28 controlled variables. Data
are sampled every 2 minutes since the MPC uses the same sampling time.
The step responses of the model used by MPC are identified by the engi-
neers and cannot be changed when operating.
Table 3.3: Model Quality Index of Results of the Industrial Air Separation
Process
Kalman Gain MQI of MPC1 MQI of MPC2 Improvement
CV3 0.6842 0.1656 0.2027 22.40%
CV9 0.1037 0.1285 0.0899 -30.04%
CV11 0.8981 0.1730 0.4148 139.77%
CV12 0.5410 0.1671 0.1484 -11.19%
CV13 0.5905 0.1745 0.2417 38.51%
CV14 0.0313 0.1701 0.2039 19.87%
CV15 0.7859 0.2281 0.3623 58.83%
CV16 0.5860 0.1592 0.1876 17.84%
CV23 0.7668 0.1883 0.3150 67.29%
CV24 0.4162 0.5259 0.4817 -8.40%
CV25 0.1014 0.5910 0.5720 -3.21%
A larger MQI value implies a better model.
Disturbance modeling algorithm is implemented first to build the
IMA(1,1) model for each CV based on the closed-loop data. The modeling
results are shown in Table 3.3 in terms of the Kalman gain. A different
MPC can be generated with the same plant model but with the identified
disturbance model.
Since the disturbance assumed here has a diagonal structure, the
MIMO process can be decomposed to 28 MISO processes. The two data
segments are collected under different operating conditions like setpoints
59
and active constraints. Fortunately, the MQI metric is able to exclude these
factors. Table 3.3 shows the model quality assessment results obtained from
the closed-loop data collected. It can be seen that the models are improved
significantly for most of the CVs after the disturbance adaptation method
is applied to the MPC. Although there exists performance degradation for
some CVs, nevertheless, the degradation is acceptable. We can conclude
that the predictions made by the MPC with adaptive disturbance model
are more accurate, and thus a better control performance is achieved.
3.5 Summary
The disturbance adaptation technique is extended to MIMO pro-
cesses in this chapter. The process and disturbance model of MPC is
introduced. Several methods identifying the disturbance are discussed.
For industrial MPCs, the IMA(1,1) disturbance model is recommended.
It can be applied to commerical MPC packages through prediction error
filtering. Then, the general steps for adapting disturbance are provided.
Simulation results show that as compared with traditional MPC the control
performance is improved by implementation of disturbance adaptation.
It is also found that process model mismatch can be compensated to
some extent by the proposed adaptive technique. The industrial example
confirms that the method is also applicable to nonlinear processes and
considerable improvements are observed.
60
Chapter 4
Improving Industrial MPC Performance with Data-
Driven Disturbance Modeling
4.1 Introduction
Model predictive control (MPC) is a widely used control technique
in process industry. One of the first MPC implementations was developed
by Cutler and Ramaker [11] as dynamic matrix control (DMC). The DMC
algorithm makes predictions of future outputs using a step response model.
The DMC approach was then extended to quadratic DMC (QDMC) [24]
to handle constraints in a quadratic programming (QP) problem. Recent
survey shows that DMC based MPC products still play an important role in
industrial applications [68].
As a model-based control method, MPC relies heavily on the accu-
racy of the model. Either plant model mismatch or unmodeled disturbances
may degrade the control performance. DMC and QDMC uses a step-like
disturbance model to achieve offset-free control. Muske and Badgwell [55]
proposed a block diagonal disturbance model to include both state and out-
put disturbances. Pannocchia and Rawlings [63] provided a more general
disturbance model without special structures. However, the filter gain of
61
the state estimator is predetermined. To avoid the need to specify the dis-
turbances at the input or the output side, the equivalence between different
disturbance models are shown in [70]. Recently, Xu et al. [88] proposed an
adaptive disturbance modeling method, but it imposes a special structure
on the Kalman filter.
Subspace identification (SID) has drawn tremendous interests in the
last two decades. Several representative algorithms have been developed,
including CVA [45], N4SID [59] and MOESP [83]. The unifying theorem by
Van Overschee and De Moor [60] provides a framework in which these algo-
rithms can be interpreted as a singular value decomposition of a weighted
matrix. These algorithms have been widely used for the open-loop systems
because of the advantage of simple parametrization. However, due to the
correlation between the current input and the past output caused by feed-
back, the application of most SID methods to the closed-loop data typically
gives biased estimation. To address this problem, The SID methods for the
closed-loop system were investigated and developed in the last decade. Qin
and Ljung developed the closed-loop SID algorithm by innovation estima-
tion [67]. In [35], the pre-estimation of the Markov parameters are proposed
to avoid the non-causal projection. The statistical consistency is studies in
[10].
Most closed-loop subspace identification methods involve three
62
stages [69]. The first step is to identify a high order ARX model from
properly collected input and output data. Providing the observer Markov
parameter estimation for high order ARX model, this step by itself is a
non-parametric procedure.
Motivated by DMC and subspace estimation of non-parametric
models, we propose a new DMC technique with a disturbance model
identified by the subspace method using closed loop data. The observer
Markov parameters of the plant and disturbance model are obtained.
Then, system Markov parameters of the disturbance model is computed
recursively. Moreover, we propose a new recursive relationship utilizing
the plant model information to better estimate the disturbance model. A
more accurate DMC prediction is then made possible by the identified
disturbance model. This scheme can be implemented adaptively as needed.
4.2 Dynamic Matrix Control
4.2.1 Plant model in step response form
DMC algorithm uses a step response model of the plant to predict
future outputs. Denote y
t
2 R
ny
and u
t
2 R
nu
as the output and input
vector of a MIMO plant respectively. The relationship inputs and outputs
can be written as [24]
y(t) =
N
X
i=1
S
i
u(ti) +y
0
+d(t) (4.1)
63
where y
0
is the output initial condition, d(t) denotes unmodelled distur-
bances u(j) =u(j)u(j 1),N is the number of steps for plant to reach
steady-state, andS
i
is theith step response matrix
S
i
=
0
B
B
B
@
s
1;1;i
s
1;2;i
s
1;nu;i
s
2;1;i
s
2;2;i
s
2;nu;i
.
.
.
.
.
.
.
.
.
.
.
.
s
ny;1;i
s
ny;2;i
s
ny;nu;i
1
C
C
C
A
(4.2)
in whichs
j;k;i
represents theith step response coefficient from thekth input
to thejth output. Thek-step-ahead prediction can be made by
^ y(t +kjt) =y(t) +
k
X
i=1
S
i
u(t +ki) +
^
d(t +kjt) (4.3)
where ^ y(t +kjt) and
^
d(t +kjt) denotes the prediction ofy(t +k) and unmod-
elled disturbances given the information up to timet , respectively,
The step response model and finite impulse response (FIR) model
can be used interchangeably:
0
B
B
B
@
S
1
S
2
.
.
.
S
P
1
C
C
C
A
=
0
B
B
B
@
1
1 1
.
.
.
.
.
.
.
.
.
1 1
1
C
C
C
A
0
B
B
B
@
G
p
1
G
p
2
.
.
.
G
p
P
1
C
C
C
A
(4.4)
whereG
p
1
;:::;G
p
P
are the impulse response coefficients, also known as sys-
tem Markov parameters, of the plant model.
64
4.2.2 DMC predictions and oset-free control
From (4.3), the predicted outputs can be rewritten as
^ y(t +kjt) =
k
X
i=1
S
i
u(t +ki) +
P
X
i=k+1
S
i
u(t +ki)
+y
0
+
^
d(t +kjt) (4.5)
whereP indicates the prediction horizon. Define
y

(t +k) =y
0
+
P
X
i=k+1
S
i
u(t +ki)
which is the effect of past inputs, and stack ^ y(t +ijt); (i = 1;:::;P ) using
(4.5):
0
B
@
^ y(t + 1jt)
.
.
.
^ y(t +Pjt)
1
C
A
=
0
B
@
y

(t + 1)
.
.
.
y

(t +P )
1
C
A
+A
M
P
0
B
@
u(t)
.
.
.
u(t +M 1)
1
C
A
+
0
B
@
^
d(t + 1)
.
.
.
^
d(t +P )
1
C
A
(4.6)
whereM is the control horizon and
A
M
P
=
0
B
B
B
B
B
B
B
B
B
B
B
B
@
S
1
0 0
S
2
S
1
0
.
.
.
.
.
.
.
.
.
S
M
S
M1
S
1
.
.
.
.
.
.
.
.
.
S
P
S
P1
S
PM+1
.
.
.
.
.
.
.
.
.
S
P
S
P
S
P
1
C
C
C
C
C
C
C
C
C
C
C
C
A
(4.7)
is called the dynamic matrix of the system.
DMC assumes a step disturbance model and estimate the step size
using the difference between the current measured and estimated output,
65
i.e.,
^
d(t +ijt) =y(t) ^ y(tjt 1); i = 1;:::;P (4.8)
This disturbance model is integrating and is able to achieve offset-free con-
trol for DMC.
4.2.3 Plant model in state-space form
We use the state-space form of plant model which will be useful in
the derivations later. Assume the plant is described by
8
<
:
x
p
(t + 1) =A
p
x
p
(t) +B
p
u(t)
y
p
(t) =C
p
x
p
(t)
(4.9)
which is considered to be deterministic and the disturbance model is writ-
ten in the following innovation form of Kalman filter:
8
<
:
^ x
d
(t + 1) =A
d
^ x
d
(t) +K
d
e(t)
y
d
(t) =C
d
^ x
d
(t) + e(t)
(4.10)
which is considered to be stochastic. When the step-like disturbance model
in DMC is assumed,A
d
= I
ny
,K
d
= I
ny
, andC
d
= I
ny
; this type of distur-
bance model corresponds to the case that disturbance is integrated white
noise.
By augmenting (4.9) and (4.10), we obtain

^ x(t + 1) =A
aug
^ x(t) +B
aug
u(t) +K
aug
e(t)
y(t) =C
aug
^ x(t) +e(t)
(4.11)
66
where
A
aug
=

A
p
0
0 A
d

; B
aug
=

B
p
0

C
aug
=

C
p
C
d

; K
aug
=

0 K
d

T
:
(4.12)
The subscript aug refers to augmented system.
Clearly, the augmented Kalman gain with the structure K
aug
=

0 K
d

T
is insufficient when there exists plant model mismatch or distur-
bance model mismatch. With the help of closed loop data and subspace
identification, it is possible to estimate a complex disturbance model to
improve the model prediction and control performance. Instead, a full
Kalman gain K
aug
=

K
p
K
d

T
will be used. We propose a method to
estimate disturbance model from closed-loop data and a revised form of
DMC that includes disturbance model.
4.3 Subspace Identication of a Disturbance Model
4.3.1 Identication of the observer Markov parameters
In this section, we develop a method for identifying the disturbance
model in the closed-loop subspace identification framework. Most closed-
loop subspace identification methods involves three stages:
Stage 1: Perform high-order ARX to give observer Markov parameter
estimates, which is a non-parametric procedure.
Stage 2: Perform SVD on the Hankel matrix or weighted Hankel ma-
trix, leading to extended observability matrix

f
and extended controllabil-
67
ity

L
p
estimates, which by themselves are reduced order (n) non-parametric
models to estimate the state and predict the output.
Stage 3: Perform least squares to estimate (

A;B;C;K) from

f
and

L
p
leading to state-space parametric models.
The Stage 1 of the closed-loop subspace identification gives the ob-
server Markov parameter estimation in the least-square sense. Collect the
sequence of Markov parameters in the following form:

G =

C
aug

A
l1
aug
B
aug
C
aug

A
l2
aug
B
aug
C
aug
B
aug

H =

C
aug

A
l1
aug
K C
aug

A
l2
aug
K C
aug
K

and
G =

C
aug
A
l1
aug
B
aug
C
aug
A
l2
aug
B
aug
C
aug
B
aug

H =

C
aug
A
l1
aug
K C
aug
A
l2
aug
K C
aug
K

where l denotes the size of the past horizon and

A
aug
= A
aug
K
aug
C
aug
.
Without loss of generality, we setH to be monic, i.e.,H
0
=I
ny
. The consis-
tent identification of G and H is to solve the following equation in by the
least-squares,
Y
t
=C
aug

A
l
aug
X
tl
+GU
tljt1
+HY
tljt1
+E
t
(4.13)
where
Y
t
=

y (t) y (t +N 1)

; (4.14)
E
t
=

e (t) e (t +N 1)

; (4.15)
68
andX
tl
=

^ x (tl) ^ x (t +Nl)

is the sequence of the initial states.
We introduce the assumption that

A
aug
is stable andl is chosen to be suffi-
ciently large,

A
l
aug
0. Then we have the following least squares estimation,

^

G
^

H

=Y
t

U
tljt1
Y
tljt1

T

U
tljt1
Y
tljt1

U
tljt1
Y
tljt1

1
(4.16)
where
U
tljt1
=
0
B
@
u (tl) u (tl + 1) u (tl +N 1)
.
.
.
.
.
.
.
.
.
u (t 1) u (t) u (t +N 2)
1
C
A
and
Y
tljt1
=
0
B
@
y (tl) y (tl + 1) y (tl +N 1)
.
.
.
.
.
.
.
.
.
y (t 1) y (t) y (t +N 2)
1
C
A
:
^

G and
^

H are asymptotically unbiased estimate of the true observer Markov
parameters since the innovations term in (4.13) is uncorrelated to the past
input and output data.
4.3.2 System Markov parameters of the disturbance model
In order to apply the estimated system Markov parameters for the
disturbance model of MPC, the estimated observer Markov parameters can
be converted to the system Markov parameters based on the recursive rela-
tionship [37]. Define the estimated observer Markov parameters
^

G =

G
l1

G
l2

G
0

; (4.17)
^

H =

H
l1

H
l2

H
0

: (4.18)
69
4.3.2.1 Using
^

H to estimate H
Knowledge of the estimated observer Markov parameters allows one
to obtain the estimated open-loop Markov parameters by using the recur-
sion [37]
H
i
=

H
i
+
i1
X
j=1

H
j
H
ij1
; i = 0;:::;l 1: (4.19)
4.3.2.2 Using G and
^

G to estimate H
To better utilize the information of plant modelG
p
in DMC, we will
show the recursive relationship betweenH
i
,G and
^

G.
Equation (4.5) provides the knowledge of the system Markov param-
eter of the process model, which areG
p
i
=C
p
A
i1
p
B
p
. Given the structure of
the augmented system (4.12),
G
i
=C
aug
A
i1
aug
B
aug
=

C
p
C
d

A
p
0
0 A
d

i1

B
p
0

=C
p
A
i1
p
B
p
=G
p
i
i = 1;:::;l 1
which implies that the system Markov parameters of the augmented model
equals that of the plant model. Then,

H can be determined by
G
i
=

G
i
+
i1
X
j=0

H
j
G
ij1
; i = 1;:::;l 1: (4.20)
70
Therefore, H can be estimated by using

H from (4.20) along with (4.19).
Note that either method has no approximation; the quality of estimatedH
solely depends on how good the estimated observer Markov parameter are.
From either (4.19) or (4.20), an FIR model of disturbance for DMC can
be built using the estimated system Markov parameters. This procedure is
non-parametric and convenient for use in MPC calculations.
4.4 DMC with Moving Average Disturbance Model
In this section, we will discuss how to implement DMC with the
information of disturbance model. Although in the last section there is no
assumption on A
d
other than that

A
d
is stable, one may still desire to set
A
d
= I
ny
because offset-free control can be achieved by augmenting the
system state with an integrating disturbance model [72].
It is convenient to differentiate input and output data in order to
identify an integrating disturbance model. Converting step response model
of the plant to impulse response model, (4.1) becomes
y(t) =
N
X
i=1
G
p
i
u(ti) +y
0
+
N
X
i=0
H
i
e(ti): (4.21)
Differentiation of (4.21) gives
y(t) =
N
X
i=1
G
p
i
u(ti) +
N
X
i=0
H
i
e(ti): (4.22)
In this model, the differentiated disturbance model H
i
is stable with its
71
integrators cancelled by the differentiation. Thus, we can directly identify
H
i
by applying the method introduced in the last section.
In order to calculate predictions, integrating (4.22), we obtain
y(t) =
N
X
i=1
S
i
u(ti) +y
0
+
N
X
i=0
H
i
e(ti) +d
0
(4.23)
whered
0
is the initial condition for disturbance. Then,k-step-ahead predic-
tion can be expressed as
^ y(t +kjt) =
k
X
i=1
S
i
u(t +ki) +
P
X
i=k+1
S
i
u(t +ki)
+y
0
+
^
d(t +kjt) (4.24)
^
d(t +kjt) =
N
X
i=0
H
i
e(ti) +d
0
: (4.25)
The optimal input moves can be optimized by QP problem with the predic-
tion model (4.24) and (4.25).
Remark 4.1. Unlike traditional DMC,
^
d(t +kjt) usually does not equal
^
d(t +
1jt) fork > 1, unless there is no model mismatch and disturbances can be
modeled as integrated moving average (IMA) (1,1). The reason is that with
the mismatch, the Kalman gain for plant model K
p
is nonzero. Thus, the
system Markov parameters of disturbance model are
H
i
=C
aug
A
i1
aug
K
aug
=

C
p
C
d

A
p
0
0 A
d

i1

K
p
K
d

=C
p
A
i1
p
K
p
+C
d
K
d
(4.26)
72
where the second term remains constant while the first term does not.
Therefore, it can be inferred that the proposed non-parametric disturbance
model is more general than the step-like one used in DMC.
4.5 Simulations
In the simulation, we use a SISO plant described by the following
system matrices.
A =

0:6 0:3
0 0:2

; B =

0
1

C =

1 1

; D = 0:
(4.27)
We generate 10000 data points using an arbitrary MPC in MATLAB.
The setpoint is selected to be 0, and a PRBS signal is added to input signal
to improve closed-loop system identifiability. The size of past horizon in
subspace identification is set as 30.
To compare the proposed control scheme with DMC, we run both
controllers with 1000 steps. The prediction and control horizon are chosen
asP = 15 andM = 4. We put zero weights on input and input move and
unit weight on output so that output variance can be served as performance
metric. The setpoint of reference signal is always 0.
73
4.5.1 Case I: no mismatch in plant model, subject to IMA(1,1) distur-
bance
In this case, we assume there is no mismatch of the plant model. The
additive output disturbance is generated by
d(t) =
1 0:3q
1
1q
1
e(t) (4.28)
where e(t) N(0; 1). Due to the integrator, the disturbance sequence is
non-stationary.
0 5 10 15 20 25 30
0
0.2
0.4
0.6
0.8
1
step
Figure 4.1: Identified system Markov parameters of the disturbance model
in case I.
The integrating disturbance model A
d
= I with K
d
is sufficient for
(4.28). Fig. 4.1 shows the estimated system Markov parametersH
i
. It can
be seen that since there is no mismatch,H
i
(i = 1;:::; 30) is flat. The values
is quite close to the theoretical oneH
i
=C
d
K
d
= 0:7.
The output variance of traditional DMC is 1.0925, while the output
variance of DMC with disturbance model is 0.9989. Since both controllers
74
shares the same model, the difference in control performance is caused by
K
d
.
4.5.2 Case II: no plant model mismatch, subject to ARIMA(1,1,1) dis-
turbance
Now consider the disturbance to be more complicated:
d(t) =
1 0:3q
1
1q
1
1
1 + 0:5q
1
e(t) (4.29)
which has an additional first-order filter compared to (4.28). In this case,K
p
will no longer be zero, andH
i
will be varying in accordance with (4.26).
0 5 10 15 20 25 30
0
0.2
0.4
0.6
0.8
1
step
Figure 4.2: Identified system Markov parameters of the disturbance model
in case II.
Estimated system Markov parameters of disturbance model is
demonstrated by Fig. 4.2. One can observe that unlike Case I, the system
Markov parameters do not reach its final value in one step. This is because
first term on the right hand side of (4.26) is nonzero.
75
The output variance of traditional DMC is 1.8684, whereas the out-
put variance of DMC with disturbance model is 1.4965. If compared to the
control results in Case I, DMC with disturbance model has a larger output
variance, which is the same as variance of prediction errors when minimum
variance control is applied. This shows that reduced order modeling of dis-
turbances tends to yield larger prediction errors.
4.5.3 Case III: with plant model mismatch, subject to ARIMA(1,1,1)
disturbance
We further add plant model mismatch to the simulation. Let the
model (4.27) be multiplied by 1:2 (1 0:4q
1
)=(1 0:5q
1
). It is then
converted to step disturbance model and fed into DMC. Then, the system
matrices of plant model can be described by
A =
0
@
0:6 0:3 0
0 0:2 0:4
0 0 0:5
1
A
; B =
0
@
0
1:2
0:3
1
A
C =

1 1 0

; D = 0:
whose order is higher than that of the plant. There is a pole at q = 0:5 in
addition to the polesq = 0:6 andq = 0:2 of the plant. The disturbance to be
added is still (4.29).
In this work, our purpose is to find the best non-parametric distur-
bance model that accommodates the existing DMC model. To achieve best
control performance, an alternative approach is to substitute both plant and
76
0 5 10 15 20 25 30
0
0.2
0.4
0.6
0.8
1
step
Figure 4.3: Identified system Markov parameters of the disturbance model
in case III a).
disturbance models from close-loop identification results. We make a com-
parison between two approaches in this case.
a) Disturbance model identified by the proposed approach is first shown in
Fig. 4.3. It is observed that due to the model mismatchH
i
oscillates less.
This implies that the mismatch has been modeled into the disturbance
part. The output variance of traditional DMC is 1.5888, compared to
1.2482 of DMC with disturbance model.
b) Fig. 4.4 provides the Markov parameters of the identified plant and dis-
turbance models. The output variance of the process controlled by the
MPC using identified process and disturbance models is 1.3950, which is
greater than the proposed method that only updates disturbance model.
After carefully checking the identified model, it is found that the MAT-
LAB n4sid method, which is capable of performing closed-loop subspace
77
0 5 10 15 20 25 30
0
1
2
step
0 5 10 15 20 25 30
0
0.5
1
step
Figure 4.4: Identified system Markov parameters of the plant and distur-
bance model in case III b). The upper one is Markov parameters of the
plant model, and the lower one is the Markov parameters of the disturbance
model.
Table 4.1: Control Performance of Different DMC Schemes
Output variance
Traditional DMC DMC w/ dist. mdl. DMC w/ proc.
& dist. mdl.
Case I 1.0925 0.9989 N/A
Case II 1.8684 1.4965 N/A
Case III 1.5888 1.2482 1.3950
identification, has problem in determining the order of augmented sys-
tem. A second-order is estimated, as one can observe in Fig. 4.4, while
the true augmented system is 5-th order. The reduced-order model iden-
tified by traditional SID methods is unable to capture some dynamics in
the true model. This problem can be even worse, when the process has
considerably higher order than the disturbance model does.
78
Table 4.1 summarizes the control performance of both controllers in
all three cases. It can be concluded from the comparison that DMC with
non-parametric disturbance model is able to compensate the performance
degradation caused by plant model mismatch or high order disturbances.
4.6 Summary
In this chapter, we proposed a subspace method for estimating an
integrated moving average disturbance model for MPC from closed-loop
data. In contrast to the traditional subspace identification method, the
proposed approach does not focus on system matrices A, B, C, which
are known in the MPC model. Instead, we try to utilize closed loop data
and the model information in MPC to obtain or update the disturbance
model. Estimation of system matrices and the Kalman gain, which is
usually an intermediate step, is then avoided. We further incorporate
this non-parametric disturbance model with DMC, which typically uses
a (non-parametric) step response plant model. The simulation results
demonstrate that DMC with the estimated disturbance model is able to
compensate model mismatch and outperforms traditional DMC that uses a
step-like disturbance model.
79
Chapter 5
Performance Monitoring of Model-Based Con-
trollers via Model Residual Assessment
5.1 Introduction
Control performance monitoring (CPM) is concerned with assess-
ing the performance of an online controller against some benchmark
estimated from closed-loop data. Originally developed by
˚
Astr¨ om [2] to
achieve minimum variance control of the controlled output of an single-
input-single-output (SISO) process, the minimum variance benchmark is
developed by Harris [27] who shows that this best achievable benchmark
can be estimated from normal closed-loop data. This control performance
monitoring method requires the knowledge of time delay. Consider a
SISO process under feedback control as shown in Fig. 5.1, where G

(q)
and H

(q) are the true but unknown process and disturbance models,
G
c
(q) is a linear time-invariant (LTI) controller, and q is the backward
shift operator. e

(k) is white noise with zero mean. The leading impulse
response coefficients of H

(q) up to the time delay of G

(q) are feedback
invariant and are used to calculate the minimum variance benchmark. For
multi-input-multi-output (MIMO) systems, the time delay information
80
of the process is described by the interactor matrix, which requires some
leading impulse response coefficients ofG

(q) as well [34]. There are many
forms of interactors, such as the unitary [34] and triangular interactors
[78], leading to different definitions of multivariate minimum variance
controllers. Interactors other than simple or diagonal structures are usually
difficult to obtain, as they require a number of leading impulse response
coefficients of the process model. To alleviate the need for more complex
interactors, Yu and Qin propose a combined left/right diagonal interactor
to ease the computation of a general one [90]. However, knowledge about
process time delays is still required. Other methods try to utilize less a
priori process model knowledge. For example, Huang et al. [31] use the
order of interactor matrix (OIM) instead of the full interactor to compute
a benchmark. Unfortunately, the use of reduced knowledge of process
models usually leads to an approximation of the benchmark, which may
not be achievable.
There has been increasing interest in assessing the control perfor-
mance of MPC from both industry and academia. The design-case method
[65] and the expectation-case method [92] take the advantage of predictions
made in MPC; they compare the actual key performance index (KPI) with
the designed KPI based on predicted outputs. The KPI refers to the control
objective function chosen at the MPC design stage. Schaefer and Cinar [76]
81
G
c
q G
o
q
y k
H
o
q
e
o
k
u k
r k
Figure 5.1: Schematic diagram of a closed-loop controlled process.
improve the method by using the ratio of the historical and actual perfor-
mance as performance metrics and using the ratio of designed and achieved
performance for diagnosis. These approaches, however, have drawbacks
in handling the case with stochastic disturbances. Furthermore, the use
of control objective as a benchmark assumes the adequacy of the existing
process model, which by itself can be a major cause of poor control perfor-
mance. Recently, Harrison and Qin [30] develop a minimum variance per-
formance map that considers the impact of constraints on the MPC perfor-
mance, showing that the minimum-variance-based method is inadequate
for MPC performance assessment. This work, rather than being a practical
assessment method, elucidates the complexity of the MPC performance.
The challenge of MPC performance monitoring is multi-facet be-
cause performance degradations can come from inappropriate control
parameters such as horizon lengths and weights in the objective functions,
82
poor model quality in either the input-output plant channel or the distur-
bance channel, inappropriate constraint setup, or inconsistency between
the upper level optimization and the dynamic MPC. Among these factors,
plant and disturbance models are the most important part of MPC that can
greatly affect the control performance. Predictions from a poor model can
result in computed inputs to be far from optimal control moves. This can
result in significant economic loss, or even worse, instability. Practitioners
usually spend up to 80% of the overall MPC design effort to obtain an
adequate model for MPC. When the control performance degrades, engi-
neers usually try to perform model maintenance by conducting a system
identification exercise due to lack of a model assessment tool, leading to
high cost of MPC maintenance. Julien et al. [38] present a performance
assessment method focusing on the model for SISO MPC. However, in
their approach, re-identification of both plant and disturbance models is
required.
It is desirable to have a performance monitoring method that can
derive the most information from closed loop operation data and does not
require a system identification exercise or other equivalent model informa-
tion. Should a system identification exercise be required, one might as well
redesign the controller instead of assessing its performance. Although CPM
methods based on KPI can be a useful means for assessing the MPC control
83
performance, it is a lumped metric of many aspects of MPC. For example,
in addition to model deficiency, tuning problems, excessive moves of set-
points and changes in constraints can also contribute to the KPI degrada-
tion. Therefore, it is desirable to be able to assess the model quality under
closed loop using routine operation data only. If the adequacy of the model
can be determined in a non-intrusive way independent of the control pa-
rameters, one can easily isolate performance degradation due to model de-
ficiency or inappropriate control design.
In this chapter, we propose a new closed-loop model assessment ap-
proach to assessing model deficiency from routine closed-loop data. Simi-
lar to the minimum variance benchmark for CPM proposed by Harris [27],
we propose a new model quality index (MQI) that is a minimum variance
benchmark for the model residuals obtainable from closed-loop data only.
To obtain the plant disturbance innovations from closed-loop data, we start
from the feedback invariant principle to show that the disturbance innova-
tions at current instance cannot be affected by the controller even if it is a
nonlinear time-varying (NTV) controller. Then it is shown that the distur-
bance innovations sequence can be estimated from closed loop data by an
orthogonal projection of the current output onto the space spanned by past
outputs, inputs or setpoints. With the disturbance innovations as the bench-
mark, a model quality index is developed by using the ratio of a quadratic
84
form of model residuals and that of the estimated disturbance innovations.
5.2 Feedback Invariance of Disturbance Innovations
5.2.1 Motivating SISO processes
We first start from a SISO process under feedback control as shown
in Fig. 5.1. The true but unknown process model is described as follows,
y(k) =G

(q)u(k) +H

(q)e

(k) (5.1)
whereG

(q) andH

(q) are the true but unknown process and disturbance
models,e

(k) is the process innovations, andq is the backward shift opera-
tor. Due to discretization sampling there should be at least one time delay
inG

(q), which can be expressed as
G

(q) =q
1
~
G

(q)
where
~
G

(q) is causal. SinceH

(q) is monic without loss of generality, it is
denoted as
H

(q) = 1 +q
1
~
H

(q)
where
~
H

(q) is causal. Assuming the setpointr(k) does not change, it fol-
lows from Fig. 5.1 that
y(k) =
H

1 +G

G
c
e

(k)
=
"
1 +q
1
~
H

~
G

G
c
1 +G

G
c
#
e

(k)
= e

(k) +L(q)e

(k 1) (5.2)
85
where the argumentq is omitted for brevity. SinceL(q) =
~
H

~
G

Gc
1+G

Gc
is causal,
e

(k) is affected by the controller G
c
(q), i.e., e

(k) is feedback invariant.
Moreover, it can be concluded from (5.2) that the closed-loop output y(k)
can be exactly modeled by an autoregressive moving average (ARMA)
process with the open-loop disturbance innovations as the driving noise,
provided that the closed-loop system is stable . This result implies that the
disturbance innovationse

(k) can be estimated from closed loop data.
5.2.2 General MIMO processes
Consider a general MIMO process still represented by (5.1), where
G

(q) and H

(q) are the plant and disturbance transfer function matrices
with appropriate dimensions, respectively. The disturbance transfer matrix
H

(q) = I + q
1
~
H

(q) is monic, where
~
H

(q) is causal. e

(k) is the dis-
turbance innovations with a covariance
e
. Consider the general case with
setpoint changes. From Fig. 5.1, The relationship fromr(k) ande

(k) toy(k)
is
y(k) = (I +G

G
c
)
1
G

G
c
r(k) + (I +G

G
c
)
1
H

e

(k) (5.3)
= (I +G

G
c
)
1
G

G
c
r(k) + (I +G

G
c
)
1
(I +q
1
~
H

)e

(k)
= (I +G

G
c
)
1
G

G
c
r(k) +e

(k) +L(q)e

(k 1) (5.4)
whereL(q) = (I +G

G
c
)
1
(
~
H

~
G

G
c
) is causal.
Since e

(k) is independent of r(k) and L(q) is causal, e

(k) is inde-
86
pendent of the controller and only the past terms ofe

(k) can be affected by
feedback. Therefore, the closed loop outputy(k) is an autoregressive mov-
ing average with exogenous inputs (ARMAX) process driven by the open-
loop disturbance innovationse

(k), provided that the closed-loop system is
stable. This feedback-invariant result implies that the disturbance innova-
tionse

(k) can be estimated from closed loop data.
5.3 Estimating Disturbance Innovations from Routine
Closed-Loop Data
To establish a model assessment benchmark, it is desirable to use
routine data only without any prior model knowledge or external perturba-
tions. In the previous section, the open-loop disturbance innovationse

(k)
are shown to be the same as closed-loop innovations, i.e., the disturbance
innovations are feedback invariant. Based on this finding, in this section,
we propose a new method to estimatee

(k) from routine closed-loop data.
The following assumption is made first.
Assumption 5.1. The closed-loop controlled outputy(k) is stationary.
In practice, Assumption 5.1 is generally satisfied if the closed-loop
system is stable and the setpoint change is stationary, even though the
disturbance H

(q)e

(k) is necessarily non-stationary. As usual we assume
87
H

(q) is integrating to be non-stationary. Denoting
T
G
= (I +G

G
c
)
1
G

G
c
T
H
= (I +G

G
c
)
1
H

as the closed-loop transfer matrices, where T
H
is integrating and monic
sinceH

(q) is integrating and monic, (5.3) becomes
y(k) =T
G
r(k) +T
H
e

(k) (5.5)
It is straightforward to derive the following one-step ahead prediction from
(5.5) using, for example, results in the test by Ljung [48],
^ y(kjk 1) = (IT
1
H
)y(k) +T
1
H
T
G
r(k) (5.6)
y(k) = ^ y(kjk 1) +e

(k) (5.7)
whereT
1
H
= H
1
(I +G

G
c
) exists and is monic becauseH

is invertible
and monic. As a consequence, (IT
1
H
) is stable and strictly causal. Further,
T
G
is strictly causal becauseG

is strictly causal. Therefore, denoting
IT
1
H
=
1
X
i=1
H
i
q
i
T
1
H
T
G
=
1
X
i=1
G
i
q
i
(5.6) and (5.7) can be combined and expanded into
y(k) =
1
X
i=1
H
i
y(ki) +
1
X
i=1
G
i
r(ki) +e

(k)

M
X
i=1
H
i
y(ki) +
N
X
i=1
G
i
r(ki) +e

(k) (5.8)
88
for sufficiently largeM andN, which is a stable high order ARX (HOARX)
model.
It is now clear that e

(k) can either be estimated from the residuals
of (5.5), an ARMAX model, or those of (5.8), a HOARX model. Taking the
route of the HOARX model and defining
y
p
(k) =

y(k) y(k 1) y(kp)

e

p
(k) =

e

(k) e

(k 1) e

(kp)

r
p
(k) =

r(k) r(k 1) r(kp)

Y
M
(k 1) =
2
6
6
6
4
y
p
(k 1)
y
p
(k 2)
.
.
.
y
p
(kM)
3
7
7
7
5
;R
N
(k 1) =
2
6
6
6
4
r
p
(k 1)
r
p
(k 2)
.
.
.
r
p
(kN)
3
7
7
7
5

Z
p
(k) =

Y
M
(k 1)
R
N
(k 1)

wherep is the window size of the data, we can obtain an extended version
of (5.8):
y
p
(k) =

L
p

Z
p
(k) +e

p
(k) (5.9)
where

L
p
= [H
1
H
2
H
M
G
1
G
2
G
N
]. The following theorem is in order
to obtain the disturbance innovationse

(k) from closed-loop data.
Theorem 5.1. Define

?

Zp(k)
=I

Z
p
(k)
T

Z
p
(k)

Z
p
(k)
T

1

Z
p
(k) (5.10)
89
as the orthogonal projection to the complement of the row space of

Z
p
(k). For a lin-
ear process controlled by an LTI controller, the sequence of disturbance innovations
e

(k) is obtained by performing the orthogonal projection on (5.9),
e

p
(k) =y
p
(k)
?

Zp(k)
asp!1: (5.11)
Proof. Using

Z
p
(k)
?

Zp(k)
= 0, post-multiplying (5.9) by
1
p

Zp(k)
leads to
1
p
y
p
(k)
?

Zp(k)
=
1
p

L
p

Z
p
(k)
?

Zp(k)
+
1
p
e

p
(k)
?

Zp(k)
=
1
p
e

p
(k)
1
p
e

p
(k)

Z
p
(k)
T

Z
p
(k)

Z
p
(k)
T

1

Z
p
(k) (5.12)
Since e

(k) is uncorrelated to past data and setpoint changes, i.e.,
E

e

(k)y(ki)
T

= 0 andE

e

(k)r(ki)
T

= 0,8i 1, we have
1
p
e

p
(k)

Z
p
(k)
T
! 0 asp!1:
Therefore, the second term on the right hand side of (5.12) approaches zero
whenp!1, and hence
e

p
(k) =y
p
(k)
?

Zp(k)
asp!1:
For nonlinear or time-varying feedback controllers, the HOARX
model (5.8) offers a linear approximation. Sun et al. [79] present an
approach to dealing with nonlinear time-varying feedback controllers.
90
While Theorem 5.1 gives the least squares estimation of the inno-
vations sequence, it is not an efficient implementation as our objective is
only to estimate the residuals. In addition, since the closed-loop data are
not necessarily well excited, the matrix inversion of

Z
p
(k)

Z
p
(k)
T
can be ill-
conditioned. To implement the projection in Theorem 5.1 efficiently, we
perform the following QR decomposition,

Z
p
(k)
y
p
(k)

=

R
11
R
21
R
22

Q
1
Q
2

(5.13)
Using these relations and the fact that the rows ofQ
1
andQ
2
are orthogonal,
y
p
(k)

Z
p
(k)
T
= R
21
R
T
11

Z
p
(k)

Z
p
(k)
T
= R
11
R
T
11
The estimated innovations sequence in (5.11) can be derived as follows.
e

p
(k) = y
p
(k)

I

Z
p
(k)
T
[

Z
p
(k)

Z
p
(k)
T
]
1

Z
p
(k)

= y
p
(k)R
21
R
T
11
[R
11
R
T
11
]
1

Z
p
(k)
= y
p
(k)R
21
R
+
11

Z
p
(k)
Note that R
21
and R
11
are required only to calculate the innovations se-
quence, whose dimensions are independent of the number of the data sam-
ples. There is no need to calculate or store theQ matrices.
91
5.4 Estimating Disturbance Innovations from Mismatched
Model Residuals
Assume that the plant and disturbance model used in the predictive
controller has the following form,
y(k) =G(q)u(k) +H(q)e(k) (5.14)
which is different from the process (5.1) butH(q) is still monic. Given past
process data y
p
(k 1) and u
p
(k 1), the one-step-ahead prediction and
prediction error are
^ y(kjk 1) = (IH
1
)y(k) +H
1
Gu(k)
e(k) = H
1
[y(k)Gu(k)] (5.15)
= H
1
[(G

G)u(k) +H

e

(k)]
The model residual e(k) can be calculated by (5.15). It is generally differ-
ent from the disturbance innovationse

(k) because of the model mismatch.
If there is no model mismatch for both plant model G(q) and disturbance
modelH(q), the predictor (5.15) is a minimum variance predictor [41], and
e(k) = e

(k). The use of e

(k) to assess e(k) is the proposal of this work
and is given in the next section. In the general case of mismatched models,
the relation betweene

(k) ande(k) is complex due to feedback, bute

(k) is
indeed the innovations ofe(k) as shown in the following theorem.
92
Theorem 5.2. Consider a MIMO process regulated by an LTI controller. Denoting
T
G
= (I +GG
c
)
1
GG
c
andT
H
= (I +GG
c
)
1
H, we have
e(k) =T
1
H
[(T
G
T
G
)r(k) +T
H
e

(k)] (5.16)
which is an ARMAX model withe

(k) as the innovations sequence fore(k).
Proof. From (5.15),
e(k) = H
1
[y(k)GG
c
(r(k)y(k))]
= H
1
[(I +GG
c
)y(k)GG
c
r(k)]
= H
1
(I +GG
c
)[y(k) (I +GG
c
)
1
GG
c
r(k)]
= T
1
H
[y(k)T
G
r(k)]
Substituting (5.5) into the above relation leads to (5.16). Since both T
H

and T
H
are monic, the relation between e

(k) and e(k) in (5.16) is monic
and therefore, (5.16) is an ARMAX model withe

(k) as the innovations se-
quence.
Theorem 5.2 shows that the disturbance innovationse

(k) can be es-
timated from the residuals e(k) from the imperfect model by projecting it
into two orthogonal parts:
e
p
(k) =e

p
(k) + ^ e
p
(k) (5.17)
93
( ) y k
0
ˆ ( ) y k
ˆ( ) y k
( ) e k
ˆ( ) e k
0
( ) e k
Figure 5.2: Relationship of measurements, predictions, simulation errors
and their orthogonal projections. The plane represents the space spanned
byZ
p
(k) or

Z
p
(k).
where
^ e
p
(k) = e
p
(k)
Zp(k)
(5.18)
e

p
(k) = e
p
(k) ^ e
p
(k) (5.19)
and
Zp(k)
=

Z
p
(k)
T
[

Z
p
(k)

Z
p
(k)
T
]
1

Z
p
(k). Note that an efficient QR imple-
mentation of (5.19) is very similar to that for (5.11), simply replacingy
p
(k)
withe
p
(k).
The orthogonal projections to estimate e

(k) from e(k) or y(k) are
depicted in Fig. 5.2. Sincey
p
(k) =L
p
Z
p
(k) +e(k),e

(k) obtained from (5.19)
is the same ase

(k) from (5.11), and hence estimation based on eithere(k)
ory(k) results in the samee

(k) sequence.
5.5 Data-Driven Model Quality Assessment in MPC
5.5.1 MPC objective and key performance index
Consider the process (5.1). A linear MPC can be designed to calculate
the control moves by optimizing the following quadratic objective function
94
min
u(k);:::;u(k+Nu1)
^
J(k) =
P
X
j=1
k(^ y(k +jjk)r(k +j))k
2
Q
+ku(k +j)k
2
S
(5.20)
while fulfilling the constraints
u
min
u(k)u
max
y
min
y(k)y
max
Here, Q and S denote the weight matrices for ^ y(k +jjk) and u(k +j),
respectively, andPN
u
is the prediction horizon.
It is well known thatQ andS are tuning parameters that can signif-
icantly affect the control performance and robustness to mode errors. As-
suming these matrices are tuned properly at the design stage, the control
performance still depends on the model quality, the prediction horizonP ,
and the control horizonN
u
. It is proposed to useQ andS to define a key
performance index (KPI) to actual assess the control performance as [66, 76],
J(k) =
k+P
X
j=k+1
ky(j)r(j)k
2
Q
+ku

(j)k
2
S
(5.21)
which is different from the designed KPI foru

(j) is the closed-loop optimal
control moves, while u(j) in the designed KPI is the closed-loop optimal
control moves. By the relationy(k) = ^ y(kjk 1) +e(k), the actual KPI (5.21)
can be rewritten as
J(k) =
k+P
X
j=k+1
k^ y(jjj 1)r(k +j) +e(j)k
2
Q
+ku

(j)k
2
S
(5.22)
95
If the model is perfect so that the prediction error is white noise, we have
J

(k) =
k+P
X
j=k+1
k^ y(jjj 1)r(k +j) +e

(j)k
2
Q
+ku

(j)k
2
S
whereu

(j) is the closed-loop optimal control move with a perfect model.
It is clear now that the model errore(j) in(5.22) affects to the control perfor-
mance index, which leads to the definition of a model quality index given
in the next subsection.
5.5.2 Model residual assessment
Previous work (e.g. Harrison and Qin [29]) points out that the model
quality is directly related to the model residual. A model residual sequence
that resembles white noise usually indicates a good model, while a model
residual sequence far from white noise indicates a poor model. Sincee

(k)
is the innovations ofe(k) and can be estimated from closed-loop data, we
propose the following model quality index (MQI) with respect to the MPC
control objective,
=
N
X
k=1
e

(k)
T
Qe

(k)
N
X
k=1
e(k)
T
Qe(k)
(5.23)
where Q is the weight matrix for outputs that are properly chosen at the
design stage. The MQI is in (0; 1]. An approaching 1 means thate(k) is
close to e

(k), thus indicating a near perfect model. In contrast, a small
value implies a poor process or disturbance model.
96
The advantage of using MQI for model quality assessment is that the
model errors are weighted by the same matrix that is used in the control
performance index, or KPI. The KPI is not a reliable benchmark unless the
model is perfect or accurate, while the MQI can always be assessed from
closed loop data without any model knowledge.
5.6 Case Studies
5.6.1 Wood-Berry distillation column
The proposed model quality assessment method is first demon-
strated by the Wood-Berry distillation column model [87] is assumed. The
process transfer matrix is described by
G(s) =
2
6
4
12:8e
s
16:7s + 1
18:9e
3s
21:0s + 1
6:6e
7s
10:9s + 1
19:4e
3s
14:4s + 1
3
7
5
:
The input firstu
1
is reflux rate in lb/min; the second inputu
2
is steam flow
rate in lb/min; the first outputy
1
is composition of top product of column
in mol %; and the second output y
2
is composition of bottom product of
column in mol %. With a sampling time of one minute, the discretized
process transfer matrix is
G

(q) =
2
6
4
q
2
0:7665
1 0:9419q
1
q
4
0:9000
1 0:9535q
1
q
8
0:6055
1 0:9123q
3
q
4
1:3472
1 0:9329q
1
3
7
5
:
97
The disturbance model is chosen to be diagonal:
H

(k) =
"
10:5q
1
1q
1
10:7q
1
1q
1
#
(5.24)
wheree

(k)2R
2
is independent white noise with covariance diagf7:6
2
; 0:14
2
g.
The process is controlled by an MPC. The prediction and control
horizons are chosen to be 100 and 10, respectively. The weight matrices
areQ = diagf1; 100g;S = 0. No constraint is enforced for simplicity, thus
implying an LTI controller. The setpoints are selected as y
1
= 90% and
y
2
= 5%.
The MPC model is converted to a state-space model with an added
disturbance model as follows.

^ x
p
(k + 1)
^ x
d
(k + 1)

=

A
p
0
0 A
d

^ x
p
(k)
^ x
d
(k)

+

B
0

u(k) +

0
K

e(k)
y(k) =

C
p
C
d

^ x
p
(k)
^ x
d
(k)

+e(k)
(5.25)
where ^ x
p
(k) andx
d
(k) are the plant and disturbance state estimates, respec-
tively, andA
p
andA
d
are the corresponding state transition matrices. To re-
produce the disturbances used in dynamic matrix control (DMC) we choose
A
d
=I
2
andC
d
=I
2
, which leads to the following disturbance model
H(q) = I
2
+C
d
(qI
2
A
d
)
1
K
=
"
1(1K
1
)q
1
1q
1
1(1K
2
)q
1
1q
1
#
which is an IMA(1,1) model. The step-like disturbance used in dynamic
matrix control (DMC) can be realized by choosing K
1
= 1 and K
2
= 1,
98
which assumes the future disturbance to be the same as the disturbance
measured at current time instant.
We first simulate the base case in which there is no mismatch in either
the plant or the disturbance models. The calculated MQI and KPI are 0.9886
and 64.64, respectively. One can conclude that the model is good enough as
the MQI close to 1. Fig. 5.3 shows more details of the result. The upper plot
is the autocorrelation coefficients of the disturbance source generated in the
simulation. In the lower plot, the innovations sequence e

(k) obtained by
the orthogonal projection is closer to white noise. It can be seen that the
autocorrelation function in the upper plot is not a strictly white noise due
to finite data length.
2 4 6 8 10 12 14 16 18 20
−0.05
0
0.05
Lag
Sample Autocorrelation
2 4 6 8 10 12 14 16 18 20
−0.05
0
0.05
Lag
Sample Autocorrelation
Figure 5.3: Autocorrelation of the first element ofe

(k), lag 1. The upper
one is based on the disturbance source sequence, and the lower one is based
on the sequence calculated from closed-loop data.
Table 5.1 summarizes the model assessment results and the cor-
99
Table 5.1: Model Residual Assessment with Perfect Plant Model but Mis-
matched Disturbance Models
MQI KPI
CV1:K
1
= 1:0, CV2:K
2
= 1:0 0.7881 81.20
CV1:K
1
= 1:0, CV2:K
2
= 0:3 0.8126 77.37
CV1:K
1
= 0:5, CV2:K
2
= 1:0 0.9333 68.47
responding KPI’s when there is no plant model mismatch but there are
disturbance model mismatches. Several disturbance models with different
Kalman gains are applied. Compared to the first row, which is based on the
step disturbance for both CVs, improvements in MQI and MPC KPI can be
observed in the second and third row where a better disturbance model
is applied. If further comparing the results in the second and third row to
the results of Case I, it can be seen that the MQI and KPI of Case I are still
superior as a result of nearly perfect model. Therefore, all the observation
is consistent with the a priori knowledge about the model deficiency.
Table 5.2 gives the results of both plant and disturbance model mis-
matches. The plant model used in the MPC is:
G(q) =
2
6
4
0:8q
2
0:7665
1 0:9419q
1
0:6q
4
0:9000
1 0:9535q
1
q
8
0:6055
1 0:9123q
3
1:2q
4
1:3472
1 0:9329q
1
3
7
5
The disturbance models are IMA(1,1) with various Kalman gains as shown
in Table 5.2. It can be seen that higher MQI corresponds to lower KPI. There-
fore, the MQI can be an effective metrics of the model quality and for MPC
100
performance.
Table 5.2: Model Residual Assessment with Mismatched Plant and Distur-
bance Models
MQI KPI
CV1:K
1
= 1:0, CV2:K
2
= 1:0 0.0250 3205.3
CV1:K
1
= 0:3, CV2:K
2
= 0:5 0.1620 1150.9
CV1:K
1
= 0:05, CV2:K
2
= 0:6 0.3330 385.5
FI
FI
FI
FI
FI
FI
FI
FI
FI
TI
TI
TI
TI
PI
LI
LI
JI
SC
C
A
D
E
CWR
CWR
CWS
CWS
Stm
Cond
Purge
Product
XA
XB
XC
XE
XF
XA
XB
XC
XD
XE
XF
XG
XH
XD
XE
XF
XG
XH
Reactor
Condenser
Stripper
Compressor
Vap/Liq
Separator
LI
XD
FI
TI
PI
PI
1
2
3
4
5
6
7
8
12
13
10
11
9
Figure 5.4: Schematic diagram of the Tennessee Eastman process [1].
5.6.2 Tennessee Eastman challenge problem
The Tennessee Eastman process (TEP) is a benchmark problem de-
veloped by the Eastman Chemical company [17]. It is widely used in pro-
cess systems engineering research to assess process control and monitoring
methods. A schematic diagram is shown in Fig. 5.4. The process consists of
eight components: two main products, four reactants, one byproduct and
101
one inert. It contains five units: a reactor, a condenser, a compressor, a sepa-
rator and a stripper. Two main products as well as the byproduct and inert
are formed in the reactor, which is fed by four gaseous reactants. The pro-
cess contains 12 manipulated variables and 41 measured variables, among
which 19 measured variables are quality variables that are sampled every 3
minutes and the rest 22 variables are measured continuously. There are six
operating modes of the process; Mode 3 is chosen in this work.
The process is open-loop unstable. A variety of control schemes has
been studied by previous researchers [73, 53, 50]. It is reported that the
process is difficult to control using a centralized nonlinear MPC [73]. There-
fore, in this work, in order to study the model quality for MPC, we choose
a two-by-two subsystem based on which the MPC is designed and applied.
The controlled variables of the subsystem are the reactor temperature (CV
9) and the product separator temperature (CV 11) since they are difficult to
control due to coupling of the CVs. Accordingly, the manipulated variables
are the reactor cooling water flow (MV 10) and the condenser cooling water
flow (MV 11). The rest of CVs are regulated by the decentralized controllers
described in the work by Ricker [73], most of which are of PID type. The
TEP will be simulated for 72 hours, and data are sampled every 36 seconds
(0.01 hour).
The model of the two-by-two subsystem is first identified. The MAT-
102
0 10 20
−0.1
−0.05
0
0 10 20
0
0.2
0.4
0.6
0.8
0 10 20
−0.1
−0.05
0
0 10 20
0
0.2
0.4
0.6
0.8
time/min
Figure 5.5: Step responses of plant model of MPC1 controlling TEP .
0 10 20
−6
−4
−2
0
0 10 20
−6
−4
−2
0
0 10 20
−0.04
−0.03
−0.02
−0.01
0
0 10 20
−0.04
−0.03
−0.02
−0.01
0
time/min
Figure 5.6: Step responses of plant model of MPC2 controlling TEP .
103
Table 5.3: MQI of the TEP Example
MQI
MPC1 0.012
MPC2 0.213
LAB function n4sid is used to identify the process model from data collected
under closed-loop control with decentralized PID [73]. The two MVs are
added with PRBS signals with a magnitude of 0.1 to improve closed-loop
identifiability. Among all the 7200 data points collected, two sets of data
are selected for identification; each set consists of 1200 samples. One set
of data is under limited excitation, while the other is excited more due to
larger changes in setpoints. The step response of two identified plant mod-
els are shown in Fig. 5.6 and Fig. 5.5, respectively. It can be seen that the
sub-process has strong interactions. These identified models are used for
the design of two MPCs, MPC1 and MPC2.
The MPC is designed using the MPC Toolbox of MATLAB. The pre-
diction and control horizons are 10 and 2, respectively, and the weights
are Q = diagf100; 1g;S = diagf5; 0:25g, since the reactor temperature has
a smaller variance. The MVs are constrained such that they are kept within
the normal range.
The control results are illustrated in Fig. 5.7. It is seen that MPC2
has much better performance for the TEP process than MPC1. Table 5.3
104
0 20 40 60
120.5
121
121.5
122
122.5
MPC1
reactor temp/
o
C
0 20 40 60
120.5
121
121.5
122
122.5
MPC2
reactor temp/
o
C
0 20 40 60
81
82
83
84
85
condensor temp/
o
C
0 20 40 60
81
82
83
84
85
condensor temp/
o
C
time/hr
Figure 5.7: Control results of TEP by MPC1 and MPC2.
presents the model residual assessment results in terms of MQI. As can been
seen, the MQI of MPC2 is greater than the MQI of MPC1, and therefore has
better control performance. The results in this TEP example demonstrate
the usefulness of MQI to indicate the quality of MPC model.
5.6.3 An air separation plant of Praxair
The same process in Chapter 3 is studied here again as an indus-
trial application of MQI. In Fig. 5.8, the process flow is provided schemat-
ically. The air is filtered, compressed and purified first before feeding into
the high pressure column, where liquid nitrogen and oxygen-enriched liq-
uid air is produced. After subcooling in the nitrogen superheaters and pass-
ing through the crude argon condenser, the liquid air is fed into low pres-
sure column for final separation.
105
Figure 5.8: Simplified diagram of an industrial air separation process [84].
The process is so complicated that it involves 22 MV and 28 CV . A
DMC-type commercial MPC is applied to regulate the process. Originally,
the disturbance model configured in MPC, denoted by MPC1, is the de-
fault one, step-like disturbance. After collected enough closed-loop data,
we apply the disturbance modeling technique introduced in Chapter 3 to
identify an IMA(1,1) disturbance model for each CV . The Kalman gain of
disturbance model is listed in Table 5.4. It is expected that the model qual-
ity indices are increased by the improved disturbance settings of the new
MPC (MPC2).
The model quality assessment results are shown in Table 5.4. It is
shown that for most of the outputs the model quality is improved signifi-
106
Table 5.4: Model Quality Index of Results of the Industrial Air Separation
Process
Kalman Gain MQI of MPC1 MQI of MPC2 Improvement
CV3 0.6842 0.1656 0.2027 22.40%
CV9 0.1037 0.1285 0.0899 -30.04%
CV11 0.8981 0.1730 0.4148 139.77%
CV12 0.5410 0.1671 0.1484 -11.19%
CV13 0.5905 0.1745 0.2417 38.51%
CV14 0.0313 0.1701 0.2039 19.87%
CV15 0.7859 0.2281 0.3623 58.83%
CV16 0.5860 0.1592 0.1876 17.84%
CV23 0.7668 0.1883 0.3150 67.29%
CV24 0.4162 0.5259 0.4817 -8.40%
CV25 0.1014 0.5910 0.5720 -3.21%
cantly as expected. Note that the process is nonlinear, which could be partly
responsible for the performance degradation for some CVs. Nevertheless,
the model quality index shows an average of 28% improvement of the MPC
model. In addition, according to the relative importance of the outputs, the
engineers made a better decision on model tuning based on MQI results.
Therefore, the proposed MQI can serve as an effective closed-loop-based
evaluation tool of MPC model.
5.7 Summary
After decades of research and practice in MPC design, the perfor-
mance monitoring of MPC using routine closed-loop data has received in-
creased attention recently. Model quality among other factors plays a key
107
role in the MPC control performance. In this chapter, a new closed-loop
model assessment approach is successfully developed to assess the model
quality from routine closed-loop data. The proposed model quality index
is a minimum variance benchmark for the current model residuals and is
estimated from closed-loop data without any plant perturbations or model
information. This is a new utilization of the feedback invariant principle,
that is, the disturbance innovations at current instance cannot be affected by
the controller. This approach is analogous to the use of minimum variance
control benchmark by Harris for control performance monitoring since it is
achievable at the optimal conditions and is estimated from routine closed-
loop data.
108
Chapter 6
Control Performance Monitoring of LP-MPC Cas-
cade Systems
6.1 Introduction
Research on control performance monitoring has been developing
quickly these years in both academia and industry. Routine operation data
are compared with some benchmark to determine the performance quality
of the controller. Pioneer work by Harris [27] shows that MVC can be used
as a SISO performance benchmark and be estimated from routing operation
data. For MIMO systems, however, it is difficult to extract time delay struc-
ture known as the interactor without additional model knowledge beyond
the input-output time delays, which plays an important role in determining
feedback invariant terms [34]. Huang et al. [32] use unitary interactors to
assess control performance based on the sum of output variances. Recently,
Yu and Qin [90] propose left and right diagonal interactor to take advan-
tage of diagonal delay structure to avoid the need for the first few Markov
parameters.
The MPC objective includes not only an output penalty term but also
input reference penalty and input move suppression term. Performance
109
benchmark based on MVC may not be achievable by MPC. Huang and
Shah [34] compare MPC with linear quadratic regulator (LQR) and propose
the LQG benchmark. To fully utilize model and constraint information,
some model based approaches, such as the design-case method [65] and
the expectation-case method [92], are proposed. This type of methods
evaluate MPC performance by comparing achieved MPC objective costs
and designed/expected costs. However, these methods ignore the impact
of constraints on the control performance. Recently, Harrison and Qin
[30] propose a minimum variance map considering effects of constraints.
Nevertheless, in the MPC practice, CVs have their priorities according to
economics and safety. The ability of closely tracking setpoints does not
necessarily mean more economic and safe operation.
Industrial MPC applications usually have a steady-state optimizer,
mostly via linear programming (LP), above MPC to optimize MPC setpoints
and achieve better economics [68]. Ying and Joseph [89] provide stabil-
ity conditions for the LP-MPC system. Recent research by Nikandrov and
Swartz [58] provides sensitivity analysis on constraints in LP-MPC. They
point out that for MIMO systems disturbance on one CV may affect other
CVs through interactions. The effects can be even amplified, thus cause
large variations on CVs with small disturbances. For CVs at their bounds,
such behavior is undesirable because the probability of violations of safety
110
limits increases. Meanwhile, if variance of other less critical CVs decreases,
the MVC based performance monitoring methods will still regard controller
as improved, failing to identify potential performance degradations of the
more critical CVs. Although we can add weights to unitary interactor, the
CVs far away from the constraints cannot be excluded completely. There-
fore, it calls for an MPC monitoring approach that considers CVs with pri-
ority.
In this chapter, we propose a method to evaluate the minimum vari-
ance based on CV priority. In the MPC formulation, constraints on ma-
nipulated variables (MV) are always ensured, but due to disturbances CVs
may exceed their bounds which are treated as soft constraints. Hence, CV
constraints of both LP and MPC actually back off from their safety limits.
The conditional MVC with group priority in the proposed method enables
one to reduce violations of safety limits, and the performance benchmark
evaluates the potential reduction of violations. Moreover, improved perfor-
mance constraints of LP-MPC may be moved closer to safety limits, which
generally leads to better economics.
6.2 Revisit of MIMO Interactor and Minimum Variance
Benchmark
Let the MIMO process take the following form:
y
t
=G(q)u
t
+N(q)a
t
(6.1)
111
whereG(q) andN(q) are process and disturbance transfer functions which
are proper and rational; y
t
, u
t
and a
t
are process output, input and noise
vectors. It is further assumed that a
t
is Gaussian with zero mean and co-
variance of
a
. N(q) is then realized in a form such that it is monic. For
nn transfer function matrixG(q), there exists a non-unique, non-singular
nn polynomial matrixD such that det(D) =q
r
and
lim
q
1
!0
DG(q) = lim
q
1
!0
~
G(q) =K (6.2)
whereK is a non-zero and finite matrix with full rank, the integerr is de-
fined as the number of infinite zeros ofG(q), and
~
G is the delay-free transfer
function matrix corresponding to G(q). The interactor matrix, D, can be
written as
D =D
0
q
d
+D
1
q
d1
+ +D
v
q
dv
(6.3)
whered denotes the order of interactor matrixD,v is defined as the relative
degree of D [33], and D
i
(i = 0;:::;v) are coefficient matrices. The inter-
actor matrix D can be categorized into three forms. The simplest form is
D = q
d
I, called simple interactor matrix. A more complicated interactor
matrix can beD = diag(q
d
1
;q
d
2
;:::;q
dn
), which is called diagonal interactor
matrix. Otherwise, for some special processes, interactor matrix can have
off-diagonal elements and is called general interactor matrix. The common
forms of general interactor matrix include lower triangular, nilpotent and
112
unitary matrix. The unitary interactor matrix satisfies
D
T
(q
1
)D(q) =I: (6.4)
For unitary interactor, by introducing interactor filtered output [34,
28]:
~ y
t
=q
d
Dy
t
=
~
F
0
a
t
+
~
F
1
a
t1
+ +
~
F
d1
a
t(d1)
+
~
F
d
a
td
+ ;
the sum of the variances of original output can be minimized as
minE(y
T
t
y
t
) =E(~ y
T
t
~ y
t
) = tr[cov(~ y
t
)] (6.5)
where the first equality is due to unitary interactor matrix. The MVC bench-
mark can be evaluated as
=
E(y
T
t
y
t
)
min
E(y
T
t
y
t
)
: (6.6)
For a lower triangular interactor matrix, the minimum variance is
conditional [18]. The minimum variance control law derived from interac-
tor filtered output minimizes:
1. the variance of first output;
2. the variance of second output on the condition that the variance of
first output is minimized;
113
3. the variance of third output on the condition that the variance of first
output is minimized first and then the variance of second output is
minimized; and so on.
6.3 Introduction of Block Lower Triangular Interactor Ma-
trix
The lower triangular interactor matrix introduced in the previous
section can be extended to block lower triangular interactor matrix.
Definition 6.1. For anynn proper and rational transfer function matrix
G(q),D(q) is a block lower triangular interactor matrix if it satisfies
G(q) =D(q)
1
~
G(q) (6.7)
lim
q
1
!0
~
G(q) = lim
q
1
!0
D(q)G(q) =K (6.8)
and
D(q) =
0
B
B
B
@
D
11
(q) 0 0
D
21
(q) D
22
(q) 0
.
.
.
.
.
.
.
.
.
.
.
.
D
m;1
(q) D
m;2
(q) D
m;m
(q)
1
C
C
C
A
(6.9)
where K denotes a full rank finite and non-zero matrix and D
i;i
(i =
1;:::;m) are unitary interactor matrices satisfying (6.4).
As an example, CVs are partitioned to two groups according to LP
results. We aim at minimizing variance of CVs in Group 1 unconditionally
and then minimize the variance of CVs in Group 2 under the condition
114
that CVs in Group 1 have minimum variance. Hence, we focus on the case
m = 2. The results can be easily extended to other cases with m > 2 in a
similar fashion.
Based on LP results,G(q) is partitioned as
G(q) =

G
1
(q)
G
2
(q)

(6.10)
where G
1
(q) and G
2
(q) are n
1
n and n
2
n transfer function matrices,
respectively. Using algorithms discussed by Huang and Shah [34], one can
obtain unitary interactorsD
11
(q) andD
22
(q) forG
1
(q) andG
2
(q). According
to property of unitary interactor, there exists full rank constant matricesK
1
andK
2
such that
K
1
= lim
q!1
D
11
(q)G
1
(q)
K
2
= lim
q!1
D
22
(q)G
1
(q)
whereK
1
2R
n
1
n
andK
2
2R
n
2
n
. If matrix
K =

K
1
K
2

(6.11)
is non-singular, then
D(q) =

D
11
(q) 0
0 D
22
(q)

(6.12)
is the block lower triangular interactor matrix forG(q).
115
Otherwise, there must be non-zero off-diagonal block in D matrix.
SinceK is non-singular, the left null spaceU
T
of matrix (K
T
2
K
T
1
)
T
can be
factorized in reduced row echelon form.
U
T

K
2
K1

= 0 (6.13)
whereU2R
nun
1. whose rows have the form (0 0 1 );
2. the number of leading zeros in each row is strictly greater than num-
ber of leading zeros in the row above it; and
3. the column containing leading one is the only nonzero component in
its column.
Since bothK
1
andK
2
are full row rank, the leading zeros in last row ofU is
less thann
2
. This means that givenU some rows ofK
2
can be represented
by linear combination ofK
1
andK
2
. Let the leading one in thei-th row of
U be in columnn
i
andU = (u
i;j
).
K
2
(n
i
; :) =
n
2
X
j=n
i
+1
u
i;j
K
2
(j; :) +
n
1
+n
2
X
j=n
2
+1
u
i;j
K
1
(j; :); i = 1;:::;n
u
where K
1
(j; :) and K
2
(j; :) denote j-th row of K
1
and K
2
, respectively, by
borrowing notations in MATLAB. Hence, delays inG
2
(q) can be increased
in a similar fashion as the method constructing lower triangular interactor
116
matrix [26]. LetD
0
2
be (0 D
22
)
T
which is the second row ofD(q). Set
D
1
2
(n
i
; :) =q
d
1
2;n
i
"
D
0
2
(n
i
; :) +
n
2
X
j=n
i
+1
u
i;j
D
2
(j; :)+
n
1
+n
2
X
j=n
2
+1
u
i;j
D
1
(j; :)
#
; i = 1;:::;n
u
(6.14)
where d
1
2;n
i
is a unique integer such that lim
q!1
D
1
2
(n
i
; :)G(q) = K
1
2
(n
i
; :)
is finite and nonzero. If K
1
2
(n
i
; :) linearly independent of rows of K
1
and
K
2
(1 :n
i
1; :), thenn
i
-th row ofD
2
is set to beD
1
2
(n
i
; :). If not, delay ofn
i
-
th row ofD
2
can be increased in the same way. The approach will eventually
terminate since detG(q)6= 0 andD
11
(q);D
22
(q) have at least delay of 1.
Remark 6.1. Although the lower triangular interactor is unique [26], the
block lower triangular interactor is non-unique because of unitary diagonal
blocks. For any unitary interactor matrix D(q) and unitary real matrix ,

D(q) = D(q) is also a unitary interactor matrix.
Example 6.1. Assume
G(q) =
0
B
B
B
B
B
B
@
4q
1
10:4q
1
q
2
10:1q
1
q
2
10:2q
1
q
2
10:3q
1
3q
1
10:1q
1
q
2
10:8q
1
q
2
10:3q
1
2q
2
10:5q
1
q
1
10:5q
1
q
1
10:6q
1
q
1
10:8q
1
q
1
10:7q
1
4q
1
10:7q
1
q
1
10:5q
1
q
1
10:4q
1
q
1
10:6q
1
1
C
C
C
C
C
C
A
=
0
@
G
1
(q)
G
2
(q)
1
A
117
in which the first two CVs are of primary interest. It can be easily verified
that the interactors forG
1
(q) andG
2
(q) are
D
11
(q) =

0:8q 0:6q
0:6q
2
0:8q
2

; (6.15)
D
22
(q) =

q 0
0 q

(6.16)
respectively. The corresponding
K = lim
q!1

D
11
(q) 0
0 D
22
(q)

G
1
(q)
G
2
(q)

=
0
B
B
@
5 0 0 0
0:72 0:2 0:2 1
1 1 1 1
4 1 1 1
1
C
C
A
=

K
1
K
2

Then, the left null space of (K
2
K
1
)
T
in reduced row echelon form is found
to be
U
T
=

1 1 1 0

which indicates thatD
2
(1; :) +D
1
(1; :)D
2
(2; :) = 0. Therefore, we need to
increase the delay for this row.
D
1
2
(1; :) =q
d
1
2;1

D
0
2
(1; :) +D
1
(1; :)D
2
(2; :)

=q
d
1
2;1
[(0 0 q 0) + (0:8q 0:6q 0 0)
(0 0 0 q)]
=q
d
1
2;1
(0:8q 0:6q q q)
118
whered
1
2;1
= 1 leads to a non-singular and finiteK = lim
q!1
DG(q):
K =
0
B
B
@
5 0 0 0
0:72 0:2 0:2 1
0 1:3 1 1:9
4 1 1 1
1
C
C
A
:
Therefore, the block lower triangular interactor matrix forG(q) is
D(q) =

D
1
(q) D
2
(q)

T
=
0
B
B
@
0:8q 0:6q 0 0
0:6q
2
0:8q
2
0 0
0:8q
2
0:6q
2
q
2
q
2
0 0 0 q
1
C
C
A
6.4 LP-based MPC Control Performance Monitoring
6.4.1 LP problem and active CV constraints
Assume that the setpoint of MPC controlling the process described
by (6.1) is obtained from the following LP problem:
min
y; u

T
y +
T
u
s.t. y =G
ss
u +d (6.17)
y
min
y y
max
u
min
u u
max
where y = ( y
1
;:::; y
n
)
T
and u = ( u
1
;:::; u
m
)
T
. Note that the constraints on
y are soft constraints in the MPC optimization problem. For simplicity, the
119
solution to (6.17) is denoted as y; u. The active constraints can be inferred
from sensitivity analysis or the dual solution.
Based on the solution to (6.17), CVs can be separated into two
groups: CVs with active (binding or more important) constraints and CVs
with inactive (non-binding or less important) constraints. The CVs in first
group are of our primary interest, whose variances along normal direction
of active constraints will be minimized. The CVs in second group are of
less importance, but one may still keep them under control.
Remark 6.2. In the case of multiple solutions to the LP problem (6.17), it
is recommended to choose the center of solution subspace as setpoint of
MPC. The reason is that the number of active constraints may be reduced.
For example, let ( y
1
; u
1
) and ( y
2
; u
2
) be two optimal vertices, and y
1
= y
max
1
and y
1
= y
min
1
are the active constraints for ( y
1
; u
1
) and ( y
2
; u
2
) respectively.
Then, ( y; u) = [( y
1
; u
1
) + ( y
2
; u
2
)]=2 is still optimal but none of y
1
constraint
is active.
Remark 6.3. Although the LP problem optimizing setpoint of MPC usually
takes the form (6.17) whose constraints are upper and lower bounds, there
may exists general constraintA
1
y +A
2
u b in LP constraints. Under this
circumstance, active constraints can involve more than one CV . However, it
can be converted to upper and lower bound types of constraints by intro-
ducing variable transformation.
120
6.4.2 Conditional MVC law and performance benchmark under block
lower triangular interactor
With CVs divided to two groups, the conditional MVC law can be
further derived. According to the CV priority,y
t
is partitioned as
y
t
=

y
1;t
y
2;t

wherey
1;t
2R
n
1
is the vector of CVs with high priority andy
2;t
2R
n
2
is the
vector of other CVs. The corresponding block lower triangular interactor
matrix is
D =

D
11
0
D
21
D
22

:
Then, the process (6.1) becomes
y
t
=D
1
~
Gu
t
+Na
t
: (6.18)
Defining interactor filtered output ~ y
t
=q
d
D(q)y
t
, we obtain
~ y
t
=

~ y
1;t
~ y
2;t

=q
d

D
11
0
D
21
D
22

y
1;t
y
2;t

(6.19)
~ y
1;t
=q
d
D
11
y
1;t
(6.20)
~ y
2;t
=q
d
(D
21
y
1;t
+D
22
y
2;t
): (6.21)
Equation (6.18) can be simplified as
~ y
t+d
=
~
Gu
t
+DNa
t
: (6.22)
Since N(q) is a rational transfer function matrix of disturbance dynamics
without delay, DN has some factors with positive power of q. It can be
121
factorized as
DN =F
d
q
d
+F
d1
q
d1
+ +F
1
q
| {z }
F (q)
+R(q):
Equation (6.22) becomes
~ y
t+d
=
~
Gu
t
+Fa
t
+Ra
t
: (6.23)
The MVC law is achieved by
u
t
=
~
G
1
Ra
t
(6.24)
since
~
G is invertible. We have the following lemma for MVC law (6.24)
under block lower triangular interactor.
Lemma 6.1. The control law (6.24) minimizes
1. the sum of variances of firstn
1
CVs; and
2. the sum of variances of lastn
2
CVs when 1) is achieved.
Proof. Consider first n
1
CVs. Multiplying e
T
1
= (I
n
1
0) to both sides of
(6.22) or (6.23), the first block row of (6.22) is obtained:
~ y
1;t+d
=e
T
1
~
Gu
t
+e
T
1
Fa
t
+e
T
1
Ra
t
: (6.25)
Because ~ y
t
is stationary anda
t
is white noise,
tr[cov(~ y
1;t
)] = tr[cov(~ y
1;t+d
)]
= tr[cov(e
T
1
~
Gu
t
+e
T
1
Ra
t
) + cov(e
T
1
Fa
t
)]
tr[cov(e
T
1
Fa
t
)] = tr(e
T
1
F
a
F
T
e
1
) (6.26)
122
On the other hand, substituting control low (6.24) leads to
~ y
1;t+d
=e
T
1
Fa
t
: (6.27)
which means the equality in (6.26) holds. Recalling D
11
is unitary, (6.20)
implies cov(~ y
1;t
) = cov(y
1;t
). Therefore, the sum of variances of firstn
1
CVs
is minimized.
Next, consider the last n
2
CVs. Similarly, by multiplying e
T
2
=
(0 I
n
2
) to both sides of (6.22), we obtain
~ y
2;t+d
=e
T
2
~
Gu
t
+e
T
2
Fa
t
+e
T
2
Ra
t
: (6.28)
Further, substitute (6.21) into (6.28),
D
21
y
1;t
+D
22
y
2;t
=e
T
2
~
Gu
t
+e
T
2
Fa
t
+e
T
2
Ra
t
D
22
y
2;t
=D
21
D
T
11
~ y
1;t
+e
T
2
~
Gu
t
+e
T
2
Fa
t
+e
T
2
Ra
t
(6.29)
where D
T
11
= D
T
11
(q
1
) is the conjugate transpose. Since the control law
(6.24) yields (6.27), (6.29) becomes
D
22
y
2;t
=D
21
D
T
11
e
T
1
Fa
t
+e
T
2
~
Gu
t
+e
T
2
Fa
t
+e
T
2
Ra
t
=e
T
2
~
Gu
t
+e
T
2
Ra
t
+ (e
T
2
D
21
D
T
11
e
T
1
)Fa
t
: (6.30)
Recall thatD
22
is unitary,
tr[cov(y
2;t
)] = tr[cov(e
T
2
~
Gu
t
+e
T
2
Ra
t
) + cov((e
T
2
D
21
D
T
11
e
T
1
)Fa
t
)]
tr[cov[(e
T
2
D
21
D
T
11
e
T
1
)Fa
t
]] (6.31)
123
where the equality holds if and only if the control law (6.24) is applied.
Therefore, the sum of variances of lastn
2
CVs is minimized by (6.24) when
1) is achieved.
Remark 6.4. Block diagonal interactor is a special case of block lower tri-
angular interactor. In this case, minimum variance of the CV groups are
independent from one another.
Based on the MVC law (6.24), the conditional minimum variance
benchmark is given by
J
MV;1
= tr(e
T
1
F
a
F
T
e
1
) (6.32)
for CVs with higher priority and
J
MV;2
= tr[(e
T
2
D
21
D
T
11
e
T
1
)F
a
F
T
(e
2
e
1
D
11
D
T
21
)] (6.33)
for CVs with lower priority whenJ
1
is minimized.
6.5 Simulation Results
In this section, an example is demonstrated.
Example 6.2. The following process is taken from [90].
G(q) =

q
1
0
q
1
q
2

and
N(q) =

1
1q
1
1
1q
1

124
where
a
= 0:01. Suppose the LP problem is
max
y; u
y
2
s.t. y =G
ss
u +d
0 y
1
; y
2
2
0 u
1
; u
2
5
where d is integrating white noise subject to (6.2) such that it does not af-
fect active set. There are two alternative solutions ( y
1
; y
2
)
T
= (0; 2)
T
and
( y
1
; y
2
) = (0; 2)
T
. The middle point of them (1; 2)
T
is selected as MPC set-
point. The only active constraint is y
2
= 2. Reversing y
1
and y
2
to get y
2
prioritized, we obtain
G
0
(q) =

q
1
q
2
q
1
0

and
N
0
(q) =

1
1q
1
1
1q
1

The unitary interactor ofG
0
(q) is
D
0
u
=

q
2
=
p
2 q
2
=
p
2
q=
p
2 q=
p
2

and the corresponding feedback invariant term is
F
0
u
=

(1 +q
1
)=
p
2 (1 +q
1
)=
p
2
q
1
=
p
2 q
1
=
p
2

:
The variance ofy
2
using MVC based on unitary interactor is
var(y
2
) = var(y
0
1
) =e
T
1
F
0
u

a
F
0
T
u
e
1
= 0:02:
125
After partitionG
0
(q) according to LP results, the (block) lower trian-
gular interactor is
D
0
tri
=

q 0
q
2
q
2

and feedback invariant term becomes
F
0
tri
=

q 0
q
2
+q q
2
q

:
The minimum variance ofy
2
based on block lower triangular interactor is
var(y
2
) = var(y
0
1
) =e
T
1
F
0
tri

a
F
0
T
tri
e
1
= 0:01:
In this case, variance ofy
2
is reduced at the expense of variance ofy
1
, as can
be observed from second row ofF
0
tri
.
Control results of both controllers are illustrated in Fig. 6.1. Sincey
2
upper limit constraint is active, MVC based on block-lower triangular inter-
actor minimizes the variance ofy
2
first and then minimizes the variance of
y
1
. From Fig. 1, we can see that the variance ofy
2
controlled by block-lower
triangular interactor based MVC is the smallest. This verifies the condi-
tional minimum variance property of the proposed method. If the safety
limit ofy
2
is chosen asy
max
2
= 2:3, the violation rate of conditional MVC and
traditional unitary interactor based MVC are 0.2% and 1%, respectively.
Next, performance monitoring will be conducted on an MPC. Since
LP results are known, one may want to emphasize more ony
2
. The weights
126
0 0.5 1 1.5 2
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
2.3
2.4
2.5
y
1
y
2
Figure 6.1: Control results of Example 6.2 in y
1
-y
2
plane. Green markers
are results by MVC using unitary interactor; blue markers are results by
conditional MVC using block lower triangular interactor; and red markers
are the dynamic control results by an MPC.
for CV are chosen to be (1; 4), and the weights for MV are (0:5; 0:5). Predic-
tion and control horizon are 5 and 2 respectively. Results are also shown in
Fig. 6.1. The variance ofy
2
is
J
y
2
= 0:0150
as compared to
J
MV;y
2
= 0:01:
The ratio isJ
MV;y
2
=J
y
2
= 0:67. It is concluded that this MPC has the potential
to improve in terms ofy
2
variance ory
2
upper limit violations.
127
6.6 Summary
For the LP-MPC system, we prioritize CVs based on importance of
their constraints for performance assessment. Since CV constraints are gen-
erally soft constraints, smaller variance along normal direction of these con-
straints results in less violations of the constraint limits. To achieve this,
we introduce block lower triangular interactor to construct the conditional
MVC law, whose output variance is minimized with priority and used as
the performance benchmark. The method is illustrated by simulation re-
sults showing that traditional unitary interactor based method is unable to
distinguish the importance of CVs and that the proposed method is.
128
Chapter 7
Conclusions
This dissertation presents serval new methods focusing on different
aspects of model predictive control, including disturbance modeling, adap-
tation and control performance monitoring. SISO run-to-run control, which
is a simple model-based controller, is studied first. It has been found that
even with an incorrect choice of process gain the disturbance adaptation
scheme will drive the controller to converge to the optimal control law. It
could be concluded that the gain mismatch is able to be compensated by
adaptive run-to-run control.
The method is further extended to MIMO systems. Several distur-
bance models are studied. A practical disturbance model for industrial
MPC applications could be the diagonal IMA(1,1) model. The model
parameter is related to the prediction error filtering constant, which can
be applied with commercial MPC software. Results from simulations
and an industrial process show that the proposed disturbance adaptation
technique can reduce MPC prediction errors considerably.
Another method for disturbance modeling and adaptation is pro-
posed. Some industrial MPC applications are based on the non-parametric
129
model. A new non-parametric disturbance modeling method is developed
under closed-loop conditions. With the HOARX model identified, it con-
verts observer Markov parameters to system Markov parameters incorpo-
rating with information of plant model. This enables the estimated distur-
bance model to compensate plant model mismatch.
It is usually unnecessary to update MPC model frequently. A control
performance assessment approach is developed for the purpose of deter-
mining when adaptation is need. Based on the feedback invariant principle,
the disturbance innovation sequence can be estimated from the closed-loop
data. The model quality of MPC is determined by the ratio of variance of
prediction errors and variance of disturbance innovations. The advantages
of this method are (a) no process knowledge is required, (b) the benchmark
is best achievable, and (c) no exogenous perturbation is needed. Results of
two cases show that the proposed model quality index is an effective tool
measuring MPC model quality.
For the LP-MPC cascade system, in which the MPC setpoints are de-
termined by an LP problem, the CVs are prioritized according to LP solu-
tions. Block lower triangular interactor is introduced to describe the de-
lay structure of systems with CV priorities. A minimum variance based
benchmark is proposed to evaluate the conditional minimum variance. A
simulation is performed to compare the proposed benchmark with other
130
conventional performance benchmarks.
Future directions of research could be further exploration of the
model quality index. The model quality for a specific input/output channel
could be investigated. It could be also an interesting topic to quantify
how much the plant model mismatch can be compensated by adapting
disturbance.
131
Bibliography
[1] C. F. Alcala, “Fault diagnosis with reconstruction-based contributions
for statistical process monitoring,” Ph.D. dissertation, University of
Southern California, Los Angeles, CA, 2011.
[2] K. J.
˚
Astr¨ om, Introduction to Stochastic Control Theory. San Diego, Cal-
ifornia: Academic Press, 1970.
[3] K.
˚
Astr¨ om and B. Wittenmark, “On self tuning regulators,” Automatica,
vol. 9, no. 2, pp. 185 – 199, 1973.
[4] A. Bemporad, “Reducing conservativeness in predictive control of con-
strained systems with disturbances,” in Decision and Control, 1998. Pro-
ceedings of the 37th IEEE Conference on, vol. 2, Dec 1998, pp. 1384 – 1389
vol.2.
[5] A. Bemporad and M. Morari, “Robust model predictive control: A
survey,” in Robustness in identification and control, ser. Lecture Notes
in Control and Information Sciences, A. Garulli and A. Tesi, Eds.
Springer Berlin / Heidelberg, 1999, vol. 245, pp. 207–226.
[6] A. Bemporad and E. Mosca, “Fulfilling hard constraints in uncertain
linear systems by reference managing,” Automatica, vol. 34, no. 4, pp.
451 – 461, 1998.
[7] G. Box and G. Jenkins, Time Series Analysis: Forecasting and Control.
Oakland, California: Holden-Day, 1976.
[8] P . J. Campo and M. Morari, “Robust model predictive control,” in
American Control Conference, June 1987, pp. 1021 – 1026.
[9] E. D. Castillo and A. Hurwitz, “Run-to-run process control: Literature
review and extensions,” Journal of Quality Technology, vol. 29, no. 1, pp.
184–196, 1997.
132
[10] A. Chiuso, “On the relation between cca and predictor-based subspace
identification,” Automatic Control, IEEE Transactions on, vol. 52, no. 10,
pp. 1795 –1812, 2007.
[11] C. Cutler and B. Ramaker, “Dynamic matrix control – a computer con-
trol algorithm,” in AIChE National Meeting, 1979.
[12] D. DeHaan, M. Guay, and V . Adetola, Adaptive Robust MPC: A
Minimally-Conservative Approach, 2nd ed., ser. Lecture Notes in Con-
trol and Information Sciences. Berlin, Germany: Springer, 2009, vol.
384/2009, pp. 55 – 67.
[13] E. Del Castillo and A. Hurwitz, “Run-to-run process control: Literature
review and extensions,” Journal of Quality Technology, vol. 29, pp. 184–
196, 1997.
[14] E. Del Castillo and R. Rajagopal, “A multivariate double EWMA pro-
cess adjustment scheme for drifting processes,” IIE Transactions, vol. 34,
pp. 1055–1068, 2002.
[15] E. Del Castillo and J.-Y. Yeh, “An adaptive run-to-run optimizing con-
troller for linear and nonlinear semiconductor processes,” Semiconduc-
tor Manufacturing, IEEE Transactions on, vol. 11, no. 2, pp. 285 –295, may
1998.
[16] D. Dougherty and D. Cooper, “A practical multiple model adaptive
strategy for single-loop mpc,” Control Engineering Practice, vol. 11,
no. 2, pp. 141 – 159, 2003.
[17] J. Downs and E. Vogel, “A plant-wide industrial process control prob-
lem,” Computers & Chemical Engineering, vol. 17, no. 3, pp. 245–255,
1993.
[18] L. Dugard, G. Goodwin, and X. Xianya, “The role of the interactor ma-
trix in multivariable stochastic adaptive control,” Automatica, vol. 20,
no. 5, pp. 701 – 709, 1984.
[19] B. Francis and W. Wonham, “The internal model principle of control
theory,” Automatica, vol. 12, no. 5, pp. 457 – 465, 1976.
133
[20] H. Fukushima, T.-H. Kim, and T. Sugie, “Adaptive model predictive
control for a class of constrained linear systems based on the compari-
son model,” Automatica, vol. 43, no. 2, pp. 301 – 308, 2007.
[21] C. E. Garcia and M. Morari, “Internal model control. 2. design pro-
cedure for multivariable systems,” Industrial & Engineering Chemistry
Process Design and Development, vol. 24, no. 2, pp. 472 – 484, 1985.
[22] C. E. Garcia and M. Morari, “Internal model control. 3. multivariable
control law computation and tuning guidelines,” Industrial & Engineer-
ing Chemistry Process Design and Development, vol. 24, no. 2, pp. 484 –
494, 1985.
[23] C. E. Garcia and M. Morari, “Internal model control. a unifying review
and some new results,” Industrial & Engineering Chemistry Process De-
sign and Development, vol. 21, no. 2, pp. 308 – 323, 1985.
[24] C. E. Garcia and A. M. Morshedi, “Quadratic programming solution
of dynamic matrix control (QDMC),” Chemical Engineering Communica-
tions, vol. 46, no. 1, pp. 73 – 87, 1986.
[25] C. E. Garcia, D. M. Prett, and M. Morari, “Model predictive control:
Theory and practice–a survey,” Automatica, vol. 25, no. 3, pp. 335–348,
1989.
[26] G. C. Goodwin and K. S. Sin, Adatpive Filtering, Prediction and Control.
Englewood Cliffs, NJ: Prentice-Hall, 1984.
[27] T. J. Harris, “Assessment of control loop performance,” The Canadian
Journal of Chemical Engineering, vol. 67, no. 10, pp. 856–861, 1989.
[28] T. J. Harris, F. Boudreau, and J. F. Macgregor, “Performance assessment
of multivariable feedback controllers,” Automatica, vol. 32, no. 11, pp.
1505–1518, 1996.
[29] C. A. Harrison and S. J. Qin, “Discriminating between disturbance and
process model mismatch in model predictive control,” Journal of Process
Control, vol. 19, no. 10, pp. 1610–1616, 2009.
134
[30] C. A. Harrison and S. J. Qin, “Minimum variance performance map
for constrained model predictive control,” Journal of Process Control,
vol. 19, no. 7, pp. 1199 – 1204, 2009.
[31] B. Huang, S. Ding, and N. Thornhill, “Practical solutions to multivari-
ate feedback control performance assessment problem: reduced a pri-
ori knowledge of interactor matrices,” Journal of Process Control, vol. 15,
pp. 573–583, 2005.
[32] B. Huang, S. L. Shah, and E. K. Kwok, “Good, bad or optimal? per-
formance assessment of multivariable processes,” Automatica, vol. 33,
no. 6, pp. 1175–1183, 1997.
[33] B. Huang, S. X. Ding, and N. Thornhill, “Practical solutions to mul-
tivariate feedback control performance assessment problem: reduced
a priori knowledge of interactor matrices,” Journal of Process Control,
vol. 15, no. 5, pp. 573 – 583, 2005.
[34] B. Huang and S. L. Shah, Performance Assessment of Control Loops: The-
ory and Applications, ser. Advances in Industrial Control. London:
Springer-Verlag, 1999.
[35] M. Jansson, “Subspace identification and arx modelling,” in Proceed-
ings of the 13th IFAC SYSID Symposium, Rotterdam, NL, 2003.
[36] G. Ji, K. Zhang, and Y. Zhu, “A method of mpc model error detection,”
Journal of Process Control, vol. 22, no. 3, pp. 635 – 642, 2012.
[37] J.-N. Juang, M. Phan, L. Horta, and R. W. Longman, “Identification of
observer/kalman filter markov parameters: Theory and experiments,”
Journal of Guidance, Control and Dynamics, vol. 16, no. 2, pp. 320–329,
1993.
[38] R. H. Julien, M. W. Foley, and W. R. Cluett, “Performance assessment
using a model predictive control benchmark,” Journal of Process Control,
vol. 14, no. 4, pp. 441–456, 2004.
135
[39] M. Kano and M. Ogawa, “The state of the art in chemical process con-
trol in japan: Good practice and questionnaire survey,” Journal of Pro-
cess Control, vol. 20, no. 9, pp. 969 – 982, 2010.
[40] S. Karra, R. Shaw, S. C. Patwardhan, and S. Noronha, “Adaptive model
predictive control of multivariable time-varying systems,” Ind. Eng.
Chem. Res., vol. 47, no. 8, pp. 2708–2720, 2008.
[41] T. Katayama, Subspace Methods for System Identification. Springer, 2005.
[42] B.-S. Ko and T. F. Edgar, “Performance assessment of constrained
model predictive control systems,” AIChE Journal, vol. 47, no. 6, pp.
1363–1371, 2001.
[43] D. Kozub, “Controller performance monitoring and diagnosis. indus-
trial perspective,” in 15th Triennial World Congress of IFAC, Barcelona,
Spain, 2002.
[44] P . Kumar and P . Varaiya, Stochastic systems: estimation, identification, and
adaptive control. Englewood Cliffs, New Jersey: Prentice Hall, 1986.
[45] W. Larimore, “Canonical variate analysis in identification, filtering,
and adaptive control,” in Decision and Control, 1990., Proceedings of the
29th IEEE Conference on, 1990, pp. 596 –604 vol.2.
[46] J. H. Lee and C. M. Kiew, “Robust forecasts and run-to-run control for
processes with linear drifts,” Journal of Process Control, vol. 19, no. 4,
pp. 636 – 643, 2009.
[47] J. H. Lee, M. Morari, and C. E. Garcia, “State-space interpretation of
model predictive control,” Automatica, vol. 30, no. 4, pp. 707 – 717, 1994.
[48] L. Ljung, System Identification: Theory for the User. Englewood Cliffs,
New Jersey: Prentice-Hall, Inc., 1999.
[49] L. Ljung and T. S¨ oderstr¨ om, Theory and Practice of Recursive Identifica-
tion. Cambridge, MA: MIT Press, 1983.
136
[50] P . Lyman and C. Georgakis, “Plant-wide control of the tennessee east-
man problem,” Computers & Chemical Engineering, vol. 19, no. 3, pp.
321–331, 1995.
[51] U. Maeder, F. Borrelli, and M. Morari, “Linear offset-free model predic-
tive control,” Automatica, vol. 45, no. 10, pp. 2214 – 2222, 2009.
[52] D. Q. Mayne, J. B. Rawlings, C. V . Rao, and P . O. M. Scokaert, “Con-
strained model predictive control: Stability and optimality,” Automat-
ica, vol. 36, no. 6, pp. 789 – 814, 2000.
[53] T. J. McAvoy and N. Ye, “Base control for the tennessee eastman prob-
lem,” Computers & Chemical Engineering, vol. 18, no. 5, pp. 383–413,
1994.
[54] M. Morari and G. Stephanopoulos, “Studies in the synthesis of control
structures for chemical processes: Part iii: Optimal selection of sec-
ondary measurements within the framework of state estimation in the
presence of persistent unknown disturbances,” AIChE Journal, vol. 26,
no. 3, pp. 247 – 260, 1980.
[55] K. R. Muske and T. A. Badgwell, “Disturbance modeling for offset-free
linear model predictive control,” Journal of Process Control, vol. 12, no. 5,
pp. 617 – 632, 2002.
[56] K. R. Muske and J. B. Rawlings, “Linear model predictive control of
unstable processes,” Journal of Process Control, vol. 3, no. 2, pp. 85 – 96,
1993.
[57] K. R. Muske and J. B. Rawlings, “Model predictive control with linear
models,” AIChE Journal, vol. 39, no. 2, pp. 267 – 287, 1993.
[58] A. Nikandrov and C. L. Swartz, “Sensitivity analysis of lp-mpc cascade
control systems,” Journal of Process Control, vol. 19, no. 1, pp. 16 – 24,
2009.
137
[59] P . V . Overschee and B. D. Moor, “N4sid: Subspace algorithms for the
identification of combined deterministic-stochastic systems,” Automat-
ica, vol. 30, no. 1, pp. 75 – 93, 1994, special issue on statistical signal
processing and control.
[60] P . V . Overschee and B. D. Moor, “A unifying theorem for three sub-
space system identification algorithms,” Automatica, vol. 31, no. 12, pp.
1853 – 1864, 1995.
[61] G. Pannocchia, “Robust disturbance modeling for model predictive
control with application to multivariable ill-conditioned processes,”
Journal of Process Control, vol. 13, no. 8, pp. 693 – 701, 2003.
[62] G. Pannocchia and A. Brambilla, “How to use simplified dynamics in
model predictive control of superfractionators,” Industrial & Engineer-
ing Chemistry Research, vol. 44, no. 8, pp. 2687–2696, 2005.
[63] G. Pannocchia and J. B. Rawlings, “Disturbance models for offset-free
model-predictive control,” AIChE Journal, vol. 49, no. 2, pp. 426 – 437,
2003.
[64] N. Patel and S. Jenkins, “Adaptive optimization of run-to-run con-
trollers: the ewma example,” Semiconductor Manufacturing, IEEE Trans-
actions on, vol. 13, no. 1, pp. 97 –107, feb 2000.
[65] R. S. Patwardhan and S. L. Shah, “Issues in performance diagnostics of
model-based controllers,” Journal of Process Control, vol. 12, no. 3, pp.
413–427, 2002.
[66] R. S. Patwardhan, S. L. Shah, and K. Z. Qi, “Assessing the performance
of model predictive controllers,” The Canadian Journal of Chemical Engi-
neering, vol. 80, no. 5, pp. 954–966, 2002.
[67] S. J. Qin and L. Ljung, “Parallel qr implementation of subspace iden-
tification with parsimonious models,” in Proceedings of the 13th IFAC
SYSID Symposium, Rotterdam, NL, 2003, pp. 1631–1636.
138
[68] S. J. Qin and T. A. Badgwell, “A survey of industrial model predictive
control technology,” Control Engineering Practice, vol. 11, no. 7, pp. 733–
764, 2003.
[69] S. J. Qin, Y. Zhao, Z. Sun, and T. Yuan, “Progressive parametrization
in subspace identification models with finite horizons,” in Decision and
Control (CDC), 49th IEEE Conference on, 2010, pp. 2819 –2824.
[70] M. R. Rajamani, J. B. Rawlings, and S. J. Qin, “Achieving state esti-
mation equivalence for misassigned disturbances in offset-free model
predictive control,” AIChE Journal, vol. 55, no. 2, pp. 396–407, 2009.
[71] C. V . Rao and J. B. Rawlings, “Steady states and constraints in model
predictive control,” AIChE Journal, vol. 45, no. 6, pp. 1266 – 1278, 1999.
[72] J. B. Rawlings and D. Q. Mayne, Model Predictive Control: Theory and
Design. Madison, WI: Nob Hill Publishing, 2009.
[73] N. L. Ricker, “Decentralized control of the tennessee eastman challenge
process,” Journal of Process Control, vol. 6, no. 4, pp. 205–221, 1996.
[74] E. Sachs, A. Hu, and A. Ingolfsson, “Run by run process control: com-
bining spc and feedback control,” Semiconductor Manufacturing, IEEE
Transactions on, vol. 8, no. 1, pp. 26 –43, feb 1995.
[75] M. Safonov and T.-C. Tsao, “The unfalsified control concept and learn-
ing,” Automatic Control, IEEE Transactions on, vol. 42, no. 6, pp. 843 –847,
1997.
[76] J. Schafer and A. Cinar, “Multivariable MPC system performance as-
sessment, monitoring, and diagnosis,” Journal of Process Control, vol. 14,
pp. 113–129, 2004.
[77] F. G. Shinskey, Feedback Controllers for the Process Industries. McGraw-
Hill, 1994.
[78] Z. Sun, S. J. Qin, A. Singhal, and L. Megan, “Control performance mon-
itoring of LP-MPC cascade systems,” in American Control Conference
(ACC), 2011, pp. 4422–4427.
139
[79] Z. Sun, S. J. Qin, A. Singhal, and L. Megan, “Control performance mon-
itoring via model residual assessment,” in American Control Conference
(ACC) 2012, 2012.
[80] M. J. Tenny and J. B. Rawlings, “Closed-loop behavior of nonlinear
model predictive control,” TWMCC technical report, 2002.
[81] S.-T. Tseng, A. Yeh, F. Tsung, and Y.-Y. Chan, “A study of variable
ewma controller,” Semiconductor Manufacturing, IEEE Transactions on,
vol. 16, no. 4, pp. 633 – 643, nov. 2003.
[82] A. Vahidi, A. Stefanopoulou, and H. Peng, “Model predictive control
for starvation prevention in a hybrid fuel cell system,” in American Con-
trol Conference, vol. 1, June-July 2004, pp. 834–839.
[83] M. Verhaegen and P . Dewilde, “Subspace model identification. part i:
the output-error state space model identification class of algorithms,”
Int. J. Control, vol. 56, 1992.
[84] D. R. Vinson, “Air separation control technology,” Computers & Chemi-
cal Engineering, vol. 30, no. 10–122, pp. 1436–1446, 2006.
[85] Y. J. Wang and J. B. Rawlings, “A new robust model predictive control
method i: theory and computation,” Journal of Process Control, vol. 14,
no. 3, pp. 231 – 247, 2004.
[86] W. C. Wong and J. H. Lee, “Realistic disturbance modeling using hid-
den markov models: Applications in model-based process control,”
Journal of Process Control, vol. 19, no. 9, pp. 1438 – 1450, 2009.
[87] R. K. Wood and M. W. Berry, “Terminal composition control of a binary
distillation column,” Chemical Engineering Science, vol. 28, no. 9, pp.
1707–1717, 1973.
[88] Z. Xu, Y. Zhu, K. Han, J. Zhao, and J. Qian, “A multi-iteration pseudo-
linear regression method and an adaptive disturbance model for mpc,”
Journal of Process Control, vol. 20, no. 4, pp. 384 – 395, 2010.
140
[89] C.-M. Ying and B. Joseph, “Performance and stability analysis of lp-
mpc and qp-mpc cascade control systems,” AIChE Journal, vol. 45,
no. 7, pp. 1521–1534, 1999.
[90] J. Yu and S. J. Qin, “Mimo control performance monitoring using
left/right diagonal interactors,” Journal of Process Control, vol. 19, no. 8,
pp. 1267 – 1276, 2009.
[91] E. Zafiriou, “Robust model predictive control of processes with hard
constraints,” Computers & Chemical Engineering, vol. 14, no. 4-5, pp. 359
– 371, 1990.
[92] Y. Zhang and M. A. Henson, “A performance measure for constrained
model predictive controllers,” in Proceedings of the European control con-
ference, Karlsruhe, Germany, 1999.

Abstract (if available)

Abstract Model predictive control (MPC) is a widely used advanced process control technique in the process industry. According to the internal model principle, the internal model of MPC has to include both exact plant and disturbance models to be optimal. However, in practice, the MPC usually assumes a step-like disturbance or a fixed disturbance model. As a result, the MPC will be suboptimal when disturbance changes slowly. Moreover, it lacks a tool assessing the optimality of control performance in terms of the MPC model. ❧ In this dissertation, a new MPC disturbance adaptation method is presented. Starting from a single-input-single-output (SISO) semiconductor manufacturing process, we replaced the conventional run-to-run controller by an adaptive EWMA controller. It is shown that the plant model mismatch can be compensated by adapting the disturbance model. Analysis has been done to show that the adaptive controller is stable and converges to the optimal controller. ❧ The proposed method is then extended to multi-input-multi-output (MIMO) systems. For the ease of practical applications, the integrated moving average (IMA) model with order (1,1) is recommended. The equivalence between the IMA(1,1) parameter and the prediction error filter constant in commercial MPC has been established. Implementation of disturbance adaptation is explained. ❧ Another disturbance modeling tool is presented. It focuses on the closed-loop identification of a nonparametric disturbance model. The method incorporates the plant model information during the conversion from observer Markov parameters to system Markov parameters. ❧ A new control performance assessment method evaluating MPC model quality is then presented. Feedback invariant principle is introduced, based on which a method estimating disturbance innovations is given. A model quality index is developed as the performance benchmark, which compares prediction errors with disturbance innovations. It is shown that the model quality index related to the MPC performance index. ❧ Most industrial processes are optimized by a linear programming (LP) problem on top of the MPC. A new control performance monitoring method for cascaded LP-MPC system is developed. The block-lower-triangular interactor matrix is introduced to form a new method that is able to determine the performance benchmark based on controlled variable (CV) priorities coming from the LP results.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Non‐steady state Kalman filter for subspace identification and predictive control

PDF

Performance prediction, state estimation and production optimization of a landfill

PDF

Economic model predictive control for building energy systems

PDF

Novel multi-stage and CTLS-based model updating methods and real-time neural network-based semiactive model predictive control algorithms

PDF

Concurrent monitoring and diagnosis of process and quality faults with canonical correlation analysis

PDF

Adaptive control: transient response analysis and related problem formulations

PDF

Learning to detect and adapt to unpredicted changes

PDF

Inferential modeling and process monitoring for applications in gas and oil production

PDF

Unfalsified adaptive control with reset and bumpless transfer

PDF

Multiple model adaptive control with mixing

PDF

High-accuracy adaptive vibrational control for uncertain systems

PDF

Adaptive control with aerospace applications

PDF

Data-driven performance and fault monitoring for oil production operations

PDF

Process data analytics and monitoring based on causality analysis techniques

PDF

LQ feedback formulation for H∞ output feedback control

PDF

Statistical modeling and process data analytics for smart manufacturing

PDF

Control strategies for the regulation of blood glucose

PDF

Bumpless transfer and fading memory for adaptive switching control

PDF

Nanostructure interaction modeling and estimation for scalable nanomanufacturing

PDF

Practical adaptive control for systems with flexible modes, disturbances, and time delays

Asset Metadata

Creator Sun, Zhijie (author)

Core Title Performance monitoring and disturbance adaptation for model predictive control

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Chemical Engineering

Publication Date 08/27/2012

Defense Date 08/15/2012

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag control performance monitoring,model predictive control,OAI-PMH Harvest,run-to-run control

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Qin, S. Joe (committee chair), Safonov, Michael G. (committee member), Wang, Pin (committee member)

Creator Email sunzhijie@gmail.com,zhijiesu@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-92439

Unique identifier UC11289349

Identifier usctheses-c3-92439 (legacy record id)

Legacy Identifier etd-SunZhijie-1167.pdf

Dmrecord 92439

Document Type Dissertation

Rights Sun, Zhijie

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA