Vladimir V Goncharov

Fátima F Pereira

Geometric Conditions for Regularity of Viscosity Solution to the Simplest Hamilton-Jacobi Equation *

Keywords: optimal time control problem, viscosity solution, eikonal equation, duality mapping, proximal normals, proximal regularity, Hölder continuity

 on well-posedness of the optimal time control problem with a constant convex dynamics (in a Hilbert space), we adapt one of the regularity conditions obtained there to a slightly more general problem, where nonaffine additive term appears. We prove existence and uniqueness of a minimizer in this problem as well as continuous differentiability of the value function (it can be seen as the viscosity solution to a Hamilton-Jacobi equation) near the boundary.

Introduction

Let us start with the first order partial differential equation in finite dimensions Γ (x, u (x) , ∇u (x)) = 0 [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF] where Γ : Ω × R × R n → R is a continuous function, nonlinear with respect to (w.r.t.) the third variable; Ω ⊂ R n is an open bounded region. Due to applications in optimal control and dynamical systems (1) is traditionally called (stationary) Hamilton-Jacobi equation. There are various notions of solutions to this equation. For instance, a function u

: Ω → R of class C Ω ∩ C 1 (Ω)
satisfying [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF] for all x ∈ Ω is said to be classical solution, while a Lipschitz continuous function u : Ω → R such that ( 1) holds for almost each (a.e.) x ∈ Ω is usually called generalized (or almost everywhere) solution. Speaking about solutions of (1) we always have in mind some prescribed boundary condition

u (x) = θ (x) , x ∈ ∂Ω, (2)
where θ : Ω → R is a (continuous) given function. Since in practice a classical solution to the boundary value problem ( 1)-( 2) often fails to exist while generalized solution may not be unique, another physically reasonable concept (so named viscosity solution) was introduced by M. Crandall and P.-L. Lions in 1983 (see [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF]) whereas similar constructions under different names were known earlier (see, e.g., [START_REF] Oleinik | On discontinuous solutions of non-linear differential equations[END_REF], [START_REF] Douglis | The continuous dependence of generalized solutions of nonlinear partial differential equations upon initial data[END_REF], [START_REF] Kružkov | Generalized solutions of the Hamilton-Jacobi equation of the Eikonal type[END_REF]). This concept was mainly based on the idea of "vanishing viscosity" in the sense that (under appropriate conditions) each viscosity solution is the uniform limit as ε → 0+ of the sequence of solutions u ε (•) to the respective boundary value problems for the nonlinear elliptic equations Γ (x, u, ∇u) -ε∆u = 0

where ∆ =

n i=1 ∂ 2 ∂x 2 i
is the Laplace operator (notice that (3) has a unique classical solution for each ε > 0 small enough due to Theorem 3.2 [START_REF] Kružkov | Generalized solutions of the Hamilton-Jacobi equation of the Eikonal type[END_REF]).

The exact definition of viscosity solution can be given either in terms of the suitable test functions (similarly as the notion of the generalized solutions of linear PDE in the sense of distrubutions), or by involving a Fréchet generalization of the gradient of a function at the point of nondifferentiability. It turned out that for each suitable boundary data θ (•) a (continuous) viscosity solution to the problem (1)-(2) exists, is unique and stable w.r.t. both θ (•) and Γ (•). Furthermore, it is consistent with other types of solutions. In particular, each viscosity solution belonging to C 1 (Ω) is classical one. For the main results of Theory of Viscosity Solutions, very developed and powerful field of the modern mathematics, we refer to [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF]- [START_REF] Barles | Solutions de Viscosité des Équations de Hamilton-Jacobi[END_REF] and to the bibliography therein. For a concise survey of viscosity solutions in finite dimensions see also the excellent tutorial lessons by A. Bressan [START_REF] Bressan | Hamilton-Jacobi equations and optimal control[END_REF].

Afterwards, the concept and the main results concerning viscosity solutions were generalized to Banach spaces with the Radon-Nikodym property (see [START_REF] Crandall | Hamilton-Jacobi equations in infinite dimensions. I: Uniqueness of viscosity solutions[END_REF], [START_REF] Crandall | Hamilton-Jacobi equations in infinite dimensions. II: Existence of viscosity solutions[END_REF]), in particular, to Hilbert spaces. Notice that although the definition based on the Fréchet sub-and superdifferentials remains the same, the interpretation of viscosity solutions via "vanishing viscosity" is no longer valid in infinite dimensions. The motivation, however, comes now from the Theory of Differential Games.

In our paper we deal only with the case when the hamiltonian Γ in ( 1) does not depend of x neither u, and is convex w.r.t. the third variable. Already S. N. Kružkov studied in [START_REF] Kružkov | Generalized solutions of the Hamilton-Jacobi equation of the Eikonal type[END_REF] such Hamilton-Jacobi equations arizing from the geometric optics. For instance, when n = 3 and Γ(x, u, ξ) = |ξ|a, with a constant a > 0, one has so called eikonal equation describing the propagation of a light wave from a point source placed at the origin in homogeneous medium with refraction index 1/a. If, instead, this medium is anisotropic and has constant coefficients of refraction of light rays parallel to the coordinate axes (say c i ) then the propagation of light can be described by the (more general) elliptic equation

n i=1 c 2 i u 2 xi -1 = 0. (4)
If, besides that, the medium moves with a constant velocity -→ v then the equation contains already a linear additive term and admits the form

n i=1 c 2 i u 2 xi + 2 c - → v , ∇u -1 = 0, (5)
where c means the speed of the light in a vacuum. In general, denoting by F the closed convex hull of the set of zeros

{ξ ∈ R n : Γ (ξ) = 0}
and assuming F to be bounded with intF = ∅ (the hamiltonians in ( 4) and ( 5) satisfy these conditions), the equation ( 1) can be reduced to

ρ F (∇u (x)) -1 = 0, (6)
where ρ F (•) is the Minkowski functional (gauge function) associated to F ,

ρ F (ξ) := inf {λ > 0 : ξ ∈ λF } .
More precisely, it was proved in [START_REF] Cardaliaguet | Geometric restrictions for the existence of viscosity solutions[END_REF] that under appropriate conditions involving a kind of geometric compatibility of F , θ (•) and the domain Ω the (unique) viscosity solution û (•), û | ∂Ω = θ, of ( 6) is the viscosity solution of the problem (1)-(2) (belonging to the space W 1,∞ (Ω)) and vice versa. Furthermore, this viscosity solution can be given by the formula

û (x) = inf y∈C {ρ F 0 (x -y) + θ (y)} (7) whenever θ : R n → R is a Lipschitz continuous function such that ∇θ (x) ∈ intF for a.e. x ∈ R n . ( 8
)
Here C := R n \Ω and F 0 is the polar set for F . Let now H be a Hilbert space with the norm • and the inner product •, • . Then the convolution [START_REF] Crandall | Hamilton-Jacobi equations in infinite dimensions. I: Uniqueness of viscosity solutions[END_REF] remains the unique viscosity solution to the equation ( 6) with the boundary data u (x) = θ (x), x ∈ C, whenever the slope condition

θ (x) -θ (y) < ρ F 0 (x -y) ∀x, y ∈ C (9)
holds. Notice that the inequality (9) follows from (8) in finite dimensions while in an arbitrary Hilbert space it can be deduced from the inclusion ∂ c θ (x) ⊂ F ,

x ∈ H. Here ∂ c is the Clarke generalized gradient of a Lipschitz continuous function. So, we are interested in regularity properties of the function [START_REF] Crandall | Hamilton-Jacobi equations in infinite dimensions. I: Uniqueness of viscosity solutions[END_REF], which was well studied when θ ≡ 0 (see [START_REF] Colombo | Variational Analysis for a class of minimal time functions in Hilbert spaces[END_REF], [START_REF] Goncharov | Neighbourhood retractions of nonconvex sets in a Hilbert space via sublinear functionals[END_REF], [START_REF] Goncharov | Geometric conditions for regularity in a time-minimum problem with constant dynamics[END_REF]). In the latter case let us notice the following:

1) existence and regularity of the (Fréchet) gradient ∇û (x) depends on uniqueness (in infinite dimensions also on existence) of a minimizer in [START_REF] Crandall | Hamilton-Jacobi equations in infinite dimensions. I: Uniqueness of viscosity solutions[END_REF];

2) it is not possible that ∇û (x) exists everywhere out of C unless some special situations;

3) the function û (x) can be interpreted as the minimal time for achieving the closed set C from x ∈ H\C by trajectories of the differential inclusion ẋ (t) ∈ -F 0 .

Taking into account 1) and 2) it is natural to study the regularity only in an (open) neighbourhood of C (target set due to 3). If F = B is the closed unit ball centred in the origin then û (•) is nothing else than the distance from C, and the minimizers in [START_REF] Crandall | Hamilton-Jacobi equations in infinite dimensions. I: Uniqueness of viscosity solutions[END_REF] are the usual metric projections onto C. In this case the (necessary and sufficient) condition guaranteeing both well-posedness of the problem and the (Lipschitz) continuity of the gradient ∇û (x) near C is so named ϕ -convexity (or proximal smoothness) of the set C well studied up to now (see survey [START_REF] Colombo | Prox-regular sets and applications[END_REF] and the bibliography therein).

As concerns an arbitrary gauge F (and θ ≡ 0) then in [START_REF] Goncharov | Neighbourhood retractions of nonconvex sets in a Hilbert space via sublinear functionals[END_REF], [START_REF] Goncharov | Geometric conditions for regularity in a time-minimum problem with constant dynamics[END_REF] two different hypotheses are given, under which both a unique minimizer in [START_REF] Crandall | Hamilton-Jacobi equations in infinite dimensions. I: Uniqueness of viscosity solutions[END_REF] (that is a point on the boundary ∂C attained from x for the minimal time) and the gradient ∇û (x) are (Hölder) continuous in a neighbourhood of C. It turns out that one of these hypotheses (based on certain ballance between external normals to the sets C and F ) can be adapted to the case of a Lipschitz continuous perturbation θ (•).

We start in Section 2 with the basic definitions and an auxiliary statement. Then, in Section 3, we study the mathematical programming problem [START_REF] Crandall | Hamilton-Jacobi equations in infinite dimensions. I: Uniqueness of viscosity solutions[END_REF] from the viewpoint of the existence, uniqueness and the (Lipschitz) regularity of minimizers near the set C. The geometric condition ensuring such well-posedness is emphasized here. Finally, in Section 4 we examine the (Fréchet) differentiability of the value function û (•) and justify the (Hölder) continuity of its gradient also under the assumption that either F 0 or the restriction θ | C is smooth.

Preliminaries

Given a convex closed bounded set F ⊂ H with 0 ∈ intF we consider the so called duality mapping J F : ∂F 0 → ∂F , which associates to each ξ * ∈ ∂F 0 the set of (normalized) linear functionals that support F 0 at ξ * ,

J F (ξ * ) := {ξ ∈ ∂F : ξ, ξ * = 1} .
In other words, J F (ξ * ) = N F 0 (ξ * ) ∩ ∂F where N F 0 (ξ * ) is the normal cone to the polar F 0 at ξ * . It can be interpreted also as the subdifferential ∂ρ F 0 (ξ * ) in the sense of Convex Analysis. For each dual pair (ξ, ξ * ), i.e., such that ξ ∈ ∂F , ξ * ∈ ∂F 0 and ξ, ξ * = 1 let us define the modulus of rotundity (see [START_REF] Goncharov | Neighbourhood retractions of nonconvex sets in a Hilbert space via sublinear functionals[END_REF])

C F (r, ξ, ξ * ) := inf { ξ -η, ξ * : η ∈ F, ξ -η ≥ r} , r > 0.
If the set F is strictly convex (rotund ) at ξ w.r.t. ξ * , i.e., C F (r, ξ, ξ * ) > 0 ∀r > 0, then ξ is an exposed point of F and, in particular, ξ is the unique element of J F (ξ * ). So, in this case ξ is well defined whenever ξ * is fixed. Furthermore, given a set U ⊂ ∂F 0 we say that F is uniformly rotund w

.r.t. U if inf C F (r, ξ, ξ * ) : ξ ∈ J F (ξ * ) , ξ * ∈ U > 0 ∀r > 0.
By [START_REF] Goncharov | Geometric conditions for regularity in a time-minimum problem with constant dynamics[END_REF]Proposition 2.1] this property is equivalent to the uniform continuity of J F (•) in the following sense

sup η∈J F (η * ) J F (ξ * ) -η → 0 as ξ * -η * → 0, ξ * ∈ U , η * ∈ ∂F 0 (10)
(we clearly identify J F (ξ * ) with its element whenever it is a singleton). Uniform rotundity implies also the existence and the uniform continuity on U of the Fréchet gradient ∇ρ F 0 (ξ * ).

Besides the concepts of Convex Analysis above we will use the following notations. For a lower semicontinuous function ϕ : H → R ∪ {+∞} we denote by ∂ p ϕ (x), ∂ l ϕ (x), ∂ -ϕ (x) and ∂ c ϕ (x) the proximal, limiting (Mordukhovich), Fréchet and Clarke subdifferential, respectively, at a point x, ϕ (x) < +∞. All the definitions and the basic facts of the calculus for non convex sets can be found, e.g., in [START_REF] Clarke | Nonsmooth Analysis and Control Theory[END_REF]. Here we observe only that the inclusions

∂ p ϕ (x) ⊂ ∂ -ϕ (x) ⊂ ∂ l ϕ (x) ⊂ ∂ c ϕ (x) (11)
always hold, while one of the reverse inclusions takes place whenever ϕ (•) is regular at x in some sense. For instance, ϕ (•) is said to be proximal (Clarke)

regular at x if ∂ p ϕ (x) = ∂ l ϕ (x) (respectively, ∂ -ϕ (x) = ∂ c ϕ (x)).
If C ⊂ H is a nonempty closed set then the notion of some kind of normal cone to C at a point x ∈ C can be given as the respective subdifferential of the indicator function I C (•) equal to 0 on C and to +∞ elsewhere. In particular, the proximal normal cone N p C (x) := ∂ p I C (x). Further on we denote by ∂ * C := {x ∈ ∂C : N p C (x) = {0}} the effective boundary, which is dense in ∂C.

Returning to the minimization problem [START_REF] Crandall | Hamilton-Jacobi equations in infinite dimensions. I: Uniqueness of viscosity solutions[END_REF] let us formulate first an approximation result, which is crucial for what follows. It can be proved similarly as Lemma 5.1 [START_REF] Goncharov | Neighbourhood retractions of nonconvex sets in a Hilbert space via sublinear functionals[END_REF] by using the Ekeland's variational principle as well as the fuzzy sum rule for the proximal subdifferentials (see [START_REF] Clarke | Nonsmooth Analysis and Control Theory[END_REF]Theorem 1.8.3]).

} such that v n ∈ ∂ p (θ | C ) (x ′ n ), ξ n ∈ ∂ρ F 0 (x -x ′′ n ) and x ′ n -x n + x ′′ n -x n → 0, v n -ξ n → 0 as n → ∞. Here θ | C := θ + I C . Notice that if the points x ′ n are such that ∂ p (θ | C ) (x ′ n ) = ∂ p θ (x ′ n )+N p C (x ′ n ) and ∂ p θ (x ′
n ) ⊂ γF for some 0 < γ < 1, then without loss of generality we can assume that x ′ n ∈ ∂ * C and v n ∈ ∂F .

Existence, uniqueness and regularity of minimizers

Our standing hypothesis in what follows will be a slightly strengthened condition than ( 9):

(H) there exists 0 < γ <

1 F F 0 such that θ (x) -θ (y) ≤ γρ F 0 (x -y) ∀x, y ∈ C,
where

F := sup { ξ : ξ ∈ F }.
Hence θ (•) is lipschitzean on C with the constant γ F . Given now an arbitrary point x 0 ∈ ∂C let us emphasize the main local assumptions, under which the well-posedness results hold:

(H 1 ) the mapping x → J F 0 (∂ p (θ | C ) (x) ∩ ∂F ) is single-valued and lipschitzean (with a constant L = L (x 0 ) > 0) on the set

C δ (x 0 ) := {x ∈ ∂ * C : x -x 0 ≤ δ} , δ > 0;
(H 2 ) in the δ-neighbourhood of x 0 the sum rule

∂ p (θ | C ) (x) = ∂ p θ (x) + N p C (x)
takes place;

(H 3 ) F 0 is uniformly rotund w.r.t. the set

U δ (x 0 ) := ∂F ∩ x∈C δ (x0) ∂ p (θ | C ) (x) . ( 12
)
For each x ∈ H we denote by π F,θ C (x) the (possibly empty) set of all minimizers of the function y → ρ F 0 (xy) + θ (y) on C.

Theorem 1 Under the hypotheses (H 1 ) -(H 3 ) there exists a neighbourhood U (x 0 ) such that the mapping x → π F,θ C (x) is single-valued and continuous on U (x 0 ).

Proof. Let us give a sketch of the proof. Taking without loss of generality δ > 0 such that δγ F < 1γ F F 0 /L, let us set

U (x 0 ) := x ∈ H : x -x 0 < 1 -γ F F 0 δ 2 F F 0 , û (x) < û (x 0 ) + 1 -γ F F 0 L -δγ F . ( 13
)
Fix x ∈ U (x 0 ) \C and a minimizing sequence

{x n } ⊂ C of y → ρ F 0 (x -y) + θ (y). Let us choose {x ′ n } ⊂ C, {x ′′ n } ⊂ H, v n ∈ ∂ p (θ | C ) (x ′ n ) and ξ n ∈ ∂ρ F 0 (x -x ′′ n ) as in Lemma 1.
Our goal is to prove that {x ′ n } (hence {x n } as well) is a Cauchy sequence.

To this end we show, first, that x ′ nx 0 ≤ δ. It follows then from (H 2 ) and from the remark after Lemma 1 that x ′ n ∈ C δ (x 0 ) and v n ∈ ∂F for all n large enough. Considering a (nonincreasing) sequence ν n → 0+ such that Lemma 1). Then using the hypothesis (H 3 ) together with the property (10) (applied to the gauge F 0 ) gives that

x ′ n -x n + x ′′ n -x n ≤ ν n ; (14)
ρ F 0 (x -x ′ n ) + θ (x ′ n ) ≤ û (x) + ν n ; (15) v n -ξ n ≤ ν n , n = 1, 2, ... (see
β n := sup ξ-η ≤νn sup { J F 0 (ξ) -J F 0 (η) : ξ ∈ ∂F , η ∈ U δ (x 0 )} → 0 as n → ∞.
Taking into account that v n ∈ U δ (x 0 ) (see [START_REF] Douglis | The continuous dependence of generalized solutions of nonlinear partial differential equations upon initial data[END_REF]) and ξ n ∈ ∂F we obtain

J F 0 (ξ n ) -J F 0 (v n ) ≤ β n ,
and, consequently, by (H 1 )

J F 0 (ξ n ) -J F 0 (ξ m ) ≤ 2β n + L x ′ n -x ′ m (16)
for all m ≥ n ≥ 1.

On the other hand, by the elementary properties of the convex subdifferentials and duality mappings we find that (x -

x ′′ n ) /ρ F 0 (x -x ′′ n ) = J F 0 (ξ n ),

and hence

x ′′ n -x ′′ m = ρ F 0 (x -x ′′ n ) J F 0 (ξ n ) -ρ F 0 (x -x ′′ m ) J F 0 (ξ m ) ≤ ≤ ρ F 0 (x -x ′′ n ) J F 0 (ξ n ) -J F 0 (ξ m ) + + |ρ F 0 (x -x ′′ n ) -ρ F 0 (x -x ′′ m )| F 0 . ( 17
)
The terms |ρ

F 0 (x -x ′′ n ) -ρ F 0 (x -x ′′ m )| and ρ F 0 (x -x ′′ n
) can be approximately estimated by x ′ nx ′ m and by û (x)û (x 0 ), respectively (we use for that the inequalities ( 14) and ( 15)). Hence, taking into account also [START_REF] Oleinik | On discontinuous solutions of non-linear differential equations[END_REF] we deduce from (17) that

x ′ n -x ′ m ≤ x ′′ n -x ′′ m + 2ν n ≤ µ n + + L (µ ′ n + û (x) -û (x 0 ) + δγ F ) + γ F F 0 x ′ n -x ′ m
for all m ≥ n ≥ 1, where {µ n } and {µ ′ n } are some sequences, converging to zero. We conclude the proof recalling the definition of the neighbourhood U (x 0 ). Thus, the limit x := lim n→∞ x ′ n will be (unique) minimizer from π F,θ C (x). The continuous dependence of this singleton on x ∈ U (x 0 ) also follows.

Remark 1. In fact, adapting the proof of Theorem 3.1 [START_REF] Goncharov | Geometric conditions for regularity in a time-minimum problem with constant dynamics[END_REF] we can show that the mapping x → π F,θ C (x) is locally lipschitzean on the same neighbourhood (13) with the Lipschitz constant tending to +∞ as x tends to the boundary ∂U (x 0 ).

Differentiability of the viscosity solution

We announce, first, a result on subdifferential regularity of the function ( 7) at a fixed point x / ∈ C, similar to Proposition 5.1 [START_REF] Goncharov | Geometric conditions for regularity in a time-minimum problem with constant dynamics[END_REF]. We see that the regularity relies upon well-posedness of the minimizers studied in the previous section.

Theorem 2 Let us fix x ∈ H\C and assume that

• π F,θ
C (y) is a singleton for each y, xy ≤ δ, δ > 0;

• the restriction θ | C is proximally regular at x := π F,θ C (x); • the following "centred" Hölder property

π F,θ C (x) -π F,θ C (y) ≤ K x -y β ∀y, x -y ≤ δ,
holds with an exponent 1/2 < β ≤ 1 and a constant K = K (x) > 0.

Then the function û (•) is Clarke regular at x. More precisely,

∂ c û (x) = ∂ -û (x) = ∂ρ F 0 (x -x) ∩ ∂ -(θ | C ) (x) = ∅. (18)
Recalling Theorem 1 and Remark 1 we immediately obtain from the statement above that under the hypotheses (H 1 ) -(H 3 ) the viscosity solution û (•) is Clarke regular and (18) holds for all x ∈ U (x 0 ) where the neighbourhood U (x 0 ) is defined by [START_REF] Goncharov | Neighbourhood retractions of nonconvex sets in a Hilbert space via sublinear functionals[END_REF]. Thus, for the (Fréchet) continuous differentiability it suffices to require that Φ(x) := ∂ρ F 0 (xx) ∩ ∂ -(θ | C ) (x) is a singleton continuously depending on x ∈ U (x 0 ). However, this is difficult to verify directly because Φ(x) (which is nothing else than the Fréchet gradient ∇û (x)) depends on the point x through a priori unknown function π F,θ C (•). On the other hand, the latter condition (single-valuedness and continuity of the mapping Φ(•)) splits into two different hypotheses regarding the smoothness either of both the function θ(•) and the set C, or of the polar gauge F 0 . Moreover, such hypotheses can be given plainly in terms of boundary points of the sets C and F 0 . Notice that C is said to be with smooth boundary near x 0 if N l C (x) ∩ ∂B is a singleton (say {n C (x)}) continuously depending on x ∈ ∂C, xx 0 ≤ δ, δ > 0. Also the smoothness of F 0 can be equivalently substituted by the rotundity assumption for F .

Thus, we arrive at the following result.

Theorem 3 Given x 0 ∈ ∂C and δ > 0 let us assume the hypotheses (H 1 ) -(H 3 ). Suppose also that at least one of the following two conditions holds:

(i) C has smooth boundary, and

θ (•) is of class C 1 near x 0 ; (ii) F is rotund w.r.t. each ξ * ∈ J F 0 (∂ p (θ | C ) (x)), x ∈ ∂C with x -x 0 ≤ δ.
Then û (•) is Fréchet continuously differentiable on the neighbourhood U (x 0 ). Furthermore, in the first case we have

∇û (x) = ∇θ (x) + λn C (x) ,
where λ = λ (x) > 0 is the unique positive root of the equation

ρ F (∇θ (x) + λn C (x)) = 1,
while in the second ∇û (x) = ∇ρ F 0 (xx) .

Here x := π F,θ C (x), x ∈ U (x 0 ), as before.

Proof. If the condition (i) is fulfilled then taking into account that π F,θ C (•) is single-valued and continuous on U (x 0 ), and that ∂ l (θ | C ) (x) = ∇θ (x) + N l C (x) = {∇θ (x) + λn C (x) : λ ≥ 0} whenever x ∈ U (x 0 ), we obtain that the intersection in (18) reduces to the singleton {∇θ (x) + λ (x) n C (x)} (see also [START_REF] Dal Maso | A Lipschitz selection from the set of minimizers of a nonconvex functional of the gradient[END_REF]). The continuity of ∇û (•) can be shown now by the standart implicit function argument. Under the alternative assumption (ii) it sufices to observe that due to a necessary condition of optimality (in the proximal form) the (unique) minimizer x = π F,θ C (x) must satisfy the relationship

∂ p (θ | C ) (x) ∩ N F 0 x - x ρ F 0 (x -x) ∩ ∂F = ∅. (19)
Then, it follows from (19) that

x-x ρ F 0 (x-x) ∈ J F 0 (ξ) for some ξ ∈ ∂ p (θ | C ) (x). Therefore, if x ∈ U (x 0 ) then x is closed to x 0 as well, and taking ξ * = x-x ρ F 0 (x-x) we deduce from (ii) that ρ F 0 (•) is (Fréchet) continuously differentiable at ξ * . So, the intersection in (18) reduces to {∇ρ F 0 (xx)}, and the continuity w.r.t.

x also follows.

Remark 2. If in addition to the hypothesis (i) in Theorem 3 we assume that both unit normal vector n C (•) and the gradient ∇θ (•) are Hölder continuous with an exponent 0 < α ≤ 1 on a δ-neighbourhood of x 0 then ∇û (•) will be also Hölder continuous near x 0 with the same exponent (we say that û (•) is of class C 1,α loc on U (x 0 )). One can derive the Hölder inequality for ∇û (•) by using Theorem 3 and the estimates for the Hausdorff distance between the polars for convex solids (see Lemma 2 [START_REF] Dal Maso | A Lipschitz selection from the set of minimizers of a nonconvex functional of the gradient[END_REF]).

Lemma 1

 1 Let C ⊂ H be a nonempty closed set, and θ : H → R ∪ {+∞} be a lower semicontinuous function, lipschitzean on C. If x ∈ H\C and {x n } ⊂ C is a minimizing sequence for the function y → ρ F 0 (xy) + θ (y) on C, then there exist another minimizing sequence {x ′ n } ⊂ C and sequences {x ′′ n }, {v n } and {ξ n

* Work is realized in framework of the project "Variational Analysis: Theory and Applications" (PTDC/MAT/111809/2009) financially supported by Fundação para Ciência e Tecnologia (FCT), the Portuguese institutions COMPETE, QREN and the European Regional Development Fund (FEDER)

† Corresponding author.


Martin Gugat
email: gugat@math.fau.de

Markus Dick
email: dick@math.fau.de

Günter Leugering ⋆⋆

Stabilization of the Gas Flow in Star-Shaped Networks by Feedback Controls with Varying Delay

Keywords: boundary feedback stabilization, Euler equations, gas network, Lyapunov function, star-shaped network, time-varying delay

We consider the subcritical gas flow through star-shaped pipe networks. The gas flow is modeled by the isothermal Euler equations with friction. We stabilize the isothermal Euler equations locally around a given stationary state on a finite time interval. For the stabilization we apply boundary feedback controls with time-varying delays. The delays are given by C 1 -functions with bounded derivatives. In order to analyze the system evolution, we introduce an L 2 -Lyapunov function with delay terms. The boundary controls guarantee the exponential decay of the Lyapunov function with time.

Introduction

Recently, there has been intense research on the system dynamics in gas networks (see e.g. [START_REF] Banda | Coupling conditions for gas networks governed by the isothermal Euler equations[END_REF][START_REF] Banda | Gas flow in pipeline networks[END_REF][START_REF] Dick | Classical solutions and feedback stabilization for the gas flow in a sequence of pipes[END_REF][START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF][START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF][START_REF] Herty | A new model for gas flow in pipe networks[END_REF][START_REF] Osiadacz | Simulation and Analysis of Gas Networks[END_REF]). Due to the pipe wall friction, there is a loss of pressure along the pipe. A common model for the gas flow in pipes is the isothermal Euler equations with friction, a 2 × 2 PDE system of balance laws (see [START_REF] Banda | Coupling conditions for gas networks governed by the isothermal Euler equations[END_REF]). We study the isothermal Euler equations on a star-shaped network of N (N ≥ 2) pipes that meet at a central node ω. The flow at the node ω is governed by the continuity of the density and conservation of mass (see [START_REF] Dick | Classical solutions and feedback stabilization for the gas flow in a sequence of pipes[END_REF]). Our main result, stated in Theorem 1, is a method to stabilize the gas flow locally around a given stationary subcritical state on a finite time interval. To do so, we use boundary feedback controls with varying delays at the pipe ends which are not at the node ω (see ( 16)). In order to measure the system evolution, we introduce a network Lyapunov function (see (22)) which is the sum of a weighted and squared L 2 -norm (see (20)) and a delay term (see (21)) for each pipe. The feedback controls guarantee the exponential decay of the Lyapunov function with time (see (36)) and, hence, the exponential stability of the system. In contrast to a previous work which studies only the case of constant delays (see [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF]), the novelty of this paper is that we consider nonconstant, i.e. timedependent, delays. This is very important for the daily operation of real gas networks. E.g., in the control of such networks via electrical and mechanical systems, nonconstant time delays often appear (see [START_REF] Bogel | Method and apparatus for control of pipeline compressors[END_REF]). This paper is structured as follows: In Sect. 2 we give the network notation, consider the isothermal Euler equations in terms of the physical and characteristic variables and present the coupling conditions at the node ω. In Sect. 3 stationary and nonstationary states are studied. The stabilization method, i.e. the feedback controls, the Lyapunov function and the exponential stability result (Theorem 1), are stated in Sect. 4. In Sect. 5 we prove Theorem 1. The weighted and squared L 2 -norm from (20) for the Euler equations has first been presented in [START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF]. It is an extension of the Lyapunov function introduced in [START_REF] Coron | A strict Lyapunov function for boundary control of hyperbolic systems of conservation laws[END_REF]. Delay terms of the form (21) have been presented in [START_REF] Nicaise | Stability of the heat and of the wave equations with boundary time-varying delays[END_REF] for the time-delayed stabilization of the wave equation. Related questions of the stabilization of the wave equation are e.g. studied in [START_REF] Gugat | Boundary feedback stabilization by time delay for one-dimensional wave equations[END_REF][START_REF] Gugat | Stars of vibrating strings: switching boundary feedback stabilization[END_REF][START_REF] Wang | Wave equation stabilization by delays equal to even multiples of the wave propagation time[END_REF].

Gas Flow in a Star-Shaped Pipe Network

In this section we consider the gas flow in a star-shaped pipe network. First, we give the network notation. Then, we present the isothermal Euler equations in terms of the physical variables (see [START_REF] Banda | Coupling conditions for gas networks governed by the isothermal Euler equations[END_REF]) and in terms of the characteristic variables (see [START_REF] Bogel | Method and apparatus for control of pipeline compressors[END_REF]). The coupling conditions at the central node of the network are stated in (5).

Network Notation

We consider a star-shaped network of N (N ≥ 2) cylindrical pipes with the same diameter δ > 0 that meet at a central node ω. We define the index set I = {1, ..., N } and number the pipes from pipe 1 to pipe N . Variables referring to pipe i (i ∈ I) are denoted with a superscript (i). We model the pipes by a one-dimensional space model and parameterize the length L (i) > 0 of pipe i by the space interval [0, L (i) ] such that the end x = 0 is at the node ω. We consider the system on a finite time interval [0, T ] with T > 0.

Isothermal Euler Equations in Physical Variables

A common model for the system dynamics in gas pipes is the isothermal Euler equations with friction, a hyperbolic 2 × 2 system of balance laws (see [START_REF] Banda | Coupling conditions for gas networks governed by the isothermal Euler equations[END_REF][START_REF] Banda | Gas flow in pipeline networks[END_REF][START_REF] Herty | A new model for gas flow in pipe networks[END_REF][START_REF] Osiadacz | Simulation and Analysis of Gas Networks[END_REF]): For pipe i (i ∈ I) they have the following form on [0, T ] × [0, L (i) ]:

   ∂ t ρ (i) (t, x) + ∂ x q (i) (t, x) = 0 , ∂ t q (i) (t, x) + ∂ x (q (i) (t,x)) 2 ρ (i) (t,x) + a 2 ρ (i) (t, x) = -θ 2 q (i) (t,x) |q (i) (t,x)| ρ (i) (t,x) (1)
where ρ (i) (t, x) > 0 is the density of the gas and q (i) (t, x) = 0 the mass flux. The sign of q (i) depends on the direction of the gas flow. It is positive if the gas in pipe i flows away from the node ω. The constant a > 0 is the sonic speed in the gas and the constant θ = ν/δ is the quotient of the friction factor ν > 0 and the pipe diameter δ > 0. The first equation in [START_REF] Banda | Coupling conditions for gas networks governed by the isothermal Euler equations[END_REF] states the conservation of mass and the second equation is the momentum equation. In this paper we study subsonic or subcritical C 1 -states, i.e. C 1 -states with |q (i) |/ρ (i) < a for all i ∈ I.

The equations ( 1) have the same form for each pipe. However, our calculations would also be true if we had different pressure laws on different pipes.

Isothermal Euler Equations in Characteristic Variables

The equations ( 1) can be transformed to the Riemann invariants (characteristic variables)

R (i) ± = -q (i) /ρ (i) ∓ a ln(ρ (i) ) . (2)
For the calculation of the Riemann invariants see [START_REF] Dick | Classical solutions and feedback stabilization for the gas flow in a sequence of pipes[END_REF][START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF]. In terms of R (i)

± the system (1) has the form

∂ t R (i) + R (i) - + λ (i) + 0 0 λ (i) - ∂ x R (i) + R (i) - = - θ 8 (R (i) + + R (i) -) |R (i) + + R (i) -| 1 1 (3)
with the eigenvalues

λ (i) ± = q (i) ρ (i) ± a = - R (i) + + R (i) - 2 ± a . (4)
In the subsonic case, λ

+ is strictly positive and λ (i) -strictly negative.

Coupling Conditions

At the node ω, i.e. at the ends x = 0 of the pipes, we have the following coupling conditions in terms of the physical variables (t ∈ [0, T ]) (see [START_REF] Banda | Coupling conditions for gas networks governed by the isothermal Euler equations[END_REF][START_REF] Marigo | Entropic solutions for irrigation networks[END_REF]):

ρ (1) (t, 0) = ρ (i) (t, 0) (i ∈ I\{1}) and i∈I q (i) (t, 0) = 0 . (5)
The first equation in [START_REF] Dick | Classical solutions and feedback stabilization for the gas flow in a sequence of pipes[END_REF] says that the density is continuous at the node ω. Due to the parameterization of the pipes, the second equation in [START_REF] Dick | Classical solutions and feedback stabilization for the gas flow in a sequence of pipes[END_REF] means that the total ingoing mass flux is equal to the total outgoing mass flux at ω. The conditions (5) can equivalently be stated in terms of R

(i) ± as (t ∈ [0, T ]) R (1)
+ (t, 0), ..., R

(N ) + (t, 0) = R (1)
-(t, 0), ..., R (N ) -(t, 0) A ω (6)
with the orthogonal, symmetric (N ×N )-matrix A ω = (a kl ) N k,l=1 with the entries

a kk = (N -2)/N (k ∈ I) and a kl = -2/N (k, l ∈ I, k = l) . ( 7
)
3 Stationary and Nonstationary States

Stationary States

The existence and behavior of stationary solutions of the isothermal Euler equations, i.e. solutions which do not explicitly depend on the time t, is studied in [START_REF] Dick | Classical solutions and feedback stabilization for the gas flow in a sequence of pipes[END_REF][START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF]. We denote the stationary variables as ρ(x), q(x), R± (x) and λ± (x). In [START_REF] Dick | Classical solutions and feedback stabilization for the gas flow in a sequence of pipes[END_REF][START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF] it is shown that q(x) is constant along a pipe and ρ(x) is strictly monotonically decreasing in the direction of the gas flow. Furthermore, a stationary subsonic C 1 -solution of the isothermal Euler equations exists on the whole pipe if the pipe length is shorter than a critical length (see [START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF]). The critical length depends on the inflow density, the mass flux, the friction factor and the pipe diameter. For typical high-pressure gas pipes the critical length is around 180km (see [START_REF] Dick | Classical solutions and feedback stabilization for the gas flow in a sequence of pipes[END_REF]).

For the system (1) on a star-shaped network with the conditions ( 5), the existence and behavior of stationary subsonic C 1 -solutions is in detail discussed in [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF]. In the following we summarize the main results from [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF]: The stationary mass fluxes q(i) = 0 are constant and have to satisfy i∈I q(i) = 0. At the node ω, we have a constant density ρω with ρ(i) (0) = ρω and |q (i) |/ρ ω < a for all i ∈ I. Furthermore, the lengths of the pipes with positive mass flux, i.e. with q(i) > 0, have to satisfy

L (i) < 1 θ a 2 ρ2 ω (q (i) ) 2 -1 + 2 ln q(i) aρ ω . ( 8
)

Nonstationary States

Assume that on the star-shaped network we have a given stationary subsonic state (ρ (i) (x), q(i) ) with the corresponding Riemann invariants ( R(i) 2), ( 4)). We define the numbers (see (4))

+ (x), R(i) -(x)) ∈ (C 1 ([0, L (i) ])) 2 and the eigenvalues λ(i) (x) (i ∈ I) (see (
σ (i) = sign(q (i) ) = -sign( R(i) + + R(i) -) ∈ {-1, 1} . (9)
Now we consider a nonstationary state ( R(i)

+ (x) + r (i) + (t, x), R(i) -(x) + r (i) -(t, x)) on [0, T ] × [0, L (i) ] in a local neighborhood of ( R(i) + , R(i) -)
. That is we assume the given stationary state R(i) ± (x) to be slightly perturbed by r (i) ± (t, x) where the

C 1 -norms ||r (i) ± || C 1 ([0, T ]×[0, L (i) ]
) are small. In particular, we suppose that the mass flux directions of the nonstationary state are the same as of the stationary state, i.e. (i ∈ I) (see [START_REF] Gugat | Stars of vibrating strings: switching boundary feedback stabilization[END_REF])

σ (i) = -sign( R(i) + + R(i) -+ r (i) + + r (i) -) . ( 10
)
The stabilization method presented in Sect. 4 guarantees that the absolute values of r (i)

± are small enough (see ( 24), ( 35)), such that (10) holds and such that the direction of the mass fluxes does not change during the stabilization process. From (3) we obtain the following quasilinear system for r (i)

± (t, x) (i ∈ I)    ∂ t r (i) + + ( λ(i) + - r (i) + +r (i) - 2 ) ∂ x r (i) + = (r (i) + + r (i) -)(-K (i) + + σ (i) θ 8 (r (i) + + r (i) -)) , ∂ t r (i) -+ ( λ(i) -- r (i) + +r (i) - 2 ) ∂ x r (i) -= (r (i) + + r (i) -)(-K (i) -+ σ (i) θ 8 (r (i) + + r (i) -)) (11) on [0, T ] × [0, L (i)
] with the strictly positive C 1 -functions

K (i) ± (x) = θ 8 | R(i) + (x) + R(i) -(x)| 4a ∓ ( R(i) + (x) + R(i) -(x)) 2a ∓ ( R(i) + (x) + R(i) -(x)) > 0 . ( 12
)
The linearity of the equation ( 6) implies that for r

-(t, 0), ..., r (N ) -(t, 0) A ω (13)
with the matrix A ω as in [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF].

Feedback Stabilization with Time-Varying Delay

In this section we present a method for the stabilization of the system (11) on [0, T ] × [0, L (i) ] with time-delayed feedbacks. The corresponding boundary controls are given in (16). In order to measure the system evolution, we define the Lyapunov function F(t) in ( 22) with the weighted and squared L 2 -norms E (i) (t) in (20) and the delay terms D (i) (t) in (21). Theorem 1 states the existence of a unique C 1 -solution of [START_REF] Marigo | Entropic solutions for irrigation networks[END_REF] with small C 1 -norm (see ( 35)) for which F(t) decays exponentially with time (see (36)). Hence, Theorem 1 gives the exponential stability of the system (11) around (r

(i) + , r (i)
-) = (0, 0).

Boundary Feedback Controls with Time-Varying Delay

Let a finite time T > 0 and a stationary subsonic state ( R(i)

+ (x), R(i) -(x)) ∈ (C 1 ([0, L (i) ]
)) 2 on the star-shaped network be given with the eigenvalues λ(i) ± (x) as in (4) (i ∈ I). Define the numbers σ (i) ∈ {-1, 1} as in ( 9) and the functions

K (i) ± (x) ∈ C 1 ([0, L (i) ]) as in (12). Let functions τ (i) (t) ∈ C 1 ([0, T ]) (i ∈ I) be given that satisfy (t ∈ [0, T ]): 0 < τ (i) (t) < T 2 and | d dt τ (i) (t)| < 1 . ( 14
)
Define the constants

τ (i) = max t∈[0, T ] τ (i) (t), τ (i) = max t∈[0, T ] | d dt τ (i) (t)| and τ max = max i∈I τ (i) .
(15) For a nonstationary state R(i)

± (x) + r (i) ± (t, x) on [0, T ] × [0, L (i)
] we consider the system [START_REF] Marigo | Entropic solutions for irrigation networks[END_REF] with the condition [START_REF] Osiadacz | Simulation and Analysis of Gas Networks[END_REF] with the matrix A ω as in [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF]. At the end x = L (i) of pipe i (i ∈ I) we apply the controls

r (i) -(t, L (i) ) =        ϑ (i) (t) (t ∈ [0, τ (i) ]) , (1 -ς (i) (t)) ϑ (i) (t) + ς (i) (t) k (i) r (i) + (t -τ (i) (t), L (i) ) (t ∈ (τ (i) , 2τ (i) ]) , k (i) r (i) + (t -τ (i) (t), L (i) ) (t ∈ (2τ (i) , T ]) (16)
with feedback constants k (i) ∈ (-1, 1) and functions

ϑ (i) (t) ∈ C 1 ([0, 2τ (i) ]) and ς (i) (t) ∈ C 1 ([τ (i) , 2τ (i) ]
) that have the following properties:

ς (i) (τ (i) ) = d dt ς (i) (τ (i) ) = d dt ς (i) (2τ (i) ) = 0 and ς (i) (2τ (i) ) = 1 . (17)

Network Lyapunov Function with Delay Terms

In order to define the network Lyapunov function, we define the numbers (i ∈ I)

µ (i) = L (i) 0 1 λ(i) + (x) + 1 | λ(i) -(x)| dx -1 > 0 (18)
and the functions (i ∈ I, x ∈ [0, L (i) ])

h (i) ± (x) = exp -µ (i) x 0 1 λ(i) ± (ξ) dξ > 0 . (19) For constants A (i)
± > 0 we define the weighted and squared L 2 -norms (i ∈ I, t ∈ [0, T ])

E (i) (t) = L (i) 0 A (i) + λ(i) + (x) h (i) + (x) (r (i) + (t, x)) 2 + A (i) - | λ(i) -(x)| h (i) -(x) (r (i) -(t, x)) 2 dx
(20) and the delay terms (i ∈ I, t ∈ [τ (i) , T ])

D (i) (t) = τ (i) (t) 0 A (i) + h (i) + (L (i) ) exp(-µ (i) s) (r (i) + (t -s, L (i) )) 2 ds . (21)
The network Lyapunov function F(t) is defined as (t ∈ [τ max , T ])

F(t) = i∈I E (i) (t) + D (i) (t) . (22)
Constants of the form µ (i) in (18) and functions of the form h (i)

± (x) and E (i) (t) in ( 19) and (20) have been introduced in [START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF]. Delay terms of the form (21) have been presented in [START_REF] Nicaise | Stability of the heat and of the wave equations with boundary time-varying delays[END_REF].

Main Result: Exponential Stability

Theorem 1 states the existence of a unique C 1 -solution of the system (11) on [0, T ] × [0, L (i) ] (i ∈ I) with the boundary conditions ( 13) and ( 16) and initial data of the form (i ∈ I, x ∈ [0, L (i) ]) (r

(i) + (0, x), r (i) -(0, x)) = (ϕ (i) + (x), ϕ (i) -(x)) (23)
with functions ϕ

(i) ± ∈ C 1 ([0, L (i) ]
). For this solution, the function F(t) decays exponentially on [2τ max , T ] (see ( 36)). The decay rate is η = min i∈I α (i) β (i) µ (i) with numbers α (i) ∈ (0, 1) and β (i) ∈ (0, 1). The C 1 -norm of the solution r (i) ± (t, x) is bounded by a constant ε 1 > 0 (see ( 35)) which has to satisfy (i ∈ I)

ε 1 < min x∈[0, L (i) ] | λ(i) ± (x)| , 2ε 1 < min x∈[0, L (i) ] | R(i) + (x) + R(i) -(x)| (24)
and

ε1 θ 4 + 1 2 3 + max exp(1) A (i) - A (i) + λ(i) + (L (i) ) | λ(i) -(L (i) )| , A (i) + A (i) - | λ(i) -(0)| λ(i) + (0) < (1 -β (i) )α (i) µ (i) . (25)
The C 1 -norms of the initial data ϕ (i) ± and the functions ϑ (i) in the controls (16) have to be sufficiently small. More precisely, there exists a number ε 2 ∈ (0, ε 1 ] such that the following inequalities have to hold:

||ϕ (i) ± || C 1 ([0, L (i) ]) ≤ ε 2 (26)
and

||ϑ (i) || C 1 ([0, τ (i) ]) ≤ ε 2 , ||(1 -ς (i) )ϑ (i) || C 1 ([τ (i) , 2τ (i) ]) ≤ ε 2 2 . ( 27
)
Note that the second inequality in (27) holds for any ς (i) if the C 1 -norm of ϑ (i) on [τ (i) , 2τ (i) ] is small enough. For Theorem 1 we define the positive real numbers (i ∈ I)

U (i) ± = max x∈[0, L (i) ] λ(i) ± (x) λ(i) ∓ (x) K (i) ∓ (x) K (i) ± (x) > 0 (28)
and

V (i) ± = min x∈[0, L (i) ] λ(i) ± (x) λ(i) ∓ (x) K (i) ∓ (x) K (i) ± (x) 1 + (1 -α (i) )µ (i) K (i) ∓ (x) > 0 . (29)
Theorem 1. Consider a star-shaped network of pipes as described in Sect. 2.1. Let a finite time T > 0 and functions τ (i) (t) ∈ C 1 ([0, T ]) with ( 14) be given (i ∈ I). Define the constants τ (i) , τ (i) and τ max as in [START_REF] Wang | Exact controllability for nonautonomous first order quasilinear hyperbolic systems[END_REF]. Let real numbers α (i) ∈ (0, 1), β (i) ∈ (0, 1) and a stationary subsonic state ( R(i)

+ (x), R(i) -(x)) ∈ (C 1 ([0, L (i) ]
)) 2 with the eigenvalues λ(i) ± (x) as in (4) be given (i ∈ I) that satisfies the coupling conditions [START_REF] Gugat | Boundary feedback stabilization by time delay for one-dimensional wave equations[END_REF]. Define the numbers σ (i) ∈ {-1, 1} and µ (i) > 0 as in ( 9) and ( 18) and the functions K (i) ± > 0 and h (i) ± > 0 as in [START_REF] Nicaise | Stability of the heat and of the wave equations with boundary time-varying delays[END_REF] and ( 19). Define the numbers U (i) ± > 0 and V (i) ± > 0 as in ( 28) and ( 29) and assume that we have exp( 1)

U (i) + ≤ V (i) + or exp(1) U (i) -≤ V (i) -. (30)
Choose constants A

± > 0 that satisfy

A (i) + ≤ 1 ≤ A (i) - (31)
and assume that we have

A (i) + /A (i) -∈ [exp(1) U (i) + , V (i) + ] or A (i) -/A (i) + ∈ [U (i) -, exp(-1) V (i) -] . ( 32
)
Choose a real number ε 1 > 0 that satisfies ( 24) and ( 25) for all i ∈ I. Then there exists ε 2 ∈ (0, ε 1 ] such that the following statements hold: Choose functions ϑ (i) (t) ∈ C 1 ([0, 2τ (i) ]), ς (i) (t) ∈ C 1 ([τ (i) , 2τ (i) ]) and ϕ (i) ± (x) ∈ C 1 ([0, L (i) ]) (i ∈ I) that satisfy ( 17), ( 26) and ( 27) and such that the C 1compatibility conditions are satisfied at the points (t, x) = (0, 0) and (t, x) = (0, L (i) ) (see Remark 1). Choose constants k (i) ∈ (-1, 1) (i ∈ I) that satisfy exp(1) (k (i) ) 2 A (i)

-≤ A (i) + exp(-µ (i) τ (i) ) (1 -τ (i) ) ( 33
)
and

|k (i) | ≤ ε 2 /(8ε 1 ||ς (i) || C 1 ([τ (i) , 2τ (i) ]) ) . ( 34
)
Then the initial-boundary value problem [START_REF] Marigo | Entropic solutions for irrigation networks[END_REF], ( 13), ( 16), ( 23) has a unique solution (r

(i) + , r (i)
-) ∈ (C 1 ([0, T ] × [0, L (i) ])) 2 that satisfies

||r (i) ± || C 1 ([0, T ]×[0, L (i) ]) ≤ ε 1 . ( 35
)
For this solution define the functions E (i) (t), D (i) (t) and F(t) as in ( 20), ( 21) and ( 22). Then the Lyapunov function F(t) satisfies the following inequality with η = min i∈I α (i) β (i) µ (i) :

F(t) ≤ F(2τ max ) exp (-η(t -2τ max )) for t ∈ [2τ max , T ] . ( 36
)
Remark 1. The C 1 -compatibility conditions guarantee that the initial data (23) and the boundary conditions ( 13) and ( 16) and their first derivatives fit together at the points (t, x) = (0, 0) and (t, x) = (0, L (i) ) (i ∈ I). They can be calculated from ϕ

(i) ± , ϑ (i)
± , A ω and the equations [START_REF] Marigo | Entropic solutions for irrigation networks[END_REF] (see [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF]).

±

  (t, 0) at the node ω we have the equation (t ∈ [0, T ])

⋆⋆ This work was supported by the DFG SPP 1253 and DAAD D/0811409 (Procope 2009/10). The authors like to thank the Institut Henri Poincaré (Paris, France) for providing a very stimulating environment during the "Control of Partial and

Remark 2. The inequalities (31) and the second assumption in (32) hold if, e.g., V (i) -is sufficiently large such that max{1, U (i)

More precisely, if (37) is satisfied, we can first choose the quotient

-}, exp(-1)

Then, without changing the quotient

± > 0 such that (31) holds. The conditions (32) and (37) are satisfied if the number µ (i) > 0 from (18) is sufficiently large which is the case if the length L (i) is not too long.

Proof of the Main Result Stated in Theorem 1

In this section we prove Theorem 1. The existence of a unique solution of [START_REF] Marigo | Entropic solutions for irrigation networks[END_REF] follows from Theorem 2.1 in [START_REF] Wang | Exact controllability for nonautonomous first order quasilinear hyperbolic systems[END_REF] where initial-boundary value problems for first order quasilinear hyperbolic systems are studied (see also [START_REF] Dick | Classical solutions and feedback stabilization for the gas flow in a sequence of pipes[END_REF][START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF]). For the proof of (36) we use the estimates (38), ( 40) and (41) for d dt E (i) (t), d dt D (i) (t) and d dt F(t). The estimate (38) is the same as in [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF] where 40) is more complicated than in [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF] where only constant delays are considered. Using integration by parts, Young's Inequality and the conditions (25) and (32), we obtain the following estimate for d dt

For a detailed calculation of (38) see [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF]. For the derivative d dt D (i) (t) we get (t ∈ [τ (i) , T ])

Using the equation ∂ s r

) and integration by parts, from (39) we obtain (t ∈ [τ (i) , T ])

Hence, the inequalities ( 14) and the definition of τ (i) and τ (i) in ( 15) imply

From ( 38) and (40) we obtain the following estimate for d dt F(t) with η = min i∈I α (i) β (i) µ (i) (t ∈ [τ max , T ]):

with the boundary terms

The equation ( 13), the orthogonality of the matrix A ω from [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF]


Janko Petereit
email: janko.petereit@iosb.fraunhofer.de

Thomas Bernard
email: thomas.bernard@iosb.fraunhofer.de

Real-Time Nonlinear Model Predictive Control of a Glass Forming Process Using a Finite Element Model

The control of complex forming processes (e.g., glass forming processes) is a challenging topic due to the mostly strongly nonlinear behavior and the spatially distributed nature of the process. In this paper a new approach for the real-time control of a spatially distributed temperature profile of an industrial glass forming process is presented. As the temperature in the forming zone cannot be measured directly, it is estimated by the numerical solution of the partial differential equation for heat transfer by a finite element scheme. The numerical solution of the optimization problem is performed by the solver HQP (Huge Quadratic Programming). In order to meet real-time requirements, in each sampling interval the full finite element discretization of the temperature profile is reduced considerably by a spline approximation. Results of the NMPC concept are compared with conventional PI control results. It is shown that NMPC stabilizes the temperature of the forming zone much better than PI control. The proposed NMPC scheme is robust against model mismatch of the disturbance model. Furthermore, the allowed parameter settings for a real-time application (i.e., control horizon, sampling period) have been determined. The approach can easily be adapted to other forming processes where the temperature profile shall be controlled.

Introduction

The control of complex forming processes (e.g., glass forming processes) is a challenging topic due to the mostly strongly nonlinear behavior and the spatially distributed nature of the process. The nonlinearity is caused on the one hand by the physics of the process (e.g., radiation) and on the other hand by nonlinear material properties (e.g., memory effect). In recent years considerable progress has been made in the field of simulation tools (e.g., finite element based simulation of complex rheological processes [START_REF] Krause | Mathematical Simulation in Glass Technology[END_REF], [START_REF]Glass Days 2005 -Analysis and Simulation of Processes in Glass Production and Processing[END_REF]) and model order reduction (MOR) methods [START_REF] Benner | Dimension Reduction of Large-Scale Systems[END_REF], [START_REF] Antoulas | Approximation of Large-Scale Dynamical Systems[END_REF]. Consequently, the use of these (reduced) models for process control and optimization is in the focus of actual research activities [START_REF]EU project PROMATCH (PROmoting and structuring Multidisciplinary Academic-industrial collaboration in research and Training trough SME teCHnology developers[END_REF]. Nonlinear Model Predictive Control (NMPC) concepts play an important role in this context, as NMPC can be applied for nonlinear and spatially distributed dynamic systems. Main advantages of NMPC approaches are that the Fig. 1. Glass forming process overview performance criteria can be designed in a transparent way and that even time dependent constraints can be applied to manipulated, output and state variables of the controlled system. Drawbacks of NMPC are that computational costs for the solution of the nonlinear programming problems often are high, hence powerful optimization solvers are needed [START_REF] Diehl | An efficient algorithm for nonlinear model predictive control of large-scale systems part 1[END_REF], [START_REF] Diehl | An efficient algorithm for nonlinear model predictive control of large-scale systems part 2[END_REF]. A suited solver for large scale dynamic optimization problems is HQP (Huge Quadratic Programming, [START_REF] Franke | The solver Omuses/HQP for structured largescale constrained optimization: algorithm, implementation and example application[END_REF]). HQP has been used e.g., for NMPC appplications in process industry [START_REF] Nagy | Real-time Implementation of Nonlinear Model Predictive Control of Batch Processes in an Industrial Framework[END_REF], power plants [START_REF] D'amato | Industrial application of a model predictive control solution for power plant startups[END_REF] as well as for optimal management of water resources [START_REF] Bernard | Optimal management of regional water supply systems using a reduced finite-element groundwater model[END_REF].

For complex forming processes in many cases an optimal spatially distributed (and sometimes time-dependent) viscosity profile of the material has to be assured. As viscosity is difficult to measure on-line, temperature often is used as an auxiliary controlled variable. Hence the problem consists in controlling a spatially distributed temperature profile. In [START_REF] Bernard | Nonlinear model predictive control of a glass forming process based on a finite element model[END_REF] a first NMPC concept has been investigated. In this work a predictive functional control scheme (PFC) has been applied in order to optimize the forming control loops, but not the temperature control. Only a very short prediction horizon is assumed in the mentioned simulation studyy. In the present paper a new approach for the real-time control of a spatially distributed temperature profile is presented. As the temperature in the forming zone cannot be measured directly, it is estimated by the numerical solution of the partial differential equation (PDE) for heat transfer by a finite element (FE) scheme. As the dimension of the state space model, which is result of the FE algorithm, is too large for real-time optimization, in each sampling interval the full finite element discretization of the temperature profile is reduced considerably by a spline approximation. The numerical solution of the optimization problem is performed by HQP. The approach can easily be adapted to other forming processes where the temperature profile shall be controlled.

The paper is organized as follows. Section 2 describes an industrial glass forming process to which the NMPC approach is applied. The details of the proposed NMPC concept are presented in section 3. In section 4 the results of the NMPC concept are discussed. Section 5 gives some conclusions.

2 Optimal temperature profile of an industrial glass forming process

The NMPC concept which will be introduced in the next section has been applied to an industrial glass forming process (Fig. 1). A vertically hung glass cylinder is fed into a ring-shaped furnace where it begins to liquefy. It then starts flowing downwards viscously and is taken up by pullers which pull the resulting glass tube with an appropriate speed. It is essential for the process to stabilize the viscosity profile along the cylinder. The viscosity can be controlled by the furnace temperature.

Figure 1 depicts the relevant sensors and actuators regarding the temperature control: As the furnace tightly encloses the glass cylinder, there are only three temperature sensors available -one each right above (T cyl ) and below (T tube ) the furnace, and one which measures the furnace temperature (T furnace ) itself. There is no possibility to measure the temperature in the forming zone (T form ) directly. Thus, a strategy is needed to reconstruct this control variable using the measurements of the only three available sensors. This is achieved by a finite element model. The forming temperature is controlled by increasing or decreasing the heating power and hence furnace temperature.

In order to ensure optimal product quality there are three major variables which must be controlled in an optimal way. The first two are the feeding and pulling speed (respectively force) which are accounted for by a common control loop. We call this one the geometry control. The third variable is the (spatially distributed) temperature of the glass in the forming zone, whose control is done by a seperate loop, which we call the temperature control. In the scope of this paper we will concentrate on the latter, i.e., the temperature control.

For glass the system's behavior is in general strongly nonlinear due to the impact of radiation. The spatio-temporal temperature distribiution can be calculated solving the heat transfer PDE. If a symmetric cylinder is assumed and radial temperature distribution the following one dimensional PDE describes the temperature distribution T (z, t) along the vertical z-axis:

ρc p (T ) ∂T ∂t = λ(T ) ∂ 2 T ∂z 2 + qconv (T, z) + qrad,oven (T, z) + qrad,dist (z, t) (1)
Radiative heat transfer exchange between the oven and the glass is described by the nonlinear term qrad,oven (Stefan-Boltzmann law). Radiation inside the glass is considered by means of an effective heat conduction coefficient λ(T ) which is nearly exponential increasing with temperature and hence introduces a second nonlinearity in the model. The term qconv describes the convective heat transfer due to the movement of the glass cylinder. A further nonlinearity is introduced by temperature dependence of the specific heat capacity c p (T ).

At the end of a batch production (i.e., the when the length of the cylinder which is fed on the oven tends to zero), disturbances arise due to radiation effects which are caused by the end of the cylinder. This radiative disturbance are described by the space and time dependent term qrad,dist (z, t). As the disturbance cannot be measured directly, it is estimated online by a special parameter estimation approach.

Due to the transport length and the corresponding dead time it is important that the disturbances are suppressed as early as possible by an optimal strategy of furnace temperature during process end phase. As the steadiness of the viscosity in the forming zone exerts the main impact upon the quality of the final product, an optimal control must stabilize the spatially distributed temperature (and hence viscosity) profile of the forming zone as good as possible.

To accommodate for the deficiencies of the PI controller we developed a control concept which explicitly takes into account the large dead time of the process and the disturbances arising during the end phase. As a PI controller can only react to changes in the control variable which are already in effect, we switched over to a predictive control strategy, which is able to take corrective action before any impact of the disturbances is visible to the temperature sensors. This behavior is necessary since the delay time of the furnace is too large to lower the temperature in sufficient time. Otherwise, overshooting of the forming temperature would be unavoidable.

As the prediction model for the NMPC scheme we exploit the finite element model of the 1D temperature distribution T (z, t) along the vertical z-axis by solving the PDE (1) online. The finite element model considers heat conduction, convection by the movement of the glass, radiation and the time dependent disturbances caused by the end of the cylinder.

Concept for real-time Optimization of a spatially distributed temperature profile

In our solution concept (Fig. 2) we propose a real-time nonlinear model predictive control of the temperature of a forming zone (here for glass forming processes). In most cases the temperature profile cannot be measured directly, hence it has to be estimated by a spatially distributed model. The one dimensional heat transfer PDE (1) has to be solved numerically for the temperature distribution T (z, t), where z denotes the height and t the time. The time-dependent disturbances at the end of a cylinder are estimated by a disturbance model using the knowledge from previous productions to determine the model's parameters in a predictive way.

The NMPC approach is formulated as follows. The predicted values of the controlled variable y(t) are collected in the vector ŷ for the timesteps from t k ...t k+np , where t k is the actual time and n p the prediction horizon. The predicted values are calculated by the model where x k denotes the actual state of the model and the vector ∆u denotes all changes of the controlled variables from the actual time step t k to the prediction horizon (we set prediction horizon n p equal control horizon n p ). The optimization

ŷ = F (x k , ∆u) (2)
J = (w -ŷ) ⊤ Q (w -ŷ) + ∆u ⊤ R ∆u (3)
∆u opt = argmin ∆u J (4)
is performed subject to

∆u min ≤ ∆u(k + j|k) ≤ ∆u max j = 0, . . . , n C -1 ŷmin ≤ ŷ(k + j|k) ≤ ŷmax j = 1, . . . , n P (5)
In this formulation the changes of the manipulated variable u are calculated. If necessary, also the absolute value u could be considered directly in the optimization problem. The proposed control concept has been implemented in a C++ application by utilizing the software suite HQP (Huge Quadratic Programming), a solver for large scale nonlinear programming problems, to solve the arising nonlinear optimal control problems. In order to meet real-time requirements, in each sampling interval the full finite element discretization of the temperature profile is reduced considerably by a spline approximation (see Fig. 4 and3). The spline approximation is admissible due to the temperature profile's smoothness. In our application the number of states could be reduced from 200 to 30.

Results

We evaluated the control performance resulting with the proposed NMPC concept and compared it to a PI controller using previously recorded production data. It is shown that NMPC stabilizes the temperature of the forming zone much better than PI control. As can be seen in Fig. 5 with very small prediction and control horizon (n p = n c = 2) the performance of NMPC and PI control do not differ very much. In subplot 1 of Fig. 5 the deviation e forming from the reference value is shown for a production time of 40 minutes. In subplot 2 the change of the furnace temperature (manipulated variable) compared to the value at starting time t = 0 min. can be seen (i.e., Tfurnace (t) = T furnace (t) -T furnace (0)). Subplot 3 shows the changes of furnace temperature in each sample time. These changes have been restricted to ±5 K during the optimization in order to guarantee a smooth change of the furnace temperature. Finally, in subplot 4 the computation time for each sample interval (30 s) is shown. It can be seen that the computation time does not differ very much (about 2 s). Fig. 6 shows the results with a much larger prediction and control horizon (n p = n c = 30). From subplot 1 it can be observerd that the control deviation with NMPC is much smaller compared to the results with PI control. Accordingly, the decrease of the furnace temperature (subplot 2) with NMPC starts much earlier. This is caused by the internal disturbance model which predicts the disturbances early enough. From subplot 3 it can be seen that the defined threshold for changes of the furnace temperature are reached for a period of about 10 minutes. In subplot 4 it can be seen that the computation is always smaller than the sample time of 30 s, hence the real-time condition is satisfied. The computation time decreases strongly in the last 15 minutes of the production due to the shortened horizon.

The allowed parameter settings for a real-time application (e.g., control horizon, sampling period) have been determined. In Fig. 7 the computation time is plotted versus the prediction and control horizon (n p , n c ). The internal time discretization dτ of the finite element model (n FDM = T s /dτ ) is used as parameter (T s : sample time). From Fig. 7 it is obvious that for a wide range of parameter values n p , n c and n FDM real-time application (computation time < sample time = 30 s) is possible.

Furthermore, simulation based application of the proposed NMPC scheme to a large number of historical runs showed that the concept is robust against model mismatch of the disturbance model. This property is very important for the industrial implementation.

Conclusion

In this paper a new approach for the real-time control of a spatially distributed temperature profile of an industrial glass forming process has been presented. A numerically efficient solution of the partial differential equation has been implemented. By means of the powerful solver HQP a real-time application could be achieved. It has been shown that with the proposed NMPC concept much better disturbance rejection and hence stabilization of the temperature profile than conventional PI control is achieved. The approach can easily be adapted to other forming processes where the temperature profile shall be controlled.

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. NMPC temperature control scheme

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Comparison NMPC and PI control (nP = nC = 2)


L Gerard Van Willigenburg
email: gerard.vanwilligenburg@wur.nl

Willem L De Koning

L G Van Willigenburg

Keywords: Temporal system properties, linear discrete-time systems, cheap LQ control problems, j-step k-step Kalman decomposition

Time-varying discrete-time linear systems may be temporarily uncontrollable and unreconstructable. This is vital knowledge to both control engineers and system scientists. Describing and detecting the temporal loss of controllability and reconstructability requires considering discretetime systems with variable dimensions and the j-step, k-step Kalman decomposition. In this note for linear discrete-time systems with variable dimensions measures of temporal and one-step stabilizability and detectability are developed. These measures indicate to what extent the temporal loss of controllability and reconstructability may lead to temporal instability of the closed loop system when designing a static state or dynamic output feedback controller. The measures are calculated by solving specific linear quadratic cheap control problems.

Introduction

Feedback control design and stability analysis of nonlinear systems along trajectories is often performed using the linearized dynamics about the trajectory [START_REF] Athans | The role and use of the Linear-Quadratic-Gaussian problem in control system design[END_REF], [START_REF] Van Willigenburg | On the synthesis of time-varying LQG weights and noises along optimal control and state trajectories[END_REF]. If the trajectory is time-varying the linearized model is time-varying. If in addition the nonlinear dynamics or the controls are non-smooth, i.e. in the case of bang-bang or digital control, the structure of the time-varying linearized system may change. Even if the nonlinear dynamics and the controls are smooth the structure of the time-varying linearized system may almost change. For control system design this is vital information since this structure reveals the temporal loss of controllability and reconstructability of the linearized system. They in turn may lead to temporal instability of a closed-loop control system [START_REF] Van Willigenburg | A Kalman decomposition to detect temporal linear system structure[END_REF], [START_REF] Van Willigenburg | Temporal linear system structure[END_REF]. Recently we investigated these issues for continuous-time systems assuming continuous-time control . This investigation lead to the introduction of the properties temporal and differential stabilizability and detectability for continuoustime linear systems [START_REF] Van Willigenburg | Temporal and differential stabilizability and detectability of piecewise constant rank systems[END_REF]. In addition measures of these properties were introduced and calculated by solving specific linear quadratic cheap control problems [START_REF] Van Willigenburg | Temporal and differential stabilizability and detectability of piecewise constant rank systems[END_REF], [START_REF] Jameson | Cheap control of the time-invariant regulator[END_REF], [START_REF] Kokotovic | Singular perturbations and order-reduction in control theory -an overview[END_REF].

Associated with computer control are digital control problems (sampled-data control problems). They concern the control of continuous-time systems by means of piecewise constant controls using sampled measurements. A common approach is to transform such control problems into equivalent discrete-time control problems [START_REF] Levis | On the behavior of optimal linear sampled-data regulators[END_REF], [START_REF] Van Willigenburg | The digital optimal regulator and tracker for stochastic time-varying systems[END_REF], [START_REF] Van Willigenburg | Digital optimal control and LQG compensation of asynchronous and aperiodically sampled nonlinear systems[END_REF]. Following this approach feedback control system design is performed in discrete-time. This motivates the discrete-time development in this paper that on the one hand parallels, but on the other is also very different from the one in continuous-time. The fact that discrete-time is not dense, as opposed to continuous-time, causes some major differences. In continuous-time our investigation required the introduction of piecewise constant rank systems and the differential Kalman decomposition [START_REF] Van Willigenburg | A Kalman decomposition to detect temporal linear system structure[END_REF], [START_REF] Van Willigenburg | Temporal linear system structure[END_REF]. In discrete-time their counterparts are discrete-time linear systems with variable state dimensions and the j-step, k-step Kalman decomposition [START_REF] Van Willigenburg | Temporal linear system structure: The discrete-time case[END_REF].

This paper develops measures of temporal stability of time-varying linear discrete-time systems over arbitrary finite time intervals, notably intervals where controllability or reconstructability is lost temporarily. Associated to this, measures of temporal and one-step stabilizability and detectability are developed. These measures can for instance be used to analyse temporal instability of a closed loop control system design using LQG output feedback.

Temporal stability may sound as a contradiction because formally stability relates to behavior when time tends to infinity. However, in one of his early seminal papers [START_REF] Kalman | Control system design via the "Second Method" of Lyapunov, I Continuous-time systems[END_REF] Kalman together with Bertram already proposed measures of stability over finite time intervals (page 386). Intuitively stability relates to growth of the system state. Intuitively over intervals where the state grows we call the system temporal unstable and over intervals where the state decays, we call the system temporal stable. This intuition is formalized by the temporal stability property proposed in this note. This property is derived from a measure of temporal stability also proposed in this note that measures the maximum growth of the state over an arbitrary interval. Our concept of stability over a finite time interval differs from what is called finite-time stability [START_REF] Amato | Finite-Time Control of Linear Systems: A Survey[END_REF], [START_REF] Amato | Finite-time stability of linear timevarying systems with jumps[END_REF]. The reason we make a different choice is that our measures, their computation and the associated control system designs, come down to solving standard LQ problems. The standard LQ problems are of a special type called cheap control LQ problems [START_REF] Jameson | Cheap control of the time-invariant regulator[END_REF], [START_REF] Kokotovic | Singular perturbations and order-reduction in control theory -an overview[END_REF] . They are characterized by a control penalty that tends to zero. Computations and control system design associated to finite-time stability concern matrix inequalities [START_REF] Amato | Finite-Time Control of Linear Systems: A Survey[END_REF], [START_REF] Amato | Finite-time stability of linear timevarying systems with jumps[END_REF]. Generally these are much more difficult to solve.

Temporal and one-step stability, stabilizability and detectability

Temporal uncontrollability/unreachability and temporal unreconstructability/ unobservability of linear time-varying systems was introduced and investigated in continuous-time [START_REF] Van Willigenburg | A Kalman decomposition to detect temporal linear system structure[END_REF], [START_REF] Van Willigenburg | Temporal linear system structure[END_REF] and in discrete-time [START_REF] Van Willigenburg | Temporal linear system structure: The discrete-time case[END_REF]. Intuitively, temporal stabilizability and temporal detectability are associated properties that apply over intervals where the system is temporal uncontrollable/unreachable and temporal unreconstructable/unobservable respectively. In continuous-time this was formalized in [START_REF] Van Willigenburg | Temporal and differential stabilizability and detectability of piecewise constant rank systems[END_REF]. In this section we formalize the discrete-time case. This requires considering variable dimension discrete-time linear systems (VDD systems) [START_REF] Van Willigenburg | Temporal linear system structure: The discrete-time case[END_REF], [START_REF] Gohberg | Minimality and realization of discrete time-varying systems[END_REF], [START_REF] Van Willigenburg | Minimal and non-minimal optimal fixed-order compensators for time-varying discrete-time systems[END_REF], [START_REF] Sandberg | Balanced truncation of linear time-varying systems[END_REF] as well as j-step controllability, j-step reachability, k-step reconstructability, k-step observability and the associated j-step, k-step Kalman decomposition. All these are introduced in [START_REF] Van Willigenburg | Temporal linear system structure: The discrete-time case[END_REF] that relies partly on [START_REF] Van Willigenburg | Linear systems theory revisited[END_REF]. In this section we consider VDD systems with a time domain [ ] 0 , N i i where 0 i may tend to -∞ and N i may tend to +∞ . Intervals where the VDD system is temporal uncontrollable/unreachable or temporal unreconstructable/unobservable are denoted by ,

s f i i     .
Definition 1

A VVD system is called j-step unreachable over the interval ,

s f i i j   +   / j-step uncontrollable over the interval , s f i j i   -   , 0 s f N i j i i i j + ≤ < ≤ -if , s f i i i j   ∀ ∈ + 
 the system is not j-step reachable at time i / not j-step controllable from time i j -.

Lemma 1

If , s f i i satisfy the conditions in Definition 1 then over the interval ,

s f i i   
 the VDD system is 1) not j- step reachable at each time and 2) not j-step controllable from each time.

Proof

Follows immediately from [START_REF] Van Willigenburg | Temporal linear system structure: The discrete-time case[END_REF] and Definition 1.

Definition 2

A VDD system that satisfies the conditions in Definition 1 is called j-step uncontrollable/unreachable over the interval ,

s f i i     .
Definition 3 (Dual of Definition 1)

A VVD system is called k-step unobservable over the interval ,

s f i k i   -   / k-step unreconstructable over the interval , s f i i k   +   , 0 s f N i k i i i k + ≤ < ≤ -if , s f i i k i   ∀ ∈ -   the system is not k-step observable at time i / not k-step reconstructable from time i k + .
Lemma 2 (dual of Lemma 1)

If , s f
i i satisfy the conditions in Definition 3 then over the interval ,

s f i i   
 the VDD system is 1) not k-step observable at each time and 2) not k-step reconstructable from each time.

Definition 4 (dual of Definition 2)

A VDD system that satisfies the conditions in Definition 3 is called k-step unreconstructable/unobservable over the interval ,

s f i i     .
Application of the j-step k-step Kalman decomposition [START_REF] Van Willigenburg | Temporal linear system structure: The discrete-time case[END_REF], [START_REF] Boley | Computing the Kalman decomposition: An optimal method[END_REF] at each time

[ ] 0 , N i i i ∈
, reveals all closed intervals (i.e. consisting of at least two consecutive discrete-time instants) where the system is jstep uncontrollable/unreachable and dually all closed intervals where the system is k-step unreconstructable/unobservable. As in Definition 2 and Definition 4 such intervals will be denoted by ,

s f i i   
 . These closed intervals are precisely the intervals where stability of the closed loop system may be lost temporarily when designing static state and dynamic output feedback controllers.

Stabilizability is a property that relates entirely to the uncontrollable part of a system. A general approach to determine stabilizability is to extract this uncontrollable part, that is autonomous, by means of a Kalman decomposition, and to determine its stability. It will become clear in this section that application of a state basis transformation changes temporal stability and stabilizability properties. To recover them we therefore need to transform back to the original state basis. As opposed to this general approach, the stabilizability analysis presented in this section is much more straightforward and simple. It does not require transformation of the state basis because it relies fully on well established standard LQ theory applied to the original system representation. Therefore the associated numerical computations are also very efficient.

The stabilizability analysis in this section is unconventional in the sense that stability, stabilizability and detectability over finite time intervals is required. Stability over an interval relates to growth of the magnitude of the state over this interval. Throughout this paper i denotes the matrix 2 norm. For vectors this amounts to the L2 norm. In the next section we will demonstrate how to compute numerically the temporal and one-step stabilizability and detectability measures presented in this section, using only evaluations of the system matrices.

Definition 5

An autonomous VVD system is called temporal stable over the interval ,

s f i i     if for any 0 s i x ≠ , 1 f s i i x x < .
Loosely speaking, according to Definition 5 an autonomous VDD system is called temporal stable over ,

s f i i   
 if for any initial state the magnitude of the associated terminal state is smaller than that of the initial state. An important difference between our definition and other finite-time stability concepts [START_REF] Amato | Finite-Time Control of Linear Systems: A Survey[END_REF], [START_REF] Amato | Finite-time stability of linear timevarying systems with jumps[END_REF] is that ours does not impose any restrictions on the magnitude of the state inside the interval. The advantage of Definition 5 is that it matches LQ control design as opposed to finite-time stability that relates to control system design using matrix inequalities [START_REF] Amato | Finite-Time Control of Linear Systems: A Survey[END_REF] that is generally much more complicated.

Definition 6

Associate to Definition 5 the following temporal stability measure, ( )

2 2 0 , max 0 f i s s i s f x i x i i x ρ ≠     = ≥       . ( 1
)
Observe that ( ) 1) are needed to achieve compatibility with LQ control computations.

Theorem 1

An autonomous VDD system is temporal stable over the time interval ,

s f i i     if and only if, ( ) , , , 1
s f s f T s f i i i i i i ρ = Φ Φ < , (2) where
,

s f i i
Φ represents the state transition matrix of the associated autonomous system from time s i to f i .

Proof

Because Theorem 1 applies to autonomous systems, ,

f s f s i i i i x x = Φ . ( 3
)
Using equation ( 3) the temporal stability measure (1) becomes, ( )

. s f s i s s s s f s f s s f s f i s s s i i i s f x i T T i i i i i i T i i i i T x i i x i i x x x x x ρ ≠ ≠   Φ   = =         Φ Φ   = Φ Φ     2 , 2 0 , , , , 0 , max max
The last equality in equation ( 4) holds because , ,

s f s f T i i i i
Φ Φ is nonnegative symmetric. Theorem 1 now follows from (4), Definition 5 and Definition 6 and,

2 2 1 1 f s f s i i i i x x x x < ⇔ < . (5)
Stabilizability over a finite time-interval relates to the ability to stabilize the system over that interval by means of control.

Definition 7

Associate to Definition 5 and Definition 6 the following temporal stabilizability measure that applies to VDD systems considered over the interval ,

s f i i , ( ) 2 | in 2 0 m min max , 0 f i i s i s s i u x s f x i x i i x ρ ≠     ≥       = , (6)
where | s i i u x indicates a control law dependent on s i

x .

Definition 8

A VDD system is called temporal stabilizable over ,

s f i i     if ( ) min , 1
s f i i ρ < .
Theorem 2

A VDD system is temporal controllable over ,

s f i i     ⇒ ( ) min , 0 s f i i ρ = ⇒ the VDD system is temporal stabilizable over , s f i i     .

Proof

If a VDD system is temporal controllable over ,

s f i i     ,
s f i i     .

Remark 1

As with ordinary controllability and stabilizability, temporal controllability is a stronger property than temporal stabilizability.

To state the main theorem in this section consider the following parameterized discrete-time LQ problem. Given the system, 1 , , 1

i i i i i s f x x u i i i +   = Φ + Γ ∈ -   , (7)
with initial state,

s i x , (8)
find the control , , 1

i s f u i i i   ∈ - 
 that minimizes the cost function, ( )

1 f f f s i T T T LQ i i i i i i i i i i J x Hx x Q x u R u ε ε - =   = + +   ∑ , (9)
with,

n H I = , 0 i Q = , , 0 1 i m R I ε ε ε = ≤ ≪ . ( 10
)
If 0 ε > the Linear Quadratic control problem ( 7), ( 8)-( 10) satisfies 0, 0, 0

i i H Q R ε ≥ ≥ > .
In this standard case it is well known that the optimal control is given by, ( )

1 1 1 , T T i i i i i i i i i i i u L x L S R S ε ε ε ε ε - + + = - = Γ Γ + Γ Φ , ( 11
)
and the minimum cost by, ( )

* s s s T LQ i i i J x S x ε ε = , ( 12
)
where , ,

i s f S i i i ε   ∈ 
 is the solution of the matrix Riccati difference equation, ( )

1 1 , f T T T i i i i i i i i i i i i S S L S R L Q S H ε ε ε ε ε + + = Φ Φ - Γ Γ + + = . ( 13
)
Theorem 3 * 0 lim i i S S ε ε ↓ = , (14)
exists, where , 0 i S ε ε ↓ satisfies the matrix Riccati difference equation ( 13) with data as specified by equation [START_REF] Van Willigenburg | Digital optimal control and LQG compensation of asynchronous and aperiodically sampled nonlinear systems[END_REF]. Furthermore,

( ) * min , , , 1
f i s f i i S i i i ρ   = ∈ -   . ( 15
)
As a special case of ( 15),

( ) * min , s s f i i i S ρ = . ( 16
)

Proof

First observe that in the parameterized LQ problem ( 7)-( 10) we may replace the initial time s i by , 1

s f i i i ′   ∈ - 
 . This also hold for the stabilizability measure min ρ . Next from equations ( 9), [START_REF] Van Willigenburg | Digital optimal control and LQG compensation of asynchronous and aperiodically sampled nonlinear systems[END_REF] observe that ( ) ( )

| | | 2 min 0 min min i i i i f f i i f LQ u x u x u x T i i i x x x J ′ ′ ′ = = (17)
Now the key to proving ( 14), ( 15) is to prove that,

( ) ( ) | * * 0 0 lim min i i u x T LQ LQ i i i J J x S x ε ε ′ ′ ′ ′ ↓ = = (18)
Suppose equation ( 18) holds. Then from equations ( 6), ( 17), (

f i i i i i u x f x i T i i i i T x i i x i i x x S x S x x ρ ′ ′ ′ ≠ ′ ′ ′ ′ ′ ′ ′     ′ = =         =     (19)
The last equality in equation ( 19) holds because * i S ′ is nonnegative symmetric. So we are left to prove equation [START_REF] Van Willigenburg | Linear systems theory revisited[END_REF]. Consider the j-step, k-step Kalman decomposition at time f i with

f j i i′ = -. According
to this decomposition the linear system (7) can be decomposed into a part that is j-step controllable from time i′ and a part that is autonomous. The contribution of the j-step controllable part to J ε is fixed and independent of ε . Because the system matrices are bounded this contribution is also finite. This proves the existence of the limit ( 16) and the equality (18).

Remark 2

There are three reasons for considering 0 1 ε < ≪ in equation [START_REF] Van Willigenburg | Digital optimal control and LQG compensation of asynchronous and aperiodically sampled nonlinear systems[END_REF], instead of 0

ε = . Taking 0 1 ε < ≪
, ε may be used to 1) keep the control within certain bounds that apply in practice and 2) as a numerical tolerance to prevent ill-conditioning of the computation of equation ( 11) when

1 T i i i S ε + Γ Γ is not full rank and i L → ∞ as 0 ε ↓ .
In practice the selection of 0 1 ε < ≪ will be a compromise and i S ε will

approximate * i S , , 1
s f i i i   ∈ -   .
As a result all computations in this paper involving i S ε will be approximations, although generally very good ones. Thirdly 0 ε = leads to a singular LQ problem that is generally much more difficult to solve and the solution of which need not be unique.

When analyzing control systems the state behavior over the entire interval ,

s f i i   
 is generally of interest, not just the behavior at the initial time s i and the final time f i . This behavior is partly considered by equation [START_REF] Gohberg | Minimality and realization of discrete time-varying systems[END_REF] of Theorem 3 that determines the stabilizability measure for each sub interval , , , 1

f s f i i i i i     ∈ -     .
The following theorem introduces a one-step stabilizability measure that applies to individual time instants.

Theorem 4

* * 1 i i S S + -
is a one-step stabilizability measure (os-stabilizability measure) at time , 1

s f i i i   ∈ -   .

Proof

From ( 15), ( ) ( )

1 * * * * min 1 , f s f s i s f i i i i i i i i S S S S ρ - + = = = + - ∑ (20) so
* * 1 i i S S + - , , 1
s f i i i   ∈ -  

Definition 9

A VDD system is called one-step stabilizable (os-stabilizable) at time , 1

s f i i i   ∈ -   if * * 1 0 i i S S + - < .
Because for a VDD system temporal and one-step detectability are dual to temporal and one-step stabilizability, the following definitions and theorems are stated without further explanation and proof.

Theorem 5 (dual of Theorem 3) * 0 lim

i i P P ε ε ↓ = , (21)
exists, where , 0 i P ε ε ↓ satisfies the matrix Riccati difference equation that is dual to (13), ( )

1 , s T T T i i i i i i i i i i i i P P L C P C R L Q P H ε ε ε ε ε ε + = Φ Φ - + + = , (22) with,
( )

1 T T i i i i i i i i L P C C P C R ε ε ε ε - = Φ + , (23)
with data as specified by equation [START_REF] Van Willigenburg | Digital optimal control and LQG compensation of asynchronous and aperiodically sampled nonlinear systems[END_REF]. Furthermore, (

i s f i i P i i i σ   = ∈ +   . ) * min , , 1, s
where

( ) min , f i i σ is a temporal detectability measure over the interval , f i i     . As a special case, ( ) * min , f s f i i i P σ = . ( 25
)
Definition 10 (dual of Definition 8)

A VDD system is called temporal detectable over ,

s f i i     if ( ) min , 1 s f i i σ < . Theorem 6 (dual of Theorem 4) * * 1 i i P P + -
is a one-step detectability measure (os-detectability measure) at time , 1

s f i i i   ∈ -   .
Definition 11 (dual of Definition 9)

A VDD system is called one-step detectable (os-detectable) at time , 1

s f i i i   ∈ -   if * * 1 0 i i P P + - < .

Conclusions

New temporal properties and associated measures for control system design concerning time-varying linear discrete-time systems were introduced in this paper. The properties and associated measures concern temporal and one-step stabilizability and detectability. They indicate to what extent control system design is problematic when discrete-time linear time-varying systems are temporal uncontrollable or temporal unreconstructable. Temporal uncontrollability and unreconstructability are detected by the j-step, k-step Kalman decomposition. As demonstrated in this paper, after introduction of a suitable, simple stability property, that applies over finite time intervals, application of ordinary standard LQ theory and algorithms enables the computation of associated temporal and one-step stabilizability and detectability measures. These determine to what extent a static or dynamic feedback control system becomes temporal unstable. A major application concerns the temporal stability analysis of digital perturbation output feedback controllers for nonlinear systems tracking control and state trajectories that may be optimal [START_REF] Athans | The role and use of the Linear-Quadratic-Gaussian problem in control system design[END_REF], [START_REF] Van Willigenburg | Digital optimal control and LQG compensation of asynchronous and aperiodically sampled nonlinear systems[END_REF].

As an alternative to LQ theory, temporal stabilizability may be determined by extracting the temporal uncontrollable or temporal unreconstructable subsystems and analyzing their temporal stability. In principle, the j-step, k-step Kalman decomposition is able to extract these subsystems. The extraction employs state basis transformations that generally change temporal stability properties. The approach presented in this paper is more simple and direct because it applies standard LQ theory to the original, untransformed system.

Although the LQ problems in this paper are singular in principle, it is advantageous to approximate them by non-singular LQ problems, as demonstrated in this paper. The interpretation of * i S as a temporal stabilizability measure is new and highly interesting. The same applies to the interpretation of * * 1 i i S S + as a one-step stabilizability measure that measures the contribution to stabilizability of each single time-step.

Along the lines of this paper we are also currently exploring temporal properties of time-varying linear systems with white stochastic parameters [START_REF] Van Willigenburg | Compensatability and optimal compensation of systems with white parameters in the delta domain[END_REF]. Among others these enable robust digital optimal perturbation feedback design for nonlinear systems.

  stability. It will become clear that the squares in equation (

  then according to Definition 1 and[START_REF] Van Willigenburg | Temporal linear system structure: The discrete-time case[END_REF], any state

  is the one-step contribution at time i to the temporal stabilizability

	measure	ρ	min	( ) , s f i i	. If this contribution is negative	ρ	min	( ) , s f i i	decreases and temporal stabilizability
	increases.