Pierre Apkarian

Dominikus Noll

Laleh Ravanbod

Jon Borwein

Non-smooth optimization for robust control of infinite-dimensional systems *

Keywords: Robust control, non-smooth optimization, non-smooth trust-region method, infinite-dimensional system, boundary control

We use a non-smooth trust-region method for H ∞ -control of infinite-dimensional systems. Our method applies in particular to distributed and boundary control of partial differential equations. It is computationally attractive as it avoids the use of system reduction or identification. For illustration the method is applied to control a reaction-convection-diffusion system, a Van de Vusse reactor, and to a cavity flow control problem.

Introduction

Feedback control of partial differential equations and other infinite-dimensional systems encounters limitations due to computational issues. In state-space, PDE models are overly complex and not directly suited for controller synthesis. System reduction is required to bring the state-space down to a size where synthesis methods are applicable. Not only is this technically demanding, it also bears the risk of producing inaccurate or oversimplified models, where H ∞ -performance can no longer be guaranteed.

Computing the system transfer function directly from the infinite dimensional model avoids this loss of information, but encounters a second difficulty. Customary strategies now try to fit a finite-dimensional state-space model to the infinite-dimensional transfer function. This uses optimization-based identification techniques, which are in conflict with the H ∞ -objective, as the two optimization procedures in series are no longer meaningful. In addition, for unstable systems the identification often uses heuristics or ad hoc approaches, which have no certificates.

The method we propose here avoids both pitfalls. We synthesize controllers directly from the pre-computed frequency response, thereby avoiding system reduction and identification. Discretization for computation is performed in frequency space on a lowdimensional object, which avoids the loss of information. Our tests demonstrate that this works fast and reliably, once the transfer function is available. It turns out that the success of our method hinges on the use of non-smooth optimization. We use a non-smooth trust-region method first proposed in [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF], which allows trial steps tailored to the specific application. We prove convergence under Kiwiel's aggregation rule, a question which had 1 remained open in [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF]. This given an affirmative answers to a question already posed in [START_REF] Ruszczyński | Nonlinear Optimization[END_REF] for the convex non-smooth trust-region method. For complementary information on bundle methods see [START_REF]Composite proximal bundle method[END_REF][START_REF] Hare | A redistributed proximal bundle method for nonconvex optimization[END_REF][START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF], a mix of bundle and trust-regions is [START_REF] Schramm | A version of the bundle idea for minimizing a nonsmooth function: conceptual ideas, convergence analysis, numerical results[END_REF].

Design of controllers in the frequency domain based on non-smooth optimization has already been performed in [START_REF] Mayne | Algorithms for the design of control systems subject to singular value inequalities[END_REF][START_REF] Mayne | Computer aided design via optimization[END_REF][START_REF] Polak | A nondifferential optimization algorithm for the design of control systems subject to singular value inequalities over the frequency range[END_REF][START_REF] Polak | On the design of linear multivariable feedback systems via constrained nondifferentiable optimization in H ∞ space[END_REF][START_REF] Kiwiel | A linearization algorithm for computing control systems subject to singular value inequalities[END_REF][START_REF] Cullum | The minimization of certain nondifferentiable sums of eigenvalues of symmetric matrices[END_REF]. Structured H -∞-control for infinitedimensional systems is addressed in [START_REF] Apkarian | Structured H ∞ -control of infinite-dimensional systems[END_REF]. H ∞ -control of a heat exchange system is discussed in [START_REF] Sano | H ∞ -control of a parallel-flow heat exchange process[END_REF]. These approaches use either unstructured controllers, are based on matrix inequalities, or differ with regard to the optimization technique.

The structure of the paper is as follows. In section 2 we outline our approach to H ∞ -control of infinite dimensional systems. In section 3 we discuss optimization and present our non-smooth trust-region method originally proposed in [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF], on which the present approach rests. Convergence of the non-smooth trust-region method is discussed in section 3.2. In section 4 we point to some particularities when applying the trust-region method to H ∞ -optimization. This concerns the choice of working model, trial step, and stability barrier, as well as the approximation error between the infinite-dimensional H ∞program and its discretization.

Numerical results for applications to infinite-dimensional control problems are presented in section 5. Subsection 5.1 shows how the method is, in general, applied to a boundary control problem, subsection 5.2 illustrates this in boundary H ∞ -control of a non-linear reaction-convection-diffusion equation, and subsection 5.3 for a non-linear Van de Vusse reactor. Subsection 5.4 discusses a cavity flow control problem.

Control strategy

We consider an abstract linear time-invariant control system of the form G :

{ ẋ = Ax + Bu y = Cx + Du ( 1
)
where A is an unbounded linear operator on a Hilbert space Z generating a strongly continuous semi-group, B a closed linear operator mapping the control input space U to Z, C a closed linear operator mapping Z to the space Y of measured outputs, and D a closed linear operator mapping U to Y . For practical reasons we assume that U ≃ R p and Y ≃ R m , which reflects the fact that the process is assessed by a finite number of sensors and actuators, rendering our control law physically implementable. This means that B, D are bounded, while C is allowed to be closed unbounded. In addition, the domains of A and C satisfy D(A) ⊂ D(C). As a consequence of the finite rank assumption on Y, U , the transfer function G(s) = C(sI -A) -1 B + D is defined on the resolvent set ρ(A) and meromorphic on C, with values G(s) ∈ C p × C m , see e.g. [START_REF] Jacobson | Linear state-space systems in infinite-dimensional space: the role and characterization of joint stabilizability/detectability[END_REF][START_REF] Engel | One-parameter semigroups for linear evolution equations[END_REF]. We consider a class K ∈ K of feedback control laws which have similar state-space realizations

K : { ẋK = A K x K + B K y u = C K x K + D K y (2)
on a Hilbert space Z K , with input space Y ≃ R m and output space U ≃ R p , so that we can put G and K ∈ K in lower feedback F ℓ (G, K) as in Figure 1. Candidate controllers K ∈ K have to stabilize F ℓ (G, K) internally in closed loop, by which we mean that the infinitesimal generator of the closed loop system generates an exponentially stable semi-group [START_REF] Jacobson | Linear state-space systems in infinite-dimensional space: the role and characterization of joint stabilizability/detectability[END_REF][START_REF] Curtain | An introduction to infinite-dimensional linear system theory[END_REF]. In H ∞ -control one has not only to assure stability in closed loop, but also to good performance and robustness of the feedback system. For that purpose G is embedded in a plant P with a similar state-space realization P :

⎧ ⎨ ⎩ ẋ = Ax + B 1 w + B 2 u z = C 1 x + D 11 w + D 12 u y = C 2 x + D 21 w + D 22 u
where the channel y → u is used for control, the channel w → z for performan the u-y-loop in (3) with (2) as in Figure 1 leaves us with the closed-loop transf T wz (K) from exogenous input w to regulated output z. The H ∞ -control proble now in minimizing the L 2 -L 2 -operator norm of T wz (K) over a suitable class admissible control laws. This operator norm is also known as the H ∞ -norm, g

∥T wz (K)∥ ∞ = max ω∈[0,∞] σ (T wz (K, jω)) ,
where σ(M ) denotes the maximum singular value of a matrix M . If candidate K ∈ K in (2) are parametrized as K(x) for a finite-dimensional vector x ∈ R n parameters, then this H ∞ -optimization program takes the form

minimize f ∞ (x) = max ω∈[0,∞] σ (T wz (K(x), jω)) subject to K(x) stabilizes G in closed-loop x ∈ R n
For finite-dimensional real-rational P, K(x) the objective (4) may be compu iterative procedure [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF][START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF], but for infinite dimensional G, K(x), this is no longe and we need to approximate (4) on a finite grid

Ω opt ⊂ [0, ∞]. Introducing the version ∥T wz (K)∥ ∞,d = max ω∈Ω opt σ (T wz (K, jω))
of the H ∞ -norm on the grid Ω opt , this leads to the discretized

H ∞ -optimizatio minimize f (x) = max ω∈Ω opt σ (T wz (K(x), jω)) subject to K(x) stabilizes G in closed-loop Figure 1: Lower feedback interconnection F ℓ (P, K) with H ∞ -performance channel w → z.
In H ∞ -control one has not only to assure stability in closed loop, but also to guarantee good performance and robustness of the feedback system. For that purpose the system G is embedded in a plant P with a similar state-space realization

P :    ẋ = Ax + B 1 w + B 2 u z = C 1 x + D 11 w + D 12 u y = C 2 x + D 21 w + D 22 u (3)
where the channel y → u is used for control, the channel w → z for performance. Closing the u-y-loop in (3) with (2) as in Figure 1 leaves us with the closed-loop transfer function T wz (K) from exogenous input w to regulated output z. The H ∞ -control problem consists now in minimizing the L 2 -L 2 -operator norm of T wz (K) over a suitable class K ∈ K of admissible control laws. This operator norm is also known as the H ∞ -norm, given as

∥T wz (K)∥ ∞ = max ω∈[0,∞] σ (T wz (K, jω)) , ( 4
)
where σ(M ) denotes the maximum singular value of a matrix M . If candidate controllers K ∈ K in (2) are parametrized as K(x) for a finite-dimensional vector x ∈ R n of tunable parameters, then this H ∞ -optimization program takes the form

minimize f ∞ (x) = max ω∈[0,∞] σ (T wz (K(x), jω)) subject to K(x) stabilizes G in closed-loop x ∈ R n (5)
For finite-dimensional real-rational P, K(x) the objective (6) may be computed by an iterative procedure [START_REF] Boyd | A bisection method for computing the H ∞norm of a transfer matrix and related problems[END_REF][START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF], but for infinite dimensional G, K(x), this is no longer possible, and we need to approximate (6) on a finite grid

Ω opt ⊂ [0, ∞]. Introducing the discretized version ∥T wz (K)∥ ∞,d = max ω∈Ωopt σ (T wz (K, jω))
of the H ∞ -norm on the grid Ω opt , this leads to the discretized

H ∞ -optimization program minimize f (x) = max ω∈Ωopt σ (T wz (K(x), jω)) subject to K(x) stabilizes G in closed-loop x ∈ R n (6)
to which our non-smooth optimization method is applied. The overall procedure for H ∞ -control is now given in algorithm 1.

Algorithm 1. Infinite-dimensional H ∞ -synthesis

Parameters: Tolerance ϑ > 0.

▷ Step 1 (Steady-state). Compute steady-state and linearize infinite-dimensional system about steady-state.

▷

Step 2 (Transfer function). Use linearized infinite-dimensional system to compute transfer function G(jω), either formally for arbitrary ω ∈ [0, ∞], or numerically at a very high precision on a very fine grid Ω fine .

▷

Step 3 (Plant). Set up plant P which defines the H ∞ -performance channel w → z.

▷

Step 4 (Grid for optimization). Find initially stabilizing K(x 0 ) for G and use it to compute grid for optimization

Ω opt such that ∥T wz (K(x 0 ))∥ ∞ ≤ f (x 0 ) + ϑ.
Either use a formal or a numerical function G(jω), or extract Ω opt from the pre-computed high precision grid Ω fine .

▷

Step 5 (Non-smooth optimization). Use non-smooth trust-region algorithm 2 to compute locally optimal solution K(x * ) of ( 6).

▷ Step 6 (Refined grid). Check whether ∥T wz (K(x * ))∥ ∞ ≤ ∥T wz (K(x * ))∥ ∞,d + ϑ.
If not then add nodes to Ω opt and go back to step 5.

Remark 1. It should be stressed that the computation of G(jω) in steps 1 and 2 of algorithm 1 is the only moment where the full infinite-dimensional, or likewise, large-scale finite-dimensional model is used. Since this step is performed prior to optimization, and typically |Ω opt | ≪ |Ω fine |, optimization is speedy. Also, the process of finding the correct P , which needs going back to step 3 of algorithm 1, is not slowed down by steps 1 and 2. In other words, once the transfer function G(jω) is available, the complexity of the process is the same as that of a finite-dimensional structured H ∞ -design procedure in the sense of [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF][START_REF] Boyd | Linear controller design. Limits of performance[END_REF]. This will be illustrated in our experimental section.

Remark 2. Stability in closed-loop in step 4, and during optimization in step 5, uses the well-known Nyquist stability test in tandem with the barrier approach to be addressed in section 4.1. For the theoretical justification in the context of infinite-dimensional systems see [START_REF] Apkarian | Structured H ∞ -control of infinite-dimensional systems[END_REF].

We now address the individual steps of algorithm 1. The central ingredient is optimization, which is needed to solve [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF], and which is discussed in the next sections 3 and 3.2. Generation of the grid Ω opt in step 4, and the certificate in step 6, are discussed in section 4.

Non-smooth trust-region method

We work with the non-smooth trust-region method introduced in [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF], which has already been successfully used in mechanical contact problems [START_REF] Dao | Nonconvex bundle method with application to a delamination problem[END_REF], and in system theory [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF] for computing the worst-case H ∞ -performance of a system, its stability margin, and its distance to instability. Here we use it for H ∞ -control of infinite-dimensional systems. We also answer a question left open in [5, Remark 16], which concerns the theoretical justification of Kiwiel's aggregation technique [START_REF] Kiwiel | An aggregate subgradient method for nonsmooth convex minimization[END_REF] in non-smooth trust-regions. This question goes back to [START_REF] Ruszczyński | Nonlinear Optimization[END_REF] for the convex trust-region method, but had until now remained open. An affirmative answer will be obtained in section 3.2.

Presentation of the algorithm

We briefly recall the essentials of the non-smooth trust-region method. For the present work it is sufficient to apply it to optimization programs of the form

min x∈R n f (x), (7)
where f : R n → R is locally Lipschitz but non-smooth and non-convex. Following [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF], a function ϕ : R n × R n → R is called a model of f if it satisfies the following properties:

(M 1 ) ϕ(•, x) is convex, ϕ(x, x) = f (x), and ∂ 1 ϕ(x, x) ⊂ ∂f (x). (M 2 ) If y k → x, then there exist ϵ k → 0 + such that f (y k ) ≤ ϕ(y k , x) + ϵ k ∥y k -x∥. (M 3 ) If x k → x, y k → y, then lim sup k→∞ ϕ(y k , x k ) ≤ ϕ(y, x).
We may interpret ϕ(•, x) as a substitute for the first-order Taylor expansion of f at x.

For convergence theory we need a slightly stronger type of model, which is given by the following:

Definition 1. A first-order model ϕ of f is called strict if it satisfies the following stronger version of axiom (M 2 ): ( M 2 ) If x k , y k → x, then there exist ϵ k → 0 + such that f (y k ) ≤ ϕ(y k , x k ) + ϵ k ∥y k -x k ∥.

□

The difference between (M 2 ) and the strict version ( M 2 ) is analogous to the difference between differentiability and strict differentiability, hence the nomenclature. For additional information on the model concept see [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF][START_REF] Apkarian | A proximity control algorithm to minimize non-smooth and non-convex semi-infinite maximum eigenvalue functions[END_REF][START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF].

Remark 3. A typical example of a strict model ϕ is obtained when f is a maximum eigenvalue function f (x) = λ 1 (F (x)), with F : R n → S m a class C 1 -mapping into the space of m × m Hermitian matrices. We take ϕ(y, x) = λ 1 (F (x) + F ′ (x)(y -x)). See [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF][START_REF] Apkarian | A proximity control algorithm to minimize non-smooth and non-convex semi-infinite maximum eigenvalue functions[END_REF].

Definition 2. Let x be the current serious iterate of the trust-region algorithm, z a trial step. Let g be a subgradient of ϕ(•, x) at z. Then the affine function m z (•, x) = ϕ(z, x) + g ⊤ (• -z) is called a cutting plane of f at serious iterate x and trial step z. □

If z = x, then due to axiom (M 1 ) a cutting plane m x (•, x) at serious iterate x and trial step z = x is just a tangent plane to f at x. Since m x (x, x) = f (x), including m x (•, x) in the working model ϕ k (•, x) at x guarantees exactness ϕ k (x, x) = f (x) of the working model at all counters k. We cast this in the following Definition 3. A cutting plane m x (•, x) at serious iterate x and trial step z = x is called an exactness plane. □

As is standard in bundle and cutting plane algorithms, by storing cutting planes at unsuccessful trial steps z k , we accumulate information, which we use to build polyhedral models of f near x. We use the notation ϕ k (•, x) for these working models of f formed by cutting planes, where k denotes the counter of the inner loop. Note that ϕ k ≤ ϕ by construction of the cutting planes. If in addition a positive semi-definite symmetric matrix Q(x) ⪰ 0 is available as a substitute for the Hessian of f at x, then we call

Φ k (•, x) = ϕ k (•, x) + 1 2 (• -x) ⊤ Q(x)(• -x)
a second-order working model of f at serious iterate x.

We are now ready to present the bundle trust-region algorithm. (See algorithm 2).

Algorithm 2. Non-smooth trust-region algorithm

Parameters: 0 < γ < γ < 1, 0 < γ < Γ ≤ 1, 0 < θ ≪ 1, M ≥ 1, q > 0.
▷ Step 1 (Initialize outer loop). Fix initial iterate x 1 and memory trust-region radius R ♯ 1 > 0. Initialize Q 1 ⪰ 0 with ∥Q 1 ∥ ≤ q. Put outer loop counter j = 1. ▷ Step 2 (Stopping test). At outer loop counter j, stop if x j is a critical point of [START_REF] Boyd | A bisection method for computing the H ∞norm of a transfer matrix and related problems[END_REF]. Otherwise go to inner loop.

▷

Step 3 (Initialize inner loop). Put inner loop counter k = 1 and initialize trustregion radius as R 1 = R ♯ j . Build polyhedral first-order working model ϕ 1 (•, x j ), where at least one exactness plane at x j is included. Possibly enrich by adding recycled planes from previous steps, or by including anticipated cutting planes. Build secondorder working model

Φ 1 (•, x j ) = ϕ 1 (•, x j ) + 1 2 (• -x j ) ⊤ Q j (• -x j ). ▷ Step 4 (Trial step generation). At inner loop counter k compute solution y k of trust-region tangent program minimize Φ k (y, x j ) subject to ∥y -x j ∥ ≤ R k Then admit any z k satisfying ∥z k -x j ∥ ≤ M ∥y k -x j ∥ and f (x j ) -Φ k (z k , x j ) ≥ θ ( f (x j ) -Φ k (y k , x j ) ) as trial step. ▷ Step 5 (Acceptance test). If ρ k = f (x j ) -f (z k ) f (x j ) -Φ k (z k , x j ) ≥ γ put x j+1 = z k (serious
Then Φ k+1 (•, x j ) = ϕ k+1 (•, x j ) + 1 2 (• -x j ) ⊤ Q j (• -x j )
is the new second-order working model. Continue with step 7.

▷

Step 7 (Update trust-region radius). Compute secondary control parameter

ρ k = f (x j ) -ϕ k+1 (z k , x j ) f (x j ) -Φ k (z k , x j )
and put

R k+1 = { R k if ρ k < γ 1 2 R k if ρ k ≥ γ
R ♯ j+1 = { R k if ρ k < Γ 2R k if ρ k ≥ Γ Update Q j → Q j+1 respecting Q j+1 ⪰ 0 and ∥Q j+1 ∥ ≤ q.
Increase outer loop counter j and go back to step 2.

Remark 4. Before we discuss convergence of algorithm 2 in the next section, we recall the form of the tangent program in step 4 from [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF]. Let the first-order working model at inner loop counter k have the form

ϕ k (•, x j ) = max i∈I k a i + g ⊤ i (• -x j )
for some finite set I k , and suppose the trust-region norm is the maximum norm. Then the tangent program at serious iterate x j and inner loop instant k is the following CQP

minimize t + 1 2 (y -x j ) ⊤ Q j (y -x j ) subject to a i + g ⊤ i (y -x j ) ≤ t, i ∈ I k -R k ≤ y i -x j i ≤ R k , i = 1, . . . , n (8)
with decision variable (t, y) ∈ R × R n , giving rise to the solution y k in step 4.

Convergence

In this section we prove convergence of the trust-region algorithm 2 toward a Clarke critical point. As is standard, we start by proving that the inner loop ends finitely if 0 ̸ ∈ ∂f (x j ), where throughout ∂f denotes the Clarke subdifferential. During this part of the proof we write x = x j and Q = Q j , as those are fixed during the inner loop at counter j. Note that by the necessary optimality condition for the tangent program in step 4 of algorithm 2, there exists a subgradient

g * k ∈ ∂ 1 ϕ k (y k , x) such that g * k +Q(y k -x)+v k = 0, where v k is in the normal cone to the trust-region norm ball B(x, R k ) at y k . Definition 4. We call g * k the aggregate subgradient. The affine function m * k (•, x) = ϕ k (y k , x) + g * ⊤ k (• -y k ) is called the aggregate plane. □ Lemma 1.
There exists σ > 0 depending only on the constants θ ∈ (0, 1) and M > 0 in algorithm 2 and on the trust-region norm ∥ • ∥, such that for every trial point z k at inner loop instant k with corresponding solution y k of the trust-region tangent program in step 4, and for the corresponding aggregate subgradient g * k ∈ ∂ 1 ϕ(y k , x), we have the estimate

f (x) -ϕ k (z k , x) ≥ σ∥g * k + Q(y k -x)∥∥z k -x∥. (9)
Proof: This is essentially the same as [5, Lemma 1]. □ Lemma 2. Suppose the inner loop at x turns infinitely, and lim inf k→∞ R k = 0. Then x is a critical point of [START_REF] Boyd | A bisection method for computing the H ∞norm of a transfer matrix and related problems[END_REF].

Proof: According to step 7 of algorithm 2 we have

ρ k ≥ γ for infinitely many k ∈ K. Since R k is never increased during the inner loop, that implies R k → 0. Hence y k , z k → x as k → ∞
, where we use the trial step generation rule of step 4 of algorithm 2. We argue that this implies

ϕ k (z k , x) → f (x). Indeed, lim sup k→∞ ϕ k (z k , x) ≤ lim sup k→∞ ϕ(z k , x) = lim k→∞ ϕ(z k , x) = f (x) is al- ways true due to ϕ k ≤ ϕ and axiom (M 1 ). On the other hand, ϕ k includes (i.e. domi- nates) an exactness plane m 0 (•, x) = f (x) + g ⊤ 0 (• -x), hence f (x) = lim k→∞ m 0 (z k , x) ≤ lim inf ϕ k (z k , x). The two together show ϕ k (z k , x) → f (x)
, and then immediately also Φ k (z k , x) → f (x). We also readily obtain ϕ k (y k , x) → f (x) from the link between z k , y k in step 4 of algorithm 2.

We now prove that lim inf k→∞ ∥g * k ∥ = 0. Assume on the contrary that ∥g

* k ∥ ≥ η > 0 for all k. Choose k large enough to have ∥g * k + Q(y k -x)∥ ≥ 1 2 ∥g * k ∥
. This is possible because the sequence g * k is bounded away from 0 and Q(y k -x) → 0. Then by estimate ( 9)

we have f (x)-ϕ k (z k , x) ≥ 1 2 ση∥z k -x∥. Therefore, for k large enough, f (x)-Φ k (z k , x) ≥ 1 4 ση∥z k -x∥, as the quadratic term in Φ k is of the order ∥z k -x∥ 2 . Since z k → x, by axiom ( M 2 ) there exist ϵ k → 0 + such that f (z k ) -ϕ(z k , x) ≤ ϵ k ∥z k -x∥. Now we estimate ρ k = ρ k + f (z k ) -ϕ(z k , x) f (x) -Φ k ( z k , x) ≤ ρ k + ϵ k ∥z k -x∥ 1 4 ση∥z k -x∥ = ρ k + 4ϵ k /(ση).
Since ϵ k → 0 and ρ k < γ, we have lim sup ρ k ≤ γ < γ, a contradiction with ρ k > γ for the infinitely many k ∈ K. That proves g * k → 0 for a subsequence k ∈ N . Next observe that by the subgradient inequality and ϕ k ≤ ϕ we have

g * ⊤ k h ≤ ϕ k (y k + h, x) -ϕ k (y k , x) ≤ ϕ(y k + h, x) -ϕ k (y k , x). Since ϕ k (y k , x) → f (x) = ϕ(x, x), passing to the limit k ∈ N and using g * k → 0, y k → x implies 0 ≤ ϕ(x + h, x) -ϕ(x, x).
Since h was arbitrary, we have 0 ∈ ∂ 1 ϕ(x, x) ⊂ ∂f (x) by (M 1 ). That proves the Lemma. □

This result needs only the fact that ϕ k ≤ ϕ, so it is not in conflict with any rule in step 6 used to taper out ϕ k → ϕ k+1 , as long as some exactness plane is present in ϕ k at all inner loop instants k. In the next two lemmas we examine the more involved case when R k is bounded away from 0. Here we require not only that an exactness plane is present at all times, but also that the latest cutting plane is added into ϕ k+1 .

However, this leads to a tangent program of size growing with k, which raises the question whether it is in principle possible to limit the number of planes included in the working model ϕ k (•, x) in step 6. For the convex bundle method this question is answered in the affirmative by Kiwiel's aggregation rule [START_REF] Kiwiel | An aggregate subgradient method for nonsmooth convex minimization[END_REF], according to which only three planes are required, an exactness plane, the latest cut, and the aggregate plane to account for the past. For the non-convex bundle method, an affirmative answer was first given in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF]. The aggregate plane is a convex combination of active cuts at the trial point y k , and can be described as follows. Had we removed from the last tangent program all active planes, and substituted instead the aggregate plane, the solution y k would have been the same. The idea of the aggregation technique is to add the aggregate plane into the working model after an unsuccessful trial step y k , which allows to remove the active planes for the next sweep. Inactive planes may leave the model in any case.

Remark 5. Ruszczyński [START_REF] Ruszczyński | Nonlinear Optimization[END_REF] stresses that the situation is more delicate for the convex non-smooth trust region method, and asks whether convergence could be proved under the aggregation rule, or under any other rule allowing to limit the number of planes included in the ϕ k . Here we address this question in the general non-convex case. In [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF] we had shown that the number of cuts in step 6 may at least be limited to n + 2 using Carathéodory's theorem, but we remarked that it would be far more attractive to have a maximum number independent of the dimension n as in Kiwiel's rule. In [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF]Remark 16] we observed that the question whether Kiwiel's rule could also be justified for the trust-region method was still open. Here we shall answer this question in the affirmative by proving convergence under the aggregation rule.

The following result justifies the use of aggregation in the first place for the special case z k = y k . Note that the trivial choice z k = y k in step 4 is always authorized (due to M ≥ 1 and θ ≤ 1), but of course we want to use the additional freedom offered by z k to improve performance of our method, so z k = y k is rather restrictive, and we will remove it later. Lemma 3. Suppose the inner loop at x turns infinitely and the trust-region radius R k stays bounded away from 0. Let Q ≻ 0 and suppose the aggregation rule is used to taper out the models in step 6. Suppose the y k are chosen as trial steps. Then x is a critical point of [START_REF] Boyd | A bisection method for computing the H ∞norm of a transfer matrix and related problems[END_REF].

Proof: Since the trust-region radius is frozen R k = R k 0 from some counter k 0 onwards, we write R := R k 0 . According to step 7 of the algorithm that means ρ k < γ for k ≥ k 0 . The only progress in the working model as we update ϕ k → ϕ k+1 is now the addition of the cutting plane and the aggregate plane. The working models ϕ k now contain at least three planes, an exactness plane, the latest cut from the last unsuccessful trial step, and the aggregate plane. They may contain more planes, but those will not be used in our argument below.

We want to prove y k → x. Since R k stays bounded away from 0, this is more involved than in the previous Lemma. Since Q ≻ 0 is fixed, we introduce the Euclidian norm

|x| 2 Q = x ⊤ Qx. With this arrangement the objective function of the tangent program becomes Φ k (•, x) = ϕ k (•, x) + 1 2 | • -x| 2 Q . We know that the cutting plane m k (•, x) at trial step y k satisfies m k (y k , x) = ϕ(y k , x), so it memorizes the value ϕ(y k , x), while the aggregate plane m * k (•, x) satisfies m * k (y k , x) = ϕ k (y k , x), so it memorizes in turn the value ϕ k (y k , x). The latter gives Φ k (y k , x) = m * k (y k , x) + 1 2 |y k -x| 2 Q . ( 10
)
Now we introduce the quadratic function

Φ * k (•, x) = m * k (•, x) + 1 2 | • -x| 2 Q ,
then from what we have just seen in ( 10)

Φ * k (y k , x) = Φ k (y k , x). (11)
Moreover, we have

Φ * k (•, x) ≤ Φ k+1 (•, x), (12)
because according to the aggregation rule we include the aggregate plane m * k (•, x) in the built of the new model ϕ k+1 , that is, we have m * k (•, x) ≤ ϕ k+1 (•, x), and hence [START_REF] Dochain | Approximation of the dynamical model of fixed bed reactors via a singular perturbation approach[END_REF]. Expanding the quadratic function

Φ * k (•, x) at y k gives Φ * k (•, x) = Φ * k (y k , x) + ∇Φ * k (y k , x) ⊤ (• -y k ) + 1 2 | • -y k | 2 Q ,
where

∇Φ * k = g * k + Q(y k -x).
From the optimality condition of the tangent program at

y k we get g * k + Q(y k -x) = -v k with v k in the normal cone to the ball B(x, R) at y k , hence Φ * k (•, x) = Φ * k (y k , x) -v ⊤ k (• -y k ) + 1 2 | • -y k | 2 Q . ( 13
)
Now we argue as follows:

Φ k (y k , x) = Φ * k (y k , x) (by (11)) ≤ Φ * k (y k , x) + 1 2 |y k+1 -y k | 2 Q = Φ * k (y k+1 , x) + v ⊤ k (y k+1 -y k ) (by (13)) ≤ Φ * k (y k+1 , x) (since v ⊤ k (y k+1 -y k ) ≤ 0) ≤ Φ k+1 (y k+1 , x) (by (12)) ≤ Φ k+1 (x, x) (y k+1 minimizer of Φ k+1 (•, x)) = ϕ(x, x) = f (x). ( 14
)
Therefore the sequence Φ k (y k , x) is increasing and bounded above, and converges to a limit Φ * ≤ f (x). Going back with this information to the estimation chain [START_REF] Fogler | Elements of Chemical Reaction Engineering, Fourth Edition[END_REF] shows

1 2 |y k+1 -y k | 2 Q → 0 and also v ⊤ k (y k+1 -y k ) → 0.
Then also

1 2 |y k+1 -x| 2 Q -1 2 |y k -x| 2 Q → 0, because | • | Q is a Euclidian norm. In consequence ϕ k+1 (y k+1 , x) -ϕ k (y k , x) = Φ k+1 (y k+1 , x) -Φ k (y k , x) -1 2 |y k+1 -x| 2 Q + 1 2 |y k -x| 2 Q → 0.
Now recall that the cutting plane m k (•, x) is an affine support function of ϕ k+1 (•, x) at y k . Hence by the subgradient inequality

g ⊤ k (• -y k ) ≤ ϕ k+1 (•, x) -ϕ k+1 (y k , x).
Since ϕ k+1 (y k , x) = ϕ(y k , x), we deduce

ϕ(y k , x) + g ⊤ k (• -y k ) ≤ ϕ k+1 (•, x). ( 15
)
Now using [START_REF] Hare | A redistributed proximal bundle method for nonconvex optimization[END_REF] we estimate as follows:

0 ≤ ϕ(y k , x) -ϕ k (y k , x) = ϕ(y k , x) + g ⊤ k (y k+1 -y k ) -ϕ k (y k , x) -g ⊤ k (y k+1 -y k ) ≤ ϕ k+1 (y k+1 , x) -ϕ k (y k , x) -g ⊤ k (y k+1 -y k ).
Since y k+1 -y k → 0 and the g k are bounded, we have

g ⊤ k (y k+1 -y k ) → 0, hence we deduce ϕ(y k , x) -ϕ k (y k , x) → 0, and also Φ(y k , x) -Φ k (y k , x) → 0. Now we claim that ϕ k (y k , x) → f (x). Since ϕ k (y k , x) ≤ Φ k (y k , x) → Φ * ≤ f (x), it remains to prove lim inf ϕ k (y k , x) ≥ f (x)
. Suppose that this is not the case, and let ϕ k (y k , x) → f (x)-η for a subsequence and some η > 0. Then also ϕ(y k , x) → f (x)-η for that subsequence. (Here we use that ϕ k (y k , x) -ϕ(y k , x) → 0 proved above). Passing to yet another subsequence, and using boundedness of the y k , we may assume 1 2

|y k -x| 2 Q → ℓ ≥ 0. Choose δ > 0 such that δ < (1 -γ)η. From what we have just seen there exists k 1 such that ϕ(y k , x) -ϕ k (y k , x) < δ for all k ≥ k 1 . Now recall that ρ k ≤ γ for every k ≥ k 0 , hence γ ( Φ k (y k , x) -f (x) ) ≤ ϕ(y k , x) -f (x) ≤ ϕ k (y k , x) -f (x) + δ.
Passing to the limit gives -γη

+ γℓ ≤ -η + δ, hence (1 -γ)η + ℓ γ ≤ δ, which contradicts the choice of δ. Hence ϕ k (y k , x) → f (x). We immediately deduce that Φ k (y k , x) → f (x) and Φ(y k , x) → f (x).
We now argue that y k → x. This follows from the estimates

ϕ k (y k , x) ≤ Φ k (y k , x) = ϕ k (y k , x) + 1 2 |y k -x| 2 Q ≤ Φ * ≤ f (x) because ϕ k (y k , x) → f (x) now
shows that all terms go to f (x), and that implies |y k -x| Q → 0. Since Q ≻ 0 we deduce y k → x, and this is where the proof no longer works if only Q ⪰ 0. Note that this shows that y k is in the interior of B(x, R) from some counter onward, so that v k = 0.

Let us now show that 0 ∈ ∂f (x). From the subgradient inequality we have

g * ⊤ k (x + h -y k ) ≤ ϕ k (x + h, x) -ϕ k (y k , x) ≤ ϕ(x + h, x) -ϕ k (y k , x). Passing to the limit using g * k → 0, ϕ k (y k , x) → f (x) = ϕ(x, x), we obtain 0 ≤ ϕ(x + h, x) -ϕ(x, x),
and since h is arbitrary and ϕ(•, x) is convex, this gives 0 ∈ ∂ 1 ϕ(x, x) ⊂ ∂f (x). □ Remark 6. For the general case Q ⪰ 0 and z k different from y k it is still not known whether aggregation is justified, but with [5, Lemma 2] we can prove convergence if we keep all cuts in the model, or if we use the Carathéodory type argument of that reference to limit the number of cutting planes in the inner loop to ≤ n + 2.

Remark 7. Note that Q ≻ 0 is not a restriction in practice, but z k = y k is. Fortunately the aggregation technique may still be justified in the case z k ̸ = y k if we proceed as follows. In the first place we allow z k as a trial point in step 4. If acceptance fails, then we perform step 7. However, if step 7 gives no reduction of R k , then we are in the difficult case. We then do the following. We fall back on y k as the trial point, i.e., we forget about z k . When y k is not accepted, we proceed with step 6 and apply aggregation. This is now justified because we are in the situation covered by Lemma 3. Note that the additional work required in steps 6 and 7 is marginal, so we do not waste time by this evasive maneuver. We could even perform this maneuver as default (i.e. checking y k if z k fails). We have proved the following Lemma 4. Suppose the inner loop turns infinitely. Suppose Q ≻ 0 and that the aggregation rule is used in step 6 to limit the size of ϕ k to any pre-defined fixed number N ≥ 3. Suppose we accept to fall back on y k if z k fails in step 5 with ρ k < γ in step 7. Then 0 ∈ ∂f (x). □ Remark 8. In the classical trust-region method failure of the trial step always leads to reduction of the trust-region radius. One occasionally sees non-smooth versions where authors do the same. As we have already shown in [5, section 5.5], that must fail. The example in that reference also shows that the Cauchy point fails in the non-smooth case.

We are now ready to state the main convergence result for our optimization method. The proof may be adapted from [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF] with minor changes, so we skip it here. Theorem 1. Suppose x 1 is such that the level set {x ∈ R n : f (x) ≤ f (x 1 )} is bounded. Let x j be the sequence of serious iterates generated by the bundle trust-region algorithm. Then every accumulation point x * of the x j is a critical point of [START_REF] Boyd | A bisection method for computing the H ∞norm of a transfer matrix and related problems[END_REF]. □

The case of H ∞ -synthesis

The general purpose algorithm 2 is readily applicable to the H ∞ -design problem [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF], as this is a special case of [START_REF] Boyd | A bisection method for computing the H ∞norm of a transfer matrix and related problems[END_REF]. Since f is the square-root of a maximum eigenvalue function, the ideal model ϕ in algorithm 2 is chosen as in remark 3. There are, however, some particularities in the application of algorithm 2 to (6), on which we comment in this section.

Stability barrier

The hidden constraint of closed-loop stability may occasionally lead to unacceptable trial points y k , but this can be avoided by complementing the objective in ( 6) by the following barrier function:

f (x) = max{∥T wz (K(x))∥ ∞,d , c∥S(K(x))∥ ∞,d }, where S = (I + GK) -1
is the closed-loop sensitivity function, and c > 0 a small constant. It is well-known that ∥S∥ -1 ∞ , also known as the modulus margin, is an indicator for the distance of the Nyquist curve to the point -1, and since the Nyquist curve crosses -1 when iterates become destabilizing, the term ∥S∥ ∞ becomes large as iterates approach the limit of the region of stability. In other words, c∥S∥ ∞ has the effect of a barrier function at the boundary of the hidden constraint. Note that the maximum of two H ∞ -norms is again a H ∞ -norm, i.e., we may still represent the modified objective as f (x) = ∥T w ′ ,z ′ (K)∥ ∞,d for a modified channel w ′ → z ′ . Note that other closed-loop transfer function, which may contribute a similar complementary barrier effect include GS = G(I + GK) -1 , (I -KG) -1 and K(I + KG) -1 , see e.g. [8, p.36].

Exploiting freedom in steps 3 and 4

An important practical aspect of program ( 6) is to use an adapted initialization of the model ϕ 1 (•, x) at the beginning of the inner loop, the idea being that in the vast majority of cases the first trial step will then be successful. This is achieved by including not only active frequencies ω from T wz (K(x), jω) at x in the model, but also branches belonging to secondary peaks, which are susceptible to become active at the next trial step. Selecting near active frequencies of the H ∞ -norm is decisive for the quality of the working models ϕ k and was already discussed in [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF] and [START_REF] Apkarian | A proximity control algorithm to minimize non-smooth and non-convex semi-infinite maximum eigenvalue functions[END_REF] in the context of the bundle method, and this is shown schematically in Figure 2. We refer to this type of affine functions as anticipated cutting planes, and their integration into the working models is covered by convergence theory as long as they are affine minorants of ϕ. Remark 9. We next exploit the freedom in step 4 of algorithm 2 in the case of program [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF]. In case of failure of the solution y k of the tangent program it is attractive to use a backtracking linesearch to generate trial steps of the form

z k = x + t(y k -x), 0 < t < 1. By convexity, Φ k (x, x) -Φ k (x + t(y k -x), x) ≥ t ( Φ k (x, x) -Φ k (y k , x) )
, hence every z k with t ≥ θ gives automatically a trial point in the sense of step 4, and we only have to check acceptance ρ k ≥ γ. Even for smaller steps t < θ it may still be possible to have

f (x j ) -Φ k (z k , x j ) ≥ θ ( f (x j ) -Φ k (y k , x j ) )
, in which case z k remains a candidate.

Performance certificate

The strategy in algorithm 1 is to perform discretization at the level of the system transfer function G(s), and not before, avoiding system reduction or identification. To justify this we emphasize the important fact that when singular values σ (G(κ, jω)) are simple on Ω e (κ), which is the rule in practice, we have Z ω = Tr Z ω so that SDP (4.21) simplifies to a much faster convex QP. A fact that can also be exploited for computing minimal norm elements.

Nonsmooth loop-shaping design

In this section we recall some key facts from the loop-shaping design technique introduced by McFarlane and Glover [98], and we discuss how it can be merged with our nonsmooth framework to arrive at a new practical PID loop-shaping design method.

Loop-shaping design

Loop-shaping design is an efficient and practical technique which has been applied successfully to a variety of difficult design problems, see [42,53,120] to cite a few. It proceeds as presented in figure 4.4. Firstly, the open-loop plant G is altered by pre-and postcompensators W 1 and W 2 , respectively, to achieve desired open-loop shapes. Roughly speaking, the shaped plant W 2 GW 1 should have large gains at low frequencies for performance and small gains at high frequencies for noise attenuation. Also, the roll-off rate should not be too large in the mid-frequency range. Secondly, an H ∞ synthesis is we have to select a discretization Ω opt for optimization [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF] such that optimal value f (x * ) = ∥T wz (K(x * ))∥ ∞,d is not too far from the true infinite-dimensional value f ∞ (x * ) = ∥T wz (K(x * ))∥ ∞ . This hinges on a suitable grid generation technique. This is crucial in our approach, but we stress that this takes place in a low-dimensional space, whereas system reduction and identification need heavy large scale linear algebra machinery, and yet remain on the level of heuristics.

Given a continuously differentiable function ϕ : [0, ∞] → R + , we want to construct a finite grid Ω opt such that max ω∈[0,∞] |ϕ(ω)-P ϕ (ω)| ≤ ϑ for a fixed tolerance ϑ, where P ϕ is the piecewise linear function interpolating {ϕ(ω) :

ω ∈ Ω opt }. We call M [•, •] a first-order bound of ϕ : R → R if |ϕ ′ (ω)| ≤ M [ω -, ω + ]
for all ω -< ω + and all ω ∈ [ω -, ω + ]. We need the following preparatory: Lemma 5. Suppose we have constructed a grid Ω on [0, ∞] such that for two consecutive nodes ω i , ω i+1 ∈ Ω and some γ * ≥ max{ϕ(ω i ), ϕ(ω i+1 )} the inequality

M [ω i , ω i+1 ](ω i+1 -ω i ) < 2γ * + 2ϑ -ϕ(ω i ) -ϕ(ω i+1 ). ( 16
)
is satisfied. Then ϕ(ω) < γ * + ϑ for every ω ∈ [ω i , ω i+1 ].
Proof: Suppose on the contrary that there exists

ω * ∈ [ω i , ω i+1 ] such that ϕ(ω * ) ≥ γ * +ϑ.
Then the polygon through ϕ(ω i ), ϕ(ω * ), ϕ(ω i+1 ) has length greater than or equal to L, where

L = √ A 2 + (ω * -ω i ) 2 + √ B 2 + (ω i+1 -ω * ) 2 with A = γ * + ϑ -ϕ(ω i ) and B = γ * + ϑ -ϕ(ω i+1 ). Now L ≥ ℓ, where ℓ := min ω∈[ω i ,ω i+1 ] √ A 2 + (ω -ω i ) 2 + √ B 2 + (ω i+1 -ω) 2 = √ (A + B) 2 + (ω i+1 -ω i ) 2 ,
the minimum being attained at ω = ω i B+ω i+1 A A+B . On the other hand the curve

{(ω, ϕ(ω)) : ω ∈ [ω i , ω i+1 ]} has length L = ∫ ω i+1 ω i √ 1 + ϕ ′ (ω) 2 dω ≤ √ 1 + M [ω i , ω i+1 ] 2 (ω i+1 -ω i )
, and L ≥ L, so in combining the two estimates we get L ≥ ℓ, which yields the estimate

√ 1 + M [ω i , ω i+1 ] 2 ≥ √ (A + B) 2 /(ω i+1 -ω i ) 2 + 1. We deduce M [ω i , ω i+1 ] ≥ (A + B)/(ω i+1 -ω i ), and since A + B = 2γ * + 2ϑ -ϕ(ω i ) -ϕ(ω i+1
), this contradicts [START_REF] Jacobson | Linear state-space systems in infinite-dimensional space: the role and characterization of joint stabilizability/detectability[END_REF].

□

We would now like to apply this to the function ϕ(ω) = σ (T wz (K * , jω)), where K * = K(x * ) is the optimal H ∞ -controller computed by algorithm 2. For that we have to prove differentiability of ϕ. We have the following: Lemma 6. [6, Theorem 2.3] The function ϕ has only a finite number of points of nonsmoothness, and in particular, is of class C 2 in the neighborhood of all primary and secondary peaks (all global and local maxima).

Proof: By [17, Thm. 6.1] the one-parameter family of Hermitian matrices ω → T (ω) = T wz (K * , jω) H T wz (K * , jω) has real analytic eigenvalue functions λ ν (ω), hence ϕ 2 (ω) is a finite maximum of real analytic functions, and then also ϕ because ϕ > 0. The rest of the argument is now as in [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF]. □ Remark 10. In consequence, ϕ is twice continuously differentiable in a neighborhood of each peak, and in particular on a set {ω ∈ [0, ∞] : ϕ(ω) > ∥T wz (K * )∥ ∞ -ϑ 0 } for some ϑ 0 > 0. This means Lemma 5 is applicable.

We use this to construct the grid Ω opt used in (6) as follows. Start with ω 0 = 0. Having constructed ω i , compute an extrapolation ω ♯ i > ω i and obtain

M = max{M [ω i , ω] : ω ∈ [ω i , ω ♯ i ]}. Then choose ω i+1 ∈ (ω i , ω ♯ i ] such that M (ω i+1 -ω i ) < 2 max{ϕ(ω i ), ϕ(ω i+1 )} + 2ϑ -ϕ(ω i ) -ϕ(ω i+1 ). ( 17
) If G(s) is available analytically, then M [•, •] is computable.
In the numerical approach we use a finite difference estimation ϕ ′ (ω) ≈ (ϕ(ω + )-ϕ(ω))/(ω + -ω). Since ϕ is continuously differentiable near the peak values, this gives excellent results. In our experience, the method rarely leads to grids with more than a few hundred of nodes, which allows an efficient solution of the optimization program. A typical example is shown in Figure 2.

The following result justifies our method theoretically. [START_REF] Kato | Perturbation theory for linear operators[END_REF] is used in step 4 of algorithm 1 to construct Ω opt , then the gain γ * achieved by the solution

Theorem 2. If 0 < ϑ ≤ ϑ 0 and if a first-order bound M [•, •] for ϕ = σ (T wz (K * , •)) in tandem with rule
K * of (6) is certified to satisfy γ * ≥ ∥T wz (K * )∥ ∞ -ϑ. Proof: Since ϕ(ω) = σ (T wz (K * , jω)) and ∥T wz (K * )∥ ∞ = max ω∈[0,∞] ϕ(ω), we have to show that γ * ≥ ϕ(ω) -ϑ for every i and all ω ∈ [ω i , ω i+1 ],
where ω i are the nodes of the grid Ω opt constructed in step 4 of algorithm 1 based on the rule [START_REF] Kato | Perturbation theory for linear operators[END_REF]. Since γ * is the gain achieved by the solution K * of ( 6) on that grid, it satisfies γ * = max i ϕ(ω i ). Hence γ * ≥ max{ϕ(ω i ), ϕ(ω i+1 )}, and so condition ( 16) of Lemma 5 is satisfied. Since by Lemma 6 we may apply Lemma 5 to ϕ, we obtain the conclusion ϕ(ω

) ≤ max{ϕ(ω i ), ϕ(ω i+1 )}+ϑ ≤ γ * + ϑ on [ω i , ω i+1 ], as claimed. □ Remark 11.
In rule [START_REF] Kato | Perturbation theory for linear operators[END_REF] we apply [START_REF] Jacobson | Linear state-space systems in infinite-dimensional space: the role and characterization of joint stabilizability/detectability[END_REF] with

γ * = max{ϕ(ω i ), ϕ(ω i+1 )} on each interval [ω i , ω i+1 ].
When it comes to just certifying the optimal value f (x * ) = ∥T wz (K(x * ))∥ ∞ = ϕ(ω * ) in step 6 of algorithm 1, then we can construct an even coarser grid by applying [START_REF] Jacobson | Linear state-space systems in infinite-dimensional space: the role and characterization of joint stabilizability/detectability[END_REF] with γ * = ϕ(ω * ) the same for all [ω i , ω i+1 ]. Namely, our grid can be very coarse at frequencies ω where ϕ(ω) ≪ ϕ(ω * ), and still capture sharp peaks, as illustrated in Figure 2. We refer to this as a verification grid Ω ver . According to our experience, the outlined method to construct Ω opt is well-adapted to discretize the controller design problem.

We can further exploit lemma 5 to obtain information on how close the values γ * of ( 6) and γ ∞ of the infinite-dimensional program (5) are. Writing as before f (x) = ∥T wz (K(x))∥ ∞,d for the discrete H ∞ -norm on Ω opt , and f ∞ (x) = ∥T wz (K(x))∥ ∞ for the true H ∞ -norm, we compare the discretized H ∞ -program min x f (x), i.e., [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF], to the underlying infinite-dimensional min x f ∞ (x), i.e. [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF].

Corollary 1. Let x ∞ be a local minimum of (5) with value γ ∞ , and x * a local minimum of [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF] with value γ * . Suppose a first-order bound in tandem with rule [START_REF] Kato | Perturbation theory for linear operators[END_REF] has been used in step 6 of algorithm 1. Then if x * , x ∞ are within neighborhoods of local optimality of each other, we have f Lemma 5, because construction of the grid uses the bound M [•, •] and rule [START_REF] Jacobson | Linear state-space systems in infinite-dimensional space: the role and characterization of joint stabilizability/detectability[END_REF]

(x ∞ ) ≥ f (x * ) ≥ f ∞ (x * ) -ϑ ≥ f ∞ (x ∞ ) -ϑ ≥ f (x ∞ ) -ϑ. Proof: Indeed, f (x ∞ ) ≥ f (x * ) because x * is a minimum of f on a neighborhood U (x * ), and x ∞ ∈ U (x * ) by hypothesis. Next f (x * ) ≥ f ∞ (x * ) -ϑ by
. Next f ∞ (x * ) ≥ f ∞ (x ∞ ), because x ∞ is a minimum of f ∞ on a neighborhood U (x ∞ ), and x * ∈ U (x ∞ ) by hypothesis. The last inequality is satisfied because f ≤ f ∞ . □
This means comparable locally optimal values of the infinite dimensional H ∞ -program (5) and its approximation ( 6) differ by at most ϑ, our apriori chosen tolerance. Since most of the time our algorithm finds the global minimum of ( 6), this is a very useful result in practice, as it determines the value of the infinite dimensional H ∞ -program (5) within a prior tolerance level ϑ.

The argument remains valid if x ∞ , x * are only approximate local minima, say up to the same tolerance ϑ in the values. Then we get the chain f

(x ∞ ) ≥ f (x * ) -ϑ ≥ f ∞ (x * ) -2ϑ ≥ f ∞ (x ∞ ) -3ϑ ≥ f (x ∞ ) -3ϑ
, so here our approximation [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF] gives the correct value up to an error of 3ϑ in the values.

Applications

In this section, we apply our method to several challenging studies in control of infinite dimensional system, and in particular, to boundary and distributed control of systems of parabolic partial differential equations.

Computation of G(s) in boundary control

We illustrate how our procedure is applied to boundary control of parabolic PDEs. Consider a boundary problem of the form

z t (x, t) = ∑ |α|,|β|≤m (-1) |α| D α (a αβ (x)D β z)(x, t) = 0 (x, t) ∈ Q × [0, ∞) D i-1 ν z(x, t) = U i (x, t) x ∈ ∂Q, i = 1, . . . , m, ( 18
)
where U i are abstract controls acting on the boundary ∂Q. Here for convenience Q is bounded open with ∂Q a compact orientable C ∞ -manifold, the coefficients are a α,β ∈ C ∞ (Q), and uniform ellipticity ∑ |α|,|β|≤m a αβ (x)ξ α ξ β ≥ c|ξ| 2 is assumed for x ∈ Q. Then by [START_REF] Salamon | Infinite dimensional linear systems with unbounded control and observation: a functional analytic approach[END_REF] problem [START_REF] Kiwiel | An aggregate subgradient method for nonsmooth convex minimization[END_REF] may be represented in the abstract form [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF], where however the input space U is potentially still infinite dimensional. To comply with our assumption that K(s) should be finite-rank, i.e., that input and output spaces U, Y should be finitedimensional, we select basis functions ϕ ik on the boundary ∂Q and replace the boundary control action U in ( 18) by a finite-dimensional version

D i-1 ν z(x, t) = N ∑ k=1 ϕ ki (x)u ik (t), x ∈ ∂Q, i = 1, . . . , m,
which now has input space U ≃ R N m . A finite-dimensional output space Y ≃ R p could be obtained by taking measurements of the form

y i (t) = ∫ Q ψ i (x)z(x, t) dx, i = 1, . . . , p,
with another set of basis functions ψ i on Ω representing sensors. For one-dimensional Q point evaluations on ∂Q are possible. This case will be used in our numerical experiments.

Referring to [START_REF] Salamon | Infinite dimensional linear systems with unbounded control and observation: a functional analytic approach[END_REF][START_REF] Curtain | An introduction to infinite-dimensional linear system theory[END_REF][START_REF] Sano | H ∞ -control of a parallel-flow heat exchange process[END_REF] for the correct setup of (1), we directly pass to the computation of G(s). Laplace transforming [START_REF] Kiwiel | An aggregate subgradient method for nonsmooth convex minimization[END_REF] with initial condition z(x, 0) = 0 leads for fixed s ∈ C to the elliptic boundary value problem

sz(x, s) = ∑ |α|,|β|≤m (-1) |α| D α (a αβ (x)D β z)(x, s) = 0 x ∈ Q D i-1 ν z(x, s) = N ∑ k=1 ϕ ki (x)u ik (t) x ∈ ∂Q, i = 1, . . . , m. (19)
Then G ikr (s) is obtained by solving [START_REF] Kiwiel | A linearization algorithm for computing control systems subject to singular value inequalities[END_REF] with

u ik = 1 u i ′ ,k ′ = 0 for (i ′ , k ′ ) ̸ = (i, k),
and by computing y r (s) = ∫ Q ψ r (x)z(x, jω) dx. For the one-dimensional case Q = [0, 1], point evaluations y r (s) = z r (0, s), y r (s) = z r (1, s), or linear combinations of those, are possible.

Remark 12. Computation of G ′ (s) can also be obtained by solving an elliptic boundary value problem, which is [START_REF] Kiwiel | A linearization algorithm for computing control systems subject to singular value inequalities[END_REF] differentiated with respect to s. Since K is known explicitly, this is useful when computing the bound M [•, •] of ϕ in Lemma 5.

Remark 13. In those cases where computations are not performed formally, a high spatial resolution is used to solve [START_REF] Kiwiel | A linearization algorithm for computing control systems subject to singular value inequalities[END_REF] accurately. One such solve can then be interpreted as a function evaluation G ikr (s). If carried out numerically, we perform the computation of G(s) for fixed s = jω with the highest spatial discretization available in our setup. Typically, this is at least as accurate spatially as our final simulation of the closed-loop system. Since the number of inputs and outputs is not very large, pre-computing G(s) will not seriously burden the overall performance of algorithm 1. Since pre-computing G(jω) is done off-line, it neither impedes the optimization phase, nor the plant modeling phase.

Reaction-convection-diffusion equation

We consider a non-linear reaction-convection-diffusion equation with Danckwaerts and von Neumann boundary conditions

∂C(z, t) ∂t = D ∂ 2 C(z, t) ∂z 2 -U (t) ∂C(z, t) ∂z -kC(z, t) (z, t) ∈ [0, L] × [0, ∞) D ∂C(0, t) ∂z -U (t)(C(0, t) -C in ) = 0, ∂C(L, t) ∂z = 0. (20)
The process represents a chemical reaction in a cylindrical plug flow reactor with timevarying flow velocity U (t), constant axial dispersion D, and constant reaction rate k. The dynamics of the reaction A k → B are described by the spatially and temporally varying concentration C(z, t) of reactant A, the concentration of product B being a dependent state. Using online measurement y(t) = C(L, t)-C ss (L) of the concentration of ingredient A at the outflow position z = L we steer the plug flow velocity U (t) to maintain the process in steady-state U ss , C ss (z), y ss , while attenuating measurement noise and a disturbance of the flow velocity, and to enable speedy tracking of set-point changes in the steady-state concentration. We refer to [START_REF] Sano | H ∞ -control of a parallel-flow heat exchange process[END_REF] or [START_REF] Curtain | An introduction to infinite-dimensional linear system theory[END_REF]Example 3.3.5] for the correct setup of this problem as a Hilbert space linear system (1). Fixing a steady-state flow velocity U ss , we compute the corresponding steady-state concentration C ss (z) by solving the one-dimensional boundary value problem

DC ′′ ss (z) -U ss C ′ ss (z) -kC ss (z) = 0 DC ′ ss (0) -U ss (C ss (0) -C in ) = 0, C ′ ss (L) = 0. (21)
Linearization about steady-state with U (t) = U ss + u(t) and C(z, t) = C ss (z) + c(z, t) leads now to the linearized boundary and distributed control problem

c t (z, t) = Dc zz (z, t) -U ss c z (z, t) -C ′ ss (z)u(t) -kc(z, t) Dc z (0, t) -U ss c(0, t) + (C in -C ss (0))u(t) = 0, c z (L, t) = 0. (22)
The linearized output is y(t) = c(L, t). In this case steady-state C ss (z) and transfer function G(s) = y(s)/u(s) of the linearized equation can be computed formally using Maple. We obtain

C ss (z) = C in b (b -f )e f (1-z/L)-bz/L 2a -(b + f )e f (z/L-1)-bz/L 2a (-bf -2ak -b 2 )e -f 2a -(bf -2ak -b 2 )e f 2a (23)
where

a = D L 2 , b = -U ss L , f = √ b 2 + 4ak. (24)
Then the transfer u(s) → c(z, s) is obtained analytically as C in P 1 (z,s) s) , where

P 2 (z,
P 1 (z, s) = [ 2La(s + k) ∫ 1 0 f 1 (x)dx -LT 2 ∫ z L 0 f 2 (x)dx -T 3 (m -1) ] e T 4 (z) + [ LT 6 ∫ z L 0 f 1 (x)dx + 2La(s + k) ∫ 1 0 f 2 (x)dx + T 8 (m -1) ] e T 9 (z) + e T 5 (z) LT 2 ∫ z L 1 f 1 (x)dx + e T 7 (z) LT 6 ∫ 1 z L f 2 (x)dx P 2 (z, s) = LT 1 (T 2 e T 1 2a + T 6 e -T 1 2a )
with

T 1 = √ b 2 + 4a(s + k), T 2 = bT 1 + 2a(-k -s) -b 2 , T 3 = bT 1 + 4a(-k -s) -b 2 , T 4 (z) = z L b + ( z L -1)T 1 -2a , T 5 (z) = -z L b + ( z L + 1)T 1 2a , T 6 = bT 1 + 2a(s + k) + b 2 , T 7 (z) = z L b + ( z L + 1)T 1 2a , T 8 = bT 1 + 4a(s + k) + b 2 , T 9 (z) = -z L b + ( z L -1)T 1 2a f 1 (z) = C ssz (z)e (b-T 1 )z 2aL , f 2 (z) = C ssz (z)e (b+T 1 )z 2aL , m = C ss (0) C in .
The transfer u → y is then obtained at z = L. Adopting numerical values D = 1.05 m 2 /min, U ss = 1.24 m/min, C in = 0.5 mol/m 3 , k = 0.25 m 3 /mol, and L = 6.36 m from a study in [START_REF] Fogler | Elements of Chemical Reaction Engineering, Fourth Edition[END_REF], we can compare the the infinitedimensional transfer function G(s) with a finite-difference approximation G fd (s). Figure 3 left shows the comparison of G(s) and G fd (s).

The scheme for synthesis is shown in Figure 4 and uses the filters in Figure 5 (right), which are defined as W e (s) = 0.00001s + 5 s + 0.25 , W n (s) = 0.00125s 2 + 0.00035s + 0.00005 0.000025s 2 + 0.007s + 1 , W u = 0.1.

The controller structure K pi includes SISO PI-controllers with two parameters

K(s) = k p + k i s , so x = (k p , k i ).
With the mixed performance-robustness channel w = (r, n) → z = (z e , z u ) = (W e e, W u u) we have now defined the objective ∥T wz (K(x))∥ ∞ of our problem, and according to section 4.1 the objective f is complemented by the barrier function c∥S(K(x))∥ ∞ . Note that in the scheme of Figure 4 the sensitivity function S equals the unfiltered closed-loop transfer function T re , so altogether with regard to the general form [START_REF] Boyd | A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L ∞ -norm[END_REF]

Van de Vusse reactor

The following more challenging study uses a parallel Van de Vusse process

A k 1 -→ B k 2 -→ C A + A k 3 -→ D
which we want to operate in continuous mode in an isothermal reactor. In [START_REF] Spingarn | Submonotone subdifferentials of Lipschitz functions[END_REF] Van de Vusse considers A = cyclopentadiene, B = cyclopentenol, C = cyclopentanediol, D = dicyclopentadiene. The desired output product is B, the inlet is the primary product A, whereas C, D are waste products. The reaction takes place in a steady-plug-flow cylindrical chemical reactor of length L, through which a fluid with velocity U is flowing, where ingredients are mixing axially with dispersion coefficient D, while the now nonlinear chemical reaction is described by the rates k i . In ODE-based models the axial dispersion is often neglected, or replaced by a singular perturbation approach [START_REF] Dao | Nonconvex bundle method with application to a delamination problem[END_REF]. Here we discuss the full non-linear model. For the functional analytic setup of the problem see

Van de Vusse reactor

The following more challenging study uses a parallel Van de Vusse process

A k 1 -→ B k 2 -→ C A + A k 3 -→ D
which we want to operate in continuous mode in an isothermal reactor. In [START_REF] Van De Vusse | Plug-flow type reactor versus tank reactor[END_REF] Van de Vusse considers A = cyclopentadiene, B = cyclopentenol, C = cyclopentanediol, D = dicyclopentadiene. The desired output product is B, the inlet is the primary product A, whereas C, D are waste products. The reaction takes place in a steady-plug-flow cylindrical chemical reactor of length L, through which a fluid with velocity U is flowing, where ingredients are mixing axially with dispersion coefficient D, while the now nonlinear chemical reaction is described by the rates k i . In ODE-based models the axial dispersion is often neglected, or replaced by a singular perturbation approach [START_REF] Dochain | Approximation of the dynamical model of fixed bed reactors via a singular perturbation approach[END_REF]. Here we discuss the full non-linear model. For the functional analytic setup of the problem see again [START_REF] Sano | H ∞ -control of a parallel-flow heat exchange process[END_REF].

Assuming radially homogeneous conditions in the tube, the system can be described by one spatial dimension z, and the reaction for ingredients A and B is governed by the following diffusion-convection-reaction system of parabolic PDEs:

∂C A (z, t) ∂t = D ∂ 2 C A (z, t) ∂z 2 -U ∂C A (z, t) ∂z -k 1 C A (z, t) -k 3 C 2 A (z, t) ∂C B (z, t) ∂t = D ∂ 2 C B (z, t) ∂z 2 -U ∂C B (z, t) ∂z + k 1 C A (z, t) -k 2 C B (z, t) (25)
for (z, t) ∈ [0, L] × [0, ∞), with Danckwaerts and von Neumann boundary conditions

D ∂C A (0, t) ∂z -U (C A (0, t) -C Ain ) = 0, D ∂C B (0, t) ∂z -U C B (0, t) = 0 ∂C A (L, t) ∂z = 0, ∂C B (L, t) ∂z = 0 (26)
for all t ∈ [0, ∞). The meaning of these boundary conditions at z = 0 is that as soon as the feed enters the reactor at z = 0, it will be diluted by the axial mixing caused by the flow. At z = L we have Neumann boundary conditions, which simply require that the concentration stops changing at the point where the flow leaves the reactor.

The goal of the study is to operate the reactor at a steady-state flow U ss leading to a steady outflow C Bss (L) of product B at the outlet of the reactor. This steady-state flow has to be controlled by feedback, where we have the possibility to act on the velocity U (t) = U ss + u(t), and where we use the deviation y(t) = C B (L, t) -C Bss (L) from the steady-state production as our online measurement at the outlet. It is assumed that changes of the axial flow velocity do not affect the axial dilution D assumed constant. Control has to maintain a stable steady-state, attenuate measurement noise and disturbances at the inflow, and enable the system to react to set-point changes in the flow velocity U .

Our procedure starts by computing the steady-state, which leads to solving the system of ODEs

DC ′′ Ass (z) -U ss C ′ Ass (z) -k 1 C Ass (z) -k 3 C 2 Ass (z) = 0 DC ′′ Bss (z) -U ss C ′ Bss (x) + k 1 C Ass (z) -k 2 C Bss (z) = 0 (27)
with steady-state boundary conditions

DC ′ Ass (0) -U ss (C Ass (0) -C Ain ) = 0, DC ′ Bss (0) -U ss C Bss (0) = 0, C ′ Ass (L) = 0, C ′ Bss (L) = 0. (28)
This defines a mapping U ss → (C Ass (L), C Bss (L)), which allows us to see what flow U ss gives the largest output. In the numerical study we fix U ss = 6.175e-3, which leads to the solutions in Figure 6.

Remark 15. According to our strategy we assume that C Ass (z), C Bss (z) are computed with a very high precision, representing a quasi-analytic solution. Indeed, it is in principle possible to solve the steady-state system formally using Taylor series expansions, but since the difference with a high precision numerical solution is marginal, we proceed with the numerical approach. Once the steady-state is computed, we linearize the system by putting

C A (z, t) = C Ass (z) + c A (z, t), C B (z, t) = C Bss (z) + c B (z, t), U (t) = U ss + u(t) with off-sets c A , c B , u, which leads to the linearized system ∂c A (z, t) ∂t = D ∂ 2 c A (z, t) ∂z 2 -U ss ∂c A (z, t) ∂z -(k 1 + 2k 3 C Ass (z))c A (z, t) - ∂C Ass (z) ∂z u(t) ∂c B (z, t) ∂t = D ∂ 2 c B (z, t) ∂z 2 -U ss ∂c B (z, t) ∂z + k 1 c A (z, t) -k 2 c B (z, t) - ∂C Bss (z) ∂z u(t) (29)
with left boundary conditions

D ∂ ∂z c A (0, t) -U ss c A (0, t) + (C Ain -C Ass (0))u(t) = 0 D ∂ ∂z c B (0, t) -U ss c B (0, t) -C Bss (0)u(t) = 0. (30)
and right boundary conditions

∂c A (L, t) ∂z = 0, ∂c B (L, t) ∂z = 0. (31)
The measured output is y(t) = c B (L, t). The transfer function G(s) = y(s)/u(s) is in principle also available analytically, but we continue with the high precision numerics strategy. We first compute the transfers c A (z, s)/u(s) and c B (z, s)/u(s), which we obtain by Laplace transforming the linearized system. This leads to the linear boundary value problem

D(c A ) zz (z, s) -U ss (c A ) z (z, s) -(k 1 + 2k 3 C Ass (z) -s)c A (z, s) -(C Ass ) z (z)u(s) = 0 D(c B ) zz (z, s) -U ss (c B ) z (z, s) + k 1 c A (z, s) -(k 2 + s)c B (z, s) -(C Bss ) z (z)u(s) = 0
with Laplace transformed boundary conditions

D(c A ) z (0, s) -U ss c A (0, s) + (C Ain -C Ass (0))u(s) = 0, D(c B ) z (0, s) -U ss c B (0, s) + (C Bin -C Bss (0))u(s) = 0 D(c A ) z (L, s) = 0 D(c B ) z (L, s) = 0
Solving this boundary value problem for fixed s = jω with u(s) ≡ 1 gives the value G(jω) of the transfer function y(s) = c B (L, s) = G(s)u(s). This is in fact the method presented in section 5.1. The magnitude of G(s) is shown in Figure 7 (right). Note that for every fixed s = jω we have to solve a static elliptic problem associated with the dynamic equation ( 29), and we perform these computations with the finest scale available, so that G(s) is essentially lossless. In Figure 3 it can be seen that in order to achieve the accuracy in G(s) with a finite-difference discretization G f d (s) we would requires at least 2000 states. So in state-space we would have to perform synthesis on a system of order 2000 to be sure that we do not lose information. This size is beyond existing synthesis techniques. In the same vein, any approach based on system reduction would run the risk of losing information in forming a transfer function based on a reduced model. In the synthesis scheme of Figure 4 we use again the channel (r, n) → (z e , z u ), now with the filters W e (s) = 10 -5 s + 1.502 s + 0.07509 , W n (s) = 0.00125s 2 + 0.00035s + 5 • 10 -5 2.5 • 10 -5 s 2 + 0.007s + 1 , W u = 0.1.

Optimization is now over the class K 3 of third order controllers, which leads to a tunable vector x ∈ R 14 , as the system matrix A K of the controller (2) is parametrized in tridiagonal form. The channel is again complemented by the barrier c∥S(K(x))∥ ∞ , where S = T re The blue curve (z = L) is the one chosen in the experiment. Right: Response of closed-loop system to a noisy reference for optimal controller computed by algorithm 1.

and c = 0.2. The optimal H ∞ -controller K(x * ) with x * ∈ R 14 computed by algorithm 1 is obtained in the form (2) as the optimal gain being f (x * ) = 0.464. The final number of nodes required for a certified result was |Ω opt | = 101, where one update of the grid in step 6 of algorithm 1 was needed. The study ends with a non-linear simulation of the optimal controller. In Figure 7 (right) the response if the nonlinear system to a noisy reference signal is displayed.

A K =   -.

Cavity flow

We consider a challenging cavity flow study from [START_REF] Yan | Experimental study of linear closed-loop control of subsonic cavity flow[END_REF], where the infinite-dimensional transfer function is available analytically and of the form G(s) = e -τ 1 s p 2 (s) + q 2 (s)e -τ 2 s + ce -τ 3 s with quadratic polynomials p 2 (s), q 2 (s) and delay parameters τ i > 0. Figure 8 (left) shows the magnitude plot of G in blue, indicating a large number of resonant peak frequencies.

As H ∞ -objective we have chosen the channel ∥(W 1 S, W 2 T )∥ ∞ with S the closed-loop sensitivity function, T the complementary sensitivity functions, and with the frequency weighing filters W 1 (s) = (0.01s + 177.4)/(s + 50.68), W 2 (s) = (100s + 500)/(s + 50000). Optimization ( 6) is over the class K 2 of order 2 controllers, which features 9 tunable parameters. The optimal controller K(x * ) with x * ∈ R 9 computed by algorithm 1 is given in transfer function form as K * (s) = (0.5069s 2 +119.2s+1419)/(s 2 +308s+1.276e04), and achieves a gain of γ * = 1.937, which improves over the value 1.948 obtained in [START_REF] Yan | Experimental study of linear closed-loop control of subsonic cavity flow[END_REF] using a coprime factorization approach. The final grid size is |Ω opt | = 382. The fact that the class K 2 with which we achieve γ * = 1.937 is much simpler than the infinite-dimensional controller structure used in [START_REF] Yan | Experimental study of linear closed-loop control of subsonic cavity flow[END_REF] explains why our novel approach represents a significant improvement in this study.

non-dimensional pressure referenced to the commonly used value of 20 µPa. Thirty two narrowband power spectra, each from 8192 points, are computed using fast Fourier transform with Hanning window, converted to Sound Pressure Level (SPL) spectra, and then averaged. The resulting spectra have a spectral resolution of about 24 Hz and are accurate within ±1 dB. A stereo particle imaging velocimetry (PIV) is

Figure 1 :

 1 Figure 1: Lower feedback interconnection F ℓ (P, K) with H ∞ -performance channe

Figure 4 . 3 :

 43 Figure 4.3: Selection of frequencies: 1)lhs active only, 2)rhs active and secondary peaks.

Figure 2 :

 2 Figure 2: Left: grid Ω opt based on the criterion (16). Right: selection of extended set of frequencies around active frequencies with primary and secondary peaks (i.e., global and local maxima of ω → σ(T wz (jω)).

10

 10

Figure 3 :

 3 Figure 3: Bode plot of infinite-dimensional transfer function G(s) compared with finite-difference based G f d (s). Left shows study 1 computed with Maple, compared to finite-differences of order N = 1000. Even with 5000 states the transmission zero at frequency ω1.6Hz, which is missed by the discretization of order N = 1000. Right shows Van de Vusse study with G(s) computed via numerical Maple solve of (29), compared to finite-differences of order N = 2000.

Figure 4 :

 4 Figure 4: Scheme for synthesis. The H ∞ -norm of the performance channel (r, n) → (z e , z u ) is minimized, which assures that the system reacts to a set-point change r, and attenuates noise with bandwidth specified by the filter W n . Tracking e is in the low-frequency range specified by W e .

Figure 5 :

 5 Figure 5: Left: Time response to a noisy reference signal of analytic closed loop. Right: Filter W e for synthesis and noise filter W n .

Figure 4 :

 4 Figure 4: Scheme for synthesis. In studies 5.2 and 5.3 the H ∞ -norm of the performance channel (r, n) → (z e , z u ) is minimized, which assures that the system reacts to a set-point change r, and attenuates noise with bandwidth specified by the filter W n . Tracking e is in the low-frequency range specified by W e . In study 5.4 the channel r → (z e , z y ) is optimized.

Figure 5 :

 5 Figure 5: Left: Time response to a noisy reference signal of analytic closed loop. Right: Filter W e for synthesis and noise filter W n .

Figure 6 :

 6 Figure 6: Steady-state flow computed by Maple for U ss = 6.175e-3. The parameters were maxmesh=1e6, abserr = 1e-7. The plots show C Ass (z i ), C Bss (z i ) for 500 equidistant nodes on [0, 1].

Figure 7 :

 7 Figure 7: Left: Comparison of transfer function magnitude for different positions of the sensor.

Fig. 1 .

 1 Fig. 1. Cutout of the wind tunnel showing the converging nozzle, the test section, the cavity, the actuator coupling, and the placement of a Kulite transducer in the cavity floor.

Fig. 2 .Figure 8 :

 28 Fig. 2. Empirical Rossiter acoustic modes, and measu Mach number.

  Increase inner loop counter k and go back to step 4.

	▷ Step 8 (Update memory radius). Store new memory trust-region radius

  Ass (z) steady-state concentration of A mol/l C Bss (z) steady-state concentration of B mol/l

			denomination	numerical value	unit
	k 1	exchange rate A → B	1.39 × 10 -2	s -1
	k 2	exchange rate B → C	2.78 × 10 -2	s -1
	k 3	exchange rate A + A → D	2.77 × 10 -4	l/mol/s
	D	axial dispersion coefficient	3.33 × 10 -4	m 2 /s
	U ss	steady-state velocity	6.175 × 10 -3	m/s
	C Ain	inlet concentration of component A		10.0	mol/l
	L	length of reactor		1.0	m
	The spatiotemporally, spatially, and temporally varying quantities are
		Quantity	denomination		unit
		C A (z, t)	concentration of reactant A	mol/l
		C B (z, t)	concentration of reactant B	mol/l
		C U (t)	superficial velocity		m/s

We have presented a novel method to compute H ∞ -controllers for infinite dimensional systems and in particular for boundary and distributed control of PDEs. At the core our approach uses a non-smooth trust-region bundle algorithm to solve a frequency discretized version of the infinite-dimensional problem. The method was justified theoretically and tested numerically on a reaction-convection-diffusion equation, on a Van de Vusse reactor, and for control of a cavity flow. A convergence certificate for the non-smooth trust-region algorithm under Kiwiel's aggregation rule was proved, allowing to limit the number of cuts in the tangent program (8) to any fixed number N ≥ 3, and answering in the affirmative a question left open in [5].