n
a
e
d
a
ila
IM
d
regression through an over-parameterisation technique. Then, a state-space realisation of the system is retrieved
IET
do
www.ietdl.orgusing a recursive subspace identification method. Simulation results are provided to demonstrate the
effectiveness of the algorithm.
1 Introduction
System identificationbased on input–output datameasurements
is a necessary step in most applications such as control design,
plant diagnosis and monitoring. A number of methods have
been developed for the identification of linear models in
input–output form as well as in state-space representation.
Some of the proposed techniques have been also extended to
slowly time-varying systems by means of recursive
implementations (see, e.g. [1–4] for references on recursive
methods for state-space models). The identification of models
for non-linear non-stationary systems, however, remains a very
challenging research area. To that purpose the identification of
Hammerstein models is considered in this paper. Such models
consist of a non-linear memoryless function f (.) followed by a
linear dynamic system. The difficulty in separating the non-
linear part from the linear one is related to the fact that their
respective terms are present only in cross-product forms in the
system equations. A typical way to proceed is to treat each of
these cross-product terms as an unknown parameter. This
procedure, which results in an increased number of unknowns,
is usually referred to as the over-parameterisation method
[5, 6]. Generally, the static non-linearity f (.) is expanded on a
basis of known elementary functions (radial basis functions
(RBFs), polynomials, B-splines, . . .) and then, the coordinates
associated to it are estimated [7–9].
Finally, while most of the proposed methods are based on
input–output models, it is interesting to point out that state-
space models are often more convenient for describing
systems, particularly in the multivariable case, and are also a
suitable setting for the formulation of the Hammerstein
identification problem (see, e.g. [10], where the use of
subspace model identification for Hammerstein systems was
first proposed). More recently, in [11] the efficiency of
component-wise least squares support vector machines (LS-
SVM) combined with subspace methods for Hammerstein
state-space model identification was first demonstrated. Note
that the proposed approach is very similar in the results to the
well-known technique of expanding the non-linearity on
gaussian functions basis. In order to take into account also the
case of slowly time-varying systems, for which an on-line
update of the estimate of the model parameters would be
useful, in this paper we extend the work presented in [11] to
time-varying systems by using LS-SVM to recursively
estimate the non-linear part of the system and ordinary least
squares for recovering the linear part in state-space form.
The paper is organised as follows. In Section 2 we state
the main assumptions about the system we are concerned
with. In Section 3 an LS-SVM regression technique is
considered to estimate the Markov parameters and the
non-linear function f (.). In Section 4 we present thePublished in IET Control Theory and Applications
Received on 7th August 2008
Revised on 12th January 2009
doi: 10.1049/iet-cta.2008.0339
Recursive subspace ide
Hammerstein models b
squares support vector
L. Bako1 G. Merce`re2 S. Leco
1De´partement Informatique et Automatique, Ecole des Mines
2Universite´ de Poitiers, Laboratoire d’Automatique et d’Inform
3Dipartimento di Elettronica e Informazione, Politecnico di M
E-mail: lovera@elet.polimi.it
Abstract: A recursive scheme for the identification of S
scheme, first the Markov parameters of the system areControl Theory Appl., 2009, Vol. 3, Iss. 9, pp. 1209–1216
i: 10.1049/iet-cta.2008.0339ISSN 1751-8644
tification of
sed on least
machines
uche1 M. Lovera3
e Douai, France
tique Industrielle, Poitiers, France
no, Italy
O Hammerstein models is presented. In the proposed
etermined, by a least squares support vector machines1209
& The Institution of Engineering and Technology 2009

12
&
www.ietdl.orgrecursive version of the method making use of an adapted
subspace fitting technique to estimate first the observability
matrix and then a state-space realisation of the system.
2 Problem statement
Consider a discrete-time SIMO Hammerstein state-space
system represented as follows
xtþ1 ¼ Axt þ Bf (ut)þ wt
yt ¼ Cxt þDf (ut)þ vt ,
(1)
where yt [ Rny , ut [ R and xt [ Rnx are the output, input
and state vectors, respectively, and wt [ Rnx , vt [ Rny are
white noise vector sequences. Moreover, the system under
consideration is assumed to be stable and of known order
nx. f (.) is a scalar-valued smooth non-linear function.
In view of these definitions, the considered system
identification problem can be formulated as follows: given N
consecutive input–output data measurements ut , yt
 N
t¼1,
derive an algorithm to recursively update estimates of the
state-space matrices (A, B, C , D) and of the static non-
linearity f of the system in (1) up to a similarity transformation.
To deal with deterministic and stochastic signals in a
compact manner, the following operator is defined
E ½  ¼ lim
N!1
1
N
XN
t¼1
E ½  (2)
where E is the expectation operator. For two signals at and bt ,
the cross-covariance matrix will be denoted as Rab ¼ E[atbTt ]
whereas estimates of signal correlations based on
measurements up to time t will be denoted by R^ab,t .
3 LS-SVM identification of
Hammerstein models
Let us write yt as a d-step ahead predictive model using (1).
The obtained equation is very similar to a finite impulse
response (FIR) model
yt ¼ CAdxtd þ
Xd
i¼0
Hif (uti)þ
Xd
i¼1
CAi1wti þ vt
|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
et
(3)
8t . d , where d is a positive integer, et is the equation error
and Hi, defined as
Hi ¼
D if i ¼ 0
CAi1B if i  1

(4)
are the so-called Markov parameters and represent the
impulse response of the system. These parameters do not
depend on the basis of the state-space representation in (1).10
The Institution of Engineering and Technology 2009Equation (3) amounts to
yt(s) ¼ CsAdxtd þ
Xd
i¼0
Hi(s)f (uti)þ et(s) (5)
s ¼ 1, . . . , ny, where yt(s) is the s-th output, Cs is the s-th
row of the observation matrix C and Hi(s) is the s-th row
of Hi .
In [11], a method for off-line identification of Hammerstein
models using LS-SVM regression techniques was proposed:
LS-SVM formulations are in terms of primal (constrained
optimisation problem) and dual problems (problem in
Lagrange multipliers) where feature maps are used in the
primal and positive definite kernels in the dual; examples of
positive definite kernels are linear, polynomial, RBF kernels.
This method produces reliable results but is computationally
heavy since it requires costly singular value decomposition
(SVD) steps (for computing the states sequences and for
extracting the non-linearity). Here, we would like to
investigate a way of avoiding the SVD computation in this
method for online application purposes. Although it is
feasible to estimate a fully parameterised MIMO model by
applying twice the LS-SVM regression procedure as in [11],
only the case of SIMO models is considered in this paper for
the sake of simplicity. More precisely, the proposed approach
can be extended to cover the case of multi-input models; the
only requirement concerns the non-linearity f (.), which has to
be a scalar-valued function.
We will first estimate the terms Hi f () and then the
extended observability matrix. In (5), we choose the integer
d at least larger than the system order and sufficiently large
so that the term CAdxtd may be considered negligible
with respect to the noise term. This approximation (which
is common in many subspace-based approaches, see [12]) is
all the more acceptable as the system has rapid dynamics or
d is relatively large. Then, (5) becomes
yt(s) 
Xd
i¼0
Hi(s)f (uti)þ et(s) (6)
for s ¼ 1, . . . , ny. In the case that the system does not
present rapid dynamics it may be safer to use the observer-
augmenting procedure described in [13] instead. This latter
algorithm introduces an observer gain to set the poles of
the resulting system as close to zero as possible so that the
accuracy of the previous equation is noticeably improved
even for small values of d.
Recall that neither the parameters Hi nor the nonlinear
function f (.) are known. The LS-SVM regression procedure
developed by Suykens [14] consists in modelling the non-
linearity as f () ¼ v`w()þ d where w : Rnu ! Rnh is a
non-linear mapping from the input space to a possibly high
dimensional feature space Rnh , d is a bias compensator
and v [ Rnh . However, direct application of this ideaIET Control Theory Appl., 2009, Vol. 3, Iss. 9, pp. 1209–1216
doi: 10.1049/iet-cta.2008.0339

IET
do
www.ietdl.orgresults here in a hard non-convex optimisation problem [6].
Fortunately, this difficulty can be circumvented by
combining over-parameterisation and component-wise LS-
SVM [15]. Letting
gi,s() ¼ Hi(s)f () ¼ v`i,sw()þ di,s (7)
we obtain from (6)
yt(s) ¼
Xd
i¼0
gi,s(uti)þ et(s) ¼
Xd
i¼0
v`i,sw(uti)þ ds þ et(s)
(8)
with ds ¼
Pd
i¼0 di,s, s ¼ 1, . . . , ny, t ¼ d þ 1, . . . , N . Since
only input–output data (ut , yt) are available, we would like
to compute the non-linear components gi,s from the sample
values of their sum along i. That corresponds to the so-
called component-wise LS-SVM problem. One can then
notice that the solution of (8) is not unique. In fact, for any
solution {gi,s()}
d
i¼0 of (8), {gi,s()þ 1i,s()}
d
i¼0 is also solution
for all arbitrary functions {1i,s()}
d
i¼0, satisfying
Pd
i¼0
1i,s(uti) ¼ 0, 8t ¼ d þ 1, . . . , N . Hence, each function
gi,s can only be estimated up to an additive undetermined
function term. A typical way to overcome that difficulty is to
introduce some additional constraint on the v`i,sw() to be
obtained by solving (8) so that the redundancy is removed.
In [11], this constraint is chosen as a centring relation over
the time window of the measurement points. An alternative
way of dealing with this problem, which is more practical in
view of a recursive implementation, will be presented in the
following.
The LS-SVM constrained optimisation problem is based
on the cost function
J wi,s, et(s)
 
¼ 1
2
Xny
s¼1
Xd
i¼0
v`i,svi,s þ
g
2
Xny
s¼1
XN
t¼dþ1
e2t (s) (9)
subject to (8), where g is a regularisation parameter and serves
as a weight in the cost function to be minimised. In order to
apply the Lagrange multipliers method, we formulate the
associated Lagrangian as
L(vi,s, ds, at,s, et,s) ¼ J (vi,s, et,s)
Xny
s¼1
XN
t¼dþ1
 at,s
Xd
i¼0
v`i,sw(uti)þ :ds þ et(s) yt(s)
( )
where at,s stands for the Lagrange multipliers. The associated
Karush–Kuhn–Tucker optimality conditions are given by
@L
@vi,s
¼ 0! vi,s ¼
XN
t¼dþ1
at,sw(uti),
8i ¼ 0, . . . , d 8s ¼ 1, . . . , ny (10)Control Theory Appl., 2009, Vol. 3, Iss. 9, pp. 1209–1216
i: 10.1049/iet-cta.2008.0339@L
@ds
¼ 0!
XN
t¼dþ1
at,s ¼ 0, 8s ¼ 1, . . . , ny (11)
@L
@at,s
¼ 0, 8t ¼ d þ 1, . . . ,N , 8s ¼ 1, . . . , ny ! yt(s)
¼
Xd
i¼0
v`i,sw(uti)þ ds þ et(s) (12)
@L
@et(s)
¼ 0! et(s) ¼ g1at,s, 8t ¼ d þ 1, . . . ,N
8s ¼ 1, . . . , ny (13)
Taking into account (10) and (13) in (12) gives
yt(s) ¼
Xd
i¼0
XN
q¼dþ1
aq,sw
`(uqi)w(uti)
0
@
1
Aþ ds þ g1at,s
(14)
The kernel function K (uq, ut) is then defined fromRnu  Rnu
to R as K (uq, ut) ¼ w`(uq)w(ut) and selected here to be a
gaussian RBF K (uq, ut) ¼ exp (kuq  utk22=s
2), with s a
constant. Equations (11) and (14) can be summarised in the
following linear equation with as unknowns the dual
parameters of the optimisation problem
0
1M
1`M
Vþ g1IM



ds
as


¼ 0Ys


8s ¼ 1, . . . , ny (15)
where M ¼ N  d , IM is the identity matrix of order M and
V ¼ Vij
 
[ RMM , 1M , as and Ys [ RM are defined as
Vij ¼
Xd
m¼0
K (uiþdm, ujþdm)
1M ¼ 1 . . . 1
`
as ¼ adþ1,s . . . aN ,s
`
Ys ¼ ydþ1(s) . . . yN (s)
`
The solution of (15) is obtained as
ds ¼ 1`MT 1 Ys=(1`MT 11M )
as ¼ T 1 Ys  1Mds
 
, s ¼ 1, . . . , ny
(16)
where T ¼ Vþ g1IM . Note that matrix T is the same for all
the components of the output vector as it depends only on the
input u through the kernel function K. Obviously, the
properties of T (e.g. its conditioning) are directly related to
those of the input.1211
& The Institution of Engineering and Technology 2009

12
&
www.ietdl.orgOnce the dual parameters at,s and ds are known, the
searched non-linear components from (7) can be expressed as
v^`i,sw() ¼
XN
t¼dþ1
at,sK (uti,  ) ¼ a`s ki() (17)
s ¼ 1, . . . , ny with
ki() ¼ K (udþ1i,  ) . . . K (uNi ,  )
`
It should be noticed that, depending on the selected kernel,
the LS-SVM solution to the non-linear regression problem
is very similar to that obtained using the standard
techniques of expanding the non-linearity on polynomials
or RBF bases. Solution (17) is known to be non-sparse,
that is, some of the terms at,sK (uti,  ) have too small a
contribution in the general error reduction and can
therefore be removed without a great impact. Some
pruning techniques have been proposed to render it sparse
(see, e.g. [16]) but our presentation does not address this
problem.
It is known (see, e.g. [11]) that the over-parameterisation
technique may result in an estimate for gi,s which does not
satisfy relation (7) any longer. This means that the
estimates for the components gi,s may not be collinear as
they are expected to be. More precisely, the estimate may
be corrupted by an unknown additive component 1i,s()
such that
Pd
i¼0 1i,s() ¼ 0, 8 t ¼ d þ 1, . . . ,N . Hence, we
have
g^i,s ¼ v^
`
i,sw()þ d^ i,s ¼ Hi(s)f ()þ 1i,s()
Summing the above equation along the subscript i results in
Xd
i¼0
v^`i,sw()þ d^ s ¼ ms f ()
wherems ¼
P
i Hi(s) is assumed to be non-zero. Note that for
any non-zero scalar l, Hi(s)f () ¼ Hi(s)l
 
1
l f ()
 
. In other
words, f (.) can only be estimated up to a scalar factor just as
the matrices (A, B, C, D) can be determined only up to a
similarity transformation. In the light of this remark, the
static non-linearity f can arbitrarily set to be f^ () ms f ()
f^ () ¼
Xd
i¼0
v^`i,sw()þ d^ s, s ¼ 1, . . . , ny (18)
Since the non-linearity f (.) does not depend on a particular
output, it is not worth repeating this estimation for every
output as that may increase the computational cost. The
non-linear function f (.) can be extracted from the
measurements of the input and the ones of a single arbitrary
output s or a combination of all the outputs.12
The Institution of Engineering and Technology 2009From now on, given that the non-linearity is known, the
system matrices can be computed using any subspace
identification method. The linear system to be considered
has input f (u) and output given by the entire output vector
y (instead of one component as above). We propose to
estimate first the Markov parameters Hi from the FIR
model (5) using least squares regression and then compute
the system matrices. The regression equation is then given by
yt ¼ H0    Hd

f (ut)
..
.
f (utd )
2
664
3
775þ et , t ¼ d þ 1, . . . , N
(19)
When this step is completed, all the Hi terms, that is, D, CB,
CAB, . . . , CAd1B and f are known. In the case of a batch
implementation, the system matrices can be directly
computed from these parameters using an SVD; in this
paper however we are interested in developing a recursive
scheme which avoids this step, which is known to be
computationally demanding.
4 Recursive implementation
In the following, a recursive version of the identification
algorithm presented in the previous section will be derived.
More precisely, we assume that not only the linear part of
the system but also the non-linearity f (.) may undergo
some changes (both in its parameters as in its form) during
time. It is then necessary to update its estimate. The
proposed approach to the recursive implementation can be
summarised as follows. Whenever a new input–output data
pair (ut , yt) is measured, the following steps are carried out:
1. A recursive update of the parameters a and d is
performed. This can be obtained by working out a recursive
form for (16).
2. The non-linear function is computed, using the updated
a and d.
3. Once the static non-linearity is known, it is possible
to update the estimate of the Markov parameters Hi,
i ¼ 0, . . . , d of the linear part from (19), using recursive
least squares.
4. Finally, the estimates A^, B^, C^, D^ of the system matrices
for the linear part are updated, using ideas from linear
recursive subspace model identification, such as the
propagator method.
In Section 4.1 a recursive approach to the estimation of the
Markov parameters and of the non-linearity (i.e. Steps 1–3
in the algorithm) will be presented, whereas the recursive
update of the system matrices (i.e. Steps 4–5) will be
discussed in Section 4.2.IET Control Theory Appl., 2009, Vol. 3, Iss. 9, pp. 1209–1216
doi: 10.1049/iet-cta.2008.0339

IE
do
www.ietdl.org4.1 Recursive estimation of the
non-linear function and of the
Markov parameters
The first step to be performed in order to derive a complete
recursive version of the algorithm is to derive an explicit
recursive update for (16). Although most recursive
identification algorithms are developed by using a forgetting
factor, in the present case we opt for a sliding window
approach (note that a similar approach has been used in [17]
to implement an online LS-SVM-based classifier). For,
assume that at a given instant N observations {ut , yt}
kþN1
t¼k
are available. When a new pair of data is acquired, the
window is slid to be {ut , yt}
kþN
t¼kþ1, that is, the oldest data
sample is removed whenever a new sample becomes
available. It is apparent from the expressions of the estimates
of a and d given in the cited equations, that in a recursive
framework the most critical issue is the propagation of the
inverse of matrix T. Note that, in view of the sliding window
implementation, matrix T changes in such a way that its
dimensions remain fixed, so both a downdate and an update
have to be performed. More precisely, considering times k
and kþ 1, matrix T in (16) is formally given by
T (k) ¼ hk w
`
k
w
k
Dk


, T (kþ 1) ¼ Dk ukþ1u`kþ1 ykþ1


where by hk we denote the element of T (k) to be removed in
the downdate step, by vkþ1 the element of T (kþ 1) to be
added in the update step and where the underlined elements
indicate column vectors. Making use of the matrix block
inversion identity we obtain
T 1(k) ¼ g1k
1 w`
k
D1k
D1k wk gkI þ D
1
k wkw
`
k
 
D1k
2
4
3
5 (20)
with gk ¼ hk  w`k D
1
k wk [ R. Assume now that T
1(k) is
known from the previous iteration of the algorithm. Then, if
we partition it as
T 1(k) ¼
ak q
`
k
q
k
Qk
" #
we have that D1k can be deduced with respect to (20) as
D1k ¼ Qk 
1
ak
q
k
q`
k
After that, T 1(kþ 1) may be computed as
T 1(kþ1)¼ h1kþ1
(hkþ1I þD1k ukþ1u`kþ1)D1k D1k ukþ1
u`kþ1D1k 1
" #
(21)
with hkþ1 ¼ ykþ1  u`kþ1D1k ukþ1 [ R.T Control Theory Appl., 2009, Vol. 3, Iss. 9, pp. 1209–1216
i: 10.1049/iet-cta.2008.0339The dual parameters as and ds follow immediately from
T 1(kþ 1) and then, the non-linear function f (.). Using
the obtained estimate of the non-linearity, the parameters
Hi are computed in a recursive least squares fashion from (19).
4.2 Online extraction of the system
matrices
In order to recursively update the estimates of the system
matrices on the basis of the computed Markov parameters
Hi, results from the field of recursive subspace identification
(see, e.g. [1, 3, 4] and the references therein) will be
exploited. In particular, the so-called propagator algorithm
for system identification will be adopted [3, 4], which is
briefly summarised in the following. Consider system (1)
and define the data vectors
y n,t ¼ y`t    y`tþn1
`
un,t ¼ u`t    u`tþn1
`
where n is chosen such that nx , n  d . From (1), it is
possible to derive the following algebraic relation between
the input and output vectors un,t and y n, t
yn,t ¼ G nxt þH n f (un,t)þ en,t (22)
where Gn is the extended observability matrix
Gn ¼ (C)` (CA)`    (CA n1)`
`
Hn¯ is the block Toeplitz matrix of the impulse responses
from f (u) to y
H n ¼
D 0    0
CB D    0
..
. . .
. . .
. ..
.
CA n2B    CB D
2
66664
3
77775
¼
H0 0    0
H1 H0    0
..
. . .
. . .
. ..
.
Hn1    H1 H0
2
66664
3
77775
(23)
and en,t is the sum of all the noise terms acting on the system
output. Equation (22) can be equivalently written as
zn,t ¼ Gnxt þ e n,t (24)
where zn,t ¼ yn,t Hnf (un,t) is the so-called observation
vector. At each time step, zn, t can be explicitly computed as
zn,t ¼ yn,t  H^ n f^ (u n,t) thanks to the estimates of f (.) and of
the elements of Hn obtained from the algorithm presented in
the previous Section. Under the assumption that the system
in (1) is observable, Gn has at least nx linearly independent
rows which can be gathered in a non-singular submatrix G1
using row permutations (see [4] for details on how to perform1213
& The Institution of Engineering and Technology 2009

12
&
www.ietdl.orgthis step), so Gn can be partitioned as
Gn ¼
G1
PTn G1


¼
Inx
PTn
" #
G1 (25)
where the matrix PTn is the so-called propagator [3, 4]. Thus,
since G1 is non-singular, a basis of the observability subspace
for system (1) can be obtained by estimating the propagator,
using (24). Indeed (possibly applying a data reorganisation to
the observation vector so that the first nx rows of Gn are
linearly independent), (24) can be partitioned as
zn,t ¼
z1,t
z2,t


¼
Inx
PTn
" #
G1xt þ
e1,t
e2,t


(26)
where z1,t [ Rnx and z2,t [ R(ny nnx) are the components of
zn,t respectively corresponding to G1 and P
T
n G1, and e1,t and
e2,t represent a consistent partition of the noise vector en,t . In
the noise free case, it is easy to show that
z2,t ¼ PTn z1,t (27)
In the presence of noise, this relation no longer holds.However,
by assuming we have collected (or are about to collect) N I/O
measurements, an unbiased estimate of the propagator PTn
can be obtained by introducing an instrumental variable
h [ Rnh , assumed to be uncorrelated with the noise but
sufficiently correlated with the state vector x (see, e.g. [3] for
further details on the choice of IVs in propagator-based
recursive subspace identification), in the sense that
Renh ¼ 0 and rank Rxh
n o
¼ nx (28)
and by defining the cost function
J (P n) ¼ R^z2h,N  P
T
n R^z1h,N


2
F
(29)
The estimate of the propagator based ondata up to time t is then
obtained as P^
T
n,t ¼ R^z2h,t R^
1
z1h,t
(see [3, 4] for details on the
recursive estimation of Pn ) and the system matrices are finally
recovered from
C^ ¼ G^n(1:l )
A^ ¼ G^n(1:l ( n 1), :)
y G^n(l þ 1:l n, :)
D^ ¼ H^ 0, B^ ¼ G^n(1:l ( n 1), :)
y
H^ 1
..
.
H^ n1
2
664
3
775
where † denotes the Moore-Penrose pseudo-inverse and
G^n ¼
Inx
P^ n
" #
.
4.3 Summary of the recursive algorithm
Finally, the complete recursive identification algorithm can
be summarised as follows. Whenever a new pair (ut , yt) of
data is measured:14
The Institution of Engineering and Technology 20091. Compute the non-linear function by computing the
parameters a and d. That is realised in a recursive way by
using (21) and (16).
2. Update the estimate of the Markov parameters of the
linear part from (19) using recursive least squares.
3. Update the estimate G^ n of extended observability matrix
by recursively adapting the propagator P^ n.
4. Update the estimate A^, B^, C^, D^ of the system matrices.
5 Simulation examples
5.1 LTI system
In order to evaluate the convergence properties of the
proposed algorithm, the SISO second-order system
xtþ1 ¼
0:15 0:325
1:0 0


xt þ
2
0


f (ut)
yt ¼ 0:1750 1:1625

xt þ f (ut)þ yt
(30)
is considered, where the static non-linearity has been chosen
as f (u) ¼ sinc(u)u2.
The identification data, consisting of 3000 samples, have
been generated by using as excitation a zero mean white
gaussian noise sequence u(t)  N (0, 1). The variance of
the output measurement noise yt , also zero mean white and
gaussian, has been chosen such that the signal-to-noise
ratio is 30 dB. The algorithm parameters have been chosen
as s ¼ 1, d ¼ 10, g1 ¼ 0:003, n ¼ 6, and the instrumental
variable h taken as a vector of past inputs [ u(t  nx)   
u(t  1)]` of length nx. The simulation is driven with
N ¼ 300 over a time interval of 3000 points. Fig. 1
illustrates the products between the non-linearity and the
two first Markov parameters on the interval [25, 5]. The
actual functions and their respective estimates are almost
undistinguishable, at least on the density interval of the
input sequence. Fig. 2 shows the eigenvalues of the system
estimated recursively.
5.2 LTV system
Similarly, the ability of the proposed recursive algorithm of
dealing with a situation in which the dynamics of the
system generating the data varies with time has been
verified by considering as example a modified version of
system (30) in which the A matrix has been made slowly
time-varying according to
A(t) ¼ 1 102
ffiffi
t
p  0:15 0:325
1:0 0


(31)IET Control Theory Appl., 2009, Vol. 3, Iss. 9, pp. 1209–1216
doi: 10.1049/iet-cta.2008.0339

IE
do
www.ietdl.orgFigure 2 Evolution of the eigenvalues of the estimated
model (dotted) and of the actual system (30) (solid)
Figure 1 Estimates of the Markov parameters of system (30)
a True (solid line) and estimated (dotted line) Markov parameter
Df(.)
b True (solid line) and estimated (dotted line) Markov parameter
CBf(.)T Control Theory Appl., 2009, Vol. 3, Iss. 9, pp. 1209–1216
i: 10.1049/iet-cta.2008.0339The results obtained in tracking the eigenvalues of the
system, under the same excitation conditions considered in
the previous case, are depicted in Fig. 3. As can be seen
from the figure, the estimated eigenvalues show a quite
satisfactory tracking of the actual system parameters.
6 Concluding remarks
A recursive algorithm combining LS-SVM regression and
propagator-based subspace identification techniques has
been presented for the identification of SIMO
Hammerstein systems. The algorithm operates in a
recursive fashion and is designed for slowly time-varying
model identification. Numerical simulations show quite
good results. Future work will focus on the recursive
identification of more general non-linear state-space models.
7 Acknowledgment
This work was supported in part by the Italian National
Research Project ‘New Techniques of Identification and
Adaptive Control for Industrial Systems’.
8 References
[1] LOVERA M., GUSTAFSSON T., VERHAEGEN M.: ‘Recursive
subspace identification of linear and non-linear Wiener
state-space models’, Automatica, 2000, 36, (11),
pp. 1639–1650
[2] OKU H., KIMURA H.: ‘Recursive 4SID algorithms using
gradient type subspace tracking’, Automatica, 2002, 38,
(6), pp. 1035–1043
[3] MERCE`RE G., LOVERA M.: ‘Convergence analysis of
instrumental variable recursive subspace identification
algorithms’, Automatica, 2007, 43, (8), pp. 1377–1386
Figure 3 Recursive estimates of the eigenvalues for the
slowly time-varying system (31)1215
& The Institution of Engineering and Technology 2009

[4] MERCE`RE G., BAKO L., LECOEUCHE S.: ‘Propagator-based
methods for recursive subspace model identification’,
Signal Process., 2008, 88, (3), pp. 468–491
[5] BAI E.-W.: ‘An optimal two-stage identification algorithm
for Hammerstein-Wiener non-linear systems’, Automatica,
1998, 34, (3), pp. 333–338
[6] GOETHALS I., PELCKMANS K., SUYKENS J., DE MOOR B.:
‘Identification of MIMO Hammerstein models using least
squares support vector machines’, Automatica, 2005, 41,
(7), pp. 1263–1272
[7] DEMPSEY E., WESTWICK D.: ‘Identification of Hammerstein
models with cubic spline non-linearities’, IEEE Trans.
Biomed. Eng., 2004, 51, (2), pp. 237–245
[8] RAMOS J., DURAND J.-F.: ‘Identification of non-linear
systems using a b-splines parametric subspace approach’.
American Control Conf., San Diego, California, 1999
[9] VO¨RO¨S J.: ‘Modeling and identification of Wiener
systems with two-segment non-linearities’, IEEE Trans.
Control Syst. Technol., 2003, 11, (2), pp. 253–257
[10] VERHAEGEN M., WESTWICK D.: ‘Identifying MIMO
Hammerstein systems in the context of subspace model
identification methods’, Int. J. Control, 1996, 63, (2),
[11] GOETHALS I., PELCKMANS K., SUYKENS J., DE MOOR B.: ‘Subspace
identification of Hammerstein systems using least squares
support vector machines’, IEEE Trans. Autom. Control,
2005, 50, (10), pp. 1509–1519
[12] VERHAEGEN M., VERDULT V.: ‘Filtering and System
Identification – a Least Squares Approach’ (Cambridge
University Press, 2007)
[13] JUANG J.N., PHANM., HORTA L.G., LONGMANR.W.: ‘Identification of
observer/Kalman filter Markov parameters: theory and
experiments’, J. Guid. Control Dyn., 1993, 16, (2), pp. 320–329
[14] SUYKENS J., VAN GESTEL T., DE BRABANTER J., DE MOOR B.,
VANDEWALLE J.: ‘Least Squares Support Vector Machines’
(World Scientific, 2002)
[15] PELCKMANS K., GOETHALS I., DE BRABANTER J., SUYKENS J.A.K.,
DE MOOR B.: ‘Support vector machines: theory and
applications’, chapter Componentwise Least Squares
Support Vector Machine (Springer, 2005), pp. 77–98
[16] DE KRUIF B., DE VRIES T.: ‘Pruning error minimization in
least squares support vector machines’, IEEE Trans. Neural
Netw., 2003, 14, (3), pp. 696–702
[17] LIU J., CHEN J., JIANG S., CHENG J.: ‘Online LS-SVM for
function estimation and classification’, J. Univ. Sci.
12
&
www.ietdl.orgpp. 331–34916
The Institution of Engineering and Technology 2009Technol. Beijing, 2003, 10, (5), pp. 73–77IET Control Theory Appl., 2009, Vol. 3, Iss. 9, pp. 1209–1216
doi: 10.1049/iet-cta.2008.0339