Minh Ngoc Dao
email: minhdn@hnue.edu.vn

Dominikus Noll
email: noll@mip.ups-tlse.fr

Minimizing memory effects of a system

Keywords: System ringing, system memory, Hankel norm, system reduction, controller design, system with tunable parameters

Given a stable linear time-invariant system with tunable parameters, we present a method to tune these parameters in such a way that undesirable responses of the system to past excitations, known as system ringing, are avoided or reduced. This problem is addressed by minimizing the Hankel norm of the system, which quantifies the influence of past inputs on future outputs. We indicate by way of examples that minimizing the Hankel norm has a wide scope for possible applications. We show that the Hankel norm minimization program may be cast as an eigenvalue optimization problem, which we solve by a nonsmooth bundle algorithm with a local convergence certificate. Numerical experiments are used to demonstrate the efficiency of our approach.

Introduction

Ringing generally designates undesired responses of a system to past excitations. In electronic systems, ringing arises under various forms of noise, such as gate ringing in converters, undesired oscillations in digital controllers, or input ring back in clock signals. In mechanical systems, ringing effects, when combined with resonance, may accelerate breakdown. In audio systems, ringing may cause echoes to occur before transients.

In more abstract terms, ringing may be understood as a tendency of the system to store energy, which is retrieved later to produce undesired effects. One way to quantify this capacity uses the Hankel norm of a system, which measures the effects of past inputs on future outputs.

This paper focuses on the problem of minimizing system ringing by casting it as a Hankel norm minimization program. This leads to an eigenvalue optimization problem, for which we propose a nonsmooth bundle algorithm which assures convergence to a critical point from an arbitrary starting point. We demonstrate that a variety of problems such as Hankel synthesis, maximizing the memory of a system, and control of flow in a graph, can be interpreted as Hankel norm minimization programs and solved efficiently using the proposed algorithm.

There is a considerable body of literature dedicated to Hankel norm system reduction, the original contribution being [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞ -error bounds[END_REF]. Our present approach is complementary to this classical line, as we focus on Hankel norm optimization problems which cannot be solved by linear algebra techniques. This makes our method closer in spirit to H 2 -or H ∞ -controller or filter design [START_REF] Zhou | Robust and optimal control[END_REF].

The structure of the paper is as follows. After presenting the problem in abstract form in Sect. 2, we show in Sect. 3 how it can be cast as a nonconvex eigenvalue optimization program. Section 4 describes how Clarke subgradients of a Hankel norm objective can be computed. In Sect. 5 we extend the Hankel norm to systems with direct transmission in a physically meaningful way. Sections 6, 7 present typical applications for the purpose of motivation of the Hankel minimization problem. Section 8 discusses a proximal bundle algorithm used to solve the Hankel norm minimization program. We propose a smooth relaxation of the Hankel norm in Sect. 9. Experiments with typical applications are given in Sect. 10.

Notation

Terminology in nonsmooth optimization is covered by [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF], system theory by [START_REF] Zhou | Robust and optimal control[END_REF]. Following the latter reference, given a transfer matrix function G(s) = C(sI -A) -1 B + D, we use the standard notations

G(s) = A B C D or G = (A, B, C, D)
to indicate that

G : ẋ = Ax + Bw z = Cx + Dw
is a state-space realization of z(s) = G(s)w(s). Similar notations apply to discrete time systems. We shall work in the set of rectangular matrices with the corresponding scalar product M, N = Tr(M N ) = Tr(N M ), where M and Tr(M ) are transpose and trace of a matrix. For symmetric matrices, M 0 means positive definite, M 0 positive semidefinite.

Hankel norm minimization

Consider a linear time-invariant system

G : ẋ = Ax + Bw z = Cx
with state x ∈ R nx , input w ∈ R m , and output z ∈ R p . Suppose G is internally stable in the sense that all eigenvalues of A have negative real part. If we think of w(t) as an excitation at the input which acts over the time period 0 t T with dynamics started at x(0) = 0, then the ring of the system after the excitation has stopped at time T is z(t) for t > T . If signals are measured in the energy norm, this leads to the definition of the Hankel norm of an internally stable system G = (A, B, C) with input w and output z = Gw as

G H = sup T >0 ∞ T z z dt 1/2 : T 0 w w dt 1, w(t) = 0 for t > T .
For the discrete time case, the Hankel norm of an internally stable system G :

x(t + 1) = Ax(t) + Bw(t)

z(t) = Cx(t)
is given by

G H = sup T >0    ∞ t=T z(t) z(t) 1/2 : T t=0 w(t) w(t) 1, w(t) = 0 for t > T    ,
where now internally stable means that all eigenvalues of A have magnitude < 1, and where it is again understood that z = Gw. A formula which works in both cases is

G H = sup T >0
z 2,[T,∞) : w 2,[0,T ] 1, w ∈ L 2 [0, T ], w(t) = 0, t > T . [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF] Note that the system G in the above definition has no direct transmission D. This accounts for the fact, proved in Lemma 2 in Sect. 5, that D causes no memory effects, and is therefore not seen by the Hankel norm [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF]. In consequence, on the space of systems G = (A, B, C, D) with direct transmission,

• H is only a semi-norm and not a norm. By definition, the Hankel norm can be interpreted as a measure of the effects of past inputs, that is, the memory of the system, on the states and future outputs. Here, we are interested in systems G(x) with tunable parameters x ∈ R n , where the matrices A(x), B(x), C(x) depend smoothly on a design parameter x varying in R n or in some constrained subset of R n . Our goal is to tune x such that system ringing is avoided or reduced while internal stability of the system is guaranteed. This leads to the following Hankel norm minimization program minimize G(x) H subject to G(x) internally stable

x ∈ R n .

(2)

We will discuss various instances, where program (2) may be of interest. Then, we present a nonsmooth optimization method based on techniques from eigenvalue optimization to solve [START_REF] Apkarian | Mixed H 2 /H∞ control via nonsmooth optimization[END_REF], and discuss a smooth relaxation motivated by a result of Nesterov in [START_REF] Nesterov | Smoothing technique and its applications in semidefinite optimization[END_REF].

Γ G : L 2 (-∞, 0] -→ L 2 [0, ∞) be the
Hankel operator associated with G, defined by

(Γ G w)(t) = 0 -∞ Ce A(t-τ ) Bw(τ )dτ, t 0.
Then, the following definitions are equivalent:

(i) G H = sup T >0 z 2,[T,∞) : w 2,[0,T ] 1, w ∈ L 2 [0, T ], w(t) = 0, t > T . (ii) G H = Γ G = sup Γ G w 2,[0,∞) : w 2,(-∞,0] 1, w ∈ L 2 (-∞, 0] . (iii) G H = λ 1 (XY )
, where λ 1 denotes the maximum eigenvalue of a matrix, and X, Y are the controllability and observability Gramians of the system.

Proof We assume x(-∞) = 0 for the Hankel operator Γ G and obtain

z(t) = t -∞
Ce A(t-τ ) Bw(τ )dτ.

If we now focus on input signals w -that live for times t 0 and vanish for t > 0, then the output restricted to t 0 is

z + (t) = 0 -∞ Ce A(t-τ ) Bw -(τ )dτ = Γ G w -, t 0.
Assuming x(0) = 0 in (i), it now follows from the time-invariance that sup

T >0 0 =w∈L 2 [0,T ] w(t)=0, t>T z 2,[T,∞) w 2,[0,T ] = sup T >0 0 =w∈L 2 [-T,0] w(t)=0, t>0 z 2,[0,∞) w 2,[-T,0] = sup 0 =w∈L 2 (-∞,0] w(t)=0, t>0 z 2,[0,∞) w 2,(-∞,0] = sup 0 =w-∈L 2 (-∞,0] z + 2,[0,∞) w -2,(-∞,0] = Γ G .
This gives the equivalence of (i) and (ii). Next, we have

w, Γ * G z L 2 (-∞,0] = Γ G w, z L 2 [0,∞) = ∞ 0 0 -∞ w(τ ) B e A (t-τ ) C dτ z(t)dt = 0 -∞ w(τ ) ∞ 0 B e A (t-τ ) C z(t)dt dτ, which implies (Γ * G z)(τ ) = ∞ 0 B e A (t-τ ) C z(t)dt, τ 0.
Note that the operator norm of Γ G is equal to its maximum singular value. Therefore, to complete the proof, we show that

σ 2 i (Γ G ) = λ i (XY )
, where σ i (•) and λ i (•) denote, respectively, the ith singular value and ith eigenvalue of an operator or matrix. Suppose σ is a nonzero singular value of Γ G , and w is an eigenvector corresponding to the eigenvalue

σ 2 of Γ * G Γ G , i.e., Γ * G Γ G w = σ 2 w. Setting z(t) = (Γ G w)(t) = Ce At x 0 with x 0 = 0 -∞ e -Aτ
Bw(τ )dτ , and noting by [START_REF] Zhou | Robust and optimal control[END_REF]Lemma 3.18] that

X = ∞ 0 e At BB e A t dt, Y = ∞ 0 e A t C Ce At dt, we have σ 2 w = Γ * G z = B e -A τ ∞ 0 e A t C z(t)dt = B e -A τ ∞ 0 e A t C Ce At x 0 dt = B e -A τ Y x 0 .
It follows that

σ 2 x 0 = 0 -∞ e -Aτ Bσ 2 w(τ )dτ = 0 -∞ e -Aτ BB e -A τ Y x 0 dτ = XY x 0 .
Moreover, x 0 = 0 since otherwise σ 2 w = 0, which is impossible. Thus, σ 2 is an eigenvalue of XY . Conversely, if σ 2 = 0 is an eigenvalue and x 0 = 0 is a corresponding eigenvector of XY , i.e., XY x 0 = σ 2 x 0 , then by setting w = B e -Aτ Y x 0 we obtain w = 0 and

Γ * G Γ G w = σ 2 w. Hence, σ 2 i (Γ G ) = λ i (XY ), and so Γ G = σ 1 (Γ G ) = λ 1 (XY ).
The lemma is proved.

Lemma 1 shows that the Hankel norm can be considered as a measure of controllability and observability of the system, and that it does not depend on the state-space representation of the system. It is now clear that problem (2) may be cast as an eigenvalue optimization program. In the sequel, we examine how this problem can be solved algorithmically.

Subgradients of the Hankel norm

In this section, we compute Clarke subgradients [8, Section 2.1] of the nonconvex composite function f (x) = G(x) 2 H . This is a fundamental tool for our optimization method.

Let G(x) be a linear time-invariant system with state-space realization (A(x), B(x), C(x)) depending smoothly on a design parameter x ∈ R n . Let X(x), Y (x) be the controllability and observability Gramians. Suppose the maximum eigenvalue λ 1 (Z(x)) of the matrix Z(x) = X(x)

1 2 Y (x)X(x)
1 2 has multiplicity r(x), and let R = R(x) be a matrix whose columns form an orthonormal basis of the eigenspace associated with λ 1 (Z(x)). For any matrix

function M (x), put M k (x) = ∂M (x) ∂x k and write M 1 2 k for (M 1 2 ) k , k = 1, . . . , n.
We have the following

Proposition 1 The function f (x) = G(x) 2
H is well defined and locally Lipschitz on the set S = {x ∈ R n : A(x) stable}. In addition, for every x in the set S 0 = {x ∈ S : (A(x), B(x)) controllable} the Clarke subgradients of f at x have the form

g U = Tr(U R Z 1 (x)R) . . . Tr(U R Z n (x)R) , ( 3
)
where U is symmetric of size r × r, U 0, Tr(U ) = 1, and where the partial derivatives Z k (x), k = 1, . . . , n are given by

Z k (x) = X 1 2 k (x)Y X 1 2 + X 1 2 Y k (x)X 1 2 + X 1 2 Y X 1 2 k (x). (4)
Here,

X k (x), Y k (x) and X 1 2
k (x) are the solutions of the following Lyapunov equations

AX k (x) + X k (x)A = -A k (x)X -XA k (x) -B k (x)B -BB k (x) , (5) A Y k (x) + Y k (x)A = -A k (x) Y -Y A k (x) -C k (x) C -C C k (x), (6)
X 1 2 X 1 2 k (x) + X 1 2 k (x)X 1 2 = X k (x). ( 7
)
Proof 1. By Lemma 1,

f (x) = G(x) 2 H = λ 1 (X(x)Y (x)),
where the Gramians X(x) and Y (x) depend on the tunable parameters x and are the solutions of the Lyapunov equations

A(x)X + XA(x) + B(x)B(x) = 0, (8)
A(x) Y + Y A(x) + C(x) C(x) = 0. ( 9
)
Note that despite the symmetry of X and Y the product XY is not necessarily symmetric, but stability of A(x) guarantees X 0, Y 0 in (8), [START_REF] Dao | Minimizing the memory of a system[END_REF], so that we can write

λ 1 (XY ) = λ 1 (X 1 2 Y X 1 2 ) = λ 1 (Y 1 2 XY 1 2 ),
which brings us back to the realm of eigenvalue theory of symmetric matrices. By positive semidefiniteness of X(x) and Y (x), the function f is now well defined on S. 2. Let us next prove that f is locally Lipschitz on S. Using the Kronecker product [START_REF] Bellman | Kronecker products and the second method of Lyapunov[END_REF], Eq. ( 8) can be written as

(I ⊗ A(x) + A(x) ⊗ I)vec(X(x)) = -vec(B(x)B(x) ),
where I is a conformable identity matrix, and where vec(•) vectorizes a matrix by stacking its columns in order. Since A(x) is smooth in x and M (x) = (I ⊗ A(x) + A(x) ⊗ I) is invertible by the stability of A(x), M (x) -1 is also smooth in x, and since B(x) depends smoothly on x, then so does vec(X(x)) = -M (x) -1 vec B(x)B(x) . A similar argument shows smooth dependence of Y (x) on x. This can also be justified based on the explicit formulas

X(x) = ∞ 0 e A(x)t B(x)B(x) e A(x) t dt, Y (x) = ∞ 0 e A(x) t C(x) C(x)e A(x)t dt
(see e.g., [START_REF] Zhou | Robust and optimal control[END_REF]Lemmas 2.7 and 3.18]), where uniform convergence of these integrals on any bounded set of x gives differentiability in x. We infer that the coefficients of the characteristic polynomial of X(x)Y (x) also depend smoothly on x. Since this characteristic polynomial is hyperbolic, that is, has only real roots, we may invoke the multi-parameter version of Bronstein's theorem [START_REF] Bronshtein | Smoothness of roots of polynomials depending on parameters[END_REF], for which an elegant proof is given in [START_REF] Parusiński | A new proof of Bronshtein's theorem[END_REF]Theorem 2], to conclude that f (x) = λ 1 (X(x)Y (x)) is locally Lipschitz on S. 3. Let us finally establish formula (3) for the subdifferential ∂f (x) at points

x ∈ S 0 . By the above argument, f (x) = λ 1 (Z(x)). Observe that controllability of (A(x), B(x)) implies that X(x) is positive definite [26, Theorem 3.1], and since the operator X → X 1 2 is smooth on the set of matrices X 0, the chain rule gives smoothness of x → X 1 2 (x), and so of Z

(x) = X 1 2 Y X 1 2
, on S 0 . Applying [START_REF] Overton | Large-scale optimization of eigenvalues[END_REF]Theorem 3], the Clarke subgradients of f at x are of the form g U = g 1 . . . g n , where

g k = U, R Z k (x)R = Tr(U R Z k (x)R)
for U symmetric of size r × r, U 0, Tr(U ) = 1. It now remains to calculate Z k (x), k = 1, . . . , n. We first have (4) by the definition of Z. Taking derivatives with respect to x on both sides of (8)-( 9), we get (5)- [START_REF] Bronshtein | Smoothness of roots of polynomials depending on parameters[END_REF], and then also X k (x), Y k (x). Finally, to compute X 1 2 k (x), we use [START_REF] Burke | Differential properties of the spectral abscissa and the spectral radius for analytic matrix-valued mappings[END_REF], which is obtained by differentiating X 1 2 X 1 2 = X. Altogether, we obtain Clarke subgradients of f at each x due to (3)- [START_REF] Dao | Minimizing the memory of a system[END_REF]. Remark 1 Formula (3) also holds if controllability of (A(x), B(x)) is replaced by observability of (A(x), C(x)) (cf. [START_REF] Zhou | Robust and optimal control[END_REF]Definition 3.4]). Here, we work with

Z = Y 1 2 XY 1 2 instead.
Remark 2 In the discrete time case, the Gramians X(x) and Y (x) are the solutions of the discrete Lyapunov equations

A(x)XA(x) -X + B(x)B(x) = 0, A(x) Y A(x) -Y + C(x) C(x) = 0,
so that X k (x) and Y k (x) are solutions, respectively, of the following equations

AX k (x)A -X k (x) = -A k (x)XA -AXA k (x) -B k (x)B -BB k (x) , A Y k (x)A -Y k (x) = -A k (x) Y A -A Y A k (x) -C k (x) C -C C k (x).
Remark 3 Subgradients of f at x ∈ S \ S 0 are no longer represented by (3), since the solution of (7) need not exist, as only X 1 2 0 is guaranteed. Nonetheless, by Clarke subdifferentiability at points x ∈ S \S 0 proved above, we can be sure that for every sequence x k ∈ S 0 converging to x ∈ S \ S 0 and g k ∈ ∂f (x k ) computed via (3), the g k stay bounded and each of their accumulation points g is an element of ∂f (x). This guarantees stability of our numerical procedure even when iterates get close to the set S \ S 0 .

Remark 4 Practical parametrizations G(x) use elementary computable operations, which can be expressed in mathematical terms by assuming that A(x), B(x), C(x) are smooth definable functions of x in the sense of [START_REF] Van Den Dries | Tame topology and o-minimal structures[END_REF]Chap. 1,Sect. 5.3]. In that case, one can say a little more about the behavior of f at points x ∈ S. Namely, it then follows from [21, Theorem 4.12] that for every smooth definable curve x(t) ∈ S the eigenvalues λ i (t) = λ i (X(x(t))Y (x(t))) are smooth functions of t, so that f (x(t)) is a finite maximum of smooth functions of t. On S 0 this property is a consequence of symmetric eigenvalue theory, which is true without the definability hypothesis. Note that this does not mean that f is a finite maximum of smooth functions of x ∈ R n , but it nonetheless indicates a favorable structure.

An extension of the Hankel norm

Lemma 1 shows why the Hankel norm is only a semi-norm on the space of internally stable systems G. It does not see a direct transmission D from w to z, as the latter does not create memory transmitted from the past to the future. This rises the question how to assess a direct transmission block in the context of (1) or [START_REF] Apkarian | Mixed H 2 /H∞ control via nonsmooth optimization[END_REF]. Namely, in some applications, attributing no cost to a block D(x) which is free to vary with the tunable parameters x bears the risk that optimization favors a solution with a high energy direct transmission.

It is well known that G H G ∞ in the case D = 0 (See e.g., [5, Sect. 5.5]), and this may guide us to define an extension. Note first that Lemma 2 (A, B, C) H (A, B, C, D) ∞ for every internally stable system G = (A, B, C, D).

Proof Let G 0 = (A, B, C) be the system where the direct transmission is ignored. Consider an input w with w(t) = 0 for t > T , and let z 0 = G 0 w, z = Gw. Then, z(t) = z 0 (t) for t > T , because the direct transmission creates no memory, and since w(t) = 0 for t > T , its influence on the output ends at T . Combining this with w 2,[0,T ] = w 2 and z 2,[T,∞) z 2 , we obtain

(A, B, C) H = sup T >0 0 =w∈L 2 [0,T ] w(t)=0, t>T z 2,[T,∞) w 2,[0,T ] sup T >0 0 =w∈L 2 [0,T ] w(t)=0, t>T z 2 w 2 sup w =0 z 2 w 2 = (A, B, C, D) ∞ .
This suggests the following extension of Hankel norm

• H to systems G = (A, B, C, D) with direct transmission D.
Definition 1 Let G = (A, B, C, D) be an internally stable linear time-invariant system. Then,

G H = max { (A, B, C) H , σ 1 (D)} (10)
is called the extended Hankel norm of the system. Here, σ 1 denotes the maximum singular value of a matrix.

This definition agrees with the usual Hankel norm for a system without direct transmission, and also preserves the inequality

G H G ∞ , since the term σ 1 (D) is part of the maximum G ∞ = max ω σ 1 (G(jω)) at ω = ∞.
As the proof of Lemma 2 shows, a direct transmission does not change the value of • H defined according to (1). In the sequel, we therefore adopt the convention that in the case D = 0, (A, B, C) H is the usual Hankel norm, where the direct transmission is ignored, while (A, B, C, D) H is the extended Hankel norm.

An advantage of (10) is that the new function is still a maximum eigenvalue function. Namely, stability of G implies positive semidefiniteness of the Gramians X and Y , and so

G 2 H = max λ 1 (X 1 2 Y X 1 2 ), λ 1 (D D) = λ 1 X 1 2 Y X 1 2 0 0 D D . (11)
Proceeding as in the proof of Proposition 1, we get immediately the following Corollary 1 Let G(x) be a linear time-invariant system depending smoothly on x ∈ S with S = {x ∈ R n : A(x) stable}. Suppose the maximum eigenvalue λ 1 (Z(x)) of the matrix

Z(x) = X(x) 1 2 Y (x)X(x) 1 2 0 0 D(x) D(x)
has multiplicity r = r(x), and R = R(x) is a matrix whose columns form an orthonormal basis of the eigenspace associated with λ 1 (Z(x)). With the notations of Proposition 1, the function f (x) = G(x) 2 H is locally Lipschitz on S and its Clarke subgradients on S 0 = {x ∈ S : (A(x), B(x)) controllable} have the form

g U = Tr(U R Z 1 (x)R) . . . Tr(U R Z n (x)R) ,
for U symmetric of size r × r, U 0, Tr(U ) = 1, where the partial derivatives Z k (x), k = 1, . . . , n are given by

Z k (x) = Z k (x) 0 0 D k (x) D(x) + D(x) D k (x)
and the Z k (x) are defined in Proposition 1.

To justify the use of (10) rigorously, we consider the extended Hankel norm minimization program (2) based on (10), and compare it to the following constraint program

minimize f (x) = (A(x), B(x), C(x)) H subject to h(x) = σ 1 (D(x)) η. ( 12
)
For the following, recall from [START_REF] Mangasarian | Nonlinear programming[END_REF]

that x * ∈ R n is called a Fritz John critical point of the constraint program min{f (x) : h(x) η} if there exist multipliers λ * 0 0, λ * 1 0, not both zero, such that 0 ∈ λ * 0 ∂f (x * ) + λ * 1 ∂h(x * ), h(x * ) η, λ * 1 (h(x * ) -η) = 0. If in addition λ * 0 > 0, then x * is called a Karush-Kuhn-Tucker point.
Remember that every local minimum x * of the constraint program is automatically a Fritz John critical point, while it will in general only be a Karush-Kuhn-Tucker point if an additional constraint qualification is satisfied [START_REF] Mangasarian | Nonlinear programming[END_REF]Chapter 7]. For later on, we call x * a critical point of constraint violation if 0 ∈ ∂h(x * ) and h(x * ) > η.

With these preparations, we have the following Proposition 2 Let x * be a critical point of the extended Hankel norm minimization program (2) with (10). Then, x * is a Fritz John critical point of program [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞ -error bounds[END_REF] for a suitable choice of η. More precisely, x * is either a Karush-Kuhn-Tucker point of (12), or a critical point of h

(x) = σ 1 (D(x)) alone. Proof Note that G(x) H = max{f (x), h(x)}. Now, if x * is a critical point of G(x) H , then we have three possibilities, f (x * ) > h(x * ), f (x * ) = h(x * ), or f (x * ) < h(x *
). In the first case, x * is a critical point of f alone, hence also a Karush-Kuhn-Tucker point of [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞ -error bounds[END_REF]. The third case corresponds to a critical point of h alone. In the case of equality, the situation is more complex. There exist multipliers

λ * 0 0, λ * 1 0, not both zero, such that 0 ∈ λ * 0 ∂f (x * ) + λ * 1 ∂h(x * ). If λ * 0 = 0 then λ * 1 = 0 and 0 ∈ ∂h(x * ), so x * is a critical point of h. In case λ * 0 = 0, we have 0 ∈ ∂f (x * ) + (λ * 1 /λ * 0 )∂h(x *
). This is the first part of the Karush-Kuhn-Tucker conditions. If we put η = f (x * ), then we also get the second half. That completes the argument.

Remark 5 Suppose we solve program min{f (x) : h(x)

η} starting at an infeasible point h(x 1 ) > η, then we will usually try to minimize h alone to find a feasible iterate. Suppose a descent method used to minimize h runs into a local minimum x * of h satisfying h(x * ) > η. Such a local minimum of constraint violation indicates a failure, since nothing better will be found in a neighborhood of x * due to local optimality, so that the search for a feasible point has to be stared anew elsewhere; cf. [20, Section 2.2] for this theme complex.

By Proposition 2 we can now interpret minimization of the extended Hankel norm ( 2) with (10) as a trade-off between minimizing the memory effects of (A(x), B(x), C(x)), subject to a constraint σ 1 (D(x))

η, or dually, as of minimizing σ 1 (D(x)) subject to a constraint on the memory effects of G(x). Since f (x) is a valid measure of the memory or ringing effects of G(x), such an interpretation is physically meaningful.

We conclude this section by showing that the Hankel norm is amenable to optimization techniques, as this will be needed later. According to Spingarn [24] a function f : U → R, where U is an open set in R n , is lower -C 1 on U , if for each x 0 ∈ U , there are a compact space K, a neighborhood V of x 0 , and a jointly continuous function F : V × K → R whose partial derivative D x F with respect to x exists and is jointly continuous, such that f (x) = max z∈K F (x, z) for all x ∈ V .

Proposition 3 Let G(x) = (A(x), B(x), C(x), D(x)
) be a linear time-invariant system depending smoothly on the set S 0 of all x ∈ R n such that A(x) is stable and

(A(x), B(x)) is controllable or (A(x), C(x)) is observable. Then, f (x) = G(x) 2 H is lower-C 1 on S 0 .
Proof For each x ∈ S 0 , according to [START_REF] Gabarrou | Design of a flight control architecture using a non-convex bundle method[END_REF] and using the Rayleigh quotient,

f (x) = λ 1 (Z(x)) = max z =1 z Z(x)z,
where Z is symmetric and depends smoothly on x. Set K = {z ∈ R m : z = 1} and F (x, z) = z Z(x)z, then K is compact, f (x) = max z∈K F (x, z), and both F and its partial derivatives F x are jointly continuous on S 0 × K and smooth in x. Therefore, f is lower-C 1 on S 0 .

Hankel synthesis

The first application of program (2) we consider is output feedback controller synthesis, where performance is assessed by the Hankel norm. Consider a linear time-invariant plant in standard form

P (s) :   ẋ z y   =   A B 1 B 2 C 1 D 11 D 12 C 2 D 21 D 22     x w u   , (13)
where x ∈ R nx is the state, u ∈ R m2 the control, w ∈ R m1 the vector of exogenous inputs, y ∈ R p2 the measurements, and z ∈ R p1 the controlled or performance vector, P (s) := P 11 (s) P 12 (s)

P 21 (s) P 22 (s) = C 1 C 2 (sI -A) -1 B 1 B 2 + D 11 D 12 D 21 D 22 .
Without loss of generality, it is assumed that D 22 = 0. Let u(s) = K(s)y(s) be an output feedback controller for the open-loop plant [START_REF] Mangasarian | Nonlinear programming[END_REF], with

K : ẋK u = A K B K C K D K x K y ,
where x K ∈ R k is the state of K. The closed-loop transfer function of the performance channel w → z is obtained as

T w→z (K, s) = P 11 (s) + P 12 (s)K(s)(I -P 22 (s)K(s)) -1 P 21 (s).
Our aim is to find an optimal controller K which stabilizes the system in closed-loop such that T w→z (K) H is minimized among all stabilizing K. By substituting u = Ky into (13), the state-space representation of the closed-loop performance channel w → z is

T w→z (K) : ξ z = A(K) B(K) C(K) D(K) ξ w ,
where ξ = (x, x K ) and

A(K) = A + B 2 D K C 2 B 2 C K B K C 2 A K , B(K) = B 1 + B 2 D K D 21 B K D 21 , C(K) = C 1 + D 12 D K C 2 D 12 C K , D(K) = D 11 + D 12 D K D 21 .
This problem is now a specific instance of [START_REF] Apkarian | Mixed H 2 /H∞ control via nonsmooth optimization[END_REF], where in agreement with our general theme we try to minimize the memory of a specific channel w → z within the plant P . If we allow structured control laws K(x) in the sense of [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF], then we obtain the following optimization program minimize T w→z (K) H subject to K stabilizes [START_REF] Mangasarian | Nonlinear programming[END_REF] internally

K = K(x), x ∈ R n . ( 14
)
Example 1 Typical examples of structured controllers are, for instance, PIDs or observer-based controllers, which in state-space have the form

K pid (x) =   0 0 r i 0 -τ r d 1 1 d K   , K obs (x) = A + B 2 K c + K f C 2 -K f K c 0 .
For a PID, the tunable parameters are x = (r i , r d , d K , τ ), while for observerbased controllers K obs (x) the vector x gathers the elements of K c , K f . Other examples are decentralized, fixed reduced order controllers, and more generally, control architectures combining basic building blocks such as PIDs with filters, feed-forward blocks, and much else (see [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF]).

Remark 6

The norm in program ( 14) is the usual Hankel norm (1) if D(K) = 0, which is the case e.g., under standard assumption as in H 2 -synthesis, where D 11 = 0 and either D 21 = 0 or D 12 = 0 or K strictly proper. In contrast, if D(K) = 0, then we should use the extended Hankel norm (10), or likewise, the constraint program [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞ -error bounds[END_REF], to control the direct transmission. It is also possible to neglect the direct transmission term D(K) and optimize the semi-norm (A(K), B(K), C(K)) H . We then exercise caution by monitoring the term σ 1 (D(K)) during optimization to check whether a large direct transmission gain σ 1 (D(K)) is favored. If that is the case, switching to the extended Hankel norm becomes mandatory.

In the sequel of this section, we discuss two particular cases of the Hankel synthesis problem (14).

System reduction

System reduction is the most widely known application of the Hankel norm minimization problem. Given a stable system

G : ẋ = Ax + Bw z = Cx + Dw
of order n x , we wish to find a stable system

G k : ẋ = A k x + B k w z = C k x + Dw
G -G k (x) H subject to G -G k (x) internally stable x = (A k , B k , C k ) ( 15
)
is a particular case of ( 14), where we define plant and controller as

P :   A B 0 C D -I 0 I 0   , K : A k B k C k D , (16)
the tunable parameters x being the elements of

A k , B k and C k . E E G G red c T c + - E e w (17)
Due to the seminal work of Glover [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞ -error bounds[END_REF], program [START_REF] Nesterov | Smoothing technique and its applications in semidefinite optimization[END_REF] has an explicit solution based on linear algebra, at least when no additional structural constraints on the matrices A k , B k , C k are imposed. This allows us to implement a blind testing of Algorithm 1 in Sect. 8, which is applied to [START_REF] Nesterov | Smoothing technique and its applications in semidefinite optimization[END_REF], considered as a particular case of ( 14) using [START_REF] Noll | Cutting plane oracles to minimize non-smooth non-convex functions[END_REF]. The value obtained by Algorithm 1 is then compared to the theoretical value obtained by an explicit Hankel system reduction.

Maximizing the memory of a system

Within the present framework, it is also possible to maximize the memory effects of a system G via feedback if a reference system G ref with desirable memory properties is used. In other words, while minimizing G(x) H leads to a system which is the least biased, we now bias G(x) as much as possible by bringing it as close as possible to G ref , and we achieve this by making G(x) -G ref as less biased as possible.

Example 2 As a motivating example, we consider a 2-DOF synthesis scheme of the following form

E q E F K E E c d E E z 2 z 1 G d E T y u v w e - + + - E G ref y ref T d E (18)
where the decentralized controller structure was chosen to challenge our method in a typical situation in practice.

Assuming that G ref has desirable memory features which do not lead to ringing, the idea is to tune the parameters in feed-forward filter F and controller K in such a way that G in closed-loop follows G ref , independently of the input w. That is, the undesirable part of the memory of G, which contributes to the mismatch z 1 = y -y ref , is reduced by minimizing T w→z1 (F, K) H . It may be beneficial to arrange this by adding a constraint z 2 2 η 2 or z 2 ∞ η ∞ , where z 2 = u + v, to avoid exceedingly large controller actions. This problem can be cast as a particular case of program [START_REF] Mifflin | A modification and extension of Lemaréchal's algorithm for nonsmooth minimization. Nondifferential and variational techniques in optimization[END_REF] if the following plant and decentralized controller structures are used

P :         A 0 0 B B 0 A ref B ref 0 0 C -C ref -D ref D D 0 0 0 I I -C 0 I -D -D 0 0 I 0 0         , K :     A F 0 B F 0 0 A K 0 B K C F 0 D F 0 0 C K 0 D K     .
Notice that

F : ẋF = A F x F + B F w v = C F x F + D F w , K : ẋK = A K x K + B K e u = C K x K + D K e
can be further structured if we wish. In our experiment, we will use this example with F a reduced-order filter, and K a PID.

Control of flow in a graph

We consider the flow in a directed graph G = (V , A ) with interior nodes, sources and sinks, V = V stay ∪ V in ∪ V out , and not excluding self-arcs. For nodes i, j ∈ V connected by an arc (i, j) ∈ A the transition probability i → j quantifies the tendency of flow going from node i towards node j. As an example consider for instance a large fairground with separated entrances and exits, with itineraries represented by the graph. By acting on the transition probabilities between nodes connected by arcs, we expect to guide the crowd in such a way that a steady flow is assured, and a safe evacuation is possible.

Assume that an individual at interior node j ∈ V stay decides with probability a jj 0 to proceed to a neighboring node j ∈ V stay , where neighboring means (j, j ) ∈ A , or with probability a jk 0 to move to a neighboring exit node k ∈ V out , where (j, k) ∈ A . The case (j, j) ∈ A of deciding to stay at stand j ∈ V stay is not excluded. Similarly, an individual entering at i ∈ V in proceeds to a neighboring interior node j ∈ V stay with probability b ij 0, where (i, j) ∈ A . We suppose for simplicity that there is no direct transmission from entrances to exits. Then, j ∈Vstay:(j,j )∈A

a jj + k∈Vout:(j,k)∈A a jk = 1, (19)
for every j ∈ V stay , and

j∈Vstay:(i,j)∈A b ij = 1 ( 20
)
for every i ∈ V in . Let x j (t) denote the number of people present at interior node j ∈ V stay and time t, and w i (t) the number of people entering the fairground through entry i ∈ V in at time t. Then, the number of people present at interior node j ∈ V stay and time t + 1 is

x j (t + 1) = j ∈Vstay:(j ,j)∈A

a j j x j (t) + i∈Vin:(i,j)∈A b ij w i (t),
while the number of people leaving the fairground at time t through exit k ∈ V out is j∈Vstay:(j,k)∈A a jk x j (t). To assess the evacuation pattern, we quantify the total number of people still inside the fairground at time t via the weighted sum

z(t) = j∈Vstay c j x j (t),
where c j > 0 are fixed weights, and where c j = 1 would correspond to simply counting the number of people inside the fairground. We let x regroup the parameters a jj , a jk , b ij , so that the discrete linear time-invariant system has the form G(x) = (A(x), B(x), C), where C is the row vector of c j 's.

Let us now consider an evacuation scenario, where at time T the inflow w(t) through the entrance gates is stopped by closing the gates, and the time until the fairground is evacuated is assessed by measuring the evacuation pattern z(t), t > T . This corresponds to computing the Hankel norm G(x) H , which identifies the worst case evacuation scenario. Minimizing z 2,[T,∞) / w 2,(0,T ] may then be understood as enhancing overall safety of the network by orienting the crowd in such a way that the worst case evacuation time is minimized. This leads to the optimization program minimize G(x) H subject to G(x) internally stable a jj 0, a jk 0, b ij 0, ( 19), (

which is a discrete version of (2) including linear constraints. Notice that these linear constraints are readily added in our algorithmic approach. In an extended model, one might consider measuring the number of people y at some selected nodes i ∈ V stay ∪ V out , and use this to react via feedback u = Ky at the entry gates. This leads to a problem where controller and parts of the plant are optimized simultaneously. Other variants include cases, where some of the probabilities a jj , b ij are imposed and cannot be modified by the designer.

Proximal bundle algorithm

In this section, we present our main algorithm to solve programs (2) and [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞ -error bounds[END_REF]. Let us consider an abstract constrained optimization program of the form minimize f (x) subject to h(x) 0 [START_REF] Rockafellar | Variational analysis[END_REF] where x ∈ R n is the decision variable, and f and h are locally Lipschitz but potentially nonsmooth and nonconvex functions, representing objective and constraints. To find solutions of the constraint program [START_REF] Rockafellar | Variational analysis[END_REF], using an idea inspired by Polak [20, Section 2.2.2], we introduce the progress function

F (y, x) = max{f (y) -f (x) -νh(x) + , h(y) -h(x) + },
where h(x) + = max{h(x), 0}, and ν > 0 is some fixed parameter (with ν = 1 a typical value). One can think of x as the current iterate, and y as the next iterate or as a candidate to become the next iterate. We need to collect a few facts about F . Note first that F (x, x) = 0. For the subdifferential, we have the useful Lemma 3 Suppose f and h are lower-C 1 functions. Then, the Clarke subdifferential of the progress function F with respect to the first variable is obtained as

∂ 1 F (x, x) =      ∂f (x) if h(x) < 0, conv{∂f (x) ∪ ∂h(x)} if h(x) = 0, ∂h(x) if h(x) > 0.
Proof Applying the formula for the Clarke subdifferential of a maximum [8, Proposition 2.3.12] we readily get

∂ 1 F (x, x) = ∂f (x) if h(x) < 0, ∂ 1 F (x, x) ⊂ conv{∂f (x) ∪ ∂h(x)} if h(x) = 0, and ∂ 1 F (x, x) = ∂h(x) if h(x) > 0.
But since f and g are lower-C 1 , according to [24, Proposition 2.4, Theorem 3.9], they are Clarke regular, so we have equality in the second case h(x) = 0.

Lemma 4 Suppose x * is a local minimum of program [START_REF] Rockafellar | Variational analysis[END_REF], then it is also a local minimum of F (•, x * ), and

0 ∈ ∂ 1 F (x * , x * ). Conversely, if 0 ∈ ∂ 1 F (x * , x * )
then x * is either a Karush-Kuhn-Tucker point of [START_REF] Rockafellar | Variational analysis[END_REF], or a critical point of constraint violation.

Proof Since x * is a local minimum of ( 22), we have feasibility h(x * ) 0, and so h(x * ) + = 0, which implies F (y, x * ) = max{f (y) -f (x * ), h(y)}. Now, there exists a neighborhood U of x * such that f (y) f (x * ) for every y ∈ U with h(y) 0. We argue that F (y, x * ) F (x * , x * ) for every y ∈ U . Namely, if h(y) > 0, then F (y, x * ) h(y) > 0 = F (x * , x * ). On the other hand, if h(y) 0, then y is feasible, and we have f (y) f (x * ) by what was said before. But then F (y, x * ) f (y) -f (x * ) 0 = F (x * , x * ). This proves x * is a local minimum of F (•, x * ), and so 0 ∈ ∂ 1 F (x * , x * ).

Next, suppose 0 ∈ ∂ 1 F (x * , x * ), then by Lemma 3, there exist non-negative constants λ * 0 , λ * 1 summing up to 1 such that 0 ∈ λ * 0 ∂f (x * ) + λ * 1 ∂h(x * ). If h(x * ) > 0, we have ∂ 1 F (x * , x * ) = ∂h(x * ), and then 0 ∈ ∂h(x * ), meaning that x * is a critical point of h. If h(x * ) < 0 then ∂ 1 F (x * , x * ) = ∂f (x * ), so λ * 1 = 0 and x * is a Karush-Kuhn-Tucker point of [START_REF] Rockafellar | Variational analysis[END_REF]. Assume that h(x * ) = 0 but x * fails to meet the Karush-Kuhn-Tucker conditions, we then obtain λ * 0 = 0 and 0 ∈ ∂h(x * ). This completes the proof of the lemma.

The consequence of this argument is that we should seek points x * with 0 ∈ ∂ 1 F (x * , x * ). We now present our method for computing solutions of program [START_REF] Rockafellar | Variational analysis[END_REF], which is based on this rationale. It generates a sequence x j of estimates which converges to a solution x * in the sense of subsequences. At the current iterate x, the inner loop of the algorithm constructs first-order working models φ k (•, x) and the corresponding second-order working models

Φ k (y, x) = φ k (y, x) + 1 2 (y -x) Q(x)(y -x),
updated with counter k. The Φ k (•, x) are approximations of F (•, x) around x, where Q(x) is symmetric, depends only on the current iterate x, and may reflect second-order information of F around x. The first-order working model

φ k (•, x) has to satisfy φ k (x, x) = F (x, x) = 0 and ∂ 1 φ k (x, x) ⊂ ∂ 1 F (x, x) at all instants k. This is guaranteed when m e (•, x) = g(x) (• -x) with g(x) ∈ ∂ 1 F (x, x
) is an affine minorant of φ k (•, x) at all times k. We refer to m e (•, x) as the exactness plane at x. For a given working model, we solve the tangent program

min y∈R n Φ k (y, x) + τ k 2 y -x 2 ,
with the so-called proximity control parameter τ k > 0. We require Q(x)+τ k I 0, which assures that the tangent program is strictly convex and has a unique solution y k , called the trial step. According to standard terminology, y k is called a serious step if it is accepted as the new iterate y k = x + , and a null step otherwise. Suppose y k is a null step, then we will have to make sure that the next working model φ k+1 (•, x) improves over φ k (•, x). This is achieved by adding cutting and aggregate planes. Let us first look at aggregation. The optimality condition for the tangent program implies

g * k := (Q(x) + τ k I)(x -y k ) ∈ ∂ 1 φ k (y k , x). We call m * k (•, x) = φ k (y k , x) + g * k (• -y k ) = a * k + g * k (• -x) with a * k = φ k (y k , x) + g * k (x -y k ) the aggregate plane. By assuring that m * k (•, x) is an affine minorant of φ k+1 (•, x), we have φ k+1 (y k , x) m * k (y k , x) = φ k (y k , x).
A central element in bundle methods is the cutting plane whose role is to cut away the unsuccessful trial step y k . For each subgradient

g k ∈ ∂ 1 F (y k , x), the affine function t k (•) = F (y k , x) + g k (• -y k ) is a tangent to F (•, x) at y k .
Without convexity, we cannot use t k (•) directly as a cutting plane. Instead, we use a technique first analyzed in [START_REF] Mifflin | A modification and extension of Lemaréchal's algorithm for nonsmooth minimization. Nondifferential and variational techniques in optimization[END_REF], which shifts the tangent down. Fixing a parameter c > 0, we define the cutting plane as

m k (•, x) = t k (•) -s = a k + g k (• -x), (23)
where a k = min{t k (x), -c y kx 2 }, and where s = [t k (x) + c y kx 2 ] + is the downshift. The detailed statement is described as Algorithm 1, while a flowchart of the algorithm is shown in Fig. 1. For more details we refer to [17, Section 3], [START_REF] Noll | Cutting plane oracles to minimize non-smooth non-convex functions[END_REF]Section 4] for unconstrained optimization case, and [2, Section 5], [START_REF] Gabarrou | Design of a flight control architecture using a non-convex bundle method[END_REF]Section 3] for the constrained case. Next, we establish the following result on the convergence of Algorithm 1.

Theorem 1 Suppose that f and h in [START_REF] Rockafellar | Variational analysis[END_REF] are lower-C 1 functions, and let {x ∈ R n : f (x) f (x 1 )} be bounded. Then, every accumulation point x * of the sequence of serious iterates x j generated by Algorithm 1 satisfies 0 ∈ ∂ 1 F (x * , x * ). In other words, x * is either a critical point of constraint violation, or a Karush-Kuhn-Tucker point of [START_REF] Rockafellar | Variational analysis[END_REF].

Proof We will adapt the proof of Theorem 6.6 and Corollary 6.7 in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF] to our needs. For that let us recall a notion from [17, Definitions 2.1 and 6.1], which we apply here to the progress function F . We call φ : R n × S → R a strict first-order model of F on the set S ⊂ R n if for every x ∈ S the function φ(•, x) is convex and the following axioms hold: Algorithm 1 Proximal bundle algorithm with downshifted tangents

Parameters: 0 < γ < γ < Γ < 1, 0 < δ 1, 0 < q < T ∞.
Step 1 (Initialize outer loop). Choose initial feasible guess x 1 , fix memory control parameter τ 1 , and put outer loop counter j = 1.

Step 2 (Stopping test). At outer loop counter j, stop if 0 ∈ ∂ 1 F (x j , x j ). Otherwise, take a symmetric matrix Q j respecting -qI Q j qI, and goto inner loop.

Step 3 (Initialize inner loop). Put inner loop counter k = 1 and initialize control parameter τ 1 = max{τ j , -λ min (Q j ) + δ}, where λ min (•) denotes the minimum eigenvalue of a symmetric matrix. Choose initial working model φ 1 (•, x j ) = g(x j ) (• -x j ) with g(x j ) ∈ ∂ 1 F (x j , x j ).

Step 4 (Tangent program). At inner loop counter k, let Φ k (y, x j ) = φ k (y, x j ) + 1 2 (yx j ) Q j (y -x j ) and find solution y k (trial step) of the tangent program

min y∈R n Φ k (y, x j ) + τ k 2 y -x j 2 .
Step 5 (Acceptance test). Compute the quotient

ρ k = F (y k , x j ) Φ k (y k , x j ) .
If ρ k γ (serious step), put x j+1 = y k and update memory element τ j+1 as τ k if ρ k < Γ , and 1 2 τ k otherwise. Reset τ j+1 = T if τ j+1 > T , increase outer loop counter j and loop back to step 2. If ρ k < γ (null step), continue inner loop with step 6.

Step 6 (Update working model). Generate a cutting plane m k (•, x j ) at null step y k and counter k using downshifted tangents. Compute aggregate plane m * k (•, x j ) at y k , and then build new working model φ k+1 (•, x j ) by adding the new cutting plane, keeping the exactness plane and using aggregation to avoid overflow.

Step 7 (Update control parameter). Compute secondary control parameter

ρ k = M k (y k , x j ) Φ k (y k , x j ) ,
with M k (y, x j ) = m k (y, x j ) + 1 2 (y -x j ) Q j (y -x j ). If ρ k < γ then keep τ k+1 = τ k , otherwise step up τ k+1 = 2τ k . Increase inner loop counter k and loop back to step 4.

(M 1 ) φ(x, x) = F (x, x) = 0 and ∂ 1 φ(x, x) ⊂ ∂ 1 F (x, x).
( M 2 ) If y j → x and x j → x then there exists ε j → 0 + such that F (y j , x j )φ(y j , x j ) ε j y j -x j . (M 3 ) φ is jointly upper semicontinuous on R n × S, i.e., if (y j , x j ) → (y, x) then lim sup j→∞ φ(y j , x j ) φ(y, x).

Representing the cutting plane in (23) as

m y + (•, x) = a+g (•-x) with g ∈ ∂ 1 F (y + , x) and a = min{t y + (x), -c y + -x 2 }, t y + (•) = F (y + , x)+g (•-y + ), we define φ(y, x) = sup{m y + (y, x) : y + ∈ B(x, r)},
where B(x, r) is a fixed ball large enough to contain all possible trial steps, and where the supremum is over all possible cases of m y + (•, x). It then follows that φ is a strict model of F in the sense of the above definition. This can be shown as in [START_REF] Noll | Cutting plane oracles to minimize non-smooth non-convex functions[END_REF]. Axiom ( M 2 ) relies on the fact that F (•, x) is lower-C 1 by the assumptions on f and h. Furthermore, the construction of φ and φ k also guarantees that the working models φ k are lower approximations of φ satisfying

start initialize x 1 , τ ♯ 1 put j = 1 outer loop τ ♯ j+1 := τk/2 j := j + 1 0 ∈ ∂1F (x j , x j ) exit ρk Γ τ ♯ j+1 := τk initialize Qj,
φ k (x, x) = φ(x, x) = F (x, x) = 0, ∂ 1 φ k (x, x) ⊂ ∂ 1 φ(x, x) and φ k (•, x) φ(•, x).
The difference with [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF] is that here the cutting planes m k (•, x) are not directly tangents of φ, but we shall argue that the essential link between φ k and φ rests the same.

The proof now follows essentially [17, Theorem 6.6, Corollary 6.7], which assures that every accumulation point x * of the iterates x j satisfies 0 ∈ ∂ 1 F (x * , x * ). Note that f (x j ) and f (y k ) used in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF] have to be replaced by F (x j , x j ) = 0 and F (y k , x). The fact that Φ(y k+1 ) in the definition of ρ k in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF] is changed to M k (y k , x) can be treated using the property that if y j → x and x j → x then there exists ε j → 0 + such that F (y j , x j ) -m yj (y j , x j ) ε j y j -x j , as follows from [16, Lemma 8], using again crucially that F (•, x) is lower-C 1 . The equality φ k+1 (y k+1 , x) = φ(y k+1 , x) used in the proof of [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF]Lemma 4.2] is now replaced by φ k+1 (y k , x) m k (y k , x). Finally, Lemma 4 completes the last statement of the theorem.

A smooth relaxation of the Hankel norm

Here, we introduce a smooth relaxation of the Hankel norm based on a result of Nesterov in [START_REF] Nesterov | Smoothing technique and its applications in semidefinite optimization[END_REF]. He provides a fine analysis of the convex bundle method in situations where the objective f (x) has the specific structure of a maxfunction, including the case of a convex maximum eigenvalue function. These findings indicate that for a given precision, such programs may be solved with lower algorithmic complexity using smooth relaxations. While these results are a priori limited to the convex case, it may be interesting to apply Nesterov's idea as a heuristic in the nonconvex situation. This leads to the following Proposition 4 Let Z be a symmetric matrix of order m depending smoothly on a parameter x ∈ R n with eigenvalues λ 1 (Z) • • • λ m (Z). Then, for a tolerance parameter µ > 0, the function

f µ (x) = µ ln m i=1 e λi(Z(x))/µ (24)
is a uniform smooth approximation of the nonsmooth function

f (x) = λ 1 (Z(x))
in the sense that f µ (x) converges uniformly to f (x) as µ → 0.

Proof Following [15, Section 4], f µ is smooth in Z and

∇f µ (Z) = m i=1 e λi(Z)/µ -1 m i=1 e λi(Z)/µ q i q i ,
where q i is the ith column of the orthogonal matrix Q(Z) from the eigendecomposition of the symmetric matrix Z = Q(Z)D(Z)Q(Z) . This implies that f µ is smooth at x with the gradient given by

∇f µ (x) = Tr(∇f µ (Z(x)) Z 1 (x)) . . . Tr(∇f µ (Z(x)) Z m (x)) .
On the other hand, we have the estimate

f (x) f µ (x) f (x) + µ ln m,
which says that f µ (x) is a uniform approximation of the function f (x). Now, we can try to solve problem ( 2) and ( 12) on replacing the function f (x) = λ 1 (Z(x)) by its smooth approximation f µ (x) in (24). Due to the estimate in the above proof, to find an ε-solution x of problem (2) and ( 12), we have to find an ε 2 -solution of the smooth problem

min{f µ (x) : h(x) 0} (25)
with µ = ε 2 ln m . Here, we use a local solution of (25) to initialize the nonsmooth Algorithm 1. The smooth problem (25) can be solved using standard NLP software.

Numerical experiments

In this section, we apply our approach to a variety of problems. Let us start by commenting on practical ways to implement the stopping test 0 ∈ ∂ 1 F (x j , x j ) in step 2 of the algorithm. In practice, this is delegated to the inner loop. If the inner loop at x j finds a new feasible serious iterate x j+1 satisfying

|f (x j+1 ) -f (x j )| 1 + |f (x j )| < tol 1 , (26)
then we accept x j+1 as optimal. This corresponds to stopping the algorithm in step 2 of the (j +1)st outer loop. In our experiments, we have used tol 1 = 10 -8 . On the other hand, if the inner loop has difficulties finding a serious step and provides three unsuccessful trial steps satisfying

x j -y k 1 + x j < tol 2 , (27)
then we interpret this in the sense that x j is already optimal. This corresponds to stopping the algorithm in step 2 of the jth outer loop. Here, we have used tol 2 = 10 -7 . Theoretically, both tests are based on the observation that 0 ∈ ∂ 1 F (x j , x j ) if and only if y k = x j is solution of the tangent program in the trial step generation (see [START_REF] Gabarrou | Design of a flight control architecture using a non-convex bundle method[END_REF] for theoretical results).

In general, our stopping strategy is similar to recommendations in smooth optimization, see e.g., [10, Chapter 7], where the goal is to obtain scale independent choices of the tolerances tol 1 and tol 2 . Nonetheless, one has to accept that a nonsmooth algorithm converges very slowly at the final stages, which makes stopping a delicate task.

Before applying Algorithm 1 to solve examples of (2), note that internal stability is not a constraint in the usual sense of mathematical programming since the set S = {x ∈ R n : G(x) internally stable} is open. The stability of the system can be formulated as a constraint α (A(x))

-ε using the spectral abscissa α(A) = max{Re(λ) : λ eigenvalue of A} in the continuous time case, and as ρ (A(x)) 1 -ε using the spectral radius ρ(A) = max{|λ| : λ eigenvalue of A} in the discrete time case, for ε > 0 some small threshold. Theoretical properties of the spectral abscissa and the spectral radius have been studied in [START_REF] Burke | Differential properties of the spectral abscissa and the spectral radius for analytic matrix-valued mappings[END_REF]. In general, before optimization can start, one has, indeed, to find a stabilizing x. Using the method in [START_REF] Bompart | Non-smooth techniques for stabilizing linear systems[END_REF], this can be achieved by an initial phase where α (A(x)) is minimized until an iterate x 1 with α A(x 1 ) -ε is found.

Hankel feedback synthesis

We introduce an application of program ( 14) to a classical 1-DOF control system design, using an example from [START_REF] Boyd | Linear controller design: limits of performance[END_REF]Section 2.4]. The open-loop system G, exogenous input w and regulated output z, are given by

G = 10 -s s 2 (10 + s) , w =   d n y r   , z = y p u .
The corresponding plant is

P :   A B 1 B 2 C 1 0 D 12 C 2 D 21 0   ,
where

A =   -10 0 0 1 0 0 0 1 0   B 1 =   1 0 0 0 0 0 0 0 0   B 2 =   1 0 0   C 1 = 0 -1 10 0 0 0 D 12 = 0 1 C 2 = 0 1 -10 D 21 = 0 -1 1 .
Inspired by a manually tuned controller

Hankel system reduction

In this section, we solve program [START_REF] Nesterov | Smoothing technique and its applications in semidefinite optimization[END_REF] with the usual Hankel norm, where our tests use the 15th order Rolls-Royce Spey gas turbine engine model described in [23, Chapter 11], with data available for download on I. Postlethwaites's homepage as aero0.mat. The goal of this study is to use the theoretical values to perform a blind testing of our algorithm. For k = 1, 2, . . . , 14, using Algorithm 1, we computed Hankel reduced-order systems G k of order k, and compared the achieved objective f (x * ) = G-G k (x * ) H of ( 15) with the theoretically known optimal Hankel norm approximation errors G-G k H = σ k+1 , the (k + 1)st Hankel singular value of G. As can be seen in columns 2 and 3 of Table 1, this error is within the limits of numerical precision. In each run, the algorithm was started from a random initial guess, and no information as to the specific structure of problem (15) was provided. On average, the algorithm needed about 103 serious steps to reach the optimal objective function value within a tolerance of < 10 -10 . See Table 1 for number of iterations and running times in seconds.

Remark 7

The results show no clear relation between running times and the order of the reduced system, as one might have expected. This is due to the fact that local optimization techniques depend very sensibly on the initial guess, which in this comparison was chosen randomly.

Remark 8

In [START_REF] Dao | Minimizing the memory of a system[END_REF], we have used the same example to give a comparison between Hankel system reduction and H ∞ -system reduction, which is compared to the H ∞ -bound (see [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞ -error bounds[END_REF]).

Maximizing the memory of a system

We use here an illustrative example for [START_REF] Overton | Large-scale optimization of eigenvalues[END_REF], where G and G ref are defined as

G(s) = 1 s -1 , G ref = 11.11 s 2 + 6s + 11.11
.

The filter F is chosen of order 2,

F (s) = as 2 + bs + c s 2 + ds + e =   -d -e 1 1 0 0 b -ad c -ae a   ,
which leads to 5 tunable parameters, whereas K is a PID

K(s) = k p + k i s + k d s T f s + 1 =    0 0 k i 0 -1 T f -k d T 2 f 1 1 k p + k d T f    ,
adding another 4 unknowns. We have added a low-pass filter W 1 (s) = 0.25s+0.6 s+0.006 to the output z 1 to asses the tracking error y-y ref in low-frequency, and a highpass filter W 2 (s) = s s+0.001 on the control output z 2 to reduce high-frequency components of the control signal u + v.

Due to the choice of the performance channel w → z = (W 1 z 1 , W 2 z 2 ), the closed-loop has a non-vanishing direct transmission term. We therefore solve problem [START_REF] Mifflin | A modification and extension of Lemaréchal's algorithm for nonsmooth minimization. Nondifferential and variational techniques in optimization[END_REF] for the setup (18) using the extended Hankel program (2) with (10), and also using the constraint program [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞ -error bounds[END_REF]. Running Algorithm 1 from the same starting point, these two methods give Hankel controllers (F eH , K eH ) and (F cH , K cH ) with F eH (s) = -3.4778s 2 -13.9996s-0.0546 s 2 +1.9202s+0.0001 , K eH (s) = 6.3078 + 3.6689 s -1.0924 0.4739s+1 , F cH (s) = -3.6552s 2 -13.6987s-0.0522 s 2 +1.9588s+0.0001 , K cH (s) = 6.1959 + 3.8435 s -0.7121 0.3644s+1 , where we used the constraint σ 1 (D) η with η = 1. For comparison, we also synthesized the usual Hankel norm controller, where the direct transmission is ignored, and the H ∞ -controller, both with the same architecture:

F H (s) = -2.2376s 2 -1.9738s-2.4161 s 2 +0.9054s+0.9836 , K H (s) = 2.4482 + 0.7883 s + 0.8023 0.7817s+1 , F ∞ (s) = -9.9366s 2 -1.5077s-0.0349 s 2 +0.9969s+0.0273
, K ∞ (s) = 11.5131 + 0.2673 s -0.5507 1.0117s+1 . Figure 5 compares step responses y and step reference responses y ref for these controllers. The evolution of the optimization method for the three Hankel controllers can be traced in Fig. 6. The achieved Hankel norms are

T w→z (F eH , K eH ) H = 0.8767 < T w→z (F cH , K cH ) H = 0.8862 < T w→z (F H , K H ) H = 1.0160 < T w→z (F ∞ , K ∞ ) H = 1.0277.
This example is again interesting in so far as the parametrization of F and K arises naturally in the frequency domain.

Control of flow in a graph

Here, we give an application of program [START_REF] Rainer | Smooth roots of hyperbolic polynomials with definable coefficients[END_REF]. Let V stay = {1, 2, . . . , n x }, V in = {1, 2, . . . , m}, V out = {1, 2, . . . , p}. Let x regroup the unknown tunable parameters a jj , b ij and set A 10) and [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞ -error bounds[END_REF] give comparable results, minimization of (A, B, C) H alone (left) gives a large direct transmission Remark that the linear constraint conditions in [START_REF] Rainer | Smooth roots of hyperbolic polynomials with definable coefficients[END_REF] can be transferred to the form A eq x = b eq , x 0, which are added in each trial step generation of Algorithm 1.

(x) = [a jj ] nx×nx , B(x) = [b ij ] m×nx , C = [c 1 , . . . , c nx ], where a jj = 0 if (j, j ) ∈ A , b ij = 0 if (i, j) ∈ A . We have a discrete linear time-invariant system G(x) : x(t + 1) = A(x)x(t) + B(x)w(t) z(t) = Cx(t).
We now take the following graph G = (V , A ) with n x = 24, m = 4 and p = 4. We next consider an example using the second graph with n x = 36, m = 2 and p = 2. Let z(t) quantify the number of people on the fairground, where the 6 central nodes are counted twice. In this example, we will directly compare our nonsmooth method to the heuristic in Sect. 9. Optimization starts again at the uniform distribution x 1 . Minimizing smooth function f µ (x) in (24) with initial x 1 leads to x † , where f (x † ) = 21.7291, G(x † ) H = 4.6614, while f (x 1 ) = 578.6875, G(x 1 ) H = 24.0559. We now use x † to initialize the nonsmooth Algorithm 1. After 44 serious steps with 168 s CPU, our algorithm returns the optimal x * with f (x * ) = 14.8353, meaning G(x * ) H = 3.8517.

For the two displayed graphs, Figs. 7 and 8 compare ringing effects in unit step and white noise responses truncated at T = 30 for the three systems G(x 1 ), G(x † ) and G(x * ). We can see that ringing for G(x † ) and G(x * ) is substantially reduced.

Tables 2 and3 show a simulated study, where we compare the effects of the transition probability distributions x 1 , x † , x * by recording the evacuation of people from the fairground. We simulate crowd entering through the gates 1, . . . , 4 for different scenarios w. We then close the entrance gates at time T = 15, when in the first study 6994 people have entered the ground, and our tests w 1 is a step signal, w 2 is a sine wave, and w 3 is a square wave. A similar approach is chosen in the second graph. Column z 1 (T ) gives the number of people still present on the fairground at time T when distribution x 1 is used, and column G(x 1 ) gives the time which then elapses until this crowd is reduced below 10% of the total number 6994. Columns 5-8 are analogous. As compared to x 1 , the optimal strategy x * reduces the evacuation time to close to 1/5 in the first graph, and to close to 1/3 in the second graph.

Conclusion

We have proposed a new methodology to reduce unwanted ringing effects in a tunable linear time-invariant system. The problem was addressed by minimizing the Hankel norm of the system, a problem which leads to an eigenvalue optimization program for the associated Hankel operator. A proximal bundle algorithm was presented to solve a variety of test problems successfully, and a smooth heuristic, based on work of Nesterov [START_REF] Nesterov | Smoothing technique and its applications in semidefinite optimization[END_REF], was added and used to initialize the algorithm with a favorable initial seed.

  of reduced order k < n x with input-output behavior as close as possible to the original system G. If the model matching error e = (G -G k )w is measured in the Hankel norm, then the program minimize

τ1, put k = 1 Fig. 1

 11 Fig. 1 Flowchart of proximal bundle algorithm. Inner loop is shown in the lower right box

  K b = 219.6s 2 + 1973.95s + 724.5 s 3 + 19.15s 2 + 105.83s + 965.95 , proposed in [5, Section 2.4], we compute the optimal Hankel controller K H with the same proposed structure and compare it to K b and also to the optimal H ∞ -controller K ∞ of that same structure K(x) = as 2 + bs + c s 3 + ms 2 + ns + p = where x = [m, n, p, a, b, c] regroups the unknown tunable parameters. Using the Matlab function hinfstruct based on [1], we obtain K ∞ = 7941.9s

Fig. 2

 2 Fig.2Hankel feedback synthesis. Bearing of the algorithm. Top left shows j → f (x j ) and j → α(A(x j )) + 10 -8 . Top right shows j → x jx j+1 . Lower left shows j → k j , lower right shows j → τ j , the evolution of the memory control parameter at serious steps

Fig. 3 Fig. 4

 34 Fig. 3 Hankel feedback synthesis. Step responses (left), impulse responses (middle), magnitude plot (right) for controllers K b (dotted), K∞ (dashed), and K H (solid)

Fig. 5 Fig. 6

 56 Fig. 5 Maximizing memory. Comparison between step responses y and y ref for H∞controller and Hankel controllers computed by programs (2) with monitoring (dotted), (12) (dashed) and (2) with (10) (solid)

  Let z(t) be the total number of people on the fairground, which corresponds to the weights c 1 = • • • = c nx = 1. We start Algorithm 1 at the uniform distribution x 1 , where f (x 1 ) = 714.8634, and G(x

Fig. 8

 8 Fig.8Ringing effects of three systems G(x 1 ) (dotted), G(x † ) (dashed) and G(x * ) (solid) for the second graph. Input: Unit step signal (top) and white noise signal (bottom)

3

  Representation of the Hankel normA representation of the Hankel norm • H amenable to computations is obtained through the observability and controllability Gramians, defined in[START_REF] Zhou | Robust and optimal control[END_REF] Section 3.8]. Based on the results in [12, Section 2.3], see also[START_REF] Zhou | Robust and optimal control[END_REF] Theorem

8.1], we have the following

Lemma 1 Let G = (A, B, C) be an internally stable linear time-invariant system with input w and output z, and let

  The interest in this example is also to show that parametrizations x may arise naturally in the frequency domain. Note also that the closed-loop has no direct transmission term since D 11 = 0 and K is strictly proper. To compute K H , ) = 10.8419, meaning T w→z (P, K H ) H = 3.2927 andK H := K(x * ) =77.0614s 2 + 255.2324s + 74.6195 s 3 + 188.0709s 2 + 133.9333s + 22.2401 .

	with f (x
	2 + 13028.4s + 3611.6 s 3 + 3206.2s 2 + 12528.3s + 11078.3	.
	we solve (14) with the standard Hankel norm (1) and start Algorithm 1 at an
	initial stabilizing controller
	x 1 = [2.1460, 12.7448, 7.4208, 1.2271, 1.8013, 0.3517]
	with f (x 1 ) = 455.2874, using the stability constraint h(x) = α(A(x)) + ε 0
	with a typical value ε = 10 -8 . The stopping tests were (26) and (27). The
	algorithm came to a halt due to (26) and returned the optimal solution
	x * = [77.0614, 255.2324, 74.6195, 188.0709, 133.9333, 22.2401]

*

Table 1

 1 Hankel system reduction. Comparison of optimal values G -G k (x * ) H with theoretical values σ k+1

	k	σ k+1	G -G red H	No of iterations Time
	1	4.046418	4.046418	26	3.5
	2	2.754623	2.754624	71	21.0
	3	1.763527	1.763529	124	47.3
	4	1.296531	1.299542	151	101.5
	5	0.629640	0.629640	88	118.0
	6	0.166886	0.166887	183	197.3
	7	0.093407	0.093408	93	185.8
	8	0.022193	0.022201	76	132.4
	9	0.015669	0.015675	162	203.7
	10	0.013621	0.013624	175	191.3
	11	0.003997	0.003997	140	380.0
	12	0.001179	0.001179	57	488.4
	13	0.000324	0.000324	24	224.2
	14	0.000033	0.000033	68	372.5

  1 ) H = 26.7369. After 2469 serious iterates with 8768 s CPU, our algorithm returns the optimal x * with f (x * ) = 8.6056, meaning G(x * ) H = 2.9335. For comparison, with the Matlab function fmincon started at x 1 , we obtain x † with f (x † ) = 12.5994 > f (x * ) = 8.6056. However, if we take x † as initial for Algorithm 1, the result is f (x * ) = 8.6056, meaning G(x * ) H = 2.9335, which is achieved very fast (29 serious iterates, 87 s CPU).

Acknowledgements The authors acknowledge helpful discussions with Dr. Armin Rainer (University of Vienna).

record the time which passes until 90% of the crowd has been evacuated. In