close

Вход

Забыли?

вход по аккаунту

?

1.4038267

код для вставкиСкачать
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
Convergence Properties of the Kalman Inverse
Filter
K. Soltani Naveh
Graduate Student
School of Mechanical
and Mining Engineering
The University of Queensland
Brisbane, Australia
Email: k.soltaninaveh@uq.edu.au
sc
rip
tN
ot
Co
P. R. McAree ?
Professor
School of Mechanical
and Mining Engineering
The University of Queensland
Brisbane, Australia
Email: p.mcaree@uq.edu.au
py
ed
ite
d
B. S. Petschel
Post-doctoral Research Fellow
School of Mechanical and
Mining Engineering
The University of Queensland
Brisbane, Australia
ABSTRACT
Ac
ce
pt
ed
Ma
nu
The Kalman filter has a long history of use in input deconvolution where it is desired to estimate structured
inputs or disturbances to a plant from noisy output measurements. However, little attention has been given to the
the convergence properties of the deconvolved signal, in particular the conditions needed to estimate inputs and
disturbances with zero bias. The paper draws on ideas from linear systems theory to understand the convergence
properties of the Kalman filter when used for input deconvolution. The main result of the paper is to show that,
in general, unbiased estimation of inputs using a Kalman filter requires both an exact model of the plant and an
internal model of the input signal. We show that for unbiased estimation, an identified sub-block of the Kalman
filter that we term the plant model input generator must span all possible inputs to the plant and that the robustness
of the estimator with respect to errors in model parameters depends on the eigen-structure of this sub-block. We
give estimates of the bias on the estimated inputs/disturbances when the model is in error. The results of this paper
provide insightful guidance in the design of Kalman filters for input deconvolution.
? Corresponding
DS-13-1302
author.
1
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
by ASME
List2017
of Figures
3
4
5
6
7
8
Kalman Inverse Filter block diagram. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Kalman Inverse Filter, including the open-loop system governed by A1 and A2 , and the closed-loop
feedback system governed by A?1 and A3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The double mass-spring-damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Response of PIG and PMIG for step disturbance in double mass-spring-damper system. . . . . . . . . . . .
Response of PIG and PMIG for combined sinusoids disturbance in double mass-spring-damper system. . .
The block diagram of the mass spring damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Response of PIG and PMIG for combined step and sinusoid disturbance in double mass-spring-damper
system with imperfect model parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Bias and bias upper bound for imperfect plant model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ed
ite
d
1
2
3
4
18
19
21
23
26
27
Co
py
List of Tables
1
Parameters of the double mass-spring-damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
2
Parameters of the real system versus the system model for the double mass-spring-damper shown in Figure 3. 24
Introduction
Inverse problems in systems theory involve reconstructing the inputs to a system from measured responses. The adaptation of the Kalman filter to inverse problems, sometimes termed the Kalman inverse filter (KIF), was first proposed by
Bayless and Brigham [1] as a framework for removing sensor reverberation from seismic signals. The idea is to augment
the states of the plant model by additional dynamics that captures known structure in the inputs to be estimated. This is
equivalent to assuming the autocorrelation function of the unknown inputs or distrubances that are to be estimated.
The formulation of the KIF in [1] was for continuous-time systems. A discrete-time formulation followed shortly after
in [2]. Most of the applications of the KIF over the last thirty years have been for geophysical signal processing, see for
example [3?5], but the method has also attracted interest for input and disturbance identification in mechanical systems.
Applications include estimation of aerodynamic loads on aircraft [6, 7], friction in machines [8?10], and the estimation of
machine-environment interaction forces [11?14].
The KIF is closely linked to other methods for input identification, e.g. the regularization method of Tikhonov [15].
Work by De Nicolao and Ferrari-Trecate [16, 17] has shown that under appropriate assumptions, Tikhonov regularization
coincides with Bayesian estimation where the autocorrelation function of the unknown input is specified and the plant is a
linear system driven by white noise. De Nicolao and Ferrari-Trecate [16] show that the inverse problem can be solved by
application of the Kalman filter, revealing an equivalence between Tikhonov regularization and the KIF.
As a method for input identification the KIF has a number of attractions including: (i) that it is recursive and therefore
suitable for real-time implementations; (ii) the underlying Kalman filter structure accommodates model and measurement
uncertainty; and (iii) the method can be adapted to non-linear systems through the extended Kalman filter. Limitations of
the method are: (i) that the structure of the inputs/disturbances must be specified, and (ii) the bandwidth to which input/disturbances can be reconstructed depends on the dynamics of the plant, and the process and measurement uncertainties, as
specified in the Kalman filter.
The structure of the KIF is illustrated in Fig. 1, and has four components: the Plant, the Plant Input Generator (PIG),
the Plant Model and the Plant Model Input Generator (PMIG). The Plant is the system whose inputs are to be estimated.
The PIG represents the process that generates these to-be-estimated inputs/disturbances acting on plant. The KIF comprises
the Plant Model and the PMIG. The PMIG assumes the inputs have structure. The idea is to design the PMIG so that its
outputs serve as estimates of the input to the Plant. This is done through the feedback structure of the Kalman filter which
provides an implicit inverse model of the Plant.
This paper establishes the conditions needed for convergence of the PMIG output to the Plant input in steady state for
linear plants. Specifically, the paper establishes the conditions that must be satisfied by the Plant Model and the PMIG for
convergence of the plant input/disturbance estimates to their actual values. Section 2 describes the structure of the Kalman
Ac
ce
pt
ed
Ma
nu
sc
rip
tN
ot
1
DS-13-1302
2
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME real system
Plant Input Generator
(PIG)
Plant
plant inputs
(unknown)
+
estimation error
ed
ite
d
estimated plant inputs
system model/estimator
-
Plant Model Input Generator
(PMIG)
+
py
Plant Model
Co
Kalman gain
ot
Fig. 1: Kalman Inverse Filter block diagram.
The Kalman Inverse Filter
The KIF is described by the block diagram in Fig. 2 and governed by the equations
Ac
ce
pt
ed
2
Ma
nu
sc
rip
tN
Inverse Filter and details the notation used throughout the paper. Section 3 summarizes key assumptions. The main work
of the paper begins in Section 4 where we give general conditions for an asymptotically unbiased KIF by formulating the
problem in terms of the solution to a constrained linear matrix equation. It is shown that if the model of the Plant and PMIG
are exact, the estimates produced by the KIF converge to their true values. Section 5 presents a more detailed investigation
of the requirements on the PMIG and shows that for unbiased estimation it must span all possible inputs to the plant and that
the robustness of the estimator with respect to errors in parameters depends on the eigen-structure of the PIG. It is shown,
moreover, that a robust PMIG can be constructed from knowledge of just the number of outputs and the minimal polynomial
of the PIG dynamics. Inter alia, the analysis also yields a simple proof of the Internal Model Principle which was first proven
in [18]. Section 6 gives error bounds on estimates given knowledge of model errors. The results of this section show that the
estimated bias is approximately proportional to the parameter error and that relative bias is bounded.
DS-13-1302
x?1 = A1 x1 + B1 u + w1
(1a)
x?2 = A2 x2 + w2
(1b)
y = C1 x1 + v
(1c)
u = C2 x2
x?? 1 = A?1 x?1 + B?1 u? + Kb e
(1d)
x?3 = A3 x3 + Kc e
(1f)
(1e)
y? = C?1 x?1
(1g)
u? = C3 x3
(1h)
e = y ? y?
(1i)
z = u ? u?
(1j)
3
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017
ASME
where
thebyvector
x is the state of the Plant, x is the state of the Plant Input Generator (PIG), u is the unknown plant input
1
2
produced by the PIG and y is the measured output signal from the plant. The vector x?1 is the state of the plant model, u? is
the estimate of the input which is produced by the Plant Model Input Generator (PMIG) with state vector x3 generally not
the same size as x2 , and e is the innovation sequence available for feedback. The matrices A1 , A2 , B1 , C1 , C2 , A?1 , B?1 , C?1 , A3 ,
C3 have appropriate sizes satisfying the Eqns. (1). In particular the PMIG and the PIG have the same number of outputs and
the Plant and the Plant Model have the same number of outputs.
Plant
PMIG
Plant Model
Ma
nu
sc
rip
tN
ot
Co
py
ed
ite
d
PIG
ed
Fig. 2: The Kalman Inverse Filter, including the open-loop system governed by A1 and A2 , and the closed-loop feedback
system governed by A?1 and A3 .
Ac
ce
pt
We note, however, that the Plant Model and PMIG may contain additional dynamics not present in the Plant and PIG
under investigation or vice-versa. Hence the model state vectors x?1 , x3 may be of different sizes to the corresponding plant
w (t) state vectors x1 , x2 . The measurement noise v and the process noise w(t) = w1 (t) are vector white noises with covariance
2
structure
E[w(s)wT (t)] = Q?(s ? t)
E[v(s)vT (t)] = R?(s ? t)
E[v(s)wT (t)] = 0.
The system can be thought of as regulating z, the difference between the input signal and the estimated input signal. When
E[z] = 0 the estimates u? track the actual inputs to the Plant u.
DS-13-1302
4
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
ASME
The by
state
vector of the combined system encompassing the Plant and the PIG is
" #
x1
x=
.
x2
" #
x?1
x? =
,
x3
which shall form the model state. The state equations (1) combine to form the system
y = CP x
(3b)
(3c)
(3d)
Co
y? = CM x?
(3a)
py
x? = AP x + w,
x?? = AL x? + BP x + Kv,
ed
ite
d
The state vector of the model of the system, encompassing the Plant Model and the PMIG, is
(3e)
ot
z = DP x ? DM x?,
tN
where
"
#
A?1 B?1C3
AM =
0 A3
sc
rip
"
#
A1 B1C2
AP =
,
0 A2
AL = AM ? KCM ,
h
i
CP = C1 0 ,
h
i
DP = 0 C2 ,
" #
Kb
K=
,
Kc
BP = KCP
h
i
CM = C?1 0 ,
h
i
DM = 0 C3
Ma
nu
(4)
T ?1
K = PCM
R
pt
ed
The so-called Kalman gain K is given by
T ?1
R CM P = 0.
AM P + PATM + Q ? PCM
Ac
ce
where P satisfies the Algebraic Riccati Equation [19]
In order to regulate z, it is necessary to have closed loop stability, that is
?(AL ) ? C? ,
where ?(AL ) is the eigenvalue spectrum of the matrix AL and C? is the open left half of the complex plane. A necessary
condition for existence of a stabilizing gain K is that the pair (AM ,CM ) is detectable. We assume that the pair (AM ,CM ) is
observable; this guarantees existence of a stabilizing gain K such that AP and AL have no common eigenvalues.
DS-13-1302
5
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
by ASME
3 2017
Assumptions
The following assumptions are made at various places and referred to by number as required. They are gathered here to
avoid repetition. As indicated, some of the assumptions may be replaced with less restrictive alternatives.
ed
ite
d
A1 The plant system is observable, i.e. (AP ,CP ) is observable.
A2a All eigenvalues of the plant system matrix AP are unstable.
A2b The plant system matrix AP and the closed loop system matrix AL have no common eigenvalues.
A3 The system model is observable, i.e. (AM ,CM ) is observable.
A4 All eigenvalues of the closed loop matrix AL are stable.
A5a The Plant Model exactly describes the Plant, i.e. (A?1 , B?1 , C?1 ) = (A1 , B1 ,C1 ).
A5b The Plant Model has the same transfer function as the Plant.
A6 The Plant matrix A1 has no eigenvalues in common with the Plant Input Generator matrix A2 .
py
Assumption A4 is essential to the results of this paper. Assumption A3 implies the existence of stabilizing gains K such
that Assumption A4 is satisfied.
Convergence properties of the KIF
This section gives general conditions for an asymptotically unbiased KIF, namely E[z(?)] = 0. The approach taken is
to reformulate the problem in terms of the solution to a constrained linear matrix equation by decomposing the trajectories
into terms arising from the state and the state estimation error. Then using Assumption A4, the conditions for E[z(?)] = 0
are derived.
tN
General conditions for unbiased input estimation
Consider the state equations (3), then x and x? are related by
sc
rip
4.1
ot
Co
4
x? = ? + Xx
(5)
Ma
nu
where X has appropriate dimensions and ? represents the state estimation error. Differentiating both sides with respect to
time and using state equations (3) we can write
??? = x?? ? X x?
= AL x? + BP x + Kv ? X(AP x + w)
= AL (?? + Xx) + BP x + Kv ? X(AP x + w)
ed
= AL ? + (AL X + BP ? XAP )x + Kv ? Xw.
XAP ? AL X = BP
(6)
Ac
ce
pt
Under Assumption A4, AL is stable and hence AL ? converges to zero as t ? ?. Therefore, in order to satisfy limt?? ? = 0,
the term due to x must be zero, so
must hold. This relation is the well-known Sylvester matrix equation and has unique solution X under Assumption A2b that
is if AP and AL have no common eigenvalues [20?22]. Assumption A2b can be satisfied in a number of ways. One sufficient
condition is that Assumptions A2a and A4 hold.
Using the above relations, the input estimation error z can then be expressed in terms of the state and state estimation
error. Specifically
z = DP x ? DM x?
= (DP ? DM X)x ? DM ? .
DS-13-1302
6
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017
ASME
Since
E[??by
]?
0 the term to due to ? converges to 0 . Therefore in order for E[z] ? 0 , the component due to x must be zero,
leading to the following results:
Lemma 1 (Sufficient conditions for zero bias). Under Assumptions A2b and A4, the bias is asymptotically zero if
DM X = DP .
(7)
where X is the solution to the matrix equation (6).
ed
ite
d
Lemma 2 (Necessary conditions for zero bias). Under Assumptions A2b and A4, the bias is asymptotically zero if and only
if
(DP ? DM X)U = 0
(8)
where X is the solution to the matrix equation (6) and U is a matrix whose columns span the unstable subspace of AP .
Co
py
Note that when all eigenvalues of AP are stable, U = 0 and (8) holds regardless of the accuracy of the model. The KIF is
asymptotically unbiased because both the input signal and the input estimate converge to zero. Conversely, if all eigenvalues
of AP are unstable (Assumption A2a), U = I in (8) which reduces to to (7).
Exact Plant Model, Exact PMIG
Before discussing the question of robustness with respect to parameters, we remark that for an exact model the bias is
zero, i.e. for
tN
ot
4.2
sc
rip
(A?1 , A3 , B?1 , C?1 ,C3 ) = (A1 , A2 , B1 ,C1 ,C2 ),
A?1 = T1 A1 T1?1 ,
A3 = T2 A2 T2?1 ,
Ma
nu
the equations (6) and (7) are satisfied by X = I.
In fact, the choice of coordinates for the model state space does not affect the bias: given invertible matrices T1 and T2 ,
the model described by
C?1 = C1 T1?1 ,
B?1 = T1 B1 ,
C3 = T1C2 T2?1 ,
"
#
T1 0
X=
.
0 T2
pt
ed
has zero bias, with equations (6) and (7) satisfied by
Ac
ce
Here the converged state estimates are related to the true states by
x?1 = T1 x1
x3 = T2 x2 .
4.3
Zero-bias equations
To show that for unbiased input estimation, the plant model must be similar to the plant, partition the matrix X into four
blocks of compatible size, with
"
#
X11 X12
X=
X21 X22
DS-13-1302
7
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
by equation
ASME (6) becomes
The2017
matrix
X11 A1 = (A?1 ? KbC?1 )X11 + B?1C3 X21 + KbC1
(9a)
X21 A1 = ?KcC?1 X11 + A3 X21 + KcC1
(9b)
X11 B1C2 + X12 A2 = (A?1 ? KbC?1 )X12 + B?1C3 X22
(9c)
X21 B1C2 + X22 A2 = ?KcC?1 X12 + A3 X22 ,
(9d)
ed
ite
d
and the zero bias equation (7) becomes
C3 X21 = 0
C3 X22 = C2 .
(10a)
(10b)
Co
py
Any given Plant and PIG has a very large class of models for which (9) and (10) can be solved for X. We limit our attention
to models for which the Plant Model has the same dimensionality as the true Plant. In this situation the Plant and the Plant
Model must necessarily obey a similarity relationship. The argument is as follows.
After substitution of (10a), (9a) can be rearranged as
ot
A?1 X11 ? X11 A1 = Kb (C?1 X11 ?C1 ).
tN
In order for this equation to hold regardless of the feedback gain Kb , the left and right sides must be zero, giving
A?1 X11 = X11 A1 ,
sc
rip
C?1 X11 = C1 .
Similarly, (9c) can be reduced to
Ma
nu
(X11 B1 ? B?1 )C2 = (A?1 ? KbC?1 )X12 ? X12 A2 ,
which in order to hold regardless of the signal generator gain C2 necessitates
ed
X11 B1 ? B?1 = 0.
pt
Since the plant (A1 ,C1 ) is observable and the model has the same number of dimensions, X11 must be square, invertible and
unique (see Sec. 5), leading to the following result:
Ac
ce
Theorem 1. Suppose the Plant (A1 , B1 ,C1 ) is observable and controllable, and that the associated zero bias Plant Model
(A?1 , B?1 , C?1 ) has the same number of state space dimensions. Then the Plant and the Plant Model must be related by a unique
similarity transform
?1
A?1 = X11 A1 X11
,
?1
C?1 = C1 X11
,
B?1 = X11 B1 .
Moreover the Plant and the Plant Model have the same transfer function
C?1 (s ? A?1 )?1 B?1 = C1 (s ? A1 )?1 B1
for all s ? C.
DS-13-1302
8
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
ASME
Thisby
result
implies a requirement for exact knowledge of the Plant transfer function: to exactly determine the plant
inputs from the plant outputs requires the exact relationship between the inputs and the outputs.
In general this requirement implies a need for a perfect Plant Model including exact parameter knowledge. An interesting
exception is where all eigen-values of A1 are zero corresponding, for example, to free rigid-body dynamics problems. In this
situation the eigen-structure is independent of the parameter values. The model used by Siegrist [13] to estimate road tyre
interaction forces for off-highway mining vehicles is in this class of problems.
Conditions for unbiased input estimation with a perfect plant model
For algebraic simplicity, we will henceforth assume a perfect Plant Model with A?1 = A1 , B?1 = B1 , C?1 = C1 so that
X11 = I. We now seek conditions for unbiased input estimation under this assumption. A sufficient condition for zero bias is
the existence of a solution to the constrained Sylvester equation
ed
ite
d
4.4
C3 X22 = C2
(11b)
py
A3 X22 = X22 A2 .
(11a)
Co
The solvability of (11) becomes a necessary condition for zero bias when the signal generator has no modes in common with
the plant, as demonstrated in the following result.
tN
ot
Theorem 2. Under Assumptions A2b, A4, A5a, suppose that the PMIG matrix A3 has no eigenvalues in common with the
Plant Model matrix A?1 and that Kb has been chosen so that (A?1 ? KbC?1 ) has no eigenvalues in common with A1 or A2 (this is
possible because (A1 ,C1 ) is observable). The KIF has zero bias if and only if there exists a solution X22 to the simultaneous
matrix equations (11).
sc
rip
Proof. To prove sufficiency, let X22 be a solution to (11). Then the equations (9), (10) are uniquely satisfied by
"
#
I 0
X=
.
0 X22
Ma
nu
To prove necessity, suppose the KIF has zero bias, i.e. (10) holds. Then (9a) reduces to
X11 A1 = (A?1 ? KbC?1 )X11 + KbC1 .
X21 A1 = A3 X21
X12 A2 = (A?1 ? KbC?1 )X12
ce
pt
ed
Since the plant model (A?1 , C?1 ) is exact and A1 has no eigenvalues in common with (A?1 ? KbC?1 ), the solution X11 = I is
unique. The next two equations (9b) and (9c) reduce to
Ac
which, by the assumption of no common eigenvalues, are uniquely satisfied by X21 = 0 and X12 = 0. The other equation (9d)
reduces to
X22 A2 = A3 X22
as required.
We now show that equations (11) imply a relationship between the observability matrices of the PMIG and the PIG.
Suppose that equations (11) hold for the matrix X22 . Multiply (11a) on the right by A2 and substitute (11b) to obtain
C3 A3 X22 = C2 A2 .
DS-13-1302
9
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 bythe
ASME
Repeating
process will yield a sequence of equations
j
j
j ? 0.
C3 A3 X22 = C2 A2 ,
Terminating the sequence at j = k ? 1,1 the system can be expressed as
Ok (A3 ,C3 )X22 = Ok (A2 ,C2 ),
(12)
ed
ite
d
where Ok (A,C) is the kth order observability matrix defined by
?
C
?
?
? CA ?
?
Ok (A,C) ? ? . ?
?.
? .. ?
CAk?1
py
?
Co
When k = n, where A is n О n, the matrix O (A,C) is the observability matrix of control theory.
The noise-free component of an input signal produced by state x2 of the signal generator satisfies
ot
u = C2 x2 ,
tN
u? = C2 x?2 = C2 A2 x2 ,
u? = C2 x?2 = C2 A22 x2 ,
sc
rip
..
.
du
dx2
= C2 k?1 = C2 Ak?1
2 x2 .
dt k?1
dt
Ma
nu
Similarly the input estimate produced by state x3 satisfies
u? = C3 x3 ,
u?? = C3 x?3 = C3 A3 x3 ,
pt
ed
u?е = C3 x?3 = C3 A23 x3 ,
..
.
d u?
dx3
= C3 k?1 = C3 Ak?1
3 x3 .
dt k?1
dt
ce
These relations can be summarized in two systems of equations
u
u?
u?
..
.
Ac
?
?
?
?
?
?
?
?
?
1 As
du
dt k?1
?
?
?
?
?
?
? = O (A2 ,C2 )x2 ,
?
?
?
?
?
?
?
?
?
?
?
u?
u??
u?е
..
.
d u?
dt k?1
?
?
?
?
?
? = O (A3 ,C3 )x3 .
?
?
?
a consequence of the Cayley-Hamilton theorem of linear algebra, the sequence can be stopped at k = n2 n3 , where A2 is n2 О n2 and A3 is n3 О n3 ;
in the next section we show that k = n3 + 1 suffices.
DS-13-1302
10
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017isby
ASMEand x = X x , the PIG produces an identical signal as shown by the following relation.
If (12)
satisfied
3
22 2
?
?
?
?
?
?
?
?
?
u
u?
u?
..
.
du
dt k?1
?
?
?
?
?
?
?
?
?
?
? = O (A3 ,C3 )X22 x2 = O (A3 ,C3 )x3 = ?
?
?
?
?
?
?
u?
u??
u?е
..
.
d u?
dt k?1
?
?
?
?
?
?.
?
?
?
tN
ot
Co
py
ed
ite
d
The interpretation we give to this result is that the signals generated by the PMIG must span all possible inputs to the
plant, so that if, for example, the plant input contains a sinusoidal component of frequency ?, the PMIG must include a mode
with this same frequency. Furthermore, the matrix X22 is a map from the states of the PIG to the states of the PMIG so that
the inputs to the model correspond to those of the plant in a way that is consistent with the dynamics of the two systems.
The above results provide a test of whether a given PMIG will produce unbiased input estimates for a given PIG.
No assumptions have been made about the size of the PMIG relative to the PIG. The PMIG may have additional signal
components that do not appear in the PIG, or multiple copies of components of the signal generator. The PMIG can therefore
contain a library of modes whose form corresponding to all expected inputs to be applied to the plant. By way of example,
if it is known that the plant can subjected to inputs/disturbances having n distinct frequencies, ?1 , . . ., ?n , applied alone or
in combination, the PMIG should be designed to have eigenvalues corresponding to those frequencies.
The next section looks in more detail at the relationship between the structure of the signal generator and the structure
of models that produce zero-bias signal estimates. The analysis shows that the robustness of the estimator with respect to
certain parameters depends largely on the eigen-structure of the PIG. This amounts to determining the general structure of
pairs of matrices (A3 ,C3 ) and (A2 ,C2 ) which satisfy Equation (11) for some X22 .
Conditions on the PMIG for zero bias estimates
In this section, we consider the abstract properties of the relationship between the PIG and the PMIG implied by the
zero bias equations (11), culminating in the result that a robust PMIG can be constructed from knowledge of just the number
of outputs and the minimal polynomial of the PIG dynamics. For ease of argument, we introduce the following notation:
sc
rip
5
Ma
nu
Definition 1. For any matrices A,C, B, D of size n О n, k О n, m О m, k О m respectively, say (A,C) models (B, D), denoted
(A,C) & (B, D), if
(13a)
AX = XB,
(13b)
ed
CX = D
pt
for some matrix X of size n О m. Say (A,C) &X1 (B, D) if (11) holds for X = X1 . In particular, the PMIG will provide zero
bias estimates for a given PIG if (A3 ,C3 ) & (A2 ,C2 ).
ce
The relation has useful ordering properties:
Ac
Reflexive: (A,C) &In (A,C), where In is the n О n identity matrix. This means any plant is a model of itself.
Transitive: (A1 ,C1 ) &X1 (A2 ,C2 ) &X2 (A3 ,C3 ) implies (A1 ,C1 ) &X1 X2 (A3 ,C3 ). This means any model of a model of the
plant is, in turn, a model of the plant.
From these properties, two equivalence relations naturally arise:
Similarity: Say (A,C) ? (B, D) if (A,C) &X (B, D) where the matrix X is square and invertible. This is similarity in the
usual linear systems sense.
Equivalence: Say (A,C) h (B, D) if (A,C) & (B, D) & (A,C). Similarity implies equivalence, but not necessarily vice-versa.
Pairs of models can be related in three fundamental ways (subject, of course, to compatibility of the size of the sub
blocks):
DS-13-1302
11
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
by ASME
R12017
(change
of basis): "
(P?1 AP,CP)
&!
(Q?1 AQ,CQ) by X"=#P?1 Q
#
h
i
A B1
I
R2 (bigger model):
, C D & (A,C) by X =
0 B2
0
"
#
!
h
i
h i
A 0
R3 (unobservable signals): (A,C) &
, C0
by X = I 0
B3 B4
ed
ite
d
Relation R1 says that the choice of coordinate system is irrelevant, relation R2 says that the PMIG can contain additional
dynamics not present in the PIG. Relation R3 says that unobservable components of the PIG can be ignored.
In fact these relations characterize &, since any matrix X that relates (A,C) and (B, D) can be reduced by row and
column operations to echelon form by P?1 XQ = 0I 00 for some invertible matrices P and Q2 . This is summarized in the
following result:
Theorem 3. (A,C) & (B, D) if and only if there are change of basis matrices P and Q that reduce (A,C) and (B, D) to upperand lower- block triangular form respectively such that
(Q?1 BQ, DQ) =
#
!
i
A12 h
, C1 C2
A21
#
!
h
i
0
, C1 0 .
B22
py
(P AP,CP) =
"
A11
0
"
A11
B21
Co
?1
Ma
nu
Lemma 3. If (A,C) &X (B, D), then
sc
rip
tN
ot
The above result characterizes the perfect plant model as containing a subsystem that is equivalent to the observable
part of the plant. In terms of structural stability, this says that the model is insensitive to perturbations that do not change the
equivalence of this subsystem that models the plant.
The aim of the remainder of this section is to show that any observable model (A,C) with suitable replication of the
eigenstructure of B is a perfect model of (B, D), i.e. contains a subsystem equivalent to (B, D). The implication of this result
is that the model is insensitive to any perturbations that do not alter the eigenstructure or compromise the observability of
the model.
In order to achieve this result, we re-examine the results relating the observability matrices of the systems. In the
previous section, we demonstrated that (11) implies (12), and hence
Ok (A,C) ? X = Ok (B, D)
for all k ? 1.
(14)
pt
ed
The reverse is not necessarily true, however note that (14) implies CAi (AX ? XB) = 0 for i ? k ? 2, leading to the
following result:
On (A,C) ? (AX ? XB) = 0.
(15)
Ac
ce
Lemma 4. If (14) holds for k = n + 1, then (AX ? XB) satisfies
Moreover if (A,C) is observable (i.e. On (A,C) has rank n), AX ? XB = 0 and hence (A,C) & (B, D).
Several results immediately follow:
Corollary 3.1.
2 Note
that some of the sub blocks may have zero size.
DS-13-1302
12
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
by ASME
If (A,C)
& (B, D) and (A,C) is observable, such an X is unique.
X
If (A,C) &X (B, D) and (A,C) is observable, any matrices E, F satisfying E ? X = F are of the form E = Q ? O (A,C), F =
Q ? O (B, D) for some matrix Q.
If (A,C) &X (B, D) and (B, D) is observable, then X has full column rank.
If (A,C) &X (B, D) and (B, D) is observable and n = m, then X is invertible (hence (A,C) ? (B, D) and (A,C) is
observable).
If (A,C) h (B, D) and both are observable, then n = m and (A,C) ? (B, D).
ed
ite
d
It is well known that any two single output systems (A, cT ), (B, dT ) are equivalent if they are observable and A and B are
similar (e.g. Casti [23]). A slightly more general result follows:
Corollary 3.2. Suppose A ? B, and cT , dT are any two row vectors such that (A, cT ) is observable (single-output case).
Then (A, cT ) &X (B, dT ) with X = (On (A, cT ))?1 On (B, dT ).
py
Proof. Since (A, cT ) is observable, On (A, cT ) is n О n and invertible. Since A and B have the same characteristic polynomial,
if (14) holds for k = n it must hold for k = n + 1 (the final equation is the same linear combination of the previous n equations)
and the result follows.
ot
Co
The following result shows that every model of an output system contains a subsystem equivalent to that system.
"
#
!
A1 , A2
Lemma 5. Suppose (A,C) &X (B, D) where (B, D) is observable. Then (A,C) ?
, [C1 ,C2 ] for some A1 , A2 , A4 ,C1 ,C2
0, A4
such that (A1 ,C1 ) ? (B, D).
tN
Proof. Let (A,C) &X (B, D). Since (B, D) is observable, X has full column rank by corollary 3.1. Then
sc
rip
" #
I
Q?1 ,
X =P
0
where P and Q are invertible matrices. Hence by Eq.(13)
Ma
nu
"
#
!
i
A1 , A2 ?1 h
P
P , C1 , C2 P?1 ,
0, A4
(A,C) =
where A1 = Q?1 BQ and C1 = DQ, and the result follows
"
#
A1 0
A1 ? A2 =
,
0 A2
ce
pt
ed
Larger systems can be built up from sums of smaller models in parallel: define the direct sum of two matrices as
and we have the following result:
Ac
Lemma 6. If (A1 ,C1 ) & (B1 , D1 ) and (A2 ,C2 ) & (B2 , D2 ), then (A1 ? A2 , [C1 ,C2 ]) & (B1 ? B2 , [D1 , D2 ]) (provided the matrices C1 and C2 have the same number of rows) and (A1 ? A2 ,C1 ?C2 ) & (B1 ? B2 , D1 ? D2 ).
A generalization of this result to k-output systems forms a basis for the internal model principle:
Lemma 7. Suppose A is in block-diagonal form with k identical blocks, i.e.
?
?
A
0
1
k
M
? .
?
.. ?,
A=
A1 = ?
?
?
0
A1
DS-13-1302
13
(16)
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017
where
A by
is ASME
n О n , and
1
1
1
?
?
cT1
0
?
? .
.. ?,
C=?
?
?
T
0
ck
?
dT11
? .
.
D=?
? .
dTk1
?
dT1k
.. ?
?
. ?,
и и и dTkk
иии
..
.
ed
ite
d
where cT1 , . . . , cTk , dT11 , . . . , dTkk are row vectors. Then (A,C) & (A, D) if and only if (A1 , cTi ) & (A1 , dTij ) for all i, j, 1 ? i ?
k, 1 ? j ? k.
Proof. Setting
Co
where the sub-blocks Xi j are each of size n1 О n1 , (13) reduces to
tN
ot
A1 Xi j = Xi j A1
cTi Xi j = dTij ,
py
?
?
X11 и и и X1k
? . .
?
. . . ... ? ,
X =?
? .
?
Xk1 и и и Xkk
and the result follows.
sc
rip
We note in Lemma 7 that observability of all subsystems (A1 , cTi ) is equivalent to observability of (A,C). Also by
Corollary 3.1, if (A, D) is observable then the systems are equivalent, so we have the following result:
Ma
nu
Lemma 8. Suppose (A,C1 ), (A,C2 ) are systems where A is of block-diagonal form (16), and C1 ,C2 are k О n matrices. If
(A,C1 ) is observable, (A,C1 ) & (A,C2 ). Moreover, if both systems are observable, they are equivalent.
Proof. All that remains is to show that if (A,C1 ) is observable, there exists a vector cT1 such that (A1 , cT1 ) is observable. Since
A consists of k copies of A1 , they have the same minimal polynomial and hence
ed
rank(On (A,C)) ? km1 ,
Ac
ce
pt
where m1 is the degree of the minimal polynomial of A1 . The matrix has full rank only if m1 = n1 , i.e. A1 is cyclic (also
termed nonderogatory [21]: each eigenvalue has geometric multiplicity 1). For some change of basis matrix P, such a matrix
will have Jordan form
?
J1
?
?1
P A1 P = J = ?
?0
0
?
0 0
.. ?
. 0?
?,
0 Jr
where the Jordan blocks
?
?
?i 1 0
? .
?
?
.
Ji = ?
?0 . 0?
0 0 ?i
DS-13-1302
14
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 by to
ASME
correspond
distinct eigenvalues ? . Set dT = [1, 0, . . . , 0]; it easily follows that (J , dT ) is observable, as is (J, dT ) where
i
i
i
i
dT = [dT1 , . . . , dTk ]. Then set cT1 = dT P?1 , so we have (A1 , c1 ) ? (J, dT ) and hence (A1 , c1 ) is observable.
With this choice of c1 , apply Lemma 7 with cTi = cT1 , 2 ? i ? k to obtain C such that (A,C) & (A,C1 ) and (A,C) & (A,C2 ).
If (A,C1 ) is observable, (A,C1 ) ? (A,C) and hence (A,C1 ) & (A,C2 ). If both (A,C1 ) and (A,C2 ) are observable, they are
both equivalent to (A,C), and hence (A,C1 ) ? (A,C2 ).
This result will be used to show that any system (B, D) can be perfectly modelled by an observable system (A,C), where
A is of the form (16). The matrix A1 is any matrix consisting of the largest cyclic component of B, i.e. a cyclic matrix
satisfying the minimal polynomial of B.
Lk
A1 where A1 is the largest cyclic component of B. Then (A,C) & (B, D)
ed
ite
d
Lemma 9. Suppose (B, D) is observable, and A =
if (A,C) is observable.
Lp
"
#
Bi Bi1
B1 ?
0 Bi2
ot
for some matrices Bi1 , Bi2 , and hence we can augment the model to obtain
Co
py
Proof. Without loss of generality, assume A1 is in Jordan form. By the cyclic decomposition theorem [20, 21], B ? i=1 Bi
where the Bi are cyclic matrices in Jordan form such that B1 = A1 and each Bi satisfies the minimal polynomial of B1 .
Each Bi can be augmented to a matrix similar to B1 , i.e.
tN
"
#
!
i
B B01 h
, D 0 & (B, D),
0 B02
sc
rip
such that
"
#
p
M
B B01
B1
?
0 B02
i=1
Ma
nu
and
B01 =
ed
B02 =
p
M
i=1
p
M
Bi1
Bi2 .
i=1
ce
pt
Since (B, D) is observable, p ? k,3 and so the model can be augmented with k ? p copies of B1 to obtain
Ac
(A,C) &
"
#
!
!
k?p
h
i
M
B B01
?
B1 , D 0 & (B, D),
0 B02
i=1
and the result follows.
As an easy corollary of the previous result and relations R1 and R2, the main result follows, namely that any observable
system that models k copies of the largest cyclic component of B, is a perfect model of (B, D), where k is the number of rows
of C and D. This gives an alternative (and arguably simpler) proof of the Internal Model Principle [18].
3 PROOF:
Suppose p > k. Some eigenvalue ? has geometric multiplicity p, so rank
B??I D
? n + k ? p < n, and hence (B, D) is not observable. Hence
p ? k.
DS-13-1302
15
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 by4 ASME
Theorem
(Internal Model Principle). Suppose D has k rows, B is the largest cyclic component of B, and
1
"L
#
B1 A1
,
0 A2
k
A?
(17)
for some matrices A1 and A2 . Then (A,C) models (B, D) for any matrix C such that (A,C) is observable.
"L
ed
ite
d
Proof. Let X be a change of basis matrix such that X ?1 AX = A0 , where A0 is the right hand side of (17), and so
#
!
i
B1 A1 h
, C1 C2 ,
0 A2
k
(A,C) ?
where [C1 ,C2 ] = CX. By relations R1 and R2, (A,C) & ( k B1 ,C1 ). Since (A,C) is observable, the subsystem ( k B1 ,C1 )
is observable, and by the previous lemma, models (B, D) for any D of compatible size. Hence, by the transitive property,
(A,C) & (B, D).
L
Co
py
L
sc
rip
tN
ot
We have proved that a robust model of an observable plant (and in particular, a robust PMIG for a given PIG) can
be constructed from knowledge of only the minimal polynomial (i.e. the largest cyclic component) of the plant dynamics,
and the number of outputs, by incorporating multiple copies of the largest cyclic component. This model is robust to any
perturbations that do not affect the eigenstructure of the model, or the observability of the model, and moreover it may
contain additional dynamics not present in the underlying plant.
We now prove the converse of the internal model principle, namely that any ?universal? model of an observable plant
necessarily contains a subsystem with multiple copies of the largest cyclic component.
Theorem 5. Suppose (A,C) & (B, D) for all matrices D having the same number of rows as C, and let B1 be the largest
cyclic component of B. Then A is equivalent to a matrix of the form (17).
ed
Ma
nu
"
#
B1 , B2
Proof. Let B1 be the largest cyclic component of B, and let B ?
. As in the proof of Lemma 8, choose a row vector
0, B4
T
dT1 such that (B
Since (A,C) & (B, D) for all D, by the transitivity property and property R2, it follows
?
?1 , d1?) is ?observable.
0
? ? ??
that (A,C) & ?B1 , ?dT1 ??. By Lemma 6, it follows that (A,C) & (?k B1 , ?k dT1 ), and by Lemma 5 it follows that (A,C)
0
contains a subsystem equivalent to (?k B1 , ?k dT1 ).
Structural stability
Any PMIG satisfying (11) (for some X22 ), that is subject to perturbations in some or all of its parameters will continue
to produce unbiased input estimates provided (11) is still solvable. Such a PMIG is said to be structurally stable with respect
to the perturbed parameters. Structurally stable parameters do not need to be known exactly for unbiased input estimates.
The degree to which the solvability of (11) is stable with respect to the parameters of PIG and PMIG depends partly
on the knowledge of the problem and, for multiple input systems, on the number of identical copies of the PIG within the
PMIG.
By way of example, consider a plant having two inputs. If the inputs are sinusoids of the same frequency but different
amplitude and phase, then for unbiased input estimates the PMIG must have two undamped modes at this frequency. This
PMIG does not, however, require a priori specification of the amplitude and phase of each mode.
Ac
5.1
ce
pt
The amount of redundancy required in the model depends entirely on the knowledge of the plant; if the parameters of
the plant dynamics and outputs are known exactly, a single copy of the exact model will suffice.
DS-13-1302
16
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
by ASME
For more
complicated multiple-output signal generators, the canonical form of the system needs to be taken into consideration. As an example, the system
?
?
#
0?0 "
? abc ?
??
?
?? 0 0 0 ? ,
de f
000
??
is equivalent to the system
?
?
#
010 "
? y01 ?
??
?
?? 0 0 0 ? ,
100
000
ed
ite
d
??
1
0
0
0
0
0
0
0
?
?
0 "
#
?
0?
? 1000 ?
?,
?,
1? 0 0 1 0 ?
0
ot
0
?? 0
??
??
?? 0
0
Co
??
py
for some y depending on the parameters ?, a, b, c, d, e, f , and hence a model of this form is not robust with respect to all of
these parameters. The model
? ? ?
0
x?1
?x? ? ? 0
? 2? ?
? ? = ? ?k1 ?k2
?x?1 ? ? m
0
0
?? ? ?
?
0 0 " #
x1
? ?x ? ? 0 0 ? F
? ? 2? ?
? 1
?
c2 ? ? ? + ? 1
?
?
?
?
? F2
0
x?
1
m1
m1
?c2 ?c3
1
x?2
0 m2
m2
? ?
x1
?x ?
? 2?
y = ? ?.
?x?1 ?
x?2
1
0
0
1
Ma
nu
k2
?c2
m1
m1
?k2 ?k3 c2
m2
m2
(18)
ed
x?2
1
k2
m2
sc
rip
tN
which is observable and contains two copies of the maximal cyclic component, is structurally stable with respect to the
parameters ?, a, b, c, d, e, f .
The examples that follow illustrate input estimation and structural stability when the Plant Model and PMIG are exact.
These examples use the forced mass-spring-damper system shown in Figure 3. The system parameters are shown in Table 1.
The dynamics of the system are described by
Ac
ce
pt
Further we assume that the measurement noise and process noise are both zero-mean Gaussians with covariance R and Q
respectively.
After substituing the parameters in Table 1 the plant is represented by A1 , B1 ,C1 matrices shown below.
?
?
0
0
1 0
? 0
0
0 1?
?
?
A1 = ?
?
??10 5 ?1.5 1.5?
5 ?10 1.5 ?3
?
?
0 0
?0 0?
?
?
B1 = ?
?
?10 0 ?
0 10
C1 = I4О4
DS-13-1302
17
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Value
k1
0.5 N/m
k2
0.5 N/m
k3
0.5 N/m
m1
0.1 kg
m2
0.1 kg
c2
0.15 N.s/m
c3
0.15 N.s/m
py
Parameter
ed
ite
d
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
sc
rip
tN
ot
Co
Table 1: Parameters of the double mass-spring-damper system.
Ma
nu
Fig. 3: The double mass-spring-damper system.
h i
A2 = 0
" #
1
C2 =
0.5
ce
pt
ed
Example 1. Step Input Consider the double mass-spring-damper system. Suppose the driving forces are unknown steps,
F2
but the ratio between them is known. Suppose
= 0.5, then PIG is represented using A2 and C2 matrices.
F1
Ac
Note that since x?2 = A2 x2 + w2 and A2 = [0], x2 (0) must be set to a non-zero value in order to generate a step input. In this
example x2 (0) = 1. Assuming exact Plant Model specifying the PMIG so that it has a single eigenvalue at the zero consistent
with the structure of a step input gives.
A?1 = A1
B?1 = B1
C?1 = C1
A3 = A2
C3 = C2 .
In PMIG the initial state is set to zero i.e. x3 = 0 so that the model is unaware of the actual input generated by PIG.
Additionally we set the Q = I5О5 О 10?3 and R = 5 О I4О4 О 10?3 .
Simulating the system response produces the results shown in Figure 4, with the key point to note being that the inputs
converge to their correct values when the system transients have decayed.
DS-13-1302
18
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
3
4
2.5
3
4
3
2
2
1.5
2
1
1
0
0.5
1
-1
0
0
3
6
9
12
15
-0.5
0
3
6
9
12
15
time (s)
time (s)
1
3
6
9
12
0
15
0
3
6
time (s)
0
0
12
15
0
3
6
Ma
nu
9
time (s)
15
0
0
3
6
9
12
time (s)
ot
0.5
6
12
tN
0.5
sc
rip
1
3
9
time (s)
1
0
12
9
12
15
time (s)
ed
Fig. 4: Response of PIG and PMIG for step disturbance in double mass-spring-damper system.
Ac
ce
pt
Example 2. Compound Sinusoid-step Input Consider the double mass-spring-damper system. Suppose the driving forces
in Figure 3 are a combination of two sinusoids function with known frequencies ?1 , ?2 and a step function. The magnitude
and phase of these sinusoids and the magnitude of the the step function are unknown to the model, but the ratio between
components in F1 and F2 are known. More specifically if the inputs are written as
F1 = a1 u(t) + b1 sin(?1t + ?1 ) + c1 sin(?2t + ?2 )
F2 = a2 u(t) + b2 sin(?1t + ?1 ) + c2 sin(?2t + ?2 )
then
a2 b2 c2
, , are known.
a1 b1 c1
DS-13-1302
15
0.5
Co
0
9
py
0.5
-1
-2
6
1
1
0
3
1.5
3
2
0
time (s)
1.5
4
-2
ed
ite
d
0
19
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
15
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
by example
ASME suppose that PIG is represented by A and C matrices shown below.
In this
2
2
ed
ite
d
?
?
0 0 0 0 0
?
?
?0 0 ?1 0 0 ?
?
?
?
A2 = ?
?0 ??1 0 0 0 ?
?
?
?0 0 0 0 ?2 ?
0 0 0 ??2 0
"
#
1 0.5 0.5 1 ?1
C2 =
0.5 0.5 ?0.5 1 1
A?1 = A1
B?1 = B1
py
h
iT
The initial state x2 (0) determines the phase and magnitude of the sinusoids in each input. We set x2 (0) = 0.5 0.25 0.25 0.5 0.5 .
Additionally we choose f1 = 1Hz, f2 = 0.2Hz and set ?1 = 2? f1 , ?2 = 2? f2 .
Assuming exact Plant Model and PMIG gives
C?1 = C1
A3 = A2
C3 = C2 .
tN
ot
Co
In the PMIG the initial state is set to zero i.e. x3 = 0 so that the model is unaware of the actual input generated by PIG.
Additionally we set the Q and R to the values of the previous example.
The the system shown in Figure 2 is simulated and the responses are presented in Figure 5. The estimated inputs are
seen to track the inputs after initial transients have decayed.
Imperfect plant model: error bounds
The results of the previous section showed that the asymptotic un-biasedness of the filter is not generally preserved
under perturbations in parameters relating to the eigenvalues of the plant and signal.
In practical applications, the parameters are never known exactly and the model is likely to capture only an approximation of the dynamics of interest, hence model error must be taken into consideration. This section shows that (i) the estimated
bias is approximately proportional to the parameter error and (ii) the relative bias is bounded.
The following result gives an estimate of the bias under the assumption that the model is ?close? to an exact model.
Ma
nu
sc
rip
6
0 , D0 ) is in error, such that
Theorem 6. Suppose that the system model (A0L ,CM
M
A0L = AL + ?AL ,
D0M = DM + ?DM ,
pt
ed
0
CM
= CM + ?CM ,
Ac
ce
where (AL ,CM , DM ) is a model producing unbiased estimates of the system (AP , BP ,CP , DP ) governed by (3). Then the bias
E[z0 ] of the imperfect system model satisfies
E[z0 ] ? ?(?DM X + (DM + ?DM )?X)E[x]
(19)
where X satisfies (6) and ?X is the unique solution to
?XAP ? AL ?X = ?AL X.
(20)
Proof. The imperfect estimate satisfies
x?? 0 = (AL + ?AL )x?0 + BP x
DS-13-1302
20
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
4
4
3
3
2
2
2
1
1
0
0
0
-1
-1
-2
-2
4
-2
0
3
6
9
12
15
-3
0
3
6
9
12
15
time (s)
time (s)
2
12
-2
-2
15
0
3
6
1
1
0.5
0.5
0
0
6
9
12
time (s)
-0.5
15
0
3
Ma
nu
3
15
0
3
6
9
12
15
time (s)
tN
1.5
12
sc
rip
1.5
0
9
time (s)
time (s)
-0.5
15
py
Co
9
ot
6
12
-1
-1
-4
3
9
0
0
0
6
1
-2
-6
3
2
1
0
0
time (s)
4
2
-6
ed
ite
d
-3
-4
6
9
12
15
time (s)
Fig. 5: Response of PIG and PMIG for combined sinusoids disturbance in double mass-spring-damper system.
z0 = DP x ? (DM + ?DM )x?0 .
pt
ed
and hence z0 satisfies
Ac
ce
For the perfect model, apply Lemma 1 to obtain
DM X = DP ,
where X satisfies
XAP ? AL X = BP .
Let ?X be the solution of
(X + ?X)AP ? (AL + ?AL )(X + ?X) = BP ,
DS-13-1302
21
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017
by ASME
which
simplifies
to (20). Defining ? 0 = x?0 ? (X + ?X)x, we have
???0 = (AL + ?AL )??0 ? (X + ?X)w.
The resulting bias E[z0 ] then satisfies
E[z0 ] = (DP ? (DM + ?DM )(X + ?X))E[x] ? (DM + ?DM )E[??0 ]
ed
ite
d
= ?(?DM X + (DM + ?DM )?X)E[x] ? (DM + ?DM )E[??0 ].
Since ? 0 is stable, E[??0 ] = 0 and the result follows.
Co
py
When the imperfect model contains multiple copies of an approximation of the maximal cyclic component of the system,
the model error can be taken to be entirely due to the error in the eigenstructure, i.e. ?CM = 0, ?DM = 0, ?AL = ?AM , where
?AM is the error of the model relative to a perfect model which contains suitably replicated copies of the true maximal cyclic
component of the system.
The above result, although it requires the solution of the matrix equation (20), does give a measure of the error of the
input estimate, modulo noise, relative to the plant state:
Corollary 6.1. The bias relative to the plant state satisfies
tN
ot
kE[u ? u?]k kE[z]k
=
? k?DM k и kXk + (kDM k + k?DM k) k?Xk .
kE[x]k
kE[x]k
sc
rip
Unfortunately the result does not give an error estimate relative to the input signal; for example if the input signal is
sinusoidal then the relative error can be infinite at the instants when the input is zero.
The next result provides an estimate of the bias when the parameter errors are large.
Ma
nu
Lemma 10. For the plant and model governed by (3) with closed loop stability, the bias E[z] is O(t n eat ) for some a ? 0 and
n ? 0. Moreover, if the plant state is bounded, with kE[x]k ? M, then E[z] is bounded with
kE[z]k ? M kDP k +C1 M kDM k max(1, kBP k /х)
for some positive constants C1 , х depending only on AL .
ed
Proof. The solution to (3) is
pt
x(t) = ?P (t)x(0) +
0 ?P (t ? ?)w(?)d?,
(21)
= exp(tAP ),
x?(t) = ?L (t)x?(0) +
?L (t)
(22)
Rt
0 ?L (t ? ?)BP x(?)d?,
(23)
= exp(tAL ).
(24)
Ac
ce
?P (t)
Rt
The matrix AL is stable, and thus
k?L (t)k ? C1 e?хt
for some positive constants C1 , х. This gives us the bound
kE[x?(t)]k ? C1 e?хt kx(0)k +C1 kBP k
DS-13-1302
22
Z t
e?х(t??) kE[x(?)]k d?.
0
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 by
ASME
It easily
follows
that if kE[x(t)]k is O(t n eat ) as t ?? ? for some a > 0 then kE[x?(t)]k is also O(t n eat ). The bias E[z] is a
linear combination of E[x?] and E[x] and thus of the same order.
For the bounded case kE[x(t)]k ? M, we have
kE[x?(t)]k ? C1 Me?хt +C1 M kBP k (1 ? e?хt )/х
= C1 M kBP k /х + e?хt C1 M(1 ? kBP k /х)
? max(C1 M,C1 M kBP k /х),
ed
ite
d
and applying the triangle inequality to
z = DP x + DM x?,
py
the result follows.
sc
rip
tN
ot
Co
The sensitivities can be difficult to calculate symbolically, so we give a numerical example here.
Ma
nu
Fig. 6: The block diagram of the mass spring damper system.
Example 3. Consider the mass spring damper system in Figure 6 with step-plus-sine and parameters w = 1, m = 1, k = 1, c =
10. Suppose that model parameters are in error, with
k? = k + ?k
c? = c + ?c,
ed
w? = w + ?w
m? = m
ce
pt
where ?w, ?k, ?c are small. For the stabilizing gains
"
Ac
Kb =
? ?
0
? ?
Kc = ?25? ,
50
#
0
100
the tracking error satisfies
ku ? u?k ? (8.77 |?w| + 0.21 |?c| + 1.41 |?k|) kxk .
Example 4. Step-Sinusoid Input with Imperfect Model Consider the double mass-spring-damper system shown in Figure 3. Suppose that the Plant Model and PMIG contain errors as shown in Table 2.
DS-13-1302
23
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
Parameter
Plant
Parameter Plant Model
0.5 N/m
k?1
0.55N/m
k2
0.5 N/m
k?2
0.55N/m
k3
0.5 N/m
k?3
0.55N/m
m1
0.1 kg
m?1
0.11kg
m2
0.1 kg
m?2
0.11kg
c2
0.15 N.s/m
c?2
0.165 N.s/m
c3
0.15 N.s/m
c?3
0.165 N.s/m
Parameter
PIG
Parameter
PMIG
f
1Hz
f?
1.1Hz
py
ed
ite
d
k1
Co
Table 2: Parameters of the real system versus the system model for the double mass-spring-damper shown in Figure 3.
B?1 6= B1
A3 6= A2
tN
A?1 6= A1
ot
Therefore
C?1 = C1 .
sc
rip
The initial state of the PIG, x2 (0), determines the real input. Specifically x2 (0) determines the phase of sinusoid component,
h
iT
the magnitude of the sinusoid component, and the size of step input component. In this example x2 (0) = 0.25 0.25 1 .
By substituting the model parameters into Equation 18 the following matrices are calculated for Plant Model.
?
0
0
1
0
?
? 0
0
0
1
?
?
A?1 = ?
?
??10 5 ?1.5273 1.5273 ?
5 ?10 1.5273 ?3.0545
?
?
0
0
? 0
?
0
?
?
B?1 = ?
?.
?9.0909
?
0
0 ?9.0909
pt
ed
Ma
nu
?
ce
C?1 = C1 = I4О4
Ac
Furthermore PMIG is given by
?
? ?
?
o 2? f? 0
0
6.9115 0
?
? ?
?
A3 = ??2? f? 0 0? = ??6.9115 0 0?
0
0 0
0
0 0
"
#
0.5 0.5 1
C3 = C2 =
.
0.5 ?0.5 0.5
Using the above matrices yields the results shown in Figure 7. To verify that the input estimation bias satisfies CorolDS-13-1302
24
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
lary2017
6.1, by ASME
kE[uu ? u?u]k
kE[xx]k
is plotted at every time step. The supremum of steady state estimate bias
kE[uu ? u?u]k
kuu ? u?uk
= lim
}
t?? kx?
xk
kE[xx]k
ed
ite
d
sup{
is compared with
k?DM k и kXk + (kDM k + k?DM k) k?Xk
kE[[y1 , y2 , y?1 , y?2 , x2 ]]k
ot
and calculate
T
F 1 , F?
F 2 ]T E[[[F 1 , F 2 ] ? [ F?
=
[y?1 , y?2 , y??1 , y??2 , x 3 ]
Co
kE[uu ? u?u]k
=
kE[xx]k
T
F 1 , F?
F 2 ]T E[[[F 1 , F 2 ] ? [ F?
py
. In this example we use
sc
rip
tN
T
F 1 , F?
F 2 ]T [[F 1 , F 2 ] ? [ F?
[y?1 , y?2 , y??1 , y??2 , x 3 ]
for every point obtained from the simulation. Figure 8 shows that as the response converges to steady state,
revealed. The figure shows
kE[uu?u?u]k
kE[xx]k
is
kE[uu ? u?u]k
} = 0.1408.
kE[xx]k
Ma
nu
sup{
Additionally by substituting the matrices into Equation 4 and solving Equation 20 the following can be obtained.
Ac
ce
pt
ed
?DM = 0 2О7
"
#
00001 1 1
DM =
0 0 0 0 1 ?1 1
?
?
?0.1288 0.1349 ?0.0198 0.0376 0.2144 0.0855 0.0348
? 0.0844 ?0.0888 0.0236 ?0.0383 ?0.0416 ?0.2628 ?0.0223?
?
?
?
?
? 0.2623 ?0.3413 ?0.0571 0.0207 0.1697 1.2240 0.0179 ?
?
?
?
?X = ?
??0.4445 0.5172 ?0.0104 0.0645 1.6985 0.0463 0.0513 ? .
?
?
??0.0415 0.0426 ?0.0088 0.0142 0.1291 1.2262 0.0124 ?
?
?
? 0.0384 ?0.0472 ?0.0049 ?0.0002 ?1.2068 0.0869 ?0.0011?
0.0167 ?0.0175 ?0.0017 ?0.0003 ?0.0522 0.0587 0.0863
Therefore k?DM k = 0, kDM k = 1.3229, kXk = 1, k?Xk = 2.1996
k?DM k и kXk + (kDM k + k?DM k) k?Xk = kDM k k?Xk = 2.9098
DS-13-1302
25
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 bythat
ASME
It follows
kE[uu ? u?u]k
? k?DM k и kXk + (kDM k + k?DM k) k?Xk .
kE[xx]k
3
3
2
2
1
1
4
3
ed
ite
d
2
1
0
0
0
-1
-1
-1
-2
3
6
9
12
15
0
3
6
time (s)
9
12
time (s)
1
1
0.8
0
0.6
-1
0.4
-2
0.2
-3
0
3
6
9
12
15
0
1.5
1
0.5
3
Ma
nu
0.5
0
6
9
12
15
12
15
time (s)
1
6
9
12
0.6
0.4
0.2
15
time (s)
time (s)
1
3
0.8
0
0
3
6
9
time (s)
sc
rip
0
0
1.2
ot
1.2
2
tN
3
-3
Co
4
15
py
0
0
-0.5
-0.5
0
3
6
9
12
0
3
6
9
12
15
time (s)
ed
time (s)
15
Conclusion
We have shown that for unbiased input or disturbance estimation using the Kalman Inverse Filter, an exact model of
the plant and a suitably replicated model describing the structure of the input signals is required. Moreover, we have shown
that the model of the input signal generator is structurally stable with respect to all parameters that do not change the eigen
structure of the PMIG. A new proof of the internal model principle has also been given.
These results provide insight to the design of Kalman Inverse Filters. In several of the recent papers using this method,
e.g. [6?13] the structural form of the PMIG is a chain of integrators without clear justification for these choices. The results
of this paper show that this approach amounts to approximating the input signal as a piecewise polynomial whose order
corresponds to the degree of the minimal polynomial of the PMIG matrix. A single integrator corresponds to a system with
DS-13-1302
Ac
7
ce
pt
Fig. 7: Response of PIG and PMIG for combined step and sinusoid disturbance in double mass-spring-damper system with
imperfect model parameters.
26
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
Bias relative to state
10
bias relative to state
upper bound
9
8
7
ed
ite
d
6
5
4
3
py
2
Co
1
0
4
6
8
time (s)
14
bias relative to state versus upper bound calculated by k?DM k и kXk + (kDM k + k?DM k) k?Xk.
tN
kE[uu?u?u]k
kE[xx]k ,
12
sc
rip
Fig. 8:
10
ot
2
Ma
nu
dynamics x? = 0 and thus would be consistent with an input having the structure of a constant (i.e. step) input. Similarly a
chain of two integrators is consistent with the input having the structure of a ramp. For example, in [12], a double integrator
is used to model the signal associated with each of the unknown forces acting on the wheels of a vehicle. Similarly a chain
of three integrators yield a piecewise parabolic fit to an input. De Nicolao [16], in his work on Tikhonov regularization has
similarly observed that chains of integrators can be associated with smoothing splines.
ed
For applications where the inputs have a higher structure that is known, e.g. harmonic inputs and/or disturbances, the
insights provided by this work informs us that the the PMIG should be modelled as a resonator having the appropriate
frequency. Where the signal is not a pure sinusoid, the PMIG might be extended to include harmonics of the base frequency.
The efficacy of this approach has been demonstrated in [14]. If the system has multiple inputs suitable replication of the
input structure is required. We note finally the possibility that where the structure of the input may be time varying, or not
precisely defined, the PMIG could include appropriately rich structure to adapt to the input signal.
8
Ac
ce
pt
From the point of view of the philosophical foundations of the design and in terms of guiding practical implementation
of input-and-disturbance estimators, we have found the principles revealed in this paper extremely powerful in applications
where we seek to estimate inputs and disturbances in systems using the KIF.
Acknowledgements
Kianoosh Soltani Naveh would like to acknowledge The University of Queensland Centennial Scholarship for the support of this work.
References
[1] Bayless, J. W., and Brigham, E. O., 1970. ?Application of the Kalman filter to continuous signal restoration?. Geophysics, 35(1), Feb, pp. 2?23.
DS-13-1302
27
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Ac
ce
pt
ed
Ma
nu
sc
rip
tN
ot
Co
py
ed
ite
d
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)[2]2017
by ASME
Crump,
N. D., 1974. ?A Kalman filter approach to the deconvolution of seismic signals?. Geophysics, 39(1), Feb,
pp. 1?13.
[3] Dimri, V., 1992. Deconvolution and Inverse Theory. Elsevier Publishing Company.
[4] Leite, L. W. B., and da. Rocha, M. P. C., 2000. ?Deconvolution of non-stationary seismic process?. Rev. Bras. Geof.,
18(1), pp. 75?89.
[5] Rocha, M. P., and Leite, L. W., 2003. ?Treatment of geophysical data as a non-stationary process.?. Mat. Apl. Comput.,
22(2), pp. 149?166.
[6] Stalford, H. L., 1981. ?High-alpha aerodynamic model identification of t-2c aircraft using the ebm method?. Journal
of Aircraft, 18(10), pp. 801?809.
[7] Sri-Jayantha, M., and Stengel, R. F., 1988. ?Determination of nonlinear aerodynamic coefficients using the estimationbefore-modeling method?. Journal of Aircraft, September, pp. 796?804.
[8] Ray, L. R., Ramasubramanian, A., and Townsend, J., 2001. ?Adaptive friction compensation using extended KalmanBucy filter friction estimation?. Control Engineering Practice, 9(2), pp. 169?179.
[9] Ray, L. R., Townsend, J., and Ramasubramanian, A., 2001. ?Optimal filtering and Bayesian detection for friction-based
diagnostics in machines?. ISA Transactions, 40(3), pp. 207?221.
[10] Austin, K. J., and McAree, P. R., 2007. ?Transmission friction in an electric mining shovel?. Journal of Field Robotics,
24(10), pp. 863?875.
[11] Ray, L. R., 1995. ?Nonlinear state and tire force estimation for advanced vehicle control?. IEEE Transactions on
Control Systems Technology, 3(1), pp. 117?124.
[12] Ray, L. R., 1997. ?Nonlinear tire force estimation and road friction identification: Simulation and experiments?.
Automatica, 33(10), pp. 1819?1833.
[13] Siegrist, P., and McAree, P. R., 2006. ?Tyre-force estimation by kalman inverse filtering: applications to off-highway
mining trucks?. Vehicle system dynamics, 44(12), pp. 921?937.
[14] Reid, A. W., McAree, P. R., Meehan, P., and Gurgenci, H., 2014. ?Longwall shearer cutting force estimation?. Journal
of Dynamic Systems, Measurement, and Control, 136(3).
[15] Tikhonov, A. N., and Arsenin, V. Y., 1977. Solutions of ill-posed problems. Washington : Winston.
[16] G., D. N., and G., F.-T., 2001. ?Regularization networks: Fast weight calculation via kalman filtering?. IEEE Transactions On Neural Networks, 12(2), March, pp. 228?235.
[17] Nicolao, G. D., and Ferrari-Trecate, G., 2003. ?Regularization networks for inverse problems: A state-space approach?.
Automatica, 39(4), April, pp. 669?676.
[18] Francis, B. A., and Wonham, W., 1976. ?The internal model principle of control theory?. Automatica, 12(5), pp. 457?
465.
[19] Gelb, A., 1974. Applied optimal estimation. MIT press.
[20] Gantmacher, F. R., 1959. The Theory of Matrices, Vol. 1. Chelsea Publishing Company.
[21] Horn, R. A., and Johnson, C. R., 1985. Matrix Analysis. Cambridge University Press.
[22] Datta, B. N., 1994. ?Linear and numerical linear algebra in control theory: Some research problems?. Linear Algebra
and Its Applications, 197, pp. 755?790.
[23] Casti, J. L., 1987. Linear Dynamical Systems. Academic Press.
DS-13-1302
28
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
ned sinusoids disturbance in double mass-spring-damper system. . .
The block diagram of the mass spring damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Response of PIG and PMIG for combined step and sinusoid disturbance in double mass-spring-damper
system with imperfect model parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Bias and bias upper bound for imperfect plant model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ed
ite
d
1
2
3
4
18
19
21
23
26
27
Co
py
List of Tables
1
Parameters of the double mass-spring-damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
2
Parameters of the real system versus the system model for the double mass-spring-damper shown in Figure 3. 24
Introduction
Inverse problems in systems theory involve reconstructing the inputs to a system from measured responses. The adaptation of the Kalman filter to inverse problems, sometimes termed the Kalman inverse filter (KIF), was first proposed by
Bayless and Brigham [1] as a framework for removing sensor reverberation from seismic signals. The idea is to augment
the states of the plant model by additional dynamics that captures known structure in the inputs to be estimated. This is
equivalent to assuming the autocorrelation function of the unknown inputs or distrubances that are to be estimated.
The formulation of the KIF in [1] was for continuous-time systems. A discrete-time formulation followed shortly after
in [2]. Most of the applications of the KIF over the last thirty years have been for geophysical signal processing, see for
example [3?5], but the method has also attracted interest for input and disturbance identification in mechanical systems.
Applications include estimation of aerodynamic loads on aircraft [6, 7], friction in machines [8?10], and the estimation of
machine-environment interaction forces [11?14].
The KIF is closely linked to other methods for input identification, e.g. the regularization method of Tikhonov [15].
Work by De Nicolao and Ferrari-Trecate [16, 17] has shown that under appropriate assumptions, Tikhonov regularization
coincides with Bayesian estimation where the autocorrelation function of the unknown input is specified and the plant is a
linear system driven by white noise. De Nicolao and Ferrari-Trecate [16] show that the inverse problem can be solved by
application of the Kalman filter, revealing an equivalence between Tikhonov regularization and the KIF.
As a method for input identification the KIF has a number of attractions including: (i) that it is recursive and therefore
suitable for real-time implementations; (ii) the underlying Kalman filter structure accommodates model and measurement
uncertainty; and (iii) the method can be adapted to non-linear systems through the extended Kalman filter. Limitations of
the method are: (i) that the structure of the inputs/disturbances must be specified, and (ii) the bandwidth to which input/disturbances can be reconstructed depends on the dynamics of the plant, and the process and measurement uncertainties, as
specified in the Kalman filter.
The structure of the KIF is illustrated in Fig. 1, and has four components: the Plant, the Plant Input Generator (PIG),
the Plant Model and the Plant Model Input Generator (PMIG). The Plant is the system whose inputs are to be estimated.
The PIG represents the process that generates these to-be-estimated inputs/disturbances acting on plant. The KIF comprises
the Plant Model and the PMIG. The PMIG assumes the inputs have structure. The idea is to design the PMIG so that its
outputs serve as estimates of the input to the Plant. This is done through the feedback structure of the Kalman filter which
provides an implicit inverse model of the Plant.
This paper establishes the conditions needed for convergence of the PMIG output to the Plant input in steady state for
linear plants. Specifically, the paper establishes the conditions that must be satisfied by the Plant Model and the PMIG for
convergence of the plant input/disturbance estimates to their actual values. Section 2 describes the structure of the Kalman
Ac
ce
pt
ed
Ma
nu
sc
rip
tN
ot
1
DS-13-1302
2
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME real system
Plant Input Generator
(PIG)
Plant
plant inputs
(unknown)
+
estimation error
ed
ite
d
estimated plant inputs
system model/estimator
-
Plant Model Input Generator
(PMIG)
+
py
Plant Model
Co
Kalman gain
ot
Fig. 1: Kalman Inverse Filter block diagram.
The Kalman Inverse Filter
The KIF is described by the block diagram in Fig. 2 and governed by the equations
Ac
ce
pt
ed
2
Ma
nu
sc
rip
tN
Inverse Filter and details the notation used throughout the paper. Section 3 summarizes key assumptions. The main work
of the paper begins in Section 4 where we give general conditions for an asymptotically unbiased KIF by formulating the
problem in terms of the solution to a constrained linear matrix equation. It is shown that if the model of the Plant and PMIG
are exact, the estimates produced by the KIF converge to their true values. Section 5 presents a more detailed investigation
of the requirements on the PMIG and shows that for unbiased estimation it must span all possible inputs to the plant and that
the robustness of the estimator with respect to errors in parameters depends on the eigen-structure of the PIG. It is shown,
moreover, that a robust PMIG can be constructed from knowledge of just the number of outputs and the minimal polynomial
of the PIG dynamics. Inter alia, the analysis also yields a simple proof of the Internal Model Principle which was first proven
in [18]. Section 6 gives error bounds on estimates given knowledge of model errors. The results of this section show that the
estimated bias is approximately proportional to the parameter error and that relative bias is bounded.
DS-13-1302
x?1 = A1 x1 + B1 u + w1
(1a)
x?2 = A2 x2 + w2
(1b)
y = C1 x1 + v
(1c)
u = C2 x2
x?? 1 = A?1 x?1 + B?1 u? + Kb e
(1d)
x?3 = A3 x3 + Kc e
(1f)
(1e)
y? = C?1 x?1
(1g)
u? = C3 x3
(1h)
e = y ? y?
(1i)
z = u ? u?
(1j)
3
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017
ASME
where
thebyvector
x is the state of the Plant, x is the state of the Plant Input Generator (PIG), u is the unknown plant input
1
2
produced by the PIG and y is the measured output signal from the plant. The vector x?1 is the state of the plant model, u? is
the estimate of the input which is produced by the Plant Model Input Generator (PMIG) with state vector x3 generally not
the same size as x2 , and e is the innovation sequence available for feedback. The matrices A1 , A2 , B1 , C1 , C2 , A?1 , B?1 , C?1 , A3 ,
C3 have appropriate sizes satisfying the Eqns. (1). In particular the PMIG and the PIG have the same number of outputs and
the Plant and the Plant Model have the same number of outputs.
Plant
PMIG
Plant Model
Ma
nu
sc
rip
tN
ot
Co
py
ed
ite
d
PIG
ed
Fig. 2: The Kalman Inverse Filter, including the open-loop system governed by A1 and A2 , and the closed-loop feedback
system governed by A?1 and A3 .
Ac
ce
pt
We note, however, that the Plant Model and PMIG may contain additional dynamics not present in the Plant and PIG
under investigation or vice-versa. Hence the model state vectors x?1 , x3 may be of different sizes to the corresponding plant
w (t) state vectors x1 , x2 . The measurement noise v and the process noise w(t) = w1 (t) are vector white noises with covariance
2
structure
E[w(s)wT (t)] = Q?(s ? t)
E[v(s)vT (t)] = R?(s ? t)
E[v(s)wT (t)] = 0.
The system can be thought of as regulating z, the difference between the input signal and the estimated input signal. When
E[z] = 0 the estimates u? track the actual inputs to the Plant u.
DS-13-1302
4
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
ASME
The by
state
vector of the combined system encompassing the Plant and the PIG is
" #
x1
x=
.
x2
" #
x?1
x? =
,
x3
which shall form the model state. The state equations (1) combine to form the system
y = CP x
(3b)
(3c)
(3d)
Co
y? = CM x?
(3a)
py
x? = AP x + w,
x?? = AL x? + BP x + Kv,
ed
ite
d
The state vector of the model of the system, encompassing the Plant Model and the PMIG, is
(3e)
ot
z = DP x ? DM x?,
tN
where
"
#
A?1 B?1C3
AM =
0 A3
sc
rip
"
#
A1 B1C2
AP =
,
0 A2
AL = AM ? KCM ,
h
i
CP = C1 0 ,
h
i
DP = 0 C2 ,
" #
Kb
K=
,
Kc
BP = KCP
h
i
CM = C?1 0 ,
h
i
DM = 0 C3
Ma
nu
(4)
T ?1
K = PCM
R
pt
ed
The so-called Kalman gain K is given by
T ?1
R CM P = 0.
AM P + PATM + Q ? PCM
Ac
ce
where P satisfies the Algebraic Riccati Equation [19]
In order to regulate z, it is necessary to have closed loop stability, that is
?(AL ) ? C? ,
where ?(AL ) is the eigenvalue spectrum of the matrix AL and C? is the open left half of the complex plane. A necessary
condition for existence of a stabilizing gain K is that the pair (AM ,CM ) is detectable. We assume that the pair (AM ,CM ) is
observable; this guarantees existence of a stabilizing gain K such that AP and AL have no common eigenvalues.
DS-13-1302
5
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
by ASME
3 2017
Assumptions
The following assumptions are made at various places and referred to by number as required. They are gathered here to
avoid repetition. As indicated, some of the assumptions may be replaced with less restrictive alternatives.
ed
ite
d
A1 The plant system is observable, i.e. (AP ,CP ) is observable.
A2a All eigenvalues of the plant system matrix AP are unstable.
A2b The plant system matrix AP and the closed loop system matrix AL have no common eigenvalues.
A3 The system model is observable, i.e. (AM ,CM ) is observable.
A4 All eigenvalues of the closed loop matrix AL are stable.
A5a The Plant Model exactly describes the Plant, i.e. (A?1 , B?1 , C?1 ) = (A1 , B1 ,C1 ).
A5b The Plant Model has the same transfer function as the Plant.
A6 The Plant matrix A1 has no eigenvalues in common with the Plant Input Generator matrix A2 .
py
Assumption A4 is essential to the results of this paper. Assumption A3 implies the existence of stabilizing gains K such
that Assumption A4 is satisfied.
Convergence properties of the KIF
This section gives general conditions for an asymptotically unbiased KIF, namely E[z(?)] = 0. The approach taken is
to reformulate the problem in terms of the solution to a constrained linear matrix equation by decomposing the trajectories
into terms arising from the state and the state estimation error. Then using Assumption A4, the conditions for E[z(?)] = 0
are derived.
tN
General conditions for unbiased input estimation
Consider the state equations (3), then x and x? are related by
sc
rip
4.1
ot
Co
4
x? = ? + Xx
(5)
Ma
nu
where X has appropriate dimensions and ? represents the state estimation error. Differentiating both sides with respect to
time and using state equations (3) we can write
??? = x?? ? X x?
= AL x? + BP x + Kv ? X(AP x + w)
= AL (?? + Xx) + BP x + Kv ? X(AP x + w)
ed
= AL ? + (AL X + BP ? XAP )x + Kv ? Xw.
XAP ? AL X = BP
(6)
Ac
ce
pt
Under Assumption A4, AL is stable and hence AL ? converges to zero as t ? ?. Therefore, in order to satisfy limt?? ? = 0,
the term due to x must be zero, so
must hold. This relation is the well-known Sylvester matrix equation and has unique solution X under Assumption A2b that
is if AP and AL have no common eigenvalues [20?22]. Assumption A2b can be satisfied in a number of ways. One sufficient
condition is that Assumptions A2a and A4 hold.
Using the above relations, the input estimation error z can then be expressed in terms of the state and state estimation
error. Specifically
z = DP x ? DM x?
= (DP ? DM X)x ? DM ? .
DS-13-1302
6
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017
ASME
Since
E[??by
]?
0 the term to due to ? converges to 0 . Therefore in order for E[z] ? 0 , the component due to x must be zero,
leading to the following results:
Lemma 1 (Sufficient conditions for zero bias). Under Assumptions A2b and A4, the bias is asymptotically zero if
DM X = DP .
(7)
where X is the solution to the matrix equation (6).
ed
ite
d
Lemma 2 (Necessary conditions for zero bias). Under Assumptions A2b and A4, the bias is asymptotically zero if and only
if
(DP ? DM X)U = 0
(8)
where X is the solution to the matrix equation (6) and U is a matrix whose columns span the unstable subspace of AP .
Co
py
Note that when all eigenvalues of AP are stable, U = 0 and (8) holds regardless of the accuracy of the model. The KIF is
asymptotically unbiased because both the input signal and the input estimate converge to zero. Conversely, if all eigenvalues
of AP are unstable (Assumption A2a), U = I in (8) which reduces to to (7).
Exact Plant Model, Exact PMIG
Before discussing the question of robustness with respect to parameters, we remark that for an exact model the bias is
zero, i.e. for
tN
ot
4.2
sc
rip
(A?1 , A3 , B?1 , C?1 ,C3 ) = (A1 , A2 , B1 ,C1 ,C2 ),
A?1 = T1 A1 T1?1 ,
A3 = T2 A2 T2?1 ,
Ma
nu
the equations (6) and (7) are satisfied by X = I.
In fact, the choice of coordinates for the model state space does not affect the bias: given invertible matrices T1 and T2 ,
the model described by
C?1 = C1 T1?1 ,
B?1 = T1 B1 ,
C3 = T1C2 T2?1 ,
"
#
T1 0
X=
.
0 T2
pt
ed
has zero bias, with equations (6) and (7) satisfied by
Ac
ce
Here the converged state estimates are related to the true states by
x?1 = T1 x1
x3 = T2 x2 .
4.3
Zero-bias equations
To show that for unbiased input estimation, the plant model must be similar to the plant, partition the matrix X into four
blocks of compatible size, with
"
#
X11 X12
X=
X21 X22
DS-13-1302
7
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
by equation
ASME (6) becomes
The2017
matrix
X11 A1 = (A?1 ? KbC?1 )X11 + B?1C3 X21 + KbC1
(9a)
X21 A1 = ?KcC?1 X11 + A3 X21 + KcC1
(9b)
X11 B1C2 + X12 A2 = (A?1 ? KbC?1 )X12 + B?1C3 X22
(9c)
X21 B1C2 + X22 A2 = ?KcC?1 X12 + A3 X22 ,
(9d)
ed
ite
d
and the zero bias equation (7) becomes
C3 X21 = 0
C3 X22 = C2 .
(10a)
(10b)
Co
py
Any given Plant and PIG has a very large class of models for which (9) and (10) can be solved for X. We limit our attention
to models for which the Plant Model has the same dimensionality as the true Plant. In this situation the Plant and the Plant
Model must necessarily obey a similarity relationship. The argument is as follows.
After substitution of (10a), (9a) can be rearranged as
ot
A?1 X11 ? X11 A1 = Kb (C?1 X11 ?C1 ).
tN
In order for this equation to hold regardless of the feedback gain Kb , the left and right sides must be zero, giving
A?1 X11 = X11 A1 ,
sc
rip
C?1 X11 = C1 .
Similarly, (9c) can be reduced to
Ma
nu
(X11 B1 ? B?1 )C2 = (A?1 ? KbC?1 )X12 ? X12 A2 ,
which in order to hold regardless of the signal generator gain C2 necessitates
ed
X11 B1 ? B?1 = 0.
pt
Since the plant (A1 ,C1 ) is observable and the model has the same number of dimensions, X11 must be square, invertible and
unique (see Sec. 5), leading to the following result:
Ac
ce
Theorem 1. Suppose the Plant (A1 , B1 ,C1 ) is observable and controllable, and that the associated zero bias Plant Model
(A?1 , B?1 , C?1 ) has the same number of state space dimensions. Then the Plant and the Plant Model must be related by a unique
similarity transform
?1
A?1 = X11 A1 X11
,
?1
C?1 = C1 X11
,
B?1 = X11 B1 .
Moreover the Plant and the Plant Model have the same transfer function
C?1 (s ? A?1 )?1 B?1 = C1 (s ? A1 )?1 B1
for all s ? C.
DS-13-1302
8
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
ASME
Thisby
result
implies a requirement for exact knowledge of the Plant transfer function: to exactly determine the plant
inputs from the plant outputs requires the exact relationship between the inputs and the outputs.
In general this requirement implies a need for a perfect Plant Model including exact parameter knowledge. An interesting
exception is where all eigen-values of A1 are zero corresponding, for example, to free rigid-body dynamics problems. In this
situation the eigen-structure is independent of the parameter values. The model used by Siegrist [13] to estimate road tyre
interaction forces for off-highway mining vehicles is in this class of problems.
Conditions for unbiased input estimation with a perfect plant model
For algebraic simplicity, we will henceforth assume a perfect Plant Model with A?1 = A1 , B?1 = B1 , C?1 = C1 so that
X11 = I. We now seek conditions for unbiased input estimation under this assumption. A sufficient condition for zero bias is
the existence of a solution to the constrained Sylvester equation
ed
ite
d
4.4
C3 X22 = C2
(11b)
py
A3 X22 = X22 A2 .
(11a)
Co
The solvability of (11) becomes a necessary condition for zero bias when the signal generator has no modes in common with
the plant, as demonstrated in the following result.
tN
ot
Theorem 2. Under Assumptions A2b, A4, A5a, suppose that the PMIG matrix A3 has no eigenvalues in common with the
Plant Model matrix A?1 and that Kb has been chosen so that (A?1 ? KbC?1 ) has no eigenvalues in common with A1 or A2 (this is
possible because (A1 ,C1 ) is observable). The KIF has zero bias if and only if there exists a solution X22 to the simultaneous
matrix equations (11).
sc
rip
Proof. To prove sufficiency, let X22 be a solution to (11). Then the equations (9), (10) are uniquely satisfied by
"
#
I 0
X=
.
0 X22
Ma
nu
To prove necessity, suppose the KIF has zero bias, i.e. (10) holds. Then (9a) reduces to
X11 A1 = (A?1 ? KbC?1 )X11 + KbC1 .
X21 A1 = A3 X21
X12 A2 = (A?1 ? KbC?1 )X12
ce
pt
ed
Since the plant model (A?1 , C?1 ) is exact and A1 has no eigenvalues in common with (A?1 ? KbC?1 ), the solution X11 = I is
unique. The next two equations (9b) and (9c) reduce to
Ac
which, by the assumption of no common eigenvalues, are uniquely satisfied by X21 = 0 and X12 = 0. The other equation (9d)
reduces to
X22 A2 = A3 X22
as required.
We now show that equations (11) imply a relationship between the observability matrices of the PMIG and the PIG.
Suppose that equations (11) hold for the matrix X22 . Multiply (11a) on the right by A2 and substitute (11b) to obtain
C3 A3 X22 = C2 A2 .
DS-13-1302
9
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 bythe
ASME
Repeating
process will yield a sequence of equations
j
j
j ? 0.
C3 A3 X22 = C2 A2 ,
Terminating the sequence at j = k ? 1,1 the system can be expressed as
Ok (A3 ,C3 )X22 = Ok (A2 ,C2 ),
(12)
ed
ite
d
where Ok (A,C) is the kth order observability matrix defined by
?
C
?
?
? CA ?
?
Ok (A,C) ? ? . ?
?.
? .. ?
CAk?1
py
?
Co
When k = n, where A is n О n, the matrix O (A,C) is the observability matrix of control theory.
The noise-free component of an input signal produced by state x2 of the signal generator satisfies
ot
u = C2 x2 ,
tN
u? = C2 x?2 = C2 A2 x2 ,
u? = C2 x?2 = C2 A22 x2 ,
sc
rip
..
.
du
dx2
= C2 k?1 = C2 Ak?1
2 x2 .
dt k?1
dt
Ma
nu
Similarly the input estimate produced by state x3 satisfies
u? = C3 x3 ,
u?? = C3 x?3 = C3 A3 x3 ,
pt
ed
u?е = C3 x?3 = C3 A23 x3 ,
..
.
d u?
dx3
= C3 k?1 = C3 Ak?1
3 x3 .
dt k?1
dt
ce
These relations can be summarized in two systems of equations
u
u?
u?
..
.
Ac
?
?
?
?
?
?
?
?
?
1 As
du
dt k?1
?
?
?
?
?
?
? = O (A2 ,C2 )x2 ,
?
?
?
?
?
?
?
?
?
?
?
u?
u??
u?е
..
.
d u?
dt k?1
?
?
?
?
?
? = O (A3 ,C3 )x3 .
?
?
?
a consequence of the Cayley-Hamilton theorem of linear algebra, the sequence can be stopped at k = n2 n3 , where A2 is n2 О n2 and A3 is n3 О n3 ;
in the next section we show that k = n3 + 1 suffices.
DS-13-1302
10
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017isby
ASMEand x = X x , the PIG produces an identical signal as shown by the following relation.
If (12)
satisfied
3
22 2
?
?
?
?
?
?
?
?
?
u
u?
u?
..
.
du
dt k?1
?
?
?
?
?
?
?
?
?
?
? = O (A3 ,C3 )X22 x2 = O (A3 ,C3 )x3 = ?
?
?
?
?
?
?
u?
u??
u?е
..
.
d u?
dt k?1
?
?
?
?
?
?.
?
?
?
tN
ot
Co
py
ed
ite
d
The interpretation we give to this result is that the signals generated by the PMIG must span all possible inputs to the
plant, so that if, for example, the plant input contains a sinusoidal component of frequency ?, the PMIG must include a mode
with this same frequency. Furthermore, the matrix X22 is a map from the states of the PIG to the states of the PMIG so that
the inputs to the model correspond to those of the plant in a way that is consistent with the dynamics of the two systems.
The above results provide a test of whether a given PMIG will produce unbiased input estimates for a given PIG.
No assumptions have been made about the size of the PMIG relative to the PIG. The PMIG may have additional signal
components that do not appear in the PIG, or multiple copies of components of the signal generator. The PMIG can therefore
contain a library of modes whose form corresponding to all expected inputs to be applied to the plant. By way of example,
if it is known that the plant can subjected to inputs/disturbances having n distinct frequencies, ?1 , . . ., ?n , applied alone or
in combination, the PMIG should be designed to have eigenvalues corresponding to those frequencies.
The next section looks in more detail at the relationship between the structure of the signal generator and the structure
of models that produce zero-bias signal estimates. The analysis shows that the robustness of the estimator with respect to
certain parameters depends largely on the eigen-structure of the PIG. This amounts to determining the general structure of
pairs of matrices (A3 ,C3 ) and (A2 ,C2 ) which satisfy Equation (11) for some X22 .
Conditions on the PMIG for zero bias estimates
In this section, we consider the abstract properties of the relationship between the PIG and the PMIG implied by the
zero bias equations (11), culminating in the result that a robust PMIG can be constructed from knowledge of just the number
of outputs and the minimal polynomial of the PIG dynamics. For ease of argument, we introduce the following notation:
sc
rip
5
Ma
nu
Definition 1. For any matrices A,C, B, D of size n О n, k О n, m О m, k О m respectively, say (A,C) models (B, D), denoted
(A,C) & (B, D), if
(13a)
AX = XB,
(13b)
ed
CX = D
pt
for some matrix X of size n О m. Say (A,C) &X1 (B, D) if (11) holds for X = X1 . In particular, the PMIG will provide zero
bias estimates for a given PIG if (A3 ,C3 ) & (A2 ,C2 ).
ce
The relation has useful ordering properties:
Ac
Reflexive: (A,C) &In (A,C), where In is the n О n identity matrix. This means any plant is a model of itself.
Transitive: (A1 ,C1 ) &X1 (A2 ,C2 ) &X2 (A3 ,C3 ) implies (A1 ,C1 ) &X1 X2 (A3 ,C3 ). This means any model of a model of the
plant is, in turn, a model of the plant.
From these properties, two equivalence relations naturally arise:
Similarity: Say (A,C) ? (B, D) if (A,C) &X (B, D) where the matrix X is square and invertible. This is similarity in the
usual linear systems sense.
Equivalence: Say (A,C) h (B, D) if (A,C) & (B, D) & (A,C). Similarity implies equivalence, but not necessarily vice-versa.
Pairs of models can be related in three fundamental ways (subject, of course, to compatibility of the size of the sub
blocks):
DS-13-1302
11
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
by ASME
R12017
(change
of basis): "
(P?1 AP,CP)
&!
(Q?1 AQ,CQ) by X"=#P?1 Q
#
h
i
A B1
I
R2 (bigger model):
, C D & (A,C) by X =
0 B2
0
"
#
!
h
i
h i
A 0
R3 (unobservable signals): (A,C) &
, C0
by X = I 0
B3 B4
ed
ite
d
Relation R1 says that the choice of coordinate system is irrelevant, relation R2 says that the PMIG can contain additional
dynamics not present in the PIG. Relation R3 says that unobservable components of the PIG can be ignored.
In fact these relations characterize &, since any matrix X that relates (A,C) and (B, D) can be reduced by row and
column operations to echelon form by P?1 XQ = 0I 00 for some invertible matrices P and Q2 . This is summarized in the
following result:
Theorem 3. (A,C) & (B, D) if and only if there are change of basis matrices P and Q that reduce (A,C) and (B, D) to upperand lower- block triangular form respectively such that
(Q?1 BQ, DQ) =
#
!
i
A12 h
, C1 C2
A21
#
!
h
i
0
, C1 0 .
B22
py
(P AP,CP) =
"
A11
0
"
A11
B21
Co
?1
Ma
nu
Lemma 3. If (A,C) &X (B, D), then
sc
rip
tN
ot
The above result characterizes the perfect plant model as containing a subsystem that is equivalent to the observable
part of the plant. In terms of structural stability, this says that the model is insensitive to perturbations that do not change the
equivalence of this subsystem that models the plant.
The aim of the remainder of this section is to show that any observable model (A,C) with suitable replication of the
eigenstructure of B is a perfect model of (B, D), i.e. contains a subsystem equivalent to (B, D). The implication of this result
is that the model is insensitive to any perturbations that do not alter the eigenstructure or compromise the observability of
the model.
In order to achieve this result, we re-examine the results relating the observability matrices of the systems. In the
previous section, we demonstrated that (11) implies (12), and hence
Ok (A,C) ? X = Ok (B, D)
for all k ? 1.
(14)
pt
ed
The reverse is not necessarily true, however note that (14) implies CAi (AX ? XB) = 0 for i ? k ? 2, leading to the
following result:
On (A,C) ? (AX ? XB) = 0.
(15)
Ac
ce
Lemma 4. If (14) holds for k = n + 1, then (AX ? XB) satisfies
Moreover if (A,C) is observable (i.e. On (A,C) has rank n), AX ? XB = 0 and hence (A,C) & (B, D).
Several results immediately follow:
Corollary 3.1.
2 Note
that some of the sub blocks may have zero size.
DS-13-1302
12
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
by ASME
If (A,C)
& (B, D) and (A,C) is observable, such an X is unique.
X
If (A,C) &X (B, D) and (A,C) is observable, any matrices E, F satisfying E ? X = F are of the form E = Q ? O (A,C), F =
Q ? O (B, D) for some matrix Q.
If (A,C) &X (B, D) and (B, D) is observable, then X has full column rank.
If (A,C) &X (B, D) and (B, D) is observable and n = m, then X is invertible (hence (A,C) ? (B, D) and (A,C) is
observable).
If (A,C) h (B, D) and both are observable, then n = m and (A,C) ? (B, D).
ed
ite
d
It is well known that any two single output systems (A, cT ), (B, dT ) are equivalent if they are observable and A and B are
similar (e.g. Casti [23]). A slightly more general result follows:
Corollary 3.2. Suppose A ? B, and cT , dT are any two row vectors such that (A, cT ) is observable (single-output case).
Then (A, cT ) &X (B, dT ) with X = (On (A, cT ))?1 On (B, dT ).
py
Proof. Since (A, cT ) is observable, On (A, cT ) is n О n and invertible. Since A and B have the same characteristic polynomial,
if (14) holds for k = n it must hold for k = n + 1 (the final equation is the same linear combination of the previous n equations)
and the result follows.
ot
Co
The following result shows that every model of an output system contains a subsystem equivalent to that system.
"
#
!
A1 , A2
Lemma 5. Suppose (A,C) &X (B, D) where (B, D) is observable. Then (A,C) ?
, [C1 ,C2 ] for some A1 , A2 , A4 ,C1 ,C2
0, A4
such that (A1 ,C1 ) ? (B, D).
tN
Proof. Let (A,C) &X (B, D). Since (B, D) is observable, X has full column rank by corollary 3.1. Then
sc
rip
" #
I
Q?1 ,
X =P
0
where P and Q are invertible matrices. Hence by Eq.(13)
Ma
nu
"
#
!
i
A1 , A2 ?1 h
P
P , C1 , C2 P?1 ,
0, A4
(A,C) =
where A1 = Q?1 BQ and C1 = DQ, and the result follows
"
#
A1 0
A1 ? A2 =
,
0 A2
ce
pt
ed
Larger systems can be built up from sums of smaller models in parallel: define the direct sum of two matrices as
and we have the following result:
Ac
Lemma 6. If (A1 ,C1 ) & (B1 , D1 ) and (A2 ,C2 ) & (B2 , D2 ), then (A1 ? A2 , [C1 ,C2 ]) & (B1 ? B2 , [D1 , D2 ]) (provided the matrices C1 and C2 have the same number of rows) and (A1 ? A2 ,C1 ?C2 ) & (B1 ? B2 , D1 ? D2 ).
A generalization of this result to k-output systems forms a basis for the internal model principle:
Lemma 7. Suppose A is in block-diagonal form with k identical blocks, i.e.
?
?
A
0
1
k
M
? .
?
.. ?,
A=
A1 = ?
?
?
0
A1
DS-13-1302
13
(16)
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017
where
A by
is ASME
n О n , and
1
1
1
?
?
cT1
0
?
? .
.. ?,
C=?
?
?
T
0
ck
?
dT11
? .
.
D=?
? .
dTk1
?
dT1k
.. ?
?
. ?,
и и и dTkk
иии
..
.
ed
ite
d
where cT1 , . . . , cTk , dT11 , . . . , dTkk are row vectors. Then (A,C) & (A, D) if and only if (A1 , cTi ) & (A1 , dTij ) for all i, j, 1 ? i ?
k, 1 ? j ? k.
Proof. Setting
Co
where the sub-blocks Xi j are each of size n1 О n1 , (13) reduces to
tN
ot
A1 Xi j = Xi j A1
cTi Xi j = dTij ,
py
?
?
X11 и и и X1k
? . .
?
. . . ... ? ,
X =?
? .
?
Xk1 и и и Xkk
and the result follows.
sc
rip
We note in Lemma 7 that observability of all subsystems (A1 , cTi ) is equivalent to observability of (A,C). Also by
Corollary 3.1, if (A, D) is observable then the systems are equivalent, so we have the following result:
Ma
nu
Lemma 8. Suppose (A,C1 ), (A,C2 ) are systems where A is of block-diagonal form (16), and C1 ,C2 are k О n matrices. If
(A,C1 ) is observable, (A,C1 ) & (A,C2 ). Moreover, if both systems are observable, they are equivalent.
Proof. All that remains is to show that if (A,C1 ) is observable, there exists a vector cT1 such that (A1 , cT1 ) is observable. Since
A consists of k copies of A1 , they have the same minimal polynomial and hence
ed
rank(On (A,C)) ? km1 ,
Ac
ce
pt
where m1 is the degree of the minimal polynomial of A1 . The matrix has full rank only if m1 = n1 , i.e. A1 is cyclic (also
termed nonderogatory [21]: each eigenvalue has geometric multiplicity 1). For some change of basis matrix P, such a matrix
will have Jordan form
?
J1
?
?1
P A1 P = J = ?
?0
0
?
0 0
.. ?
. 0?
?,
0 Jr
where the Jordan blocks
?
?
?i 1 0
? .
?
?
.
Ji = ?
?0 . 0?
0 0 ?i
DS-13-1302
14
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 by to
ASME
correspond
distinct eigenvalues ? . Set dT = [1, 0, . . . , 0]; it easily follows that (J , dT ) is observable, as is (J, dT ) where
i
i
i
i
dT = [dT1 , . . . , dTk ]. Then set cT1 = dT P?1 , so we have (A1 , c1 ) ? (J, dT ) and hence (A1 , c1 ) is observable.
With this choice of c1 , apply Lemma 7 with cTi = cT1 , 2 ? i ? k to obtain C such that (A,C) & (A,C1 ) and (A,C) & (A,C2 ).
If (A,C1 ) is observable, (A,C1 ) ? (A,C) and hence (A,C1 ) & (A,C2 ). If both (A,C1 ) and (A,C2 ) are observable, they are
both equivalent to (A,C), and hence (A,C1 ) ? (A,C2 ).
This result will be used to show that any system (B, D) can be perfectly modelled by an observable system (A,C), where
A is of the form (16). The matrix A1 is any matrix consisting of the largest cyclic component of B, i.e. a cyclic matrix
satisfying the minimal polynomial of B.
Lk
A1 where A1 is the largest cyclic component of B. Then (A,C) & (B, D)
ed
ite
d
Lemma 9. Suppose (B, D) is observable, and A =
if (A,C) is observable.
Lp
"
#
Bi Bi1
B1 ?
0 Bi2
ot
for some matrices Bi1 , Bi2 , and hence we can augment the model to obtain
Co
py
Proof. Without loss of generality, assume A1 is in Jordan form. By the cyclic decomposition theorem [20, 21], B ? i=1 Bi
where the Bi are cyclic matrices in Jordan form such that B1 = A1 and each Bi satisfies the minimal polynomial of B1 .
Each Bi can be augmented to a matrix similar to B1 , i.e.
tN
"
#
!
i
B B01 h
, D 0 & (B, D),
0 B02
sc
rip
such that
"
#
p
M
B B01
B1
?
0 B02
i=1
Ma
nu
and
B01 =
ed
B02 =
p
M
i=1
p
M
Bi1
Bi2 .
i=1
ce
pt
Since (B, D) is observable, p ? k,3 and so the model can be augmented with k ? p copies of B1 to obtain
Ac
(A,C) &
"
#
!
!
k?p
h
i
M
B B01
?
B1 , D 0 & (B, D),
0 B02
i=1
and the result follows.
As an easy corollary of the previous result and relations R1 and R2, the main result follows, namely that any observable
system that models k copies of the largest cyclic component of B, is a perfect model of (B, D), where k is the number of rows
of C and D. This gives an alternative (and arguably simpler) proof of the Internal Model Principle [18].
3 PROOF:
Suppose p > k. Some eigenvalue ? has geometric multiplicity p, so rank
B??I D
? n + k ? p < n, and hence (B, D) is not observable. Hence
p ? k.
DS-13-1302
15
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 by4 ASME
Theorem
(Internal Model Principle). Suppose D has k rows, B is the largest cyclic component of B, and
1
"L
#
B1 A1
,
0 A2
k
A?
(17)
for some matrices A1 and A2 . Then (A,C) models (B, D) for any matrix C such that (A,C) is observable.
"L
ed
ite
d
Proof. Let X be a change of basis matrix such that X ?1 AX = A0 , where A0 is the right hand side of (17), and so
#
!
i
B1 A1 h
, C1 C2 ,
0 A2
k
(A,C) ?
where [C1 ,C2 ] = CX. By relations R1 and R2, (A,C) & ( k B1 ,C1 ). Since (A,C) is observable, the subsystem ( k B1 ,C1 )
is observable, and by the previous lemma, models (B, D) for any D of compatible size. Hence, by the transitive property,
(A,C) & (B, D).
L
Co
py
L
sc
rip
tN
ot
We have proved that a robust model of an observable plant (and in particular, a robust PMIG for a given PIG) can
be constructed from knowledge of only the minimal polynomial (i.e. the largest cyclic component) of the plant dynamics,
and the number of outputs, by incorporating multiple copies of the largest cyclic component. This model is robust to any
perturbations that do not affect the eigenstructure of the model, or the observability of the model, and moreover it may
contain additional dynamics not present in the underlying plant.
We now prove the converse of the internal model principle, namely that any ?universal? model of an observable plant
necessarily contains a subsystem with multiple copies of the largest cyclic component.
Theorem 5. Suppose (A,C) & (B, D) for all matrices D having the same number of rows as C, and let B1 be the largest
cyclic component of B. Then A is equivalent to a matrix of the form (17).
ed
Ma
nu
"
#
B1 , B2
Proof. Let B1 be the largest cyclic component of B, and let B ?
. As in the proof of Lemma 8, choose a row vector
0, B4
T
dT1 such that (B
Since (A,C) & (B, D) for all D, by the transitivity property and property R2, it follows
?
?1 , d1?) is ?observable.
0
? ? ??
that (A,C) & ?B1 , ?dT1 ??. By Lemma 6, it follows that (A,C) & (?k B1 , ?k dT1 ), and by Lemma 5 it follows that (A,C)
0
contains a subsystem equivalent to (?k B1 , ?k dT1 ).
Structural stability
Any PMIG satisfying (11) (for some X22 ), that is subject to perturbations in some or all of its parameters will continue
to produce unbiased input estimates provided (11) is still solvable. Such a PMIG is said to be structurally stable with respect
to the perturbed parameters. Structurally stable parameters do not need to be known exactly for unbiased input estimates.
The degree to which the solvability of (11) is stable with respect to the parameters of PIG and PMIG depends partly
on the knowledge of the problem and, for multiple input systems, on the number of identical copies of the PIG within the
PMIG.
By way of example, consider a plant having two inputs. If the inputs are sinusoids of the same frequency but different
amplitude and phase, then for unbiased input estimates the PMIG must have two undamped modes at this frequency. This
PMIG does not, however, require a priori specification of the amplitude and phase of each mode.
Ac
5.1
ce
pt
The amount of redundancy required in the model depends entirely on the knowledge of the plant; if the parameters of
the plant dynamics and outputs are known exactly, a single copy of the exact model will suffice.
DS-13-1302
16
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
by ASME
For more
complicated multiple-output signal generators, the canonical form of the system needs to be taken into consideration. As an example, the system
?
?
#
0?0 "
? abc ?
??
?
?? 0 0 0 ? ,
de f
000
??
is equivalent to the system
?
?
#
010 "
? y01 ?
??
?
?? 0 0 0 ? ,
100
000
ed
ite
d
??
1
0
0
0
0
0
0
0
?
?
0 "
#
?
0?
? 1000 ?
?,
?,
1? 0 0 1 0 ?
0
ot
0
?? 0
??
??
?? 0
0
Co
??
py
for some y depending on the parameters ?, a, b, c, d, e, f , and hence a model of this form is not robust with respect to all of
these parameters. The model
? ? ?
0
x?1
?x? ? ? 0
? 2? ?
? ? = ? ?k1 ?k2
?x?1 ? ? m
0
0
?? ? ?
?
0 0 " #
x1
? ?x ? ? 0 0 ? F
? ? 2? ?
? 1
?
c2 ? ? ? + ? 1
?
?
?
?
? F2
0
x?
1
m1
m1
?c2 ?c3
1
x?2
0 m2
m2
? ?
x1
?x ?
? 2?
y = ? ?.
?x?1 ?
x?2
1
0
0
1
Ma
nu
k2
?c2
m1
m1
?k2 ?k3 c2
m2
m2
(18)
ed
x?2
1
k2
m2
sc
rip
tN
which is observable and contains two copies of the maximal cyclic component, is structurally stable with respect to the
parameters ?, a, b, c, d, e, f .
The examples that follow illustrate input estimation and structural stability when the Plant Model and PMIG are exact.
These examples use the forced mass-spring-damper system shown in Figure 3. The system parameters are shown in Table 1.
The dynamics of the system are described by
Ac
ce
pt
Further we assume that the measurement noise and process noise are both zero-mean Gaussians with covariance R and Q
respectively.
After substituing the parameters in Table 1 the plant is represented by A1 , B1 ,C1 matrices shown below.
?
?
0
0
1 0
? 0
0
0 1?
?
?
A1 = ?
?
??10 5 ?1.5 1.5?
5 ?10 1.5 ?3
?
?
0 0
?0 0?
?
?
B1 = ?
?
?10 0 ?
0 10
C1 = I4О4
DS-13-1302
17
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Value
k1
0.5 N/m
k2
0.5 N/m
k3
0.5 N/m
m1
0.1 kg
m2
0.1 kg
c2
0.15 N.s/m
c3
0.15 N.s/m
py
Parameter
ed
ite
d
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
sc
rip
tN
ot
Co
Table 1: Parameters of the double mass-spring-damper system.
Ma
nu
Fig. 3: The double mass-spring-damper system.
h i
A2 = 0
" #
1
C2 =
0.5
ce
pt
ed
Example 1. Step Input Consider the double mass-spring-damper system. Suppose the driving forces are unknown steps,
F2
but the ratio between them is known. Suppose
= 0.5, then PIG is represented using A2 and C2 matrices.
F1
Ac
Note that since x?2 = A2 x2 + w2 and A2 = [0], x2 (0) must be set to a non-zero value in order to generate a step input. In this
example x2 (0) = 1. Assuming exact Plant Model specifying the PMIG so that it has a single eigenvalue at the zero consistent
with the structure of a step input gives.
A?1 = A1
B?1 = B1
C?1 = C1
A3 = A2
C3 = C2 .
In PMIG the initial state is set to zero i.e. x3 = 0 so that the model is unaware of the actual input generated by PIG.
Additionally we set the Q = I5О5 О 10?3 and R = 5 О I4О4 О 10?3 .
Simulating the system response produces the results shown in Figure 4, with the key point to note being that the inputs
converge to their correct values when the system transients have decayed.
DS-13-1302
18
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
3
4
2.5
3
4
3
2
2
1.5
2
1
1
0
0.5
1
-1
0
0
3
6
9
12
15
-0.5
0
3
6
9
12
15
time (s)
time (s)
1
3
6
9
12
0
15
0
3
6
time (s)
0
0
12
15
0
3
6
Ma
nu
9
time (s)
15
0
0
3
6
9
12
time (s)
ot
0.5
6
12
tN
0.5
sc
rip
1
3
9
time (s)
1
0
12
9
12
15
time (s)
ed
Fig. 4: Response of PIG and PMIG for step disturbance in double mass-spring-damper system.
Ac
ce
pt
Example 2. Compound Sinusoid-step Input Consider the double mass-spring-damper system. Suppose the driving forces
in Figure 3 are a combination of two sinusoids function with known frequencies ?1 , ?2 and a step function. The magnitude
and phase of these sinusoids and the magnitude of the the step function are unknown to the model, but the ratio between
components in F1 and F2 are known. More specifically if the inputs are written as
F1 = a1 u(t) + b1 sin(?1t + ?1 ) + c1 sin(?2t + ?2 )
F2 = a2 u(t) + b2 sin(?1t + ?1 ) + c2 sin(?2t + ?2 )
then
a2 b2 c2
, , are known.
a1 b1 c1
DS-13-1302
15
0.5
Co
0
9
py
0.5
-1
-2
6
1
1
0
3
1.5
3
2
0
time (s)
1.5
4
-2
ed
ite
d
0
19
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
15
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017
by example
ASME suppose that PIG is represented by A and C matrices shown below.
In this
2
2
ed
ite
d
?
?
0 0 0 0 0
?
?
?0 0 ?1 0 0 ?
?
?
?
A2 = ?
?0 ??1 0 0 0 ?
?
?
?0 0 0 0 ?2 ?
0 0 0 ??2 0
"
#
1 0.5 0.5 1 ?1
C2 =
0.5 0.5 ?0.5 1 1
A?1 = A1
B?1 = B1
py
h
iT
The initial state x2 (0) determines the phase and magnitude of the sinusoids in each input. We set x2 (0) = 0.5 0.25 0.25 0.5 0.5 .
Additionally we choose f1 = 1Hz, f2 = 0.2Hz and set ?1 = 2? f1 , ?2 = 2? f2 .
Assuming exact Plant Model and PMIG gives
C?1 = C1
A3 = A2
C3 = C2 .
tN
ot
Co
In the PMIG the initial state is set to zero i.e. x3 = 0 so that the model is unaware of the actual input generated by PIG.
Additionally we set the Q and R to the values of the previous example.
The the system shown in Figure 2 is simulated and the responses are presented in Figure 5. The estimated inputs are
seen to track the inputs after initial transients have decayed.
Imperfect plant model: error bounds
The results of the previous section showed that the asymptotic un-biasedness of the filter is not generally preserved
under perturbations in parameters relating to the eigenvalues of the plant and signal.
In practical applications, the parameters are never known exactly and the model is likely to capture only an approximation of the dynamics of interest, hence model error must be taken into consideration. This section shows that (i) the estimated
bias is approximately proportional to the parameter error and (ii) the relative bias is bounded.
The following result gives an estimate of the bias under the assumption that the model is ?close? to an exact model.
Ma
nu
sc
rip
6
0 , D0 ) is in error, such that
Theorem 6. Suppose that the system model (A0L ,CM
M
A0L = AL + ?AL ,
D0M = DM + ?DM ,
pt
ed
0
CM
= CM + ?CM ,
Ac
ce
where (AL ,CM , DM ) is a model producing unbiased estimates of the system (AP , BP ,CP , DP ) governed by (3). Then the bias
E[z0 ] of the imperfect system model satisfies
E[z0 ] ? ?(?DM X + (DM + ?DM )?X)E[x]
(19)
where X satisfies (6) and ?X is the unique solution to
?XAP ? AL ?X = ?AL X.
(20)
Proof. The imperfect estimate satisfies
x?? 0 = (AL + ?AL )x?0 + BP x
DS-13-1302
20
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
4
4
3
3
2
2
2
1
1
0
0
0
-1
-1
-2
-2
4
-2
0
3
6
9
12
15
-3
0
3
6
9
12
15
time (s)
time (s)
2
12
-2
-2
15
0
3
6
1
1
0.5
0.5
0
0
6
9
12
time (s)
-0.5
15
0
3
Ma
nu
3
15
0
3
6
9
12
15
time (s)
tN
1.5
12
sc
rip
1.5
0
9
time (s)
time (s)
-0.5
15
py
Co
9
ot
6
12
-1
-1
-4
3
9
0
0
0
6
1
-2
-6
3
2
1
0
0
time (s)
4
2
-6
ed
ite
d
-3
-4
6
9
12
15
time (s)
Fig. 5: Response of PIG and PMIG for combined sinusoids disturbance in double mass-spring-damper system.
z0 = DP x ? (DM + ?DM )x?0 .
pt
ed
and hence z0 satisfies
Ac
ce
For the perfect model, apply Lemma 1 to obtain
DM X = DP ,
where X satisfies
XAP ? AL X = BP .
Let ?X be the solution of
(X + ?X)AP ? (AL + ?AL )(X + ?X) = BP ,
DS-13-1302
21
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017
by ASME
which
simplifies
to (20). Defining ? 0 = x?0 ? (X + ?X)x, we have
???0 = (AL + ?AL )??0 ? (X + ?X)w.
The resulting bias E[z0 ] then satisfies
E[z0 ] = (DP ? (DM + ?DM )(X + ?X))E[x] ? (DM + ?DM )E[??0 ]
ed
ite
d
= ?(?DM X + (DM + ?DM )?X)E[x] ? (DM + ?DM )E[??0 ].
Since ? 0 is stable, E[??0 ] = 0 and the result follows.
Co
py
When the imperfect model contains multiple copies of an approximation of the maximal cyclic component of the system,
the model error can be taken to be entirely due to the error in the eigenstructure, i.e. ?CM = 0, ?DM = 0, ?AL = ?AM , where
?AM is the error of the model relative to a perfect model which contains suitably replicated copies of the true maximal cyclic
component of the system.
The above result, although it requires the solution of the matrix equation (20), does give a measure of the error of the
input estimate, modulo noise, relative to the plant state:
Corollary 6.1. The bias relative to the plant state satisfies
tN
ot
kE[u ? u?]k kE[z]k
=
? k?DM k и kXk + (kDM k + k?DM k) k?Xk .
kE[x]k
kE[x]k
sc
rip
Unfortunately the result does not give an error estimate relative to the input signal; for example if the input signal is
sinusoidal then the relative error can be infinite at the instants when the input is zero.
The next result provides an estimate of the bias when the parameter errors are large.
Ma
nu
Lemma 10. For the plant and model governed by (3) with closed loop stability, the bias E[z] is O(t n eat ) for some a ? 0 and
n ? 0. Moreover, if the plant state is bounded, with kE[x]k ? M, then E[z] is bounded with
kE[z]k ? M kDP k +C1 M kDM k max(1, kBP k /х)
for some positive constants C1 , х depending only on AL .
ed
Proof. The solution to (3) is
pt
x(t) = ?P (t)x(0) +
0 ?P (t ? ?)w(?)d?,
(21)
= exp(tAP ),
x?(t) = ?L (t)x?(0) +
?L (t)
(22)
Rt
0 ?L (t ? ?)BP x(?)d?,
(23)
= exp(tAL ).
(24)
Ac
ce
?P (t)
Rt
The matrix AL is stable, and thus
k?L (t)k ? C1 e?хt
for some positive constants C1 , х. This gives us the bound
kE[x?(t)]k ? C1 e?хt kx(0)k +C1 kBP k
DS-13-1302
22
Z t
e?х(t??) kE[x(?)]k d?.
0
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 by
ASME
It easily
follows
that if kE[x(t)]k is O(t n eat ) as t ?? ? for some a > 0 then kE[x?(t)]k is also O(t n eat ). The bias E[z] is a
linear combination of E[x?] and E[x] and thus of the same order.
For the bounded case kE[x(t)]k ? M, we have
kE[x?(t)]k ? C1 Me?хt +C1 M kBP k (1 ? e?хt )/х
= C1 M kBP k /х + e?хt C1 M(1 ? kBP k /х)
? max(C1 M,C1 M kBP k /х),
ed
ite
d
and applying the triangle inequality to
z = DP x + DM x?,
py
the result follows.
sc
rip
tN
ot
Co
The sensitivities can be difficult to calculate symbolically, so we give a numerical example here.
Ma
nu
Fig. 6: The block diagram of the mass spring damper system.
Example 3. Consider the mass spring damper system in Figure 6 with step-plus-sine and parameters w = 1, m = 1, k = 1, c =
10. Suppose that model parameters are in error, with
k? = k + ?k
c? = c + ?c,
ed
w? = w + ?w
m? = m
ce
pt
where ?w, ?k, ?c are small. For the stabilizing gains
"
Ac
Kb =
? ?
0
? ?
Kc = ?25? ,
50
#
0
100
the tracking error satisfies
ku ? u?k ? (8.77 |?w| + 0.21 |?c| + 1.41 |?k|) kxk .
Example 4. Step-Sinusoid Input with Imperfect Model Consider the double mass-spring-damper system shown in Figure 3. Suppose that the Plant Model and PMIG contain errors as shown in Table 2.
DS-13-1302
23
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
Parameter
Plant
Parameter Plant Model
0.5 N/m
k?1
0.55N/m
k2
0.5 N/m
k?2
0.55N/m
k3
0.5 N/m
k?3
0.55N/m
m1
0.1 kg
m?1
0.11kg
m2
0.1 kg
m?2
0.11kg
c2
0.15 N.s/m
c?2
0.165 N.s/m
c3
0.15 N.s/m
c?3
0.165 N.s/m
Parameter
PIG
Parameter
PMIG
f
1Hz
f?
1.1Hz
py
ed
ite
d
k1
Co
Table 2: Parameters of the real system versus the system model for the double mass-spring-damper shown in Figure 3.
B?1 6= B1
A3 6= A2
tN
A?1 6= A1
ot
Therefore
C?1 = C1 .
sc
rip
The initial state of the PIG, x2 (0), determines the real input. Specifically x2 (0) determines the phase of sinusoid component,
h
iT
the magnitude of the sinusoid component, and the size of step input component. In this example x2 (0) = 0.25 0.25 1 .
By substituting the model parameters into Equation 18 the following matrices are calculated for Plant Model.
?
0
0
1
0
?
? 0
0
0
1
?
?
A?1 = ?
?
??10 5 ?1.5273 1.5273 ?
5 ?10 1.5273 ?3.0545
?
?
0
0
? 0
?
0
?
?
B?1 = ?
?.
?9.0909
?
0
0 ?9.0909
pt
ed
Ma
nu
?
ce
C?1 = C1 = I4О4
Ac
Furthermore PMIG is given by
?
? ?
?
o 2? f? 0
0
6.9115 0
?
? ?
?
A3 = ??2? f? 0 0? = ??6.9115 0 0?
0
0 0
0
0 0
"
#
0.5 0.5 1
C3 = C2 =
.
0.5 ?0.5 0.5
Using the above matrices yields the results shown in Figure 7. To verify that the input estimation bias satisfies CorolDS-13-1302
24
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
lary2017
6.1, by ASME
kE[uu ? u?u]k
kE[xx]k
is plotted at every time step. The supremum of steady state estimate bias
kE[uu ? u?u]k
kuu ? u?uk
= lim
}
t?? kx?
xk
kE[xx]k
ed
ite
d
sup{
is compared with
k?DM k и kXk + (kDM k + k?DM k) k?Xk
kE[[y1 , y2 , y?1 , y?2 , x2 ]]k
ot
and calculate
T
F 1 , F?
F 2 ]T E[[[F 1 , F 2 ] ? [ F?
=
[y?1 , y?2 , y??1 , y??2 , x 3 ]
Co
kE[uu ? u?u]k
=
kE[xx]k
T
F 1 , F?
F 2 ]T E[[[F 1 , F 2 ] ? [ F?
py
. In this example we use
sc
rip
tN
T
F 1 , F?
F 2 ]T [[F 1 , F 2 ] ? [ F?
[y?1 , y?2 , y??1 , y??2 , x 3 ]
for every point obtained from the simulation. Figure 8 shows that as the response converges to steady state,
revealed. The figure shows
kE[uu?u?u]k
kE[xx]k
is
kE[uu ? u?u]k
} = 0.1408.
kE[xx]k
Ma
nu
sup{
Additionally by substituting the matrices into Equation 4 and solving Equation 20 the following can be obtained.
Ac
ce
pt
ed
?DM = 0 2О7
"
#
00001 1 1
DM =
0 0 0 0 1 ?1 1
?
?
?0.1288 0.1349 ?0.0198 0.0376 0.2144 0.0855 0.0348
? 0.0844 ?0.0888 0.0236 ?0.0383 ?0.0416 ?0.2628 ?0.0223?
?
?
?
?
? 0.2623 ?0.3413 ?0.0571 0.0207 0.1697 1.2240 0.0179 ?
?
?
?
?X = ?
??0.4445 0.5172 ?0.0104 0.0645 1.6985 0.0463 0.0513 ? .
?
?
??0.0415 0.0426 ?0.0088 0.0142 0.1291 1.2262 0.0124 ?
?
?
? 0.0384 ?0.0472 ?0.0049 ?0.0002 ?1.2068 0.0869 ?0.0011?
0.0167 ?0.0175 ?0.0017 ?0.0003 ?0.0522 0.0587 0.0863
Therefore k?DM k = 0, kDM k = 1.3229, kXk = 1, k?Xk = 2.1996
k?DM k и kXk + (kDM k + k?DM k) k?Xk = kDM k k?Xk = 2.9098
DS-13-1302
25
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)
2017 bythat
ASME
It follows
kE[uu ? u?u]k
? k?DM k и kXk + (kDM k + k?DM k) k?Xk .
kE[xx]k
3
3
2
2
1
1
4
3
ed
ite
d
2
1
0
0
0
-1
-1
-1
-2
3
6
9
12
15
0
3
6
time (s)
9
12
time (s)
1
1
0.8
0
0.6
-1
0.4
-2
0.2
-3
0
3
6
9
12
15
0
1.5
1
0.5
3
Ma
nu
0.5
0
6
9
12
15
12
15
time (s)
1
6
9
12
0.6
0.4
0.2
15
time (s)
time (s)
1
3
0.8
0
0
3
6
9
time (s)
sc
rip
0
0
1.2
ot
1.2
2
tN
3
-3
Co
4
15
py
0
0
-0.5
-0.5
0
3
6
9
12
0
3
6
9
12
15
time (s)
ed
time (s)
15
Conclusion
We have shown that for unbiased input or disturbance estimation using the Kalman Inverse Filter, an exact model of
the plant and a suitably replicated model describing the structure of the input signals is required. Moreover, we have shown
that the model of the input signal generator is structurally stable with respect to all parameters that do not change the eigen
structure of the PMIG. A new proof of the internal model principle has also been given.
These results provide insight to the design of Kalman Inverse Filters. In several of the recent papers using this method,
e.g. [6?13] the structural form of the PMIG is a chain of integrators without clear justification for these choices. The results
of this paper show that this approach amounts to approximating the input signal as a piecewise polynomial whose order
corresponds to the degree of the minimal polynomial of the PMIG matrix. A single integrator corresponds to a system with
DS-13-1302
Ac
7
ce
pt
Fig. 7: Response of PIG and PMIG for combined step and sinusoid disturbance in double mass-spring-damper system with
imperfect model parameters.
26
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c) 2017 by ASME
Bias relative to state
10
bias relative to state
upper bound
9
8
7
ed
ite
d
6
5
4
3
py
2
Co
1
0
4
6
8
time (s)
14
bias relative to state versus upper bound calculated by k?DM k и kXk + (kDM k + k?DM k) k?Xk.
tN
kE[uu?u?u]k
kE[xx]k ,
12
sc
rip
Fig. 8:
10
ot
2
Ma
nu
dynamics x? = 0 and thus would be consistent with an input having the structure of a constant (i.e. step) input. Similarly a
chain of two integrators is consistent with the input having the structure of a ramp. For example, in [12], a double integrator
is used to model the signal associated with each of the unknown forces acting on the wheels of a vehicle. Similarly a chain
of three integrators yield a piecewise parabolic fit to an input. De Nicolao [16], in his work on Tikhonov regularization has
similarly observed that chains of integrators can be associated with smoothing splines.
ed
For applications where the inputs have a higher structure that is known, e.g. harmonic inputs and/or disturbances, the
insights provided by this work informs us that the the PMIG should be modelled as a resonator having the appropriate
frequency. Where the signal is not a pure sinusoid, the PMIG might be extended to include harmonics of the base frequency.
The efficacy of this approach has been demonstrated in [14]. If the system has multiple inputs suitable replication of the
input structure is required. We note finally the possibility that where the structure of the input may be time varying, or not
precisely defined, the PMIG could include appropriately rich structure to adapt to the input signal.
8
Ac
ce
pt
From the point of view of the philosophical foundations of the design and in terms of guiding practical implementation
of input-and-disturbance estimators, we have found the principles revealed in this paper extremely powerful in applications
where we seek to estimate inputs and disturbances in systems using the KIF.
Acknowledgements
Kianoosh Soltani Naveh would like to acknowledge The University of Queensland Centennial Scholarship for the support of this work.
References
[1] Bayless, J. W., and Brigham, E. O., 1970. ?Application of the Kalman filter to continuous signal restoration?. Geophysics, 35(1), Feb, pp. 2?23.
DS-13-1302
27
McAree
Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Ac
ce
pt
ed
Ma
nu
sc
rip
tN
ot
Co
py
ed
ite
d
Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013;
Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267
Copyright (c)[2]2017
by ASME
Crump,
N. D., 1974. ?A Kalman filter approach to the deconvolution of seismic signals?. Geophysics, 39(1), Feb,
pp. 1?13.
[3] Dimri, V., 1992. Deconvolution and Inverse Theory. Elsevier Publishing Company.
[4] Leite, L. W. B., and da. Rocha, M. P. C., 2000. ?Deconvolution of non-stationary seismic process?. Rev. Bras. Geof.,
18(1), pp. 75?89.
[5] Rocha, M. P., and Leite, L. W., 2003. ?Treatme
Документ
Категория
Без категории
Просмотров
0
Размер файла
1 315 Кб
Теги
4038267
1/--страниц
Пожаловаться на содержимое документа