Pierre Apkarian
email: pierre.apkarian@onera.fr

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

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

Parametric robust structured control design

Keywords: Real uncertain parameters, structured H∞synthesis, parametric robust control, nonsmooth optimization, local optimality, inner approximation

come

I. INTRODUCTION

P ARAMETRIC uncertainty is among the most challenging problems in control system design due to its NP-hardness. Albeit, being able to provide solutions to this fundamental problem is a must for any practical design tool worthy of this attribute. Not surprisingly, therefore, parametric uncertainty has remained high up on the agenda of unsolved problems in control for the past three decades.

It is of avail to distinguish between analysis and synthesis techniques for parametric robustness. Analysis refers to assessing robustness of a closed-loop system when the controller is already given. If the question whether this given controller renders the closed loop parametrically robustly stable is solved exhaustively, then it is already an NP-hard problem [START_REF] Özbay | On the N P-hardness of the purely complex µ computation, analysis/synthesis, and related problems in multi-dimensional systems[END_REF]. Parametric robust synthesis, that is, computing a controller which is robust against uncertain parameters, is even harder, because it essentially involves an iterative procedure where at every step an analysis problem is solved. Roughly, we could say that in parametric robust synthesis we have to optimize a criterion, a single evaluation of which is already NP-hard.

For the analysis of parametric robustness, theoretical and practical tools with only mild conservatism and acceptable CPUs have been proposed over the years [START_REF] Packard | Linear, multivariable robust control with a µ perspective[END_REF]. In contrast, no tools with comparable merits in terms of quality and CPU are currently available for synthesis. It is fair to say that the parametric robust synthesis problem has remained open. The best currently available techniques for synthesis are the µ tools going back to [START_REF] Balas | µ-Analysis and synthesis toolbox: User's Guide[END_REF], made available to designers through the MATLAB Robust Control Toolbox. These rely on upper bound approximations of µ and follow a heuristic which alternates between analysis and synthesis steps. When this works, it gives performance and stability certificates, but the approach may turn out conservative, and the computed controllers are often too complicated for practice.

The principal obstruction to efficient robust synthesis is the inherent nonconvexity and nonsmoothness of the mathematical program underlying the design. These obstacles have to some extent been overcome by the nonsmooth approaches [START_REF] Polak | Nondifferentiable optimization algorithm for designing control systems having singular value inequalities[END_REF], [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF], [START_REF] Burke | Stabilization via nonsmooth, nonconvex optimization[END_REF] for control. These techniques allow to address multi-model and multi-objective structured control design as discussed in [START_REF] Gahinet | Automated tuning of gain-scheduled control systems[END_REF], [START_REF] Apkarian | Tuning controllers against multiple design requirements[END_REF], [START_REF] Apkarian | Optimization-based control design techniques and tools[END_REF], [START_REF] Burke | HIFOO -A Matlab package for fixed-order controller design and H∞ optimization[END_REF], [START_REF] Apkarian | Nonsmooth optimization for multidisk H∞ synthesis[END_REF]. These have become available to designers through synthesis tools like HINFSTRUCT or SYSTUNE from [START_REF]Robust Control Toolbox 5.0[END_REF]. However, these methods can no longer be used for the substantially harder parametric robust synthesis problem, and this has been an incentive for our present work, where we investigate new classes of optimization programs involving upper-C 1 stability and performance criteria.

In order to understand our approach, it is helpful to distinguish between inner and outer approximations of the robust control problem on a set ∆ of uncertain parameters. Outer approximations relax the problem over ∆ by choosing a larger, but more convenient, set ∆ ⊃ ∆, the idea being that the problem on ∆ becomes accessible to computations. If solved successfully on ∆, this provides performance and robustness certificates for ∆. Typical tools in this class are the upper bound approximation µ of the structured singular value µ developed in [START_REF] Fan | Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics[END_REF], the DK-iteration function DKSYN of [START_REF]Robust Control Toolbox 5.0[END_REF], or LMI-based approaches like [START_REF] Balakrishnan | Linear matrix inequalities in robustness analysis with multipliers[END_REF]. The principal drawback of outer approximations is the inherent conservatism, which increases significantly with the number of uncertainties and their repetitions, and the fact that failures occur more often.

Inner approximations are preferred in practice and relax the problem by solving it on a smaller typically finite subset ∆ a ⊂ ∆. This avoids conservatism and leads to acceptable CPUs, but has the disadvantage that no immediate stability or performance certificate for ∆ is obtained. Our principal contribution here is to show a way how this shortcoming can be avoided. We present an efficient technique to compute an inner approximation with structured controllers with a local optimality certificate such that robust stability and performance are achieved over ∆ in the majority of cases. We then also show how this can be certified a posteriori over ∆, when combined with outer approximation for analysis. The new method we propose is termed dynamic inner approximation, as it generates the inner approximating set ∆ a dynamically. The idea of using inner approximations, and thus multiple models, to solve robust synthesis problems is not new and was employed in different contexts, see e.g. [START_REF] Nyström | Multimodel robust control of nonlinear plants: a case study[END_REF], [START_REF] Magni | A multimodel-based approach to robust and self-scheduled control design[END_REF], [START_REF] Ackermann | Robust control. Systems with Uncertain Physical Parameters[END_REF].

To address the parametric robust synthesis problem we use a nonsmooth optimization method tailored to minimizing a cost function, which is itself a semi-infinite minimum of smooth functions. This is in contrast with previously discussed nonsmooth optimization problems, where a semi-infinite maximum of smooth functions is minimized, and which have been dealt with successfully in [START_REF] Apkarian | Optimization-based control design techniques and tools[END_REF]. At the core of our new approach is therefore understanding the principled difference between a min-max and a min-min problem, and the algorithmic strategies required to solve them successfully. Along with the new synthesis approach, our key contributions are

• an in-depth and rigorous analysis of worst-case stability and worst-case performance problems over a compact parameter range.

• the description of a new resolution algorithm for worstcase programs along with a proof of convergence in the general nonsmooth case. Note that convergence to local minima from an arbitrary, even remote, starting point is proved. Convergence to a global minimum is as a rule unrealistic in terms of CPU due to the NP-hardness of the problems.

The paper is organized as follows. Section II states the problem formally, and Subsection II-B presents our novel dynamic inner approximation technique and the elements needed to carry it out. Section III highlights the principal differences between nonsmooth min-min and min-max problems. Sections IV-A and IV-B examine the criteria which arise in the optimization programs, the H ∞ -norm, and the spectral abscissa. Section V presents the optimization method we designed for min-min problems and Subsections V-B, V-C are dedicated to convergence analysis. Section VI-A presents an assessment and a comparison of our algorithm on a bench of test examples. Section VI-B gives a more refined study of a challenging missile control problem.

NOTATION

For complex matrices X H denotes conjugate transpose. For Hermitian matrices, X 0 means positive definite, X 0 positive semi-definite. We use concepts from nonsmooth analysis covered by [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]. For a locally Lipschitz function f : R n → R, ∂f (x) denotes its (compact and convex) Clarke subdifferential at x ∈ R n . The Clarke directional derivative at x in direction d ∈ R n can be computed as

f • (x, d) = max g∈∂f (x) g T d .
The symbols F l , F u denote lower and upper Linear Fractional Transformations (LFT) [START_REF] Zhou | Robust and Optimal Control[END_REF]. For partitioned 2 × 2 block matrices, stands for the Redheffer star product [START_REF] Redheffer | On a certain linear fractional transformation[END_REF].

II. PARAMETRIC ROBUSTNESS

A. Setup

We consider an LFT plant as in Fig. 1 with real parametric uncertainties F u (P, ∆) where

P (s) :        ẋ = Ax + B p p + B w w + Bu q = C q x + D qp p + D qw w + D qu u z = C z x + D zp p + D zw w + D zu u y = Cx + D yp p + D yw w + Du (1)
and x ∈ R nx is the state, u ∈ R m2 the control, w ∈ R m1 the vector of exogenous inputs, y ∈ R p2 the output, and z ∈ R p1 the regulated output. The uncertainty channel is defined as p = ∆q where the uncertain matrix ∆ is without loss assumed to have the block-diagonal form

∆ = diag [δ 1 I r1 , . . . , δ m I rm ] (2)
with δ 1 , . . . , δ m representing real uncertain parameters, and r i giving the number of repetitions of δ i . We assume without loss that δ = 0 represents the nominal parameter value. Moreover, we consider δ ∈ ∆ in one-to-one correspondence with the matrix ∆ in [START_REF] Packard | Linear, multivariable robust control with a µ perspective[END_REF].

We consider an LFT plant with real parametric uncertainties F u (P, ) where

(1) P (s) :

8 > > < > > : ẋ = Ax + B p p + B w w + Bu q = C q x + D qp p + D qw w + D qu u z = C z x + D zp p + D zw w + D zu u y = Cx + D yp p + D yw w + Du
and x 2 R nx is the state, u 2 R m2 the control, w 2 R m1 the vector of exogenous inpu y 2 R p2 the output, and z 2 R p1 the regulated output. The uncertainty channel defined as p = q where the uncertain matrix is without loss assumed to have t block-diagonal form

= diag [ 1 I r1 , . . . , m I rm ] (2)
with 1 , . . . , m representing real uncertain parameters, and r i giving the number of rep titions of i . We assume without loss that = 0 represents the nominal parameter val Moreover, we consider 2 in one-to-one correspondence with the matrix in (2). Given a compact convex set ⇢ R m containing = 0, the parametric robust structur H 1 control problem consists in computing a structured output-feedback controller u K( ⇤ )y with the following properties:

(i) Robust stability. K( ⇤ ) stabilizes F u (P, ) internally for every 2 .

(ii) Robust performance. Given any other robustly stabilizing controller K() w the same structure, the optimal controller satisfies max

2 kT zw ( ,  ⇤ ) k 1  max 2 kT zw ( , ) k 1 .
Here T zw ( , ) := F l (F u (P, ( )), K()) denotes the closed-loop transfer function of t performance channel w ! z of (1) when the control loop with controller K() and t uncertainty loop with uncertainty are closed.

Fig. 1: Robust synthesis interconnection

Given a compact convex set ∆ ⊂ R m containing δ = 0, the parametric robust structured H ∞ control problem consists in computing a structured output-feedback controller u = K(κ * )y with the following properties:

(i) Robust stability. K(κ * ) stabilizes F u (P, ∆) internally for every δ ∈ ∆. (ii) Robust performance. Given any other robustly stabilizing controller K(κ) with the same structure, the optimal controller satisfies

max δ∈∆ T zw (δ, κ * ) ∞ ≤ max δ∈∆ T zw (δ, κ) ∞ .
Here T zw (δ, κ) := F l (F u (P, ∆(δ)), K(κ)) is the closed-loop transfer function of the performance channel w → z of (1) when the control loop with K(κ) and the uncertainty loop with ∆ are both closed. We recall that according to [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF] a controller

K(κ) : ẋK = A K (κ)x K + B K (κ)y u = C K (κ)x K + D K (κ)y (3)
in state-space form is called structured if A K (κ), B K (κ), . . . depend smoothly on a design parameter κ varying in a design space R n or in some constrained subset of R n . Typical examples of structure include PIDs, reduced-order controllers, observer-based controllers, or complex control architectures combining controller blocks such as set-point filters, feedforward, washout or notch filters, and much else [START_REF] Apkarian | Optimization-based control design techniques and tools[END_REF]. In contrast, full-order controllers are state-space representations with the same order as P (s) without particular structure, and are sometimes referred to as unstructured, or as black-box controllers.

Parametric robust control is among the most challenging problems in linear feedback control. The structured singular value µ developed in [START_REF] Zhou | Robust and Optimal Control[END_REF] is the principled theoretical tool to describe problem (i), (ii) formally. In the same vein, based on the spectral abscissa α(A) = max{Re(λ) : λ eigenvalue of A} of a square matrix A, criterion (i) may be written as

max δ∈∆ α (A(δ, κ * )) < 0, (4)
where A(δ, κ) is the A-matrix of the closed-loop transfer function T zw (δ, κ).

If the uncertain parameter set is a cube ∆ = [-1, 1] m , which is general enough for applications, then the same information is obtained from the distance to instability in the maximum-norm

d * = min{ δ ∞ : α (A(δ, κ * )) ≥ 0}, (5)
because criterion (i) is now equivalent to d * > 1. It is known that the computation of any of these elements, µ, (4), or ( 5) is NP-complete, so that their practical use is limited to analysis of small problems, or to the synthesis of tiny ones. Practical approaches have to rely on intelligent relaxations, or heuristics, which use either inner or outer approximations.

In the next chapters we will develop our dynamic inner approximation method to address problem (i), (ii). We solve the problem on a relatively small set ∆ a ⊂ ∆, which we construct iteratively.

B. Dynamic inner approximation

The following static inner approximation to (i), (ii) is near at hand. After fixing a sufficiently fine approximating static grid ∆ s ⊂ ∆, one solves the multi-model H ∞ -problem

min κ∈R n max δ∈∆s T zw (δ, κ) ∞ . (6)
This may be addressed with the methods of [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF], [START_REF] Apkarian | Nonsmooth optimization for multidisk H∞ synthesis[END_REF], [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF], and the software tools HINFSTRUCT and SYSTUNE [START_REF] Apkarian | Multi-model, multi-objective tuning of fixed-structure controllers[END_REF], [START_REF] Gahinet | Automated tuning of gain-scheduled control systems[END_REF], [START_REF] Apkarian | Tuning controllers against multiple design requirements[END_REF], [START_REF]Robust Control Toolbox 5.0[END_REF], but has a high computational burden due to the large number of scenarios in ∆ s , which makes it prone to failure. Straightforward gridding becomes very quickly intractable for sizable dim(δ).

Here we advocate a different strategy, which we call dynamic inner approximation, because it operates on a substantially smaller set ∆ a ⊂ ∆ generated dynamically, whose elements are called the active scenarios, which we update a couple of times by applying a search procedure locating problematic parameter scenarios in ∆. This leads to a rapidly converging procedure, much less prone to failure than [START_REF] Burke | Stabilization via nonsmooth, nonconvex optimization[END_REF]. The method can be summarized as shown in Algorithm 1.

The principal elements of Algorithm 1 will be analyzed in the following sections. We will focus on the optimization programs v * in step 4, α * in step 3, and d * , h * in step 6, which represent a relatively unexplored type of nonsmooth programs, with some common features which we shall put into evidence here. In contrast, program v * in step 2 is accessible to numerical methods through [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF] and can be addressed with Algorithm 1. Dynamic inner approximation for parametric robust synthesis over ∆ Parameters: ε > 0.

Step 1 (Nominal synthesis). Initialize the set of active scenarios as ∆ a = {0}.

Step 2 (Multi-model synthesis). Given the current finite set ∆ a ⊂ ∆ of active scenarios, compute a structured multi-model H ∞ -controller by solving

v * = min κ∈R n max δ∈∆a T zw (δ, κ) ∞ .
The solution is the structured H ∞ -controller K(κ * ).

Step 3 (Destabilization). Try to destabilize the closed-loop system T zw (δ, κ * ) by solving the destabilization problem

α * = max δ∈∆ α (A(δ, κ * )) .
If α * ≥ 0, then the solution δ * ∈ ∆ destabilizes the loop. Include δ * in the active scenarios ∆ a and go back to step 2. If no destabilizing δ was found then go to step 4.

Step 4 (Degrade performance). Try to degrade the robust H ∞ -performance by solving

v * = max δ∈∆ T zw (δ, κ * ) ∞ .
The solution is δ * .

Step 5 (Stopping test). If v * < (1 + ε)v * degradation of performance is only marginal. Then exit, or optionally, go to step 6 for post-processing. Otherwise include δ * among the active scenarios ∆ a and go back to step 2.

Step 6 (Post-processing). Check robust stability (i) and performance (ii) of K(κ * ) over ∆ by computing the distance d * to instability [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF], and its analogue

h * = min{ δ ∞ : T zw (δ, κ * ) ∞ ≥ v * }.
Possibly use µtools from [START_REF]Robust Control Toolbox 5.0[END_REF] to assess d * , h * approximately. If all δ * obtained satisfy δ * ∈ ∆, then terminate successfully.

tools like HINFSTRUCT or SYSTUNE available through [START_REF]Robust Control Toolbox 5.0[END_REF], or HIFOO available through [START_REF] Burke | HIFOO -A Matlab package for fixed-order controller design and H∞ optimization[END_REF]. Note that our approach is heuristic in so far as we have relaxed (i) and (ii) by computing locally optimal solutions, so that a global stability/performance certificate is only provided in the end as a result of step 6.

III. NONSMOOTH MIN-MAX VERSUS MIN-MIN PROGRAMS

A. Classification of the programs in Algorithm 1

Introducing the functions a ± (δ) = ±α (A(δ)), the problem of step 3 can be equivalently written in the form

minimize a -(δ) = -α (A(δ)) subject to δ ∈ ∆ (7)
for a matrix A(δ) depending smoothly on the parameter δ ∈ R m . Here the dependence of the matrix on controller K(κ * ) is omitted for simplicity, as the latter is fixed in step 3 of the algorithm. Similarly, if we introduce h

± (δ) = ± G(δ) ∞ ,
with G(s, δ) a transfer function depending smoothly on δ ∈ R m , then problem of step 4 has the abstract form

minimize h -(δ) = -G(δ) ∞ subject to δ ∈ ∆ (8)
where again controller K(κ * ) is fixed in step 4, and therefore suppressed in the notation. In contrast, the

H ∞ -program v * in step 2 of Algorithm 1 has the form minimize h + (κ) = G(κ) ∞ subject to κ ∈ R n (9)
which is of the more familiar min-max type. Here we use the well-known fact that the H ∞ -norm may be written as a semi-infinite maximum function h

+ (κ) = max ω∈[0,∞] σ (G(κ, jω)).
The maximum over the finitely many δ ∈ ∆ a in step 2 complies with this structure and may in principle be condensed into the form ( 9), featuring only a single transfer G(s, κ). In practice this is treated as in [START_REF] Apkarian | Nonsmooth optimization for multidisk H∞ synthesis[END_REF].

Due to the minus sign, programs ( 7) and ( 8), written in the minimization form, are now of the novel min-min type, which is given special attention here. This difference is made precise by the following Definition 1 (Spingarn [START_REF] Spingarn | Submonotone subdifferentials of Lipschitz functions[END_REF]). A locally Lipschitz function f :

R n → R is lower-C 1 at x 0 ∈ R n if there exist a compact space K, a neighborhood U of x 0 , and a mapping F : R n × K → R such that f (x) = max y∈K F (x, y) (10)
for all x ∈ U , and F and ∂F/∂x are jointly continuous. The function f is said to be upper-

C 1 if -f is lower-C 1 .
We expect upper-and lower-C 1 functions to behave quite differently in descent algorithms. Minimization of lower-C 1 functions, as required in [START_REF] Apkarian | Optimization-based control design techniques and tools[END_REF], should lead to a genuinely nonsmooth problem, because iterates of a descent method move toward the points of nonsmoothness. In contrast, minimization of upper-C 1 functions as required in [START_REF] Gahinet | Automated tuning of gain-scheduled control systems[END_REF] and ( 8) is expected to be better behaved, because iterates move away from the nonsmoothness. Accordingly, we will want to minimize upper-C 1 functions in [START_REF] Gahinet | Automated tuning of gain-scheduled control systems[END_REF] and [START_REF] Apkarian | Tuning controllers against multiple design requirements[END_REF] in much the same way as we optimize smooth functions in classical nonlinear programming, whereas the minimization of lower-C 1 functions in (9) requires specific techniques like nonconvex bundle methods [START_REF] Apkarian | A proximity control algorithm to minimize nonsmooth and nonconvex semi-infinite maximum eigenvalue functions[END_REF], [START_REF]A trust region spectral bundle method for nonconvex eigenvalue optimization[END_REF], [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF]. See Fig. 2 for an illustration.

Remark 1 (Distance to instability). Note that the computation of the distance to instability d * defined in (5) for step 6 of Algorithm 1 has also the features of a min-min optimization program. Namely, when written in the form

minimize t subject to -t ≤ δ i ≤ t, i = 1, . . . , m -α (A(δ)) ≤ 0 (11)
with variable (δ, t) ∈ R m+1 , the Lagrangian of ( 5) is

L(δ, t, λ, µ ± ) = t+ m i=1 µ i-(-t -δ i )+µ i+ (δ i -t)-λα (A(δ))
for Lagrange multipliers λ ≥ 0 and µ ± ≥ 0. In particular, if (δ * , t * , λ * , µ * ± ) is a Karush-Kuhn-Tucker point of ( 11), (see [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]Prop. 6.6.4], or [START_REF] Bertsekas | Nonlinear Programming[END_REF]Sect. 3.3.1]), then the local minimum (δ * , t * ) we are looking for is also a critical point of the unconstrained program

min δ∈R m ,t∈R L(δ, t, λ * , µ * ± ),
which features the function a -and is therefore of min-min type. Therefore, in solving (5), we expect phenomena of minmin type to surface rather than those of a min-max program.

A similar comment applies to the computation of h * in step 6 of the algorithm.

Remark 2 (Well-posedness). Yet another aspect of Algorithm 1 is that in order to be robustly stable over the parameter set ∆, the LFTs must be well-posed in the sense that (I -∆D) -1 exists for every δ ∈ ∆, where D is the closed-loop D-matrix.

Questioning well-posedness could therefore be included in step 3 of the algorithm, or added as posterior testing in step 6. It can be formulated as yet another min-min program

minimize -σ((I -∆D) -1 ) subject to δ ∈ ∆ ( 12
)
where one would diagnose the solution δ * to represent an illposed scenario as soon as it achieves a large negative value. Program ( 12) exhibits the same properties as minimizing h -in Section IV-A and is handled with the same novel techniques.

For programs v * in step 4, α * in step 3, and d * , h * in step 6 of Algorithm 1, well-posedness ( 12) is a prerequisite. However, we have observed that it may not be necessary to question well-posedness over ∆ at every step, since questioning stability over ∆ has a similar effect. Since the posterior certificate in step 6 of the algorithm covers also well-posedness, this is theoretically justified. Remark 3. Our notation makes it easy for the reader to distinguish between min-min and min-max programs. Namely, minimizations over the controller variable κ turn out the minmax ones, while minimizations over the uncertain parameters δ lead to the min-min type.

Remark 4. Note that min-min programs [START_REF] Gahinet | Automated tuning of gain-scheduled control systems[END_REF], [START_REF] Apkarian | Tuning controllers against multiple design requirements[END_REF] require only smoothness of A(δ), respectively, G(δ). This is practical in problems where LFT representations are not readily available.

B. Highlighting the difference between min-max and min-min

In this section we look at the typical difficulties which surface in min-max and min-min programs. This is crucial for the understanding of our algorithmic approach. Consider first a min-max program of the form

min κ∈R n max i∈I f i (κ), (13)
where the f i are smooth. When the set I is finite, we may simply dissolve this into a classical nonlinear programming (NLP) using one additional dummy variable t ∈ R:

minimize t subject to f i (κ) ≤ t, i ∈ I.
The situation becomes more complicated as soon as the set I is infinite, as is for instance the case in program v * in step 2 of Algorithm 1. The typical difficulty in min-max programs is to deal with this semi-infinite character, and one is beholden to use a tailored solution, as for instance developed in [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF], [START_REF]A trust region spectral bundle method for nonconvex eigenvalue optimization[END_REF], [START_REF] Apkarian | A proximity control algorithm to minimize nonsmooth and nonconvex semi-infinite maximum eigenvalue functions[END_REF]. Altogether this type of difficulty is well-known and has been thoroughly studied.

In contrast, a min-min program

min δ∈R n min i∈I f i (δ) (14)
cannot be converted into an NLP even when I is finite. The problem has disjunctive character, and if solved to global optimality, min-min programs lead to combinatorial explosion. On the other hand, a min-min problem has some favorable features when it comes to solely finding a good local minimum. Namely, when meeting a nonsmooth iterate δ j , where several branches f i are active, we can simply pick one of those branches and continue optimization as if the objective function were smooth. In the subsequent sections we prove that this intuitive understanding is indeed correct. Our experimental section will show that good results are obtained if a good heuristic is used.

The above considerations lead us to introduce the notion of active indices and branches for functions f (δ) defined by the inner max and min in ( 13) and ( 14). Definition 2. The set of active indices for f at δ is defined as

I(δ) := {i ∈ I : f i (δ) = f (δ)} .
Active branches of f at δ are those corresponding to active indices, i.e, f i , i ∈ I(δ). along which branch to proceed, hence the potentially combinatorial character of [START_REF] Balakrishnan | Linear matrix inequalities in robustness analysis with multipliers[END_REF]. Minimizing a maximum (left) requires simultaneous minimization of active branches, hence the non smoothness and the potentially semi-infinite character of (13).

IV. COMPUTING SUBGRADIENTS

In this section we briefly discuss how the subgradient information needed to minimize h -and a -is computed.

A. Case of the H ∞ -norm

We start by investigating the case of the H ∞ -norm h ± . We recall that function evaluation is based on the Hamiltonian algorithm of [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] and its further developments [START_REF] Benner | L∞-norm computation for continuous-time descriptor systems using structured matrix pencils[END_REF].

Computation of subgradients of h -in the sense of Clarke can be adapted from [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF], see also [START_REF] Boyd | Linear Controller Design: Limits of Performance[END_REF]. We assume the controller is fixed in this section and investigate the properties of h - as a function of δ. To this aim, the controller loop is closed by substituting the structured controller (3) in ( 1), and we obtain the transfer function M (κ) := F l (P, K(κ)). Substantial simplification in Clarke subdifferential computation is then obtained by defining the 2 × 2-block transfer function * T qw (δ)

T zp (δ) T zw (δ) := 0 I I ∆ M , (15)
where the dependence on κ has now been suppressed, as the controller will be fixed to κ * after step 2. It is readily seen that T zw coincides with the closed-loop transfer function where both controller and uncertainty loops are closed. Now consider the function h -(δ) := -T zw (δ) ∞ , which is well defined on its domain D := {δ ∈ R m : T zw (δ) is internally stable}. We have the following Proposition 1. The function h -is everywhere Clarke subdifferentiable on D. The Clarke subdifferential at δ ∈ D is the compact and convex set

∂h -(δ) = φ Y : Y = (Y ω ), ω ∈ Ω(δ), Y ω 0, ω∈Ω(δ) Trace(Y ω ) = 1 ,
where the i-th entry of φ Y is Trace ∆ T i Φ Y with ∆ i = ∂∆/∂δ i , and

Φ Y = - ω∈Ω(δ) Re T qw (δ, jω)P ω Y ω Q H ω T zp (δ, jω) T .
Here Ω(δ) is the set of active frequencies at δ, Q ω is a matrix whose columns are the left singular vectors associated with the maximum singular value of T zw (δ, jω), P ω is the corresponding matrix of right singular vectors, and Y ω is an Hermitian matrix of appropriate size.

Proof: Computation of the Clarke subdifferential of h - can be obtained from the general rule ∂(-h) = -∂h, and knowledge of ∂h + , see [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF]. Note that in that reference the Clarke subdifferential is with respect to the controller and relies therefore on the Redheffer star product

P K(κ) I I 0 .
Here we apply this in the upper loop in ∆, so we have to use the analogue expression (15) instead.

Remark 5. In the case where a single frequency ω 0 is active at δ and the maximum singular value σ of T zw (δ, jω 0 ) has multiplicity 1, h -is differentiable at δ and the gradient is

∂h -(δ) ∂δ i = -Trace Re T qw (δ, jω 0 )p ω0 q H ω0 T zp (δ, jω 0 ) T ∆ i ,
where p ω0 and q ω0 are the unique right and left singular vectors of T zw (δ, jω 0 ) associated with σ(T zw (δ, jω 0 )) = h + (δ).

Proposition 2. Let D = {δ : T zw (δ) is internally stable}. Then h + : δ → T zw (δ) ∞ is lower-C 1 on D, so that h -: δ → -T zw (δ) ∞ is upper-C 1 there.
Proof: Recall that the maximum singular value has the variational representation

σ(G) = sup u =1 sup v =1 u T Gv .
Now observe that z → |z|, being convex, is lower-C 1 as a mapping R 2 → R, so we may write it as

|z| = sup l∈L Ψ(z, l)
for Ψ jointly of class C 1 and L compact. Then

h + (δ) = sup jω∈S 1 sup u =1 sup v =1 sup l∈L Ψ u T T zw (δ, jω)v, l , (16)
where S 1 = {jω : ω ∈ R ∪ {∞}} is homeomorphic with the 1-sphere. This is the desired representation [START_REF] Burke | HIFOO -A Matlab package for fixed-order controller design and H∞ optimization[END_REF], where the compact space K is obtained as K := S 1 ×{u : u = 1}×{v : v = 1} × L, F as F (δ, jω, u, v, l) := Ψ u T T zw (δ, jω)v, l and y as y := (jω, u, v, l).

B. Case of the spectral abscissa

For the spectral abscissa the situation is more complicated, as a ± is not locally Lipschitz everywhere. Recall that an eigenvalue

λ i of A(δ) is called active at δ if Re(λ i ) = α (A(δ)).
We use I(δ) for the indices of active eigenvalues. Let us write the LFT describing A(δ) as A(δ) = A + C∆(I -D∆) -1 B, where dependence on controller parameters κ is again omitted and considered absorbed into the state-space data A, B, etc. Proposition 3. Suppose all active eigenvalues λ i , i ∈ I(δ) of A(δ) at δ are semi-simple. Then a ± (δ) = ±α (A(δ)) is Clarke subdifferentiable in a neighborhood of δ. The Clarke subdifferential of a -at δ is ∂a

-(δ) = {φ Y : Y = (Y i ) i∈I(δ) , Y i 0, i∈I(δ) Trace(Y i ) = 1}, where the i-th entry of φ Y is -Trace ∆ i T Φ Y with ∆ i = ∂∆/∂δ i , and
Φ Y = i∈I(δ) Re (I -D∆) -1 CV i Y i U H i B(I -∆D) -1 T .
Here V i is a column matrix of right eigenvectors, U H i a row matrix of left eigenvectors of A(δ) associated with the eigenvalue λ i , and such that U H i V i = I. Proof: This follows from [START_REF] Burke | Differential properties of the spectral abscissa and the spectral radius for analytic matrix-valued mappings[END_REF]. See also [START_REF] Bompart | Non-smooth techniques for stabilizing linear systems[END_REF]. A very concise proof that semi-simple eigenvalue functions are locally Lipschitz could also be found in [START_REF] Lui | Pseudospectral mapping theorem II[END_REF].

When every active eigenvalue is simple, Y i reduces to a scalar y i and a fast implementation is possible. We use the LU-decomposition to solve for u i and v i in the linear systems

u H i (I -∆D) := u H i B, (I -D∆) v i := Cv i .
Given the particular structure (2) of ∆, subgradients with respect to the kth entry are readily obtained as a sum over i ∈ I(δ) of inner products of the form

y i Re u i (J(k)) H v i (J(k)),
where J(k) is a selection of indices associated with the rows/columns of δ k in ∆(δ). Similar inner products arise in the computation of H ∞ norm subgradients. It was observed in [START_REF] Burke | Differential properties of the spectral abscissa and the spectral radius for analytic matrix-valued mappings[END_REF] that a ± may fail to be locally Lipschitz at δ if A(δ) has a derogatory active eigenvalue. Proposition 4. Suppose every active eigenvalue of A(δ) is simple. Then a -is upper-C 1 in a neighborhood of δ.

Proof: If active eigenvalues are simple, then a + is the maximum of C 1 functions in a neighborhood of δ. The result follows from a -= -a + .

V. ALGORITHM FOR MIN-MIN PROGRAMS

In this section we present our descent algorithm to solve programs [START_REF] Gahinet | Automated tuning of gain-scheduled control systems[END_REF] and [START_REF] Apkarian | Tuning controllers against multiple design requirements[END_REF]. We consider an abstract form of the min-min program with f a general objective function of this type:

minimize

f (δ) subject to δ ∈ ∆ (17)
As we already pointed out, the crucial point is that we want to stay as close as possible to a standard algorithm for smooth optimization, while assuring convergence under the specific form of upper nonsmoothness in these programs.

Algorithm 2. Descent method for min-min programs

Parameters: 0 < γ < Γ < 1, 0 < θ < Θ < 1.
Step 1 (Initialize). Put outer loop counter j = 1, choose initial guess δ 1 ∈ ∆, fix memory stepsize t 1 > 0.

Step 2 (Stopping). If δ j is a Karush-Kuhn-Tucker point of (17) then exit, otherwise go to inner loop.

Step 3 (Inner loop). At current iterate δ j call the step finding Subroutine (Subroutine 1) started with last memorized stepsize t j to find a step t k > 0 and a new serious iterate δ j+1 such that

ρ k = f (δ j ) -f (δ j+1 ) f (δ j ) -φ k (δ j+1 , δ j ) ≥ γ.
Step 4 (Stepsize update). If ρ k ≥ Γ then update memory stepsize as t j+1 = θ -1 t k , otherwise update memory stepsize as t j+1 = t k . Increase counter j and go back to step 2.

In order to understand Algorithm 2 and its step finding subroutine (Subroutine 1), we recall from [START_REF] Noll | Convergence of non-smooth descent methods using the Kurdyka-Łojasiewicz inequality[END_REF], [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF] that

φ (η, δ) = f (δ) + f • (δ, η -δ)
is the standard model of f at δ, where f • (δ, d) is the Clarke directional derivative of f at δ in direction d [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]. This model can be thought of as a substitute for a first-order Taylor expansion at δ and can also be represented as

φ (η, δ) = f (δ) + max g∈∂f (δ) g T (η -δ), (18)
where ∂f (δ) is the Clarke subdifferential of f at δ. In the subroutine we generate lower approximations φ k of φ using finite subsets G k ⊂ ∂f (δ), putting

φ k (η, δ) = f (δ) + max g∈G k g T (η -δ).
We call φ k the working model at inner loop counter k.

Subroutine 1. Descent step finding for min-min programs

Input: Current serious iterate δ, last memorized stepsize t > 0. Flag. Output: Next serious iterate δ + .

Step 1 (Initialize). Put linesearch counter k = 1, and initialize search at t 1 = t . Choose subgradient g 0 ∈ ∂f (δ). Put G 1 = {g 0 }.

Step 2 (Tangent program). Given t k > 0, a finite set of Clarke subgradients G k ⊂ ∂f (δ), and the corresponding working model

φ k (•, δ) = f (δ)+max g∈G k g T (•-δ), compute solution η k ∈ ∆ of the convex quadratic tangent program (TP) min η∈∆ φ k (η, δ)+ 1 2t k η-δ 2 .
Step 3 (Armijo test). Compute

ρ k = f (δ) -f (η k ) f (δ) -φ k (η k , δ)
If ρ k ≥ γ then return δ + = η k successfully to Algorithm 2. Otherwise go to step 4

Step 4 (If Flag = strict. Cutting and aggregate plane). Pick a subgradient

g k ∈ ∂f (δ) such that f (δ) + g T k (η k -δ) = φ (η k , δ), or equivalently, f • (δ, η k -δ) = g T k (η k -δ).
Include g k into the new set G k+1 for the next sweep. Add the aggregate subgradient g * k into the set G k+1 to limit its size.

Step 5 (Step management). Compute the test quotient

ρ k = f (δ) -φ k+1 (η k , δ) f (δ) -φ k (η k , δ) . If ρ k ≥ γ then select t k+1 ∈ [θt k , Θt k ],

A. Practical aspects of Algorithm 2

The subroutine of the descent Algorithm 2 looks complicated due to step 4, but as we now argue, it reduces to a standard backtracking linesearch in the majority of cases. Note that if f is certified upper-C 1 , then we completely dispense with step 4 and keep G k = {g 0 }, which by force reduces the subroutine to a linesearch along a projected gradient direction. This is what we indicate by flag = upper in step 4 of the subroutine.

If f is known to have a strict standard model φ in [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF], without being certified upper-C 1 , which corresponds to flag = strict, then step 4 of the subroutine is in principle needed. However, even then we expect the subroutine to reduce to a standard linesearch. This is clearly the case when the Clarke subdifferential ∂f (δ) at the current iterate δ is singleton,

because φ k (η, δ) = f (δ)+∇f (x) T (η -δ) is then independent of k, so ρ k ≥ γ reads f (η k ) ≤ f (δ j ) + γ∇f (δ j ) T (η k -δ j ),
which is the usual Armijo test [START_REF] Bertsekas | Nonlinear Programming[END_REF]. Moreover, η k is then a step along the projected gradient P ∆-δ (-∇f (δ)), which is easy to compute due to the simple structure of ∆. More precisely, for ∆ = [-1, 1] m and stepsize t k > 0, the solution η of tangent program (TP) in step 2 can be computed coordinatewise as

min γ i η + (2t k ) -1 η 2 + γ i -δ i t -1 k η : -1 ≤ η ≤ 1 ,
where γ i := ∂f (δ)/∂δ i . Cutting plane and aggregate plane in step 4 become redundant, and the quotient ρ k in step 5 is also redundant as it is always equal to 1.

Remark 7. There is only one case in which step 4 is fully executed, and that is when f is not certified upper-C 1 , and in addition the subgradient g 0 ∈ ∂f (δ) in step 1 of Subroutine 1 does not satisfy f (δ) + g T 0 (η kδ) = φ (η k , δ). In that rare event step 4 requires computation of a new subgradient

g k ∈ ∂f (δ) which does satisfy f (δ) + g T k (η k -δ) = φ (δ, η k -δ).
From here on the procedure changes. The sets G k+1 may now grow, because we add g k into G k+1 . This corresponds to what happens in a bundle method. The tangent program (TP) is now solved numerically using a QP-solver. Fortunately we may limit the number of elements of G k+1 using the aggregate subgradient of Kiwiel [START_REF] Kiwiel | An aggregate subgradient method for nonsmooth convex minimization[END_REF], so this is still very fast.

Remark 8. For the spectral abscissa f (δ) = a -(δ), which is not certified upper-C 1 , we use this cautious variant, where the computation of g k in step 4 may be required. For f = a -this leads to a low-dimensional semidefinite program.

Remark 9. The stopping test in step 2 of Algorithm 2 can be delegated to Subroutine 1. Namely, if δ j is a Karush-Kuhn-Tucker point of ( 17), then η k = δ j is solution of the tangent program (TP). This means we can use the following practical stopping tests: If the inner loop at iterate δ j finds δ j+1 ∈ ∆ such that

δ j+1 -δ j 1 + δ j < tol 1 , |f (δ j+1 ) -f (δ j )| 1 + |f (δ j )| < tol 2 ,
then we decide that δ j+1 is optimal and stop. That is, the (j + 1)st inner loop is not started. On the other hand, if the inner loop at δ j has difficulties finding a new iterate and provides five consecutive unsuccessful backtracks η k such that

η k -δ j 1 + δ j < tol 1 , |f (η k ) -f (δ j )| 1 + |f (δ j )| < tol 2 ,
or if a maximum k max of linesearch steps k is exceeded, then we decide that δ j was already optimal and stop. In our experiments we use tol 1 = 10 -4 , tol 2 = 10 -4 , k max = 50.

B. Convergence analysis for the negative H ∞ -norm

Algorithm 2 was studied in much detail in [START_REF] Noll | Convergence of non-smooth descent methods using the Kurdyka-Łojasiewicz inequality[END_REF], and we review the convergence result here, applying them directly to the functions a -and h -. The significance of the class of upper-C 1 functions for convergence lies in the following Proposition 5. Suppose f is upper-C 1 at δ. Then its standard model φ is strict at δ in the following sense: For every ε > 0 there exists r > 0 such that the one-sided Taylor type estimate

f (η) ≤ φ (η, δ) + ε η -δ (19)
is satisfied for all δ, η ∈ B( δ, r).

Proof:

The following, even stronger property of upper-C 1 functions was proved in [START_REF] Dao | Bundle method for nonconvex nonsmooth constrained optimization[END_REF], see also [START_REF] Noll | Cutting plane oracles to minimize non-smooth non-convex functions[END_REF], [START_REF] Noll | Convergence of non-smooth descent methods using the Kurdyka-Łojasiewicz inequality[END_REF]. Suppose δ k → δ and η k → δ, and let g k ∈ ∂f (δ k ) be arbitrary. Then there exist ε k → 0 such that

f (η k ) ≤ f (δ k ) + g T k (η k -δ k ) + ε k η k -δ k (20)
is satisfied.

Theorem 1 (Worst-case H ∞ norm on ∆). Let δ j ∈ ∆ be the sequence generated by Algorithm 2 with standard linesearch for minimizing program [START_REF] Apkarian | Tuning controllers against multiple design requirements[END_REF]. Then the sequence δ j converges to a Karush-Kuhn-Tucker point δ * of (8).

Proof: The method of proof of [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF]Theorem 6.6] shows that every accumulation point of the δ j ∈ ∆ is a critical point of [START_REF] Apkarian | Tuning controllers against multiple design requirements[END_REF], because of the one-sided Taylor estimate [START_REF] Zhou | Robust and Optimal Control[END_REF]. However, the proof in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF] still needs the extended version of Subroutine 1, which we want to avoid. Now, using the stronger hypothesis that f is upper-C 1 , we find indeed that the full construction of the step finding subroutine is not needed, because regardless how the cutting planes are chosen, the working models always satisfy a onesided Taylor estimate

f (η k ) ≤ φ k (η k , δ k ) + ε k η k -δ k with ε k → 0 as η k , δ k → δ.
This freedom allows us to simplify the step finding Subroutine 1 to a standard linesearch, while maintaining the convergence argument of [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF]. More details are given in [START_REF] Dao | Bundle method for nonconvex nonsmooth constrained optimization[END_REF].

Since h -is upper-C 1 on D by Proposition 2, the above applies to h -. Note that hypotheses assuring boundedness of the sequence δ j used in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF], [START_REF] Noll | Convergence of non-smooth descent methods using the Kurdyka-Łojasiewicz inequality[END_REF], [START_REF] Dao | Bundle method for nonconvex nonsmooth constrained optimization[END_REF] are not needed here, since ∆ is bounded.

Convergence to a single critical point is now assured through [START_REF] Noll | Convergence of non-smooth descent methods using the Kurdyka-Łojasiewicz inequality[END_REF]Cor. 1], because G in (8) depends analytically on δ, so h - is a subanalytic function, and hence satisfies the Łojasiewicz inequality [START_REF] Bolte | The Łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems[END_REF]. Subanalyticity of h -can be derived from the following fact [START_REF] Bierstone | Semianalytic and subanalytic sets[END_REF]: If F : R n × K → R is subanalytic, and K is subanalytic and compact, then f (δ) = min y∈K F (δ, y) is subanalytic. We apply this to the negative of ( 16).

Remark 10. The lightning function f : R → R in [START_REF] Klatte | Nonsmooth Equations in Optimization[END_REF] is an example which has a strict standard model but is not upper-C 1 . It is Lipschitz with constant 1 and has ∂f (x) = [-1, 1] for every x. The standard model of f is strict, because for all x, y there exists ρ

= ρ(x, y) ∈ [-1, 1] such that f (y) = f (x) + ρ|y -x| ≤ f (x) + sign(y -x)(y -x) ≤ f (x) + f • (x, y -x) = φ (x, y -x),
using the fact that sign(yx) ∈ ∂f (x). At the same time f is certainly not upper-C 1 , because it is not semi-smooth in the sense of [START_REF] Mifflin | Semismooth and semiconvex functions in constrained optimization[END_REF]. This shows that the class of functions f with a strict standard model offers a scope of its own, justifying the effort made in the step finding subroutine.

C. Convergence analysis for the negative spectral abscissa While we obtained an ironclad convergence certificate for the H ∞ -programs (8), and similarly, for [START_REF]Robust Control Toolbox 5.0[END_REF], theory is more complicated with program [START_REF] Gahinet | Automated tuning of gain-scheduled control systems[END_REF]. In our numerical testing a -(δ) = -α (A(δ)) behaves consistently like an upper-C 1 function, and we expect this to be true at least if all active eigenvalues of A(δ * ) are semi-simple. We now show that a - has a strict standard model as a rule.

Since A(δ) depends analytically on δ, the eigenvalues are roots of a characteristic polynomial p δ (λ) = λ m + a 1 (δ)λ m-1 + • • • + a m (δ) with coefficients a i (δ) depending analytically on δ. For fixed d ∈ R m , every eigenvalue λ ν (t) of A(δ * + td) has therefore a Newton-Puiseux expansion of the form

λ ν (t) = λ ν (0) + ∞ i=k λ ν,i-k+1 t i/p (21)
for certain k, p ∈ N, where the coefficients λ ν,i = λ ν,i (d) and leading exponent k/p can be determined by the Newton polygon [START_REF] Moro | On the Lidskii-Vishik-Lyusternik perturbation theory for eigenvalues of matrices with arbitrary Jordan structure[END_REF]. If all active eigenvalues of a -(δ) = -α(A(δ)) are semi-simple, then a -is Lipschitz around δ * by Proposition 3, so that necessarily k/p ≥ 1 in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF]. It then follows that either a -(δ

* , d) = 0 when k/p > 1 for all active ν, or a -(δ * , d) = -Re λ ν,1 ≤ a • -(δ * , d) for the active ν ∈ I(δ * ) if k/p = 1.
In both cases a -satisfies the strictness estimate [START_REF] Zhou | Robust and Optimal Control[END_REF] directionally, and we expect a -to have a strict standard model. Indeed, for k/p = 1 we have a

-(δ * + td) ≤ a -(δ * ) + a • -(δ * , d)t -Re λ ν,2 t (p+1)/p + o(t (p+1)/p ), while the case k/p > 1 gives a -(δ * , d) = 0, hence a • -(δ * , d) ≥ 0, and so a -(δ * + td) ≤ a -(δ * ) -Re λ ν,1 t k/p + o(t k/p ) ≤ a -(δ * ) + a • -(δ * , d)t -Re λ ν,1 t k/p + o(t k/p ).
As soon as these estimates hold uniformly over d ≤ 1, a -has indeed a strict standard model, i.e., we have the following Lemma 1. Suppose every active eigenvalue of A(δ * ) is semisimple, and suppose the following two conditions are satisfied:

lim t→0 sup d ≤1 sup ν∈I(δ * ),k/p=1 ∞ i=k+1 Re λ ν,i-k+1 (d)t i/p-1 ≥ 0 lim t→0 sup d ≤1 sup ν∈I(δ * ),k/p>1 ∞ i=k Re λ ν,i-k+1 (d)t i/p-1 ≥ 0. (22)
Then the standard model of a -is strict at δ * , i.e., a one-sided Taylor estimate of the form (19) is satisfied.

Even though these conditions are not easy to check, they seem to be verified most of the time, so that the following result reflects what we observe in practice for the min-min program of the negative spectral abscissa a -.

Theorem 2 (Worst-case spectral abscissa on ∆). Let δ j ∈ ∆ be the sequence generated by Algorithm 2 for program [START_REF] Gahinet | Automated tuning of gain-scheduled control systems[END_REF], where the step finding subroutine is carried out with step 4 activated. Suppose every accumulation point δ * of the sequence δ j is simple or semi-simple and satisfies condition [START_REF] Apkarian | Multi-model, multi-objective tuning of fixed-structure controllers[END_REF]. Then the sequence converges to a unique Karush-Kuhn-Tucker point of program [START_REF] Gahinet | Automated tuning of gain-scheduled control systems[END_REF].

Proof: We apply [34, Cor. 1], using that a -satisfies the Łojasiewicz inequality at all accumulation points. Remark 11. The results of this section give a fair theoretical explanation why the step finding subroutine reduces to a standard linesearch in almost all cases. Indeed, as ∂a -(δ) is singleton almost everywhere in the neighborhood of a semisimple eigenvalue due to Proposition 3 and Rademacher's theorem, step 4 becomes redundant almost everywhere.

D. Multiple performance measures

Practical applications often feature several design requirements combining H ∞ and H 2 performances with spectral constraints related to pole locations. The results in Section V-B easily extend to this case upon defining H(κ, δ) := max i∈I h i (T zi,wi (κ, δ)), where several performance channels w i → z i are assessed against various requirements h i , as in [START_REF] Apkarian | Multi-model, multi-objective tuning of fixed-structure controllers[END_REF], [START_REF] Apkarian | Tuning controllers against multiple design requirements[END_REF]. All results developed so far carry over to multiple requirements, because the worst-case multi-objective performance in step 4 of Algorithm 1 involves H -= -H which has the same min-min structure as before.

VI. EXPERIMENTS A. Algorithm testing

In this section our dynamic inner approximation technique (Algorithm 1) is tested on a bench of 14 examples of various sizes and structures. All tests have been adapted from the literature and are described in Table I. Some tests have been made more challenging by adding uncertain parameters in order to illustrate the potential of the technique for higher-dimensional parametric domains ∆. The notation [r 1 r 2 . . . r m ] in the rightmost column of the table stands for the block sizes in ∆ = diag [δ 1 I r1 , . . . , δ m I rm ]. Uncertain parameters have been normalized so that ∆ = [-1, 1] m , with nominal value δ = 0.

The dynamic inner approximation of Algorithm 1 is first compared to static inner approximation [START_REF] Burke | Stabilization via nonsmooth, nonconvex optimization[END_REF], which uses a dense enough static grid ∆ s of ∆ to perform a multi-model synthesis for a large number card(∆ s ) of models where card stands for set cardinality. In consequence, static approximation cannot be considered practical. Namely,

• Dense grids become quickly intractable for highdimensional ∆. • Static approximation may lead to overly optimistic answers in terms of worst-case performance if critical parametric configurations are missed by gridding. This is what is observed in columns 7 -9 of Table II, where we used a 5 m -point grid with m = dim(δ) the number of uncertain parameters. Worst-case performance is missed in tests 6, 9, 12 and 14, as we verified by Algorithm 1. Running times may rise to hours or even days for cases 1, 2, 5 and 10. On the other hand, when the worst case δ * found by Algorithm 1 is close to the grid ∆ s , then static approximation and Algorithm 1 find the same result. Dynamic inner approximation (Algorithm 1) can therefore be considered a cheap, and therefore very successful, way to cover the uncertainty box. The number of scenarios in ∆ a rarely exceeds 10 in our testing. Computations were performed with MATLAB R2013b on OS Windows 7 Home Premium with CPU Intel Core i5-2410M running at 2.30 Ghz and 4 GB of RAM.

The results h ∞ achieved by Algorithm 1 shown in column 4 of Table II are underestimates of the worst-case H ∞ performance on the unit cube. We therefore certify them a posteriori through the mixed µ upper bound [START_REF] Fan | Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics[END_REF], which can be computed using the routine wcgain of [START_REF]Robust Control Toolbox 5.0[END_REF]. This is Step 6 in Algorithm 1, which gives the final validation of the design. This gives an overestimate γA1 of the worst-case performance on the unit cube ∆, which we report in column 5 of Table II. Since h ∞ ≤ µ ≤ γA1 , we know that as soon as lower and upper worst-case estimates nearly coalesce, h ∞ ≈ γA1 , controllers designed using Algorithm 1 are certified on the unit cube.

Our last comparison is between Algorithm 1 and DKSYN for complex-and mixed-µ syntheses. To put comparison on an equal basis, we scale the performance channels by 1/h ∞ , where h ∞ are the values computed via Algorithm 1 (column 4 of Table II). If we now validate our method using wcgain as before, this gives the values γA1 ≥ 1 reported in column 2 of Table III. This allows to compare results to one and to decide whether a technique improves over Algorithm 1.

For the purpose of comparison we now run complex-and mixed-µ syntheses for solving min

K(s) max δ∈∆ 1 h ∞ T zw (δ, K(s)) ∞ . (23)
The optimal controllers K(s) are then again analyzed in closed-loop via wcgain, which gives the estimates γC and γR reported in columns 3 and 4 of Table III.

Table III shows that all synthesis techniques perform equally well in tests 4 and 6, as indicated by γA1 ≈ γC ≈ γR . Mixedµ synthesis performs better in case 8, but with a very high order controller, 61 states versus a simple PID for Algorithm 1. In all other cases, Algorithm 1 achieves better worst-case performance γA1 γC , γR with simpler structured controllers K(κ). It also proves competitive in terms of CPU.

Going a step further, it can be shown using the min-min program in Step 4 in Algorithm 1, or alternatively the routine wcgain, that the upper estimates γA1 , γC and γR are tight. We therefore conclude that the conservatism observed in µ synthesis is not an artifact of the analysis routine, but reflects the conservatism of the synthesis approach. The results of Tables II and III can be checked by downloading the entire benchmark at http://pierre.apkarian.free.fr/BENCH.zip

B. Tail fin controlled missile

We now illustrate our robust synthesis technique in more depth for a tail fin controlled missile. This problem is adapted from [52, Chapter IV] and has been made more challenging by adding parametric uncertainties in the most critical parameters. The linearized rigid body dynamics of the missile are

α q = Z α 1 M α M q α q + Z d M d u η q = V /kG Z α 0 0 1 α q + V /kG Z d 0 u
where α is the angle of attack, q the pitch rate, η the vertical acceleration and u the fin deflection. Both η and q are measured through appropriate devices as described below.

A more realistic model also includes bending modes of the missile structure. In this application, we have 3 bending modes whose contribution to η and q is additive and described as:

η i (s) q i (s) = 1 s 2 + 2ζω i s + ω 2 i s 2 Ξ ηi sΞ qi , i = 1, 2, 3 .
It is also important to account for actuator and detector dynamics. The actuator is modeled as a 2nd-order trans-fer function with damping 0.7 and natural frequency 188.5 rad./sec. Similarly, the accelerometer and pitch rate gyrometer are 2nd-order transfer functions with damping 0.7 and natural frequencies 377 rad./sec. and 500 rad./sec., respectively.

Uncertainties affect both rigid and flexible dynamics and the deviations from nominal are 30% for Z α , 15% for M α , 30% for M q , and 10% for each ω i . This leads to an uncertain model with uncertainty structure given as ∆ = diag δ Zα , δ Mα , δ Mq , δ ω1 I 6 , δ ω2 I 6 , δ ω3 I 6 , which corresponds to δ ∈ R 6 and repetitions [1 1 1 6 6 6] in the terminology of Table I. The controller structure includes both feed-forward K ff (s) and feedback K fb (s) actions

u c = K ff (s)η r + K fb (s) η r -η m -q m = K(s)   η r -η m q m η r   ,
where η r is the acceleration set-point and η m , q m are the detectors outputs. The total number of design parameters κ in K(κ, s) is 85, as a tridiagonal state space representation of a 12-th order controller was used. The missile autopilot is optimized over κ ∈ R 85 to meet the following requirements: • Penalization of the high-frequency rate of variation of the control signal and roll-off are captured through the constraint ||W u (s)T uηr || ∞ ≤ 1, where W u (s) is a highpass weighting W u (s) := (s/100(0.001s + 1)) 2 .

• Stability margins at the plant input are specified through the H ∞ constraint W o (s)S(s)W i (s) ∞ ≤ 1, where S is the input sensitivity function S := (I + K fb G) -1 and with static weights W o = W i = 0.4. Finally, stability and performance requirements must hold for the entire range of parametric uncertainties, where ∆ is the R 6 -hyperbox with limits in percentage given above. The resulting nonsmooth program v * to be solved in step 2 of Algorithm 1 takes the form min κ∈R 85 max δ∈∆a⊂R 6 T zw (δ, κ) ∞ . We have observed experimentally that controllers K(s) of order greater than 12 do not improve much. The order of the augmented plant including flexible modes, detector and actuator dynamics, and weighting filters is n x = 23.

The evolution of the worst-case H ∞ performance vs. iterations in Algorithm 2 (and its Subroutine 1) is problemdependent. For the missile example, a destabilizing uncertainty is found at the 1st iteration. The algorithm then settles very quickly in 5 iterations on a final set ∆ a consisting of 6 scenarios. The number of scenarios in the final ∆ a coincides with the number of iterations in Algorithm 1 plus the nominal scenario, and can be seen in column 3 of Table II. Note that the evolution of the worst-case H ∞ performance is not always monotonic. Typically the curve may bounce back when a bad parametric configuration δ is discovered by the algorithm. This is the case e.g. for the mass-spring example.

The achieved values of the H ∞ norm and corresponding running times are given in Table II. Responses to a step reference input for 100 models from the uncertainty set ∆ are shown in Fig. 3 to validate the robust design. Good tracking is obtained over the entire parameter range. The magnitude of the 3 controller gains of K(s) is plotted in Fig. 4. Robust roll-off and notching of flexible modes are clearly achieved. Potential issues due to pole-zero cancellations are avoided as a consequence of allowing parameter variations in the model. Finally, Fig. 5 displays the Nichols plots for 100 models sampled in the uncertainty set. We observe that good "rigid" margins as well as attenuation of the flexible modes over ∆ has been achieved. Remark 12. Mixed-µ synthesis turned out time-consuming, exceeding two hours in the missile example. The controller order inflates to 385 and conservatism is still present as compared to dynamic approximation via Algorithm 1 as shown in Table III. Resorting to interpreting uncertain parameters as complex cannot be considered an acceptable workaround. Even when it delivers a result, this approach as a rule leads to high-order controllers (149 states in the missile example) and is also conservative, as we expected. It appears that scalingor multiplier-based approaches outer approximations [START_REF] Young | Controller design with real parametric uncertainty[END_REF], [START_REF] Packard | Linear, multivariable robust control with a µ perspective[END_REF] encounter difficulties in terms of conservatism and controller complexity for repeated parametric uncertainties whereas our approach is not affected by these issues.

Remark 13. Static inner approximation remains intractable even for a coarse grid of 5 points in each dimension. See Table II.

VII. CONCLUSION

We have presented a novel algorithmic approach to parametric robust H ∞ control with structured controllers. A new inner approximation technique termed dynamic inner approximation, adapting a set of parameter scenarios ∆ a iteratively, was developed and shown to work rapidly without introducing conservatism. Global robustness and performance certificates are then best obtained a posteriori by applying analysis tools based on outer approximations. At the core our new method is leveraged by sophisticated nonsmooth optimization techniques tailored to the class of upper-C 1 stability and performance functions. The approach was tested on a bench of challenging examples, and within a case study. The results indicate that the proposed technique is a valid practical tool, capable of solving challenging design problems with parametric uncertainty. The new method discussed in this paper will be available in the R2015b release of MATLAB's Robust Control Toolbox.

Figure 1 .Figure 2 .

 12 Figure 1. Robust synthesis interconnection

Fig. 2 :

 2 Fig. 2: Minimization of a minimum (right) may require deciding

Remark 6 .

 6 else keep t k+1 = t k . Increase counter k and go back to step 2. Typical values are γ = 0.0001, γ = 0.0002, and Γ = 0.1. For backtracking we use θ = 1 4 and Θ = 3 4 .

•

  The acceleration η m should track the reference input η r with a rise time of about 0.5 seconds. In terms of the transfer function from η r to the tracking error e := η rη m this is expressed as ||W e (s)T eηr || ∞ ≤ 1, where the weighting function W e (s) is W e (s) := 1/M s/ω B + M s/ω B + A , A = 0.05, M = 1.5, ω B = 10 .

Fig. 3 :FFig. 4 :Fig. 5 :

 345 Fig. 3: Step responses of controlled missile for 100 sampled models in uncertainty range

TABLE I :

 I Test cases

	No.	Benchmark name	Ref. States Uncertainty block structure
	1	Flexible Beam	[43]	8	[1 1 1 3 1]
	2	Mass-Spring-Dashpot	[44] 12	[1 1 1 1 1 1]
	3	DC Motor	[45]	5	[1 2 2]
	4	DVD Drive	[46]	5	[1 3 3 3 1 3]
	5	Four Disk	[47] 10 [1 3 3 3 3 3 1 1 1 1]
	6	Four Tank	[48]	6	[1 1 1 1]
	7	Hard Disk Drive	[49] 18 [1 1 1 2 2 2 2 1 1 1 1]
	8	Hydraulic Servo	[50]	7	[1 1 1 1 1 1 1 1]
	9	Mass-Spring System	[51]	4	[1 1]
	10 Tail Fin Controlled Missile [52] 23	[1 1 1 6 6 6]
	11	Robust Filter Design 1 [53]	4	[1]
	12	Robust Filter Design 2 [54]	2	[1 1]
	13	Satellite	[55]	5	[1 6 1]
	14	Mass-Spring-Damper	[12]	8	[1]

TABLE II :

 II Comparisons of Algorithm 1 with static approximation on unit box, running times in sec.

	No. order	Algorithm 1 # scen. h∞ γA1 time # scenarios H∞ norm time Static relaxation
	1 3 2 5 3 PID 2 4 16	1.290 1.290 25 3125 2.960 2.990 261 15625 0.500 0.510 6 125	I I 0.500	∞ ∞ 128
	4 5	1	45.455 45.455 2	15625	45.454	4909
	5 6 6 6	6 4	0.682 0.682 68 9765625 I 5.568 5.568 42 625 5.564	∞ 3872
	7 4 8 PID 3 4 9 4 4	0.026 0.026 35 48828125 I 0.701 0.708 10 390625 I 0.745 0.745 23 25 0.759	∞ ∞ 67
	10 12 11 4	6 4	1.810 1.828 159 15625 2.636 2.636 17 5	I 2.636	∞ 7
	12 1	3	2.793 2.793 8	25	2.660	23
	13 6	5	0.154 48 125	0.156	876
	14 5	3	1.643 1.643 39 5	1.644	27
			I = intractable

TABLE III

 III

		: Comparisons of Algorithm 1 (A1) & DKSYN
		complex-(C) and mixed-µ (R)
	Worst-case H ∞ performance on unit box computed using µ
			upper bound
	No.	worst-case H∞ Controller order Running times γA1 γC γR A1 C R A1 C R
	1 1.00 1.36 U	3	52 102	25 80 86
	2 1.01 1.56 1.22 5	42 82	261 123 141
	3 1.02 24.86 26.10 PID 47 47	6 76 27
	4 1.00 1.00 1.00 5	5 5	2 27 53
	5 1.00 4.04 3.40 6	54 10	68 131 316
	6 1.00 1.00 1.00 6	6 22	42 18 29
	7 1.01 63.65 F	4	18 F	35 159 F
	8 1.01 1.20 0.95 PID 45 47	10 101 134
	9 1.00 1.10 1.21 4	24 28	23 48 113
	10 1.01 2.83 2.05 12 149 385	159 1412 7612
	11 1.00 1.04 1.49 4	14 16	17 14 22
	12 1.00 1.02 1.07 1	10 14	8 16 21
	13 1.00 2.76 8.41 6	269 259	48 183 257
	14 1.00 1.32 1.27 5	14 16	39 17 47
		U: unstable, F: failure, times in seconds


Adel Maamar
email: adel.maamar@chu-rennes.fr

Stéphanie Chevalier

Pierre Fillâtre

Vlad Botoc

Yves Le Tulzo

Arnaud Gacouin

Jean-Marc Tadié

Le Tulzo

Tadié Jm

COPD is independently associated with six-month survival in patients who have life support withheld in intensive care

Keywords: Withholding treatment, prognosis, patients outcome, ethics, acute respiratory failure, COPD, critically ill patients

COPD is independently associated with 6-month survival in patients who have life support withheld in intensive care

Introduction:

During the past 20 years, guidelines for the withholding (Wh) and withdrawal (Wd) of life support have been published in the United States 1 and France [START_REF] Ferrand | Société de Réanimation de Langue Française. End-of-life care in the intensive care unit. Recommendations of the Société de Réanimation de Langue Française[END_REF] . In 2005, a law in France has been voted for [START_REF]Du 22 Avril[END_REF] (n° 2005-370 of April 22, 2005), providing a legal framework for decisions to Wh or Wd life-sustaining treatments in critically ill patients. This law authorizes the withdrawal or withholding of life-sustaining therapies when they are "useless, disproportionate or have no other effect than solely the artificial preservation of life". It introduces the essential concept of futility and supports avoiding situations that may be considered as "irrational obstinacy".

Several studies have shown that the withholding and withdrawal of life support are common practices in Europe. For instance, Ferrand et al [START_REF] Ferrand | Withholding and withdrawal of life support in intensive-care units in France: a prospective survey[END_REF] reported in 2001 that, in 113 French ICUs, 53% of ICU deaths were preceded by a decision to withhold or withdraw life support and that such decisions affected approximately 10% of all patients. These findings have been confirmed in more recent studies [START_REF] Esteban | Withdrawing and withholding life support in the intensive care unit: a Spanish prospective multi-centre observational study[END_REF][START_REF] Lesieur | Withholding or withdrawal of treatment under French rules: a study performed in 43 intensive care units[END_REF][START_REF] Orban | Causes and Characteristics of Death in Intensive Care Units: A Prospective Multicenter Study[END_REF][START_REF] Azoulay | End-of-life practices in 282 intensive care units: data from the SAPS 3 database[END_REF][START_REF] Sprung | End-of-life practices in European intensive care units: the Ethicus Study[END_REF] and it appears that the withholding of life support in critically-ill patients does not necessarily lead to death in the ICU. Ferrand et al. [START_REF] Ferrand | Withholding and withdrawal of life support in intensive-care units in France: a prospective survey[END_REF] reported an ICU discharge rate of 43% after Wh, and a recent study [START_REF] Lesieur | Withholding or withdrawal of treatment under French rules: a study performed in 43 intensive care units[END_REF] reported an ICU discharge rate of 37% and a hospital discharge rate of 17%. Although in-hospital outcomes in ICU patients following decisions to Wh or Wd life support have been previously assessed, the outcomes in these patients several months after hospital discharge remain unknown. For this purpose, we conducted a retrospective study to assess the six-month outcomes of patients discharged from the ICU after life-sustaining treatments were limited and to identify prognostic factors for survival.

Methods:

Study design and measurements:

This study was conducted in a 15-bed ICU at the General Public Hospital of Saint-Malo, France. All patients who had been discharged from the ICU and experienced withheld life-support from Jan 1, 2009 to Dec 31, 2011 were included in the study. Patients with incomplete data such as comorbidities or vital status at six months were excluded from the study.

Withholding was defined in accordance with the French guidelines 2 edited by the French Intensive Care Society as a planned decision not to institute or optimize therapies that were otherwise warranted: endotracheal intubation, non-invasive ventilation, renal replacement therapy, increased doses of vasopressor infusions beyond a defined threshold, transfusion of blood products, antimicrobial therapy, surgery, increased fraction of inspired oxygen beyond a defined threshold, external cardiac massage, and ICU readmission. Data collected for all patients were as follows: age, sex, severity of illness at admission according to the simplified acute physiological score (SAPS II) [START_REF] Gall | A new Simplified Acute Physiology Score (SAPS II) based on a European/North American multicenter study[END_REF] and the sequential organ failure assessment score (SOFA) [START_REF] Vincent | Use of the SOFA score to assess the incidence of organ dysfunction/failure in intensive care units: results of a multicenter, prospective study. Working group on "sepsis-related problems" of the European Society of Intensive Care Medicine[END_REF] at admission. Life expectancy was estimated using the McCabe classification [START_REF] Mccabe | Gram-negative bacteremia: I. etiology and ecology[END_REF] : class A for no disease or a nonfatal underlying disease, class B for an ultimately fatal disease (death expected in a 5-year period), and class C for a rapidly fatal disease (death expected within 1 year). Diagnoses were coded according to the International Classification of Diseases, 10th revision (ICD-10). The primary diagnosis was defined according to the major organ dysfunction and classified into one of seven categories: acute respiratory failure without cardiac failure, coma including toxic causes, renal failure and metabolic disorders, acute pulmonary cardiogenic edema or shock, acute liver failure or digestive failure, surgical history, or sepsis [START_REF] Boumendil | Prognosis of patients aged 80 years and over admitted in medical intensive care unit[END_REF] . We also recorded the following reasons for Wd or Wh based on items proposed by the French intensive care society: futility and/or poor expected quality of life, a severe underlying condition, or family or patient request 2 . As the reason for Wd and Wh is required by French law, this item was systematically recorded in the medical records.

Additionally, we recorded the involvement of the referring physician, ICU and in-hospital length of stay (LOS), and history of previous chronic diseases including chronic respiratory insufficiency, chronic obstructive pulmonary disease (COPD) according to the American Thoracic Society criteria [START_REF]Standards for the diagnosis and care of patients with chronic obstructive pulmonary disease. American Thoracic Society[END_REF] (smoking, other environmental risk factors, chronic cough, acute chest illnesses, dyspnea, physical examination, and laboratory investigations such as chest radiography, lung function tests, or arterial blood gases); malignant disorders; class IV heart failure according to the New York Heart Association criteria; chronic neurological disease; proven cirrhosis; pre-existing renal insufficiency and severe psychiatric disorder (schizophrenia, major affective disorder, institutionalized or under protective supervision, such as tutorship and curatorship). We assessed the vital status at six months after ICU discharge for all patients by calling patients, their relatives, or their family doctor. If the patient was not alive, we recorded their length of survival. The primary aim was to assess the survival status at six months post-ICU discharge in patients with at least one treatment limitation. The secondary aim was to identify prognostic factors for six-month survival.

Statistical analysis:

Values are presented as the median (Inter Quartile Range) for continuous variables or number (percentages) for categorical variables. The two groups (survivors vs. non-survivors) were compared using Student's t-test for continuous variables and a chi-square or Fisher's exact test for categorical variables. First, univariate analysis was performed for each potential factor. Factors with a p value of less than 0.2 in the univariate analysis were then introduced as part of a Cox proportional hazard model, excluding redundant variables associated with long-term survival. Hazard ratios (HR) and their 95% confidence intervals (CI) were estimated from the estimated parameters of the final regression model. Survival curves with a 95% confidence interval were computed using the Kaplan-Meier method. Patients were censored at six months. Statistical analyses were performed using the StatView statistical software version 5.0 (SAS Institute, Cary, NC, USA); p values of less than 0.05 were considered significant.

Ethical considerations:

The study protocol was approved by the local ethical committee (15.108). Due to the nature of the study, patient informed consent was waived by the ethical committee.

Results: Overall Population:

During the study period, 1483 patients were admitted to the ICU. The mortality rate was 17% (255 deaths, 1228 discharged from the ICU) during the ICU stay. Of the 1483 patients admitted, 264 had life support Wh or Wd during their ICU stay (17.8% of all patients). Among them, 119 patients (45% of this group and 9.7% of all patients discharged from the ICU) were discharged from the ICU. Five patients were excluded because of incomplete data, specifically their comorbidities (2 patients) and their vital status at six months (3 foreign patients). Thus, 114 patients were included in the study group. The baseline characteristics of the study group are shown in Table 1. The Wh measures implemented are shown in Table 2. A total of 328 measures were taken. The previous quality of life was the most frequently reported reason for Wh life support. Endotracheal intubation was the most frequently cited type of limitation (81 patients, 71%). Of note, NIV was withheld in 5%.

Long-term survival analysis after ICU discharge:

After ICU discharge, 67 patients for whom life-sustaining therapies were withheld were alive (58.8%), and 47 (41%) patients had died within 6 months.

For the study group, the mean survival time was 124±72 days.

Characteristics of the patients in the survivor and non-survivor groups at six months posthospital discharge are shown in Table 3. Survival was significantly associated with a longer length of ICU and hospital stay, an ICU admission for acute respiratory failure, a medical history of chronic respiratory disease or chronic obstructive pulmonary disease, and a lower SAPS II score at admission. To assess the impact of life expectancy before ICU admission, we compared the Mc Cabe score of patients according to each main comorbidity (Table 4).

The multivariate analysis identified a previous history of COPD (HR=2.11; IC 95% 1.02-4.36) and a longer length of hospital stay (HR=1.03; IC 95% 1.004-1.5) as factors associated with better survival rates (Table 3).

Long-term outcomes of COPD patients after ICU discharge for whom life-sustaining therapies were withheld:

Forty-one patients were admitted for acute respiratory failure. The mean survival duration of these patients was 144 days versus 95 days for patients with other diagnoses. Eleven of these 41 patients (26.8%) had been previously admitted to the ICU, and life support was withheld (no intubation) in the previous stays of 4 patients. We compared the survival rate of patients by previous history of COPD. As shown in Figure 2, the survival rate at 6 months postdischarge was significantly higher for patients with a previous history of COPD. Among the 37 patients with COPD, 11 (29%) underwent invasive mechanical ventilation and 1 died. The remaining 26 COPD patients (71%) received non-invasive ventilation and 9 died within 6 months. There was no significant difference in the survival rates of patients who did or did not receive invasive ventilation (p=0.11).

Discussion:

The results of our study showed that more than half of patients for whom life support was withheld in the ICU survived for at least 6 months after discharge. Among these patients with Wh decision, we also found that COPD patients had a significantly higher survival rate at 6 months post-discharge than did non-COPD patients.

Withholding life-supporting therapies does not predict certain death in the ICU. In our study, 45% of patients for whom life support was withheld in the ICU were discharged. In accordance with the findings of previous studies, we found that half of the deaths in the ICU occurred after the decisions to withhold or withdraw life support [START_REF] Ferrand | Withholding and withdrawal of life support in intensive-care units in France: a prospective survey[END_REF][START_REF] Lesieur | Withholding or withdrawal of treatment under French rules: a study performed in 43 intensive care units[END_REF][START_REF] Orban | Causes and Characteristics of Death in Intensive Care Units: A Prospective Multicenter Study[END_REF][START_REF] Azoulay | End-of-life practices in 282 intensive care units: data from the SAPS 3 database[END_REF] . Among those patients for whom life support was withheld, 58.8% were still alive six months post-discharge.

Lautrette et al. [START_REF] Lautrette | Respective impact of no escalation of treatment, withholding and withdrawal of life-sustaining treatment on ICU patients' prognosis: a multicenter study of the Outcomerea Research Group[END_REF] found in 2015 a Day-30 survival of 25% after withholding life-sustaining treatment. This difference may have 2 explanations. First, they distinguished patients with withholding and with no-escalation of life-sustaining therapies for whom Day-30 survival was 65%. Second, the population of patients with chronic respiratory insufficiency was lower.

Additionally, another recent French study found an ICU discharge rate of 37% and a hospital discharge rate of 17% 6 . These authors 6 reported that withholding (not withdrawing) lifesupport therapies was associated with a mortality rate of only 56% after ICU discharge as we found in our study. Of note, the six-month survival of these patients was not assessed, however these results highlight the important survival difference between Wh and Wd, as all authors report that survival rate of Wd patients is almost nil.

In our study, we found that the survival rate in these patients was significantly higher in patients with a previous history of COPD. COPD patients having survived an acute hypercapnic respiratory failure have a poor prognosis as they are highly susceptible to readmissions, life-threatening events, and death [START_REF] Chu | Readmission rates and life threatening events in COPD survivors treated with non-invasive ventilation for acute hypercapnic respiratory failure[END_REF] , with or without Wh of life-sustaining therapies. Despite this poor prognosis, it seems better in our study than for other patients. This finding may have several explanations. First, we believe that this result arises from a better management of these patients in the COPD-treating ward, thus improving their prognosis during the post-ICU period. The death rate of COPD patients tends to stabilize [START_REF] Fuhrman | Deaths from chronic obstructive pulmonary disease in France, 1979-2002: a multiple cause analysis[END_REF] , probably because of a better management of COPD [START_REF] Rycroft | Epidemiology of chronic obstructive pulmonary disease: a literature review[END_REF] and an increased availability of hospital and home-based non-invasive ventilation (NIV) [START_REF] Plant | Cost effectiveness of ward based noninvasive ventilation for acute exacerbations of chronic obstructive pulmonary disease: economic analysis of randomised controlled trial[END_REF] . It has been shown recently that home based NIV improves the survival of COPD patients, even among severe, stable patients [START_REF] Köhnlein | Non-invasive positive pressure ventilation for the treatment of severe stable chronic obstructive pulmonary disease: a prospective, multicentre, randomised, controlled clinical trial[END_REF] . This improvement was previously shown in patients with chronic respiratory failure secondary to neuromuscular or chest wall disorders [START_REF] Annane | Nocturnal mechanical ventilation for chronic hypoventilation in patients with neuromuscular and chest wall disorders[END_REF] . A retrospective study suggested that long-term NIV could reduce the rates of hospital readmissions and acute hypercapnic respiratory failure in severe unstable COPD patients [START_REF] Ankjaergaard | Long-term non-invasive ventilation reduces readmissions in COPD patients with two or more episodes of acute hypercapnic respiratory failure[END_REF] . A recent multicentric randomized clinical trial of patients with persistent hypercapnia and hypoxemia showed that adding noninvasive ventilation to home oxygen therapy prolonged the time to readmission or death within 12 months [START_REF] Murphy | Effect of Home Noninvasive Ventilation With Oxygen Therapy vs Oxygen Therapy Alone on Hospital Readmission or Death After an Acute COPD Exacerbation: A Randomized Clinical Trial[END_REF] . Unlike other life support therapies that are only available in the ICU (such as vasopressor infusion), NIV is the only organ support therapy that is available daily at the same intensity at home, in the ward and in the ICU. Azoulay et al. showed that hospital mortality rates for COPD patients admitted for acute respiratory failure were lower in comparison with other patients in both do-not-intubate group and no treatment limitation group [START_REF] Azoulay | Noninvasive mechanical ventilation in patients having declined tracheal intubation[END_REF] . They also showed that there was no significant decline in their quality of life at day-90, and that 56% of patients with a donot-intubate status were alive at hospital discharge. Second, according to recent studies, invasive mechanical ventilation worsens the short-and long-term prognoses of patients with acute-on-chronic respiratory insufficiency [START_REF] Gacouin | Trends in Prevalence and Prognosis in Subjects With Acute Chronic Respiratory Failure Treated With Noninvasive and/or Invasive Ventilation[END_REF][START_REF] Gacouin | Long term survival in chronic obstructive pulmonary disease after ventilation in intensive care[END_REF] . It also has been shown to be of no benefit in the most severe patients [START_REF]Good care of the dying patient[END_REF][START_REF] Singer | Rationing, patient preferences, and cost of care at the end of life[END_REF] . As endotracheal intubation was the most frequent life-support therapies withheld, we hypothesize that the COPD patients who had withheld life-support therapy during their ICU stay still benefited from NIV in ward and ICU during the following acute hypercapnic respiratory failure. Of note, although other studies used day-30 or day-90 mortality, we decided to study 6 months mortality because of the high mortality rate expected and to avoid loss to follow-up Our study highlights the importance of individualization of post-ICU follow-up care to improve patient survival and to avoid new exacerbations. Furthermore, we think that the decision to withhold life-support should not lead to not readmit these patients and should not refrain from treating the COPD patients for whom life-support was withheld.

We also found in our study that the survivor patients had a significantly longer hospital stay after ICU discharge probably due to a better general condition at discharge.In addition to the retrospective nature of the study, some limitations should be acknowledged. First, we can't exclude uncontrolled confounders or case-mix. Second, this is a single-center study, thereby limiting its external validity. Furthermore, we only considered COPD patients in our study since its represent the major etiology of chronic respiratory disease in our studied population.

However, other chronic respiratory diseases such as idiopathic pulmonary fibrosis should be evaluated. Third, management of patients after ICU discharge should has been assessed, especially in COPD patients for whom home noninvasive ventilation prolonged the time to readmission or death within 12 months after an acute exacerbation [START_REF] Murphy | Effect of Home Noninvasive Ventilation With Oxygen Therapy vs Oxygen Therapy Alone on Hospital Readmission or Death After an Acute COPD Exacerbation: A Randomized Clinical Trial[END_REF] . Fourth, we can't exclude confounding factors that have not been assessed in this specific population such as sources of stress which could be responsible for early readmission [START_REF] Edmondson | Psychological stress and 30day all-cause hospital readmission in acute coronary syndrome patients: an observational cohort study[END_REF] . Lastly, we do not have quality of life data for surviving patients, which would have been interesting. Recent studies have emphasized that the quality of life of critically-ill survivor patients was impaired when compared to that of the general population [START_REF] Angus | Surviving intensive care: a report from the 2002 Brussels Roundtable[END_REF] , especially after septic shock [START_REF] Nesseler | Long-term mortality and quality of life after septic shock: a follow-up observational study[END_REF] , and that improving quality of life following ICU discharge is an important issue.

Conclusion:

In conclusion, our study found that a substantial proportion of patients discharged from the ICU after life-sustaining therapies were withheld survived for at least six months. Among them, patients with COPD appeared to have a significantly higher survival rate. These results highlight the notion that the decision to withhold life-sustaining therapies should not lead to a cessation of post-ICU follow-up care and that COPD patients with a decision to withhold endotracheal intubation may still benefit of ICU-readmission.

Figure 1

 1 shows the flowchart of the selection of included study participants. All patients from whom life support was Wd died in the ICU (145 patients). Of the 1483 patients admitted, 1219 had no treatment limitations and 110 died within six months (9%) after ICU discharge. Studied population: patients discharged from the ICU after a decision to withhold lifesustaining therapies.

This study was partly presented at the