University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content NON-STEADYSTATEKALMANFILTERFORSUBSPACE
IDENTIFICATIONANDPREDICTIVECONTROL
by
YuZhao
ADissertationPresentedtothe
FACULTYOFTHEUSCGRADUATESCHOOL
UNIVERSITYOFSOUTHERNCALIFORNIA
InPartialFulfillmentofthe
RequirementsfortheDegree
DOCTOROFPHILOSOPHY
(ELECTRICALENGINEERING)
May2014
Copyright2014 YuZhao
ii
Dedication
Tomyfamily.
iii
Acknowledgments
First of all, I would like to express my sincerest gratitude to my ad-
visor, Professor S. Joe Qin. I am deeply grateful for his continuous encour-
agement,intellectualguidanceandfinancialsupportforthepastfewyears.
Mydoctoralstudycannotbecompletedwithouthisingeniousdirectionand
constructivecriticism. Ithasbeenagreatpleasuretoworkwithhimandhe
has influenced me in many aspects of my life. He guided me to acquire
knowledgebackgroundsandtodevelopnewideastosolvepracticalprob-
lems. More importantly, he taught me the innovative thinking and profes-
sionalanalyticalskillswhichwillbenefitmethroughmycareer.
IappreciatetheeffortofDr. M.SafonovfromElectricalEngineering
Department and Dr. K. Shing from the Chemical Engineering Department
forservinginmycommittee. Iwouldliketothankthemfortheirinsightful
commentsandadvices.
IthasbeenanhonorandprivilegetobeamemberofProfessorQin’s
researchgroup. Itisacollectionoftalentedanddedicatedindividualswho
granted me enjoyable and memorable experience. I am grateful to the for-
mergroupmemberswhoofferedmevaluableadvicesandpatienthelp,Car-
los, Quan, Yamin, Bo, Zhijie, Jingran and Hu. I am also grateful to all the
iv
current members for the friendship and consistent assistence, Tao, Johnny,
Yining,Zhaohui,Alisha,Yuan,QinqinandvisitingscholarProfessorLi. E-
specially, Zhijie deserves sincere appreciation for all our interactions and
collaborations. He shared his broad knowledge of MPC and enlightened
meoftheinterdisciplinaryapplicationvalueofmyresearch.
I would like to thank the USC Provost’s Fellowship, and the Center
forInteractiveSmartOilfieldTechnologies(CiSoft)forthefinancialsupport
duringmygraduatestudies.
Last but not the least, I want to thank my parents for their endless
love,careandencouragement.
v
Table of Contents
Dedication ii
Acknowledgments iii
ListofTables viii
ListofFigures ix
Abstract xii
Chapter1. Introduction 1
1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
Chapter2. SubspaceIdentificationOverview 11
2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 MathematicalPreliminaries . . . . . . . . . . . . . . . . . . . . 15
2.2.1 OrthogonalProjection . . . . . . . . . . . . . . . . . . . . 16
2.2.2 ObliqueProjection . . . . . . . . . . . . . . . . . . . . . . 18
2.3 ProblemFormulation . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3.2 AssumptionsofSubspaceIdentification . . . . . . . . . . 22
2.4 Open-loopSubspaceIdentificationMethods . . . . . . . . . . . 23
2.4.1 MOESP . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4.2 N4SID . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.3 CVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.4 UnifyingTheorem . . . . . . . . . . . . . . . . . . . . . . 29
2.5 ProblemsinSubspaceIdentification . . . . . . . . . . . . . . . . 30
Chapter 3. Progressive Parameterization in SIMs with Finite Hori-
zons 32
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
vi
3.2 ProgressiveParametrizationofSIDModels . . . . . . . . . . . 34
3.2.1 High-orderARXModels . . . . . . . . . . . . . . . . . . 34
3.2.2 SubspaceIdentificationModelswithaPredictionHori-
zon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.3 SubspaceIdentificationModelswithSpecifiedOrder . . 38
3.2.4 SubspaceIdentificationModelswithRecursiveStateS-
paceParametrization . . . . . . . . . . . . . . . . . . . . 40
3.3 ComparisonofSIDModelsatVariousStages . . . . . . . . . . 41
3.4 Numerical Comparison of SID Models via Monte-Carlo Sim-
ulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
Chapter4. NKFParameterizationbasedDisturbanceModelingand
MPCApplication: PartI 51
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.2 DynamicMatrixControl . . . . . . . . . . . . . . . . . . . . . . 53
4.2.1 Plantmodelinstepresponseform . . . . . . . . . . . . . 53
4.2.2 DMCpredictionsandoffset-freecontrol . . . . . . . . . 55
4.2.3 Plantmodelinstate-spaceform . . . . . . . . . . . . . . 56
4.3 EquivalenceBetweenLSEstimationandNKF . . . . . . . . . . 57
4.4 IdentificationofNKFDisturbanceModels . . . . . . . . . . . . 60
4.4.1 IdentificationofPMPs . . . . . . . . . . . . . . . . . . . . 60
4.4.2 IdentificationofDisturbanceModels . . . . . . . . . . . 64
4.5 DMCwithNKFDisturbanceModels . . . . . . . . . . . . . . . 68
4.6 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Chapter5. NKFParameterizationbasedDisturbanceModelingand
MPCApplications: PartII 75
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.2 DerivingtheProcessandDisturbanceModels . . . . . . . . . . 80
5.2.1 EstimatingPMPs . . . . . . . . . . . . . . . . . . . . . . . 80
5.2.2 ParameterizingProcessandDisturbanceModels . . . . 83
5.3 NKFParameterization . . . . . . . . . . . . . . . . . . . . . . . 84
5.3.1 EstimatingPMPs . . . . . . . . . . . . . . . . . . . . . . . 84
vii
5.3.2 ParameterizingSystemMarkovParameters . . . . . . . 86
5.4 MPCwithImprovedDisturbanceModels . . . . . . . . . . . . 87
5.5 CaseStudies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.5.1 Wood-BerryDistillationColumn . . . . . . . . . . . . . . 89
5.5.2 TennesseeEastmanprocess . . . . . . . . . . . . . . . . . 94
5.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Chapter6. NKFParameterizationbasedSubspaceIdentification 106
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.2 SteadyStateKalmanFilterParameterization . . . . . . . . . . . 110
6.2.1 EstimatingPredictorMarkovParameters . . . . . . . . . 110
6.2.2 ObtainingSystemMarkovParameters . . . . . . . . . . 115
6.2.3 SubspaceIdentificationwithSKFParameterization . . . 116
6.3 NKFParameterizationandEstimation . . . . . . . . . . . . . . 117
6.3.1 NKFMarkovParameters . . . . . . . . . . . . . . . . . . 117
6.3.2 FromNKFtoSystemMarkovParameters . . . . . . . . . 120
6.4 SubspaceIdentificationwithNKFParameterization . . . . . . 121
6.5 SimulationExample . . . . . . . . . . . . . . . . . . . . . . . . . 124
6.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
Chapter7. Conclusions 134
Bibliography 137
viii
List of Tables
Table3.1 EigenvaluesofHankelMatrixorWeightedHankelMatrix . 47
Table4.1 ControlPerformanceofDifferentMPCs . . . . . . . . . . . . 72
Table 5.1 Control Performance of MPCs with Different Disturbance
ModelsandVariousOrders . . . . . . . . . . . . . . . . . . . 94
Table 5.2 Control Performance of MPCs with Different Disturbance
Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Table6.1 MSEforpoleestimation . . . . . . . . . . . . . . . . . . . . . 128
Table6.2 Varianceforpoleestimation . . . . . . . . . . . . . . . . . . . 128
Table6.3 MSEandVarianceforvaryingp . . . . . . . . . . . . . . . . . 131
ix
List of Figures
Figure 2.1 The left hand side shows the subspace identification ap-
proach;therighthandsideistheclassicalapproach. . . . . 13
Figure2.2 Orthogonalprojection. . . . . . . . . . . . . . . . . . . . . . 17
Figure 3.1 Various stages of subspace identification leading to the
estimatesofMarkovparameters. . . . . . . . . . . . . . . . . 48
Figure 3.2 Estimated PMPs after stage 1 and 2. Red curves are true
PMPs and dots are estimated PMPs. a) estimated PMPs
of
¯
G from HOARX after stage 1; b) estimated PMPs of
¯
H
from HOARX after stage 1; c) estimated PMPs of
¯
G from
HOARXafterstage2(reducedrank);d)estimatedPMPsof
¯
H fromHOARXafterstage2(reducedrank); e)estimated
PMPsof
¯
Gfrom
¯
H
LS
fp
;f)estimatedPMPsof
¯
H from
¯
H
LS
fp
;g)
estimated PMPs of
¯
G from
¯
H
CVA
fp
; h) estimated PMPs of
¯
H
from
¯
H
CVA
fp
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
Figure 3.3 Estimated PMPs from the Hankel matrices after stage 2.
Thefirstf PMPsareextractedfromthefirstcolumnofthe
Hankel matrices; the last (p−1) PMPs are extracted from
the last row of the Hankel matrices. a) estimated PMPs of
¯
G without CVA weighting; b) estimated PMPs of
¯
G with
CVAweighting.. . . . . . . . . . . . . . . . . . . . . . . . . . 50
x
Figure 3.4 Estimated PMPs of
¯
G after stage 3. a) estimated PMP-
s of
¯
G using SSARX without weighting, system matrices
estimated using
¯
Γ
f
and
¯
L
p
; b) estimated PMPs of
¯
G us-
ingSSARXwithCVAweighting,systemmatricesestimat-
ed using
¯
Γ
f
and
¯
L
p
; c) estimated PMPs of
¯
G using SSARX
without weighting, system matrices estimated using state
sequence;d)estimatedPMPsof
¯
GusingSSARXwithCVA
weighting,systemmatricesestimatedusingstatesequence. 50
Figure 4.1 Identified system Markov parameters of the disturbance
modelinCaseI. . . . . . . . . . . . . . . . . . . . . . . . . . 71
Figure4.2 SystemMarkovparametersofthedisturbancemodel. . . . 71
Figure 4.3 Identified system Markov parameters of the disturbance
modelinCaseII . . . . . . . . . . . . . . . . . . . . . . . . . 73
Figure5.1 Impulseresponseofthenominaldisturbancemodel. . . . . 91
Figure5.2 Impulseresponseofstep-like,SKFandNKFbaseddistur-
bance: order=8. . . . . . . . . . . . . . . . . . . . . . . . . . 92
Figure5.3 Impulseresponseofstep-like,SKFandNKFbaseddistur-
bance: order=30. . . . . . . . . . . . . . . . . . . . . . . . . 93
Figure5.4 SchematicdiagramoftheTennesseeEastmanprocess. . . . 96
Figure5.5 Open-loopresponseoftheselectedMV/CVpair. . . . . . . 98
Figure5.6 ImpulseresponseofplantmodelofMPCcontrollingTEP. . 99
Figure 5.7 Impulse response of disturbance model: step-like distur-
bancemodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
Figure 5.8 Impulse response of disturbance model: SKF parameter-
izeddisturbancemodel. . . . . . . . . . . . . . . . . . . . . . 100
xi
Figure 5.9 Impulse response of disturbance model: NKF parameter-
izeddisturbancemodel. . . . . . . . . . . . . . . . . . . . . . 100
Figure5.10ControlresultsofTEPbyMPCwithstep-likedisturbance
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
Figure 5.11Control results of TEP by MPC with SKF parameterized
disturbancemodel. . . . . . . . . . . . . . . . . . . . . . . . . 103
Figure 5.12Control results of TEP by MPC with NKF parameterized
disturbancemodel. . . . . . . . . . . . . . . . . . . . . . . . . 104
Figure6.1 ScatteredplotoftheestimatedpolesfromNKFID,p =f =
20. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
Figure6.2 ScatteredplotoftheestimatedpolesfromOKID,p =f = 20.126
Figure 6.3 Bode plot of NKFID and OKID. The red solid line is the
true values. The blue and green dashed lines are the esti-
matedvaluesaveragedfrom50runsbasedonNKFIDand
OKID,respectively. . . . . . . . . . . . . . . . . . . . . . . . . 127
Figure 6.4 Identification result from NKFID with varyingp. The left
columnistheBodemagnitudeplots;Thesolidlinearethe
true values and the dashed line are the NKFID estimated
values averaged from 50 runs. The right column is the
scatterplotsoftheeigenvaluesofestimatedA . . . . . . . . 129
Figure 6.5 Identification result from OKID with varying p. The left
column is the Bode magnitude plots for different p; The
solid line are the true values and the dashed line are the
OKID estimated values averaged from 50 runs. The right
columnisthescatterplotsoftheeigenvaluesofestimatedA130
xii
Abstract
Subspaceidentificationmethods(SIMs)havedrawntremendousre-
search interest and have been widely recognized for its computational ef-
ficiency and strong capacity in identifying MIMO systems. Most existing
SIMsarebasedonsteady-stateKalmanfilter(SKF)modelswhichimplyin-
finite horizons to be used in forming the dynamic time lags for the model.
However,inpractice,finitedatalengthistheonlyoptionthatdictateseven
shorter time horizons for the number of time lags. As a result, the SKF pa-
rameterizationforSIMscanleadtosuboptimalsolutionwithfinitenumber
ofdata.
Inthisdissertation,anovelnonsteady-stateKalmanfilter(NKF)pa-
rameterization based SIM (NKFID) is proposed to address two aspects: fi-
nitedatahorizonandclosed-loopidentification. Aprogressiveparameteri-
zationframeworkisfirstbuilttointerpretthemodelsineachstageofSIM.
The model of the first stage is shown as non-parametric high order ARX
model,whichhasappliedvalueinmodelpredictivecontrol(MPC).
In practice, MPC usually assumes a step-like disturbance model
which is insufficient to describe complex dynamics. The finite impulse
response (FIR) model of the disturbance is recommended in this work.
xiii
Two data based parameterization methods are proposed to parameterize
the FIR form disturbance model. The disturbance modeling schemes with
SKF parameterization is first developed to improve MPC performance.
Considering finite data length, the NKF parameterization based modeling
scheme is proposed and improved, and its suitability for finite data is ver-
ified theoretically. The general conversion relation between the predictor
Markov parameters (PMPs) and the system Markov parameters (SMPs)
is developed under the NKF framework, based on which the disturbance
modelisderivedandincorporatedtoMPC.
The NKFID is proposed as a novel SIM under the NKF framework.
The algorithm uses the NKF parameterization hence optimal solution is
guaranteed even with limited data horizons. A routinely used projection
that removes the impact of the future input is avoided to make NKFID
applicable to closed-loop identification. The performance of the proposed
disturbance modeling and subspace identification methods are demon-
strated by several simulation examples and the benefits are verified by the
comparisonwithexistingmethods.
1
Chapter 1
Introduction
1.1 Overview
Mathematical models of dynamic systems are built to capture the
system dynamic behavior for use in a wide range of industries including
chemical, biological, petroleum, automotive and aerospace. As a vital part
of the engineering applications, mathematical models of the dynamic sys-
tem are the bases for control system design, process analysis, optimization
and process monitoring. There are generally two approaches for building
mathematical models. The first approach is the first principles method, or
the white-box method. It requires very solid understanding of the system
and the knowledge of all related physical and chemical parameters. The
secondapproachisthesystemidentification,ortheblack-boxmethod. Sys-
temidentificationdealswiththeproblemofbuildingmathematicalmodels
ofdynamicsystemsbasedontheobserveddata. Itispartofthebasicscien-
tificmethodologyandithasabroadrangeofapplications[44].
System identification has a long history of development. It has
been thriving after the success of the model-based control design in 1960s
[20]. Major achievements for the past decades have been summarized
2
in the classic book by Ljung [44]. The traditional system identification
can be interpreted as nonlinear optimization problems requiring iterative
nonlinear search for the solution. The two best-known traditional system
identification methods are Prediction Error Method (PEM) and Maximum
LikelihoodMethod(ML)[62,6,44]. PEMdeterminesthemodelparameters
by minimizing the prediction error, while ML determines the parameters
by maximizing the likelihood function. Although those methods are ca-
pable of handling SISO systems, they are difficult to be directly applied to
MIMOsystems. Thedifficultycomesfromthenumericallyill-conditioning,
nonconvexity, and possible local minima existing in the identification of
multivariablepolynominalmodels.
Subspaceidentificationmethods,asanotherclassofsystemidentifi-
cationmethods,emergedinthe1990sandhavegainedgreatresearchinter-
estforitscomputationalefficiency,numericalreliability,andstrongcapacity
inidentifyingMIMOsystems. TheSIMsoriginatedfromtheclassicaldeter-
ministic realization theory by Ho and Kalman [25] and the stochastic real-
izationtheorybyAkaike[1]. TheearlydevelopmentsofSIMsareclassified
tothesocalledrealizationbasedSIMs[72]. TherealizationbasedSIMsrely
on the impulse response to recover the system matrices. The applications
arelimitedduetothedifficultyinobtainingsystemimpulseresponse.
On the other hand, the direct SIMs, based on the input-output data
3
instead,havedevelopedrapidlyforthelasttwodecadesandhaveachieved
remarkable success in applications. The best-known methods in this cat-
egory are multivariate output-error state space (MOESP) [69], numerical
algorithm of subspace state space system identification (N4SID) [65] and
canonicalvariate analysis (CVA)[39]. Despite thedifferencein algorithms,
theyallcontaintwomajorsteps: inthefirststepthedataisprojectedtocer-
tainsubspacestogeneratetheestimateoftheextendedobservabilitymatrix
and/or the Kalman state series. In the second step the system matrices are
retrievedfromeithertheextendedobservabilitymatrixortheKalmanstate
estimates. The unifying theorem by Van Overschee and DeMoor pointed
out that the difference among the three methods can be interpreted as em-
ploying different weighting matrices in the singular value decomposition
[66].
Comparedtothetraditionalsystemidentificationmethods,theSIM-
s are suitable for MIMO system identification. The algorithm provides a
numerically reliable way to identify state space models without nonlinear
optimization. AnotheradvantageofSIMsisitslessstrictrequirementofthe
prior knowledge. In fact, the system orderis the only necessary parameter
to choose for SIMs. This feature provides significant flexibility in practice
since the other model information, such as the model structure, is usually
unavailablewithoutsystemknowledge.
4
Although the traditional SIMs (e.g. CVA, N4SID, MOESP) are
suitable for the open-loop data, they are biased for the closed-loop system
identification. The closed-loop identification has attracted special interest
from engineers since carrying out the identification experiment under
closed-loop condition is often desirable for safety and quality reasons [53].
The difficulty in applying the traditional SIMs to the closed-loop data
comesfromthecorrelationbetweenthenoiseandinputunderclosed-loop
condition[44]. MuchefforthasbeenmadetodeveloptheSIMssuitablefor
the closed-loop data. Verhaegen [68] suggested to identify the closed-loop
model via an overall open-loop state space model identification followed
by model reduction. Ljung and McKelvey [45] proposed the recursive ap-
proachbasedonautoregressivemovingaverage(ARX)modelsasafeasible
way to handle closed-loop data. The recent development of closed-loop
SIMs includes: the use of parity space and principal component analysis
(PCA) and its further enhancement [73, 28]; row-wise regression on the
extended state space model to bypass the projection [7]; the removal of
futureinputbymanipulatingMarkovparameterestimates[29].
A complete subspace identification procedure goes through several
stages[53]:
1. PerformleastsquaresestimationonthehighorderARXmodeltoob-
tain predictor Markov parameter estimates, which by themselves are
5
coefficients of the FIR model. This step is non-parametric and has its
ownapplicationsinmodelpredictivecontrol.
2. PerformSVDontheHankelmatrixorweightedHankelmatrix,lead-
ingtoextendedobservabilityandcontrollabilitymatrixestimates. The
reducedordernon-parametricmodelsisobtainedinthisstepforstate
estimationandoutputprediction.
3. Perform least squares estimation to retrieve system matrices, leading
tothestate-spaceparametricmodels.
Theestimationofthefiniteimpulseresponse(FIR)modelofthepro-
cess and disturbance, as the first step of most SIMs, has important applied
value in the process control area. Model predictive control, being the most
popularadvancedprocesscontroltechnique[36],caneasilyincorporatethe
FIR formed process and disturbance model for optimization and predic-
tion. MPC is one of the most successful process control schemes and it is
well-known for being capable of handling multivariable constrained con-
trol problems. As reported by Qin and Badgwell [56], more than 4,500 in-
dustrial linear MPCs emerged in a broad range of industries from refining
toaerospace.
TheMPCmodelaccuracyisoneofthekeyfactorsaffectingthecon-
trol performance. The MPC model includes both the process and distur-
6
bance models. The process model describes the dynamics between manip-
ulatedvariables(MV)andcontrolledvariables(CV),whilethedisturbance
model shapes the unmeasured disturbances. The optimal control signals
are computed based on the predictions of the future outputs by the MPC
model. As a widely used MPC technique, dynamic matrix control (DMC)
[10, 19] adopts a step-like disturbance model which was recommended as
thefavorablechoice[16,17]. Unfortunately,thestep-likedisturbancemodel
is not always sufficiently general to give desirable predictions. The neces-
sity of more complete disturbance dynamics was proposed by Francis and
Wonham [14] and more complex disturbance models have been suggested
inordertoachievesuperiorcontrolperformance[48,15,43]. Therecentre-
searchonthegeneralformofthedisturbancemodelderivedthedetectabil-
ity requirement as the criterion for disturbance modeling [49, 51]. There is
still plenty space for improvement of the disturbance modeling, especially
forthesystemwithcomplicateddisturbancedynamics.
When identifying the disturbance model, the first step of most esti-
mation approaches is to identify the high-order ARX (HOARX) model re-
gardless of the model structure, e.g., FIR model or state-space model. In
this step, most existing methods assume the order of HOARX model as
infinite so that the predictor Markov parameters for high order terms can
approach zero. Based on this assumption, the identified HOARX model
7
is further parameterized by the steady-state Kalman filter (SKF). Unfortu-
nately, it is inconsistent to assume infinite order while estimating with the
finite length of data. Therefore, how to parameterize the HOARX mod-
el with limited order using closed-loop data still needs further investiga-
tion. Under open-loop condition, previous work by Sorenson [63] showed
that the nonsteady-state Kalman filter (NKF) is a recursive solution to the
finite order least-squares (LS) problems. Under the subspace framework,
Van Overschee and De Moor [65] proved that the LS regression of the fu-
tureoutputonthepastinputandoutputdataactuallyrepresentsabankof
NKF.Huang,DingandQin[28]alsointerpretedthefuturestateastheNKF
state. However,theseresultsareonlyapplicabletotheopen-loopdata. Al-
though NKF parameterization is suggested, how to use it to parameterize
theprocessanddisturbancemodelisstillanopenquestion. Thismotivates
ustobuildtheNKFframeworkfortheprocessanddisturbancemodeling.
Sincethepracticalsubspaceidentificationproblemmustdealwithfi-
nitelengthofdata,theSKFbasedparameterizationisbiased. Unfortunate-
ly,mostexistingSIMsemploytheSKFintheparameterizationofthepredic-
tor Markov parameters for derivation simplicity. The representative SIMs
based on SKF framework include the OKID [34], which directly retrieves
processanddisturbancemodelsfromthepredictorMarkovparameteresti-
mates;theSSARXmethod[29],whicheliminatestheinputnoisecorrelation
8
problem by removing the future input via pre-estimated Markov parame-
ters. However,theapplicationofthoseSIMstothefinitedatasetmaylead
tosuboptimalestimates[57]. ThisinspiresustodevelopthenovelSIMun-
der NKF parameterization framework to address the practical finite data
lengthidentificationproblem.
1.2 Outline
This dissertation focuses on developing the NKF parameterization
based SIM. In this work, a novel NKF parameterization based process and
disturbance modeling scheme is first developed. Based on this, the NK-
F parameterization based subspace identification (NKFID) is proposed to
address the issues and drawbacks of the existing SIMs. In addition, the
proposed disturbance modeling method is applied to MPC to improve the
controlperformance. Thisdissertationisorganizedasfollows:
Chapter 2 is dedicated to the review of the subspace identification
developmentduringthelastfewdecades. Thenecessarybackgroundsand
preliminaries are introduced. Three most representative open-loop SIMs
arediscussedwithgeometricalinterpretation,followedbythesummaryof
theadvantagesandcurrentproblemsinsubspaceidentification.
In Chapter 3, the progressive parameterization framework is intro-
duced to interpret the models used in each step of SIM. It also discussed
how the progressively parameterized models lead to the recursive state s-
9
pace models. In addition, the intermediate non-recursive model is shown
tobeusefulforthepurposeofstateestimation,faultdetectionandcontrol.
Chapter 4 presents the FIR formed disturbance models with both
SKF and NKF parameterizations. The existing step-like disturbance for
DMC is compared with to verify the control performance improvement
achieved by the proposed disturbance models. A general conversion
relation between system and predictor Markov parameters is developed
in order to obtain the process and disturbance FIR models. The proposed
disturbance modeling scheme is compatible with the NKF structure but
requiresprocessanddisturbancemodelsharingpoles.
In Chapter 5, a general modeling scheme based on NKF parameter-
ization is studied and further incorporated to DMC. It is shown that the
Markov parameter estimates from least squares can be parameterized by
NKF under closed-loop condition. The conversion from predictor Markov
parameterstosystemMarkovparametersisthendevelopedundertheNK-
Fframework. TheNKFparameterizationschemeisappliedtodisturbance
modelingforDMC.Controlperformanceimprovementisdemonstratedby
the simulation on Wood-Berry distillation column and Tennessee Eastman
process.
TheNKFparameterizationbasedsubspaceidentification(NKFID)is
proposed in Chapter 6. The NKF parameterization framework is built for
subspaceidentification. TheproposedNKFIDiscomparedwiththeexisting
10
OKID to demonstrate the effect of the data horizon on the identification.
TheadvantageofNKFIDisdemonstratedthroughabenchmarkproblem.
Chapter7givestheconcludingremarksforthedissertation.
11
Chapter 2
Subspace Identification Overview
2.1 Overview
In general, SIMs can be classified into two major categories: the re-
alization based SIMs and direct SIMs ([72]). The origin of the realization
basedsubspaceidentificationmethodscanbetracedbacktotherealization
algorithmofHoandKalman[25],whichprovidedthesolutiontothedeter-
ministic realization problem. Furthermore, Akaike [1] interpreted various
realizationalgorithmsfromastochasticperspective. Basedonthoseresult-
s, the subspace identification has developed from deterministic realization
theory to the so-called combined deterministic stochastic realization theo-
ry. However, due to the difficulty in obtaining accurate impulse response,
therealizationbasedmethodscannotbeeasilyimplementedinpracticeand
thustheyattractrelativelylittleattentionfromresearchers.
Over the past decade, much research effort has been made towards
thedirectSIMs([8],[12],[39],[69],[71]). Unliketherealizationbasedmeth-
odswhichrelyonimpulseresponses,thedirectSIMsestimatetheextended
observability matrix or the state sequence directly from the input and out-
put data. There are a number of representative algorithms, including the
12
multivariateoutput-errorstatespace(MOESP)[69],numericalalgorithmof
subspacestatespacesystemidentification(N4SID)[65]andcanonicalvari-
ate analysis (CVA) [39]. Although those methods vary in algorithms, the
difference can be interpreted as using different weighting matrices in the
singular value decomposition ([66]). Despite the differences among them,
those subspace identification methods consist two basic steps. In the first
step,theextendedobservabilitymatrixandablocktriangularToeplitzma-
trix containing Markov parameters are estimated by projecting the data to
certain subspaces. The second step is to retrieve the system matrices and
disturbancecharacteristicsfromtheestimation.
Therearemanydirectionsintheresearchofsubspaceidentification.
The important and ongoing topics includes asymptotic analysis and
consistency condition analysis ([2], [3], [11], [31], [38], [52]); analysis of the
weightingmatrices([22],[4]);exploringnewSIMsandapplyingtheSIMsto
closed-loopdata([45],[68],[67],[34],[29],[73],[28],[7],[54]);extendingthe
methods to unstable, linear time variant (LTI) systems, nonlinear systems
anderror-in-variablesystems([8],[42]).
The SIM aims at the state-space model estimation and it has strong
capacityformultivariableidentification. ThedifferencesbetweenSIMsand
theclassicalmethodsareanalyzedin[37]andaredemonstratedinFig. 2.1.
Intheclassicalmethodssuchasthepredictionerrormethod(PEM)[44],the
13
Input output data
Kalman state
sequence
System matrices
System matrices
Kalman state
sequence
Orthogonal or
oblique projection
Classical
identification
Least squares Kalman filter
Figure 2.1: The left hand side shows the subspace identification approach;
therighthandsideistheclassicalapproach.
transfer function model is identified, a state space model, if desired, could
then be derived through realization schemes. On contrary, in the SIMs the
Kalman state estimates are first constructed by projecting the input output
data to certain subspaces. A state space model is further retrieved and the
transferfunctionmodelcouldbeobtainedthroughtransformation.
The SIMs have notable advantages compared to the PEM, which is
mostly based on transfer function models. First of all, SIMs are simple for
MIMO processes. The SIMs provide numerically reliable ways to identify
the state-space models with little prior knowledge, thereby removing the
main limitations of PEM in terms of how to parameterize a multivariable
14
system with a suitable canonical form. For SIMs, the only required knowl-
edgeistheorderofthesystem. Inthisway, SIMscompletelyavoidthethe
canonical parameterizations. Furthermore, based on reliable numerical al-
gorithmsofthesingularvaluedecompositionortheQRfactorization,SIM-
s do not need nonlinear optimization techniques nor do they impose any
canonical form on the system. This indicates that SIMs can be convenient-
ly applied to the MIMO as well as the SISO systems. In addition, the state
space formulation is convenient for estimation, filtering and prediction in
industrialprocessapplications.
These advantages make SIMs powerful tools for MIMO system i-
dentificationandon-lineapplication, such asmodelpredictivecontroland
adaptivecontrol.
Although SIMs are regarded as one of the most significant achieve-
ment in system identification, some drawbacks are also discovered ([54]).
Threeofthemarelistedbelow:
1. Most existing SIMs are developed under the framework of SKF pa-
rameterization, which inherently assumes that the data horizon is in-
finite. Thisassumptionisnotvalidinpractice.
2. Assume the system dynamics is represented by the following state
15
spacemodel,
x
k+1
= Ax
k
+Bu
k
+w
k
(2.1)
y
k
= Cx
k
+Du
k
+v
k
(2.2)
The estimation of B and D is more problematic than that of A and
C, which is reflected in the poor estimation of zeros and steady state
gains.
3. The application of SIMs to closed-loop data is still challenging, even
for the data satisfying identifiability conditions for traditional meth-
odssuchasPEMs.
This chapter is devoted to general formulation and background in-
formationoftheSIMs. InSection2.1,somemathematicalpreliminariesare
given. In Section 2.2, the useful vectors and matrices of the subspace i-
dentification are introduced. Section 2.3 reviews three most representative
open-loop SIMs. The unifying theorem is also introduced to interpret the
differencesamongthem. Section2.4discussesthechallengesinclosed-loop
subspaceidentification.
2.2 Mathematical Preliminaries
Mostsubspaceidentificationalgorithmsarebasedongeometriccon-
cepts. Some system characteristics can be interpreted as geometric manip-
16
ulation of the row spaces of certain matrices. In this section, the basic con-
cepts of orthogonal and oblique projection are reviewed. In the following
weassumethatthematricesA∈R
pN
,B ∈R
qN
andC ∈R
rN
aregiven.
2.2.1 Orthogonal Projection
WedefineΠ
B
astheopertatorwhichprojectstherowspaceofama-
trixontotherowspaceofthematrixB,
Π
B
,B
T
(
BB
T
)
y
B (2.3)
where (•)
y
denotes the Moore-Penrose pseudo-inverse of the matrix of (•).
A=B is the shorthand for the projection of the row space of A on the row
spaceofthematrixB,
A=B , AΠ
B
(2.4)
= AB
T
(
BB
T
)
y
B (2.5)
As indicated in the Fig. 2.2, it is obvious that the result after the projection
operator Π
B
lies in the row space of B. The orthogonal projection can be
considered as the geometric interpretation of the least squares estimation
andtheQRdecompositionistheidealnumericaltoolforit.
AssumethematrixAandB arelinearlyrelatedwithV representing
thematrixofnoise,
A = ΘB +V: (2.6)
17
% $
!

% $ $ Figure2.2: Orthogonalprojection.
Byminimizingtheobjectivefunction
J =∥A−ΘB∥
2
F
(2.7)
where∥•∥
F
istheF-norm,wehavetheleastsquaressolution
ˆ
Θ =AB
T
(
BB
T
)
y
(2.8)
andthemodelpredictionis
ˆ
A =
ˆ
ΘB
= AB
T
(
BB
T
)
y
B
= AΠ
B
: (2.9)
ThismeansthatthepreditionofmatrixAbasedontheleastsquaresestima-
tionistheprojectionoftherowspaceofmatrixAtotherowspaceofmatrix
18
B. Theleastsquaresresidualis
˜
A = A−
ˆ
A
= A(I−Π
B
)
= AΠ
?
B
(2.10)
where
Π
?
B
= I−Π
B
= I−B
T
(
BB
T
)
y
B (2.11)
is defined as the projection to the orthogonal complement of B. Similar to
A=B, we can useA=B
?
as the shorthand ofAΠ
?
B
. It is easy to verify that
ˆ
A
and
˜
Aareorthogonal, andthematrixAcanbedecomposedintothosetwo
orthogonalcomponents
A =AΠ
B
+AΠ
?
B
:
2.2.2 Oblique Projection
Instead of decomposing the matrix A into two orthogonal com-
ponents, it can be decomposed as the linear combination of two non-
orthogonal matrices. To illustrate this we need to introduce the definition
of oblique projection. The oblique projection of the row space of matrix A
alongtherowspaceofmatrixB ontherowspaceofmatrixC isdefinedas
follows
A=
B
C,
[
A=B
?
][
C=B
?
]
y
C: (2.12)
19
It is worth noting that when B = 0 or when the row space of matrix B
is orthogonal to the row space of matrix C, the oblique projection can be
reducedtotheorthogonalprojection
A=
B
C =A=C: (2.13)
Somepropertiescanbederivedstraightforwardly,
B=
B
C = 0; (2.14)
C=
B
C = C: (2.15)
2.3 Problem Formulation
In this section, we start from the basic state space model to derive
the important equations and matrices for the subspace identification. As
three basic state space models for all subspace methods, the process form,
the innovation form and the predictor form are reviewed. Based on them,
the stacked input output variables are introduced, followed by the deriva-
tion of the extended state space models. Furthermore, the Hankel matrix
and extended observability matrix which play very important role in the
subspaceidentificationareformed.
2.3.1 Preliminaries
The purpose of system identification is to estimate an appropriate
dynamic model for a properly collected time series of system input and
20
output data{u
k
;y
k
} which can be collected from open-loop or closed-loop
conditions. Assumingthedatasequence{u
k
;y
k
}isgeneratedbythefollow-
ingstatespacemodel,
x
k+1
= Ax
k
+Bu
k
+w
k
(2.16)
y
k
= Cx
k
+Du
k
+v
k
(2.17)
wherey
k
∈R
ny
,x
k
∈R
n
,u
k
∈R
nu
, are the system output, state and input,
respectively. w
k
∈R
n
representsthestatenoisewhilev
k
∈R
ny
istheoutput
measurementnoise. Thestatespacemodelin(2.16)and(2.17)istheprocess
form.
AKalmanfiltercanbedesignedforthissystemtoestimatethestate
variables. Ifthesystemisobservable,wecandesignastablestateobserver
asfollows:
ˆ x
k+1
=Aˆ x
k
+Bu
k
+K(y
k
−Cˆ x
k
−Du
k
) (2.18)
whereK isthesteadystateKalmangainwhichcanbeobtainedthroughthe
algebraicRicattiequation. TheinnovationoftheKalmanfilterisdefinedas
e
k
=y
k
−Cˆ x
k
−Du
k
(2.19)
Intherestofthisdissertation,weusex
k
insteadof ˆ x
k
forconvenience,and
wehavethefollowinginnovationform,
x
k+1
= Ax
k
+Bu
k
+Ke
k
(2.20)
y
k
= Cx
k
+Du
k
+e
k
(2.21)
21
The system in (2.20) and (2.21) can also be represented by the following
equivalentpredictorforrm,
x
k+1
=
¯
Aˆ x
k
+Bu
k
+Ky
k
(2.22)
y
k
= Cx
k
+Du
k
+e
k
(2.23)
where
¯
A =A−KC.
Althoughallthreebasicstatespacemodelformscanmodelthesys-
tem dynamics accurately, different subspace methods use different forms.
Therearesomedifferencesamongthem. Theprocessformandtheinnova-
tion form use the system matricesA, B, C andD while the predictor form
usesthe
¯
A,B,K andC,D.
Based on the process form in (2.16) and (2.17), the extended state
spacemodelcanbederivedbyiteration,
y
f
(k) = Γ
f
x
k
+Hu
f
(k)+n
f
(k) (2.24)
wherethestackedvectory
f
(k) =
[
y
T
k
y
T
k+1
::: y
T
k+f1
]
T
representsthe
future output. The future input u
f
(k) and noise n
f
(k) are defined in the
same way. The contributions from the process and output noise are both
included in n
f
(k). Similar to y
f
(k), the past output y
p
(k) is defined as
y
p
(k) =
[
y
T
kp
::: y
T
k2
y
T
k1
]
T
.
22
TheextendedobservabilitymatrixΓ
f
in(2.24)isformedasfollows
Γ
f
=





C
CA
.
.
.
CA
f1





; (2.25)
andthematrixH hastheToeplitzstructure
H =





D 0 ··· 0
CB D ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
CA
f2
B ··· CB D





; (2.26)
In practice, the data matrices are used instead of the stacked vector
variables. Thefutureinputandoutputdatamatricesareasfollowing
Y
f
=
[
y
f
(k) y
f
(k+1) ··· y
f
(k+N −1)
]
; (2.27)
U
f
=
[
u
f
(k) u
f
(k+1) ··· u
f
(k+N −1)
]
: (2.28)
Thecorrespondingextendedstatespacemodelwithdatamatricesis
Y
f
= Γ
f
X
k
+HU
f
+N
f
(2.29)
wheretheextendedstateX
k
=
[
x
k
x
k+1
::: x
k+N1
]
:
2.3.2 Assumptions of Subspace Identification
Severalgeneralassumptionsforthesubspaceidentificationweredis-
cussed in [22] and [31]. In fact, those assumptions are considered to be
standard for most SIMs. We introduce those assumptions based on the in-
novationformin(2.20)and(2.21).
23
1. Theinnovationprocesse
k
isastationary,zero-mean,whitenoisepro-
cesswithsecondordermoments
E
{
e
t
e
T
t
}
=R
e

t;s
where
t;s
istheKroneckerdelta.
2. The input u
k
is a quasi-stationary signal. The correlation function is
definedby
R
u
() = lim
N!1
1
N
N
∑
t=1
E
{
u
t+
u
T
t
}
The inputu
k
and the innovation processe
k
are assumed to be uncor-
relatedforopen-loopdata.
3. Thesystemisasymptoticallystable,i.e.,theeigenvaluesofAliestrick-
lyinsidetheunitcircle.
4. (A;C)isobservableandthepair
(
A;
[
B K
])
isreachable.
2.4 Open-loop Subspace Identification Methods
In this section, we introduce three most representative open-loop
SIMs: MOESP,N4SIDandCVA.Thefocuswillbegiventotheestimationof
theextendedobservabilitymatrixandthestate,whilethedetailsofretriev-
ing system matrices will be saved for later chapters. The unifying theorem
is also introduced to provide insights about the unified formulation and
therelationsamongtheSIMs.
24
2.4.1 MOESP
The original MOESP is based on the estimation of the extended ob-
servability matrix. As the first step, the future output is projected onto the
orthogonalspaceofthefutureinput,andtheextendedobservabilitymatrix
is extracted from the SVD of the projection. The system matrices A and C
are extracted directly from the extended observability matrix, whileB and
D areobtainedbytheleastsquaresestimation.
Based on the extended state space model in (2.29), the first step of
MOESPistoeliminatethefutureinputfromoutputbyprojectingdataonto
theorthogonalspaceofU
f
Y
f
Π
?
U
f
= Γ
f
X
k
Π
?
U
f
+N
f
Π
?
U
f
(2.30)
Oneimportantfeatureoftheopen-loopdataisthatthefuturenoiseandthe
future input are uncorrelated, and thusN
f
is not affected by the projection
ontotheorthogonalspaceofU
f
. Therefore,wecansimplify(2.30)to
Y
f
Π
?
U
f
= Γ
f
X
k
Π
?
U
f
+N
f
(2.31)
TheprojectioncanbesimplycalculatedthroughtheQRfactorization,
[
U
f
Y
f
]
=
[
R
11
0
R
21
R
22
][
Q
T
1
Q
T
2
]
(2.32)
and
Y
f
Π
?
U
f
=R
22
Q
T
2
(2.33)
25
TheextendedobservabilitymatrixΓ
f
isextractedbytheSVDofR
22
andthe
numberofdominantsingularvaluesofR
22
determinesthesystemorder
R
22
= USV
T
(2.34)
=
[
U
1
U
2
]
[
S
1
0
0 S
2
][
V
T
1
V
T
2
]
(2.35)
Theextendedobservabilitymatrixisestimatedas
ˆ
Γ
f
=U
1
(2.36)
The original MOESP is biased when the process and output nois-
es are not independent. To address this, a more practical algorithm (PO-
MOESP) is proposed ([69]). The PO-MOESP method constructs the instru-
mentalvariablewiththepastinputoutputdataandremovethenoiseterm
byprojectingontotherowspaceofZ
p
whereZ
p
=
[
U
T
p
Y
T
p
]
T
.
Y
f
Π
?
U
f
Π
Zp
= Γ
f
X
k
Π
?
U
f
Π
Zp
(2.37)
TheprojectioncanalsobeobtainedfromtheQRfactorization,


U
f
Z
p
Y
f


=


R
11
0 0
R
21
R
22
0
R
31
R
32
R
33




Q
T
1
Q
T
2
Q
T
3


Similartotheequationin(2.33),itcanbeobtainedthat
Y
f
Π
?
U
f
Π
Zp
=R
32
Q
T
2
(2.38)
TheextendedobservabilitymatrixΓ
f
canbeestimatedfromtheSVDofR
32
,
R
32
=
[
U
1
U
2
]
[
S
1
0
0 S
2
][
V
T
1
V
T
2
]
26
and
ˆ
Γ
f
=U
1
S
1=2
1
(2.39)
2.4.2 N4SID
The N4SID algorithm is based on the innovation form in (2.20) and
(2.21). The extended state space model of the innovation form can be ob-
tainedbyiteration
Y
f
= Γ
f
X
k
+HU
f
+GE
f
(2.40)
whereGisformedinthesamewayasH exceptthatthematrixDisreplaced
by K. Similar to the MOESP, N4SID projects Y
f
onto the orthogonal space
ofU
f
toremovetheU
f
term,
Y
f
Π
?
U
f
= Γ
f
X
k
Π
?
U
f
+GE
f
Π
?
U
f
(2.41)
= Γ
f
X
k
Π
?
U
f
+GE
f
(2.42)
andZ
T
p
isusedastheinstrumentalvariabletoremovethenoiseterm
Y
f
Π
?
U
f
Z
T
p
= Γ
f
X
k
Π
?
U
f
Z
T
p
(2.43)
The SVD is performed on Y
f
Π
?
U
f
Z
T
p
(
Z
p
Π
?
U
f
Z
T
p
)
1
Z
p
to obtain the
estimateofΓ
f
Y
f
Π
?
U
f
Z
T
p
(
Z
p
Π
?
U
f
Z
T
p
)
1
Z
p
=
[
U
1
U
2
]
[
S
1
0
0 S
2
][
V
T
1
V
T
2
]
andthebalancedrealizationofextendedobservabilitymatrixis
ˆ
Γ
f
=U
1
S
1=2
1
(2.44)
27
It is worth noting that N4SID has the geometric interpretation as
oblique projection. It can be interpreted as an oblique projection of the
future output Y
f
along the future input U
f
onto the past input output Z
p
.
The SVD of the oblique projection is further performed to extract the state
estimates.
Y
f
=
[
Z
p
U
f
]
=
[
L
1
L
2
]
[
Z
p
U
f
]
(2.45)
andthus
Y
f
=
U
f
Z
p
=L
1
Z
p
(2.46)
Basedontheobliqueprojectionpropertyin(2.12),
Y
f
=
U
f
Z
p
= Y
f
Π
?
U
f
[
Z
p
Π
?
U
f
]
y
Z
p
(2.47)
= Y
f
Π
?
U
f
Z
T
p
(
Z
p
Π
?
U
f
Z
T
p
)
1
Z
p
(2.48)
This means in the N4SID, the SVD is performed on the oblique projection
Y
f
=
U
f
Z
p
.
2.4.3 CVA
CVA was originally developed by Larimore ([39]) and the software
Adaptix based on this algorithm is available for model identification. The
key step of CVA is to estimate the state variables by Canonical Correlation
Analysis(CCA),whichwasdevelopedbyHotelling([26]). Theideabehind
CCA is to find the most correlated linear combinations from two sets of
variables. InCVAitisusedtoanalyzetherelationbetweentheprojections
28
ofY
f
andZ
p
ontotheorthogonalspaceofU
f
. CVAisprovedtobeoptimal
insomecasescomparedtoN4SIDandMOESP,andthestatisticaloptimality
wasprovidedbyLarimore([40]).
In the equation (2.42) the state sequence X
k
are usually represented
bythepastinputandoutput.
X
k
=
¯
L
p
Z
p
(2.49)
where
¯
L
p
is the observer extended controllability matrix. If H
fp
is defined
as
H
fp
= Γ
f
¯
L
p
(2.50)
then CVA can be interpreted as estimating H
fp
first and perform SVD on
ˆ
H
fp
Z
p
to get the estimate of extended observability matrix. From this per-
spective,theideaofCVAistoperformtheCCAanalysison(2.42)toextract
thesmallestanglesbetweenY
f
Π
?
U
f
andZ
p
Π
?
U
f
. Twoweightingmatricesare
usedbeforetheSVD
W
r
Y
f
Π
?
U
f
Z
T
p
W
c
=
[
U
1
U
2
]
[
S
1
0
0 S
2
][
V
T
1
V
T
2
]
where
W
r
=
(
Y
f
Π
?
U
f
Y
T
f
)
1=2
W
c
=
(
Z
p
Π
?
U
f
Z
T
p
)
1=2
ThebalancedrealizationforΓ
f
is
ˆ
Γ
f
=W
1
r
U
1
S
1=2
1
29
2.4.4 Unifying Theorem
van Overschee and de Moor ([66]) developed an unifying theorem
which provided insights into the relations among the representative open-
loop SIMs. According to the unifying theorem, MOESP, N4SID and CVA
can all be interpreted as oblique projection with different choices of the
weightingmatricesintherankreduction.
Theorem2.1. UnifyingTheoremBesidestheassumptionsinprevioussection,if
thefollowingconditionsaresatisfied,
1. Theinputu
k
ispersistentlyexcitingoftheorder2f.
2. Aninfinitenumberofmeasurementsareavailable.
3. W
1
isoffullrankandW
2
issuchthatrank(Z
p
) =rank(Z
p
W
2
).
Also,Θisdefinedastheobliqueprojection
Θ,Y
f
=
U
f
Z
p
withtheSVDasfollowing
W
1
ΘW
2
=
[
U
1
U
2
]
[
S
1
0
0 0
][
V
T
1
V
T
2
]
(2.51)
thenwehavethefollowing
1. Theorderofthesystemisequaltothenumberofnon-zerosingularvaluesin
(2.51).
30
2. TheextendedobservabilitymatrixΓ
f
canbetakenequalto
Γ
f
=W
1
1
U
1
S
1=2
1
3. TheKalmanstatesequencecanberecoveredfrom
X
k
= Γ
y
f
Θ
ForN4SID,theweightingmatricesare
W
1
= I
W
2
= I
ForMOESP,theweightingmatricesare
W
1
= I
W
2
= Π
?
U
f
andforCVA,
W
1
=
[
Y
f
Π
?
U
f
(
Y
f
Π
?
U
f
)
T
]
1=2
W
2
= Π
?
U
f
2.5 Problems in Subspace Identification
Inthissection,webrieflyreviewtheissuesofthesubspaceidentifica-
tion. AlthoughSIMsarewidelyrecognizedforitsadvantagesincomputing
efficiency and strength of handling multivariable systems, some problems
arealsowellrecognized.
31
1. Thefirststepofmostopen-loopSIMsisprojectingY
f
totheorthogonal
subspaceofU
f
. Thefutureinputtermisremovedinthiswaywithout
affectingthe noise. However, thisrequiresthefutureinput andnoise
are uncorrelated, which is only valid under open-loop condition. For
the closed-loop data, approaches need to be developed to handle the
correlationbetweeninputandnoise.
2. Another important assumption for SIMs is that the number of data
points is infinite. This condition cannot be satisfied in the practical
systemidentificationsinceinpracticethedatalengthisalwaysfinite.
InthiscasetheSteady-stateKalmanfiltermodelin(2.20)and(2.21)is
no longer valid. To address this, the non-steady Kalman filter should
beemployedandthecorrespondingsystemidentificationframework
needstobeestablished.
32
Chapter 3
Progressive Parameterization in SIMs with Finite
Horizons
3.1 Introduction
Subspace identification (SID) has experienced fast development
in the last twenty years in both theory and practice due to its simple
parametrization for MIMO systems and the non-iterative numerical
solutions. For safety reasons or quality restrictions, it is desirable that
identification experiments are carried out under the closed-loop or partial
closed-loop condition. As pointed out by many researchers [62, 44], the
fundamental problem with closed-loop data is the correlation between
the unmeasurable noise and the input. Because of this correlation early
formulations of SID algorithms cannot work on closed loop data. Ljung
and McKelvey [45] investigated the SID formulation through the classical
realization path and proposed a recursive prediction approach based on
an ARX model as a feasible approach to closed-loop SID. Formulated in
an errors-in-variables (EIV) framework, Chou and Verhaegen [8] proposed
a new method that can be applied to closed-loop data. The algorithm has
to treat white input from non-white input differently. Wang and Qin [73]
33
proposed the use of parity space and principal component analysis (PCA)
for EIV and closed-loop identification. Recent work of Qin and Ljung [54],
Jansson[29,30],ChiusoandPicci[7]analyzedsubspaceidentificationwith
feedback, proposed several new closed-loop SID methods and provided
consistencyanalysistothesemethods.
Peternell[52]investigatedthemainstepsoftheSIDforopenloopda-
ta. For the closed loop data, SID also involves several stages to realize the
subspace in terms of state space model parameters. These stages include
i)ahighorderARXsteptoestimateMarkovparametersaspre-estimation;
ii) an SVD step to extract a reduced rank expression of the Hankel matrix
toobtaintheobservabilitymatrixandcontrollabilitymatrix;andiii)tofur-
ther parametrize the observability matrix and controllability matrix by the
state space model parameters. In many of the recent SID formulations, the
steady-stateKalmanfilterframeworkisusedtointerprettheSIDstepsand
toanalyzeconsistency. Toachievethissteadystateconditionthefutureand
past horizons used in SID have to approach infinity. In practice, however,
these horizons must be finite due to finite data length and numerical cost.
Therefore,thesteadystateKalmanfilterframeworkdoesnotapply.
In this chapter, we introduce a progressive parametrization frame-
worktointerpretthemodelsusedineachstepofSIDmethodsanddiscuss
how the progressively parametrized models lead to the recursive state s-
34
pacemodels,whenadditionalassumptionsaremade. Wealsostatethatthe
intermediate non-recursive models can be useful for the purpose of state
estimation,faultdetection,andcontrol.
3.2 Progressive Parametrization of SID Models
3.2.1 High-order ARX Models
The task of identification is to estimate an appropriate dynamic
model for a properly collected time series of system input and output data
{u
k
;y
k
};k = 1;2;:::. The identified models will be used for additional
purposes, such as state estimation, output prediction, fault detection,
and control. The experimental data sequence {u
k
;y
k
} can be collected
from open-loop or closed-loop conditions. When causality is required for
the model, a generic choice for such a model is the auto-regressive with
exogenousinput(ARX)model,
y
k
=
Nu
∑
i=1
¯ g
i
u
ki
+
Ny
∑
i=1
¯
h
i
y
ki
+e
k
(3.1)
where e
k
is the model error term. This model indicates that the current
outputcanbepredictedfromafinitewindowofpastinputandoutputdata.
IntheARXmodelonlycausalityandlinearityareenforced. Iftheda-
ta sequences{u
k
;y
k
} are collected from low-order models with a recursive
state space model relation, the ARX model is not an efficient parametriza-
tion. ForthisreasonwerefertotheARXmodelasnonparametric.
35
Assuming the data sequence {u
k
;y
k
} is generated by the following
statespacemodel,
x
k+1
=Ax
k
+Bu
k
+w
k
y
k
=Cx
k
+v
k
(3.2)
wherex
k
∈R
n
. We demonstrate that this state equation can be adequately
represented by the ARX with a high but finite order. Assuming (C;A) is
observable,wecandesignastablestateobserverasfollows:
ˆ x
k+1
=
¯
Aˆ x
k
+Bu
k
+Ky
k
y
k
=Cˆ x
k
+e
k
(3.3)
where
¯
A =A−KC. Theabovestateobserversequencewillingenerallead
toaninfiniteorderARX(IOARX)model,
y
k
=
1
∑
i=1
¯ g
i
u
ki
+
1
∑
i=1
¯
h
i
y
ki
+e
k
(3.4)
with
¯ g
i
=C
¯
A
i1
B (3.5)
¯
h
i
=C
¯
A
i1
K: (3.6)
The question now is whether there exists K such that the IOARX model
becomes finite order. This is indeed possible if
¯
A
p
= 0, which corresponds
toadeadbeatobserver. TheIOARXmodel(3.4)becomesexactly
y
k
=
p
∑
i=1
¯ g
i
u
ki
+
p
∑
i=1
¯
h
i
y
ki
+e
k
(3.7)
36
Thiscoefficientsequences{¯ g
i
}and{
¯
h
i
}areknownasthepredictorMarkov
parameters(PMP)foru
k
andy
k
,respectively. However,itisdifficulttofind
fromdatathetime-invariant
¯
Athatisnilpotent. Eventhoughsuchamatrix
¯
Aexists,itbynomeansminimizesthemodelerrore
k
. Therefore,oneoften
assigns a high order p to find an approximate estimate. To estimate these
PMPs it is convenient to use least squares, which makes the model error
{e
k
}uncorrelatedtothedata{u
ki
}
i=1:p
and{y
ki
}
i=1:p
,thatis
R
eu
(i) =
1
N −1
N
∑
k=1
e
k
u
T
ki
= 0; i = 1;:::;p (3.8)
R
ey
(i) =
1
N −1
N
∑
k=1
e
k
y
T
ki
= 0; i = 1;:::;p (3.9)
3.2.2 Subspace Identification Models with a Prediction Horizon
Assume now we are interested in a model that can predict a fu-
ture horizon of output y
k
based on an observed window of past data
{u
ki
;y
ki
}
i=1:p
. Defining
y
f
(k) =





y
k
y
k+1
.
.
.
y
k+f1





;u
p
(k) =





u
k1
u
k2
.
.
.
u
kp





;y
p
(k) =





y
k1
y
k2
.
.
.
y
kp





and
z
p
(k) =
[
u
p
(k)
y
p
(k)
]
37
a non-parametric, generic model to represent the window of future output
willbe
y
f
(k) =
¯
H
fp
z
p
(k)+
¯
G
f
u
f
(k)+
¯
H
f
y
f
(k)+e
f
(k): (3.10)
Toenforcecausalitythe(i;j)
th
blockof
¯
G
f
and
¯
H
f
mustbe0fori≤j. Infact
we can parametrize these two matrices using the estimates of the HOARX
modelasfollows:
¯
G
f
=







0 0 0 ··· 0
¯ g
1
0 0 ··· 0
¯ g
2
¯ g
1
0 ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
¯ g
f1
¯ g
f2
··· ¯ g
1
0







(3.11)
¯
H
f
=







0 0 0 ··· 0
¯
h
1
0 0 ··· 0
¯
h
2
¯
h
1
0 ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
¯
h
f1
¯
h
f2
···
¯
h
1
0







(3.12)
Defining
˜ y
f
(k) =y
f
(k)−
¯
G
f
u
f
(k)−
¯
H
f
y
f
(k)
=
¯
H
fp
z
p
(k)+e
f
(k) (3.13)
which is the prediction of the future window from the data of the past
window alone, the non-parametric matrix
¯
H
fp
can be estimated using least
squares similar to the HOARX estimates. The resulting non-parametric
model can predict the future window of output based on the past window
ofdatawithoutspecifyingthesystemorobserverorder.
38
3.2.3 Subspace Identification Models with Specified Order
If we assume the system order isn (n < f;p), then the past window
of system input/output data can be memorized in a state vector ˆ x
k
∈ R
n
,
whichcanthenbeusedtopredictitseffectonthefuturewindowofoutput.
Torepresenttheabovestatementmathematically,wehave
ˆ x
k
=
¯
L
p
z
p
(k) (3.14)
˜ y
f
(k) =
¯
Γ
f
ˆ x
k
+e
f
(k) =
¯
Γ
f
¯
L
p
z
p
(k)+e
f
(k): (3.15)
Usingsimilarargumentthatthestatespacemodelcanbeexactlyrepresent-
edbyafinitebuthighorderARXmodelif
¯
A
p
= 0,Equation(3.14)isanexact
representationofthestateobserver(6.2). However,anilpotent
¯
Acannotbe
easily found and it would by no means give a reasonably good prediction
oftheoutputwhennoiseexists.
Tointerpretthemodelof(3.14)and(3.15)withfinitef andp,were-
sort to the time-varying Kalman filter, which has a time varying observer
gain in (6.2). In this case,
¯
A will be time-varying. However, (3.14) is still a
valid state prediction with
¯
L
p
of rank no more than n, and (3.15) is still a
valid output prediction with
¯
Γ
f
of rank no more than n. For convenience,
we refer to (3.14) as the batch-wise state observer with finite data window
and(3.15)thecorrespondingbatch-wiseoutputprediction. Itisworthnot-
ing that this model, albeit non-recursive, is useful for fault detection, state
estimation,andcontroloffinitehorizons(e.g.,modelpredictivecontrol).
39
Comparing (3.15) to (3.13), it is clear
¯
H
fp
=
¯
Γ
f
¯
L
p
except that
¯
Γ
f
¯
L
p
has a reduced rank of n while
¯
H
fp
in (3.13) can be full rank. To estimate
¯
Γ
f
and
¯
L
p
of (3.14) and (3.15), canonical correlation analysis (CCA) gives
an optimal estimate of
¯
Γ
f
and
¯
L
p
. This approach has been suggested by
Larimore [39] as canonical variate analysis (CVA) so we inherit the same
terminology. Defining
R
fp
=
1
N
N
∑
k=1
˜ y
f
(k)z
T
p
(k) (3.16)
R
ff
=
1
N
N
∑
k=1
˜ y
f
(k)˜ y
T
f
(k) (3.17)
R
pp
=
1
N
N
∑
k=1
z
p
(k)z
T
p
(k) (3.18)
andperformingsingularvaluedecomposition(SVD)on
R

1
2
ff
R
fp
R

1
2
pp
=U
n
D
n
V
T
n
+Σ
n
(3.19)
whereU
n
representsthefirstnsingularvectors,theCVAestimatesare
¯
Γ
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
:
(3.20)
Themodelcoefficientmatrix
¯
H
CVA
fp
=
¯
Γ
CVA
f
¯
L
CVA
p
(3.21)
hasnowreducedorderrankofn. Thisrankreductionisnecessarytoleadto
astatevectoroflengthn. WerefertothisSIDmethodastheSSARXmethod
40
withCVAweightingduetoJansson[30]. WithouttheCVAweightingaleast
squaressolutiontotheHankelmatrixcanbefound. PerformingSVDonthis
Hankel matrix leads to an SID method to be referred to as SSARX method
withLSweighting.
Remark 1. The model (3.14) and (3.15) has finite state length n, but
isstillanon-parametricmodel. Thismodelestimatesthecurrentstatefrom
past window of data in one batch (i.e., non-recursive) and calculate the fu-
ture window of output in one bach (i.e., non-recursive). This model is in-
deedusefulinestimatingthecurrentstateorpredictingawindowoffuture
output.
Remark2. If the CVA calculation does not truncate, then the result-
ingestimatein(3.21)willapproachtheleastsquares,
¯
H
LS
fp
=R
fp
R
1
pp
.
3.2.4 Subspace Identification Models with Recursive State Space
Parametrization
If the data sequence {u
k
;y
k
} comes from a truly recursive state s-
pace representation as in (6.2), then the most efficient parametrization of
the model is the recursive representation of (6.2), which means we need to
furtherparametrizethemodelintermsof
¯
A;B;C andK. AsanoptimalK
intheKalmanfiltersensewouldnotbetime-invariantifthehorizonsarefi-
nite,thisfurtherparametrizationwithatime-invariantK isanapproximate
representation,andtherefore,itissuboptimal. However,ifpischosenlong
41
enough,suchthat
¯
A
p
≈ 0,theapproximationisadequate.
There are many ways to estimate the state space parameters from
the
¯
Γ
f
and
¯
L
p
estimates. We discuss the following two methods. For more
details,seeQin[53].
1. Use
¯
L
p
toreconstructthestatesequence{ˆ x
k
},andestimate{
¯
A;B;C;K}
from (6.2) by least squares. We refer to this method as the state se-
quence(SS)method.
2. EstimateC and
¯
Afrom
¯
Γ
f
,thenestimateB andK usinglinearregres-
sionbasedontheestimateofC and
¯
A[53].
3.3 Comparison of SID Models at Various Stages
ItisclearfromtheaboveanalysisthatacompleteSIDproceduregoes
throughseveralstages.
• Stage1: PerformhighorderARXtogiveMarkovparameterestimates,
whichisanon-parametricandhasitsownapplicationinpractice(for
example, in dynamic matrix control where output-error formulation
isused).
• Stage 2: Perform SVD on the Hankel matrix or weighted Hankel ma-
trix,leadingto
¯
Γ
f
and
¯
L
p
estimates,whichbythemselvesarereduced
order(n)non-parametricmodelstoestimatethestateandpredictthe
output.
42
• Stage 3: Performleastsquarestoestimate(
¯
A;B;C;K)from
¯
Γ
f
and
¯
L
p
leadingtostate-spaceparametricmodels.
Sinceateachoftheabovestagessomeformofmodelsareobtained,
a legitimate question is, what does each stage contribute to the estimation
of the system model? To answer this question, we compare the predictor
Markovparametersestimatedaftereachstageforthefollowingreasons:
1. TheMarkovparametersareinvarianttothechoiceofstatecoordinate,
and therefore, comparable to the real Markov parameters of the sys-
tem;
2. The predictor Markov parameters allow us to examine whether the
subspace identification improves on these estimates from stage to
stage,astheyareestimatedateachstage;and
3. The predictor Markov parameters can be easily converted to system
Markov parameters, which can be further converted to step response
modelspopularinmodelpredictivecontrolpractice[56].
Defining
g
i
=CA
i1
B (3.22)
h
i
=CA
i1
K; (3.23)
43
these system Markov parameters can be obtained from the predictor
Markovparametersasfollows[35]:
g
i
= ¯ g
i
+
i1
∑
j=1
¯
h
j
g
ij
(3.24)
h
i
=
¯
h
i
+
i1
∑
j=1
¯
h
j
h
ij
(3.25)
To summarize the various SID methods and the resulting models
aftereachstage,wegivetheflowchartinFig. 3.1.
The predictor Markov parameters after each stage have increased
restrictions:
• Stage1: high-order,non-parametricmodel;
• Stage 2: reduced-order, non-parametric model, leading to batch-wise
stateestimationandoutputprediction;and
• Stage3: reduced-order,parametricstatespacemodel,leadingtorecur-
sivestateestimationandoutputprediction.
3.4 Numerical Comparison of SID Models via Monte-Carlo
Simulation
TorevealthecontributionofeachstageofSIDproceduretothemod-
el estimates, we use the following closed-loop example adopted from Ver-
44
haegen[68].
A =






4:40 1 0 0 0
−8:09 0 1 0 0
7:83 0 0 1 0
−4:00 0 0 0 1
0:86 0 0 0 0






;B =






:00098
0:01299
0:01859
0:0033
−0:00002






;
C
T
=






1
0
0
0
0






;K =






2:3
−6:64
7:515
−4:0146
0:86336






;D = 0
The controller used has the state space model composed by following ma-
trices:
A
c
=




2:65 −3:11 1:75 −0:39
1 0 0 0
0 1 0 0
0 0 1 0




;B
c
=




1
0
0
0




;
C
c
=
[
−0:4135 0:8629 −0:7625 0:2521
]
;D
c
= 0:61:
For each method we generate 10 data sets using Monte-Carlo simulation,
eachwith4000datapointsusingtheabovesimulationmodelinMATLAB.
Thesetpointofthereferencesignalis1. TheinnovationsequenceeisGaus-
sianwhitenoisewithvariance1=9. WeaddaPRBSsignaltoinputchannel
uinordertoimproveclosed-loopsystemidentifiability. Ineachsimulation,
pandf aresetas20whicharethepastandfuturehorizon,respectively.
Theinitialorderischosentobef+p−1 = 39tohaveenoughMarkov
parametersformingtheHankelmatrix
¯
H
ARX
fp
.
45
The predictor Markov parameters after stages 1 and 2 are shown in
Fig. 3.2 with three different Hankel matrices to perform SVD, i) HOARX
estimates,ii)SSARXwithleastsquaresestimates,andiii)SSARXwithCVA
weighting. Theeigenvaluesatstage2isshowninTable3.1. FromFig. 3.2it
isobservedthat
1. weighting(byCVA)iscritical,
2. thestage2achievesdead-beatobserverswherethepredictorMarkov
parameters from SSARX tend to zero after the horizonp. Comparing
totheHOARXestimatesofthePMPs,biasisclearlyobservedinstage
2.
Toexaminetheimpactofstage3ontheMarkovparameterestimates,
we show in Fig. 3.3 the results using SSARX with CVA weighting and S-
SARXwithoutweighting.
TheHankelmatricesinFig. 3.2afterSVDdonotexactlyconformto
aHankelstructure,althoughtherankisrestrictedton. ThefirstpPMPsin
Fig. 3.2areextractedfromthefirstrowoftheHankelmatrix,whilethelast
(f−1)PMPsareextractedfromthelastcolumnoftheHankelmatrix. While
thesePMPsfromSSARXwithorwithoutCVAweightingdonotshowsig-
nificant difference, the eigenvalues shown in Table 3.1 are very different.
The HOARX or SSARX without CVA weighting have only two significant
46
eigenvalues, while the SSARX with CVA weighting has clearly five signif-
icant eigenvalues. It is concluded that the results with CVA weighing is
consistent with the real system that has five state variables. To further ex-
ploit the difference between SSARX with and without CVA weighting, we
show in Fig. 3.4 another set of PMPs from the Hankel matrices after Stage
2. The firstf PMPs in this figure are extracted from the first column of the
Hankel matrices, while the last (p− 1) PMPs are from the last row of the
Hankel matrices. It is clear that the SSARX with CVA weighting has much
betterestimatesatthisstage.
3.5 Summary
For subspace identification with finite horizons, which is the only
choiceinpractice,weintroduceaprogressiveparameterizationframework
for the models used at various stages. It is stated that the models from
the SVD stage of SID are batch-wise state observers and output predictors,
whichbythemselvesareusefulforestimation,control,andfaultdetection.
To further parameterize the batch-wise models, a sub-optimal recursive s-
tateobserverisintroduced, whichisanapproximationofthetime-varying
Kalmanfilter. Forboththebatch-wiseandrecursivestatespaceparameter-
izations, CVA weighting is shown to be critical in the SVD step of the SID
procedure.
47
Table3.1: EigenvaluesofHankelMatrixorWeightedHankelMatrix
OKID SSARX SSARXwith
(noweighting) CVAweighting
ratioofthelargest
94.15% 94.89% 93.69%
5eigenvalues
1steigenvalue 1.0000 12.1444 12.2807
2ndeigenvalue 1.1997 1.1629 0.9999
3rdeigenvalue 0.6502 0.6318 0.9998
4theigenvalue 0.4611 0.4333 0.9903
5theigenvalue 0.4374 0.3978 0.9857
6theigenvalue 0.3155 0.2482 0.0443
7theigenvalue 0.1857 0.1535 0.0419
8theigenvalue 0.1094 0.0951 0.0376
9theigenvalue 0.0759 0.0869 0.0305
10theigenvalue 0.0714 0.0754 0.0302
48
Data
{ , }
k k
u y
HOARX
, f f G H
Least Squares
LS
1
fp
fp pp
H R R

!
CVA weighting
LS
1
1 2
fp W H W

ARX
fp H
SVD , f p L " Markov parameters
ˆ ( ) p
k p
x L z k ! f "
, , , A B C K Markov parameters
Markov parameters
Figure3.1: Variousstagesofsubspaceidentificationleadingtotheestimates
ofMarkovparameters.
49
0 10 20 30 40
0
0.05
0.1
0.15
0 10 20 30 40
-2
-1
0
1
2
3
a) b)
0 10 20 30 40
0
0.05
0.1
0.15
c)
0 10 20 30 40
-2
-1
0
1
2
3
d)
0 10 20 30 40
0
0.05
0.1
0.15
e)
0 10 20 30 40
-2
-1
0
1
2
3
f)
0 10 20 30 40
0
0.05
0.1
0.15
g)
0 10 20 30 40
-2
-1
0
1
2
3
h)
Figure 3.2: Estimated PMPs after stage 1 and 2. Red curves are true PMPs
and dots are estimated PMPs. a) estimated PMPs of
¯
G from HOARX after
stage 1; b) estimated PMPs of
¯
H from HOARX after stage 1; c) estimated
PMPsof
¯
GfromHOARXafterstage2(reducedrank);d)estimatedPMPsof
¯
H fromHOARXafterstage2(reducedrank); e)estimatedPMPsof
¯
Gfrom
¯
H
LS
fp
;f)estimatedPMPsof
¯
H from
¯
H
LS
fp
;g)estimatedPMPsof
¯
Gfrom
¯
H
CVA
fp
;
h)estimatedPMPsof
¯
H from
¯
H
CVA
fp
.
50
0 10 20 30 40
0
0.05
0.1
0.15
0 10 20 30 40
0
0.05
0.1
0.15
a)
b)
Figure 3.3: Estimated PMPs from the Hankel matrices after stage 2. The
firstf PMPsareextractedfromthefirstcolumnoftheHankelmatrices;the
last(p−1)PMPsareextractedfromthelastrowoftheHankelmatrices. a)
estimatedPMPsof
¯
GwithoutCVAweighting;b)estimatedPMPsof
¯
Gwith
CVAweighting.
0 10 20 30 40
-0.05
0
0.05
0.1
0.15
0 10 20 30 40
-0.05
0
0.05
0.1
0.15
0 10 20 30 40
-0.05
0
0.05
0.1
0.15
0 10 20 30 40
-0.05
0
0.05
0.1
0.15
a) b)
d)
c)
Figure3.4: EstimatedPMPsof
¯
Gafterstage3. a)estimatedPMPsof
¯
Gusing
SSARX without weighting, system matrices estimated using
¯
Γ
f
and
¯
L
p
; b)
estimated PMPs of
¯
G using SSARX with CVA weighting, system matrices
estimated using
¯
Γ
f
and
¯
L
p
; c) estimated PMPs of
¯
G using SSARX without
weighting, system matrices estimated using state sequence; d) estimated
PMPs of
¯
G using SSARX with CVA weighting, system matrices estimated
usingstatesequence.
51
Chapter 4
NKF Parameterization based Disturbance Model-
ing and MPC Application: Part I
4.1 Introduction
Subspaceidentificationmethods(SIMs)havegainedtremendousde-
velopment during the last few decades. The traditional subspace identifi-
cation methods such as N4SID [65], CVA [39] and MOESP [70] are biased
under closed-loop conditions due to the correlation between the input and
thenoise. LjungandMcKelvey[45]proposedtheutilizationofahighorder
ARX (HOARX) model to address this problem. Jansson [29] estimated the
predictor Markov parameters from the HOARX model. These Markov pa-
rameterswereusedtoenforcemorestructureintothedatamodel. Qinand
Ljung [55] performed least squares (LS) over f HOARX models to remove
the non-causal terms. The consistency of those algorithms has been inves-
tigated in [7]. In [53] and [57], the main steps of the closed-loop subspace
identification were investigated. The first step of most closed-loop SIMs is
the pre-estimation [29, 41, 61]. In this step, the predictor Markov param-
eter estimates are obtained from the HOARX model through least squares
estimation.
52
Although the LS yields the unstructured Markov parameter esti-
mates, most subspace identification and fault detection methods interpret
the Markov parameter estimates with the steady state Kalman filter (SKF).
However, by showing the similarity between the objective functions
of the LS estimation and the non-stationary Kalman filtering, Sorenson
[63] showed that the non-stationary Kalman filter (NKF) represented a
recursive solution to the LS problems. Under the subspace framework,
Van Overschee and De Moor [65] proved that least squares regression of
the future output on the past input output data actually represents a bank
of NKF. Huang, Ding and Qin [28] also interpreted the future state as the
NKF state. However, those results were obtained under the open-loop
condition. This motivates us to parameterize HOARX using NKF, which
hasadvantagesindealingwithpracticalfinitedatahorizons.
Model predictive control (MPC) is one of the most widely used ad-
vancedprocesscontroltechniqueinindustry. Dynamicmatrixcontrol(DM-
C)[10],asthepioneerMPCalgorithmstillplaysanimportantroleinmany
areas, as reported by Qin and Badgwell [56]. As a model-based control
method, DMC control performance is affected by model accuracy greatly.
The model includes both plant and disturbance model. Compared with
plant model, disturbance model, however, receives little attention in prac-
tice. Only a simple step-like disturbance model is used in DMC to achieve
53
offset-freecontrol. Tobettermodelstateandstatedisturbances,Muskeand
Badgwell[49]proposedablockdiagonaldisturbancemodel. Itisextended
to a more general disturbance model without special structure by Pannoc-
chiaandRawlings[51]. Themethodrequiresapredeterminedfiltergainof
thestateestimator.
In this chapter we propose an NKF parameterization of the predic-
tor Markov parameters (PMPs) which are estimated by least squares. It is
shownthattheLSestimationonthefinitehorizonisequivalenttoabankof
NKF. The transformation from PMPs to system Markov parameters (SMP-
s) is developed under NKF framework. Based on the recursive relation,
the SMPs are obtained and further implemented to the disturbance model
inDMC.Comparedtothetraditionalstep-likedisturbancemodel,thepro-
posed SKF and NKF parameterized disturbance models show advantages
in improving DMC control performance. In addition, NKF parameteriza-
tionoutperformstheSKFunderfinitedatalengthcondition.
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
∈ R
ny
and u
t
∈ R
nu
as the output and input
vector of a MIMO plant respectively. The relationship inputs and outputs
54
canbewrittenas[18]
y(t) =
N
∑
i=1
S
i
∆u(t−i)+y
0
+d(t) (4.1)
where y
0
is the output initial condition, d(t) denotes unmodelled distur-
bances∆u(j) = u(j)−u(j−1),N isthenumberofstepsforplanttoreach
steady-state,andS
i
istheithstepresponsematrix
S
i
=





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





(4.2)
inwhichs
j;k;i
representstheithstepresponsecoefficientfromthekthinput
tothejthoutput. Thek-step-aheadpredictioncanbemadeby
ˆ y(t+k|t) =y(t)+
k
∑
i=1
S
i
∆u(t+k−i)+
ˆ
d(t+k|t) (4.3)
where ˆ y(t+k|t)and
ˆ
d(t+k|t)denotesthepredictionofy(t+k)andunmod-
elleddisturbancesgiventheinformationuptotimet,respectively,
The step response model and finite impulse response (FIR) model
canbeusedinterchangeably:





S
1
S
2
.
.
.
S
P





=





1
1 1
.
.
.
.
.
.
.
.
.
1 ··· 1










G
p
1
G
p
2
.
.
.
G
p
P





(4.4)
whereG
p
1
;:::;G
p
P
aretheimpulseresponsecoefficients, alsoknownassys-
temMarkovparameters,oftheplantmodel.
55
4.2.2 DMC predictions and offset-free control
From(4.3),thepredictedoutputscanberewrittenas
ˆ y(t+k|t) =
k
∑
i=1
S
i
∆u(t+k−i)+
P
∑
i=k+1
S
i
∆u(t+k−i)
+y
0
+
ˆ
d(t+k|t) (4.5)
whereP indicatesthepredictionhorizon. Define
y

(t+k) =y
0
+
P
∑
i=k+1
S
i
∆u(t+k−i)
which is the effect of past inputs, and stack ˆ y(t +i|t);(i = 1;:::;P) using
(4.5):



ˆ y(t+1|t)
.
.
.
ˆ y(t+P|t)



=



y

(t+1)
.
.
.
y

(t+P)



+A
M
P



∆u(t)
.
.
.
∆u(t+M −1)



+



ˆ
d(t+1)
.
.
.
ˆ
d(t+P)



(4.6)
whereM isthecontrolhorizonand
A
M
P
=














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














(4.7)
iscalledthedynamicmatrixofthesystem.
DMC assumes a step disturbance model and estimate the step size
using the difference between the current measured and estimated output,
56
i.e.,
ˆ
d(t+i|t) =y(t)− ˆ y(t|t−1); i = 1;:::;P (4.8)
Thisdisturbancemodelisintegratingandisabletoachieveoffset-freecon-
trolforDMC.
4.2.3 Plant model in state-space form
We use the state-space form of plant model which will be useful in
thederivationslater. Assumetheplantisdescribedby
x
p
(t+1) =A
p
x
p
(t)+B
p
u(t)
y
p
(t) =C
p
x
p
(t)
(4.9)
whichisconsideredtobedeterministicandthedisturbancemodelis
writteninthefollowinginnovationformofKalmanfilter:
ˆ 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)
whichisconsideredtobestochastic. Whenthestep-likedisturbancemodel
in DMC is assumed, A
d
= I
ny
, K
d
= I
ny
, and C
d
= I
ny
; this type of distur-
bance model corresponds to the case that disturbance is integrated white
noise.
Byaugmenting(4.9)and(4.10),weobtain
{
ˆ 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)
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)
57
Thesubscriptaugreferstoaugmentedsystem.
Clearly, the augmented Kalman gain with the structure K
aug
=
(
0 K
d
)
T
is insufficient when there exists plant model mismatch or distur-
bancemodelmismatch. Withthehelpofclosedloopdataandthefirststage
of subspace identification, it is possible to estimate a general 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
toestimatedisturbancemodelfromclosed-loopdataandarevisedformof
DMCthatincludesdisturbancemodel.
4.3 Equivalence Between LS Estimation and NKF
Beforedevelopingthedatabaseddisturbancemodelingscheme,we
need to prove the suitability of interpreting data with finite length of data
withNKF.Assumethedatasequence{u
k
;y
k
}isgeneratedbythestatespace
modelofthefollowingprocessform,
x
k+1
=Ax
k
+Bu
k
+w
k
y
k
=Cx
k
+v
k
whereu
k
∈R
nu
,y
k
∈R
ny
,x
k
∈R
n
,w
k
∈R
n
,v
k
∈R
ny
standsforthesystem
input,output,state,processnoise,andmeasurementnoise,respectively. A
Kalmanfilter(KF)canbedesignedtoestimatethestate,sothatthefollow-
inginnovationformcanbederived
x
k+1
=Ax
k
+Bu
k
+Ke
k
(4.13)
58
y
k
=Cx
k
+e
k
(4.14)
where the innovation e
k
is white noise which is independent of the past
input output and K is the steady-state Kalman filter. We use x
k
instead of
ˆ x
k
for convenience. Through the iterative substitution of (4.13) and (4.14),
thefollowingsubspaceequationscanbederived[65],
Y
p
= Γ
i
X
p
+H
d
i
U
p
+H
s
i
E
p
(4.15)
Y
f
= Γ
i
X
f
+H
d
i
U
f
+H
s
i
E
f
(4.16)
X
f
=A
i
X
p
+∆
d
i
U
p
+∆
s
i
E
p
(4.17)
wherethesubscriptpandf standforthe“past”and“future”,respectively.
Thepastandfutureinputblock-Hankelmatricesaredefinedasfollowing,
U
p
=U
0ji1
=





u
0
u
1
··· u
j1
u
1
u
2
··· u
j
.
.
.
.
.
.
.
.
.
.
.
.
u
i1
u
i
··· u
i+j2





U
f
=U
ij2i1
=





u
i
u
i+1
··· u
i+j1
u
i+1
u
i+2
··· u
i+j
.
.
.
.
.
.
.
.
.
.
.
.
u
2i1
u
2i
··· u
2i+j2





Theoutputandinnovationblock-HankelmatricesY
p
,Y
f
,E
p
andE
f
areformedsimilartoU
p
andU
f
. Wealsodenote
Z
p
=
[
U
p
Y
p
]
59
Thestatesaredefinedas
X
p
=X
0
=
[
x
0
x
1
··· x
j1
]
(4.18)
X
f
=X
i
=
[
x
i
x
i+1
··· x
i+j1
]
(4.19)
Theextendedobservermatrixisdefinedas
Γ
i
=







C
CA
CA
2
.
.
.
CA
i1







Thereversedextendedcontrollabilitymatrices∆
d
i
and∆
s
i
aregivenas
∆
d
i
=
[
A
i1
B A
i2
B ··· AB B
]
(4.20)
∆
s
i
=
[
A
i1
K A
i2
K ··· AK K
]
(4.21)
Theblock-ToeplitzmatricesH
d
i
andH
s
i
havethefollowingstructure
H
d
i
=





0 0 ··· 0
CB 0 ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
CA
i2
B CA
i3
B ··· 0





(4.22)
H
s
i
=





I 0 ··· 0
CK I ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
CA
i2
K CA
i3
K ··· I





(4.23)
Itisworthnotingthatbyimposing∆
d
i
,∆
s
i
,H
d
i
andH
s
i
withthestruc-
tureas(4.20)-(4.23),weparameterizethesubspaceequationsand(4.14)with
a steady-state Kalman gain K. However, by inspecting (4.15)-(4.17) more
60
closer,wefindthat(4.15)-(4.17)doesnotrequireasteadystateKalmangain
in(4.13). Withoutspecifyingthesub-blocksof∆
d
i
,∆
s
i
,H
d
i
andH
s
i
,weshow
that for finite data set,X
f
is the state estimation of a bank of NKF withX
p
as the initial state and Z
p
as the input/output. Therefore, parameterizing
∆
d
i
,∆
s
i
,H
d
i
andH
s
i
withNKFistheaccurateinterpretationof(4.15)-(4.17).
From(4.16),itcanbeobtainedthat
E
p
=H
s
i
Y
p
−H
s
i
Γ
i
X
p
−H
s
i
H
d
i
U
p
(4.24)
Bysubstituting(4.24)into(4.17),wehave
X
f
=
(
A
i
−∆
s
i
H
s
i
Γ
i
)
X
p
+
[
∆
d
i
−∆
s
i
H
s
i
H
d
i
| ∆
s
i
H
s
i
]
Z
p
(4.25)
Obviously, X
f
is state solution to a bank of non-stationary Kalman filter
withX
p
astheinitialstatevalue[65,21]. Thismotivatesustoparameterize
∆
d
i
,∆
s
i
,H
d
i
andH
s
i
withNKFratherthanthesteadyKalmangainK.
4.4 Identification of NKF Disturbance Models
4.4.1 Identification of PMPs
TheNKFparameterizationleadstoarepresentationofthestatespace
modelinthenon-steady-stateinnovationformasfollowing
x
k+1
=Ax
k
+Bu
k
+K
k
e
k
y
k
=Cx
k
+e
k
61
orequivalently,
x
k+1
=
¯
A
k
x
k
+Bu
k
+K
k
y
k
(4.26)
y
k
=Cx
k
+e
k
(4.27)
where
¯
A
k
=A−K
k
C.
Focusingonthefirstn
y
rowsof(4.16),wehave
Y
i
=CX
f
+E
i
(4.28)
where
Y
i
=Y
iji
=
[
y
i
y
i+1
··· y
i+j1
]
(4.29)
andE
i
hasthesimilarstructureasY
i
. Substituting(4.25)into(4.28)yields
Y
i
=C
(
A
i
−∆
s
i
H
s
i
Γ
i
)
X
p
+C
[
∆
d
i
−∆
s
i
H
s
i
H
d
i
| ∆
s
i
H
s
i
]
Z
p
+E
i
(4.30)
Ontheotherhand,theiterativesubstitutionof(4.26)and(4.27)yields
Y
i
=C
i1
∏
k=0
¯
A
k
X
p
+
[
¯
G |
¯
H
]
Z
p
+E
i
(4.31)
where
¯
G =
[
¯ g
NKF
i
¯ g
NKF
i1
··· ¯ g
NKF
1
]
¯
H =
[
¯
h
NKF
i
¯
h
NKF
i1
···
¯
h
NKF
1
]
(4.32)
and
¯ g
NKF
k
=
{
CB; k = 1
C
¯
A
i1
¯
A
i2
···
¯
A
ik+1
B; k = 2;:::;i
¯
h
NKF
k
=
{
CK
i1
; k = 1
C
¯
A
i1
¯
A
i2
···
¯
A
ik+1
K
ik
; k = 2;:::;i
62
where
{
¯ g
NKF
k
}
,
{
¯
h
NKF
k
}
are the PMPs. (·)
NKF
indicates (·) is parameterized
usingnon-steady-stateKalmangainsequence.
Acomparisonbetween(4.30)and(4.31)indicatesthat
A
i
−∆
s
i
H
s
i
Γ
i
=
i1
∏
k=0
¯
A
k
(4.33)
and
C
[
∆
d
i
−∆
s
i
H
s
i
H
d
i
| ∆
s
i
H
s
i
]
=
[
¯
G |
¯
H
]
(4.34)
Obviously, when i is large enough, the effect of the initial state can be re-
movedsincewehave
i1
∏
k=0
¯
A
k
≃ 0 (4.35)
and(4.30)becomes
Y
i
=C
[
∆
d
i
−∆
s
i
H
s
i
H
d
i
| ∆
s
i
H
s
i
]
Z
p
+E
i
(4.36)
whichisaHOARXmodel.
Weperformtheorthogonalprojectionof(4.25)ontotherowspaceof
Z
p
,
Y
i
=Z
p
=C
(
A
i
−∆
s
i
H
s
i
Γ
i
)
X
p
=Z
p
+C
[
∆
d
i
−∆
s
i
H
s
i
H
d
i
|∆
s
i
H
s
i
]
Z
p
=Z
p
+E
i
=Z
p
(4.37)
Since the future noise is independent of the past input/output, the noise
63
termof(4.37)canberemovedand(4.37)issimplifiedas
Y
i
=Z
p
=C
(
A
i
−∆
s
i
H
s
i
Γ
i
)
X
p
=Z
p
+C
[
∆
d
i
−∆
s
i
H
s
i
H
d
i
| ∆
s
i
H
s
i
]
Z
p
(4.38)
The left hand side of (4.38) is the LS estimation under the HOARX struc-
ture. Therighthandsideof(4.38)istheoutputofanon-stationaryKalman
filtering process with the initial state X
p
=Z
p
. From the result of (4.35), the
effect of the initial state can be eliminated with a sufficiently large i, (4.38)
thusbecomes
Y
i
=Z
p
=C
[
∆
d
i
−∆
s
i
H
s
i
H
d
i
| ∆
s
i
H
s
i
]
Z
p
=
[
¯
G |
¯
H
]
Z
p
(4.39)
Weuse(4.39)toestimatethePMPswithNKFparameterization.
Remark4.1. Equation(4.39)holdsforbothopen-loopandclosed-loopdata.
In [65], an oblique projection of (4.16) along the row space of U
f
onto the
rowspaceofZ
p
isperformed,
Y
f
=
U
f
Z
p
= Γ
i
X
f
=
U
f
Z
p
+H
d
i
U
f
=
U
f
Z
p
+H
s
i
E
f
=
U
f
Z
p
(4.40)
Under the open-loop condition, the noise term can be removed since
E
f
=
U
f
Z
p
= 0. However, E
f
=
U
f
Z
p
is not zero under the closed-loop condi-
tionduetothefeedback,thustheN4SIDproposedin[65]isnotapplicable
to the closed-loop data. The advantage of the HOARX model (4.36) is that
64
we do not need to deal with the future input term U
f
, and the noise term
canbeeasilyremovedbytheorthogonalprojection.
Remark 4.2. When the PMPs are parameterized with steady state Kalman
gain,wehave
¯
G =
[
¯ g
i
¯ g
i1
··· ¯ g
1
]
(4.41)
¯
H =
[
¯
h
i
¯
h
i1
···
¯
h
1
]
(4.42)
where ¯ g
k
= C
¯
A
k1
B,
¯
h
k
= C
¯
A
k1
K, and
¯
A = A−KC. In this case, (4.39)
still holds and it is often used for system identification and fault detection.
As the first step of most closed-loop subspace identification methods, the
estimated PMPs of the form (4.41) and (4.42) are further used to construct
the lower triangular Toeplitz matrices in order to remove the effect of the
futureinputandoutput[29].
4.4.2 Identification of Disturbance Models
In this section, we convert the estimated predictor Markov param-
eters into the system Markov parameters for the disturbance modeling in
MPC. Regardlessof whether the steady-state Kalman gain or time-varying
Kalmangainisusedforparameterization,thepredictorMarkovparameter
estimatescanbeobtainedfrom(4.39)bytheLS
[
ˆ
¯
G
ˆ
¯
H
]
=Y
i
Z
T
p
[
Z
p
Z
T
p
]
1
(4.43)
65
The traditional parameterization using SKF interprets the system
Markovparametersequenceasfollowing
G =
[
g
i
g
i1
··· g
1
]
H =
[
h
i
h
i1
··· h
1
]
whereg
k
=CA
k1
B andh
k
=CA
k1
K. ThevaluesofthesesystemMarkov
parameters can be computed recursively by the following conversion rela-
tion[34]
g
k
= ¯ g
k
+
k1
∑
m=1
¯
h
m
g
km
; k = 1;:::;i: (4.44)
h
k
=
¯
h
k
+
k1
∑
m=1
¯
h
m
h
km
; k = 1;:::;i: (4.45)
Wenowderivetherecursiverelationbetweenthepredictorandsys-
temMarkovparameterswhentheyarebasedonNKF.Inthiscase,wehave
G =
[
g
NKF
i
g
NKF
i1
··· g
NKF
1
]
H =
[
h
NKF
i
h
NKF
i1
··· h
NKF
1
]
whereg
NKF
k
=CA
k1
B andh
NKF
k
= CA
k1
K
ik
. Similar to [34], the conver-
sion relation between the predictor and system Markov parameters based
onNKFisgiveninarecursiveformasfollows
g
NKF
k
= ¯ g
NKF
k
+
k1
∑
m=1
¯
h
NKF
m
g
NKF
km
; k = 1;:::;i: (4.46)
h
NKF
k
=
¯
h
NKF
k
+
[
¯
h
NKF
k1
¯
h
NKF
k2
···
¯
h
NKF
1
]
Γ
k1
K
ik
;
k = 1;:::;i: (4.47)
66
Remark4.3. Itisworthnotingthat(4.46)hasexactlythesameformas(4.44),
meaningthattheconversionrelationoftheprocesssystemMarkovparam-
eters still holds although the Markov parameters are now parameterized
with NKF. On the other hand, (4.47) is a more general form of (4.45) since
(4.47) can be reduced to (4.45) when all the time-varying Kalman gain is
replacedbysteadyKalmangainK.
With the predictor Markov parameter estimates, the process system
Markov parameters {g
k
} can be obtained directly from (4.46). Inspired by
subspace identification, we use {g
k
} to construct the Hankel matrix. The
extendedobservabilitymatrixcanbeestimatedthroughtheSingularValue
Decomposition(SVD).
TheHankelmatrixH
sl
isconstructedasfollows
H
sl
=





CB CAB ··· CA
l1
B
CAB CA
2
B ··· CA
l
B
.
.
.
.
.
.
.
.
.
CA
s1
B CA
s
B ··· CA
s+j2
B





(4.48)
wheresandf satisfys+l−2 =i. PerformingSVDon
H
sl
=USV
T
≃U
n
S
n
V
T
n
whereS
n
contains then largest singular values. Γ
s
= U
n
S
1=2
n
is the estima-
tionoftheextendedobservabilitymatrix.
67
Defining =
[
K
T
i1
K
T
i2
··· K
T
is
]
T
,wehave





h
NKF
1
h
NKF
2
.
.
.
h
NKF
s





= Λ· (4.49)
whereΛ = diag{C;CA;:::;CA
s1
}.
Combined with (4.49), (4.45) can be written in the stacked form as
follows
Λ· =





¯
h
NKF
1
¯
h
NKF
2
.
.
.
¯
h
NKF
s





+Θ· (4.50)
whereΘ = diag{Θ
1
;Θ
2
;:::Θ
s
},and
Θ
k
=
{
0; k = 1
[
¯
h
NKF
k1
¯
h
NKF
k2
···
¯
h
NKF
1
]
Γ
k1
; k = 2;:::;s
(4.51)
It is obvious that Λ and Θ can be generated from the estimated predictor
Markov parameters
{
h
NKF
k
}
and the extended observability matrix Γ
s
. The
estimationofcanbeobtainedbysolvingthematrixequation(4.50)
= (Λ−Θ)
y





¯
h
NKF
1
¯
h
NKF
2
.
.
.
¯
h
NKF
s





(4.52)
andtheestimationofthedisturbancesystemMarkovparameteris





h
NKF
1
h
NKF
2
.
.
.
h
NKF
s





= Λ·(Λ−Θ)
y





¯
h
NKF
1
¯
h
NKF
2
.
.
.
¯
h
NKF
s





(4.53)
68
Remark4.4. TheNKFbaseddisturbancePMPandSMPconversionrelation
in(4.53)takesadvantageofthesysteminformationretrievedfromthepro-
cess PMPs. The PMPs of the process are converted to SMPs as a first step.
The process SMPs are used to estimate the extended observability matrix,
basedonwhichthedisturbanceSMPscanbeobtainedrecursively. Theidea
of this algorithm originated from the deterministic realization theory [25].
The estimation of the extended observability matrix, as the second stage
modeldiscussedinChapter3,isastandardstepformostSIMs.
It is important to determine the condition for the proposed distur-
bance modeling scheme. Since the disturbance SMPs are obtained based
on the knowledge of the extended observability matrix, which by itself is
estimated from the process, it is obvious that the process and disturbance
have to share poles. This condition needs to be satisfied in the practical
implementationoftheproposeddisturbancemodelingscheme.
4.5 DMC with NKF Disturbance Models
ThepredictedoutputsofDMCcanberewrittenas
ˆ y(t+k|t) =
k
∑
i=1
s
i
∆u(t+k−i)+
P
∑
i=k+1
s
i
∆u(t+k−i)
+y
0
+
ˆ
d(t+k|t) (4.54)
whereP isthepredictionhorizon. Optimalinputsequencearedetermined
fromanoptimizationproblembasedontheprediction(4.54).
69
It is convenient to differentiate input and output data in order to
identifyanintegratingdisturbancemodel. Convertingstepresponsemodel
oftheplanttoimpulseresponsemodelandthendifferentiatingit,
∆y(t) =
N
∑
i=1
g
p
i
∆u(t−i)+
N
∑
i=0
∆h
i
e(t−i): (4.55)
In this model, the differentiated disturbance ∆h
i
is stable as its integrating
factoriscancelledbythedifferentiation. Thus,wecandirectlyidentify∆h
i
byapplyingtheNKF.
Inordertocalculatepredictions,weintegrate(4.55):
y(t) =
N
∑
i=1
s
i
∆u(t−i)+y
0
+
N
∑
i=0
h
i
e(t−i)+d
0
(4.56)
where d
0
is the initial condition for disturbance. Therefore, k-step-ahead
predictioncanbewrittenas
ˆ y(t+k|t) =
k
∑
i=1
s
i
∆u(t+k−i)+
P
∑
i=k+1
s
i
∆u(t+k−i)
+y
0
+
ˆ
d(t+k|t) (4.57)
ˆ
d(t+k|t) =
N
∑
i=0
h
i
e(t−i)+d
0
: (4.58)
TheoptimalinputmovescanbeoptimizedbyQPproblemwiththepredic-
tionmodel(4.57)and(4.58).
70
4.6 Simulations
In the simulation, we use a SISO plant described by the following
systemmatrices.
A =
[
0:6 0:3
0 0:2
]
; B =
[
0
1
]
C =
[
1 1
]
; D = 0:
(4.59)
Wegenerate10000datapointsusinganarbitrarilydesignedMPCin
MATLAB. The setpoint is 0, and a PRBS signal is added to the setpoint to
improve closed-loop system identifiability. The order of HOARX model is
settobei = 40,andthesizeofΓ
s
inNKFparameterizationiss = 20.
TocomparetheDMCwithNKFtoDMCwithsteady-stateKF,werun
both controllers with 1000 steps. The steady-state KF algorithm is adopted
from Sun et al. [64], which identifies a non-parametric disturbance model.
The prediction and control horizon are chosen as P = 15 and M = 4, re-
spectively. We put zero weights on input and input move and unit weight
onoutputsothatoutputvariancecanbeservedasperformancemetric. The
setpointisalways0throughoutthesimulation.
4.6.1 Case I: no plant model mismatch, subject to ARIMA(1,1,1) dis-
turbance
Nowconsiderthefollowingadditiveoutputdisturbance
d(t) =
1−0:3q
1
1−q
1
1
1−0:95q
1
e(t) (4.60)
wheree(t)∼N(0;1).
71
2 4 6 8 10 12 14 16 18 20
0
2
4
6
8
10
2 4 6 8 10 12 14 16 18 20
0
0.5
1
1.5
Figure 4.1: Identified system Markov parameters of the disturbance model
in Case I. The upper one ish
i
of steady-state KF, and the lower one ish
i
of
NKF.
Estimated system Markov parameters of disturbance model using
steady-stateKFandNKFisdemonstratedbyFig. 4.1. Fig. 4.2indicatesthat
we need the order to be at least 80 to capture the dynamics of the distur-
bance. TheKFwithorder40isfarfromconvergingtothesteadystategain.
Therefore,theparameterizationwithsteady-stateKFyieldsfailureinMPC
20 40 60 80 100
0
1
2
3
4
5
6
Figure4.2: SystemMarkovparametersofthedisturbancemodel.
72
Table4.1: ControlPerformanceofDifferentMPCs
Outputvariance
DMCw/st. st. KF DMCw/NKF
CaseI unstable 16.02
CaseII unstable 8.97
whiletheparameterizationwithNKFleadstogoodcontrolperformance.
4.6.2 Case II: with plant model mismatch, subject to ARIMA(1,1,1)
disturbance
We further add plant model mismatch to the simulation. Let the
model(4.59)bemultipliedby×(1−0:4q
1
)=(1−0:3q
1
). Itisthenconverted
tostepdisturbancemodelandfedintoDMC.Therefore,thesystemmatrices
ofplantmodelcanbedescribedby
A =


0:6 0:3 0
0 0:2 0:4
0 0 0:5


; B =


0
1:2
0:3


C =
[
1 1 0
]
; D = 0:
whose order is higher than that of the plant. There is a pole at q = 0:5 in
additiontothepolesq = 0:6andq = 0:2oftheplant. Thedisturbancetobe
addedisstill(4.60).
Estimated system Markov parameters of disturbance model using
steady-state KF and NKF is demonstrated by Fig. 4.3. Similar to Case I,
steady-stateKFparameterizationfailstostabilizetheoutputwhiletheNKF
73
2 4 6 8 10 12 14 16 18 20
0
2
4
6
8
10
2 4 6 8 10 12 14 16 18 20
0
0.5
1
1.5
Figure 4.3: Identified system Markov parameters of the disturbance model
inCaseII.Theupperoneish
i
ofsteady-stateKF,andtheloweroneish
i
of
NKF.
parameterizationresultsinastablecontrolwiththeoutputvariance8.97.
4.7 Summary
In this chapter, we show that the least squares estimation of the
Markov parameters under the HOARX model represents a non-steady
state Kalman filtering process. The NKF based parameterization of the
Markov parameters is proposed. We further develop the algorithm to
convert the PMPs to the SMPs. The obtained SMPs are incorporated to
DMC for disturbance modeling. The simulation results indicate that the
NKFparameterizationoutperformsthetraditionalSKFparameterizationin
capturingthesystemdynamicswithfinitelengthofdata. Asthedrawback,
the proposed NKF based disturbance modeling scheme requires shared
74
poles of disturbance and process. This condition is conservative and to a
large extent limits the scope of application. This motivates us to develop
universalalgorithmwhichisintroducedinChapter5.
75
Chapter 5
NKF Parameterization based Disturbance Model-
ing and MPC Applications: Part II
5.1 Introduction
InrecentyearsMPCtechnologyhasgainedincreasingpopularityin
applications and has become the most frequently used advanced process
control technology in industry ([5]). Being well known for its capability in
controlling multivariable process and handling constraints, MPC has been
adopted by a wide range of industries including refining, chemical, auto-
motiveandaerospace([56]).
Asamodel-basedcontrolmethod, MPCreliesheavilyonthemodel
accuracy. Eitherplantmodelmismatchorunmodeleddisturbancemayde-
grade the control performance ([24]). However, the modeling and identifi-
cationoftheprocessanddisturbancearetypicallythemosttime-consuming
steps for MPC implementation. Much effort has been made to achieve a
suitableprocessmodelinordertoobtainaccurateoutputpredictions. Cut-
ler and Hawkins [9] used data from step tests to derive the finite impulse
response (FIR) model. Zang and Bitmead ([76]) developed a control rele-
vantidentificationschemeinordertoachieveoptimalcontrolperformance.
76
This approach uses a proper data prefilter before identification to benefit
thecontrollerdesign.
Compared with the plant model, disturbance models receive
relatively little attention in the past. Dynamic matrix control (DMC)
and quadratic dynamic matrix control (QDMC) use a step-like output
disturbance model to achieve offset-free control. Muske and Badgwell
[49] proposed a block diagonal step disturbance model to include both
state and output disturbances. Pannocchia and Rawlings [51] provided a
more general disturbance model without special structures. However, the
filter gain of the state estimator is pre-determined. To avoid specifying
the disturbances at the input or the output side, the equivalence between
input and output step-like disturbance models is shown in [59], where a
method for estimating the disturbance gain is established. Recently, Xu et
al. [75]proposedanadaptivedisturbancemodelingmethod,butitimposes
a special structure on the Kalman filter. Although these methods aim to
achieveimprovedpredictions,thedisturbancemodelsaresimplisticwhich
are usually inadequate in capturing the complex disturbance dynamics in
practicalindustrialprocesses.
ThewiderangeoftheMPCimplementationdemandsmorecomplex
andaccuratedisturbancemodelingwithlittleparametertuning,whichcalls
for the contribution from the recent advancement of system identification.
77
Overthelastfewdecades,subspaceidentificationmethods(SIMs),basedon
the state space representation, have become a dominating research stream
in system identification. Compared to the input-output system identifica-
tion methods, SIMs completely bypass the canonical parameterization and
are suitable for multivariable process identification [53]. Unlike the pre-
diction error methods (PEMs), SIMs do not require nonlinear optimization
and have better numerical reliability and modest computational complexi-
ty. Inaddition,tremendousefforthasbeenmadetodeveloptheclosed-loop
SIMssothattheidentificationcanbecarriedoutundertheclosed-loopcon-
dition. Several representative closed-loop SIMs were proposed during the
lastdecade([45],[68],[67],[34],[29],[73],[28],[7],[54]).
The first step of most closed-loop SIMs is the pre-estimation ([29],
[34],[41],[61]). Inthisstep,thepredictorMarkovparameters(PMPs)arees-
timatedvialeastsquaresregression. ThesystemMarkovparameters(SMP-
s) can be recovered based on PMP estimates, with which the FIR process
and disturbance models can be established. Although the process and dis-
turbance models of the FIR form are the intermediate results in subspace
identification, they are convenient models for the MPC implementations.
Hale and Qin [23] and Huang and Kadali [27] incorporated the interme-
diate models of the FIR form from SIMs with MPC to formulate subspace
model based MPC. Sun et al. [64] use the disturbance models from SIMs
78
to enhance the prediction capability of DMC, which typically uses a step-
like disturbance model. The enhanced disturbance model has the benefit
indescribingcomplexdisturbancedynamicsandcompensatingtheprocess
model mismatch, while the process model is unchanged to avoid process
identificationtesting.
Most industrial MPC implementations use the FIR models for com-
putation efficiency, which will be shown later in this chapter to be exactly
related to PMPs as intermediate results of subspace identification. Since
the PMP estimates are non-parametric, most existing SIM methods ([34])
interpret and further parameterize them with SKF. However, the SKF pa-
rameterization inherently assumes infinite data horizon which cannot be
satisfied in practice. Zhao, Sun and Qin ([77], [57]) show that interpreting
finite data with SKF can lead to suboptimal solution and degrade the con-
trolperformance. Ontheotherhand,NKFwithtime-varyingKalmangain
works with any prediction horizon length. To avoid the suboptimal effec-
t with SKF, we develop the NKF based parameterization for the practical
finitedatalength.
In this chapter, the SKF based process and disturbance modeling
approach is introduced first. We show that due to the finite data horizon
in disturbance estimation, the SKF parameterization based modeling is
biased in practice. To address this, we propose the NKF based modeling
79
approach. Both SKF and NKF parameterized disturbance models are
incorporated with MPC to improve its control performance. The proposed
novel MPC schemes are compared with the traditional DMC which uses
stepdisturbancemodels.
It is worth noting that once the process and disturbance models are
estimatedorgiven,onecanchoosetouseinfinitehorizonMPCformulation
forthecontrolimplementation. Thiscontrolhorizonshouldnotbeconfused
withthedatahorizonusedindisturbancemodelestimationinthischapter.
One may argue that for continuous plant operations there are eventually
as many data as one would like to have. However, in practice, the distur-
bancedynamicsarehardlytime-invariantoveranextendedperiodoftime.
Therefore, it is impractical to assume that distant past data share the same
disturbance dynamics with the current plant disturbances. For these rea-
sonsitisbeneficialtousefinalhorizondatafordisturbanceestimationand
prediction. It is also beneficial to re-estimate the disturbance models from
timetotimetocapturetime-varyingdynamicdisturbances.
Thereminderofthischapterisorganizedasfollows. InSection2,the
process and disturbance modeling scheme is developed in the framework
of SKF. In Section 3, the NKF framework is established and it is used to
parameterze the process and disturbance models. The disturbance model
is applied to the MPC in Section 4. As the case studies, the control perfor-
80
mances of the proposed MPC scheme on the Wood-Berry distillation col-
umnandTennesseeEastmanProcessarereportedinSection5,followedby
theconclusionsinSection6.
5.2 Deriving the Process and Disturbance Models
5.2.1 Estimating PMPs
It is well known that for observable systems, the input-output rela-
tionscanbedescribedintheinnovationform,
ˆ x
k+1
=Aˆ x
k
+Bu
k
+Ke
k
y
k
=Cˆ x
k
+e
k
(5.1)
where K is the steady-state Kalman filter obtained from an algebraic Ri-
catti equation and the innovation e
k
is stationary, zero mean white noise
independentofthepastinputandoutput. Weusex
k
insteadof ˆ x
k+1
inthe
remainder of the chapter for simplicity. By substitution, the system in (6.2)
canberepresentedbytheequivalentpredictorform,
x
k+1
=
¯
Ax
k
+
[
B K
]
[
u
k
y
k
]
(5.2)
y
k
=Cx
k
+e
k
(5.3)
where
¯
A = A− KC. By iterating (6.3) it is straightforward to derive the
followingrelation,
x
k
=
¯
A
p
x
kp
+
¯
L
p
z
p
(k);
where p is the past horizon, the predictor controllability matrix is
¯
L
p
=
[
¯
A
p1
B :::
¯
AB B
¯
A
p1
K :::
¯
AK K
]
. Thepastinputisdefinedas
81
u
p
(k) =
[
u
T
kp
::: u
T
k2
u
T
k1
]
T
and y
p
(k) is formed similarly. z
p
(k) =
[
u
p
(k)
T
;y
p
(k)
T
]
isthecombinedpastinputandoutput.
Ifp is sufficiently large,x
k
can be represented by the past input and
outputas
x
k
=
¯
L
p
z
p
(k): (5.4)
Based on (6.3), (6.4) and (6.5), the predictor form extended state space can
beformedas
y
f
(k) =
¯
Γ
f
¯
L
p
z
p
(k)+
¯
Gu
f
(k)+
¯
Hy
f
(k)+e
f
(k)
=
¯
H
fp
z
p
(k)+
¯
Gu
f
(k)+
¯
Hy
f
(k)+e
f
(k) (5.5)
wherethepredictorobservabilitymatrix
¯
Γ
f
is
¯
Γ
f
=





C
C
¯
A
.
.
.
C
¯
A
f1





;
The stacked future output is y
f
(k) =
[
y
T
k
y
T
k+1
::: y
T
k+f1
]
T
, where f
is the future horizon. u
f
(k) ande
f
(k) are formed similarly. Interpreted by
SKF,thePMPsfortheprocessanddisturbancecanbeparameterizedas¯ g
i
=
C
¯
A
i1
B and
¯
h
i
=C
¯
A
i1
K,respectively.
¯
Ghasthefollowingform
¯
G =





0 0 ··· 0
¯ g
1
0 ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
¯ g
f1
··· ¯ g
1
0





;
82
¯
H has the same format as
¯
G except that all the B matrices are replaced by
K. TheHankelmatrix
¯
H
fp
istheproductofthepredictorobservabilityand
controllabilitymatrices,
¯
H
fp
,
¯
Γ
f
¯
L
p
=





¯ g
p
··· ¯ g
2
¯ g
1
¯ g
p+1
··· ¯ g
3
¯ g
2
.
.
.
.
.
.
.
.
.
.
.
.
¯ g
p+f1
··· ¯ g
f+1
¯ g
f

¯
h
p
···
¯
h
2
¯
h
1
¯
h
p+1
···
¯
h
3
¯
h
2
.
.
.
.
.
.
.
.
.
.
.
.
¯
h
p+f1
···
¯
h
f+1
¯
h
f





: (5.6)
Theextendedstatespacemodelin(6.6)isthebasisforpre-estimation
formostclosed-loopSIMs. Theextendedstatespacemodelin(6.6)contains
f ARX models with each one increasing the order by one. Since the PMPs
in the lower order models reappear in the higher order models, Qin and
Ljung [58] proposed that one only needs to perform one least squares esti-
mationonthelastrowtoextractallthePMPs. However,manyresearchers
reported that better model estimates are obtained through estimating mul-
tiple rows simultaneously or separately. Therefore, NKF parameterization
will be used later in this chapter as a precise representation for finite data
horizons.
83
5.2.2 Parameterizing Process and Disturbance Models
Although PMPs can capture the process and disturbance dynamics,
they are difficult to be directly applied to the MPC. According to Qin and
Badgwell[56],mostMPCtechnologiesrelyontheFIRform,whichinstead
usesthesystemMarkovparametersoritsintegralascoefficients. Specifical-
ly, for DMC implementation the pre-estimated PMPs need to be converted
to the SMPs to form the FIR models. Juang, et al., [34] established the rela-
tionshipbetweenPMPsandSMPsbasedontheSKFparameterization.
DenotethesystemMarkovparametersas
{
g
0
= 0
nynu
;
g
i
=CA
i1
B; fori> 1,
and
{
h
0
=I
nyny
;
h
i
=CA
i1
K; fori> 1.
AfterestimatingthePMPsfromthelastrowof(6.6),thesystemMarkovpa-
rameterscanberecoveredrecursivelybasedonthefollowingrelationships:



g
1
= ¯ g
1
;
g
i
= ¯ g
i
+
i1
∑
m=1
¯
h
m
g
im
; fori> 2,
(5.7)



h
1
=
¯
h
1
;
h
i
=
¯
h
i
+
i1
∑
m=1
¯
h
m
h
im
; fori> 2,
(5.8)
The system Markov parameter sequences, which by themselves are
the impulse responses of the process and disturbance, can be easily imple-
84
mentedonDMCafterspecialtreatment. Moredetailswillbedemonstrated
inthenextsection.
It is obvious that the above relationships are only valid under the
SKFparameterization. TheydonotholdunderNKFparameterizationdue
tothetime-varyingstructureofbothPMPsandSMPs.
5.3 NKF Parameterization
The SKF parameterization inherently requires the past horizon p to
beinfinitetoreachsteadyKalmangain. However,thisconditioncannotbe
satisfiedinpracticesincetheindustrialMPCusesthetruncatedorder. NKF
cancompletelyremovethisconstraint. Inthissection,wedeveloptheNKF
basedprocessanddisturbancemodelingscheme.
5.3.1 Estimating PMPs
The predictor form state space model with NKF can be represented
as,
x
k+1
=
¯
A
k
x
k
+
[
B K
k
]
[
u
k
y
k
]
(5.9)
y
k
=Cx
k
+e
k
(5.10)
where
¯
A
k
=A−K
k
C andK
k
isthenon-steadyKalmangain. SimilartoSKF
model,thestatecanberepresentedas,
x
k
=
¯
L
NKF
k;p
z
p
(k); (5.11)
85
where (·)
NKF
indicates (·) is parameterized by NKF. By denoting the tran-
sition matrix
¯
Φ
i;j
=
¯
A
i1
¯
A
i2
:::
¯
A
j
with
¯
Φ
i;i
= I, the NKF parameterized
predictorcontrollabilitymatrixhasthefollowingform,
¯
L
NKF
k;p
=
[
¯
Φ
k;kp+1
B :::
¯
Φ
k;k1
B B

¯
Φ
k;kp+1
K
kp
:::
¯
Φ
k;k1
K
k2
K
k1
]
:
The predictor form extended state space model can be derived as follows
similarto(6.6)
y
f
(k) =
¯
Γ
NKF
f
¯
L
NKF
k;p
z
p
(k)+
¯
G
NKF
u
f
(k)+
¯
H
NKF
y
f
(k)+e
f
(k)
=
¯
H
NKF
fp
z
p
(k)+
¯
G
NKF
u
f
(k)+
¯
H
NKF
y
f
(k)+e
f
(k) (5.12)
WiththeNKFparameterization,thePMPsoftheprocessanddistur-
bance are formed as g
NKF
i;j
= CΦ
i;j+1
B and h
NKF
i;j
= CΦ
i;j+1
K
j
, respectively.
Wehave
¯
H
NKF
fp
,
¯
Γ
NKF
f
¯
L
NKF
k;p
=





g
NKF
k;kp
··· g
NKF
k;k2
g
NKF
k;k1
g
NKF
k+1;kp
··· g
NKF
k+1;k2
g
NKF
k+1;k1
.
.
.
.
.
.
.
.
.
.
.
.
g
NKF
k+f1;kp
··· g
NKF
k+f1;k2
g
NKF
k+f1;k1

h
NKF
k;kp
··· h
NKF
k;k2
h
NKF
k;k1
h
NKF
k+1;kp
··· h
NKF
k+1;k2
h
NKF
k+1;k1
.
.
.
.
.
.
.
.
.
.
.
.
h
NKF
k+f1;kp
··· h
NKF
k+f1;k2
h
NKF
k+f1;k1






: (5.13)
86
The NKF parameterization reveals the fact that a bank of HOARX,
ratherthanonlyoneHOARXmodel,shouldbeestimatedwhenthehorizon
isfinite. Comparedtothe
¯
H
fp
in(6.7),the
¯
H
NKF
fp
doesnotpreservetheHan-
kel structure. Although the row-wise partition of (6.14) still yields f ARX
models, the PMPs in an upper row do not reappear in the lower row. This
is a major difference between the NKF and SKF parameterization. Com-
paredtoSKF,theNKFparameterizedextendedstatespacemodelhasmore
complicated structure. As a result, one needs to estimatef ARX models of
increasingorders. Also,theexistingconversionrelationbetweenPMPsand
SMPs is not valid due to the time-varying structure, and thus additional
effortisneededtoretrieveprocessanddisturbancemodels.
5.3.2 Parameterizing System Markov Parameters
SimilartotheSKFparameterization,thesystemMarkovparameters
canbedenotedas
{
g
NKF
k;k
= 0
nynu
;
g
NKF
k;ki
=CA
i1
B; fori> 1,
and
{
h
NKF
k;k
=I
nyny
;
h
NKF
k;ki
=CA
i1
K
ki
; fori> 1.
Nowthequestionis,withNKFparameterization,howtoconvertthe
PMPstoSMPs. ByexamingthestructureoftheNKFparameterizedpredic-
tor and system Markov parameters, we reveal the general relationship as
87
follows,



g
NKF
k;k1
=g
NKF
k;k1
;
g
NKF
k;ki
=g
NKF
k;ki
+
i1
∑
m=1
h
NKF
k;km
g
NKF
km;ki
; fori> 2,
(5.14)



h
NKF
k;k1
=h
NKF
k;k1
;
h
NKF
k;ki
=h
NKF
k;ki
+
i1
∑
m=1
h
NKF
k;km
h
NKF
km;ki
;
fori> 2,
(5.15)
Therelationshipsin(6.16)and(6.17)enableonetocalculatetheSMPs
recursivelyfromthePMPestimates. Itsworthnotingthatwhenallthetime-
varyingK
i
is replaced by the steady Kalman gainK, the relations in (6.16)
and (6.17) are reduced to (6.8) and (6.9). This means the relations in (6.16)
and(6.17)holdforbothNKFandSKFandtheyaremoregeneralconversion
relationsofMarkovparameters.
5.4 MPC with Improved Disturbance Models
AstepresponseplantmodelareusedinDMCtodopredictions. We
denotey
t
∈R
ny
andu
t
∈R
nu
as the output and input of an MIMO process
respectively. Theinput-outputrelationshipcanbewrittenas[18]
y
t
=
N
∑
i=1
s
i
∆u
ti
+y
0
+d
t
(5.16)
where y
0
is the initial condition for output, d
t
denotes unmeasured distur-
bances, ∆u
j
= u
j
− u
j1
, N is a number large enough for plant to reach
steady-state,ands
i
istheithstepresponsematrix. From(5.16),thepredict-
88
edoutputscanberewrittenas
ˆ y
t+kjt
=
k
∑
i=1
s
i
∆u
t+ki
+
P
∑
i=k+1
s
i
∆u
t+ki
+y
0
+
ˆ
d
t+kjt
(5.17)
whereP isthepredictionhorizon. Optimalinputsequencearedetermined
fromanoptimizationproblembasedontheprediction(5.17).
It is convenient to differentiate input and output data in order to
identifyanintegratingdisturbancemodel. Convertingstepresponsemodel
oftheplanttoimpulseresponsemodelandthendifferentiatingit,
∆y
t
=
N
∑
i=1
g
i
∆u
ti
+
N
∑
i=0
∆h
i
e
ti
: (5.18)
In this model, the differentiated disturbance ∆h
i
is stable as its integrating
factoriscancelledbythedifferentiation. Thus,wecandirectlyidentify∆h
i
byapplyingtheSKForNKFparameterizationbaseddisturbancemodeling
scheme.
Inordertocalculatepredictions,weintegrate(5.18):
y
t
=
N
∑
i=1
s
i
∆u
ti
+y
0
+
N
∑
i=0
h
i
e
ti
+d
0
(5.19)
where d
0
is the initial condition for disturbance. Therefore, k-step-ahead
89
predictioncanbewrittenas
ˆ y
t+kjt
=
k
∑
i=1
s
i
∆u
t+ki
+
P
∑
i=k+1
s
i
∆u
t+ki
+y
0
+
ˆ
d
t+kjt
(5.20)
ˆ
d
t+kjt
=
N
∑
i=0
h
i
e
ti
+d
0
: (5.21)
TheoptimalinputmovescanbeoptimizedbyQPproblemwiththepredic-
tionmodel(5.20)and(5.21).
5.5 Case Studies
5.5.1 Wood-Berry Distillation Column
The advantages of the proposed disturbance modeling methods are
firstdemonstratedontheWood-Berrydistillationcolumnmodel[74]. With
asamplingrateofoneminute,thediscretizedprocesstransferfunctionma-
trixisdescribedby
G(q) =



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



:
Thefirstinputu
1
isrefluxrateinlb/min;thesecondinputu
2
issteam
flowrateinlb/min;thefirstoutputy
1
iscompositionoftopproductofcol-
umninmol %; and thesecond outputy
2
iscompositionof bottomproduct
ofcolumninmol%.
90
Thedisturbancemodelischosentobediagonal
H(q) =




1−0:5q
1
1−q
1
1
1−0:95q
1
0
0
1−0:7q
1
1−q
1
1
1−0:9q
1




(5.22)
The innovation e
k
is independent white noise with covariance matrix as
diag{0:19
2
;0:32
2
}. Thesetpointsareselectedas ¯ y
1
= 90%and ¯ y
2
= 5%.
ThecolumniscontrolledbytheMPCwiththepredictionandcontrol
horizon as 40 and 3, respectively. The weighting matrices of the MPC are
Q = diag{1;100};S = 0. No constraint is enforced for simplicity, thus the
MPCis an LTIcontroller. The MPCs with step-like disturbance, SKF based
disturbancemodelandNKFbaseddisturbancemodelwillbeimplemented
forcontrolperformancecomparison.
Fig. 5.1 shows the impulse response of the true disturbance mod-
el. It takes 100 time steps for the disturbance subsystem u
1
to y
1
to reach
steady state and 20 time steps for the subsystem u
2
to y
2
to reach steady
state. Theother2subsystemsshouldhavezerodynamics. Theimpulsere-
sponse of three types of disturbance models with order chosen to be 8 and
30aredemonstratedinFig. 5.2and5.3,respectively. Inthiscaseweanalyze
theeffectoforderreduction,orequivalentlythedatahorizonreduction,on
thedisturbancemodelingandtheMPCcontrolperformance. Whentheor-
deris30,bothSKFandNKFprovidesthedisturbancemodelwithdynamics
91
0
2
4
6
8
10
From: In(1)
To: Out(1)
0 20 40 60 80 100 120
−1
−0.5
0
0.5
1
1.5
2
2.5
3
To: Out(2)
From: In(2)
0 20 40 60 80 100 120
Disturbance model (true)
Time (sec)
Amplitude
Figure5.1: Impulseresponseofthenominaldisturbancemodel.
close to the nominal. When the order is reduced to 8, the SKF parameter-
ization based disturbance model degrades significantly, especially for the
subsystem u
2
to y
1
. On the contrary, NKF parameterization still provides
robustdisturbancedynamics.
The control performances are demonstrated in Table 5.1. In all the
cases,theMPCswithSKFandNKFparameterizeddisturbancemodelsout-
performtheMPCwithstep-likedisturbancemodel. Itisconsistentwiththe
fact that the SKF and NKF modeled disturbances capture more dynamics
92
0
0.5
1
1.5
2
2.5
3
3.5
4
From: In(1)
To: Out(1)
0 5 10 15 20
−0.5
0
0.5
1
1.5
2
To: Out(2)
From: In(2)
0 5 10 15 20

Disturbance model for order = 8
Time (sec)
Amplitude
Step−like
SKF
NKF
Figure5.2: Impulseresponseofstep-like,SKFandNKFbaseddisturbance:
order=8.
93
0
2
4
6
8
10
From: In(1)
To: Out(1)
0 5 10 15 20 25 30 35
0
0.5
1
1.5
2
2.5
3
To: Out(2)
From: In(2)
0 5 10 15 20 25 30 35

Disturbance model for order = 30
Time (sec)
Amplitude
Step−like
SKF
NKF
Figure5.3: Impulseresponseofstep-like,SKFandNKFbaseddisturbance:
order=30.
94
information,hencetheMPChasbettermodelsforpredictionandoptimiza-
tion. The result also confirms that for the MPC with SKF or NKF parame-
terizeddisturbance,thecontrolperformanceimproveswiththeincreaseof
horizon.
The difference between the control performance of the MPCs with
SKF and NKF parameterized disturbances becomes obvious with low or-
ders. The result is consistent with the Kalman filtering theory in the sense
thatforhighorder,NKFapproachesSKFasymptotically,butwhenthehori-
zonisnotsufficientlylarge,theMPCwithNKFparameterizeddisturbance
showssignificantadvantage.
Table5.1: ControlPerformanceofMPCswithDifferentDisturbanceModels
andVariousOrders
KPI
Order Step-like SKF NKF
11 103.1823 92.0376 84.9618
12 103.1823 89.7329 83.2468
15 103.1823 86.3617 82.2018
20 103.1823 82.5650 81.2663
25 103.1823 81.9324 80.7905
30 103.1823 80.2507 79.5972
40 103.1823 79.8911 79.6470
5.5.2 Tennessee Eastman process
BeingabenchmarkproblemoftheEastmanChemicalcompany[13],
theTennesseeEastmanprocess(TEP)iswidelyusedbymanyprocesscon-
95
trolresearcherstoevaluateprocesscontrolandsystemidentificationmeth-
ods. In Fig. 5.4, a brief diagram is shown. The TEP is composed of eight
components: twomainproducts,fourreactants,onebyproductandonein-
ert. Therearefiveunitsintheprocess: onereactor,onecondenser,onecom-
pressor, one separator and one stripper. In the reactor, fed by four gaseous
reactants,twomainproductsareformedalongwiththebyproductandthe
inert. In this process, there are 12 manipulated variables and 41 measured
variables. AmongtheMVs,19ofthemarequalityvariablessampledevery
3minutes; therest22MVsaremeasuredcontinuously. Theprocesshassix
operatingmodes;Mode3isselectedinthisexample.
Due to the open-loop instability of the TEP, it is hard to design con-
trollers to regulate this process. Some previous research work successfully
carriedoutavarietyofdecentralizedcontrolstrategies,suchas[60,47,46].
Itisreportedthattheprocessisdifficulttocontrolusingacentralizednon-
linear MPC [60]. Therefore, in this example, only one PID control loop is
substitutedbyanMPCtoillustratethecontrolperformanceimprovements
benefitedfromNKFparameterizedmodel. Thecontrolledandmanipulated
variables chosen are separator temperature (CV 11) and condenser cooling
water flow (MV 11). In the control settings by Ricker [60], this pair of MV
and CV works together with MV 10 (reactor cooling water flow) and CV
9 (reactor temperature) to form a cascade control loop. The MV/CV pair
96
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
Figure5.4: SchematicdiagramoftheTennesseeEastmanprocess.
selected is the outer loop, which usually has slow plant and disturbance
dynamics. Ontheotherhand,MV10/CV9innerloop,whichistheplantof
outer loop, responds quickly as a result of the inner PID controller. There-
fore, this unique combination has all the desired features of NKF parame-
terization: fast plant dynamics, slow disturbance dynamics. An MPC will
replace the PID controller for the MV 11/CV 11 pair to show the advan-
tages of NKF model in industrial applications. The rest of CVs are under
the control of the decentralized PID controllers designed by Ricker [60].
The prediction and control horizon are chosen to be 6 and 2, respectively.
97
Thepasthorizonpandfuturehorizonf arebothsetas6forSKFandNKF
parameterization. In this example, a total of 72 hours will be simulated in
MATLAB,anddataaresampledevery36seconds(0.01hour).
Fig. 5.5 illustrates the open-loop response of the subsystem, which
confirms that the disturbance is almost integrating and has very slow dy-
namics. The proposed NKF based modeling method is then applied to the
closed-loop data to estimate the plant model and the NKF based distur-
bancemodel. AsshowninFig. 5.6,theplantmodelestimatedsuggeststhat
the output responds very quickly to the input; it only takes two time steps
to reach the steady state. This plant model is shared by all the MPCs to be
simulated. Fig. 5.7, 5.8 and 5.9 show the impulse response of disturbance
modelofMPCusingstep-likedisturbance,SKFparameterizeddisturbance
andNKFparameterizeddisturbance,respectively.
The simulation results are demonstrated in Fig. 5.10, 5.11 and 5.12.
To make a fair comparison, the MPCs are designed to be unconstrained to
avoid the impact of constraints, which allows the MV to exceed 100% in
Fig. 5.11. It is as expected that the MPC with a SKF parameterized distur-
bancemodeloutperformstheMPCwithstep-likedisturbancemodel,since
the latter one assumes no knowledge about unmodeled disturbances. The
MPC with NKF disturbance model achieves an even better performance:
thevariancesofbothMVandCVaresmallerwhenencounteringlargedis-
98
0 10 20 30 40 50 60 70
82.5
83
83.5
84
84.5
85
condenser cooling water flow/%
0 10 20 30 40 50 60 70
65
70
75
time/hour
separator temp/
o
C
Figure5.5: Open-loopresponseoftheselectedMV/CVpair.
99
0 0.05 0.1 0.15 0.2
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
Figure5.6: ImpulseresponseofplantmodelofMPCcontrollingTEP.
0 1 2 3 4 5 6
0
0.2
0.4
0.6
0.8
1
Figure 5.7: Impulse response of disturbance model: step-like disturbance
model.
100
0 1 2 3 4 5 6
0
0.2
0.4
0.6
0.8
1
Figure5.8: Impulseresponseofdisturbancemodel: SKFparameterizeddis-
turbancemodel.
0 1 2 3 4 5 6
0
0.2
0.4
0.6
0.8
1
Figure 5.9: Impulse response of disturbance model: NKF parameterized
disturbancemodel.
101
turbance at about t = 30hr. The statistical results presented in Table 5.2
confirm the observations of MV/CV data above. Since the MPCs are un-
constrained, the average KPI of MPC, defined as the mean of MPC cost
function, can be adopted as a performance metric. The output variance is
also listed in the table for reference. The statistics suggests that the NK-
F parameterization based MPC is able to improve the control performance
considerablyascomparedtotheconventionalSKFparameterizationbased
and step-disturbance based MPC. Therefore, NKF modeling technique can
beanidealcandidateforMPCdesignunderthecircumstanceoffinitedata
length.
Table5.2: ControlPerformanceofMPCswithDifferentDisturbanceModels
KPI OutputVariance
MPCw. step-likedisturbancemodel 0.7903 0.3414
MPCw. SKFbaseddisturbancemodel 0.6378 0.3412
MPCw. NKFbaseddisturbancemodel 0.4642 0.2962
5.6 Summary
The industrially popular DMC technique uses a step disturbance
model, which is insufficient to capture complex disturbance dynamics. In
this paper, we proposed integrated moving average disturbance model
parameterizedbySKForNKF.Usingrecentresearchresultsfromsubspace
identification, time-varying Kalman gains are estimated to form the FIR
102
0 10 20 30 40 50 60 70
50
60
70
80
90
100
condenser cooling water flow/%
0 10 20 30 40 50 60 70
62
63
64
65
66
67
68
time/hour
separator temp/
o
C
Figure5.10: ControlresultsofTEPbyMPCwithstep-likedisturbancemod-
el.
103
0 10 20 30 40 50 60 70
50
60
70
80
90
100
condenser cooling water flow/%
0 10 20 30 40 50 60 70
62
63
64
65
66
67
68
time/hour
separator temp/
o
C
Figure5.11: ControlresultsofTEPbyMPCwithSKFparameterizeddistur-
bancemodel.
104
0 10 20 30 40 50 60 70
50
60
70
80
90
100
condenser cooling water flow/%
0 10 20 30 40 50 60 70
62
63
64
65
66
67
68
separator temp/
o
C
time/hour
Figure 5.12: Control results of TEP by MPC with NKF parameterized dis-
turbancemodel.
105
model of the process and disturbance. The obtained disturbance model
wasincorporatedtoDMCtoimprovethecontrolperformance.
We interpret the predictor Markov parameters with NKF to extrac-
t the system parameters for use in MPC. The advantages of the proposed
DMC with SKF and NKF parameterized disturbance model over the tradi-
tional DMC are successfully demonstrated on the Wood-Berry Distillation
Column and the Tennessee Eastman Process. Since industrial disturbance
dynamics are hardly time-invariant over an extended period of time, it is
important not to use distant past data to estimate the current plant distur-
bances. For this reason, it is beneficial to use finite horizon data for distur-
bance estimation and to re-estimate the disturbance models from time to
timetocapturetime-varyingdynamicdisturbances.
106
Chapter 6
NKF Parameterization based Subspace Identifica-
tion
6.1 Introduction
Over the last forty years, subspace identification methods (SIMs),
based on projection techniques in the Euclidean space, have become a
dominating research stream in system identification. Originated in the
Ho-Kalman paper [25], subspace identification went through the develop-
mentofstochasticrealizationtheory([1]),andthesolutionofthecombined
deterministic-stochastic realization problem ( [8], [12], [39], [50], [69], [71]).
Several representative algorithms have been developed for open-loop
identification during the last twenty years, including canonical variate
analysis (CVA, [39]), numerical algorithm of subspace state space system
identification (N4SID, [65]) and multivariate output-error state space
(MOESP,[69]). TheunifyingtheorembyVanOverscheeandDeMoor([66])
provides a framework in which these algorithms can be interpreted as a
singularvaluedecompositionofvariousweightmatrices.
However,theaforementionedSIMsarebiasedforclosedloopidenti-
ficationunlessspecialtreatmentsareintroduced. Closed-loopidentification
107
is of special interest for engineering applications. Due to safety and qual-
ity requirements, it is preferred that identification experiments are carried
out under closed-loop condition. Unlike prediction error methods (PEM-
s), the open-loop SIMs (e.g., CVA, N4SID and MOESP) are biased under
closed-loopconditionduetothecorrelationbetweenthenoiseandthecon-
trol input. To solve the closed-loop identification problem, several closed-
loopsubspaceidentificationmethodswereproposedduringthelastdecade
([45], [68], [67]). The recent developments are presented in ([34], [29], [73],
[28], [7], [54]). The OKID method ([34]) estimates the Markov parameter-
s from the high order ARX (HOARX) model and forms the Hankel ma-
trix,basedonwhichthesystemmatricesarerecoveredviaERA([33]). The
SSARX method ([29]) used the pre-estimated Markov parameters to con-
struct
¯
G
f
sothat thefutureinputeffectwasremovedfromtheoutput. Qin
and Ljung ([54]) developed the innovation estimation method which pre-
estimates the E
f
from the recursive regression on the ARMAX model. In
the whitening filter approach (WFA, [7]), Chiuso and Picci ([7]) used the
row-wise regression on the predictor form to avoid the problem of corre-
lation between input and noise. Wang and Qin ([73]) proposed the parity
spaceandprincipalcomponentanalysis(PCA)forerror-in-variableclosed-
loop identification. Chiuso and Picci ([7]) provided theoretical analysis of
these methods. Qin, et al. ([57]) proposed a progressive parameterization
108
frameworktointerprettheproceduresofclosed-loopSIMs.
Ageneralclosed-loopsubspaceidentificationalgorithmusuallygoes
throughthefollowingsteps:
1. Pre-estimation: obtain Markov parameter estimates from HOARX,
which is non-parametric and has its own application in practice (for
example,indynamicmatrixcontrol,[64]).
2. Model reduction: perform SVD on the Hankel matrix or weighted
Hankel matrix, leading to extended observability matrix
¯
Γ
f
and ex-
tended controllability matrix
¯
L
p
estimates, which by themselves are
reducedordernon-parametricmodels.
3. Systemmatricesidentification: recoverthesystemmatricesA;B;C;K
fromtheextendedobservabilityandcontrollabilitymatricesordirect-
lyfromtheestimatedstatesequence.
In the pre-estimation ([29], [34], [41], [61]) step, the predictor
Markov parameters are estimated via least squares regression. These
non-parametricestimatesareusuallyinterpretedandfurtherparametrized
with steady-state Kalman filter (SKF) as the predictor Markov parameters.
However, inspired by the similarity between the objective functions of
the least squares and the non-steady state Kalman filter (NKF), Sorenson
([63]) showed that NKF represents a recursive solution to least squares
109
problems. VanOverscheeandDeMoor([65])provedthattheleastsquares
regression of the future output on a finite horizon past input and output
datacanbeinterpretedasabankofNKF.Recently,Zhao,SunandQin([77],
[57]) pointed out that steady-state Kalman filter parameterization is not
consistentwiththefinitedatalengthandthusdoesnotleadtotheoptimal
solution. The employment of the NKF structure is necessary to obtain the
optimalestimationinpracticalcaseswhenthedatahorizonisfinite.
Inthischapter,weshowthattheSKFbasedSIMsinherentlyrequires
thedatahorizontobeinfinitewhichcannotbesatisfiedinpractice,andthus
the implementation on finite data yields sub-optimal estimation. Comply-
ing with the finite data horizon, a SIM with NKF parameterization is pro-
posed to address the practical closed-loop identification problem. In ad-
dition, the relationship between predictor Markov parameters and system
MarkovparametersisestablishedundertheNKFframework.
Thereminderofthischapterisorganizedasfollows. Thetraditional
SKFparameterizationisfirstgivenfortheeaseofdefiningnotationandthe
conversionrelationshipbetweenpredictorandsystemMarkovparameters
is introduced in Section 2. The NKF parameterization and the correspond-
ing subspace models are discussed in Section 3. An NKF parameterization
based closed-loop subspace identification method is developed in Section
4, followed by the simulation results in Section 5 and the conclusions in
110
Section6.
6.2 Steady State Kalman Filter Parameterization
6.2.1 Estimating Predictor Markov Parameters
The purpose of system identification is to estimate an appropriate
dynamic model for a properly collected time series of system input and
output data {u
k
;y
k
} which can be generated under open-loop or closed-
loop conditions. Assuming the data sequence {u
k
;y
k
} is generated by the
followingstatespaceprocessmodel,
x
k+1
=Ax
k
+Bu
k
+w
k
y
k
=Cx
k
+v
k
(6.1)
wherex
k
∈R
n
. Assuming(C;A)isobservable,wecandesignastablestate
observerasfollows:
ˆ x
k+1
=Aˆ x
k
+Bu
k
+Ke
k
y
k
=Cˆ x
k
+e
k
(6.2)
whereK isthesteady-stateKalmanfilterobtainedfromanalgebraicRicatti
equation. The innovation e
k
is stationary, zero mean white noise which is
independent of the past input and output. In the rest of this chapter, we
usex
k
instead of ˆ x
k
for convenience. (6.2) is also known as the innovation
form. By substitution, the system described by (6.2) can be represented by
thefollowingequivalentpredictorform,
x
k+1
=
¯
Ax
k
+
[
B K
]
[
u
k
y
k
]
(6.3)
y
k
=Cx
k
+e
k
(6.4)
111
where
¯
A = A−KC is stable from the Kalman filter theory. The following
assumptionsareintroducedtoestablishthefoundationsforSIMs:
1. (A;C)isobservableand
(
A;
[
B K
])
iscontrollable.
2. Theinputu
k
andinnovatione
k
areuncorrelatedundertheopen-loop
condition,buttheyarecorrelatedundertheclosed-loopcondition.
3. The input u
k
is quasi-stationary and is persistently exciting of order
f + p, where f and p are the future and past horizons, respectively.
Thedetaileddefinitionoff andpwillbeintroducedlater.
Most existing subspace identification methods iterate (6.3) to derive
thefollowingrelation
x
k
=
¯
A
p
x
kp
+
¯
L
p
z
p
(k);
where x
kp
is an initial state, the predictor controllability matrix has the
form
¯
L
p
=
[
¯
A
p1
B :::
¯
AB B
¯
A
p1
K :::
¯
AK K
]
. Thepastinputis
definedasu
p
(k) =
[
u
T
kp
::: u
T
k2
u
T
k1
]
T
andy
p
(k)isformedsimilarly
asu
p
(k). z
p
(k) =
[
u
p
(k)
T
;y
p
(k)
T
]
T
isthecombinedpastinputandoutput.
However,mostexistingSIMworkneglectsthefactthat,basedonthe
initialx
kp
topredictthestateforwardstepbystep,theoptimalprediction
requiresthe non-steady state Kalman gain sequence rather than the steady
state Kalman gain. The suboptimal effect of using a static Kalman gain is
112
more pronounced when the horizon length p is relatively short. Precisely
speaking the predictor controllability matrix
¯
L
p
must be parameterized by
the sequence of non-steady state Kalman gains to be optimal. In the next
section of this chapter this NKF parametrization is developed with associ-
atedidentificationprocedures.
Ifpissufficientlylargesothattheeffectofx
kp
canbeneglected,the
currentstatex
k
isapproximatelyrepresentedbythepastinputandoutput
x
k
=
¯
L
p
z
p
(k): (6.5)
Basedonthepredictorform(6.3)and(6.4)andthestaterepresentation(6.5),
the extended predictor form can be formulated by by iterating from k to
k+f −1,whichiscompactlyrepresentedasfollows.
y
f
(k) =
¯
Γ
f
x
k
+
¯
Gu
f
(k)+
¯
Hy
f
(k)+e
f
(k)
=
¯
Γ
f
¯
L
p
z
p
(k)+
¯
Gu
f
(k)+
¯
Hy
f
(k)+e
f
(k)
=
¯
H
fp
z
p
(k)+
¯
Gu
f
(k)+
¯
Hy
f
(k)+e
f
(k) (6.6)
wherethepredictorobservabilitymatrixis
¯
Γ
f
=





C
C
¯
A
.
.
.
C
¯
A
f1





;
The future output is y
f
(k) =
[
y
T
k
y
T
k+1
::: y
T
k+f1
]
T
, where f is the fu-
ture horizon. u
f
(k) and e
f
(k) are formed similar to y
f
(k). With the SKF
113
parameterization,thepredictorMarkovparametersfortheprocessanddis-
turbancecanbedenotedas ¯ g
i
=C
¯
A
i1
B and
¯
h
i
= C
¯
A
i1
K, respectively.
¯
G
hasthefollowingform
¯
G =





0 0 ··· 0
¯ g
1
0 ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
¯ g
f1
··· ¯ g
1
0





;
¯
H hasthesameformatas
¯
GexceptthatalltheBmatricesarereplacedbyK.
The Hankel matrix is defined as the product of the predictor observability
andcontrollabilitymatrices,
¯
H
fp
,
¯
Γ
f
¯
L
p
=





¯ g
p
··· ¯ g
2
¯ g
1
¯ g
p+1
··· ¯ g
3
¯ g
2
.
.
.
.
.
.
.
.
.
.
.
.
¯ g
p+f1
··· ¯ g
f+1
¯ g
f

¯
h
p
···
¯
h
2
¯
h
1
¯
h
p+1
···
¯
h
3
¯
h
2
.
.
.
.
.
.
.
.
.
.
.
.
¯
h
p+f1
···
¯
h
f+1
¯
h
f





: (6.7)
Theextendedstatespacemodelin(6.6)hasthefollowingimportant
features:
1. Thecoefficient
¯
H
fp
hastheHankelstructure.
¯
Gand
¯
H arelowertrian-
gularblock-Toeplitz.
2. Each row of (6.6) represents a HOARX model. Overall, the extended
114
state space model (6.6) contains f HOARX models with each one of
increasingorders.
3. As the coefficients of the HOARX, the predictor Markov parameters
in the upper row reappear in the lower row. Based on this, Qin and
Ljung ([58]) proposed to use the parameters of the last row only to
estimateallpredictorMarkovparametersfortheHankelmatrix.
The extended state space model (6.6) is the basis for pre-estimation
for most closed-loop SIMs. For example, both OKID and SSARX focus on
thelastrowof(6.6),whichbyitselfisanARXmodelwithorderp+f−1,to
estimate the predictor Markov parameters. The extended subspace matri-
ces are built with the predictor Markov parameter estimates ([41], [61]). It
willbemadeclearinthenextsectionthatestimatingpredictorMarkovpa-
rametersdoesnotrequirethesteady-stateKalmanfilterstructure,although
most methods inherently employ the steady-state Kalman filter and use it
forparameterization.
The steady-state Kalman filter parameterization suggests that one
only needs to perform one least squares estimation on the last row to ex-
tract all the predictor Markov parameters. However, many researchers ex-
periencedthatitoftenleadstobettermodelestimatesbyestimatingmulti-
ple rows simultaneously or separately. These observations implicitly raise
115
questions about the validity of the steady state Kalman filter parametriza-
tion.
6.2.2 Obtaining System Markov Parameters
Inthemodelreductionstep,theHankelmatrixisformedtocalculate
the subspace matrices, for example, process observability and controllabil-
ity matrices, and/or state sequences. It is usually desirable to estimate the
system Markov parameters to recover the system matrices. An intuitive
ideaafterobtainingpredictorMarkovparametersistoconvertthemtosys-
temMarkovparameters. Juang,etal. ([34])establishedtherelationshipbe-
tweenpredictorMarkovparametersandsystemMarkovparametersbased
ontheSKFparameterization.
DenotethesystemMarkovparametersas
{
g
0
= 0
nynu
;
g
i
=CA
i1
B; fori> 1,
and
{
h
0
=I
nyny
;
h
i
=CA
i1
K; fori> 1.
AfterestimatingthepredictorMarkovparametersfromthelastrowof(6.6),
the system Markov parameters can be calculated recursively based on the
followingrelationships:



g
1
= ¯ g
1
;
g
i
= ¯ g
i
+
i1
∑
m=1
¯
h
m
g
im
; fori> 2,
(6.8)
116



h
1
=
¯
h
1
;
h
i
=
¯
h
i
+
i1
∑
m=1
¯
h
m
h
im
; fori> 2,
(6.9)
The system Markov parameter sequences, which by themselves are
theimpulseresponsesoftheprocessanddisturbancemodel,havetheirown
applicationinthemodelpredictivecontrol.
ItisobviousthattheaboverelationshipsareonlyvalidundertheSK-
F parameterization condition. They do not hold under NKF parameteriza-
tionduetothetime-varyingstructureofbothpredictorandsystemMarkov
parameters.
6.2.3 Subspace Identification with SKF Parameterization
TheOKIDmethod([34])isarepresentativeclosed-loopidentification
method based on the SKF parameterization. Like most closed-loop identi-
fication methods, the first step of OKID is the pre-estimation, in which the
predictorMarkovparametersareestimatedvialeastsquaresestimation. In
thesecondstep,thepredictorMarkovparametersareconvertedtotheSKF
parameterizedsystemMarkovparametersviatherelationshipsin(6.8)and
(6.9). AHankelmatrixisthenbuiltasfollowing:
H
fp
=





g
f1
··· g
2
g
1
g
f
··· g
3
g
2
.
.
.
.
.
.
.
.
.
.
.
.
g
p+f1
··· g
f
g
f1

117
h
f1
··· h
2
h
1
h
f
··· h
3
h
2
.
.
.
.
.
.
.
.
.
.
.
.
h
p+f1
··· h
f
h
f1





: (6.10)
With the SKF parameterization, H
fp
preserves the Hankel structure and
thus it can be formed in this way. The system matricesA,B,C and steady
state Kalman gainK can be identified through the ERA ([33]) or ERA/DC
([32]).
6.3 NKF Parameterization and Estimation
By employing the SKF, the state space model (6.6) inherently as-
sumes an infinite past data horizon for Kalman filter to reach the steady
gain. However, in practice, the length of the data set is always finite hence
theinfinitehorizonconditioncannotbesatisfied. However, thestatespace
model with NKF puts no constraints on the data horizon and provides an
optimalmodelforfinitedatalength.
6.3.1 NKF Markov Parameters
ThestatespacepredictormodelwithNKFcanberepresentedasfol-
lows,
x
k+1
=
¯
A
k
x
k
+
[
B K
k
]
[
u
k
y
k
]
(6.11)
y
k
=Cx
k
+e
k
(6.12)
118
where
¯
A
k
=A−K
k
C andK
k
isthenon-steadyKalmangain. Similartothe
SKF based model, the state can be approximately represented by the past
inputandoutput,
x
k
=
¯
L
NKF
k;p
z
p
(k); (6.13)
where (·)
NKF
indicates (·) is parameterized by NKF. By denoting the tran-
sition matrix
¯
Φ
i;j
=
¯
A
i1
¯
A
i2
:::
¯
A
j
with
¯
Φ
i;i
= I, the NKF parameterized
predictorcontrollabilitymatrixhasthefollowingform,
¯
L
NKF
k;p
=
[
¯
Φ
k;kp+1
B :::
¯
Φ
k;k1
B B

¯
Φ
k;kp+1
K
kp
:::
¯
Φ
k;k1
K
k2
K
k1
]
:
Thecorrespondingpredictorformextendedstatespacemodelis
y
f
(k) =
¯
Γ
NKF
f
¯
L
NKF
k;p
z
p
(k)+
¯
G
NKF
u
f
(k)+
¯
H
NKF
y
f
(k)+e
f
(k)
=
¯
H
NKF
fp
z
p
(k)+
¯
G
NKF
z
f1
(k)+
¯
H
NKF
y
f
(k)+e
f
(k) (6.14)
ThepredictorMarkovparametersoftheprocessanddisturbancecan
bedenotedasg
NKF
i;j
=CΦ
i;j+1
Bandh
NKF
i;j
=CΦ
i;j+1
K
j
,respectively. Wehave
¯
H
NKF
fp
,
¯
Γ
NKF
f
¯
L
NKF
k;p
=





g
NKF
k;kp
··· g
NKF
k;k2
g
NKF
k;k1
g
NKF
k+1;kp
··· g
NKF
k+1;k2
g
NKF
k+1;k1
.
.
.
.
.
.
.
.
.
.
.
.
g
NKF
k+f1;kp
··· g
NKF
k+f1;k2
g
NKF
k+f1;k1

119
h
NKF
k;kp
··· h
NKF
k;k2
h
NKF
k;k1
h
NKF
k+1;kp
··· h
NKF
k+1;k2
h
NKF
k+1;k1
.
.
.
.
.
.
.
.
.
.
.
.
h
NKF
k+f1;kp
··· h
NKF
k+f1;k2
h
NKF
k+f1;k1






: (6.15)
Compared to the SKF parameterized extended state space model
(6.6), the extended state space model with NKF parameterization has the
followingfeatures,
1.
¯
H
NKF
fp
does not have the Hankel structure.
¯
G
NKF
and
¯
H
NKF
are not
lowertriangularblock-Toeplitz.
2. WithNKFparameterization,
¯
H
NKF
fp
isstilltheproductofpredictorcon-
trollabilityandobservabilitymatrices.
3. Similar to the SKF parameterization, a row-wise partition of (6.14)
yields f HOARX models. However, unlike the SKF model, the
predictor Markov parameters in the upper row do not appear in the
lowerrow.
The NKF parameterization reveals the fact that a bank of HOARX
models, rather than only one HOARX model, should be estimated when
the horizon is finite. In contrast with the SKF parameterization, the pre-
dictor Markov parameter sequence in the upper row does not reappear in
thelowerrow,thereforetheknowledgeoftheMarkovparameterestimates
fromthelastrowisnotsufficienttorecovertheextendedsubspacematrices,
120
¯
H
NKF
fp
,
¯
G
NKF
and
¯
H
NKF
. Asaresult,oneneedstoestimatef HOARXmodels
ofincreasingorders.
6.3.2 From NKF to System Markov Parameters
SimilartotheSKFparameterization,thesystemMarkovparameters
canbedenotedas
{
g
NKF
k;k
= 0
nynu
;
g
NKF
k;ki
=CA
i1
B; fori> 1,
and
{
h
NKF
k;k
=I
nyny
;
h
NKF
k;ki
=CA
i1
K
ki
; fori> 1.
Now the question is, with NKF parameterization, how to convert the pre-
dictorMarkovparameterstosystemMarkovparameters. Byexaminingthe
structureoftheNKFparameterizedpredictorandsystemMarkovparame-
ters,werevealthegeneralrelationshipasfollows,



g
NKF
k;k1
=g
NKF
k;k1
;
g
NKF
k;ki
=g
NKF
k;ki
+
i1
∑
m=1
h
NKF
k;km
g
NKF
km;ki
; fori> 2,
(6.16)



h
NKF
k;k1
=h
NKF
k;k1
;
h
NKF
k;ki
=h
NKF
k;ki
+
i1
∑
m=1
h
NKF
k;km
h
NKF
km;ki
;
fori> 2,
(6.17)
Remark 1: The relationships in (6.16) and (6.17) enable one to obtain
the NKF parameterized system Markov parameters recursively. When all
the time-varyingK
i
is replacedby the steady Kalman gainK, the relation-
ships in (6.16) and (6.17) are reduced to (6.8) and (6.9). This means the re-
121
lationships in (6.16) and (6.17) hold for both SKF and NKF, hence they are
generalconversionrelationshipsofMarkovparameters.
6.4 Subspace Identification with NKF Parameterization
Inthissection,weshalldevelopaNKFparameterizationbasedsub-
spaceidentificationmethod(abbreviatedasNKFID)tosolvethefinitedata
horizonclosed-loopidentificationproblem.
NKFparameterizationcomplicatesthestructureofthepredictorob-
servability matrix
¯
Γ
NKF
f
and predictor controllability matrix
¯
L
NKF
k;p
, making
the recovery of the system matrices more difficult than the SKF parame-
terization. As an alternative, the process observability matrix Γ
f
, which is
definedas,
Γ
f
=





C
CA
.
.
.
CA
f1





;
preserves the structure in the NKF parameterization. To take advantage of
this, we resort to the innovation form of the extended state space which
contains Γ
f
. Throughiterativesubstitution, theextendedstatespaceinthe
innovationformcanbederivedas
y
f
(k) = Γ
f
¯
L
NKF
k;p
z
p
(k)+G
NKF
u
f
(k)+H
NKF
e
f
(k)
= H
NKF
fp
z
p
(k)+G
NKF
u
f
(k)+H
NKF
e
f
(k) (6.18)
ComposedofthesystemMarkovparameters,G
NKF
andH
NKF
havethefol-
122
lowingforms,
G
NKF
=





0 0 ··· 0
g
NKF
k+1;k
0 ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
g
NKF
k+f1;k
··· g
NKF
k+f1;k+f2
0





;
H
NKF
=





I 0 ··· 0
h
NKF
k+1;k
I ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
h
NKF
k+f1;k
··· h
NKF
k+f1;k+f2
I





:
Analogoustothe
¯
H
NKF
fp
in(6.14),H
NKF
fp
istheproductoftheprocessobserv-
abilitymatrixandthepredictorcontrollabilitymatrix.
The extended state space models (6.14) and (6.18) are both bases for
subspace identification. With many SIMs developed based on them, they
deserve a thorough analysis and comparison. The predictor form in (6.14)
is preferred by most closed-loop identification methods due to its flexibil-
ity in handling input noise correlation. The Kalman filtered noise term in
(6.14) can be easily got rid of through least squares regression. However,
the predictor observability matrix
¯
Γ
NKF
f
in (6.14) is inconvenient for system
matricesrecovery.
NKFIDtakesadvantageofthemodelsinboth(6.14)and(6.18). Based
on the predictor form in (6.14), NKFID estimates the predictor Markov pa-
rameters through the least squares regression on f HOARX models and
builds
¯
G
NKF
accordingly. The relationship in (6.16) allows one to convert
¯
G
NKF
toG
NKF
inwhichthesystemMarkovparametersarecalculatedrecur-
123
sively. Theconversionbridgestheextendedstatespacemodelsin(6.14)and
(6.18)sothatthefuturestepscanbecarriedonthemodelin(6.18).
The knowledge ofG
NKF
enables the removal ofu
f
(k) from the right
handside,resultingintheresidualfutureoutputas
˜ y
f
(k) , y
f
(k)−G
NKF
u
f
(k)
= H
NKF
fp
z
p
(k)+H
NKF
e
f
(k): (6.19)
Thenoisetermcanthenberemovedbymultiplyingbothsidesof(6.19)with
z
p
(k)
T
,
˜ y
f
(k)z
p
(k)
T
=H
NKF
fp
z
p
(k)z
p
(k)
T
:
The H
NKF
fp
estimate can be obtained directly from the project or the CVA
weightingmaybeapplied. Systemmatricescanthenberecoveredfromthe
classicalmethodssuchasERA.
Remark 2: The novelty of NKFID is that it is completely based on
the NKF parameterization, and thus it is a more appropriate solution to
thepracticalfinitedataclosed-loopidentificationproblem. Inaddition,the
conversion relationship in (6.16) givesG
NKF
in NKFID, and the knowledge
ofG
NKF
grantstheinnovationformin(6.18)theabilitytohandleclosed-loop
data.
124
6.5 Simulation Example
In this section, we use a benchmark problem to evaluate the pro-
posedNKFIDmethodandcompareitwiththeOKID,whichisregardedas
a representative SKF parameterization based closed-loop identification ap-
proach. The closed-loop control system is adopted from Verhaegen ([68])
whofirstusedittoevaluatedifferentclosed-loopSIMs. Thesystemmodel
intheinnovationformispresentedasfollowing,
A =






4:40 1 0 0 0
−8:09 0 1 0 0
7:83 0 0 1 0
−4:00 0 0 0 1
0:86 0 0 0 0






;B =






:00098
0:01299
0:01859
0:0033
−0:00002






;
C
T
=






1
0
0
0
0






;K =






2:3
−6:64
7:515
−4:0146
0:86336






;D = 0
The controller used has the state space model composed by following ma-
trices:
A
c
=




2:65 −3:11 1:75 −0:39
1 0 0 0
0 1 0 0
0 0 1 0




;B
c
=




1
0
0
0




;
C
c
=
[
−0:4135 0:8629 −0:7625 0:2521
]
;D
c
= 0:61:
Monte-Carlo simulation is used to guarantee the statistical signifi-
cance: for each experiment 50 data sets are generated, each with 1200 data
125
points. The innovation sequence e is the Gaussian white noise with zero
mean and variance 1=9. Identification accuracy and consistency are two
important aspects in the evaluation. Both the averaged Bode plot and the
scatter plot of the estimated poles will be given to demonstrate the iden-
tification accuracy. In addition, the mean square error (MSE) will also be
calculated. Thepoleestimationvariancewillbeprovidedsinceitindicates
identificationconsistency.
Wefirstsimulatethebasecaseinwhichthepasthorizonpandfuture
horizonf arebothsetas20. Fig.6.1andFig.6.2showthescatteredplotofthe
estimatedpolesofAfromNKFIDandOKID,respectively. Theresultshows
that with a finite horizon p = 20, although both successfully identify the
poles of the system, the identification variance of NKFID is much smaller,
indicatingNKFIDoutperformsOKIDintermsofconsistency.
InFig.6.3theBodeplotshowstheresultcomparisoninthefrequency
domain. Thephaseistheresultaftermod360toremaininbetween[0;360].
ItisobviousthatNKFIDprovidestheidentificationwithbettermagnitude
fittinginlowfrequencyaswellasbetterphasefittinginhighfrequency.
Table6.1andTable6.2summarizetheMSEandestimatevariancefor
eachpole. MSEcanbeinterpretedasadirectreflectionoftheidentification
accuracy. Notsurprisingly,NKFIDprovidestheidentificationwithsmaller
MSE and smaller variance for all the poles. The result shows that NKFID
outperforms SKF parameterization based SIMs when the data horizon is
126
0 0.5 1 1.5
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
Real
Imaginary
Pole Estimation from NKFID
Figure6.1: ScatteredplotoftheestimatedpolesfromNKFID,p =f = 20.
0 0.5 1 1.5
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
Real
Imaginary
Pole Estimation from OKID
Figure6.2: ScatteredplotoftheestimatedpolesfromOKID,p =f = 20.
127
10
−2
10
−1
10
0
−40
−20
0
20
40
60
Bode Plot
Magnitude (dB)

True
NKFID
OKID
10
−2
10
−1
10
0
−400
−300
−200
−100
0
Frequency (rad/s)
Phase (deg)
Figure 6.3: Bode plot of NKFID and OKID. The red solid line is the true
values. Theblueandgreendashedlinesaretheestimatedvaluesaveraged
from50runsbasedonNKFIDandOKID,respectively.
128
Table6.1: MSEforpoleestimation
MSE(×10
4
)
NKFID OKID
1.0000 8.41 108.33
0.9681±0.1486i 7.42 62.18
0.7319±0.6007i 32.17 318.95
Table6.2: Varianceforpoleestimation
Variance(×10
4
)
NKFID OKID
1.0000 0.04 4.56
0.9681±0.1486i 7.19 12.21
0.7319±0.6007i 4.97 119.02
finite.
We vary the value of p to see the effect of the data horizon on the
identificationresult. Fig.6.4andFig.6.5showtheBodemagnitudeplotsand
the scatter plots of the pole for the cases when we keepf = 20 and reduce
pto15, 10and5. Asexpected, theshrinkofthehorizondoesnotaffectthe
NKFIDwhileitdeterioratetheidentificationfromOKIDsignificantly.
With the data horizon p reduced to 15, 10 and 5, Table 6.3 demon-
stratesastraightcomparisonoftheMSEandvariancefromtheestimationof
thepolelocatedat1.0000. Theresultsverifyourargument: astheNKFpa-
129
10
−2
10
−1
10
0
−100
0
100
Mag (dB)
NKFID: p=20
0 0.5 1 1.5
−0.5
0
0.5
Real
Imag
10
−2
10
−1
10
0
−100
0
100
Mag (dB)
NKFID: p=15
0 0.5 1 1.5
−0.5
0
0.5
Real
Imag
10
−2
10
−1
10
0
−100
0
100
Mag (dB)
NKFID: p=10
0 0.5 1 1.5
−0.5
0
0.5
Real
Imag
10
−2
10
−1
10
0
−100
0
100
Mag (dB)
NKFID: p=5
0 0.5 1 1.5
−0.5
0
0.5
Real
Imag
Figure 6.4: Identification result from NKFID with varying p. The left col-
umnistheBodemagnitudeplots;Thesolidlinearethetruevaluesandthe
dashed line are the NKFID estimated values averaged from 50 runs. The
rightcolumnisthescatterplotsoftheeigenvaluesofestimatedA
130
10
−2
10
−1
10
0
−50
0
50
100
Mag (dB)
OKID: p=20
0 0.5 1 1.5
−0.5
0
0.5
Real
Imag
10
−2
10
−1
10
0
−50
0
50
100
Mag (dB)
OKID: p=15
0 0.5 1 1.5
−0.5
0
0.5
Real
Imag
10
−2
10
−1
10
0
−50
0
50
100
Mag (dB)
OKID: p=10
0 0.5 1 1.5
−0.5
0
0.5
Real
Imag
10
−2
10
−1
10
0
−50
0
50
100
Mag (dB)
OKID: p=5
0 0.5 1 1.5
−0.5
0
0.5
Real
Imag
Figure6.5: IdentificationresultfromOKIDwithvaryingp. Theleftcolumn
istheBodemagnitudeplotsfordifferentp;Thesolidlinearethetruevalues
andthedashedlinearetheOKIDestimatedvaluesaveragedfrom50runs.
TherightcolumnisthescatterplotsoftheeigenvaluesofestimatedA
131
Table6.3: MSEandVarianceforvaryingp
MSE(×10
4
) Var(×10
4
)
NKFID OKID NKFID OKID
p=20 8.41 108.33 0.04 4.56
p=15 8.43 157.13 0.04 14.56
p=10 9.17 1013.71 0.04 74.51
p=5 8.11 9069.79 0.05 3890.52
rameterization based method, NKFID is robust to the data horizon length.
With same data horizon, NKFID provides better identification result than
OKIDinbothaccuracyandconsistencyaspects. Ontheotherhand,theSK-
FparameterizationbasedmethodOKIDissensitivetothedatahorizonand
maysufferfrombiaswhenhorizonisdeficient. Theresultisconsistentwith
the Kalman filter theory in the sense that the Kalman gain is time-varying
exceptforsufficientlylonghorizonandtheSKFbasedSIMsaresensitiveto
thelengthofthedatahorizons.
6.6 Summary
Datasetsinpracticalidentificationproblemsalwayshavefinitedata
length. This condition challenges the existing SKF parameterization popu-
larly and conveniently used in most existing SIMs since SKF imposes in-
finite data horizons for Kalman filter to become steady. In this chapter,
132
weproposeanovelsubspaceidentificationalgorithmwithnon-steadystate
Kalmanfilterparameterization. Theproposedmethodprovidesoptimalpa-
rameterizationandsolutionregardlessofthethehorizonlengthsofthepast
and future windows, and it is applicable to closed-loop identification. As
a by-product, the NKF framework can explain why it is often beneficial to
estimateabankofHOARXmodelsinSIMsoverasingleHOARXmodel.
A critical step of the proposed identification method is the conver-
sion from the NKF Markov parameters to the system Markov parameters.
Asanothercontribution, therelationshipestablishedinthischapterallows
one to obtain system Markov parameters from NKF parameters. In addi-
tion,itconnectstheinnovationformwiththepredictorformandenablethe
innovationformtohandleclosed-loopidentification.
BycomparingwiththeOKID,whichisbasedonSKF,theadvantage
of the proposed method in the finite data case are demonstrated through a
benchmarkproblem. TheadvantageofNKFoverSKFbasedmethodsisin-
creasingly noticeable with shrinking horizons which is consistent with the
non-steady state Kalman filter. While it is naturally desirable to compare
the proposed method to other existing SIMs that are applicable to closed
loop identification, the comparison with OKID is a head-to-head compari-
son where the effect of using NKF parameterization is highlighted. Future
workwillfocusonhowtheNKFparametrizationimprovesotherSKFbased
133
subspaceidentificationmethods.
134
Chapter 7
Conclusions
This dissertation presents a novel process and disturbance model-
ing approach and its application to MPC and subspace identification. The
proposedmodelingschemeaddressesthemajordrawbackassociatedwith
interpretingfinitelengthofdatawithSKF.
As a newly emerging and fast developing system identification
method, subspace identification is reviewed first. Three representative
open-loopSIMsarestudiedandtheiradvantagesovertraditionalmethods
are discussed. The main challenge of subspace identification can be
summarizedastheidentificationwithclosed-loopdataoffinitelength.
Toanalyzethefinitehorizonidentificationproblem, theprogressive
parameterization framework is introduced to provide insight of various
stages of subspace identification. It could be concluded that the first stage
of most subspace identification is estimating PMPs from high order ARX
models. It has been found that interpreting PMP estimates with SKF may
leadtosuboptimalsolutionwhendatalengthisfinite.
Based on the equivalence between NKF and least squares, the NKF
parameterization based process and disturbance modeling scheme is pro-
135
posedtohandlefinitedatasystemidentificationproblem. Itparameterizes
the HOARX model using NKF theory, leading to the best state-space mod-
el for finite data. It requires process and disturbance sharing poles. The
proposeddatabasedmodelingschemeisincorporatedtoDMCtoimprove
the disturbance model. Simulation results indicate that the proposed dis-
turbancemodelcanreducethepredictionerror.
The method is further extended to a more general NKF based mod-
elingmethodwhichdoesnotrequiresharedpoles. Thegeneraltransforma-
tion relation between PMPs and SMPs is developed. The estimated PMP-
s can be converted to SMPs under NKF structure hence the obtained FIR
modelisconsistentwithfinitedatalength. Theproposedmethodisalsoin-
corporated to DMC, and the benefits are verified by the simulation results
ofWood-BerrydistillationcolumnandTEP.
The practical system identification problem always deals with finite
data. BasedontheNKFparameterization,NKFIDisdevelopedasasuitable
SIM for the data of finite horizon. NKFID employs the NKF parameteriza-
tioninthefirststage. Thefutureinputisremovedfromoutputsothatitis
unbiasedundertheclosed-loopcondition. Theproposedalgorithmiscom-
pared to OKID to show its significant advantages in the limited horizon
cases.
Future directions of research could be further exploration of the es-
timation of B and D. The derivation of the stochastic characteristics with
136
NKF based SIM could be an interesting research topic. How to incorporat-
ingtheidentifiedstatespacemodeltoindustrialMPCalsodeservesfurther
investigation.
137
Bibliography
[1] H. Akaike, “A new look at the statistical model identification,” IEEE
Trans.Auto.Cont.,vol.19,pp.716–723,1974.
[2] D. Bauer, M. Deistler, and W. Scherrer, “Consistency and asymptotic
normality of some subspace algorithms for systems without observed
inputs,”Automatica,vol.35,pp.1243–1254,1999.
[3] D.BauerandM.Jansson,“Analysisoftheasymptoticpropertiesofthe
MOESPtypeofsubspacealgorithms,”Automatica,vol.36,pp.497–509,
2000.
[4] D. Bauer and L. Ljung, “Some facts about the choice of the weighting
matricesinLarimoretypeofsubspacealgorithms,”Automatica,vol.38,
pp.763–773,2002.
[5] M. Bauer and I. K. Craig, “Economic assessment of advanced process
control - a survey and framework,” Journal of Process Control, vol. 18,
pp.2–18,2008.
[6] G. E. P. Box, G. M. Jenkins, and G. Reinsel, Time Series Analysis: Fore-
castingandControl. EnglewoodCliffs,NJ:Prentice-Hall,1994.
[7] A.ChiusoandG.Picci,“Consistencyanalysisofsomeclosed-loopsub-
spaceidentificationmethods,” Automatica,vol.41,no.3,pp.377–391,
2005.
[8] C.T.ChouandM.Verhaegen,“Subspacealgorithmsfortheidentifica-
tionofmultivariabledynamicerrors-in-variablesmodels,”Automatica,
vol.33,no.10,pp.1857–1869,1997.
[9] C.R.CutlerandR.B.Hawkins,“Applicationofalargepredictivemul-
tivariablecontrollertoahydrocrachersecondstagereactor,”inProceed-
ingsofACC,1988,pp.284–291.
138
[10] C.R.CutlerandB.L.Ramaker,“Dynamicmatrixcontrol–acomputer
controlalgorithm,”inAIChENationalMeeting,1979.
[11] M. Deistler and K. Peternell, “Consistency and relative efficiency of
subspacemethods,”Automatica,vol.31,pp.1865–1875,1995.
[12] B. DeMoor, J. Vandewalle, L. Vandenberghe, and P. VanMieghem, “A
geometrical strategy for the identification of state space models of
linear multivariable systems with singular value decomposition,” in
Proc.8th IFAC Symp. on Identification and System Parameter Estimation,
1988.
[13] J. J. Downs and E. F. Vogel, “A plant-wide industrial process control
problem,” Computers & Chemical Engineering, vol. 17, no. 3, pp. 245–
255,1993.
[14] B. A. Francis and W. M. Wonham, “The internal model principle of
controltheory,”Automatica,vol.12,no.5,pp.457–465,1976.
[15] C. E. Garcia and M. Morari, “Internal model control. 2. design pro-
cedure for multivariable systems,” Industrial & Engineering Chemistry
ProcessDesignandDevelopment,vol.24,no.2,pp.472–484,1985.
[16] C. E. Garcia and M. Morari, “Internal model control. 3. multivariable
controllawcomputationandtuningguidelines,”Industrial&Engineer-
ing Chemistry Process Design and Development, vol. 24, no. 2, pp. 484 –
494,1985.
[17] C.E.GarciaandM.Morari,“Internalmodelcontrol.aunifyingreview
and some new results,” Industrial & Engineering Chemistry Process De-
signandDevelopment,vol.21,no.2,pp.308–323,1985.
[18] C. E. Garcia and A. M. Morshedi, “Quadratic programming solution
ofdynamicmatrixcontrol(QDMC),”ChemicalEngineeringCommunica-
tions,vol.46,no.1,pp.73–87,1986.
[19] C. E. Garcia, D. M. Prett, and M. Morari, “Model predictive control:
Theoryandpractice–asurvey,” Automatica, vol.25, no.3, pp.335–348,
1989.
139
[20] M.Gevers,“Apersonalviewonthedevelopmentofsystemidentifica-
tion,”in Proceedings of the 13th IFAC symposium on System Identification,
Rotterdam,Netherlands,2003,pp.773–784.
[21] G. C. Goodwin and K. S. Sin, Adaptive Filtering Prediction and Control.
EnglewoodCliffs,NewJersey: Prentice-Hall,1984.
[22] T. Gustafsson, “Subspace-based system identification: weighting and
pre-filtering of instruments,” Automatica, vol. 38, no. 3, pp. 433–443,
2002.
[23] E.T.HaleandS.J.Qin,“Subspacemodelpredictivecontrolandacase
study,”inProceedingsofACC,vol.6,2002,pp.4758–4763.
[24] C.A.HarrisonandS.J.Qin,“Discriminatingbetweendisturbanceand
processmodelmismatchinmodelpredictivecontrol,”JournalofProcess
Control,vol.19,pp.1610–1616,2009.
[25] B. L. Ho and R. E. Kalman, “Effective construction of linear
state-variable models from input-output functions,” Regelungstechnik,
vol.12,pp.545–548,1965.
[26] H. Hotelling, “Relations between two sets of variants,” Biometrika,
vol.28,pp.321–377,1936.
[27] B. Huang and R. Kadali, Dynamic Modeling, Predictive Control and Per-
formanceMonitoring. Springer,2008.
[28] B. Huang, S. X. Ding, and S. J. Qin, “Closed-loop subspace identifi-
cation: an orthogonal projection approach,” Journal of Process Control,
vol.15,pp.53–66,2005.
[29] M. Jansson, “Subspace identification and ARX modelling,” in Proceed-
ingsofthe13thIFACSYSIDSymposium,Rotterdam,NL,Aug2003.
[30] M. Jansson, “A new subspace identification method for open and
closed loop data,” in 16th IFAC World Congress, Prague, Czech Repub-
lic,2005.
140
[31] M.JanssonandB.Wahlberg,“Onconsistencyofsubspacemethodsfor
systemidentification,”Automatica,vol.34,no.12,pp.1507–1519,1998.
[32] J.N.Juang,J.E.Cooper,andJ.R.Wright,“Aneigensystemrealization
algorithm using data correlations (ERA/DC) for modal parameter i-
dentification,”ControlTheoryandAdvancedTechnology,vol.4,no.1,pp.
5–14,1988.
[33] J. N. Juang and R. S. Pappa, “An eigensystem realization algorithm
for modal parameter identification and modal reduction,” Journal of
Guidance,ControlandDynamics,vol.8,no.5,pp.620–627,1985.
[34] J. N. Juang, M. Phan, L. G. Horta, and R. W. Longman, “Identifica-
tion of Observer/Kalman filter Markov parameters: Theory and ex-
periments,”JournalofGuidance,ControlandDynamics,vol.16,no.2,pp.
320–329,1993.
[35] J. N. Juang and M. Q. Phan, “Identification of system, observer, and
controller from closed-loop experimental data,” in Proceedings of the
AIAA Guidance, Navigation, and Control Conference, Hilton Head, SC,
1992.
[36] M.KanoandM.Ogawa,“Thestateoftheartinchemicalprocesscon-
trol in Japan: Good practice and questionnaire survey,” Journal of Pro-
cessControl,vol.20,no.9,pp.969–982,2010.
[37] T. Katayama, Subspace methods for system identification. Berlin, Ger-
many: Springer,2005.
[38] T.Knudsen,“Consistencyanalysisofsubspaceidentificationmethods
basedonalinearregressionapproach,” Automatica, vol.37, pp.81–89,
2001.
[39] W. E. Larimore, “Canonical variate analysis in identification, filtering,
and adaptive control,” in Proceedings of the 29th IEEE Conference onDe-
cisionandControl,Dec.1990,pp.596–604vol.2.
[40] W. E. Larimore, “Statistical optimality and canonical variate analysis
systemidentification,”SignalProcessing,vol.52,pp.131–144,1996.
141
[41] W. E. Larimore, “Large sample efficiency for ADAPTX subspace sys-
temidentificationwithunknownfeedback,”inProc.IFACDYCOPS’04,
2004.
[42] W. E. Larimore and J. Baillieul, “Identification and filtering of nonlin-
earsystemsusingcanonicalvariateanalysis,”in Proceedingsof the 29th
IEEEConf.onDecisionandControl,Honolulu,HI,1990,pp.635–640.
[43] J. H. Lee, M. Morari, and C. E. Garcia, “State-space interpretation of
modelpredictivecontrol,”Automatica,vol.30,no.4,pp.707–717,1994.
[44] L. Ljung, System identification: Theory for the user. Englewood Cliffs,
NewJersey: Prentice-Hall,1999.
[45] L.LjungandT.McKelvey,“Subspaceidentification,”SignalProcessing,
vol.52,pp.209–215,1996.
[46] P. R. Lyman and C. Georgakis, “Plant-wide control of the Tennessee
Eastman problem,” Computers & Chemical Engineering, vol. 19, no. 3,
pp.321–331,1995.
[47] T.J.McAvoyandN.Ye,“BasecontrolfortheTennesseeEastmanprob-
lem,” Computers & Chemical Engineering, vol. 18, no. 5, pp. 383–413,
1994.
[48] M.MorariandG.Stephanopoulos,“Studiesinthesynthesisofcontrol
structures for chemical processes: Part III: Optimal selection of sec-
ondarymeasurementswithintheframeworkofstateestimationinthe
presence of persistent unknown disturbances,” AIChE Journal, vol. 26,
no.3,pp.247–260,1980.
[49] K.R.MuskeandT.A.Badgwell,“Disturbancemodelingforoffset-free
linearmodelpredictivecontrol,”JournalofProcessControl,vol.12,no.5,
pp.617–632,2002.
[50] A. Negiz and A. Cinar, “PLS, balanced and canonical variate realiza-
tion techniques for identifying varma models in state space,” Chemo-
metricsandIntelligentLaboratorySystems,vol.38,pp.209–221,1997.
142
[51] G. Pannocchia and J. B. Rawlings, “Disturbance models for offset-free
model-predictive control,” AIChE Journal, vol. 49, no. 2, pp. 426 – 437,
2003.
[52] K.Peternell,W.Scherrer,andM.Deistler,“Statisticalanalysisofnovel
subspace identification methods,” Signal Processing, vol. 52, pp. 161–
177,1996.
[53] S.J.Qin,“Anoverviewofsubspaceidentification,” Computer & Chem-
icalEngineering,vol.17,no.10-12,pp.1502–1513,2006.
[54] S.J.QinandL.Ljung,“Closed-loopsubspaceidentificationwithinno-
vation estimation,” in Proceedings of the 13th IFAC SYSID Symposium,
Rotterdam,NL,Aug2003,pp.887–892.
[55] S.J.QinandL.Ljung, “ParallelQRimplementationofsubspaceiden-
tification with parsimonious models,” in Proceedings of the 13th IFAC
SYSIDSymposium,Rotterdam,NL,Aug2003,pp.1631–1636.
[56] S.J.QinandT.A.Badgwell, “Asurveyofindustrialmodelpredictive
controltechnology,”ControlEngineeringPractice,vol.11,no.7,pp.733–
764,2003.
[57] S. J. Qin, Y. Zhao, Z. Sun, and T. Yuan, “Progressive parametrization
in subspace identification models with finite horizons,” in 49th IEEE
Conference on Decision and Control (CDC), 2010, dec. 2010, pp. 2819 –
2824.
[58] S.QinandL.Ljung,“Ontheroleoffuturehorizoninclosed-loopsub-
spaceidentification,” in 14th IFACProc. 51th IEEE Conf. on Decision and
Control,Newcastle,Australia,2006,pp.1080–1084.
[59] M. R. Rajamani, J. B. Rawlings, and S. J. Qin, “Achieving state esti-
mation equivalence for misassigned disturbances in offset-free model
predictivecontrol,”AIChEJournal,vol.55,no.2,pp.396–407,2009.
[60] N. L. Ricker, “Decentralized control of the Tennessee Eastman chal-
lengeprocess,”JournalofProcessControl,vol.6,no.4,pp.205–221,1996.
143
[61] R. Shi and J. F. MacGregor, “A framework for subspace identification
methods,”inProceedingsofACC,2001,pp.3678–3683.
[62] T. S. Soderstrom and P. G. Stoic, System Identification. Englewood
Cliffs,NJ:Prentice-Hall,1989.
[63] H. W. Sorenson, “Least-squares estimation: from Gauss to Kalman,”
Spectrum,IEEE,vol.7,no.7,pp.63–68,July1970.
[64] Z.Sun,Y.Zhao,andS.J.Qin,“ImprovingindustrialMPCperformance
with data-driven disturbance modeling,” in 50th IEEE Conference on
DecisionandControl(CDC),2011,dec.2011,pp.2819–2824.
[65] P.VanOverscheeandB.DeMoor,“N4SID:Subspacealgorithmsforthe
identificationofcombineddeterministic-stochasticsystems,”Automat-
ica, vol. 30, no. 1, pp. 75 – 93, 1994, special issue on statistical signal
processingandcontrol.
[66] P. VanOverschee and B. DeMoor, “A unifying theorem for three sub-
spacesystemidentificationalgorithms,”Automatica,vol.31,no.12,pp.
1853–1864,1995.
[67] P.VanOverscheeandB.DeMoor, “Closed-loopsubspacesystemiden-
tification,” in Proceedings of the 36th Conference on Decision and Control,
SanDiego,CA,Dec1997.
[68] M. Verhaegen, “Application of a subspace model identification tech-
nique to identify LTI systems operating in closed-loop,” Automatica,
vol.29,pp.1027–1040,1993.
[69] M. Verhaegen, “Identification of the deterministic part of MIMO state
spacemodelsgivenininnovationsformfrominput-outputdata,” Au-
tomatica,vol.30,pp.61–74,1994.
[70] 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,pp.1187–1210,1992.
144
[71] M. Viberg, “Subspace methods in system identification,” in Proc. 10th
IFACSymp.onIdentificationandSystemParameterEstimation,1994.
[72] M. Viberg, “Subspace-based methods for the identification of linear
time-invariantsystem,”Automatica,vol.31,no.12,pp.1835–1851,1995.
[73] J. Wang and S. J. Qin, “A new subspace identification approach based
onprincipalcomponentanalysis,”JournalofProcessControl,vol.12,pp.
841–855,2002.
[74] R.K.WoodandM.W.Berry,“Terminalcompositioncontrolofabinary
distillation column,” Chemical Engineering Science, vol. 28, no. 9, pp.
1707–1717,1973.
[75] Z.Xu,Y.Zhu,K.Han,J.Zhao,andJ.Qian,“Amulti-iterationpseudo-
linear regression method and an adaptive disturbance model for M-
PC,”JournalofProcessControl,vol.20,no.4,pp.384–395,2010.
[76] Z. Zang and R. R. Bitmead, “Iterative weighted least-squares identifi-
cation and weighted LQG control design,” Automatica, vol. 31, no. 11,
pp.1577–1594,1995.
[77] Y. Zhao, Z. Sun, and S. J. Qin, “Non-stationary Kalman filter parame-
terization with applications to MPC,” in Proc. 51th IEEE Conf. on Deci-
sionandControl,2012.

Abstract (if available)

Abstract Subspace identification methods (SIMs) have drawn tremendous research interest and have been widely recognized for its computational efficiency and strong capacity in identifying MIMO systems. Most existing SIMs are based on steady‐state Kalman filter (SKF) models which imply infinite horizons to be used in forming the dynamic time lags for the model. However, in practice, finite data length is the only option that dictates even shorter time horizons for the number of time lags. As a result, the SKF parameterization for SIMs can lead to suboptimal solution with finite number of data. ❧ In this dissertation, a novel nonsteady‐state Kalman filter (NKF) parameterization based SIM (NKFID) is proposed to address two aspects: finite data horizon and closed-loop identification. A progressive parameterization framework is first built to interpret the models in each stage of SIM. The model of the first stage is shown as non‐parametric high order ARX model, which has applied value in model predictive control (MPC). ❧ In practice, MPC usually assumes a step‐like disturbance model which is insufficient to describe complex dynamics. The finite impulse response (FIR) model of the disturbance is recommended in this work. Two data based parameterization methods are proposed to parameterize the FIR form disturbance model. The disturbance modeling schemes with SKF parameterization is first developed to improve MPC performance. Considering finite data length, the NKF parameterization based modeling scheme is proposed and improved, and its suitability for finite data is verified theoretically. The general conversion relation between the predictor Markov parameters (PMPs) and the system Markov parameters (SMPs) is developed under the NKF framework, based on which the disturbance model is derived and incorporated to MPC. ❧ The NKFID is proposed as a novel SIM under the NKF framework. The algorithm uses the NKF parameterization hence optimal solution is guaranteed even with limited data horizons. A routinely used projection that removes the impact of the future input is avoided to make NKFID applicable to closed‐loop identification. The performance of the proposed disturbance modeling and subspace identification methods are demonstrated by several simulation examples and the benefits are verified by the comparison with existing methods.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Performance monitoring and disturbance adaptation for model predictive control

PDF

Economic model predictive control for building energy systems

PDF

Performance prediction, state estimation and production optimization of a landfill

PDF

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

PDF

Process data analytics and monitoring based on causality analysis techniques

PDF

High-accuracy adaptive vibrational control for uncertain systems

PDF

Active state tracking in heterogeneous sensor networks

PDF

Data-driven H∞ loop-shaping controller design and stability of switched nonlinear feedback systems with average time-variation rate

PDF

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

PDF

Provable reinforcement learning for constrained and multi-agent control systems

PDF

Digital signal processing techniques for music structure analysis

PDF

Structure learning for manifolds and multivariate time series

PDF

Reinforcement learning with generative model for non-parametric MDPs

PDF

Elements of robustness and optimal control for infrastructure networks

PDF

Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures

PDF

Damage detection using substructure identification

PDF

Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction

PDF

Representation, classification and information fusion for robust and efficient multimodal human states recognition

PDF

Uncertainty quantification and data assimilation via transform process for strongly nonlinear problems

PDF

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

Asset Metadata

Creator Zhao, Yu (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 02/19/2014

Defense Date 01/24/2014

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

Tag ARX model,disturbance modeling,Markov parameter,non‐steady state Kalman filter,OAI-PMH Harvest,predictive control,subspace identification

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Qin, S. Joe (committee chair), Safonov, Michael G. (committee member), Shing, Katherine (committee member)

Creator Email yu.zhao1985@gmail.com,yuzhao@usc.edu

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

Unique identifier UC11297245

Identifier etd-ZhaoYu-2254.pdf (filename),usctheses-c3-363949 (legacy record id)

Legacy Identifier etd-ZhaoYu-2254.pdf

Dmrecord 363949

Document Type Dissertation

Format application/pdf (imt)

Rights Zhao, Yu

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