Minimax state estimation for linear stationary differential-algebraic equations

Preprints of the 16th IFAC Symposium on System Identification
The International Federation of Automatic Control
Brussels, Belgium. July 11-13, 2012
Minimax state estimation for linear
stationary differential-algebraic equations ?
Sergiy Zhuk ∗
∗
IBM Research, Dublin, Ireland (Tel: 35 31 826 94 51; e-mail:
[email protected]).
Abstract: This paper presents a generalization of the minimax state estimation approach for
singular linear Differential-Algebraic Equations (DAE) with uncertain but bounded input and
observation’s noise. We apply generalized Kalman Duality principle to DAE in order to represent
the minimax estimate as a solution of a dual control problem for adjoint DAE. The latter is
then solved converting the adjoint DAE into ODE by means of a projection algorithm. Finally,
we represent the minimax estimate in the form of a linear recursive filter.
Keywords: Minimax; Robust estimation; Descriptor systems; Optimization under uncertainties.
1. INTRODUCTION
This paper presents a generalization of the minimax state
estimation approach to linear Differential-Algebraic Equation (DAE) in the form
d(F x)
= Ax(t) + f (t), F x(t0 ) = x0 ,
(1)
dt
where F, A ∈ Rm×n . State estimation theory is one of the
main tools in mathematical modeling. It is used, in particular, for parameter identification in oceanography [Heitz
et al., 2010] or operational forecasts in meteorology and air
pollution [Wu et al., 2008]. Mathematically, geophysical
models are represented by systems of Partial Differential
Equations (PDE) that make state estimation for those
models very expensive from the computational point of
view. In practice, different model reduction techniques are
applied in order to get a reasonable computational time.
For instance, the classical Galerkin projection method
allows one to project the state of an infinite dimensional
state equation, described by a PDE, onto a finite dimensional manifold and focus on the dynamics of the projection coefficients only. In the closed form this dynamics
may be represented by the DAE (1): the main idea is to
model the impact of the unresolved part of the state vector
onto projection coefficients introducing the model error f
and restricting it to belong to a certain subspace [Mallet
and Zhuk, 2011]. Apart from the model reduction DAEs
are applied in robotics [Schiehlen, 2000]. Pros and cons of
using DAEs for modeling were discussed by Müller [2000].
The common way of deriving a state estimate for DAE
is to 1) convert the DAE into an Ordinary Differential
Equation (ODE) transforming the matrix pencil F − λA
to Weierstrass canonical form [Gantmacher, 1960] and
2) apply the classical Kalman Duality (KD) principle to
the resulting ODE. For the detailed description of 1),
we refer the reader to [Darouach et al., 1997], [Gerdin
? This work was carried out during the author’s tenure of an ERCIM
”Alain Bensoussan” Fellowship Programme. This Programme is supported by the Marie Curie Co-funding of Regional, National and International Programmes (COFUND) of the European Commission.
© IFAC, 2012. All rights reserved.
et al., 2007] and citations there. We note that applying the
matrix pencils theory to DAEs one restricts coefficients
F and A (det (F − λA) 6= 0 for some λ) and might
need to differentiate f . Although smoothness of f may be
appropriate for control problems, it turns out to be a very
restrictive assumption in the context of state estimation
where f is the model error which is often represented by
a measurable function with bounded L2 -norm. In what
follows we derive a minimax state estimate for DAE (1)
without these limitations.
The main contribution of this paper is an optimal recursive
state estimation algorithm for stationary DAEs in the
form (1). The algorithm is applied to DAE directly that
allows to avoid the restrictions imposed by the matrix
pencils theory. To achieve this we transform the state estimation problem for DAE (1) to the dual control problem
by means of generalized Kalman Duality (KD) principle
proposed by Zhuk [2012] and the latter is then solved
using a projection method. For the case of ellipsoidal G
(see formula (7)) generalized KD states that the optimal
Rt
estimate û(y) = t01 ûT ydt of the linear function `T F x(t1 )
may be found as a solution of quadratic control problem
for adjoint DAE:
d(F 0 z)
= −A0 z(t) + H 0 u(t), F 0 z(t1 ) = F 0 ` . (2)
dt
If F = I then this statement reduces to the classical KD
principle (see [Åström, 2006]). In the general case, the
set of all ` ∈ Rm such that (2) has a solution — this
set will be referred to as a minimax observable subspace
L (t1 ) — describes all possible directions in the state space
of (1) that have the following property: if ` ∈ L (t1 ) then
the worst-case estimation error (see formula (9)) is finite.
In other words, the estimate û(y) with finite worst-case
estimation error may be constructed for the projection
of F x(t1 ) on L (t1 ) only. It is worth noting that for all
` ∈ L (t1 ) the estimate of `F x(t1 ) corresponds, in fact, to
the solution of (1) which has the smallest eucledian norm
and is, therefore, unique (see formula (5) and the following
discussion). Therefore, there is no need to restrict F and
A in order to force unique solvability for (1). Instead, it
143
16th IFAC Symposium on System Identification
Brussels, Belgium. July 11-13, 2012
is sufficient to take ` from L (t1 ) and so constraints are
moved from F, A onto `.
We present a projection algorithm allowing one to convert (2) into equivalent ODE. It efficiently computes the
minimax observable subspace and determines algebraic
constraints for u imposed by the structure of (2). The
algorithm splits DAE (2) into differential and algebraic
parts (see (13)-(14)) and constructs an ODE (see (22))
such that all its solutions belong to the linear manifold
defined by the algebraic part of DAE. This ODE is obtained applying a finite number of linear transformations
(see (15)- (18)) to the differential part of DAE. We stress
that in the general case the number of required transformations would be infinite if the DAE was non-stationary.
An alternative way of constructing the set of all solutions
for (2) is based on the theory of singular matrix pencils. We
stress that this approach constrains u and its derivatives
as well (see [Gantmacher, 1960, p.336,form.(79)]) that, in
turn, complicates the procedure of computing the optimal
control.
Transforming the cost function of the dual control problem
we find an optimal estimate û using Pontryagin maximum
principle. We note that necessary optimality conditions for
linear quadratic control problem with stationary regular
DAE as a constraint were derived in Bender and Laub
[1987]. The general case of a singular non-stationary DAE
was considered by Zhuk [2012] where Tikhonov regularization is used to derive optimality conditions. Although the
approach of this paper is applicable to the class of DAEs
with constant coefficients only, it gives exact necessary and
sufficient optimality conditions in contrast to the general
approach presented in [Zhuk, 2012].
Finally, we derive the minimax estimate in the form of
the linear recursive filter. We refer the reader to [Milanese
and Tempo, 1985], [Chernousko, 1994] and [Kurzhanski
and Vályi, 1997] for the basic information on the minimax
framework. Minimax estimates for linear singular DAEs
with discrete time were obtained by Zhuk [2010]. Let
us note that our approach is applicable for the case of
unbounded input f : indeed, one can always introduce
an
extended state x̃ = (x, f )0 , a new matrix F1 = F 0 and
add a dummy bounded input f1 in (1). In this regard let us
mention H∞ framework [Başar and Bernhard, 1995] which
allows one to deal with unbounded f . Connections between
minimax and H∞ frameworks were revealed in [Baras and
Kurzhanski, 1995].
This paper is organized as follows. At the beginning of
section 2 we describe the formal problem statement and
definitions, and discuss ill-posedness of DAE. Then we
present the dual control problem (Proposition 2) and projection algorithm (Proposition 3). At the end of the section
the state estimation algorithm is derived (Theorem 4).
In section 3 we consider a simple ODE with unbounded
inputs and describe its minimax observable subspace and
derive the minimax estimate. Section 4 contains conclusions. All proofs and technical statements are located in
Appendix.
Notation: Rn denotes the n-dimensional Euclidean space;
L2 (t0 , t1 , Rm ) denotes a space of square-integrable functions with values in Rm (in what follows we will often
write L2 (t0 , t1 ) referring L2 (t0 , t1 , Rk ) where the dimension
k will be defined by the context); H1 (t0 , t1 , Rm ) denotes a
space of absolutely continuous functions with L2 (t0 , t1 )-
derivative; the prime 0 denotes the operation of taking
the adjoint: L0 denotes adjoint operator, F 0 denotes the
transposed matrix; R(L), N (L) denote the range and
null-space of the operator L, PR (F )(PN (F )) denotes the
orthogonal projection matrix onto the range (null-space)
of the matrix F ; c(G, ·) denotes the support function of
a set G; δ(G, ·) denotes the indicator of G: δ(G, f ) = 0 if
f ∈ G and +∞ otherwise; xT y denotes the inner product
of vectors x, y ∈ Rn , kxk2 := xT x; S > 0 means xT Sx > 0
for all x ∈ Rn ; F + denotes the pseudoinverse matrix; In
denotes n × n-identity matrix, 0n×m denotes n × m-zero
matrix, I0 := 0;
2. LINEAR MINIMAX ESTIMATION FOR DAES
Consider a pair of systems
d(F x)
= Ax(t) + f (t) ,
dt
y(t) = Hx(t) + η(t) ,
m×n
F x(t0 ) = x0 ,
(3)
(4)
p×n
where F, A ∈ R
,H∈R
, t0 , t1 ∈ R and x(t) ∈ Rn ,
f (t) ∈ Rm , y(t) ∈ Rp , η(t) ∈ Rp represent the state, input,
observation and observation’s noise respectively.
In what follows we assume that initial condition x0 , input
f and observation’s noise η are unknown and belong to
the given bounding set G . Our aim is to estimate the
state F x(t1 ) of (3) given observations y(t), t ∈ (t0 , t1 ).
Let us briefly describe the main points of our strategy. In
the case F = I the classical state estimation procedure
(see [Åström, 2006]) suggests to look for an estimate of
a linear function `(x) := `T F x(t1 ) in the class of linear
Rt
functions u(y) := t01 uT (t)y(t)dt, provided x solves (3)
and the output y is given in the form (4). One computes an
estimate û which is optimal with respect to some criteria
σ̃. Assuming that G is bounded, σ̃ may be defined (see
[Nakonechny, 1978] for details) as a worst-case estimation
error
σ̃(t1 , `, u) :=
sup (`(x) − u(y))2
(x0 ,f,η)∈G
where x is a unique solution of ODE (3) corresponding
to data x0 , f . In this case, the optimal estimate û solves
the equation σ̃(t1 , `, û) = inf u σ̃(`, u). In general case
F ∈ Rm×n (3) is solvable not for all x0 and f , and so
one needs to assume that there is at least one x0 , f, η in G
such that (3) has a solution. Further, if (3) is solvable for a
given x0 , f , it may have non-unique solution so σ̃(t1 , `, u) is
not well-defined. In fact, any solution of (3) corresponding
to the data x0 , f admits the following representation:
x = x0 + x1 , x0 ∈ X0 , x1 ⊥ X0 ,
(5)
where X0 contains all x such that x solves (3) with x0 = 0
and f = 0. We stress that the data x0 , f determines x1
uniquely and so we may define the worst-case error as
follows:
σ(t1 , `, u) := sup sup (`(x) − u(y))2 ,
x∈X0 x0 ,f,η
where an additional sup is taken over the set X0 . This
allows one to eliminate the part of the solution x belonging
to X0 . Notice that, in contrast with ODEs, σ(t1 , `, u)
may be infinite for some `, u as X0 is a linear subspace.
This observation suggests to introduce a so called minimax
observable subspace L (t1 ).
144
16th IFAC Symposium on System Identification
Brussels, Belgium. July 11-13, 2012
Definition 1. Given t1 < +∞, u ∈ L2 (t0 , t1 , Rp ) and
` ∈ Rm define a worst-case estimation error
Z t1
σ(t1 , `, u) :=
sup
(`T F x(t1 ) −
uT (t)y(t)dt)2
x∈X0 ,(x0 ,f,η)∈G
t0
Rt
A function û(y) = t01 ûT (t)y(t)dt is called a minimax
estimate in the direction ` (`-estimate) if inf u σ(t1 , `, u) =
σ(t1 , `, û). The number σ̂(t1 , `) := σ(t1 , `, û) is called a
minimax error in the direction ` at time-instant t1 (`error). The set L (t1 ) := {` ∈ Rm : σ̂(t1 , `) < +∞} is
called a minimax observable subspace.
We impose the following uncertainty description. Let
Q0 , Q(t) ∈ Rm×m , Q0 = Q00 > 0, Q = Q0 > 0, R(t) ∈
Rp×p , R0 = R > 0; Q(t), R(t), R−1 (t) and Q−1 (t) are
continuous functions of t on [t0 , t1 ]. Define
Z t1
E(x0 , f, η) := xT0 Q0 x0 +
f T Q(t)f + η T R(t)ηdt (6)
t0
and let
G = {(x0 , f, η) : E(x0 , f, η) ≤ 1}
(7)
It is easy to see that there is at least one (x0 , f, η) ∈ G
such that (3) has a solution x.
A geometrical representation of the `-estimate and error. Define a linear function Fu (x, η) := `T F x(t1 ) −
R t1 T
u (t)y(t)dt and operator
t0
d(F x)
(Lx)(t) = F x(t0 ),
− Ax(t) , x ∈ D(L)
dt
(8)
D(L) := {x ∈ L2 (t0 , t1 ) : F x ∈ H1 (t0 , t1 , Rm )}
and note that (3) is equivalent to Lx(t) = (x0 , f (t)) and
X0 = N (L). Since (x0 , f, η) ∈ G we can write that:
1
2
σ (t1 , `, u) =
sup {Fu (x, η) : (Lx, η) ∈ G }
Proposition 3. Take A01,2 ∈ Rn×n , A03,4 ∈ Rm×n
consider DAE
dp
= A01 p + A02 q, p(t1 ) = w,
dt
0 = A03 p + A04 q .
Let us set by definition
Ak+1
= (Ak1 − Ak2 (Ak4 )+ Ak3 )PN (I − PR (Ak4 ))Ak3 ,
1
Ak+1
= Ak2 PN (Ak4 ),
2
Ak+1
= (I − PN (I − PR (Ak4 ))Ak3 )Ak+1
,
3
1k+1
k+1
k
k
A4 = (I − PN (I − PR (A4 ))A3 )A2 ,
Pis := usk=i PN (I − PR (Ak4 ))Ak3 ,
pk = PN (I − PR (Ak4 ))Ak3 pk+1 , pj = p∗ ,
q k = −(Ak4 )+ Ak3 pk + PN (Ak4 )q k+1 , q j = q ∗ .
Then 1) there exists s ≤ n such that As3 = As4 =
DAE (13)-(14) has a solution p, q iff w = P0s w;
0
0
0 0
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
0m×n ; 2)
3) if p, q
solve (13)-(14) then p = p , q = q where p , q are defined
from (20)-(21) with p∗ := p, q ∗ := q, k ∈ 0, s − 1 and
dp
= As1 p(t) + As2 q(t), p(t1 ) = w, q ∈ L2 (t0 , t1 ). (22)
dt
Im
and set:
Let us define I˜ = 0p×m
0
+ 0 0
0
A1 = −(F ) A PR (F ) , A2 = −(F + )0 A0 PN (F 0 ) (F + )0 H 0 ,
A03 = −PN (F )A0 PR (F ), A04 = −PN (F )A0 PN (F 0 ) PN (F )H 0
s
0
−1
i
0
Q̃ = Q
:= PN 0(F ) I0p × us−1
−1 , β
i=0 PN (A4 ) ,
0
R
γ := β s 0 Q̃β s , αs := −
s−1
X
−1
i + i s
PN (Ai−1
4 )(A4 ) A3 Pi , A4 = 0,
i=0
−1
0 +
0
D := (I − (PN (F 0 )Q−1
0 PN (F )) PN (F )Q0 ),
˜ 0s + P0s 0 I˜0 Q̃αs ,
˜ 0s + αs 0 Q̃IP
Q := αs 0 Q̃αs + P0s 0 I˜0 Q̃IP
(9)
x∈D(L),η
so that `-error is a support function of the convex set
X = {(x, η) : (Lx, η) ∈ G }
(10)
The support function of X describes the distance of the
supporting hyperplane H (u) := {(x, η) : Fu (x, η) ≤
1
σ 2 (t1 , `, u)} from the origin. Then, the `-estimate û defines
a direction in a functional space L2 (t0 , t1 ) such that the
distance of H (û) from the origin is minimal and so in this
direction the maximal deviation of elements of X from
the origin is minimal.
Proposition 2. (dual control problem). Take ` ∈ Rm . The
`-error σ̂(t1 , `) is finite iff there are z ∈ L2 (0, t1 , Rm ) and
u ∈ L2 (0, t1 , Rp ) such that:
d(F 0 z)
= −A0 z(t) + H 0 u(t), F 0 z(t1 ) = F 0 ` .
(11)
dt
If σ̂(t1 , `) < +∞ then `-estimate û is a unique solution of
the following optimal control problem:
Z t1
I (z, d, u) :=
(uT R−1 u + z T Q−1 z)dt(:= I1 (z, u))
and
˜ s , Q̃0 := D0 Q−1 D , Q0 := P s 0 Q̃0 P s .
B s := αs + IP
0
0
0
0
The next theorem describes the minimax observable subspace, `-error and represents `-estimate in the form of a
linear recursive filter on a subspace.
Theorem 4. Assume that s ≤ n is defined as in Proposition 3 and define
dK
= (−As1 0 + B s 0 Q̃β s γ + As2 0 )K
dt
+ K(−As1 + As2 γ + β s 0 Q̃B s ) − KAs2 γ + As2 0 K
+
+ Q − B s 0 Q̃β s γ + β s 0 Q̃B s , K(t0 ) = Q0 ,
dx̂
= (−As1 0 + B s 0 Q̃β s γ + As2 0 − KAs2 γ + As2 0 )x̂
dt
+ (αs 0 + (KAs2 − B s 0 Q̃β s )γ + β s 0 ) y0 , x̂(t0 ) = 0 .
Then û(y) = `T F F + x̂(t1 ) and
σ̂(t1 , `) = `T F F + K(t1 )F F + ` , ∀` ∈ L (t1 ) ,
L (t1 ) = {` ∈ Rn : P0s (F + )0 F 0 ` = (F + )0 F 0 `} .
(23)
t0
+ (F
0+
3. EXAMPLE
+
0
F z(t0 ) − d)T Q−1
F 0 z(t0 ) − d) → min
0 (F
0
d(F 0 z)
= −A0 z + H 0 u,
dt
z,d,u
F 0 d = 0, F 0 z(t1 ) = F 0 ` .
(12)
Next proposition presents the projection algorithm allowing to convert DAE to ODE.
This example illustrates application of Theorem 4 to an
ill-posed DAE which represents an ODE with unbounded
inputs. As it was pointed out in Zhuk [2012] application of
Pontryagin maximum principle to the dual control problem for this DAE would lead to the restriction of the minimax observable subspace. We apply a projection algorithm
145
16th IFAC Symposium on System Identification
Brussels, Belgium. July 11-13, 2012
(Proposition 3) in order to preserve the structure
h 1of0L0(t01i).
0 0 1 0 1
0
0
0
Let F = [ 0 1 0 0 ], A = −1 0 0 −1 , H = 0 0 0 1 .
0 1 0 0
dx
):
Then (3)-(4) has the following form:(ẋ :=
dt
ẋ1 = x3 + f1 , ẋ2 = −x1 − x4 + f2 , x1,2 (t0 ) = x1,2
0 ,
y1 = x1 + η1 , y2 = x4 + η2 , y3 = x2 + η3 .
where x1,2
0 , f1,2 and η1,2,3 belong to (7) and E in (6) is
defined with Q0 = I2 , Q = I2 , R = I3 . We see that
x3,4 ∈ L2 (t0 , t1 ) are arbitrary functions. In this case, by
Proposition 2 the `-estimate û of F x(t1 ) may be obtained
minimizing the worst-case error
Z t1 X
Z t1
3
2
X
2
2
u2j dt (24)
σ(t1 , `, u) =
zi dt +
kzi (t0 )k +
i=1
t0
t0
j=1
over solutions to the adjoint DAE:
ż1 − z2 − u1 = 0, z1 (t1 ) = `1 , −z1 ≡ 0 ,
(25)
ż2 − u3 = 0, z2 (t1 ) = `2 , −z2 − u2 = 0 .
The structure of (25) is simple and so we can explicitly
compute û minimizing (25). This will allow to compare
the corresponding `-estimate with the one derived from
Theorem 4. (25) implies that û1,2 = −z2 . Hence, û3 is a
unique solution to the following control problem
Z t1
σ(t1 , `, u) = kz2 (t0 )k2 +
3z22 + u23 dt,
Fig. 1. `-estimate ûε (y) = x̂2 (t) (Dotted), error σ̂(t, `)
(DotDashed) and simulated trajectory x2 (t) (Solid),
t ∈ [0, 1]. It is easy to see that the minimax estimate is
a “central” point of the set of all possible realizations
of the “true” state x2 . The estimate is “biased” due
to the presence of unknown functions x2 , x3 . Nevertheless, the minimax error stays bounded at infinity
that proves robustness of the `-estimate with respect
to unknown and possibly unbounded disturbances in
the state equation.
t0
ż2 = u3 , z2 (t1 ) = `2 .
The optimality condition takes the classical form: û3 = p
where
ż2 = p, z2 (t1 ) = `2 ,
ṗ = 3z2 , p(t0 ) = z2 (t0 ).
Introducing k as a solution of the Riccati equation k̇ = 3−
k 2 , k(0) = 1 we find that û3 = kz2 where ż2 = kz2 ,
z2 (t1 ) = `2 . Let x̂ be a solution to
x̂˙ = −kx̂ − y1 − y2 + ky3 , x̂(t0 ) = 0 .
(26)
R t1 T
Then it is easy to see that û(y) = t0 û ydt = `2 x̂(t1 ) and
so x̂ represents a minimax state estimate.
Now, let us use Proposition 3. We find that s = 2,
L (t1 ) = {0} × R and
0 0 k1 k2 k̇1 k̇2
0 0 − k1 k2
=
0
3
0 1
k
k
k2 k4 ,
2
4
k̇2 k̇4
(27)
x̂˙ 1 = k2 y3 , x̂˙ 2 = −k2 x̂1 − k4 x̂2 − y1 − y2 + k4 y3 ,
0 ) k2 (t0 )
with x̂1,2 (t0 ) = 0 and kk12 (t
= 00 01 . From the
(t0 ) k4 (t0 )
first equation of (27) we see that k1 = k2 = 0 so that
x̂1 = 0 and the equation for x̂2 coincides with (26).
In order to generate y1,2,3 we take t0 = 0, t1 = 1,
x3 = cos(t) and x4 = sin(t), x1 (0) = 0.1, x2 (0) = −0.1,
f1 = f2 = 0, η1 = −0.1, η2 = −0.2, η3 = 0.3. In
Figure 1 the `-estimate û(y), and `-error are presented.
As L(t) ≡ {0} × R we see that x1 is not observable in
the minimax sense although y1 = x1 + η1 is observed.
This can be explained as follows: the derivative x3 of x1
may be any element of L2 . It is not hard to compute
RT
that the expression for σ(t1 , `, u) contains t0 x3 z1 dt and
so σ(t1 , `, u) < +∞ if and only if z1 (t) ≡ 0 for any
t ∈ (t0 , t1 ). As z1 is absolutely continuous it follows that
z1 (t1 ) = 0. Hence, if `1 = 0 then `-error is finite. If
`1 6= 0 then the only candidate for û1 is the impulse
control δ(t1 −t)`1 switching z1 from 0 to `1 at time-instant
t1 . However, δ 6∈ L2 so that `-error is infinite. In other
words, the algebraic structure of the adjoint DAE (25)
can not “compensate” the unbounded derivative of x1
though it compensates the unbounded derivative of x2 and
so the minimax observable subspace is non-trivial in that
direction. It is curious to note that minimax error in the
direction (0, 1) stays bounded with time as the solution of
Riccati equation is bounded at infinity.
4. CONCLUSION
In this paper we discuss the minimax state estimation
approach for linear stationary DAEs with unknown but
bounded initial condition, input and observation noise.
Our approach is based on the generalized Kalman duality principle that allows to avoid constraints on DAE
coefficients and input which are usually imposed in order
to apply the theory of matrix pencils. The presented approach is relevant to high-dimensional numerical models
that require state estimation algorithms in reduced form
due to computational burden. In this case F might encapsulate the vectors which span a low dimensional reduced
sub-space and the state of the resulting DAE describes
the evolution of the projection coefficients. This reduction approach addresses stability issues of the classical
Galerkin projection method (see [Mallet and Zhuk, 2010]
for details). In perspective, we plan to apply the presented
duality concept to stochastic DAEs. Another promising
direction is an ensemble state estimation where one state
is shared by different models and each model describes
the dynamics of the same process reflecting various physical/chemical/numerical parameterizations for the same
process (see [Garaud and Mallet, 2011] for an example in
air quality modeling).
146
16th IFAC Symposium on System Identification
Brussels, Belgium. July 11-13, 2012
REFERENCES
A. Albert. Regression and the Moor-Penrose pseudoinverse. Acad. press, N.Y., 1972.
K. Åström. Introduction to Stochastic Control. Dover,
2006.
T. Başar and P. Bernhard. H∞ -optimal Control and
Related Minimax Design Problems. Springer, 1995.
J. S. Baras and A.B. Kurzhanski. Nonlinear filtering: The
set-membership and the H∞ techniques. In Proc. 3rd
IFAC Symp.Nonlinear Control Sys.Design. Pergamon,
1995.
D. Bender and A. Laub. The linear-quadratic optimal
regulator for descriptor systems. IEEE TAC, 32(8):672–
688, 1987.
F. L. Chernousko. State Estimation for Dynamic Systems
. Boca Raton, FL: CRC, 1994.
M. Darouach, M. Boutayeb, and M. Zasadzinski. Kalman
filtering for continuous descriptor systems. In ACC,
pages 2108–2112, Albuquerque, June 1997. AACC.
W. Fleming and R. Rishel. Deterministic and Stochastic
Optimal control. Springer, 1975.
F. Gantmacher. The theory of matrices. Chelsea Publish.Comp., N.-Y., 1960.
Damien Garaud and Vivien Mallet. Automatic calibration
of an ensemble for uncertainty estimation and probabilistic forecast: Application to air quality. J. Geophys.
Res., 116, 2011.
Markus Gerdin, Thomas B. Schön, Torkel Glad, Fredrik
Gustafsson, and Lennart Ljung. On parameter and
state estimation for linear differential-algebraic equations. Automatica J. IFAC, 43(3):416–425, 2007.
D. Heitz, E. Mémin, and C. Schnörr. Variational fluid
flow measurements from image sequences: synopsis and
perspectives. Exp. Fluids, 48(3):369–393, 2010.
A. Ioffe and V. Tikhomirov. Theory of extremal problems.
North-Holland, Amsterdam, 1974.
Alexander Kurzhanski and István Vályi. Ellipsoidal calculus for estimation and control. Systems & Control:
Foundations & Applications. Birkhäuser Boston Inc.,
Boston, MA, 1997. ISBN 0-8176-3699-4.
V. Mallet and S. Zhuk.
Reduced minimax
state estimation.
Technical Report RR7500,
INRIA,
Paris-Rocquencourt,
2010.
http://hal.archives-ouvertes.fr.
V. Mallet and S. Zhuk. Reduced minimax filtering
by means of differential-algebraic equations. In Proc.
of 5th Int. Conf. on Physics and Control (PhysCon
2011). IPACS Electronic Library, 2011. Available at:
lib.physcon.ru.
Mario Milanese and Roberto Tempo. Optimal algorithms
theory for robust estimation and prediction. IEEE
Trans. Automat. Control, 30(8):730–738, 1985.
P. Müller. Descriptor systems: pros and cons of system
modelling by differential-algebraic equations. Mathematics and Computers in Simulation, 53(4-6):273–279,
2000.
A. Nakonechny. A minimax estimate for functionals of the
solutions of operator equations. Arch. Math. (Brno), 14
(1):55–59, 1978. ISSN 0044-8753.
W. Schiehlen. Force coupling versus differential algebraic
description of constrained multibody systems. Multibody
System Dynamics, 4:317–340, 2000.
Lin Wu, Vivien Mallet, Marc Bocquet, and Bruno
Sportisse. A comparison study of data assimilation
algorithms for ozone forecasts. J. Geophys. Res., 113
(D20310), 2008.
S. Zhuk.
Minimax state estimation for linear
discrete-time differential-algebraic equations.
Automatica J. IFAC, 46(11):1785–1789, 2010.
doi:
10.1016/j.automatica.2010.06.040.
S. Zhuk. Kalman duality principle for a class of illposed minimax control problems with linear differentialalgebraic constraints. Applied Mathematics & Optimisation, 2012. submitted.
APPENDIX
Proof of Proposition 2. Generalized Kalman duality
for DAEs in the form (3) with zero initial condition was
presented in [Zhuk, 2012, Th.2.4]. The proof given there
is based on the operator interpretation of DAEs (see
geometrical representation of `-estimate above). Applying
the method of Zhuk [2012] (see the first part of the proof
of Theorem 2.4) to the definition of operator L given by
formula (8) it is easy to see that (11) has a solution and
(see [Zhuk, 2012, frml.(2.10)]):
1
σ 2 (t1 , `, u)
=
min
+
(z0 ,v)∈N (L0 )
c(G , F 0 F 0 z(t0 ) − z0 , z − v, u) ,
where L0 is defined as follows:
dF 0 v
− A0 v, b ∈ D(L0 ),
(L0 b)(t) = −
dt
D(L0 ) := {b = (z0 , v) : F 0 v ∈ H1 (t0 , t1 , Rm ),
F 0 v(t1 ) = 0, v0 = F
+0
(.1)
(.2)
F 0 v(t0 ) + d, F 0 d = 0} .
Since (z0 , v) ∈ N (L ) and z solves (11), it follows that z+v
solves (11), and so taking the min in (.1) with respect to
(z0 , v) ∈ N (L0 ) is equivalent to minimize with respect to
all possible d and z such that d : F 0 d = 0 and z solves (11)
for the fixed u:
1
σ 2 (t1 , `, u)
(.3)
+
=
min
c(G , F 0 F 0 z(t0 ) − d, z, u)
0
0
z solves (11) and F d=0
˜ Noting that the support
and min is attained at z̃, d.
function of the ellipsoid G coincides with I , namely
+
c2 (G , F 0 F 0 z(t0 ) − d, z, u) = I (z, d, u) and recalling (.3)
we obtain
σ(t1 , `, u) =
min
I (z, d, u)
(.4)
0
z solves (11) and F d=0
This completes the proof.
Proof of Proposition 3. Consider (14). It is solvable iff
(I − PR (A04 ))A03 p(t) ≡ 0, and its general solution has the
following form (see [Albert, 1972]): q(t) = −(A04 )+ A03 p(t)+
PN (A04 )v(t), v ∈ L2 (t0 , t1 ). Using this we transform (13)(14) into the equivalent form (k = 0):
dp
= (Ak1 − Ak2 (Ak4 )+ Ak3 )p(t) + Ak2 PN (Ak4 )v(t),
(.5)
dt
(I − PR (Ak4 ))Ak3 p(t) ≡ 0, p(t1 ) = w .
(.6)
Define r1 := dimN ((I − PR (A04 ))A03 ) and r0 := n. There
are three different possibilities: r1 = r0 , r1 = 0 and
0 < r1 < r0 . Consider the case 0 < r1 < r0 . Define
147
16th IFAC Symposium on System Identification
Brussels, Belgium. July 11-13, 2012
dpk
= Ak1 pk (t) + Ak2 q k (t),
dt
0 = Ak3 pk (t) + Ak4 q k (t) .
k
p (t1 ) = w,
(.7)
(.8)
We claim that for k = 1 any solution of (.7)-(.8) corresponds to some solution of (13)-(14) through (20)-(21),
provided P01 w = w. Let us prove this. Assume P01 w = w
and p1 , q 1 solve (.7)-(.8). Let p0 , q 0 be defined by (20)(21) with p∗ = p1 , q ∗ = q 1 , k = 0. Then, by direct
substitution one checks that p0 , q 1 solve (.5)-(.6) so p0 , q 0
solve (13)-(14). On the other hand, if p, v solve (.5)-(.6)
then P01 w = w and p, v verify (.7) (this may be checked by
direct substitution noting that PN ((I − PR (A04 ))A03 )p = p
by (.6)) and (.8) (one checks this differentiating (.6)).
Now, as above we transform (.7)-(.8) into (.5)-(.6) but with
k = 1 and define r2 := dim N ((I −PR (A14 ))A13 ). Assuming
0 < r2 < r1 we continue constructing (.7)-(.8) and transforming it into (.5)-(.6). At some point for ri := dim N (I−
i−1
PR (Ai−1
should be either ri = ri−1 or ri = 0.
4 ))A3
i−1
Assume ri = 0. Then N ((I − PR (Ai−1
4 ))A3 ) = {0}
i−1
i−1
and, thus, PN (I − PR (A4 ))A3
= 0n so by (15),(17):
Ak1 ≡ 0n×n and Ak3 ≡ 0m×n for k ≥ i. We also have Ai2 =
i−1
i+1
i
i
Ai−1
= Ai2 PN (Ai4 ) = 0n
2 PN (A4 ), A4 = A2 so that A2
i+1
and thus A4 = 0m×n by (18). We see that the statement
1) of the proposition holds for s = i + 1. If ri = ri−1
i−1
n
then N ((I − PR (Ai−1
so that PN (I −
4 ))A3 ) = R
i−1
i−1
PR (A4 ))A3
= In and thus Ai3 = 0m×n , Ai4 = 0m×n
by (17)-(18) so 1) holds with s = i.
Assume 1) is verified for k = s. By construction, for any
solution p, q of (13)-(14) we have that p = p0 and q = q 0
where p0 , q 0 are defined by (20) and (21) with p∗ = ps
and q ∗ = q s , provided ps , q s solve (.7)-(.8) and P0s w = w.
Since As3 = As4 = 0m×n it follows that (22) is equivalent
to (.7)-(.8) for k = s and so statement 3) is verified. By
3) w = p(t1 ) = p0 (t1 ) = P0s p(t1 ) = P0s w so equality
w = P0s w follows from the solvability of (13)-(14). On
the other hand, if P0s w = w and p, q solve (22) then p0 , q 0
defined by (20) and (21) with p∗ = p and q ∗ = q solve (13)(14) by construction. This completes the proof.
Proof of Theorem 4. Define p := PR (F )z, q̃ :=
PN (F 0 )z, q := (q̃, u)0 and let A01,2,3,4 be defined as
above. We note that F 0 z = F 0 p and so (11) is equivalent to DAE (13)-(14). By statement 2) of Proposition 3
DAE (13)-(14) is solvable iff P0s (F + )0 F ` = (F + )0 F `. This
and Proposition 2 prove the second line in (23).
As we saw above (11) is equivalent to DAE (13)-(14). Let
us represent I in (12) as a function of p and q. Recalling
0
that p(t0 ) = F + F 0 z(t0 ) we write:
−1
I (z, d, u) = kQ0 2 (p(t0 ) − d)k2 + I1 (p + q̃,
0 I
q) .
Define I˜1 (z, u) := mind:F 0 d=0 I (z, d, u). It is easy to com−1
pute that mind:F 0 d=0 kQ0 2 (p(t0 ) − d)k2 = pT (t0 )Q̃0 p(t0 ) .
and so
I˜1 (z, u) = pT (t0 )Q̃0 p(t0 ) + I1 (p + q̃, 0 I q) .
Z
t1
I2 (p, q) = p (t0 )Q0 p(t0 ) +
pT Qpdt
t0
Z t1
(q T γq + 2q T β s 0 Q̃B s p)dt .
+
T
(.9)
t0
Since q uniquely defines p through (22), it follows that
min 0
I (z, d, u) = min I2 (p, q)
(.10)
z,u solve (11),F d=0
q
Note that I2 is convex and smooth in q. Thus, by
[Ioffe and Tikhomirov, 1974, p.81,Prop.1] q̂ is a minimizer
for I2 iff the Gateaux
derivative of I2 at q̂ vanishes:
d
I2 (pτ , q̂ + τ v)
= 0, where pτ denotes the solution
dτ
τ =0
of (22) with q = q̂ + τ v. Define
dp̂
= As1 p̂ + As2 q̂, p̂(t1 ) = F 0+ F 0 `,
dt
(.11)
dŵ
= −As1 0 ŵ + Qp̂ + B s 0 Q̃β s q̂, ŵ(t0 ) = Q0 p̂(t0 ) .
dt
It is not difficult to check that
Z t1
d
=2
v T γ q̂dt
I2 (pτ , q̂ + τ v)
dτ
t0
τ =0
Z t1
v T β s 0 Q̃B s p̂ − As2 0 ŵ dt .
+2
t0
So the Gateaux derivative of I2 at q̂ vanishes iff q̂ solves
the linear algebraic equation:
γ q̂ = −β s 0 Q̃B s p̂ + As2 0 ŵ .
(.12)
We note that, detγ = 0 in the general case. So (.12) is solvable iff the vector on the right hand side of (.12) belongs to
the range R(γ) of matrix γ. As Q̃ is a symmetric positive
definite matrix, then R(γ) = R(β s 0 ). It is clear that
−β s 0 Q̃B s p̂ ∈ R(β s 0 ). By definition, β s 0 = PN (As−1
)×
4
0
0
· · · × PN (A04 ) and As2 0 = PN (As−1
)
×
·
·
·
×
P
(A
)A
N
4
2
4
s0
s0
so A2 ŵ ∈ R(β ). Thus, the minimizer for I2 admits
representation:
q̂ = γ + As2 0 ŵ − β s 0 Q̃B s p̂ + v,
(.13)
where γv = 0. Let us plug (.13) into (.11). Then γv = 0
implies As2 v = 0. Thus p̂, ŵ do not depend on v and
so does I2 . Therefore, we can take v = 0 in (.13).
Now we split (.11) using the standard argument of the
linear regulator problem (see [Fleming and Rishel, 1975,
p.23,Exmpl.2.3,p.88]): plugging (.13) with v = 0 into (.11)
we find that ŵ = K p̂ and q̂ = γ + As2 0 K − β s 0 Q̃B s p̂
where K solves Riccati equation formulated in the theorem’s statement. Thus, q̂ is in the feed-back form. In
order to compute `-estimate
û which solves (11), we
recall that u = 0 I q and q = αs p + β s q so that
û = 0 I (αs p̂ + β s q̂). Using this representation and
definition of x̂ we find integrating by parts that: û(y) =
R t1 T 0 q̂ y dt = `T F F + x̂(t1 ). Let us compute `-error. We
t0
find plugging (.11) into (.9) and using (.4), (.10) that
σ̂(t1 , `) = I2 (p̂, q̂) = p̂T (t1 )F F + ` = `T F F + K(t1 )F F + `.
This completes the proof.
By statement 3) of Proposition 3 any solution p, q of (13)(14) may be represented as p = P0s p and q = αs p + β s q
where p, q solve (22). Using this it is straightforward to
check that I˜1 (z, u) = I2 (p, q) where:
148