A Characterization of Underspecified Anaphora and its Consequences on the Annotation of Anaphoric Relations

When resolving anaphora, the precise identification of the antecedent is sometimes difficult. In a first set of cases linked to abstract anaphora ("in that context", "it happens"), the reason lies in the exact delimitation of the antecedent (previous verbal phrase? previous sentence? whole paragraph?). In a second set of cases, several concrete antecedents are possible and the hearer or reader cannot decide between them. Moreover, deciding is not mandatory to the comprehension of the utterance or discourse, and the antecedent can stay underspecified. This phenomenon can be called 'underspecified anaphora', or anaphora with 'fuzzy' or 'vague' antecedent. We illustrate it with French data, and we propose a first classification of types of potential underspecified antecedent. We then deduce principles for their annotation. A characterisation of underspecified anaphora in French We focus on the classical French personal pronouns "il" and "elle" [START_REF] Kleiber | Anaphores et pronoms[END_REF] when they are ambiguous. To simulate the ambiguity in the English version, two different pronouns are sometimes used.

(1) « Je lui ai dit sur un ton de plaisanterie que "Do you want to take John Dupont to be your husband? -Yes, I want him/it." (6) « Jean a coupé tout le bois. Il l'a fait en petits morceaux ».

"John cut wood. He did/made it in small parts." These examples put forward nominal phrases like "le N 1 de le N 2 " as well as possessives, coordinations and juxtapositions. Attributing a precise antecedent is each time not really crucial for the comprehension of the text. In (1), the person as well as her idea can show things in an unexpected way. Even if the two semantic representations differ, they convey the same global signification, since the person can be fused to her idea in the context of this sentence. In (2) and (3), a deep analysis can be made to identify the referent the writer or speaker had in mind, but the result is not guaranteed and the underspecified solution may be better. In (4), two possible interpretations are the precision and the correction. The juxtaposition gives no indication to choose, and both stay relevant. In some other unambiguous examples, the antecedent is the first noun: "le rêve, l'évasion est nécessaire. Il permet de mieux supporter les tracas de la vie quotidienne" ("dreaming, daydreaming is necessary. It allows feeling better"). In ( 5) and ( 6), the ambiguity implies a concrete antecedent and an abstract one. Following the point of view of [START_REF] Prandi | Sémantique du contresens[END_REF] concerning some syntactic ambiguities, we can say that underspecified anaphora is more frequent than we think (in French in any case), since there is often no simple linguistic way to avoid the ambiguity. With a natural language processing point of view, we can notice that classical algorithms for anaphora resolution [START_REF] Mitkov | Anaphora resolution: to what extent does it help NLP applications?[END_REF], and for instance salience-based methods [START_REF] Miltsakaki | A rethink of the relationship between salience and anaphora resolution[END_REF] are useless face to this phenomenon.

Confronting these examples to others from scientific studies or Frantext corpus, we propose the following characterization of potential underspecified antecedents:

1. Possessives, when the possessor is animated and the possession is a personality trait to which it can be assimilated, or when the possessor is an object and the possession is its main function.

The key is that both of them have to be linked (one may be a facet of the other, or may belong to its set of properties). 2. Complex nominal phrases "le N 1 de le N 2 ", with the same participants than in the previous situation.

3. Coordinations, in the same case and in cases where the sentence elements are in plural: "the N 1 and the N 2 […] they […]". 4. Juxtapositions, in particular when it is hard to distinguish a simple enumeration from a reformulation, and, in the latter case, to distinguish a precision from a correction (example 4). Another example implies several persons and involves a possible precision: in "the First Secretary, Mister Smith, and his wife", how many persons are mentioned? two or three? 5. Evolutive referents: in examples such as "as a child, Marcello […], but adult Marcello is […] he […]", does "he" refer to "Marcello" or to "adult Marcello" (or to "child Marcello" if we consider not only text spans but also semantic referents)? It is sometimes impossible and useless to make a choice. To these cases we have to add the situations where the antecedent is not linguistically mentioned. This is the case with collective pronouns: "ils ont encore augmenté les impôts", "they have increased taxes again" [START_REF] Kleiber | Anaphores et pronoms[END_REF]. This is also the case with antecedentless anaphora [START_REF] Cornish | Antecedentless' anaphors: deixis, anaphora, or what? Some evidence from English and French[END_REF]. Consequences on the annotation of anaphoric relations A way to apprehend the identification of antecedents consists of annotating (as a means to represent semantic concerns and not to evaluate computational systems). Here is a first list of annotation principles linked to underspecified anaphora:

1. The determination of an underspecified antecedent is necessary if we don't want the interpretation to be reduced to a particular or irrelevant choice. 2. Alternatives principle: in most of cases, several alternatives confront themselves and the choice does not matter (the person or her idea, the whole documents or only the reports, etc.).

The aim is to determine what are the alternatives, and to group them as the antecedent. The initial intended antecedent from the writer or speaker is in the set of alternatives, but all of them are relevant for the reader or hearer (and for the semantic representation determination). The underspecification cannot be the source of a mistake or a misunderstanding. This principle has to be distinguished from the principle consisting of grouping some antecedents into a composite one, as in "John was sleeping. Mary was reading. They were happy". 3. Feasibles principle: in some other cases, alternatives may be hard to determine, for instance with juxtapositions. The aim here is to identify all the possibilities and to label them as feasible.

The intended antecedent from the speaker is in the set of feasibles, but only some of them are relevant for the semantic representation determination, i.e., underspecification may be a cause of mistake. This is the case with the example "the First Secretary, Mister Smith, and his wife". 4. Double mark-up principle: for abstract anaphora such as "it happens", alternatives and feasibles are not sufficient. The best solution would be to manage fuzzy or progressive limits for the antecedent, which is not compatible with mark-up constraints. We propose to manage two sets of tags, a first one corresponding to the limit inferior and the second one to the limit superior. To these tags correspond two discourse entities, and every entity in-between is acceptable. This principle differs from the 'MIN' attribute from MUC-7 (Message Understanding Conferences). The 'MIN' attribute has been designed to reduce a referential expression to a minimal mention like a proper noun, and emphasizes in this way a characteristic of a particular discourse entity. In the case of underspecified anaphora, there are several discourse entities in presence. Each of them can be characterised by a minimal mention and therefore with a 'MIN' attribute. The double mark-up principle is then independent from the 'MIN' attribute.


Member, IEEE Pierre Apkarian

Member, SIAM Dominikus Noll

Nonsmooth H ∞ Synthesis

Keywords: H∞-synthesis, static output feedback, fixedorder synthesis, simultaneous stabilization, NP-hard problems, nonsmooth optimization, bundle methods, Clarke subdifferential, bilinear matrix inequality (BMI), linear matrix inequality (LMI)

We develop nonsmooth optimization techniques to solve H∞ synthesis problems under additional structural constraints on the controller. Our approach avoids the use of Lyapunov variables and therefore leads to moderate size optimization programs even for very large systems. The proposed framework is versatile and can accommodate a number of challenging design problems including static, fixed-order, fixed-structure, decentralized control, design of PID controllers and simultaneous design and stabilization problems.

Our algorithmic strategy uses generalized gradients and bundling techniques suited for the H∞ norm and other nonsmooth performance criteria. We compute descent directions by solving quadratic programs and generate steps via line search. Convergence to a critical point from an arbitrary starting point is proved and numerical tests are included to validate our methods. The propose approach proves to be efficient even for systems with several hundreds of states.

I. INTRODUCTION

In this paper we consider H ∞ synthesis problems with additional structural constraints on the controller. This includes static and reduced-order H ∞ output feedback control, structured, sparse or decentralized synthesis, simultaneous stabilization problems, multiple performance channels, and much else. We propose to solve these problems with a nonsmooth optimization method exploiting the structure of the H ∞ norm.

In nominal H ∞ synthesis, feedback controllers are computed via semidefinite programming (SDP) [START_REF] Gahinet | A linear matrix inequality approach to H∞ control[END_REF] or algebraic Riccati equations [START_REF] Doyle | State-Space Solutions to Standard H 2 and H∞ Control Problems[END_REF]. When structural constraints on the controller are added, the H ∞ synthesis problem is no longer convex. Some of the problems above have even been recognized as NP-hard [START_REF] Nemirovskii | Several NP-Hard problems arising in robust stability analysis[END_REF] or as rationally undecidable [START_REF] Blondel | Simultaneous stabilizability question of three linear systems is rationally undecidable[END_REF]. These mathematical concepts indicate the inherent difficulty of H ∞ synthesis under constraints on the controller.

Even with structural constraints, the bounded real lemma may still be brought into play. The difference with customary H ∞ synthesis is that it no longer produces LMIs, but bilinear matrix inequalities, BMIs, which are genuinely non-convex. Optimization code for BMI problems is currently developed by several groups, see e.g. [START_REF] Henrion | Solving simultaneous stabilization BMI problems with PENNON[END_REF]- [START_REF] Leibfritz | Trust region methods for solving the optimal output feedback design problem[END_REF], but it appears that the BMI approach runs into numerical difficulties even for problems of moderate size. This is mainly due to the presence of Lyapunov variables, whose number grows quadratically with the number of states. Out present approach does not use the bounded real lemma and thereby avoids Lyapunov variables. This leads to moderate size optimization programs even for very large systems. In exchange, our cost functions are nonsmooth and require special optimization techniques, which we develop here. We evaluate the H ∞ norm via the Hamiltonian bisection algorithm [START_REF] Boyd | A bisection method for computing the H∞ norm of a transfer matrix and related problems[END_REF]- [START_REF] Gahinet | Numerical computation of the L∞ norm revisited[END_REF] and exploit it further to compute subgradients, which are then used to compute descent steps. Notice, however, that our method is not a pure frequency domain method. In fact, it allows both frequency domain and state space domain parameterizations of the unknown controller. This makes it a very flexible tool in a number of situations of practical interest.

Several iterative methods for reduced-order control have been proposed over recent years, see for instance [START_REF] De Oliveira | Numerical comparison of output feedback design methods[END_REF]- [START_REF] Grigoriadis | Low-order control design for LMI problems using alternating projection methods[END_REF]. In [START_REF] De Oliveira | Numerical comparison of output feedback design methods[END_REF], a comparison among four of these methods on a large set of test problems is arranged, with the result that successive linearization [START_REF] Elghaoui | An algorithm for static outputfeedback and related problems[END_REF], also known as the Frank and Wolfe (FW) algorithm [START_REF] Frank | An algorithm for quadratic programming[END_REF], performed best. Whenever possible, we have therefore compared our new nonsmooth methods and the augmented Lagrangian algorithm in [START_REF] Apkarian | Fixed-order H∞ control design via an augmented Lagrangian method[END_REF], [START_REF] Noll | Partially augmented Lagrangian method for matrix inequality constraints[END_REF] with the Frank and Wolfe method. The results are presented in the experimental section.

As far as comparison with existing methods is concerned, let us mention that for specific classes of plants, it is possible to compute reduced-order controllers without the use of optimization techniques. This has for instance been investigated in [START_REF] Xin | Reduced-order controllers for the H∞ control problem with unstable invariant zeros[END_REF]- [START_REF] Stoorvogel | A reduced-order oberver based controller design for H∞-optimization[END_REF]. These approaches usually make strong additional assumptions like singularity, or hypotheses about unstable invariant zeros. In such cases it may then even be possible to assure global optimality of the computed controllers. Unfortunately, in these approaches, the order of the controller is not a priori known, and in particular, it is not possible to compute static controllers with this type of technique. In the absence of these additional assumptions, and in particular when structural constraints are imposed, synthesis via nonlinear optimization appears to be the most general approach to H ∞ synthesis.

The structure of the paper is as follows. In section II we present the H ∞ synthesis problem and give several motivating examples. Section III computes subgradients of the H ∞ norm, which are then applied to closed-loop scenarios in section IV. In section V we start to develop our first-order descent method, which is completed in section VI. Section VI-G discusses practical aspects of the method, and the final section VII presents a number of experiments to validate our approach.

NOTATION

Let M n,m be the space of n×m matrices, equipped with the corresponding scalar product X, Y = Tr(X H Y ), where X H is the transconjugate of the matrix X, Tr X its trace. The space of m × m Hermitian matrices is denoted S m . For Hermitian or symmetric matrices, X ≻ Y means that X -Y is positive definite, X Y that X -Y is positive semi-definite. For ease of notations, we define the following sets of Hermitian matrices: B m := {X ∈ S m : X 0, Tr (X) = 1}. Consider q-tuples of Hermitian matrices (Y 1 , . . . , Y q ), we define the set B q m := {(Y 1 , . . . , Y q ) : Y i ∈ S m , Y i 0, q i=1 Tr (Y i ) = 1}. For short, we shall use B or B q when the dimension needs not be specified.

We use the symbol λ 1 to denote the maximum eigenvalue of a symmetric or Hermitian matrix. We shall use general notions from nonsmooth analysis covered by [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]. Notions on ǫ-subdifferentials and ǫ-enlarged subdifferentials for spectral functions and their relationships are discussed at length in [START_REF] Cullum | The minimization of certain nondifferentiable sums of eigenvalues of symmetric matrices[END_REF]- [START_REF] Oustry | A second-order bundle method to minimize the maximum eigenvalue function[END_REF]. Unless stated otherwise, the symbol x designates a vector gathering (controller) decision variables and must not be confused with the plant state in Section II. In the notation x k , the subscript k refers to the iteration index.

II. H ∞ SYNTHESIS

The general setting of the H ∞ synthesis problem is as follows. We consider a linear time-invariant plant described in standard form by the state-space equations:

P (s) :   ẋ z y   =   A B 1 B 2 C 1 D 11 D 12 C 2 D 21 D 22     x w u   , (1)
where x ∈ R n is the state vector, u ∈ R m2 the vector of control inputs, w ∈ R m1 the vector of exogenous inputs, y ∈ R p2 the vector of measurements and z ∈ R p1 the controlled or performance vector. Without loss of generality, it is assumed throughout that D 22 = 0.

Let u = K(s)y be a dynamic output feedback control law for the open-loop plant [START_REF] Gahinet | A linear matrix inequality approach to H∞ control[END_REF], and let T w→z (K) denote the closed-loop transfer function of the performance channel mapping w into z. Our aim is to compute K(s) such that the following design requirements are met:

• Internal stability: For w = 0 the state vector of the closed-loop system (1) and (2) tends to zero as time goes to infinity.

• Performance: The H ∞ norm T w→z (K) ∞ is minimized
among all stabilizing K. We assume that the controller K has the following frequency domain representation:

K(s) = C K (sI -A K ) -1 B K + D K , A K ∈ R k×k , ( 2
)
where k is the order of the controller, and where the case k = 0 of a static controller K(s) = D K is included. Often practical considerations dictate additional challenging structural constraints. For instance it may be desired to design low-order controllers (0 ≤ k ≪ n) or controllers with prescribed-pattern, sparse controllers, decentralized controllers, observed-based controllers, PID control structures, synthesis on a finite set of transfer functions, and much else. Formally, the synthesis problem may then be represented as minimize

T w→z (K) ∞ subject to K stabilizes (1) K ∈ K (3)
where K ∈ K represents a structural constraint on the controller (2) like one of the above. Without the restriction K ∈ K, and under standard stabilizability and detectability conditions, it has become customary to synthesize K(s) as follows. After substituting (2) into (1), the H ∞ synthesis problem is transformed into a matrix inequality condition using the bounded real lemma [START_REF] Anderson | Network analysis and synthesis: a modern systems theory approach[END_REF]. Then the projection lemma from [START_REF] Gahinet | A linear matrix inequality approach to H∞ control[END_REF] is used to eliminate the unknown controller data A K , B K , C K , D K from the cast, leaving an LMI problem, which may be solved by SDP. In a third step the controller state-space representation (2) is recovered.

This scenario changes dramatically as soon as constraints K ∈ K are added. Then the problem may no longer be transformed into an LMI or any other convex program, and alternative algorithmic strategies are required. The aim of this paper is to present and analyze one such alternative.

Example 1. Pure stabilization. Often the first important step in controller synthesis (3) is to find a stabilizing controller K. Already at this stage the H ∞ norm plays a prominent role, because of the well-known fact that under stabilizability and detectability, a linear-time invariant system is Lyapunov stable if and only if its H ∞ norm is finite [START_REF] Desoer | Feedback Systems: Input-Output Properties[END_REF]. More specifically, under stabilizability and detectability assumptions, the static control law u = Ky stabilizes the plant

G(s) : ẋ = Ax + B 2 u y = C 2 x , (4)
if and only if the closed-loop transfer matrix C 2 (sI -

(A + B 2 KC 2 )) -1 B 2 has finite H ∞ norm.
In order to construct a static stabilizing controller for an unstable open-loop system (4), the following procedure appears fairly natural. Suppose we are given an initial guess K 0 , which leaves the closed-loop system unstable. Then we pick a 0 > 0 such that the a 0 -shifted H ∞ -norm of the closedloop system is finite:

C 2 (sI -(A + B 2 K 0 C 2 )) -1 )B 2 -a0,∞ < +∞,
where the shifted H ∞ norm is given in [START_REF] Boyd | Linear Controller Design: Limits of Performance[END_REF]. The problem of finding a stabilizing controller K may now be addressed by an optimization program minimize

K C 2 (sI -(A + B 2 KC 2 )) -1 B 2 -a,∞ , (5)
where the shift a is either kept fixed at the initial a 0 , or is gradually decreased after each minimization step to accelerate the procedure. A stabilizing controller K is obviously obtained when the shift reaches a ≤ 0, but very often this happens already with the initial value a 0 , so that shifting is not even necessary as a rule. Numerical tests for this method will be presented in section VII. While we stop the optimization (5) as soon as a stabilizing K is reached, it may happen that for a fixed shift a, the method converges to a local minimum K of [START_REF] Henrion | Solving simultaneous stabilization BMI problems with PENNON[END_REF], which fails to stabilize the closed-loop system. This is explained by the fact that [START_REF] Henrion | Solving simultaneous stabilization BMI problems with PENNON[END_REF], just like all the other methods in this paper, are local optimization methods in the sense that they guarantee convergence to a local minimum (or a critical point). If an unsatisfactory local minimum is reached, the only possibility is to do a restart with a new initial guess K 0 , or switch to another method. Such a local convergence certificate may appear weak at first sight, but experience shows that local methods perform much better than global optimization techniques. Those may have stronger certificates, but run into numerical problems even for small problems. And indeed, our present approach is almost always successful even without restart. A similar comment applies to the pure stabilization method in [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF].

Notice that pure stabilization is just a special case of the more general H ∞ synthesis problem in section II when we specialize the standard form to

P (s) :   ẋ z y   =   A B 2 B 2 C 2 0 0 C 2 0 0     x w u   . ( 6
)
Example 2. Spectral abscissa. In [START_REF] Burke | Two numerical methods for optimizing matrix stability[END_REF], [START_REF]Robust stability and a criss-cross algorithm for pseudospectra[END_REF]

K α(A + B 2 KC 2 ), (7)
where α is the spectral abscissa of a matrix M ∈ M n,n defined as

α(M ) = max{Re λ : λ eigenvalue of M }.
Unfortunately, this function is not even locally Lipschitz, which renders application of existing nonsmooth algorithms impossible. The authors have therefore developed a probabilistic algorithm which allows to treat problems like [START_REF] Thevenet | Nonlinear spectral SDP method for BMI-constrained problems: Applications to control design[END_REF]; see [START_REF]Robust stability and a criss-cross algorithm for pseudospectra[END_REF], [START_REF]A robust gradient sampling algorithm for nonsmooth, nonconvex optimization[END_REF]. Since optimality tests and approximate subgradients for α are difficult to compute (see [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF]), we prefer the use of ( 5) over [START_REF] Thevenet | Nonlinear spectral SDP method for BMI-constrained problems: Applications to control design[END_REF]. Numerical tests for [START_REF] Henrion | Solving simultaneous stabilization BMI problems with PENNON[END_REF] are presented in section VII and show that this method is successful as a rule. Our own numerical tests with program [START_REF] Thevenet | Nonlinear spectral SDP method for BMI-constrained problems: Applications to control design[END_REF] are published in [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF].

The above scenario covers the case of a static stabilizing controller K, but it is clear that stabilization problems including structural constraints K ∈ K can be treated in exactly the same way. Several examples of classes K will be presented in the sequel.

γ i C i (sI -(A i + B i KC i ) -a i I) -1 B i ∞
where γ i > 0 are appropriate weights, and the a i are chosen so that the initial guess K renders the ith system stable after shifting by a i . A lower bound for a i is therefore the spectral abscissa (7) of the ith system. As before, the shifts are decreased after each H ∞ norm minimization step and a simultaneously stabilizing K is obtained e.g. if a i ≤ 0 for all i = 1, . . . , r. But even when some a i > 0, the solution may produce a simultaneously stabilizing K.

G stab (s) -G stab (s) , (8)
where Gstab (s) ranges over a prespecified class K of stable system of reduced order, n ≪ N , and where some norm criterion is used to evaluate the mismatch between nominal and reduced systems. If • represents the Hankel norm, an explicit expression for G stab is available [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L∞-error bounds[END_REF]. But it may be preferable to use other criteria like the H ∞ norm, a problem which then falls within the class of problems considered in this work. Once a solution G stab (s) to ( 8) is obtained, the new system G(s) = G instab (s) + G stab (s), while easier to control, may be expected to have characteristics similar to those of the original system.

III. SUBDIFFERENTIAL OF THE H ∞ NORM

In this section, we start characterizing the subdifferential of the H ∞ norm, and derive expressions for the Clarke subdifferential of several nonconvex composite functions f (x) = G(x) ∞ , where G is a smooth operator defined on some R n with values in the space of stable matrix transfer functions H ∞ .

Consider the H ∞ -norm of a nonzero transfer matrix function G(s):

G ∞ = sup ω∈R σ (G(jω)) ,
where G is stable and σ(X) is the maximum singular value of X. Suppose G ∞ = σ (G(jω)) is attained at some frequency ω, where the case ω = ∞ is allowed. Let G(jω) = U ΣV H be a singular value decomposition. Pick u the first column of U , v the first column of V , that is, u = G(jω)v/ G ∞ . Then the linear functional φ = φ u,v,ω defined as

φ(H) = Re u H H(jω)v = G -1 ∞ Re Tr vv H G(jω) H H(jω) = G -1 ∞ Re Tr G(jω) H uu H H(jω)
is continuous on the space H ∞ of stable transfer functions and is a subgradient of • ∞ at G [START_REF] Boyd | Linear Controller Design: Limits of Performance[END_REF]. More generally, assume that the columns of Q u form an orthonormal basis of the eigenspace of G(jω)G(jω) H associated with the largest eigenvalue λ 1 G(jω)G(jω) H = σ(G(jω)) 2 , and that the columns of Q v form an orthonormal basis of the eigenspace of G(jω) H G(jω) associated with the same eigenvalue. Then for all complex Hermitian matrices Y v , Y u ∈ B,

φ(H) = G -1 ∞ Re Tr Q v Y v Q H v G(jω) H H(jω) (9) = G -1 ∞ Re Tr G(jω) H Q u Y u Q H u H(jω) is a subgradient of • ∞ at G.
Finally, with G(s) rational and assuming that there exist finitely many frequencies ω 1 , . . . , ω p where the supremum G ∞ = σ(G(jω ν )) is attained, all subgradients of • ∞ at G are precisely of the form

φ(H) = G -1 ∞ Re p ν=1 Tr G(jω ν ) H Q ν Y ν Q H ν H(jω ν ),
where the columns of Q ν form an orthonormal basis of the eigenspace of G(jω ν )G(jω ν ) H associated with the leading eigenvalue G 2 ∞ , and where (Y 1 , . . . , Y p ) ∈ B p . See [22, Prop. 2.3.12 and Thm. 2.8.2] and [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF] for this.

Suppose now we have a smooth operator G, mapping R n onto the space H ∞ of stable transfer functions G. Then the composite function

f (x) = G(x) ∞ is Clarke subdifferen- tiable at x with ∂f (x) = G ′ (x) ⋆ [∂ • ∞ (G(x))], (10)
where ∂ • ∞ is the subdifferential of the H ∞ -norm obtained above, and where G ′ (x) ⋆ is the adjoint of G ′ (x), mapping the dual of H ∞ into R n , where R n is identified with its dual here.

In the sequel, we will compute this adjoint G ′ (x) ⋆ for special classes of closed-loop transfer functions. Suitable chain rules covering this case are for instance given in [22,

(K) of f := • ∞ •T w→z at K. We first notice that the derivative T ′ w→z (K) of T w→z at K is T ′ w→z (K)δK := G 12 δKG 21 ,
where δK is an element of the same matrix space as K and with the definitions

T w→z (K, s) G 12 (K, s) G 21 (K, s) ⋆ := P 11 + P 12 K(I -P 22 K) -1 P 21 P 12 (I -KP 22 ) -1 (I -P 22 K) -1 P 21 ⋆ = C(K) C 2 (sI -A(K)) -1 [ B(K) B 2 ] + D(K) D 12 D 21 ⋆
and the closed-loop state-space data

A(K) := A + B 2 KC 2 , B(K) := B 1 + B 2 KD 21 , C(K) := C 1 + D 12 KC 2 , D(K) := D 11 + D 12 KD 21 .
Note that numerical stability requires that transfer functions T w→z , G 12 and G 21 be computed through state-space realizations. [START_REF] Leibfritz | Trust region methods for solving the optimal output feedback design problem[END_REF], specified by Y ∈ B and with T w→z (K) ∞ attained at frequency ω. According to the chain rule, the subgradients

Now let φ = φ Y be a subgradient of • ∞ at T w→z (K) of the form
Φ Y of f at K are of the form Φ Y := T ′ w→z (K) ⋆ φ Y ∈ M m2,p2 , where the adjoint T ′ w→z (K) ⋆ acts on φ Y through T ′ w→z (K) ⋆ φ Y , δK = T ′ w→z (K)δK, φ Y = T w→z (K) -1 ∞ Re Tr T w→z (K, jω) H QY Q H T ′ w→z (K)δK(jω) ) = T w→z (K) -1 ∞ Re Tr T w→z (K, jω) H QY Q H G 12 (K, jω)δK(jω)G 21 (K, jω) ) = T w→z (K) -1 ∞ Re Tr G 21 (K, jω)T w→z (K, jω) H QY Q H G 12 (K, jω)δK(jω) ) . (11)
In consequence, for a static K, the Clarke subdifferential of

f (K) := T w→z (K) ∞ at K consists of all subgradients Φ Y of the form Φ Y = T w→z (K) -1 ∞ Re G 21 (K, jω)T w→z (K, jω) H QY Q H G 12 (K, jω) ) T , (12)
where Y ∈ B. Recall that Φ Y is now an element of the same matrix space as K and acts on test vectors δK through Φ Y , δK = Tr(Φ T Y δK). This formula is easily adapted if the H ∞ norm is attained at a finite number of frequencies ω 1 , . . . , ω q . In this more general situation, subgradients of f at K are of the form

Φ Y = T w→z (K) -1 ∞ q ν=1 Re (G 21 (K, jω ν ) T w→z (K, jω ν ) H QY ν Q H G 12 (K, jω ν ) ) T , ( 13
) where Y = (Y 1 , . . . , Y q ) ∈ B q .
At this stage, it is important to stress that expressions (11), ( 12) and ( 13) are general and can accommodate any problem discussed in previous sections. Below we resume and expand this list by considering more examples of practical interest.

Example 5. Dynamic controllers. Assume now that the controller is dynamic as in [START_REF] Doyle | State-Space Solutions to Standard H 2 and H∞ Control Problems[END_REF]. The subgradient set is again obtained via formula (13) by performing the substitutions:

K → A K B K C K D K , A → A 0 0 0 k , B 1 → B 1 0 , C 1 → [ C 1 0 ] , B 2 → 0 B 2 I k 0 , C 2 → 0 I k C 2 0 , D 12 → [ 0 D 12 ] , D 21 → 0 D 21 . (14)
The entire Clarke subdifferential is then described by the set of subgradients in M k+m2,k+p2

Φ Y := Φ Y,AK Φ Y,BK Φ Y,CK Φ Y,DK ,
where Y = (Y 1 , . . . , Y q ) ∈ B q and Φ Y is derived through formulas ( 12) or [START_REF] De Oliveira | Numerical comparison of output feedback design methods[END_REF].

Example 6. Structured controllers. In practice, it is sometimes required that some entries in the controller gain be put to zero, while the others may be freely assigned. This is the case in decentralized control, where the controller must enjoy a block-diagonal structure. Consider a pattern matrix W with entries W ij ∈ {0, 1}, where W ij = 0 means that the controller gain K ij = 0 must be zero, whereas W ℓk = 1 means that K ℓk can be freely assigned. The Clarke subdifferential of [START_REF] De Oliveira | Numerical comparison of output feedback design methods[END_REF] and where ⊙ denotes the entry-wise Hadamard or Schur product [START_REF] Horn | Matrix Analysis[END_REF].

f = • ∞ • F at K is then of the form W ⊙ Φ Y , where Φ Y ∈ ∂ • ∞ (F(K)) is as in
Example 7. PID controllers. PID control is one of the most classical approaches in control system design. The controller is generally written as

K(s) = K P + 1 s K I + s 1 + τ s K D ,
where K P , K I and K D are matrix static gains to be computed, and τ is a small positive scalar. Using the general formula ( 11), the subgradients with respect to κ := [

K P ; K I ; K D ] at K(s) are obtained through Φ Y,κ = T w→z (K) -1 ∞ q ν=1 Re (G 21 (K, jω ν ) T w→z (K, jω ν ) H QY ν Q H G 12 (K, jω ν )[I 1 jων I jων 1+τ jων I] T ,
where as before Y = (Y 1 , . . . , Y q ) ∈ B q . The above approach could be generalized by making τ an additional design parameter. Then a constraint τ ≥ 0 should be added. Notice also that the above formula readily extends to arbitrary basis functions {Q j (s)} j=1,...,r

K(s) := r j=1 K j Q j (s) ,
where the K j 's are the design variables.

Example 8. Matrix fraction representations. An alternative representation of controllers is via matrix fraction descriptions. For instance, the left matrix fraction representation is given as

K(s) = N (s)D(s) -1 = (N J N (s))(DJ D (s)) -1 , with J N (s) := [I sI . . . s n I] T , J D (s) := [I sI . . . s d I] T and N := [ N 0 . . . N n ] , D := [ D 0 . . . D n ]
Now, N and D are the design variables. As before, it is immediate to show that partial subgradients with respect to N and D are given as

Φ Y,N = T w→z (K) -1 ∞ q ν=1 Re (J N (jω ν ) (D(jω ν )J D (jω ν )) -1 G 21 (K, jω ν )T w→z (K, jω ν ) H QY ν Q H G 12 (K, jω ν ) ) T Φ Y,D = -T w→z (K) -1 ∞ q ν=1 Re (J D (jω ν ) (D(jω ν )J D (jω ν )) -1 G 21 (K, jω ν )T w→z (K, jω ν ) H QY ν Q H G12(K, jω ν )K(jω ν ) ) T ,
respectively.

Example 9. Multiple performance channels. Practical specifications often impose that several closed-loop channels i = 1, . . . , r be minimized simultaneously. One way to address multi-objective optimization of this type is to solve a program of the form minimize

K max γ i T i w→z (K) ∞ : i = 1, . . . , r ,
where T i w→z is the ith performance specification to be optimized. Since the maximum of a finite number of maximum eigenvalue functions is itself a maximum eigenvalue function of a block diagonal operator T = diag(T 1 w→z , . . . , T r w→z ), the Clarke subgradients could be obtained directly from [START_REF] De Oliveira | Numerical comparison of output feedback design methods[END_REF]. When the usual max formula is used, the result is the same, i.e., subgradients are of the form

φ (Y,τ ) = i∈I(K) τ i γ i φ Yi ,
where I(K) are the i = 1, . . . , r which are active at K, τ i ≥ 0,

i∈I(K) τ i = 1 and φ Yi ∈ (T i w→z ) ⋆ ∂ • ∞ (K)
as specified in [START_REF] De Oliveira | Numerical comparison of output feedback design methods[END_REF].

Before going further, it is worth mentioning that our methodology carries over to a wide range of controller structures of practical interest. This is in particular the case when the structural constraint is of the form K = {K : K = S(ℓ), ℓ ∈ L}, where (S, L) is a suitable differentiable parametrization of the class K. This includes for instance observed-based controllers, feed-forward compensators, controllers defined through Youla parameterizations and much else.

V. STEEPEST DESCENT METHOD

Nonsmooth techniques have been used before in algorithms for controller synthesis. For instance, E. Polak and co-workers have proposed a variety of techniques suited for eigenvalue or singular-value optimization and for extensions to the semiinfinite case, covering in particular the H ∞ -norm (see [START_REF] Polak | On the mathematical foundations of nondifferentiable optimization in engineering design[END_REF], [START_REF]Optimization : Algorithms and Consistent Approximations[END_REF] and the citations given there). Another reference is [START_REF] Boyd | Linear Controller Design: Limits of Performance[END_REF], where the authors exploit the Youla parameterization via convex nondifferentiable analysis to derive the cutting plane and ellipsoid algorithms.

Let us consider the problem of minimizing f (x) = G(x) ∞ , where x regroups the controller data, referred to as K in the previous section, and where G maps R n smoothly into a space H ∞ of stable transfer functions. We write G(x, s) or G(x, jω) when the complex argument of G(x) ∈ H ∞ needs to be specified.

A necessary condition for optimality is 0

∈ ∂f (x) = G ′ (x) ⋆ ∂ • ∞ (G(x)). It is therefore reasonable to consider the program d = - g g , g = argmin{ φ Y : Y = (Y 1 , . . . , Y q ) ∈ B q } (15)
which either shows 0 ∈ ∂f (x), or produces the direction d of steepest descent at x if 0 ∈ ∂f (x), and where the φ Y are as in [START_REF] De Oliveira | Numerical comparison of output feedback design methods[END_REF]. If we vectorize y = vec(Y ), Y = (Y 1 , . . . , Y q ), then we may represent φ Y by a matrix vector product, φ Y = Φy, with a suitable matrix Φ. Program ( 15) is then equivalent to the following SDP: minimize t subject to t y T Φ T Φy tI 0

Y i 0, i = 1, . . . , q e T y = 1 (16)
where e T y = 1 encodes the constraint i Tr(Y i ) = 1. The direction d of steepest descent at x is then obtained as d = -Φ y/ Φ y , where (t, y) is solution of ( 16) with y = 0. This suggests the following:

Steepest descent method for the H ∞ -norm 1) If 0 ∈ ∂f (x) stop. Otherwise:
2) Solve ( 16) and compute the direction d of steepest descent at x. 3) Perform a line search and find a descent step x + = x + t d. 4) Replace x by x + and go back to step 1. The drawback of this approach is that it may fail to converge due to the nonsmoothness of f . We believe that a descent method should at least give the weak convergence certificate that accumulation points of the sequence of iterates are critical. This is not guaranteed by the above scheme. The reason is that the steepest descent direction at x does not depend continuously on x. This is why modifications of the steepest descent scheme are discussed in the next section.

Remark. Spectral abscissa versus H ∞ norm. Notice that the tangent program ( 16) is convenient because it leads to relatively small size SDPs. Indeed, the matrices Y ν , ν = 1, . . . , q, are of the size of the multiplicity of λ 1 (F(x, jω ν )), and our experiments indicate that dim(y) in ( 16) rarely exceeds 30. The situation is very different for the spectral abscissa [START_REF] Thevenet | Nonlinear spectral SDP method for BMI-constrained problems: Applications to control design[END_REF]. In [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF] we have derived a tangent program for [START_REF] Thevenet | Nonlinear spectral SDP method for BMI-constrained problems: Applications to control design[END_REF]. The difficulty is that Lyapunov variables re-enter the scene. Indeed, x * is a local minimum of the composite function α • F like in [START_REF] Thevenet | Nonlinear spectral SDP method for BMI-constrained problems: Applications to control design[END_REF] with value α (F(x * )) = t * if and only if (x * , t * , X * ) is a local minimum of the optimization program:

(P ) minimize t subject to X σI, X I F(x) T X + XF(x) -2tX 0
for some small fixed 0 < σ ≪ 1. An optimality test for α • F is therefore derived from an optimality test for program (P ), as shown in [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF]. This leads to a SDP with unknown variable of dim(x, X), which may be prohibitively large. This is one of the reasons why our present approach privileges the use of the H ∞ norm (5) over [START_REF] Thevenet | Nonlinear spectral SDP method for BMI-constrained problems: Applications to control design[END_REF].

VI. FIRST-ORDER DESCENT METHOD In this section we devise a first-order algorithm for composite functions of the H ∞ norm. Along with G : R n → H ∞ we consider the symmetrized operators

F(x, s) = G(x, s)G(x, s) H respectively F(x, s) = G(x, s) H G(x, s). We represent f (x) as f (x) = sup ω∈R f (x, ω), f (x, ω) = λ 1 (F(x, jω)) , (17)
and solve the optimization program

min x∈R n f (x), f (x) = G(x) 2 ∞ = sup ω∈R λ 1 (F(x, jω))
Notice that for fixed ω ∈ R, x → F (x, jω) is a smooth operator into the space of Hermitian m × m matrices, S m , while λ 1 : S m → R is the maximum eigenvalue function. Similar techniques could be applied to broader classes with a structure like [START_REF] Apkarian | Fixed-order H∞ control design via an augmented Lagrangian method[END_REF]. Deriving the method will require three steps, which we regroup into subsections. We start with the important special case f = λ 1 •F, where F maps R n smoothly into S m .

A. Preparation

The function ( 17) is subject to two sources of nonsmoothness. The nonsmooth character of the maximum eigenvalue function, and the nonsmoothness introduced by the operator sup, which in the case of • ∞ is even infinite. Each individual function f (•, ω) will be analytic at x if the multiplicity of λ 1 (F(x, jω)) is one, but nonsmoothness needs to be taken into account as soon as eigenvalues coalesce.

From a practical point of view it is reasonable to make the following additional hypothesis. (H) The maximum f (x) is always attained on a finite set of frequencies. This set is denoted by Ω(x) and may contain ω = ±∞. Assumption (H) is for instance satisfied when the multiplicity of λ 1 is 1, as ω → f (x, ω) is then analytic in typical control applications.

Let us introduce some more notation. For

Ω ⊂ R ∪ {-∞, +∞} we define f Ω ≤ f as f Ω (x) = max ω∈Ω f (x, ω).
Notice that f Ω (x) = f (x) as soon as Ω(x) ⊂ Ω. Next recall the definition of the ǫ-subdifferential of the maximum eigenvalue function [START_REF] Cullum | The minimization of certain nondifferentiable sums of eigenvalues of symmetric matrices[END_REF] ∂

ǫ λ 1 (X) = {Z ∈ S m : Tr (Z X) ≥ λ 1 (X) -ǫ}
which is an important analytical tool in nonsmooth analysis. Since ∂ ǫ λ 1 (X) is difficult to compute, we follow Cullum et al. [START_REF] Cullum | The minimization of certain nondifferentiable sums of eigenvalues of symmetric matrices[END_REF] and Oustry [START_REF] Oustry | A second-order bundle method to minimize the maximum eigenvalue function[END_REF] and introduce a modification δ ǫ λ 1 (X) of ∂ ǫ λ 1 (X), called the ǫ-enlarged subdifferential for the maximum eigenvalue function. For ǫ > 0 and X ∈ S m let r(ǫ, X) the index such that

λ 1 ≥ . . . ≥ λ r(ǫ,X) ≥ λ 1 -ǫ > λ r(ǫ,X)+1 ≥ . . . ≥ λ m .
The index r(ǫ, X) is also called the ǫ-multiplicity of λ 1 (X). Let Q ǫ be a r(ǫ, X) × m-matrix whose columns form an orthonormal basis of the invariant subspace of X associated with the first r(ǫ, X) eigenvalues. Then we define

δ ǫ λ 1 (X) = {Q ǫ Y Q H ǫ : Y ∈ B r(ǫ,X) }. By construction ∂λ 1 (X) ⊂ δ ǫ λ 1 (X) ⊂ ∂ ǫ λ 1 (X), so δ ǫ λ 1 (X)
is an enlargement of ∂λ 1 (X) and an inner approximation of the ǫ-subdifferential (see also [START_REF] Oustry | A second-order bundle method to minimize the maximum eigenvalue function[END_REF]). The gap associated with the choice

ǫ is ∆(ǫ, X) = λ r(ǫ,X) (X) -λ r(ǫ,X)+1 (X) > 0. If
r is the multiplicity of λ 1 (X), then choosing ǫ small enough gives r(ǫ, X) = r. In this case we have δ ǫ λ 1 (X) = ∂λ 1 (X).

The following is an important step toward the analysis of [START_REF] Apkarian | Fixed-order H∞ control design via an augmented Lagrangian method[END_REF]. Consider a differentiable mapping F : R n → S m . We extend ∂ ǫ λ 1 (X) and the enlarged subdifferential δ ǫ λ 1 (X) to the composite function f = λ 1 • F by setting

∂ ǫ f (x) = F ′ (x) ⋆ [∂ ǫ λ 1 (X)] , X = F(x),
and similarly

δ ǫ f (x) = F ′ (x) ⋆ [δ ǫ λ 1 (X)] , X = F(x).
Here F ′ (x) ⋆ is the adjoint of the linear operator F ′ (x). Finally, going one step further, this approach allows us to consider ∂ ǫ f (x, ω) and δ ǫ f (x, ω), which are applied to the variable x.

B. Descent Step Generator

In this section we discuss a very simple mechanism which generates descent steps in such a way that the following weak form of convergence can be guaranteed: every accumulation point of the sequence of iterates is a stationary point.

Let f : R n → R be a locally Lipschitz function, and let ∂f (x) denote its Clarke subdifferential [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]. Suppose we can exhibit a mechanism s : R n → R n , the descent step generator, such that the following rules are satisfied:

(i) Whenever 0 ∈ ∂f (x), then f (s(x)) < f (x).
(ii) When 0 ∈ ∂f (x), then there exists a neighborhood B(x, ǫ) of x and some δ > 0 such that f (s(x ′ )) ≤ f (x ′ )δ for every x ′ ∈ B(x, ǫ).

While (i) simply means that s(x) is a descent step away from x, we can interpret (ii) as some weak form of continuity of s(•). Indeed, when the mapping s(•) describing descent is continuous, then (i) implies (ii) without further work. Axiom (ii) is weaker than asking s(•) to be continuous. Clearly, (ii) always implies (i).

Example. In order to understand the idea behind axiom (ii), consider a C 1 -function f and let s(x) denote the steepest descent step at x, obtained by an Armijo line search. If we define formally s(x) = xt(x)f ′ (x), where t(x) is the smallest step satisfying an Armijo rule, then s(x) turns out continuous, but it is clear that in practice we would accept any step t(x) satisfying the Armijo condition, without insisting on a continuous dependence of t(x) on x. In that case s(x) will not be continuous, but property (ii) will still hold true.

Clearly the situation we have in mind is when f is nonsmooth, so that a steepest descent step in tandem with Armijo search would typically fail even when t(x) was continuous. This is explained by the fact that under nonsmoothness, defining s(•) along the lines above would miss axiom (ii), because the steepest descent direction -f ′ (x) behaves very discontinuously. Indeed, examples where this happens are easily produced.

Proposition 1: . Suppose the descent step generator s(•) for f satisfies axioms (i) and (ii). Let x k be a sequence of iterates defined as x k+1 = s(x k ). Then every accumulation point x of x k is a critical point, that is, satisfies 0 ∈ ∂f (x).

Proof: Let N ⊂ N be an infinite sequence such that

x k → x, k ∈ N . Then by monotonicity, f (x k+1 ) → f (x). That means f (x k+1 ) -f (x k ) → 0 (k ∈ N ).
Now use axiom (ii) at the limit point x. There exist ǫ, δ > 0 such that f (s(x))f (x) ≤ -δ < 0 for all x ∈ B(x, ǫ). Since x k ∈ B(x, ǫ) for k ∈ N large enough, and since

x k+1 = s(x k ), we should have f (x k+1 ) -f (x k ) ≤ -δ for k ∈ N large enough, a contradiction.

C. Eigenvalue Optimization

How can we define a descent step generator s(•) with properties (i) and (ii) for a maximum eigenvalue function f (x) = λ 1 (F(x))? Suppose we are at a point x where 0 ∈ ∂f (x). Let the eigenvalues of X = F(x) ∈ S m be arranged into groups:

λ 1 = . . . = λ k2-1 > λ k2 = . . . = λ k3-1 > λ k3 = . . .
where k 1 = 1 and where k i are the group leaders. Consequently, eigenvalue gaps occur between k i -1 and k i . Let Q 1 be an orthonormal bases of the eigenspace associated with the first block λ 1 (X), . . . , λ k2-1 (X), Q 2 an orthonormal basis containing Q 1 associated with the first two blocks λ 1 (X), . . . λ k3-1 (X), and so on. At X = F(x) ∈ S m we compute the quantities

∆ i (X) = λ ki+1-1 (X) -λ ki+1 (X) > 0, d i (x) = min{ F ′ (x) ⋆ Q i Y i Q H i : Y i ∈ B ki+1-1 }
and keep those i = 1, . . . , r where d i (x) > 0. Notice that

d i (x) = dist(0, δ ǫi(x) f (x)), where ǫ i (x) > 0 cuts into the ith gap, that is λ 1 (X) -ǫ i (x) ∈ [λ ki+1-1 (X), λ ki+1 (X)). Put differently, the ǫ i (x)-multiplicity of λ 1 (X) is k i+1 -1. Moreover, d 1 (x) ≥ d 2 (x) ≥ . . . ≥ d r (x) and d r+1 (x) = d r+2 (x) = .
. . = 0 eventually. We compute the quantity

M (x) = max i=1...r ∆ i (X)d i (x) 2 (18)
Now we use the following Lemma 1:. Let 0 ∈ ∂f (x) and R > 0. Let k i be the leader of a group of eigenvalues of X = F(x) ∈ S m such that

d i (x) = dist(0, δ ǫi(x) f (x)) > 0. Let h i (x) be the direction of steepest ǫ i (x) enlarged descent, that is h i (x) = -gi gi with g i = argmin{ g : g = F ′ (x) ⋆ Q i Y Q H i , Y ∈ B ki+1-1 },
where the columns of Q i are an orthonormal basis of the invariant subspace of X associated with the eigenvalues up to k i+1 -1. Then there exists a descent step s i (x) away from x in direction h i (x), which decreases the value of f by at least

f (s i (x)) -f (x) ≤ -κ(x)∆ i (X)d i (x) 2 < 0.
Here κ(x) > 0 depends on sup{ F ′ (x ′ ) : xx ′ ≤ R} and in particular continuously on x. The line search required to compute this step is finite.

Proof: In the case of an affine operator A : R n → S m , Oustry [START_REF] Oustry | A second-order bundle method to minimize the maximum eigenvalue function[END_REF]Theorem 5] shows that f = λ 1 • A may be decreased by the quantity

f (x + t(x)h i (x)) -f (x) ≤ - 1 4 A * ∆ ǫi(x) (X)d i (x) 2
where A * is the linear part of A. Here κ = 1/4 A * is even independent of x. Moreover, t(x) > 0 is computed by a line search which terminates after a finite number of steps, see [25, sect. 3.2].

In [START_REF] Noll | Spectral bundle methods for nonconvex maximum eigenvalue functions: first-order methods[END_REF], [START_REF]Spectral bundle methods for nonconvex maximum eigenvalue functions: second-order methods[END_REF], this result is generalized to nonconvex maximum eigenvalue functions f = λ 1 • F, with a constant κ(x) now depending on the Lipschitz constant of F on a bounded region around x, like for instance B(x, R) with some fixed R > 0. In the nonconvex case, the line search procedure locating t(x) is more complicated than for f = λ 1 • A, but finite termination is still guaranteed (cf. [37, sect. 3.7]). The constant κ(x) may be computed via the formulae ( 15), ( 19) and ( 20) of that reference.

Following the lead of [START_REF] Noll | Spectral bundle methods for nonconvex maximum eigenvalue functions: first-order methods[END_REF], it is now clear how to obtain our step s(x). Choose i so that the maximum ( 18) is attained and take s(x) = s i (x). Then by Lemma 1, s(x) gives a guaranteed decrease of

f (s(x)) -f (x) ≤ -κ(x)∆ ǫi(x) (X)d i (x) 2 = -κ(x)M (x),
where the constant κ(x) > 0 depends continuously on x as argued above. What remains to be checked is Lemma 2:. This choice of s(x) and M (x) guarantees property (ii).

Proof: Consider a sequence x k → x. We need to compare M (x k ) to M (x). Now observe that due to the continuity of the eigenvalue functions λ j (X), every eigenvalue gap at X is also an eigenvalue gap at X k as soon as X k is sufficiently close to X. On the other hand, X k may (and will) have many more gaps than X. But notice that the maximum ( 18) is over all eigenvalue gaps. So each gap in M (x) will occur in the computation of M (x k ). More precisely, it will be approximated by some of the gaps considered in M (x k ). Put differently, for the ith eigenvalue gap of X we have 2 , where i k is the index of the eigenvalue gap of X k corresponding to the ith gap of X.

∆ i k (X k )d i k (x k ) 2 → ∆ i (X)d i (x)
Here ∆ i k (X k ) → ∆ i (X) and d i k (x k ) → d i (x)
rely on the fact that gaps at X k remain gaps at X. Since d i k (x k ) calls for all the Q 1 , Q 2 , . . . up to the i k th gap at X k , and since these converge to the corresponding basis Q regrouping the gap i, the result follows.

That means lim sup M (x k ) ≥ M (x), so we can assume M (x k ) ≥ 1 2 M (x) from some index k on. This in turn implies

f (s(x k )) -f (x k ) ≤ -κ(x k )M (x k ) ≤ - 1 4 κ(x)M (x)
as soon as κ(x k ) is close enough to κ(x), proving property (ii).

The method outlined in this section gives a convergence certificate because all eigenvalue gaps are included in the computation of [START_REF] Noll | Partially augmented Lagrangian method for matrix inequality constraints[END_REF]. This may seem inconvenient for very large size matrices F(x). If we decide to truncate and consider only some of the eigenvalue gaps among the largest eigenvalues, the theoretical convergence properties of the method are weaker, even though convergence may still be guaranteed e.g. when f is convex (see [START_REF] Noll | Spectral bundle methods for nonconvex maximum eigenvalue functions: first-order methods[END_REF]). Notice however that even for large m the quantity M (x) in ( 18) may be computed fairly reliably by considering ∆ ǫi(x) (X)d i (x) 2 for some of the first i only. As d i (x) → 0 rather quickly, the higher ∆ ǫi(x) (X)d i (x) 2 as a rule do not contribute to the computation of M (x). Also notice that since the sequence d i (x) is monotone, computing d i (x) may often be avoided, for instance when ∆ ǫi-1(x) (X) ≥ ∆ ǫi(x) (X) or ∆ ǫi(x) (X) ≤ ∆ ǫi-ν (x) (X) when the current best value is located at index iν.

D. Descent by a local model

In this section we present an alternative way to obtain a descent step generator s(•) for the maximum eigenvalue function f (x) = λ 1 (F(x)). We start by constructing an intermediate function θ(x), which serves as an optimality test, and secondly, will allow us to quantify decrease. In this section our method follows the line of [START_REF]Optimization : Algorithms and Consistent Approximations[END_REF]Thm. 2.1.6].

Let X = F(x) ∈ S m . Let µ 1 (X) > µ 2 (X) > . . . > µ r (X) be the eigenvalues of X without repetitions. That means µ i = λ ki in our old terminology, where we agree that there are r ≤ m distinct eigenvalues. For some fixed σ > 0 we define the criticality measure

θ(x) = inf h∈R n sup i=1,...,r sup Yi∈B {-f (x) + µ i (X) +Tr (Y i Q H i [F ′ (x)h]Q i ) + 1 2 σ h 2 (19)
Here Q i has k i+1 -1 columns which form an orthonormal basis of the invariant subspace of X associated with the first k i+1 -1 eigenvalues. It is immediately clear that θ(x) ≤ 0, because putting h = 0 gives the upper bound θ(x) ≤ sup i=1,...,r -f (x) + µ i (X) = -f (x) + µ 1 (X) = 0. Lemma 3:. We have θ(x) = 0 if and only if 0 ∈ ∂f (x).

Proof: The easiest way to see this is to swap max and min in [START_REF] Xin | Reduced-order controllers for the H∞ control problem with unstable invariant zeros[END_REF]. This requires that we first replace the inner double supremum in [START_REF] Xin | Reduced-order controllers for the H∞ control problem with unstable invariant zeros[END_REF] by a double supremum over the convex hull of Y i 0, Tr(Y i ) = 1, i = 1, . . . , r, a manoeuvre which does not change the value θ(x). Then we use Fenchel duality to interchange the inner (double) supremum and the outer infimum, which goes again without changing the value. The now inner infimum is unconstrained and may be computed explicitly. For fixed Y i and convex coefficients τ it is attained at h

(x) = -σ -1 r i=1 τ i F ′ (x) ⋆ Q i Y i Q H i .
Substituting this back into [START_REF] Xin | Reduced-order controllers for the H∞ control problem with unstable invariant zeros[END_REF] leaves the dual expression

θ(x) = sup τi≥0, P i τi=1 sup Yi∈B -f (x) + r i=1 τ i µ i (X) - 1 2σ r i=1 τ i F ′ (x) ⋆ Q i Y i Q H i 2    ( 20
)
which the reader recognizes as a semidefinite program. Since µ i (X) < f (x) for i ≥ 2, equality θ(x) = 0 is only possible when τ 2 = . . . = τ r = 0 and hence

0 = θ(x) = sup Y1∈B - 1 2σ F ′ (x) ⋆ Q 1 Y 1 Q H 1
But the quantity on the right hand side is only zero when

F ′ (x) ⋆ Q 1 Y 1 Q H 1 = 0, i.e.
, when 0 ∈ ∂f (x). The latter follows readily from the representation

∂f (x) = {F ′ (x) ⋆ Q 1 Y 1 Q H 1 : Y 1 ∈ B} of the subdifferential ∂f (x).
As a byproduct of the proof via duality we have the following Corollary 1: The infimum in ( 19) is attained at

h(x) = - 1 σ r i=1 τ i F ′ (x) ⋆ Q i Y i Q H i , (21)
where (τ, Y ) is the solution of the dual program [START_REF] Xin | Reduced-order controllers for continuous and discrete-time singular H∞ control problems based on LMI[END_REF]. As soon as θ(x) < 0, h(x) is a direction of descent of f at x.

In order to construct our descent step generator, we need to establish two additional properties of θ. Firstly, we need to show that decrease at x may be quantified with the help of θ. Secondly, we have to establish continuity of θ.

Lemma 4: . The function θ is continuous.

Proof: Notice that by the dual formula [START_REF] Xin | Reduced-order controllers for continuous and discrete-time singular H∞ control problems based on LMI[END_REF], θ(x) is the supremum of an infinite family of functions of the form

-f (x) + r i=1 τ i µ i (X) + c(x) = -f (x) + m j=1 σ j λ j (X) + c(x) for σ j = τ i /(k i+1 -k i ), where k i ≤ j < k i+1 , and where c(x) = -1/2σ i τ i F ′ (x) ⋆ Q i Y i Q H i 2 de- pends continuously on x. Notice that µ i (X) = (k i+1 - k i ) -1 ki+1-1 j=ki λ j (X)
. This shows that θ(x) is the supremum of a family of continuous functions, indexed by (τ, Y ). θ is therefore lower semi-continuous. It remains to prove that θ is also upper semi-continuous. Let x j → x such that θ(x j ) → θ. We have to show θ ≤ θ(x). We use the representation [START_REF] Xin | Reduced-order controllers for continuous and discrete-time singular H∞ control problems based on LMI[END_REF]. Let ǫ j > 0, ǫ j → 0 and choose τ j i and Y j i such that

θ(x j ) ≤ -f (x j ) + r i=1 τ j i µ i (X j ) - 1 2σ r i=1 τ j i g j i 2 + ǫ j ,
where

g j i = F ′ (x j ) ⋆ Q j i Y j i Q jH i .
Passing to a subsequence if necessary, we may assume that each X j has exactly r eigenvalue gaps, which remain in the same places k i , i = 1, . . . , r. That is µ i (X j ) = λ ki (X j ). Passing to yet another subsequence, we may assume

τ j → τ , Y j i → Y i and Q j i → Q i ,
where the limiting elements are all of the same types and dimensions as the elements at stage j.

However, λ ki (X) are no longer the distinct eigenvalues of X, because some of the distinct µ i (X j ) = λ ki (X j ) may coalesce in the limit j → ∞. Suppose for instance that lim j→∞ µ i (X j ) = . . . = lim j→∞ µ i+t (X j ), so that the blocks i, i + 1, . . . , i + t coalesce in the limit, forming a new larger eigenvalue block of X:

λ ki-1 (X) > λ ki (X) = . . . = λ ki+t+1-1 (X) > λ ki+t (X),
which is represented by a certain µ ν (X). Suppose there are N block leaders at X. For each of these µ ν (X) we define

σ ν = lim j→∞ i+t s=i τ j s , then lim j→∞ r i=1 τ j i µ i (X j ) = N ν=1 σ ν µ ν (X),
where N ν=1 σ ν = 1. This represents the limiting linear term in [START_REF] Xin | Reduced-order controllers for continuous and discrete-time singular H∞ control problems based on LMI[END_REF] in a new form suited for the dual representation of θ(x).

We next have to treat the norm square term arising in [START_REF] Xin | Reduced-order controllers for continuous and discrete-time singular H∞ control problems based on LMI[END_REF] in much the same way. Notice first that the Q j i , i = 1, . . . , r together form a nested sequence of basis vectors adding up to an orthonormal basis of eigenvectors of X j . Passing to the limit j → ∞ gives an orthonormal basis of eigenvectors of X. We regroup it according to the eigenvalue gaps of X and rename the corresponding parts P 1 ⊂ P 2 ⊂ . . . P N . All that remains to do is to re-write the limit lim j→∞

r i=1 τ j i Q j i Y j i Q jH i in the form N ν=1 σ ν P ν Z ν P H ν for certain Z ν ∈ B.
This is done by writing

Q i Y i Q H i = Q j Y i 0 0 0 Q H j = P ν Y i 0 0 0 P H ν
whenever i < j have to be regrouped in the same P ν , and where j is the last among the old indices subsumed into the new index ν, so that Q j = P ν . Then

Z ν = i+t s=i τ i σ ν Y i 0 0 0 is as required, because i+t s=i τ s = σ ν , hence Tr(Z ν ) = 1, while Z ν 0 is clear.
The argument shows that θ = -f (x) + N ν=1 σ ν µ ν (X) -

1 2σ N ν=1 σ ν F ′ (x) ⋆ P ν Z ν P H ν 2
, which is to say that θ is now of the form required to be admitted to the supremum [START_REF] Xin | Reduced-order controllers for continuous and discrete-time singular H∞ control problems based on LMI[END_REF] defining θ(x). In other words, θ ≤ θ(x), and this is what we had to prove.

Notice that [36, Thm. 2.1.6 (e)] is obtained as a special case of Lemma 4 if the operator F(x) is specialized to a diagonal matrix. The extension to multiple eigenvalues is possible because all the eigenvalues are taken into account simultaneously. For large matrices, this may again seem inconvenient since it will lead to large SDPs in [START_REF] Xin | Reduced-order controllers for continuous and discrete-time singular H∞ control problems based on LMI[END_REF].

Lemma 5: . The mapping h(x) defined by ( 21) is continuous.

The proof follows along the lines of the previous Lemma and is therefore omitted.

Let us now see how θ(x) may be used to quantify descent of f = λ 1 • F at x. Assume 0 ∈ ∂f (x), so that θ(x) < 0. Using the directional derivative of f at x in direction h(x), we obtain

f ′ (x; h(x)) = sup Y ∈B h(x) T F ′ (x) ⋆ Q 1 Y Q H 1 = sup Y ∈B Tr (Y Q H 1 [F ′ (x)h(x)]Q 1 ) ≤ θ(x) - 1 2σ h(x) 2 < θ(x) < 0
which follows readily from the primal formula [START_REF] Xin | Reduced-order controllers for the H∞ control problem with unstable invariant zeros[END_REF] for θ if we use the fact that -f (x) + µ 1 (X) = 0. In consequence we have the following Lemma 6: . Let 0 < τ < 1 be fixed. There exists ǫ > 0 such that f (x ′ + th(x ′ ))f (x ′ ) < tτ θ(x ′ ) for every x ′ ∈ B(x, ǫ) and all t > 0 sufficiently small.

Proof: By the definition of the directional derivative and the fact that 0 < τ < 1 we have

f (x + th(x)) -f (x) ≤ tτ f ′ (x; h(x)) ≤ tτ θ(x) - 1 2σ h(x) 2 (22)
< tτ θ(x)

for some t 0 > 0 and all 0 < t < t 0 . Since h(•) and θ(•) are continuous, we can find a neighborhood B(x, ǫ) of x such that f (x ′ + th(x ′ ))f (x ′ ) < tτ θ(x ′ ) for all x ′ ∈ B(x, ǫ) and every 0 < t < t 0 . This proves the claim.

In order to construct s(•), we follow [36, 3c, p. 223] and define s(x) = x + t(x)h(x), where

t(x) := sup{2 -k : k ∈ N, f (x+2 -k h(x))-f (x) < 2 -k τ θ(x)}.
The supremum is over a nonempty set because of ( 22), hence 0 < t(x) < +∞. Let k(x) be the integer where the supremum t(x) is attained. Let us check property (ii) with ǫ > 0 as in Lemma 6 and δ

= -2 -k(x)-1 τ θ(x) > 0. Let x ′ ∈ B(x, ǫ). Since t(x) satisfies (22), we have f (x ′ + t(x)h(x ′ )) -f (x ′ ) < t(x)τ θ(x ′ ) by Lemma 6. Therefore t(x) = 2 -k(x) is admitted for the supremum t(x ′ ), which implies t(x ′ ) ≥ t(x). Therefore f (s(x ′ )) -f (x ′ ) = f (x ′ + t(x ′ )h(x ′ )) -f (x ′ ) < t(x ′ )τ θ(x ′ ) ≤ t(x)τ θ(x ′ ) ≤ t(x)τ θ(x)/2 = -2 -k(x)-1 τ θ(x) = -δ,
when we assume that ǫ is chosen sufficiently small to assure θ(x ′ ) ≤ θ(x)/2 for every x ′ ∈ B(x, ǫ). This proves (ii) for the above choice of t(x). Other choices of t(x) like for instance [START_REF] Pironneau | On the rate of convergence of a class of methods of centers[END_REF] are possible.

E. The semi-infinite case

Our last step is now to address the semi-infinite case. We are in the situation [START_REF] Apkarian | Fixed-order H∞ control design via an augmented Lagrangian method[END_REF]. Suppose that for finite Ω ⊂ R ∪ {-∞, +∞} we already dispose of a descent step generator s Ω for f Ω , satisfying axioms (i) and (ii). This is naturally the case when the f (x, ω) are maximum eigenvalue functions, because a finite maximum of maximum eigenvalue functions is itself a maximum eigenvalue function. So here we obtain s Ω as in sections VI-C or VI-D. Suppose now that we can specify a sequence of finite sets Ω 1 ⊂ Ω 2 ⊂ . . . such that the following conditions are satisfied:

(iii) If 0 ∈ ∂f (x), then lim sup k→∞ f (s Ω k (x)) -f (x) < 0. (iv) For every x ∈ R n and ǫ > 0, lim k→∞ max x ′ ∈B(x,ǫ) inf ω ′ ∈Ω k |f (x ′ , ω) -f (x ′ , ω ′ )| = 0 .
Notice that both axioms guarantee that the approximation Ω k improves with growing k. Axiom (iii) tells that as soon as 0 ∈ ∂f (x), descent steps may be eventually generated by using approximations f Ω k of f . Axiom (iv) is simply saying that approximations get better as k increases, and that this happens uniformly on small neighborhoods of each x.

Using these axioms, let us construct a descent step generator s(•) for f . We proceed as follows:

Semi-infinite descent step generator 1) If 0 ∈ ∂f (x) then s(x) = x and return. Otherwise put counter k = 1 and continue.

2) At counter k, if 0 ∈ ∂f Ω k (x), then increase k until 0 ∈ ∂f Ω k (x). 3) At counter k with 0 ∈ ∂f Ω k (x), compute the descent step s Ω k (x) for f Ω k at x. Let ǫ k , δ k > 0 be such that f Ω k (s k (x ′ )) -f Ω k (x ′ ) ≤ -δ k < 0 for every x ′ ∈ B(x, ǫ k ), as guaranteed by axiom (ii) for s Ω k . 4) Compute η k := inf{|f (x ′ , ω) -f (x ′ , ω ′ )| : ω ′ ∈ Ω k , x ′ ∈ B(x, ǫ k ) ∪ s (B(x, ǫ k ))} 5) If 3η k < δ k , then let s(x) = s Ω k (x) and stop.
Otherwise increase k by one and go back to step 3. We have to make sure that this scheme is well-defined and introduces a step generator s(•) for the infinite maximum function f . Lemma 7: . The descent step generator s(•) for f (x) = sup ω∈R f (x, ω) is well-defined and satisfies axioms (i) and (ii). Moreover, each of the above loops ends after a finite number of iterations.

Proof: Notice first that if 0 ∈ ∂f (x), then descent around x is possible. Since f Ω k (x) → f (x) as k → ∞, we can also decrease f Ω k around x for k sufficiently large. So step 2 ends with a descent step of f Ω k at x after a finite number of trials. Moreover, this remains so for the following counters k, because Ω k ⊂ Ω k+1 .

Next observe that η k → 0 by axiom (iv), while lim sup -δ k < 0 by axiom (iii). That means 3η k < δ k for k sufficiently large, i.e., the procedure ends in step 4 after a finite number of updates k → k + 1. It remains to check that s is as required.

Suppose

f Ω k (s Ω k (x ′ )) -f Ω k (x ′ ) ≤ -δ k for every x ′ ∈ B(x, ǫ k ). Then |f (s Ω k (x ′ )) -f Ω k (s Ω k (x ′ ))| ≤ η k , |f (x ′ ) - f Ω k (x ′ )| ≤ η k , hence f (s Ω k (x ′ )) -f (x ′ ) ≤ -δ k + 2η k ≤ -1 3 δ k . This proves axiom (ii).
The procedure is sufficiently flexible to accommodate a problem oriented step generation. We may adapt the choice of Ω k to the structure of f , and what is more important, to the local behavior of f around the current x.

F. First-order algorithm for the H ∞ norm

We are now ready to present our algorithm, which follows the lines of the previous sections. We start with a version based on sections VI-C and VI-E.

First-order algorithm for the H ∞ norm: variant I 0) Fix 0 < τ < 1 and choose initial point x 0 . 1) Given x k choose a finite set Ω k containing Ω(x k ).

2) Compute M (x k ) for f Ω k according to [START_REF] Noll | Partially augmented Lagrangian method for matrix inequality constraints[END_REF]. If M (x k ) = 0 then stop, because 0 ∈ ∂f (x). Otherwise 3) Use a line search to find a step t k such that the predicted decrease satisfies

π k = f Ω k (x k + t k h k ) -f (x k ) ≤ -κ(x k )M (x k ) < 0. Compare with the actual decrease α k = f (x k + t k h k ) -f (x k ). If α k ≤ τ π k , accept t k , put x k+1 = x k + t k h k and goto step 5.
4) If α k > τ π k then reject t k and add nodes to Ω k to obtain the finer mesh Ω k+1 . Increase counter k by one and go back to step 2. 5) Increase counter k by one and go back to step 1. In steps 3 and 4 of the algorithm we recognize the mechanism of the previous section, which obtains a descent step generator for the semi-infinite function by using those of the finite models f Ω k . We accept the step for f if it exceeds a small fraction of the descent predicted by the finite model f Ω k . This is essentially the same procedure as in VI-E.

Next comes the version based on section VI-D in tandem with section VI-E. The semi-infinite case is handled in exactly the same fashion, but the descent step generators s Ω k are different.

First-order algorithm for the H ∞ norm: variant II 0) Fix 0 < τ 1 , τ 2 < 1 and δ > 0 and choose initial point x 0 . 1) Given x k choose a finite set Ω k containing Ω(x k ).

2) Compute the value θ k and the solution (τ k , Y k ) of the SDP [START_REF] Xin | Reduced-order controllers for the H∞ control problem with unstable invariant zeros[END_REF]. If θ k = 0, then stop because 0 ∈ ∂f (x k ).

Otherwise compute the descent direction h k for f Ω k at x k according to [START_REF] Stoorvogel | A reduced-order oberver based controller design for H∞-optimization[END_REF]. 3) Using a line search, find a step t k such that the predicted decrease satisfies

π k = f Ω k (x k + t k h k ) -f (x k ) < τ 1 t k θ k < 0. Compare with the actual decrease α k = f (x k + t k h k ) -f (x k ). If α k ≤ τ 2 π k , accept t k , put x k+1 = x k + t k h k

G. Practical aspects

In this section we comment on the salient features of the nonsmooth first-order descent algorithm and address some of the practical aspects.

In our testing we have observed that often the leading eigenvalues λ 1 (F(x k , jω ν )) have multiplicity 1 at all frequencies ω ν ∈ Ω(x k ). In this situation step 2 of the descent algorithm variant II could be simplified. Suppose at the current iterate x k we have selected a finite set of frequencies

Ω k = {ω 1 , . . . , ω q } containing Ω(x k ). Suppose all λ 1 (F(x k , jω ν )) have multi- plicity 1. Let g ν = F ′ (x k , jω ν ) ⋆ e ν e T ν = f ′ (x k ; ω ν )
, where e ν is the normalized eigenvector associated with λ 1 (F(x k , jω ν )). Then the semidefinite program ( 16) respectively [START_REF] Xin | Reduced-order controllers for continuous and discrete-time singular H∞ control problems based on LMI[END_REF] simplifies to a convex quadratic program θ(x k ) = sup τν ≥0,

P ν τν =1 -f Ω (x k ) + q ν=1 τ ν f (x k , ω ν ) - 1 2σ q ν=1 τ ν g ν 2
and the associated direction of descent

h k = h(x k ) is h(x k ) = - 1 σ q ν=1 τ ν g ν ,
where τ is the optimal solution of the quadratic program.

Observe that for Ω k = Ω(x k ), h(x k ) coincides with the steepest descent direction for a finite max function.

Let us now specify in which way we select the frequency set Ω k at each step. The finite set of frequencies Ω(x k ) where the H ∞ norm is attained is computed via the Hamiltonian technique [START_REF] Boyd | A bisection method for computing the H∞ norm of a transfer matrix and related problems[END_REF]. We then form an enriched set Ω k of frequencies by adding to the peak frequencies a collection of logarithmically spaced frequencies ω ν such that

T w→z (K) ∞ -σ(T w→z (K, jω ν )) ≤ ǫ ω T w→z (K) ∞ ,
where ǫ ω is a user-specified tolerance. We usually limit the set to the first 50 frequencies with largest singular values, as this appears to work well on a broad range of numerical tests. Typical values for ǫ ω range from 0.05 to 0.5. The algorithm requires that this set be iteratively refined when descent steps cannot be computed, but in practice our choice is usually satisfactory, and numerical problems due to exceedingly fine Ω k can be avoided.

VII. NUMERICAL EXPERIMENTS

In this section we test our nonsmooth algorithms on a variety of synthesis problems from the COM P l e ib collection by F. Leibfritz [START_REF] Leibfritz | COMPLeIB, COnstraint Matrix-optimization Problem LIbrary -a collection of test examples for nonlinear semidefinite programs, control system design and related problems[END_REF]. Computations were performed on a (low-level) SUN-Blade Sparc with 256 RAM and a 650 MHz sparcv9 processor. LMI-related computations for search directions used the LMI Control Toolbox [START_REF] Gahinet | LMI Control Toolbox[END_REF] or our home made SDP code [START_REF] Apkarian | A spectral quadratic-SDP method with applications to fixed-order H 2 and H∞ synthesis[END_REF] while QP computations are based on Schittkowski's code [START_REF] Schittkowski | QLD: A Fortran code for quadratic programming[END_REF].

Our algorithm is a first-order method. Not surprisingly, it may be slow in the neighborhood of a local solution. We have implemented various stopping criteria to ensure that an adequate approximation of a solution has been found and to avoid unwarranted computational efforts as is often the case with a first-order algorithm. The first of these termination criteria is an absolute stopping test, which provides a criticality assessment inf{ g :

g ∈ ∂f (x)} < ǫ 1 , (23)
This is reasonable, as 0 ∈ ∂f (x) indicates a critical point. It is also mandatory to use relative stopping criteria to reduce the dependence on the problem scaling. The test

T w→z (K) ∞ -T w→z (K + ) ∞ < ǫ 2 (1 + T w→z (K) ∞ ) , (24)
compares the progress achieved relatively to the current H ∞ performance, while

K + -K < ǫ 3 (1 + K ) (25)
compares the step-length to the controller gains. The tolerances

ǫ 1 = 1e-5, ǫ 2 = 1e-3, ǫ 3 = 1e-3
have been used in our numerical testing. For stopping we required that either the first two tests or the third one are satisfied. For the enriched set, the number of frequencies has been limited to 50. They are selected according to our discussion in section VI-G. It is sometimes possible to employ fewer frequencies, but generally better steps are performed when richer sets are used. Our choice appears reasonable and has been validated on numerous experiments. It does not restrict efficiency since QP codes are very efficient up to 500 variables.

A. Stabilization

Pure stabilization may be regarded as H ∞ synthesis under the special form [START_REF] Apkarian | A spectral quadratic-SDP method with applications to fixed-order H 2 and H∞ synthesis[END_REF]. The optimization program is ( 5), but we stop the algorithm as soon as a stabilizing controller is obtained. Iterating until a local optimum of ( 5) is reached does not seem to improve any of the usual performance specifications of the stabilizing feedback controller. Such questions have to be addressed in a second procedure, where the stabilizing controller will serve as a starting point. The spectral abscissa ( 7) is used to check stability of K, but is not used as a cost function. Previous experiments where [START_REF] Thevenet | Nonlinear spectral SDP method for BMI-constrained problems: Applications to control design[END_REF] was used to find stabilizing controllers are reported in [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF].

Table I displays results obtained for static stabilization problems borrowed from the literature. The triple (n, m, p) gives the number of states, inputs and outputs. Column 'iter' corresponds to the number of iterations required to meet the combined stopping tests. Column α displays the final closed-loop spectral abscissa [START_REF] Thevenet | Nonlinear spectral SDP method for BMI-constrained problems: Applications to control design[END_REF]. Negative α < 0 indicates a stable closed-loop system. Column 'cpu' gives the cputime in seconds. Note that the reported cputimes are only indicative as the α-version of our code includes a number of extra tests and even graphics that have been exploited for fine tuning of the various algorithm parameters.

The initial shift for the H ∞ norm minimization was chosen according to the following rule

a 0 := max (α(A) * (1 + 10%), 1e-2) , (26)
with a zero initial static controller K 0 = 0.

From the table, we observe that in all cases but one, the original choice a 0 of the shift was successful. As indicated, there is no need to run the algorithm until convergence. A stabilizing solution is therefore obtained very early. A single SDP or convex QP with line search suffices in most cases.

In the second example, we had to reduce the shift three times before stability was reached. This was done according to rule [START_REF] Anderson | Network analysis and synthesis: a modern systems theory approach[END_REF], applied to the closed-loop dynamics. In this example the H ∞ criterion is flat about a (global) optimum, which allows longer steps and renders the step-length-based termination criterion more stringent. Note that in this case the solution is globally optimal because a zero H ∞ norm is reached. [START_REF] Chen | H∞ Control and Its Applications, ser. Lectures Notes in Control and Information Sciences[END_REF] (5, 1, 3) 1 -9.95e-1 2.88 Boeing 767 [START_REF]Benchmark problems for control system design[END_REF] (55, 2, 2) 1 -2.33e-2 6.47

steepest descent method -ǫω = 0.05

For test reasons, the pure stabilization problems have been solved with the steepest descent method in section V. As expected, this technique is less stable than algorithmic variants I and II. Specifically, the choice of ǫ ω is rather critical and finding a general selection rule appears difficult. Together with our analysis in section VI, this encouraged us to use algorithmic variant II in the syntheses problems below.

B. H ∞ synthesis problems

This section is devoted to H ∞ synthesis problems. The synthesis procedure is based on the scheme (3) and must be initialized with a stabilizing controller. This initial phase I is described in the previous section, but alternative techniques to find an initial stabilizing controller K may prove useful [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF]. All examples are extracted from the COM P l e ib collection [START_REF] Leibfritz | COMPLeIB, COnstraint Matrix-optimization Problem LIbrary -a collection of test examples for nonlinear semidefinite programs, control system design and related problems[END_REF]. Column 'problem' now indicates the COM P l e ib acronym attached to each example.

Notice that for a static K, our program could formally be solved as minimize

T w→z (K) ∞ , K ∈ M k+m2,k+p2 subject to C 2 (sI -A(K)) -1 B 2 ∞ ≤ M (27)
where M > 0 is a constant e.g. such that the initial stabilizing controller K 0 satisfies

C 2 (sI -A(K 0 )) -1 B 2 ∞ ≪ M < +∞
and where A(K) refers to the closed-loop system matrix. Note that the extra constraint in ( 27) maintains asymptotic stability of the closed-loop system during the optimization of K.

It often happens that stability of the performance channel T w→z (K) alone already implies stability of the closed-loop system, so that the (internal) stability constraint in ( 27) is redundant. This observation is significant, because this constraint can be dropped, and the problem becomes unconstrained. In those cases where this does not hold true, we may if we want to avoid the big-M constraint in [START_REF] Desoer | Feedback Systems: Input-Output Properties[END_REF], introduce the composite problem minimize

K T w→z (K) 0 0 ǫ s C 2 (sI -A(K)) -1 B 2 ∞
where the lower-right term enforces internal stability, with ǫ s > 0 a small enough parameter. The new problem can now be handled by the techniques discussed so far without change. Similar modifications apply to design problem with additional structural constraints on the controller. We compare the results of our nonsmooth algorithm variant II in columns 'nonsmooth H ∞ ' to older results obtained with the specialized augmented Lagrangian (AL) algorithm described in [START_REF] Apkarian | Fixed-order H∞ control design via an augmented Lagrangian method[END_REF], displayed in columns 'H ∞ AL', and to results obtained with the Frank & Wolfe (FW) algorithm described in [START_REF] Elghaoui | An algorithm for static outputfeedback and related problems[END_REF], column 'FW'. In column 'H ∞ full', we display the gain obtained with a full-order feedback controller, synthesized via the usual Riccati-based DGKF technique [START_REF] Doyle | State-Space Solutions to Standard H 2 and H∞ Control Problems[END_REF]. Note that it gives the best achievable H ∞ performance and thus provides a lower bound for other techniques.

The results which we obtain with our nonsmooth technique are usually close to those previously obtained with the augmented Lagrangian method [START_REF] Apkarian | Fixed-order H∞ control design via an augmented Lagrangian method[END_REF], except for problems with large state dimension as 'AC10' (55 states), 'BDT2' (82 states), 'HF1' (130 states) and 'CM4' (240 states), where the augmented Lagrangian method fails, while the present nonsmooth method is still functional. For these large systems, 'AC10', 'BDT2', 'HF1' and 'CM4', we have observed that even Riccati or LMI solvers encounter serious difficulties or even break down.

Notice that our present nonsmooth technique (NS) and the AL method are rigorous in the sense that they converge to local minima (critical points). It is therefore not surprising that NS and AL often achieve the same H ∞ performance at the same K, for which optimality was established in [START_REF] Noll | Partially augmented Lagrangian method for matrix inequality constraints[END_REF].

Note that in its original form, the FW method cannot solve problems where performance is optimized under constraints on the order of the controller. We have therefore encapsulated FW into a dichotomy search in order to assure the best possible performance. According to our numerical experiments, AL and FW, which solve SDPs at every iteration, are no longer functional even for medium size systems such as the Boeing problem 'AC10'. Higher-order problems like 'BDT2', 'HF1', and 'CM4' are completely intractable with these techniques, while the nonsmooth method continues to produce valid solutions. We observed that FW, when functional, is outperformed both by AL and NS. We attribute this to the fact that, as opposed to AL and NS, FW is not supported by a sound convergence theory, and therefore often stops at iterates K which are not even locally optimal. Examples of K where this is the case are easily identified. It suffices to start AL or NS at iterates K where FW stops which almost always leads to further improvement. As an example, when we initialized NS with the FW solution in example 'AC8', the H ∞ performance was improved from 2.612 to 2.005. A phenomenon that we have also observed in examples 'HE1' and 'REA2'. This is a strong argument in favor of those optimization techniques, which generate steps based on local convergence theory. It means that NS and AL can be used to certify criticality of controllers obtained with alternative methods.

As an illustration of the nonsmooth technique, Fig. 1 shows the evolution of the maximum singular value of T w→z (K) for example 'AC8' during the first 5 iterations. The stars indicate the frequencies ω i that were regrouped in the set Ω k to construct a bundle of Clarke subgradients at iterate x k . Also, Fig. 2 depicts the evolution of the absolute value of the criticality measure of algorithm variant II [START_REF] Xin | Reduced-order controllers for the H∞ control problem with unstable invariant zeros[END_REF] and [START_REF] Xin | Reduced-order controllers for continuous and discrete-time singular H∞ control problems based on LMI[END_REF] versus iterations. As theoretically expected, θ(x k ) gradually tends to zero until a local minimum is reached. Finally, controllers for large problems 'AC10', 'BDT2', 'HF1' and 'CM4' are given in Table III.

VIII. CONCLUSION

We have proposed several new algorithms to minimize the H ∞ norm subject to structural constraints on the controller dynamics. The proposed method uses nonsmooth techniques suited for H ∞ synthesis and for semi-infinite eigenvalue or singular value optimization programs. Variant I and variant II of our algorithm are supported by global convergence theory, a crucial parameter for the reliability of an algorithm in practice. Note that the proposed tools and techniques easily extend to multidisk problems and synthesis problems on prescribed frequency intervals [START_REF] Apkarian | Nonsmooth optimization for multidisk H∞ synthesis[END_REF]. More importantly, they pave the way for investigating an even larger scope of synthesis problems, characterized through frequency domain inequalities of the form λ 1 (H(x, ω)) ≤ 0, ω ≥ 0, where H(x, ω) is Hermitianvalued and x stands for controller parameters and possibly multiplier variables, as is the case when IQC formulations are used [START_REF] Megretski | System analysis via integral quadratic constraints[END_REF]. This is a strong incentive for further development and research. Also, a second-order version of our technique with enhanced asymptotic convergence is currently under investigation [START_REF] Bompart | Second-order nonsmooth optimization for H∞ and H 2 syntheses[END_REF].

  Manuscript received January 20, 2002; revised November 18, 2002. M. Apkarian is with ONERA and Universit Paul Sabatier, Maths. Dept. M. Noll is with Universit Paul Sabatier, Maths. Dept

Example 3 .

 3 Simultaneous stabilization. Another instance of interest is the simultaneous stabilization problem, which can be cast as minimizing a finite family of closed-loop transfer functions. Formally, given the open-loop plants, G i , i = 1, . . . , r, we consider the problem minimize K max i=1,...,r

Example

and goto step 5 . 4 )

 54 If α k > τ 2 π k then reject t k and add nodes to Ω k to obtain the finer mesh Ω k+1 . Increase counter k by one and go back to step 2. 5) Increase counter k by one and go back to step 1.

  Variant II has been shown to perform satisfactorily on a number of difficult examples. In particular, four examples with large state dimension (n = 55 n = 82, n = 130 and n = 240) have been solved.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Max singular values (transport airplane 'AC8') versus frequencyfirst 5 iterations -'*' selected frequencies

  Burke et al. present an alternative approach to computing static stabilizing controllers K. The authors propose to solve the nonsmooth optimization program minimize

  P 11 + P 12 K(I -P 22 K) -1 P 21 ,where the state-space data of P 11 , P 12 , P 21 and P 22 are given in[START_REF] Gahinet | A linear matrix inequality approach to H∞ control[END_REF] and the dependence on s is omitted for brevity. Our aim is to compute the subdifferential ∂f

			section 2.3].
	IV. CLARKE SUBDIFFERENTIALS IN CLOSED-LOOP
	Given a stabilizing controller K(s) and a plant with the
	usual partition
	P (s) :=	P 11 (s) P 12 (s) P 21 (s) P 22 (s)	,
	the closed-loop transfer function is obtained as
	T w→z (K) :=

TABLE II H∞

 II SYNTHESIS WITH NONSMOOTH ALGORITHMIC VARIANT II

	problem	(n, m, p)	order	iter	cpu (sec.)	nonsmooth H∞	H∞ AL	FW	H∞ full
	AC8	(9, 1, 5)	0	20	45	2.005	2.02	2.612	1.62
	HE1	(4, 2, 1)	0	4	7	0.154	0.157	0.215	0.073
	REA2	(4, 2, 2)	0	31	51	1.192	1.155	1.263	1.141
	AC10	(55, 2, 2)	0	15	294	13.11	*	*	3.23
	AC10	(55, 2, 2)	1	46	408	10.21	*	*	3.23
	BDT2	(82, 4, 4)	0	44	1501	0.8364	*	*	0.2340
	HF1	(130, 1, 2)	0	11	1112	0.447	*	*	0.447
	CM4	(240, 1, 2)	0	2	3052	0.816	*	*	*
				ǫω = 0.05, ' * ': problem is intractable

TABLE III H∞

 III CONTROLLERS FOR LARGE PROBLEMS WITH ALGORITHM II

	problem	order							K(s)
	AC10	0					K =	»	4.3042 -9.0747e-1 2.1249e-5 2.2467e-5	-
	AC10	1	»	-6.110e-2 -2.361e-2	-	(s + 2.633) -1 [ -6.476e-2 1.245e-2 ] +	»	-3.285e-01 5.226e-5 2.698e-1 1.389e-5	-
						2	-1.0402 1.2997	1.4684 3.0555	3
	BDT2	0				6 4	-3.7306 -3.5572 -1.9894 2.2705 -2.8162 2.1944 1.8146 9.6202	7 5
							-9.8493 -5.5376 -2.4885 5.1699
	HF1	0						[ -0.1907 -1.3093 ]
	CM4	0						[ -4.5684 -0.6667 ]

ACKNOWLEDGMENTS

Thanks to Friedemann Leibfritz, Universität Trier, for providing the COM P l e ib collection and for fruitful discussions.