Alternative Form of Predictor Based Identification of
LPV-SS Models with Innovation Noise

arXiv:1609.07536v1 [cs.SY] 23 Sep 2016

Pepijn Cox† and Roland Tóth†
Abstract— In this paper, we present an approach to identify
linear parameter-varying (LPV) systems with a state-space (SS)
model structure in an innovation form where the coefficient
functions have static and affine dependency on the scheduling
signal. With this scheme, the curse of dimensionality problem
is reduced, compared to existing predictor based LPV subspace identification schemes. The investigated LPV-SS model is
reformulated into an equivalent impulse response form, which
turns out to be a moving average with exogenous inputs (MAX)
system. The Markov coefficient functions of the LPV-MAX
representation are multi-linear in the scheduling signal and
its time-shifts, contrary to the predictor based schemes where
the corresponding LPV auto-regressive with exogenous inputs
system is multi-quadratic in the scheduling signal and its timeshifts. In this paper, we will prove that under certain conditions
on the input and scheduling signals, the ℓ2 loss function of
the one-step-ahead prediction error for the LPV-MAX model
has only one unique minimum, corresponding to the original
underlying system. Hence, identifying the LPV-MAX model in
the prediction error minimization framework will be consistent
and unbiased. The LPV-SS model is realized by applying an
efficient basis reduced Ho-Kalman realization on the identified
LPV-MAX model. The performance of the proposed scheme is
assessed on a Monte Carlo simulation study.

I. I NTRODUCTION
In recent years, identification of linear parameter-varying
state-space (LPV-SS) models has received considerable attention, e.g., [1]–[4], resulting in a succesful extention of
subspace identification (SID) and prediction error minimization (PEM) methods to the LPV case. The LPV-SS
identification framework has been applied to data-driven
modelling of wind turbines [5], forced Lorenz attractors [6],
power management of Web service systems [7], or to capture
traffic flow [8], to mention a few. A popular choice is to
capture the underlying system as a discrete-time LPV-SS
representation with an innovation noise model:
qx = A(p)x + B(p)u + K(p)e,
(1a)
y = C(p) x + D(p)u + e,
(1b)
where x : Z → X = Rnx is the state variable, y : Z →
Y = Rny is the measured output signal, u : Z → U =
Rnu denotes the input signal, p : Z → P ⊆ Rnp is the
scheduling variable, q is the forward time-shift operator, e.g.,
qx(t) = x(t + 1) where t ∈ Z is the discrete time, e : Z →
Rny is a sampled path of a zero-mean i.i.d. stationary noise
processes with Gaussian distribution, i.e., e(t) ∼ N (0, Σe )
with a nonsingular covariance Σe ∈ Rny ×ny . The matrix
functions A(·), ..., K(·) defining the SS representation (1)
† Pepijn Cox and Roland Tóth are with the Control Systems
Group, Department of Electrical Engineering, Eindhoven University of
Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands,

{p.b.cox,r.toth}@tue.nl

are usually taken to be affine combinations of p, to coincide
with the majority of LPV control synthesis methods, e.g., [9].
Hence, these matrix functions are defined as
np
X
Ai pi ,
A(p) = A0 +
(2)
i=1

where B(p), . . . , K(p) are equivalently parametrized as (2)
np
being constant matrices of apwith {Ai , Bi , Ci , Di , Ki }i=0
propriate dimensions.
Identification of LPV-SS models can be formulated in
the PEM framework, where the PEM methods can roughly
be categorized as: i) methodologies which assume that
full state measurements are available, e.g., [10]; ii) grey
box schemes, where only a small subset of all parameters
is estimated, e.g., [11]; iii) set-membership approaches,
e.g., [12]; or iv) direct PEM methods using gradient based
methodologies, e.g., [2], [3]. However, assuming that a full
measurement of the state or detailed model of the system is
available, is unrealistic in many applications. Therefore, only
set-membership and direct PEM methods classify as ‘blackbox’ identification schemes. In general, the set-membership
approaches have a significant higher computational load
compared to the direct PEM methods, as the set-membership
schemes estimate a feasible model set for which only heuristic techniques exist. In addition, these schemes have often
unrealistic noise assumptions and consistency of the overall
scheme is difficult to show. For the direct PEM methods,
the nonlinear and nonunique optimization problem is usually
solved by applying gradient-based search strategies or an
expectation maximization strategy. Regardless of the strategy,
direct PEM methods are solved in an iterative way, are prone
to local minima, and their convergence depends heavily on
a proper initial guess.
On the contrary, SID methods use convex optimization,
to identify a specific LPV input-output (IO) structure, from
which an LPV-SS model is realized by using matrix decomposition methods. Subspace schemes have extensively been
applied in the linear time invariant (LTI) case and extensions
to the LPV framework can be found, e.g., [1], [4], [13], [14].
Unfortunately, SID methods usually depend on approximations to get a convex problem, and, as a consequence, in the
LPV case, these approaches suffer heavily from the curse
of dimensionality and/or result in ill-conditioned estimation
problems with high parameter variances.
Based on these considerations, efficient estimation of LPVSS models remains a central problem to be solved.
In this paper, we would like to reduce the dimensionality
problems associated with the predictor based SID methods

by restating the LPV-IO identification setting to avoid multiquadratic dependency on the scheduling signal and its timeshifts. First, the LPV-SS representation (1) is formulated
in terms of its impulse response representation (IIR). As
Section II highlights, the corresponding IIR turns out to be
a moving average with exogenous inputs (MAX) representation, where the Markov coefficient functions are multi-linear
in the scheduling signal and its time-shifts. In Section III, the
identification of the multivariable MAX model is given. We
will proof that minimizing the ℓ2 loss function of the onestep-ahead prediction error has only one unique minimum
under some mild conditions on the input and scheduling
signals, which is the main contribution of this paper. To this
end, we extent the uniqueness proof of the LTI multivariable
moving average model [15] to the LPV case, which is
essential in showing uniqueness of the overall LPV-MAX
PEM identification method. Hence, applying a pseudo linear
regression to identify the LPV-MAX model results in a
consistent and unbiased estimate. After identifying the LPVMAX model, the LPV-SS model is estimated by utilizing
the bases reduced Ho-Kalman realization scheme of [16],
which, in this paper, is modified to also realize the noise
model. In Section V, the performance of the identification
method is assessed by a Monte-Carlo study, followed by
some conclusions in Section VI.
II. I MPULSE R ESPONSE R EPRESENTATION
To be able to solve our identification problem, we use an
alternative representation of (1):
Lemma 1 (IIR representation [17]): Any asymptotically
stable1 LPV system given in terms of representation (1) has
a convergent series expansion in terms of the pulse-basis
{q −i }∞
i=0 given by
∞
∞
X
X
(3)
(hj ⋄ p)q −j e,
(gi ⋄ p)q −i u +
y=
j=0

i=0

|

{z

process model

}

|

{z

noise model

}

where gi ∈ R ny ×nu and hj ∈ R ny ×ny are the expansion
coefficient functions, i.e., Markov coefficients, of the process
and noise dynamics, respectively. The ring of all real polynomial functions with finite dimensional domain is defined
by R, the operator ⋄ : (R, PZ ) → (Rny ×nu )Z denotes
(gi ⋄ p) = gi (p(t), . . . , p(t − τ )) with τ ∈ Z.

The IIR coefficients {gi , hi }∞
of
(3)
are
given
by
i=0
y = D(p) u + C(p)B(q −1 p) q −1 u+
| {z }
|
{z
}
g0 ⋄p

g1 ⋄p
−2

C(p)A(q −1 p)B(q
{z
|

C(p)K(q
|
{z

−1

h1 ⋄p

p) q
}

g2 ⋄p
−1

I e+
p) q −2 u + . . . + |{z}
}
h0 ⋄p

−1

e+C(p)A(q p)K(q
|
{z
h2 ⋄p

−2

p) q −2 e+. . . (4)
}

where gi , hi converges to the zero function in de ℓ∞ sense
as i → ∞. The Markov coefficients of the process part can
be written as
1 An

LPV system, represented in terms of (1), is called asymptotically stable in the deterministic sense, if, for all trajectories of
(u(t), e(t), p(t), y(t)) satisfying (1), with u(t) = 0, e(t) = 0 for t ≥ 0
and p(t) ∈ P, it holds that limt→∞ |y(t)|= 0.

gm ⋄ p = C(p)A(q −1 p) · · · A(q −(m−1) p)B(q −m p) =
np np
np np
X
X X
X
Ci Aj · · · Ak Bl pi (q −1 pj)· · ·(q −m pl), (5)
···
i=0 j=0

k=0 l=0

and the Markov coefficients of the noise part hm are similar
to (5), however, B(q −m p) is exchanged for K(q −m p). The
individual products Ci Aj · · · Ak Bl or Ci Aj · · · Ak Kl are
called the sub-Markov parameters of the process model and
noise model, respectively. Due to the convergence of gi and
hj , it is often sufficient to truncate (3):
nb
nc
X
X
(hj ⋄ p)q −j e,
(6)
(gi ⋄ p)q −i u +
y≈
j=0

i=0

where nb > 0 and nc > 0 are the order of the resulting finite
impulse response (FIR) models for the process and noise part.
The sub-Markov parameters in (5) have a multi-linear
dependency on the elements of p and its time-shifts, contrary
to the LPV auto-regressive with exogenous inputs (ARX)
formulation of other predictor based subspace schemes,
e.g. [1], which have multi-quadratic dependency. The LPVARX model2 uses Ã(p) = A(p) − K(p)C(p), therefore,
increasing the complexity of the LPV-IO
Pmodel to be identified.
To
see
this,
substitute
A(p)
=
Ai pi by Ã(p) =
PP
Ai pi − Ki Cj pi pj (with p0 = 1) in (5) and, similarly,
B by B̃. To compare the additional parameters, truncate the
LPV-ARX model (7) by orders na and nd for the u and y
polynomial,
then the
model (7) has
P dLPV-ARX 2j
P a respectively;
(1 + np ) ) parameters
(1 + np )2i+1 + ny nj=1
ny (nu ni=0
Pnb +1
(1 + np )i +
and the LPV model (6) has ny (nu i=1
Pnc +1
j
ny j=2 (1 + np ) ). For example, if nu=ny=np=na=nb=
nc =nd =2 then the LPV-ARX model has 1452 parameters
and the proposed LPV model has 300, hence, a significant
reduction is achieved. Therefore, a common assumption by
other predictor based subspace schemes is to take either
K or C,D to be constant matrices, avoiding the additional
modelling complexity. However, by identifying (4), we can
keep K, C, and D to be parameter dependent.
III. I DENTIFICATION

OF A

MAX

MODEL

A. Problem Setting
The Markov coefficients (4) of the impulse response (6)
representing the LPV-SS system (1) are functions in the
scheduling signal and its time shifts, i.e., (5). Hence (6) is a
moving average with exogenous inputs (MAX) system:
y = B(p, q −1 )u + C(p, q −1 )e,
(8)
with the process and noise filter given by
nc
nb
X
X
Cj (p)q −j , (9)
Bi (p)q −i , C(p, q −1 ) =
B(p, q −1 ) =
j=0

i=0

2 The LPV-ARX representation is found by substituting e of (1b) into (1a)
and, by using this modified state equation, writing out (1b) as:
∞
i
h i−1
X
Y
Ã(q −j p) B̃(q −i p)q −i u
y = D(p)u +
C(p)
i=1

+

∞
X
i=1

j=1

h i−1
i
Y
C(p)
Ã(q −j p) K(q −i p)q −i y + e

(7)

j=1

where Ã(p) = A(p)-K(p)C(p), B̃(p) = B(p)-K(p)D(p), and

Q0

j=1 Ã = I.

where Bi (p) = (gi ⋄ p) and Cj (p) = (hj ⋄ p). Hence,
Bi (p) and Ci (p) are multi-linear matrix functions in the
scheduling signal and its times-shifts from t, . . . , t−i, similar
to the Markov coefficients (5). In this paper, the LPV-MAX
model (9) is identified by using the PEM setting. To avoid
identifiability issues, (1) with functional dependencies (2)
is assumed to be structurally state observable and state
reachable w.r.t. both u and e in the dertiministic sense3 ,
which implies joint minimality of (1).
In order to apply the PEM framework, we select the model
structure similar to (9):
y = B̂(θ, p, q −1 )u + Ĉ(θ, p, q −1 )ε(θ),
(10)
where ε(θ) : Z → Y is the one-step-ahead prediction error
and the process B̂ and noise Ĉ models are considered to
be polynomials of q −1 with p-dependent coefficients similar
to B and C in (9) with orders n̂b , n̂c , respectively. These
polynomials are parametrized in terms of the sub-Markov
parameters resulting in an overall parameter vector θ.
The data-generating system (1) is aimed to be identified using a given identification dataset DN
=
{u(t), p(t), y(t)}N
t=1 generated by (1) and minimizing the
following ℓ2 loss function
!
N
1 X
⊤
ε(θ, t)ε (θ, t)
(11)
VN (θ) = Tr
N t=1
where N is the number of data points collected in DN . Under
weak regularity conditions it is well known that
N


1X
Tr E ε(θ, t)ε⊤(θ, t) , (12)
VN (θ) → V∞ (θ) = lim
N→∞ N
t=1

with probability one and uniformly in θ [18], which removes
the effect of sample based realizations of u, p, and e in the
analysis. Due to the uniform convergence of (12), the results
of Section III-C and III-D give information about the shape
of VN (θ) for large N .
First, let us define a notation to indicate which sub-Markov
parameters (5) of the filters B, C are selected. Denote with
Ivs the set {s, s + 1, . . . , v}. Then, [Ivs ]n defines the set of
all sequences of the form (i1 , . . . , in ) with i1 , . . . , in ∈ Ivs .
The elements of Ivs will be viewed as characters and the
finite sequences of elements of Ivs will be referred to as
strings. Then [Ivs ]n is the set of all strings containing exactly
n characters.
n a selection
 n Then
 n n with n ≥n 0 is constructed
 n n
from η ∈ I0 p 0 with I0 p 0 = {ǫ} ∪ I0 p ∪ . . . ∪ I0 p
 1 2
and ǫ denoting the empty string. As an example, I0 0 =
{ǫ, 0, 1, 00, 01, 10, 11}. Define by #(η) the amount of characters of a single string in the set. With this notation, we
will simplify the notation of the sub-Markov parameters as
C[ηb ]1 A[ηb ]2 · · · A[ηb ]i−1 B[ηb ]i = gηb ,
(13)

where i = #(ηb ), [η]j denotes the j-th character of the string
 n n +1
 n n +1
is the set of all
η and ηb ∈ I0 p 2 b . Hence, I0 p 2 b
possible combinations of the sub-Markov parameters for the
process filter. Similarly, we indicate a sub-Markov parameter
of the noise filter C, process model B̂, and noise model Ĉ
3 There exists at least one p ∈ PZ such that n -step observability and
x
reachability holds for all time moments on the support of the signals [17].

 n n̂c +1
 n n̂ +1
 n nc +1
,
, ηb̂ ∈ I0 p 2 b , and ηĉ ∈ I0 p 2
by ηc ∈ I0 p 2
4
respectively. Then define

(14a)
Q(θ) = Ē ε(θ, t)ε(θ, t)⊤
∞
X
Γ̂i (θ, p, t)q −i
Γ̂(θ, p, q −1 ) =
i=0

= Ĉ−1 (θ, p, q −1 ), (Γ0 = I)
⊤

Rk (θ, t) = ε(θ, t − k)ε (θ, t), k ≥ 0,
where

(14b)
(14c)

N
1 X
E {·} .
N →∞ N
t=1

Ē {·} = lim

Remark, we use a simplified notation Γ̂i (θ, p, t) = (Γ̂i (θ) ⋄
p)(t) where each Γ̂i (θ, p, t) depends on p(t), . . . , p(t − i).
B. Pseudo linear regression
The advantage of the given problem setting of Section III-A is that it can easly be solved by applying global
pseudo linear regression methods, e.g., [19, Algorithm 3].
Pseudo linear regression has a relatively low computational
load. Due to space limitations, the algorithm is not presented
here. If we take n̂b and n̂c to be finite, the corresponding
LPV-MAX model can be seen as an extended LTI model.
In the LTI case, pseudo linear regression has been studied
extensively and it is known that it will converge if no
auto-regressive part is present [20]. However, a key part
of the convergence is to have only one unique solution
of the minimization problem (12), which is analysed in
Sections III-C and III-D.
C. Global minima
To show that a global minimum exists, the following
assumption is taken:
A1 The input signal u is uncorrelated to the noise signal e.
With the aforementioned identification setting, we get:
Theorem 2: Given the data generating system (8) with
functional dependencies (9), model structure (10), n̂b ≥ nb ,
n̂c ≥ nc , and u satisfying A1; then VN (θ) has a global
minimum
min V∞ (θ) = Tr(Σe ),
(15)
θ

which can be obtained if the elements of θ satisfy
(16a)
ĝηb = gηb , ĥηc = hηc ,
 np nc+1
 np nb+1
and ηc ∈ I0 2 . In addition,
for all indices ηb ∈ I0 2
if n̂b > nb , n̂c > nc ,
ĝηb̂ = 0, ĥηĉ = 0,
(16b)
 np n̂b+1
 np n̂c+1
for all ηb̂ ∈ I0 n +2 and ηĉ ∈ I0 n +2 .

c
b
Proof:
For notational ease, denote B(p, q −1 ),
B̂(θ, p, q −1 ), C(p, q −1 ), Γ̂(θ, p, q −1 ), and ε(θ) as B, B̂, C,
Γ̂, and ε, respectively. Then, using (8) and (10), let us rewrite
ε in terms of u and e as
h

i
ε = Γ̂

B − B̂ u + Ce .

Hence, as e and u are uncorrelated (A1), the loss function (12) separates as

4 The inverse filter Γ̂(θ, p,q −1 ) of Ĉ(θ, p, q −1 ) in (14b) is generated by
P
−1 ))i [19].
Γ̂(θ, p,q −1 ) = ∞
i=0 (I − Ĉ(θ, p, q

V∞ (θ) = V∞,1 (θ) + V∞,2 (θ),

(17)

with


 oi
h n 
V∞,1 (θ) = Tr Ē Γ̂ B− B̂ u(t)u⊤(t) B⊤− B̂⊤ Γ̂⊤ , (18a)
h n
oi
V∞,2 (θ) = Tr Ē Γ̂Ce(t)e⊤(t)C⊤ Γ̂⊤ .
(18b)

Obviously, V∞,1 (θ) ≥ 0 with equivalence if B ≡ B̂. For
V∞,2 (θ), remark that the filters C and Ĉ are monic,
ε = Γ̂Ce = e + v,
(19)
where the random signal v(t) depends linearly on the past
samples of e, i.e., {e(τ )}t−1
τ =−∞ , but is independent of e(t).
Hence, it follows that 

V∞,2 (θ) = Tr Ē (e(t) + v(t))(e(t) + v(t))⊤



= Tr E{e(t)e⊤(t)} +Tr Ē v(t)v⊤(t)

≥ Tr E{e(t)e⊤(t)} = Tr(Σe ).
(20)

Eq. (20) holds with equality if v ≡ 0, which follows if Γ̂C =
I implying C ≡ Ĉ. Hence, V∞ (θ) = Tr(Σe ) if (16) holds,
which is a global minimum of minθ V∞ (θ).
Theorem 2 does not imply uniqueness of the solution, only
that a global minimum exists if the sub-Markov parameters
of the model (10) are equal to those of the original system (9). Uniqueness of the solution is proven next.
D. Stationary point
Next, to show uniqueness of the of minimum in Theorem 2, we prove that there exists only one stationary point
of (15), hence, the global minimum of (15) is unique and it
can always be obtained by optimization methods. To prove
this, the following additional assumptions are taken:
A2 Each signal pi is assumed to be a zero-mean white noise
process with finite variance and independent of e. The
processes pi are mutually independent and the higher
order moments are bounded, given by E{pki plj } =
np
ck+j
i,j < ∞ for i, j ∈ I1 and 2 ≥ k + j ≥ 6 + n̂c .
A3 The given input u is chosen such that the signal




1
1
ũ(t) =
⊗ ···⊗
⊗ u(t) (21)
p(t − n̂b )
p(t)
Pn̂b
is persistently exciting of order ny i=1 (1+np )i nu [18].
The first assumption, the choice of the scheduling signal, can
be made less restrictive, however, it will make the analysis
more involved. In this case, we restrict the scheduling signal
to be stationary, hence, Tr(Ē{·}) = Tr(E{·}) for V∞,2
in (18b). The second assumption is needed to be able to
uniquely distinguish all parameters in B̂, which is discussed
in the proof of Theorem 5.
As the loss function V∞ can be separated as given in (17),
we will first focus on the analysis of the loss function w.r.t.
the noise model, i.e., V∞,2 (18b). The stationary points
of (18b) are solutions of the following equations
 n n̂c+1
∂V∞,2 (θ)
= 0, 1 ≤ i, j ≤ ny , ∀ηĉ ∈ I0 p 2 , (22)
∂θij,ηĉ
where θij,ηĉ only contains the parameters of the noise model.
Using the notation of (14), the partial derivative (22) is


∂ε(θ, t)
∂TrQ(θ) ∂Q(θ)
⊤
. (23)
= 2Ē ε (θ, t)
0 = Tr
∂Q(θ) ∂θij,ηĉ
∂θij,ηĉ

Taking the partial derivative of (10) w.r.t. the model parameters gives
∂ε(θ, t)
+ si s⊤
(24)
0 = Ĉ(θ, p, q −1 )
j pηĉ (t)ε(θ, t − k)
∂θij,ηĉ
for each t ∈ Z, where k = #(ηĉ ) − 1, si is a selector
vector for which only the i-th element is non-zero, and pηĉ (t)
defines the product Q
of different scheduling signals and its
ĉ )−1
time-shifts pηĉ (t) = #(η
p[ηĉ ]i (t−i) and, for notational
i=0
simplicity, p0 (t) = 1, ∀t ∈ Z. See that (24) is equivalent to
∂ε(θ, t)
= −Γ̂(θ, p, q −1 )si s⊤
(25)
j pηĉ (t)ε(θ, t − k),
∂θij,ηĉ
and substituting (25) in (23) gives
)
(
∞
X
Γ̂r (θ, p, t)si s⊤
0 = Ē ε⊤(θ, t)
j pηĉ (t − r)ε(θ, t − k − r)
=

s⊤
j Ē

(

r=0

ε(θ, t − k)

∞
X

⊤

)

ε (θ, t + r)Γ̂r (θ, p, t)pηĉ (t − r) si ,

r=0

Hence, we find that (22) is equivalent to:
(∞
)
X
 n n̂c+1
Ē
Rk+r (θ, t)Γ̂r (θ, p, t)pηĉ (t − r) = 0, ∀ηĉ ∈ I0 p 2 ,
r=0

(26)
with k = #(ηĉ ) − 1. We would like to highlight that
the analysis given next differs from the LTI case [15] as
ε is not independent of p, e.g., Γ̂r and Rk+r are not
independent in (26). Hence, this dependency makes the
analysis more involved. Note that the time-shift operator q
is non-communicative if applied on a coefficient function,
e.g., q −1 Ci (p) = Ci (q −1 p)q −1 [21]. Hence, for the following
analysis, the time-instance will be indicated as Cj (p, t − i)
for Cj (p(t − i), . . . , p(t − i − j)). First, decompose Rk+r as:
Lemma 3: If A2 holds, then are the matrices Rr for r ≥ 1
defined in (14c) given by

Ē{Rr (θ, t)} = Ē e(t−r)e⊤(t−r)Ω⊤(θ, p, t−r)Γ̄⊤
r (θ, p, t)
∞
nX
Γ̄i (θ, p, t − r)Ω(θ, p, t − i − r)e(t − r − i)·
+ Ē
i=1

o
e (t − r − i)Ω⊤(θ, p, t − i − r)Γ̄⊤
i+r (θ, p, t)
⊤

(27)

with


Γ̄k (θ, p, t) = Γ̂k−1 (θ, p, t) . . . Γ̂k−n̂c (θ, p, t) , (28a)


C1 (p, t + 1) − Ĉ1 (θ, p, t + 1)


..


.




 Cnc (p, t + nc ) − Ĉnc (θ, p, t + nc ) 
Ω(θ, p, t) = 
, (28b)


−Ĉnc +1 (θ, p, t + nc + 1)


..




.
−Ĉn̂c (θ, p, t + n̂c )
where Γl (·) = 0 for l < 0.

Proof: See Appendix.
Now we return to the uniqueness of V∞,2 :
Lemma 4: Given the data generating system (8) with
functional dependencies (9), model structure (10) with n̂c ≥
nc , B(p, q −1 ) = B̂(θ, p, q −1 ) ≡ 0, then, under assumption A2 the global minimum for

min V∞,2 (θ) = Tr(Σe ),
θ

(29)

is the unique stationary point of V∞,2 .

Proof: See that (27) can be substituted in (26) as e is
independent from p. Hence,
∞ h
nX
⊤
e(t−r−k)e⊤(t−r−k)Ω⊤(θ, p, t−r−k)Γ̄r−k
(θ, p, t)
Ē

Theorem 5: Given the data generating system (8) with
functional dependencies (9), model structure (10), n̂b ≥ nb ,
n̂c ≥ nc , then under assumptions A1-A3 the global minimum
min V∞ (θ) = Tr(Σe ),
θ

(33)

is the only stationary point of V∞ and is found if (16) is
satisfied.

∞
X
Proof:
Recall
from
(17)
that
V
=
V
+
V
∞
∞,1
∞,2 .
+ Γ̄i (θ, p, t−r−k)Ω(θ, p, t−i−r−k)e(t−r−i)e⊤(t−i−r−k) Lemma 4 proofs that for V
there
is
only
one
stationary
∞,2
i=1
i
o
point. Hence, C(p, q −1 ) = Ĉ(θ, p, q −1 ) at this stationary
Ω⊤(θ, p, t−i−r−k)Γ̄⊤
point. The underlying filter Ĉ is a polynomial in q −1 with
i+r−k (θ, p, t) Γ̂r (θ, p, t)pηĉ (t−r) = 0,
(30) parameter-varying coefficients and it is monic; hence it is
 np n̂c+1
for all ηĉ ∈ I0 2
and k = #(ηĉ ) − 1. Remark that full rank in the functional sense and, therefore, its inverse Γ̂
ĥηĉ = 0 cannot be a solution of (30) as Γ̂0 = I. Also see is also full rank. Next, we need to prove that V∞,1 has only
that the maximal amount of products of p with equivalent one stationary point. If we consider ũ as the input signal,
time-shift is 6 + n̂c in (30). Therefore, in A2, we have that then the filters B and B̂ can be written as a multi-input LTI
input signal ũ is persistent
Pn̂b
all moments of p up to the 6 + n̂c moment are needed to be filter, similar to (32). The extended
(1
+
np )i nu , thus [B − B̂]u
of
excitation
of
order
n
y
i=1
bounded. To treat (30) in the rest of the proof would require
−1
additional technical details, hence, instead we present just the cannot be zero for B 6= B̂. Recall, Γ̂(p, q ) also cannot lose
concept of the proof w.r.t. the restricted case of ηĉ ∈ [0]n̂2 c+1 rank, hence, there exists only one stationary point of V∞,1
such that pηĉ (t − r) = 1, then, due to the construction of Γ̂i , satisfying (16), e.g., see [18]. Therefore, the only stationary
it is impossible that summations and products of Γ̂i cancel point of V∞ is found for (16).
Theorem 5 proofs that under the conditions A1-A3 there
each other, which could lead to the zero solution of (30).
can
only exists one stationary point of (12). Hence, identificaHence, the solution can only be found for
tion
of the LPV-MAX model (10) is consistent and unbiased.

Ē Ω(θ, p, t − j)e(t − j)e⊤(t − j)Ω⊤(θ, p, t − j) = 0, (31)
IV. LPV-SS REALIZATION
for j > 0. We can rewrite Ω as
r=0

Ω(θ, p, t − j) = (θh − θ̂ĥ )P (t − j),

(32)

with
h[n] = [ h0...00 h0...01 . . . h0...0nc h0...10 . . . hnc ...nc nc ] ,
 [1]

h
0
... 0
 ..


... 0 
θh =  .
,
 h[nc ]
0 ... 0 
0
... 0
n+1

where h[n] ∈ Rny ×ny (1+np )
denotes the matrix of all
sub-Markov parameters associated with the hn -th Markov
coefficient. For θh , the last nc − n̂c rows and columns are
zero, such that θh and θ̂ĥ have equivalent dimensions. The
matrix θ̂ĥ is similarly parametrized as θh , however, with ĥ,
n̂c in stead of
 nc , h, and



1
1
P (t − j) =
⊗ ··· ⊗
⊗ Iny ,
p(t − j + n̂c )
p(t − j)
where ⊗ is the Kronecker product. Combing gives

(θh − θ̂ĥ )E P (t − j)e(t − j)e⊤(t − j)P⊤(t − j) (θh − θ̂ĥ )⊤= 0,
{z
}
|
Λ

where Λ is clearly positive definite if A2 holds, therefore,
θh = θ̂ĥ is the only solution for pηĉ (t − r) = 1. If
pηĉ (t − r) 6= 1, then the proof is more involved, but relies on
the above repeated concept. Hence, it is not presented here.
Concluding, θh = θ̂ĥ is the only solution for which the set
of equations (26) is satisfied. Therefore, it is the only unique
stationary point of V∞,2 .
Lemma 4 is the LPV extension of the LTI result of [15].
Using this preliminary analysis the following result holds:

The next step in the proposed identification scheme is to
realize the LPV-SS model from the identified LPV-MAX
model. By treating
the noise
as an additional input, i.e.,


extending B̃i = Bi Ki for i = 0, . . . , np , an isomorphic
to the original LPV-SS representation (1)5 is obtained by
employing the bases reduced Ho-Kalman realization scheme
of [16]. This bases reduced realization can considerably
decrease the size of the Hankel matrix by selecting only
its non-repetitive elements and, therefore, reducing the computational load, compared to realization on the full Hankel
matrix [4, Eq. (48)]. In the basis reduced realization, the SVD
is only applied on a no ×
x instead
Pnir matrix with no , nr ≥ nP
j
of a matrix with size ny l=1 (1+np)l ×(nu +ny ) l=1 (1+
np )l , for i, j ≥ 2 in the full realization case. The realization
scheme allows to reconstruct the state bases of the LPV-SS
model based upon only the process B or noise model C,
simply by only selecting their corresponding sub-Markov to
fill the Hν,ς matrix in [16, Eq. (20)]. This can be useful, if,
for example, one of both models is poorly estimated then the
other model can be used to more accurately reconstruct the
reachability and observability matrix.
V. S IMULATION E XAMPLE
The proposed identification scheme is tested on the benchmark model used in [22]. The benchmark contains an MIMO
LPV-SS model with input dimension nu = 2, scheduling
dimension np = 2, state dimension nx = 2, and output
5 There exists a constant, nonsingular transformation matrix T ∈ Rnx ×nx
n
such that: T Âi = Ai T , T [B̂i K̂i ] = [Bi Ki ], and Ĉi = Ci T , ∀i ∈ I0 p .

dimension ny = 2. We added the following parameter
independent function 

0.32 0.16
K=
.
0.64 0.24
The simulation output or one-step-ahead predicted output ŷ
of the estimated model is compared to measured output or the
one-step-ahead predicted output y of the oracle, respectively,
by means of the best fit rate (BFR)6
)
(
1 PN
t=1kyt − ŷt k2
N
, 0 · 100%,
(34)
BFR = max 1− 1 PN
t=1kyt − ȳk2
N
using a validation dataset Dval as in [22]. In (34), ȳ
defines the mean of the simulation output or one-stepahead predicted output y of the oracle. In the realization
step, the basis reduced scheme uses nr = 10, no = 8
bases, where the controllability matrix is spanned by ς =
{(ǫ, 0, 1),(ǫ, 0, 2),(ǫ, 1, 1),(0, 0, 1),(0, 0, 2),(1, 0, 1),(1, 0, 2),
(2, 0, 1),(2, 0, 2),(1, 1, 1)} and the observability is spanned
by ν = {(1, 0, ǫ),(2, 0, ǫ),(1, 0, 0),(2, 0, 0),(1, 0, 1),(2, 0, 1),
(1, 0, 2),(2, 0, 2)}. The truncation order is n̂b = 4 and n̂c =
2. To evaluate the statistical properties of the identification
scheme, NMC = 100 Monte Carlo runs are carried out.
In each run, a new realization of the input and scheduling
signal is used. We use two different methods of finding
the search direction in each of the pseudo linear regression
(PLR) iterations: i) ℓ2 regularized least squares estimate,
or ii) enhanced Gaus-Newton optimization method [23]. The
extension is trivial and, therefore, is not given here. For
the ℓ2 regularized least squares estimate, we optimized the
regularization parameters by a line-search and found λ1 = 1,
λ2 = 100 for SNR = ∞dB and λ1 = 0.1, λ2 = 1 for
SNR = {40, 10}dB (λ1 is w.r.t. the process parameters and
λ2 w.r.t. the noise parameters).
Table I displays the mean and standard deviation of the
BFR of the identification algorithm for different SNRy =
{∞, 40, 10}dB. The table indicates that all approaches are
capable of identifying the underlying dynamics. However,
our approach is mildly outperformed by [1]. The approache
of [1] has numerically efficient implementation, in terms of
a kernel based approach, which can also be used in our
case. In addition, we noticed that the iterations of the PLR
are not very robust, as reflected by the increased standard
deviation. On the other hand, with this paper, we would
like to show how to conceptually reduce the amount of
parameters used to identify the input-output model if K,
C, and D are parameter varying. Hence, evolving to an
numerical efficient implementation is for future research.
VI. C ONCLUSIONS
In this paper, an identification scheme to estimate LPVSS models in innovation form with static and affine dependence on the scheduling signal has been presented. The
LPV-SS representation has been reformulated into its corresponding LPV-MAX system, which significantly reduced
the amount of parameters of the LPV-IO model compared
6 Usually

the BFR are defined per channel. Eq. (34) are the average
performance criteria over all channels.

to other predictor based schemes. We had proven that the
LPV-MAX model could consistently be identified using the
prediction error minimization framework, under some mild
conditions on the input and scheduling signals. The LPVSS model has been realized from the LPV-MAX model
by using the basis reduced Ho-Kalman realization scheme.
This realization scheme significantly reduces the amount of
parameters needed to realize the LPV-SS model. Hence, the
overall scheme decrease the computational load significantly,
especially with moderate to large LPV-SS models where K,
C, and D are parameter varying. The proposed identification
scheme has been illustrated with a simulation example.
A PPENDIX
P ROOF OF L EMMA 3
First, we show how the sequence {Γ̂i (·)} is generated.
Note that Γ̂(θ, p, t)Ĉ(θ, p, t) = I. Hence, collecting all terms
in the polynomial with equivalent time-shift q −i , it follows
that for i ≥ 1
Γ̂i (θ, p, t) + Γ̂i−1 (θ, p, t)Ĉ1 (θ, p, t − i + 1)
+ . . . + Γ̂i−n̂c (θ, p, t)Ĉn̂c (θ, p, t − i + n̂c ) = 0, (35)
with Γ̂i (·) = 0 for i < 0 and Γ̂0 (·) = I. For example,
Γ̂1 (θ, p, t) = −Ĉ1 (θ, p, t),
(36a)
Γ̂2 (θ, p, t) = Ĉ1 (θ, p, t)Ĉ1 (θ, p, t − 1) − Ĉ2 (θ, p, t). (36b)
Second, rewrite Rr (θ, t) = ε(θ, t − r)ε⊤(θ, t) as
∞
nX
Γ̂i (θ, p, t−r)C(p, t−r−i)e(t−r−i)·
E {Rr (θ, t)} = E
i=0

∞
hX
j=0

i⊤ o
Γ̂j (θ, p, t)C(p, t − j)e(t − j)
. (37)

As p is uncorrelated with e, it is clear that all terms in (37)
for r − k 6= j are zero. Hence,
∞
nX
E {Rr (θ, t)} = E
Γ̃i (θ, p, t−r)C̃(p, t−r−i)e(t−r−i)·
i=0

o
e (t − r − i)C̃⊤(p, t − i − r)Γ̃⊤
(θ,
p,
t)
, (38)
i+r
⊤

with

i
h
Γ̃i (θ, p, t) = Γ̂i (θ, p, t) . . . Γ̂i−nc (θ, p, t) ,
⊤

⊤
C̃(p, t) = C⊤
.
0 (p, t) . . . Cnc (p, t + nc )

However, using the construction of Γ̂ (35) and the fact that
Γ̂ and C are monic polynomials, the product between Γ̂ and
C can be simplified. For example,
Γ̃1 (θ, p, t)C̃⊤ (p, t−2) = Ĉ1 (θ, p, t)Ĉ1 (θ, p, t−1)−Ĉ2 (θ, p, t)
− Ĉ1 (θ, p, t)C1 (p, t − 1) + C2 (p, t) = C2 (p, t) − Ĉ2 (θ, p, t)
h
i
+Γ1(θ, p, t) C1 (p, t − 1)− Ĉ1 (θ, p, t − 1) = Γ̄2 (θ, p, t)Ω(θ, p, t−2)

Hence, by using (28), the above example generalizes to

Γ̃i+r (θ, p, t)C̃(p, t − i − r) = Γ̄i+r (θ, p, t)Ω(θ, p, t − i − r),
(39a)
and, for i > 0
Γ̃i (θ, p, t − r)C̃(p, t − i − r) = Γ̄i (θ, p, t − r)Ω(θ, p, t − i − r).
(39b)
Substituting (39) in (38) results in (27).

TABLE I: Mean and standard deviation (std) of the BFR of the identification algorithm per Monte-Carlo run for different SNRy = {∞, 40, 10}dB. The
BFR is based on the simulated output and one-step-ahead predicted output of the estimated model on the Dval for NMC = 100 Monte-Carlo simulations.
GN indicates the enhanced Gaus-Newton optimization method [23]. The PBSID has a past and future window equal to 3.
∞dB
40dB
10dB

LPV-MAX and
Simulation
96.46 (0.5843)
96.52 (0.5954)
90.78 (1.570)

Regularized LS
Prediction
96.46 (0.5840)
97.08 (0.4496)
84.98 (1.780)

LPV-MAX and GN
Simulation
Prediction
96.69 (0.5511) 97.21 (0.4354)
96.87 (0.5429) 97.17 (0.4163)
90.76 (1.662)
84.29 (1.804)

R EFERENCES
[1] J. W. van Wingerden and M. Verhaegen, “Subspace identification of
bilinear and LPV systems for open- and closed-loop data,” Automatica,
vol. 45, no. 2, pp. 372–381, 2009.
[2] A. Wills and B. Ninness, “System identification of linear parameter varying state-space models,” in Linear parameter-varying system
identification: new developments and trends, P. Lopes dos Santos, T. P.
Azevedo Perdicoúlis, C. Novara, J. A. Ramos, and D. E. Rivera, Eds.
World Scientific, 2011, ch. 11, pp. 295–316.
[3] V. Verdult, N. Bergboer, and M. Verhaegen, “Identification of fully
parameterized linear and nonlinear state-space systems by projected
gradient search,” in Proc. of the 13th IFAC Symposium on System
Identification, Rotterdam, The Netherlands, Aug. 2003, pp. 737–742.
[4] R. Tóth, H. S. Abbas, and H. Werner, “On the state-space realization
of LPV input-output models: practical approaches,” IEEE Trans. on
Control Systems Technology, vol. 20, no. 1, pp. 139–153, Jan. 2012.
[5] J. W. van Wingerden, I. Houtzager, F. Felici, and M. Verhaegen,
“Closed-loop identification of the time-varying dynamics of variablespeed wind turbines,” Int. J. of Robust and Nonlinear Control, vol. 19,
no. 1, pp. 4–21, 2009.
[6] W. E. Larimore, P. B. Cox, and R. Tóth, “CVA identification of
nonlinear systems with LPV state-space models of affine dependence,”
in Proc. of the American Control Conf., Chicago, IL, USA, Jul. 2015,
pp. 831–837.
[7] M. Tanelli, D. Ardagna, and M. Lovera, “On-and off-line model
identification for power management of Web service systems,” in Proc.
of the 47th IEEE Conf. on Decision and Control, Cancun, Mexico,
2008, pp. 4497–4502.
[8] T. Luspay, B. Kulcsár, J. W. Van Wingerden, and M. Verhaegen, “On
the identification of LPV traffic flow model,” in Proc. of the European
Control Conf., Budapest, Hungary, Aug. 2009, pp. 1752–1757.
[9] J. Mohammadpour and C. Scherer, Eds., Control of linear parameter
varying systems with applications. Springer, 2012.
[10] M. Nemani, R. Ravikanth, and B. A. Bamieh, “Identification of linear
parametrically varying systems,” in Proc. of the 34th IEEE Conf. on
Decision and Control, New Orleans, LA, USA, 1995, pp. 2990–2995.
[11] P. Gáspár, Z. Szabó, and J. Bokor, “A grey-box identification of an
LPV vehicle model for observer-based side-slip angle estimation,” in
Proc. of the American Control Conf., New York City, USA, Jul. 2007,
pp. 2961–2965.
[12] C. Novara, “SM Identification of State-Space LPV systems,” in Linear
parameter-varying system identification: new developments and trends,
P. Lopes dos Santos, T. P. Azevedo Perdicoúlis, C. Novara, J. A.
Ramos, and D. E. Rivera, Eds. World Scientific, 2011, ch. 11, pp.
65–93.
[13] W. E. Larimore, “Identification of nonlinear parameter-varying systems via canonical variate analysis,” in Proc. of the American Control
Conf., Washington, DC, USA, Jun 2013, pp. 2247–2262.
[14] P. Lopes dos Santos, J. A. Ramos, and J. L. M. de Carvalho,
“Identification of linear parameter varying systems using an iterative
deterministic-stochastic subspace approach,” in Proc. of the European
Control Conf., Kos, Greece, Jul 2007, pp. 4867–4873.
[15] P. Stoica and T. Söderström, “Uniqueness of prediction error estimates
of multivariable moving average models,” Automatica, vol. 18, no. 5,
pp. 617–620, 1982.
[16] P. B. Cox, R. Tóth, and M. Petreczky, “Estimation of LPV-SS models
with static dependency using correlation analysis,” in Proc. of the
1st IFAC Workshop on Linear Parameter Varying Systems, Grenoble,
France, Oct. 2015, pp. 91–96.
[17] R. Tóth, Modeling and identification of linear parameter-varying
systems. Springer, 2010.
[18] L. Ljung, System identification: theory for the user. Springer, 1989.
[19] R. Tóth, P. S. C. Heuberger, and P. M. J. Van den Hof, “Predictionerror identification of LPV systems: present and beyond,” in Control

[20]

[21]
[22]
[23]

Simulation
99.50 (0.1164)
99.48 (0.1306)
93.06 (0.8757)

[1]
Prediction
99.64 (0.0834)
99.04 (0.0553)
85.25 (0.5100)

of Linear Parameter Varying Systems with Applications, J. Mohammadpour and C. W. Scherer, Eds. Springer, 2012, ch. 2, pp. 27–58.
L. Ljung, T. Söderström, and I. Gustavsson, “Counterexamples to
general convergence of a commonly used recursive identification
method,” IEEE Trans. on Automatic Control, vol. 20, no. 5, pp. 643–
652, 1975.
R. Tóth, “Identification problems in the context of compressive sensing,” Tech. Rep., 2010.
V. Verdult and M. Verhaegen, “Kernel methods for subspace identification of multivariable LPV and bilinear systems,” Automatica, vol. 41,
pp. 1557–1565, 2005.
A. Wills and B. Ninness, “On gradient-based search for multivariable
system estimates,” IEEE Trans. on Automatic Control, vol. 53, no. 1,
pp. 298–306, Feb. 2008.

