H H Bauschke

D Noll

ON THE LOCAL CONVERGENCE OF THE DOUGLAS-RACHFORD ALGORITHM

Keywords: nonconvex feasibility problem, fixed-point, discrete dynamical system, convergence, stability

We discuss the Douglas-Rachford algorithm to solve the feasibility problem for two closed sets A, B in R d . We prove its local convergence to a fixed point when A, B are finite unions of convex sets. We also show that for more general nonconvex sets the scheme may fail to converge and start to cycle, and may then even fail to solve the feasibility problem.

Introduction

The Douglas-Rachford projection algorithm, originally introduced in [START_REF] Douglas | On the numerical solution of heat conduction problems in two or three space variables[END_REF] to solve nonlinear heat flow problems, aims to find a point x * in the intersection of two closed constraint sets A, B in R d or in Hilbert space. As a consequence of more general results in monotone operator theory by Lions and Mercier [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF], it is known that the scheme converges weakly for two closed convex sets A, B in Hilbert space with non-empty intersection. A rather comprehensive analysis of the convex case is given in [START_REF] Bauschke | Finding best approximation pairs relative to two closed convex sets in Hilbert spaces[END_REF].

Due to its success in applications, the Douglas-Rachford scheme is frequently used in the nonconvex setting despite the lack of a satisfactory convergence theory. Recently Hesse and Luke [START_REF] Hesse | Nonconvex notions of regularity and convergence of fundamental algorithms for feasibility problems[END_REF] made progress by proving local convergence of the scheme for B an affine subspace intersecting the set A transversally, where A is no longer convex, but satisfies a regularity hypothesis called superregularity. Numerical experiments in the nonconvex case indicate, however, that the Douglas-Rachford scheme should converge in much more general situations. Very frequently one observes that the iterates settle for convergence after a chaotic transitory phase; see [START_REF] Aragón Artacho | Recent results on Douglas-Rachford methods[END_REF][START_REF] Elser | Searching with iterated maps[END_REF] and the references therein. Here we prove local convergence of the Douglas-Rachford scheme when A, B are finite unions of convex sets. Our result is complementary to [START_REF] Hesse | Nonconvex notions of regularity and convergence of fundamental algorithms for feasibility problems[END_REF], because no transversality hypothesis is required. This result is proved in section 3.

We will also show that for nonconvex sets A, B the Douglas-Rachford scheme may fail to converge and start to cycle without solving the feasibility problem. We show that this may even lead to continuous limiting cycles. These are more delicate to construct, because in that case the Douglas-Rachford sequence x n+1 ∈ T (x n ) is bounded and satisfies x n+1 -x n → 0, but fails to converge. Our construction is given in section 4.

Preparation

Given a nonempty closed subset A of R d , the projection onto A is the set-valued mapping P A associating with x ∈ R d the nonempty set P A (x) = a ∈ A : x -a = dist(x, A) , Date: June 2, 2014. † University of British Columbia, Kelowna, Canada. * Institut de Mathématiques de Toulouse, France.

where • is the Euclidean norm, and dist(x, A) = min{ x -a : a ∈ A}. The reflection of x in A is then the set-valued operator

R A = 2P A -I.
Given two nonempty closed sets A, B in R d , the Douglas-Rachford iterative scheme, starting at x 0 , generates a sequence x n by the recursion

x n+1 ∈ T (x n ), T := 1 2 (R B R A + I) , n ∈ N.
We call T the Douglas-Rachford operator, or shortly, DR operator. Suppose x + ∈ T (x) is one step of the Douglas-Rachford scheme. Then x + is obtained as

x + = x + b -a,
where a ∈ P A (x), y = 2a -x, and b ∈ P B (y). We call a the shadow of iterate x on A, b the reflected shadow of x in B, both used to produce x + . We write x + = T (x) if the DR-operator is single-valued, and similarly for projectors P A , P B and reflectors

R A , R B .
The fixed point set of T is defined as

F(T ) = {x ∈ R d : x ∈ T (x)}. Note that if x * ∈ F(T )
, and if a * ∈ A, b * ∈ B are the shadow and reflected shadow of x * used to produce x * ∈ T (x * ), then a * = b * ∈ A ∩ B, so every fixed point gives rise to a solution a * ∈ A ∩ B of the feasibility problem. However, in the set-valued case, it may happen that x * ∈ F(T ) has other shadow-reflected shadow pairs (ã, b) leading away from x * , i.e., where ã = b, so that x = x * + b -ã ∈ T (x * ) \ {x * }. We therefore introduce the set of strong fixed-points of T as

F (T ) = x ∈ R d : T (x) = {x} . Note that A ∩ B ⊂ F (T ) ⊂ F(T ). If T is single-valued, then F(T ) = F (T ).
These concepts are linked to discrete dynamical system theory, where fixed points are steady states, and where any sequence x n+1 ∈ T (x n ) is called an orbit or a trajectory. We recall that a steady state x * is stable in the sense of Lyapunov if for every > 0 there exists δ > 0 such that every trajectory x n+1 ∈ T (x n ) with starting point x 0 ∈ B(x * , δ) satisfies x n ∈ B(x * , ) for all n. Here and throughout B(x, r) means the closed Euclidian ball with centre x and radius r. It is clear that x * ∈ F(T ) \ F (T ) can never be stable, because x 0 = x * produces trajectories going away from x * .

Unions of convex sets

In this section A = i∈I A i and B = j∈J B j are finite unions of closed convex sets, a case which is of interest in a number of practical applications like rank or sparsity optimization [START_REF] Hesse | Projection methods for sparse affine feasibility: results and counterexamples[END_REF], or even road design [START_REF] Bauschke | Projection methods: Swiss Army knives for solving feasibility and best approximation problems with halfspaces[END_REF], where finite unions of linear or affine subspaces are used. For every i ∈ I and j ∈ J let T ij be the Douglas-Rachford operator associated with the sets A i , B j . By convexity of A i , B j , the operators T ij are single-valued, and

T (x) ⊂ {T ij (x) : i ∈ I, j ∈ J}.
Since A, B are finite unions of convex sets, every DR step is realized as the DR step of one of the operators T ij , and in that case, we say that this operator is active. To make this precise, we define the set of active indices at x as

K(x) = (i, j) ∈ I × J : P A i (x) ∈ P A (x), P B j (R A i (x)) ∈ P B (R A i (x)) . (1)
Note that if (i, j) ∈ K(x), then T ij (x) ∈ T (x). Conversely, for every x + ∈ T (x), there exists (i, j) ∈ K(x) such that x + = T ij (x). However, be aware that T ij (x) ∈ T (x) may be true without (i, j) being active at x.

Theorem 1. (Stable local attractor)

. Let A = i∈I A i and B = j∈J B j be finite unions of closed convex sets, and let x * ∈ F (T ) be a strong fixed point. Then x * has a radius of attraction r > 0 with the following property: For arbitrary fixed ∈ (0, r), suppose a Douglas-Rachford trajectory x n+1 ∈ T (x n ) enters the ball B(x * , ). Then it stays there and converges to some fixed point x ∈ F(T ). Moreover, every accumulation point of the shadow sequence a n ∈ P A (x n ) is a solution of the feasibility problem. The radius of attraction can be computed as

r = sup > 0 : K (B(x * , )) ⊂ K(x * ) . (2)
Proof. 1) The fact that x * ∈ F (T ) is a strong fixed point has the following consequence. Whenever (i, j) ∈ K(x * ) is active, then P A i (x * ) = P B j (R A i (x * )) ∈ A ∩ B. Therefore, for every (i, j) ∈ K(x * ), the operator T ij has x * as a fixed point.

2) We now show the following. There exists > 0 such that every x ∈ B(x * , ) has K(x) ⊂ K(x * ).

Let us consider the set I(x) = {i ∈ I : there exists j ∈ J such that (i, j) ∈ K(x)} of active indices i ∈ I at x. Then by definition

δ 1 := min dist(x * , A i ) : i ∈ I(x * ) -dist(x * , A) > 0. (3)
Similarly, we have

δ 2 := min dist (R A i (x * ), B j ) -dist (R A i (x * ), B) : i ∈ I(x * ), (i, j) ∈ K(x * ) > 0. ( 4
)
Choose > 0 such that 2 < min{δ 1 , δ 2 }. We show that is as claimed. Pick (i, j) ∈ K(x * ). We have to show that (i, j) ∈ K(x). First consider the case where i ∈ I \ I(x * ). We show that i ∈ I \ I(x). Indeed,

dist(x, A) ≤ x -x * + dist(x * , A) ≤ + dist(x * , A i ) -δ 1 (using (3)) ≤ 2 + dist(x, A i ) -δ 1 < dist(x, A i ),
showing that i ∈ I(x). We now discuss the case where i ∈ I(x * ), but (i, j)

∈ K(x * ). That means dist(R A i (x * ), B j ) -dist(R A i (x * ), B) ≥ δ 2 . Therefore dist(R A i (x), B) ≤ R A i (x) -R A i (x * ) + dist(R A i (x * ), B) ≤ + dist(R A i (x * ), B j ) -δ 2 (using (4)) ≤ 2 + dist(R A i (x), B j ) -δ 2 < dist(R A i (x), B j ), proving (i, j) ∈ K(x).
3) As an immediate consequence of 2) we have the following: If x ∈ B(x * , ) and x + ∈ T (x) is realized as x + = T ij (x) for some active operator T ij , that is, for some (i, j) ∈ K(x), then this operator T ij has x * as a fixed point. Namely, by 2) x satisfies K(x) ⊂ K(x * ), hence (i, j) ∈ K(x * ), and therefore

P A i (x * ) = P B j (R A i (x * )), which proves what we claimed.
5) We next show that as soon as a DR sequence x n+1 ∈ T (x n ) enters the ball B(x * , ), then it stays there and converges.

This can be seen as follows. Suppose the trajectory enters B(x * , ) at stage n. Then the active operator T injn used to produce [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Space[END_REF]Prop. 4.21], this operator satisfies

x n+1 = T injn (x n ) ∈ T (x n ) has x * as a fixed point, because (i n , j n ) ∈ K(x n ) ⊂ K(x * ). Therefore, by
x n+1 -x * = T injn (x n ) -x * ≤ x n -x * ≤ .
The conclusion is that from index n onward the sequence x n stays in the ball B(x * , ), and all operators T imjm used from here on have x * as a common fixed point. Now we invoke Elsner et al. [9, Thm. 1], who show that x n converges to a common fixed point x of the operators T imjm , m ≥ n. But x is then also a fixed point of T , as follows from the continuity of the distance functions. One has x ∈ B(x * , ), and moreover, if a n = P in (x n ) ∈ A in ⊂ A, then every accumulation point of the sequence a n is a solution of the feasibility problem. Namely, if we consider b n = P B j (R A i (x n )) ∈ B, and if we take accumulation points a * of a n and b * of b n , then a * ∈ P

A (x * ), b * ∈ P B (2a -x * )), hence a * = b * , because x * is a strong fixed point.
6) To conclude let us now define the radius of attraction r as in formula (2). In 1) -5) above we have shown that there exists > 0 such that K (B(x * , )) ⊂ K(x * ), and that this inclusion alone already implies convergence of every trajectory entering B(x * , ). This means that the supremum in ( 2) is over a nonempty set, and that is all that we need.

Remark 1. (Stable steady state.) The dynamic system interpretation of Theorem 1 is that a strong fixed point x * ∈ F (T ) is a stable steady state of the Douglas-Rachford dynamical system x + ∈ T (x) when A, B are finite unions of convex sets. Note also that we do not claim that x ∈ F(T ) is strong, nor do we claim that the iterates converge to x * itself.

The following observation is also of the essence.

Remark 2. (Strong fixed point needed). Theorem 1 is not true if x * ∈ F(T ) \ F (T ), that is, if x * is not a strong fixed point. Indeed, let A = {-1, 1} and B = {-2, 1}. Then x * =
0 is a fixed-point of T , but not a strong one. Now there exist arbitrarily small ∈ (0, 1) such that trajectories starting in B(0, ) = (-, ) will not stay in that ball. Indeed, for x ∈ (-, 0), x + will move away from 0 and will not stay in B(0, ), while for x ∈ (0, ), x + stays in that interval.

Remark 3. Note that if A, B are convex sets, then all K(x) are identical singleton sets, so formula (2) gives r = ∞, which means the DR scheme converges globally. x * x + y * } a line, where x * > 0, y * > 0. Then (0, y * ) and (x * , 0) are the two only fixed points of T , both strong, and one easily finds r(x * , 0) = x * / √ 2 and r(0, y * ) = y * / √ 2.

Remark 5. (Asymptotic stability). Let x * ∈ F (T ) be a strong fixed point, and suppose there exists δ > 0 such that B(x * , δ) contains no further fixed point of T . Then it follows from Theorem 1 that we can find ∈ (0, δ] such that every trajectory x n+1 ∈ T (x n ) entering B(x * , ) stays there and converges to x * . In the dynamical system terminology, x * is then asymptotically stable in the sense of Lyapunov. Note that this still fails for an isolated fixed point x * ∈ F(T ) \ F (T ).

Theorem 2. (Local convergence)

. Let A = i∈I A i and B = j∈J B j be finite unions of convex sets. Let x n+1 ∈ T (x n ) be a bounded Douglas-Rachford sequence satisfying x n+1 -x n → 0. Then x n converges to a fixed-point x ∈ F(T ). Moreover, every accumulation point of the shadow sequence a n ∈ P A (x n ) is a solution to the feasibility problem.

Proof. 1) For every n ∈ N let us choose an active index pair

(i n , j n ) ∈ K(x n ) such that x n+1 = T injn (x n ). Put a n = P A in (x n ) and b n = P B jn R A in (x n ) = P B jn (2a n -x n ), so that x n+1 = x n + b n -a n . Note that a n -b n = x n -x n+1 → 0 by hypothesis.
2) Let x * be any accumulation point of the sequence x n . We define a subset K 0 (x * ) of the active set K(x * ) as

K 0 (x * ) = (i, j) ∈ K(x * ) : P A i (x * ) = P B j (R A i (x * )) .
Note that every T ij with (i, j) ∈ K 0 (x * ) has x * as a fixed point.

3) We now claim that for every accumulation point x * of the sequence x n there exists > 0 and an index n 0 such that for every x n with n ≥ n 0 and x n ∈ B(x * , ), we have

(i n , j n ) ∈ K 0 (x * ).
To prove this, assume on the contrary that for every = 1 k there exists n k such that

x n k ∈ B(x * , 1 k ), but (i n k , j n k ) ∈ K 0 (x * ). Moreover, let n k < n k+1 → ∞. Then x n k → x * .
Passing to another subsequence which we also denote by x n k , we may assume that

i n k = i, j n k = j. Then a n k = P A i (x n k ) → a * ∈ A, b n k = P B j (R A i (x n k )) → b * ∈ B.
Since a n -b n → 0 by part 1), we deduce that a * = b * ∈ A ∩ B. Since we also have

a n k = P A i (x n k ) ∈ P A (x n k ) and b n k ∈ P B (R A i (x n k )) ⊂ P B (R A (x n k ))
, we get a * ∈ P A (x * ) and b * ∈ P B (R A (x * )), hence (i, j) ∈ K(x * ). Since a * = b * , we have (i, j) ∈ K 0 (x * ). This contradiction proves the claim.

4) Since x * is an accumulation point of the sequence x n , there exist infinitely many indices with x n ∈ B(x * , ). Choose one with n ≥ n 0 . Then x n+1 = T injn (x n ), and by part 3) we have (i n , j n ) ∈ K 0 (x * ). By part 2), T injn has therefore x * as a fixed point. Using [2, Prop. [START_REF] Bauschke | Restricted normal cones and the method of alternating projections: theory[END_REF].21], we deduce x n+1 -x * = T injn (x n ) -x * ≤ x n -x * ≤ , hence x n+1 stays in the ball B(x * , ). This means we can repeat the argument, showing that the entire sequence x m , m ≥ n, stays in B(x * , ). By part 3) the operators T imjm , m ≥ n, have the common fixed point x * , hence we conclude again using [9, Thm. 1] that x m converges to some x, which must be a fixed point of T . The second part of the statement follows now from a n -b n → 0. Remark 6. (Discrete limit cycle). Let A = {(x, y) : y = 0} ⊂ R 2 be the x-axis, fix 1 ≥ η > 0, and put B = {(0, 0), (7 + η, η), (7, -η)}. When started at x 1 = (7, η), the method cycles between the four points x 1 , x 2 = (7 + η, 0), x 3 = (7 + η, -η), x 4 = (7, 0). Note that B is a finite union of bounded convex sets and A is convex, the iterates x 2 and x 4 reach A, but the method fails to converge, and it also fails to solve the feasibility problem.

Remark 7. (Several shadows).

Let B be a circle in R 2 , and let A consist in the union of two circles which touch B from outside in a * 1 , a * 2 . Then the centre x * of B is a fixed-point of the Douglas-Rachford operator T , and the two points a * i ∈ A ∩ B are both shadows of x * . This shows that even in the case of convergence of x n we do not expect the shadows a n to converge. i=1 C i = ∅, the Douglas-Rachford algorithm in product space can be used to compute a point in this intersection. The interesting observation is that if each C i is a finite union of convex sets, then this remains true for the set B, hence our convergence theory applies.

Existence of a continuous limit cycle

In this section we construct two closed bounded sets A, B such that the Douglas-Rachford iteration x n+1 ∈ T (x n ) with T = 1 2 (R B R A + I) fails to converge and produces a continuum of accumulation points F ⊂ F (T ) forming a continuous limit cycle. We let A be the cylinder mantle

A = (cos t, sin t, h) : 0 ≤ t ≤ 2π, 0 ≤ h ≤ 1 ,
and B a double spiral consisting of two logarithmic spirals in 3D winding down against the cylinder, one from inside, one from outside. That is, B = ((1 ± e -t ) cos t, (1 ± e -t ) sin t, e -t/2 ) : 0 ≤ t ∪ F, where F = {(cos α, sin α, 0) : α ∈ [0, 2π]}. Note that A ∩ B = F . We will construct a Douglas-Rachford sequence x n+1 = T (x n ), whose set of accumulation points is the entire set F . It will be useful to divide the spiral in its outer and inner part B ± = ((1 ± e -t ) cos t, (1 ± e -t ) sin t, e -t/2 ) : 0 ≤ t ∪ F so that B = B -∪ B + and B -∩ B + = F with these notations. Theorem 3. (Continuous limit cycle). Let x n+1 = T (x n ) be any Douglas-Rachford sequence between A and B with starting point x 1 ∈ B -\ F . Then the sequence x n is bounded, satisfies, x n+1 -x n → 0, but fails to converge. Its set of accumulation points is

F = A ∩ B ⊂ F (T ).
Proof. 1) For t ≥ 0 let us introduce the notations a(t) = cos t, sin t, e -t/2 ∈ A, and b ± (t) = (1 ± e -t ) cos t, (1 ± e -t ) sin t, e -t/2 ∈ B ± \ F. The set {a(t) : t ≥ 0} ⊂ A is the shadow of the spiral on the cylinder mantle. Namely, it is clear that for t > 0,

P A (b + (t)) = P A (b -(t)) = a(t). (5) In particular, b ± (t) -P A (b ± (t)) = e -t . (6) In consequence R A (b + (t)) = b -(t), R A (b -(t)) = b + (t).
In words, the two branches B ± of the double spiral are the reflections of each others in the cylinder mantle.

2) Let us now analyze the projection of a(t) on the double spiral B. We consider the projections of a(t) on each of the branches B ± . We start with the analysis of b ∈ P B + (a(t)). We first claim that b ∈ F . Indeed, the point v ∈ F closest to a(t) is v = (cos t, sin t, 0) = P F (a(t)), so v -a(t) = e -t/2 . But b + (t) -a(t) = e -t < e -t/2 , hence there are points on B + \ F closer to a(t) than v. This shows that any projected point b ∈ P B + (a(t)) has to be of the form b + (τ ) for some τ ≥ 0. Now consider some such b + (τ ) ∈ P B + (a(t)), then

e -t = a(t) -b + (t) ≥ a(t) -b + (τ ) ≥ a(τ ) -b + (τ ) = e -τ , (7)
which shows τ ≥ t. Here the second estimate follows from a(τ ) = P A (b + (τ )).

3) Let us further observe that τ > t. Namely, if we had τ = t, then we would have a fixed point pair for the method of alternating projections between A and B + in the sense that a(t)

= P A (b + (t)), b + (t) ∈ P B + (a(t)). That would mean the distance squared τ → 1 2 a(t) -b + (τ ) 2 had a local minimum at τ = t.
But the derivative of this function at τ = t is -e -2t < 0, so τ = t is impossible, and we deduce τ > t.

4) Using b + (τ ) ∈ P B + (a(t)) and ( 7), we find

e -t ≥ a(t) -b + (τ ) ≥ |e -t/2 -e -τ /2 | = e -t/
2 -e -τ /2 = e -t e t/2 -e t-τ /2 , which shows 0 < e t/2 -e t-τ /2 ≤ 1. This can be re-arranged as

0 < 1 -e t/2-τ /2 ≤ e -t/2 . (8)
In particular, for t → ∞ we must have τ -t → 0.

5) Let us next show that the projection b + (τ ) = P B + (a(t)) is unique for t sufficiently large. Indeed, suppose we find t < τ 1 < τ 2 such that b + (τ 1 ), b + (τ 2 ) ∈ P B + (a(t)). Then by ( 8) we have t < τ 1 < τ 2 < t -2 log 1 -e -t/2 . Define the function d

+ (x) = 1 2 a(t) - b + (x) 2 , then d + (t) = 1 2 e -2t , d + (t) = -e -2t < 0. Since τ 1 , τ 2 are local minima, we have d + (τ 1 ) = d + (τ 2 ) = 0. But d + (x) = cos(t -x) + 2e -2x + 2e -x sin(t -x) + 3
2 e -x -1 4 e -x/2-t/2 . In consequence, for t large and x moving in the interval x ∈ t, t -2 log 1 -e -t/2 , we have d + (x) ≈ cos(t-x) ≈ 1, so certainly d + (x) > 0 for these x, and since the local minima τ i are in that interval for t large, d + (τ 2 ) = 0 is impossible. This proves b + (τ ) = P B + (a(t)) for t sufficiently large. 2 e -x (1 -cos(x -t)), so f ≥ 0, and f = 0 for x = t + 2kπ. Since t < τ < t -2 log(1 -e -t/2 ) t + 2π, this proves f (τ ) > 0.

7) Let b ∈ P B (a(t)) = P B -(a(t)) a projected point of a(t) in the inner spiral. We know already that b ∈ F , hence b = b -(σ) for some σ ≥ 0. Repeating the argument in part 2), it follows that σ > t. Indeed, like in [START_REF] Douglas | On the numerical solution of heat conduction problems in two or three space variables[END_REF] we have

e -t = a(t) -b -(t) ≥ a(t) -b -(σ) ≥ a(σ) -b -(σ) = e -σ ,
and the same argument as in part 2) shows σ > t. But then again e -t ≥ a(t) -b -(σ) ≥ |e -t/2 -e -σ/2 | = e -t/2 -e -σ/2 = e -t e t/2 -e t-σ/2 , which shows 0 < e t/2 -e t-σ/2 ≤ 1. This can be re-arranged to 0 < 1 -e t/2-σ/2 ≤ e -t/2 . (9)

In particular, for t → ∞ we must have σ -t → 0, and in particular 0 < σ -t 2π. Therefore projected points b -(σ) ∈ P B (a(t)) lie on the same tour of the spiral as a(t), b -(t), and one does not take shortcuts by jumping down a full turn of the spiral B - or more. Repeating the argument of 5), we also see that b -(σ) = P B -(a(t)) is unique for t > 0 large enough.

8) Let us now generate our Douglas-Rachford sequence x n , starting at x 1 = b -(t 1 ) ∈ B - with t 1 > 0, excluding t 1 = 0 for simplicity to have a unique projection on the cylinder mantle at the start.

We get a 1 = P A (x 1 ) = a(t 1 ), hence R A (x 1 ) = b + (t 1 ). Since b + (t 1 ) ∈ B, it is its own reflection in B, and we get b 1 = b + (t 1 ). Averaging then gives x 2 = (b 1 +x 1 )/2 = a 1 = a(t 1 ), which concludes the first DR-step.

The second DR-step proceeds as follows. Since x 2 = a(t 1 ) ∈ A, it is its own reflection in A, so a 2 = x 2 = a 1 = R A (x 1 ). Now let b 2 = P B (R A (x 2 )), then b 2 = P B (a(t 1 )) = P B -(a(t 1 )), so b 2 = b -(t 2 ) for some t 2 > t 1 , where t 2 is for t 1 what σ was for t in part 7). So the reflected point is 2b -(t 2 ) -a(t 2 ) and averaging then gives

x 3 = b -(t 2 ) ∈ B -.
Proceeding in this way, we generate a strictly increasing sequence t n such that

x 2k-1 = b -(t k ) ∈ B -, x 2k = a(t k ) ∈ A. (10)
Moreover, the shadow and reflected shadow are

P A (x 2k-1 ) = a 2k-1 = a(t k ), P B (R A (x 2k-1 )) = b 2k-1 = b + (t k ) respectively, P A (x 2k ) = a 2k = a(t k ), P B (R A (x 2k )) = b 2k = b -(t k )
. Furthermore, note that we generate a sequence of alternating projections between A and B -. Namely

b 0 := x 1 P A → a 1 = a 2 P B - → b 2 P A → a 3 = a 4 P B - → b 4 . . . (11)
9) We now argue that t n so constructed tends to ∞. Suppose on the contrary that t n < t n+1 → t * < ∞. Then from the construction we see that we create a pair a(t * ) ∈ A, b -(t * ) ∈ B -such that a(t * ) ∈ P A (b -(t * )) and b -(t * ) ∈ P B -(a(t * )). Arguing as before, this would imply that τ → 1 2 a(t * ) -b -(τ ) 2 had a local minimum at τ = t * , which it does not because its derivative at t * is -e -2t * < 0. Hence t * < ∞ is impossible, and we have t n → ∞. As a consequence, the statements about uniqueness of the operators and the estimate [START_REF] Douglas | On the numerical solution of heat conduction problems in two or three space variables[END_REF] are now satisfied from some counter n 0 onward. 10) To conclude, observe that x n -x n+1 → 0 by [START_REF] Elsner | Convergence of sequential and asynchronous nonlinear paracontractions[END_REF], and that the a(t n ) are 2e -tn -dense in the interval [t n , t n + 2π], because of

a(t n ) -a(t n+1 ) ≤ a(t n ) -b -(t n+1 ) + b -(t n+1 ) -a(t n+1 ) ≤ e -tn + e -t n+1 ≤ 2e -tn .
Using [START_REF] Bauschke | Projection methods: Swiss Army knives for solving feasibility and best approximation problems with halfspaces[END_REF] and the fact that every a(t) with t ∈ [t n , t n + 2π] is at distance ≤ e -tn/2 to the set F , we deduce that every point in F is an accumulation point of both the DR-sequence (10) and the MAP sequence [START_REF] Hesse | Nonconvex notions of regularity and convergence of fundamental algorithms for feasibility problems[END_REF].

Remark 9. (Strong fixed points need not be stable). The system-theoretic interpretation of this result is that a strong fixed-point x * ∈ F ⊂ F (T ) need not be a stable steady state. This is in contrast with Theorem 1, where this was shown to be true when A, B are finite unions of convex sets. A second interpretation is that F is a stable attractor for the dynamical system x + = T (x).

Remark 10. (Shadows need not converge). We note that not only does the DR sequence x n fail to converge in Theorem 3, also the sequences a n = P A (x n ) ∈ A, b n = P B (R A (x n )) ∈ B fail to converge and have the same continuum set of accumulation points F . Presently no example of failure of convergence of a bounded DR-sequence x n satisfying x n -x n+1 → 0 is known, where the shadow sequence a n converges to a single limit a * ∈ A ∩ B. It is clear that this could only happen when F ⊂ {x : x -a * = } for some > 0. Proof. In part 2) of the proof of Theorem 3 we analyzed this sequence, which is generated by the building blocks a(t) → b + (τ ) → a(τ ).

Remark 11. Here we have an example of a semi-algebraic convex set conv(A), and the spiral B + , which is the projection of a semi-analytic set in R 4 , where the MAP sequence fails to converge and leads to a continuous limit cycle. The first example with a continuous limit cycle appears in [START_REF] Bauschke | On cluster points of alternating projections[END_REF], but with more pathological sets A, B. The fact that B + is not subanalytic can be deduced from [START_REF] Noll | On local convergence of the method of alternating projections[END_REF]Cor. 7]. Currently we do not have an example where the DR-algorithm fails to converge and creates a continuous limit cycle with one of the sets convex.

Concluding remarks

Without convexity convergence of the DR-projection algorithm and of its companion, the method of alternating projections (MAP), are usually based on transversality assumptions, see [START_REF] Bauschke | Restricted normal cones and the method of alternating projections: theory[END_REF] for MAP, and [START_REF] Hesse | Nonconvex notions of regularity and convergence of fundamental algorithms for feasibility problems[END_REF] for DR. The first results for MAP which do not rely on transversality are the recent [START_REF] Bauschke | On cluster points of alternating projections[END_REF][START_REF] Noll | On local convergence of the method of alternating projections[END_REF]. For the DR-scheme the present work gives a first convergence result without transversality.

The difficulty to prove convergence for the DR-scheme may to some extent be highlighted by inspecting the axiomatic global convergence method of Zangwill [START_REF] Zangwill | Nonlinear Programming[END_REF]. It assures convergence under fairly general hypotheses as soon as a Lyapunov function exists, but no such Lyapunov function is currently known for the DR-method.

Remark 4 .

 4 Formula (2) allows to compute the radius of attraction of a strong fixed-point x * ∈ F (T ) in certain cases. For illustration, consider A = {(x, 0) : x ∈ R} ∪ {(0, y) : y ∈ R} a union of two lines and B = {(x, y) : y = -y *

Remark 8 .

 8 (More than two sets.) It is a standard procedure in applications to extend the Douglas-Rachford scheme to solve the feasibility problem for a finite number of constraint sets C 1 , . . . , C m in R d . One defines A to be the diagonal inR d × • • • × R d (mtimes), and chooses as B = C 1 × • • • × C m in the product space. Then if m

6 )

 6 Let us now consider the point b -(τ ) ∈ B -\ F on the inner spiral. We claim that b -(τ ) is closer to a(t) than b + (τ ). Indeed, the set of points w having equal distance to b -(τ ) and b + (τ ) is the tangent plane to the cylinder at the point a(τ ) = 1 2 (b -(τ ) + b + (τ )). But the cylinder lies in one of the half-spaces associated with this plane, namely the one containing b -(τ ), b -(τ ) -a(t) < b + (τ ) -a(t) . Since b + (τ ) is the nearest point to a(t) in B + , we deduce that P B (a(t)) ⊂ P B -(a(t)) for all t. In other words, projections from the shadow of the spiral onto the double spiral always go to the inner spiral. We could also use an analytic argument to prove this. Let d -(x) = 1 2 a(t) -b -(x) 2 and consider the function f (x) = d + (x) -d -(x), then f (x) = 1

Corollary 1 .

 1 (Continuous limit cycle for MAP). Let x n be the Douglas-Rachford sequence constructed above, and let a n ∈ A, b n ∈ B be the shadows associated with x n . Then b 2n , a 2n+1 is a sequence of alternating projections between the cylinder A and the inner spiral B -. This sequence also fails to converge and has the same set of accumulation points F = A ∩ B.Corollary 2. (Limit cycle for MAP with one set convex). Every sequence of alternating projections a n , b n between the outer spiral B + and the solid cylinder conv(A) started at b 1 ∈ B + \ F is bounded, satisfies a n -b n → 0, but fails to converge and has the set F = B + ∩ conv(A) as its set of accumulation points.

Acknowledgements H.H. Bauschke was supported by the Natural Sciences and Engineering Research Council of Canada and the Canada Research Chair Program. D. Noll acknowledges hospitality of the University of British Columbia Kelowna and support by the Pacific Institute of the Mathematical Sciences during the preparation of this paper. The visualization was made possible by GeoGebra [10].


Pierre Apkarian

Dominikus Noll

Optimization-Based Control Design Techniques and Tools

published or not. The documents may come

Optimization-Based Control Design Techniques and Tools

Pierre Apkarian Dominikus Noll

Abstract-Structured output feedback controller synthesis is an exciting new concept in modern control design, which bridges between theory and practice in so far as it allows for the first time to apply sophisticated mathematical design paradigms like H∞or H2-control within control architectures preferred by practitioners. The new approach to structured H∞-control, developed during the past decade, is rooted in a change of paradigm in the synthesis algorithms. Structured design may no longer be based on solving algebraic Riccati equations or matrix inequalities. Instead, optimization-based design techniques are required. In this essay we indicate why structured controller synthesis is central in modern control engineering. We explain why non-smooth optimization techniques are needed to compute structured control laws, and we point to software tools which enable practitioners to use these new tools in high technology applications.

I. MOTIVATIONS

In the modern high technology field control engineers usually face a large variety of concurring design specifications such as noise or gain attenuation in prescribed frequency bands, damping, decoupling, constraints on settling-or risetime, and much else. In addition, as plant models are generally only approximations of the true system dynamics, control laws have to be robust with respect to uncertainty in physical parameters or with regard to un-modeled high frequency phenomena. Not surprisingly, such a plethora of constraints presents a major challenge for controller tuning, due not only to the ever growing number of such constraints, but also because of their very different provenience.

The dramatic increase in plant complexity is exacerbated by the desire that regulators should be as simple as possible, easy to understand and to tune by practitioners, convenient to hardware implement, and generally available at low cost. Such practical constraints explain the limited use of blackbox controllers, and they are the driving force for the implementation of structured control architectures, as well as for the tendency to replace hand tuning methods by rigorous algorithmic optimization tools.

II. STRUCTURED CONTROLLERS

Before addressing specific optimization techniques, we introduce some basic terminology for control design problems with structured controllers. Given a plant P in state-space Pierre Apkarian is with ONERA, 2 Av. Ed. Belin, 31055, Toulouse, France Pierre.Apkarian@onera.fr Dominikus Noll is with Université de Toulouse, Institut de Mathmatiques, 118, route de Narbonne, 31062 Toulouse, France dominikus.noll@math.univ-toulouse.fr form P :

   ẋ = Ax + B 1 w + B 2 u z = C 1 x + D 11 w + D 12 u y = C 2 x + D 21 w + D 22 u (1)
where x ∈ R n P is the state, u ∈ R nu the control, y ∈ R ny the measured output, w ∈ R nw the exogenous input, and z ∈ R nz the regulated output, an output feedback controller K in state-space form

K : ẋK = A K x K + B K y u = C K x K + D K y (2)
with

A K ∈ R n K ×n K , is called structured if the matrices A K , B K , C K , D K depend
smoothly on a design parameter x ∈ R n , referred to as the vector of tunable parameters. Formally, we have differentiable mappings

A K = A K (x), B K = B K (x), C K = C K (x), D K = D K (x),
and we abbreviate these by the notation K = K(x). A structured controller synthesis problem is then an optimization problem of the form minimize T wz (P, K(x)) subject to K(x) closed-loop stabilizing

K = K(x), x ∈ R n (3)
where T wz (P, K) = F (P, K) is the lower feedback connection of (1) with (2) as in Fig. 1 (left), also called the Linear Fractional Transformation [START_REF] Varga | Symbolic and numerical software tools for lft-based low order uncertainty modeling[END_REF]. The norm • stands for the H ∞ -norm, the H 2 -norm, or any other system norm, while the optimization variable x ∈ R n regroups the tunable parameters in the design. Standard examples of structured controllers K(x) include realizable PIDs, observer-based, reduced-order, or decentralized controllers, which in state-space are expressed as:

  0 0 1 0 -1/τ -k D /τ k I 1/τ k P + k D /τ   , A -B 2 K c -K f C 2 K f -K c 0 , A K B K C K D K ,    q diag i=1 A Ki q diag i=1 B Ki q diag i=1 C Ki q diag i=1 D Ki    .
In the case of a PID the tunable parameters are x = (τ, k P , k I , k D ), for observer-based controllers x regroups the estimator and state-feedback gains (K f , K c ), for reduced order controllers n K < n P the tunable parameters x are the n 2 K + n K n y + n K n u + n y n u unknown entries in (A K , B K , C K , D K ), and in the decentralized form x regroups the unknown entries in A K1 , . . . , D Kq . In contrast, full-order controllers have the maximum number N = n 2 P + n P n y + n P n u + n y n u of degrees of freedom and are referred to as unstructured or as black-box controllers. contrast, full-order controllers have the maximum number N = n 2 P + n P n y + n P n u + n y n u of degrees of freedom and are referred to as unstructured or as black-box controllers. More sophisticated controllers structures K(x) arise form architectures like for instance a 2-DOF control arrangement with feedback block K 2 and a set-point filter K 1 as in Fig. 1 (right). Suppose K 1 is the 1st-order filter K 1 (s) = a/(s + a) and K 2 the PI feedback K 2 (s) = k P + k I /s. Then the transfer T ry from r to y can be represented as the feedback connection of P and K(x) with

z w P y K u r K 1 K 2 G y e u + -
P := 2 6 6 4 A 0 0 B C 0 0 D 0 I 0 0 C 0 I D 3 7 7 5 , K(x) :=  K 1 (s) 0 0 K 2 (s) ,
where K(x, s) takes a typical block-diagonal structure featuring the tunable elements x = (a, k P , k I ).

In much the same way arbitrary multi-loop interconnections of fixed-model elements with tunable controller blocks K i (x) can be re-arranged as in Fig. 2, so that K(x) captures all tunable blocks in a decentralized structure general enough to cover most engineering applications. Synthesis of K = block-diag(K 1 , . . . , K N ) against multiple requirements or models P (1) , . . . , P (M ) . Each K i (x) can be structured.

The structure concept is equally useful to deal with the second central challenge in control design: system uncertainty. The latter may be handled with µ-synthesis techniques [START_REF] Stein | Beyond singular values and loopshapes[END_REF] if a parametric uncertain model is available. A less ambitious but often more practical alternative consists in optimizing the structured controller K(x) against a finite set of plants P (1) , . . . , P (M ) representing model variations due to uncertainty, aging, sensor and actuator breakdown, un-modeled dynamics, in tandem with the robustness and performance specifications. This is again formally covered by Fig. 2 and leads to a multi-objective constrained optimization problem of the form

minimize f (x) = max k2SOFT,i2I k kT (k) wizi (K(x))k subject to g(x) = max k2HARD,j2J k kT (k) wj zj (K(x))k  1 K = K(x) closed-loop stabilizing x 2 R n (4)
where

T (k)
wizi denotes the ith closed-loop robustness or performance channel w i ! z i for the k-th plant model P (k) (s). The rationale of ( 4) is to minimize the worst-case cost of the soft constraints kT

III. OPTIMIZATION TECHNIQUES OVER THE YEARS

During the late 1990s the necessity to develop design techniques for structured regulators K(x) was recognized [START_REF] Fares | An Augmented Lagrangian Method for a Class of LMI-Constrained Problems in Robust Control Theory[END_REF], and the limitations of synthesis methods based on algebraic Riccati equations (AREs) or linear matrix inequalities (LMIs) became evident, as these techniques can only provide black-box controllers. Unfortunately, the lack of appropriate synthesis techniques for structured K(x) led to the unsatisfying situation, where sophisticated approaches like the H 1 paradigm developed by academia since the 1980s could not be brought to work for the design of those controller structures K(x) preferred by practitioners. Design engineers had to continue to rely on heuristic and ad-hoc tuning techniques, with only limited scope and reliability. As an example: post-processing to reduce a black-box controller to a practical size is prone to failure. It may at best be considered fill-in for a rigorous design method which directly computes a reduced-order controller. Similarly, hand-tuning of the parameter x remains a puzzling task because of the loop interactions, and fails as soon as complexity increases.

In the late 1990s and early 2000s, a change of methods was observed. Structured H 2 -and H 1 -synthesis problems (3) were addressed by bilinear matrix inequality (BMI) optimization, which used local optimization techniques based on augmented Lagrangian [START_REF]An Augmented Lagrangian Method for a Class of LMI-Constrained Problems in Robust Control Theory[END_REF]- [START_REF] Kocvara | A code for convex nonlinear and semidefinite programming[END_REF], sequential semidefinite programming methods [START_REF] Fares | Robust Control via Sequential Semidefinite Programming[END_REF], [START_REF] Apkarian | A Spectral Quadratic-SDP Method with Applications to Fixed-Order H 2 and H∞ Synthesis[END_REF], and non-smooth methods for BMIs [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex semi-infinite maximum eigenvalue functions[END_REF], [START_REF] Lemaréchal | Nonsmooth algorithms to solve semidefinite programs[END_REF]. However, these techniques were based on the bounded real lemma or similar matrix inequalities, and were therefore of limited success due to the presence of Lyapunov variables, whose number grows quadratically in n P +n K and represents the bottleneck of this approach.

The epoch-making change occurs with the introduction of non-smooth optimization techniques [START_REF] Noll | Spectral bundle methods for nonconvex maximum eigenvalue functions: first-order methods[END_REF]- [START_REF]Nonsmooth optimization for multidisk H∞ synthesis[END_REF] to programs (3) and ( 4). Today non-smooth methods have superseded matrix inequality-based techniques and may be considered the state-of-art as far as realistic applications are concerned. The transition took almost a decade.

Alternative control-related local optimization techniques and heuristics include the gradient sampling technique of [START_REF] Burke | Stabilization via nonsmooth, nonconvex optimization[END_REF], derivative-free optimization discussed in [START_REF] Kolda | Optimization by direct search: new perspectives on some classical and modern methods[END_REF], [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF], particle swarm optimization, see [START_REF] Oi | Development of PSObased PID tuning method[END_REF] and references therein, and also evolutionary computation techniques [START_REF] Lieslehto | PID controller tuning using evolutionary programming[END_REF]. The last three classes do not take advantage of derivative information More sophisticated controllers structures K(x) arise form architectures like for instance a 2-DOF control arrangement feedback block K 2 and a set-point filter K 1 as in Fig. 1 (right). Suppose K 1 is the 1st-order filter K 1 (s) = a/(s + a) and K 2 the PI feedback K 2 (s) = k P + k I /s. Then the transfer T ry from r to y can be represented as the feedback connection of P and K(x) with

P :=     A 0 0 B C 0 0 D 0 I 0 0 -C 0 I -D     , K(x) := K 1 (s) 0 0 K 2 (s) ,
where K(x, s) takes a typical block-diagonal structure featuring the tunable elements x = (a, k P , k I ).

In much the same way arbitrary multi-loop interconnections of fixed-model elements with tunable controller blocks K i (x) can be re-arranged as in Fig. 2, so that K(x) captures all tunable blocks in a decentralized structure general enough to cover most engineering applications.

contrast, full-order controllers have the maximum number N = n 2 P + n P n y + n P n u + n y n u of degrees of freedom and are referred to as unstructured or as black-box controllers. More sophisticated controllers structures K(x) arise form architectures like for instance a 2-DOF control arrangement with feedback block K 2 and a set-point filter K 1 as in Fig. 1 (right). Suppose K 1 is the 1st-order filter K 1 (s) = a/(s + a) and K 2 the PI feedback K 2 (s) = k P + k I /s. Then the transfer T ry from r to y can be represented as the feedback connection of P and K(x) with

P := 2 6 6 4 A 0 0 B C 0 0 D 0 I 0 0 C 0 I D 3 7 7 5 , K(x) :=  K 1 (s) 0 0 K 2 (s) ,
where K(x, s) takes a typical block-diagonal structure featuring the tunable elements x = (a, k P , k I ).

In much the same way arbitrary multi-loop interconnections of fixed-model elements with tunable controller blocks K i (x) can be re-arranged as in Fig. 2, so that K(x) captures all tunable blocks in a decentralized structure general enough to cover most engineering applications.

z 1 w 1 z 2 w 2 . . . . . . P (1) K 1 . . . K N y u . . . z 1 w 1 z 2 w 2 . . . . . . P (M ) K 1 . . . K N y u Fig. 2.
Synthesis of K = block-diag(K 1 , . . . , K N ) against multiple requirements or models P (1) , . . . , P (M ) . Each K i (x) can be structured.

The structure concept is equally useful to deal with the second central challenge in control design: system uncertainty. The latter may be handled with µ-synthesis techniques [START_REF] Stein | Beyond singular values and loopshapes[END_REF] if a parametric uncertain model is available. A less ambitious but often more practical alternative consists in optimizing the structured controller K(x) against a finite set of plants P (1) , . . . , P (M ) representing model variations due to uncertainty, aging, sensor and actuator breakdown, un-modeled dynamics, in tandem with the robustness and performance specifications. This is again formally covered by Fig. 2 and leads to a multi-objective constrained optimization problem of the form

minimize f (x) = max k2SOFT,i2I k kT (k) w i z i (K(x))k subject to g(x) = max k2HARD,j2J k kT (k) w j z j (K(x))k  1 K = K(x) closed-loop stabilizing x 2 R n (4)
where T (k) w i z i denotes the ith closed-loop robustness or performance channel w i ! z i for the k-th plant model P (k) (s). The rationale of ( 4) is to minimize the worst-case cost of the soft constraints kT (k) w i z i k, k 2 SOFT, while enforcing the hard constraints kT (k)

w j z j k  1, k 2 HARD.

III. OPTIMIZATION TECHNIQUES OVER THE YEARS

During the late 1990s the necessity to develop design techniques for structured regulators K(x) was recognized [START_REF] Fares | An Augmented Lagrangian Method for a Class of LMI-Constrained Problems in Robust Control Theory[END_REF], and the limitations of synthesis methods based on algebraic Riccati equations (AREs) or linear matrix inequalities (LMIs) became evident, as these techniques can only provide black-box controllers. Unfortunately, the lack of appropriate synthesis techniques for structured K(x) led to the unsatisfying situation, where sophisticated approaches like the H 1 paradigm developed by academia since the 1980s could not be brought to work for the design of those controller structures K(x) preferred by practitioners. Design engineers had to continue to rely on heuristic and ad-hoc tuning techniques, with only limited scope and reliability. As an example: post-processing to reduce a black-box controller to a practical size is prone to failure. It may at best be considered fill-in for a rigorous design method which directly computes a reduced-order controller. Similarly, hand-tuning of the parameter x remains a puzzling task because of the loop interactions, and fails as soon as complexity increases.

In the late 1990s and early 2000s, a change of methods was observed. Structured H 2 -and H 1 -synthesis problems (3) were addressed by bilinear matrix inequality (BMI) optimization, which used local optimization techniques based on augmented Lagrangian [START_REF]An Augmented Lagrangian Method for a Class of LMI-Constrained Problems in Robust Control Theory[END_REF]- [START_REF] Kocvara | A code for convex nonlinear and semidefinite programming[END_REF], sequential semidefinite programming methods [START_REF] Fares | Robust Control via Sequential Semidefinite Programming[END_REF], [START_REF] Apkarian | A Spectral Quadratic-SDP Method with Applications to Fixed-Order H 2 and H∞ Synthesis[END_REF], and non-smooth methods for BMIs [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex semi-infinite maximum eigenvalue functions[END_REF], [START_REF] Lemaréchal | Nonsmooth algorithms to solve semidefinite programs[END_REF]. However, these techniques were based on the bounded real lemma or similar matrix inequalities, and were therefore of limited success due to the presence of Lyapunov variables, whose number grows quadratically in n P +n K and represents the bottleneck of this approach.

The epoch-making change occurs with the introduction of non-smooth optimization techniques [START_REF] Noll | Spectral bundle methods for nonconvex maximum eigenvalue functions: first-order methods[END_REF]- [START_REF]Nonsmooth optimization for multidisk H∞ synthesis[END_REF] to programs (3) and (4). Today non-smooth methods have superseded matrix inequality-based techniques and may be considered the state-of-art as far as realistic applications are concerned. The transition took almost a decade.

Alternative control-related local optimization techniques and heuristics include the gradient sampling technique of [START_REF] Burke | Stabilization via nonsmooth, nonconvex optimization[END_REF], derivative-free optimization discussed in [START_REF] Kolda | Optimization by direct search: new perspectives on some classical and modern methods[END_REF], [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF], particle swarm optimization, see [START_REF] Oi | Development of PSObased PID tuning method[END_REF] and references therein, and also evolutionary computation techniques [START_REF] Lieslehto | PID controller tuning using evolutionary programming[END_REF]. The last three classes do not take advantage of derivative information Fig. 2.

Synthesis of K = block-diag(K 1 , . . . , K N ) against multiple requirements or models P (1) , . . . , P (M ) . Each K i (x) can be structured.

The structure concept is equally useful to deal with the second central challenge in control design: system uncertainty. The latter may be handled with µ-synthesis techniques [START_REF] Stein | Beyond singular values and loopshapes[END_REF] if a parametric uncertain model is available. A less ambitious but often more practical alternative consists in optimizing the structured controller K(x) against a finite set of plants P (1) , . . . , P (M ) representing model variations due to uncertainty, aging, sensor and actuator breakdown, un-modeled dynamics, in tandem with the robustness and performance specifications. This is again formally covered by Fig. 2 and leads to a multi-objective constrained optimization problem of the form

minimize f (x) = max k∈SOFT,i∈I k T (k) wizi (K(x)) subject to g(x) = max k∈HARD,j∈J k T (k) wj zj (K(x)) ≤ 1 K = K(x) closed-loop stabilizing x ∈ R n (4)
where

T (k)
wizi denotes the ith closed-loop robustness or performance channel w i → z i for the k-th plant model P (k) (s). The rationale of ( 4) is to minimize the worst-case cost of the soft constraints T

III. OPTIMIZATION TECHNIQUES OVER THE YEARS

During the late 1990s the necessity to develop design techniques for structured regulators K(x) was recognized [START_REF] Fares | An Augmented Lagrangian Method for a Class of LMI-Constrained Problems in Robust Control Theory[END_REF], and the limitations of synthesis methods based on algebraic Riccati equations (AREs) or linear matrix inequalities (LMIs) became evident, as these techniques can only provide black-box controllers. Unfortunately, the lack of appropriate synthesis techniques for structured K(x) led to the unsatisfying situation, where sophisticated approaches like the H ∞ paradigm developed by academia since the 1980s could not be brought to work for the design of those controller structures K(x) preferred by practitioners. Design engineers had to continue to rely on heuristic and ad-hoc tuning techniques, with only limited scope and reliability. As an example: post-processing to reduce a black-box controller to a practical size is prone to failure. It may at best be considered fill-in for a rigorous design method which directly computes a reduced-order controller. Similarly, hand-tuning of the parameter x remains a puzzling task because of the loop interactions, and fails as soon as complexity increases.

In the late 1990s and early 2000s, a change of methods was observed. Structured H 2 -and H ∞ -synthesis problems (3) were addressed by bilinear matrix inequality (BMI) optimization, which used local optimization techniques based on augmented Lagrangian [START_REF]An Augmented Lagrangian Method for a Class of LMI-Constrained Problems in Robust Control Theory[END_REF]- [START_REF] Kocvara | A code for convex nonlinear and semidefinite programming[END_REF], sequential semidefinite programming methods [START_REF] Fares | Robust Control via Sequential Semidefinite Programming[END_REF], [START_REF] Apkarian | A Spectral Quadratic-SDP Method with Applications to Fixed-Order H 2 and H∞ Synthesis[END_REF], and non-smooth methods for BMIs [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex semi-infinite maximum eigenvalue functions[END_REF], [START_REF] Lemaréchal | Nonsmooth algorithms to solve semidefinite programs[END_REF]. However, these techniques were based on the bounded real lemma or similar matrix inequalities, and were therefore of limited success due to the presence of Lyapunov variables, whose number grows quadratically in n P +n K and represents the bottleneck of this approach.

The epoch-making change occurs with the introduction of non-smooth optimization techniques [START_REF] Noll | Spectral bundle methods for nonconvex maximum eigenvalue functions: first-order methods[END_REF]- [START_REF]Nonsmooth optimization for multidisk H∞ synthesis[END_REF] to programs (3) and ( 4). Today non-smooth methods have superseded matrix inequality-based techniques and may be considered the state-of-art as far as realistic applications are concerned. The transition took almost a decade.

Alternative control-related local optimization techniques and heuristics include the gradient sampling technique of [START_REF] Burke | Stabilization via nonsmooth, nonconvex optimization[END_REF], derivative-free optimization discussed in [START_REF] Kolda | Optimization by direct search: new perspectives on some classical and modern methods[END_REF], [START_REF] Apkarian | Controller design via nonsmooth multidirectional search[END_REF], particle swarm optimization, see [START_REF] Oi | Development of PSObased PID tuning method[END_REF] and references therein, and also evolutionary computation techniques [START_REF] Lieslehto | PID controller tuning using evolutionary programming[END_REF]. The last three classes do not take advantage of derivative information and rely on function evaluations only. They are therefore applicable to a broad variety of problems including those where function values arise from complex numerical simulations. The combinatorial nature of these techniques, however, limits their use to small problems with a few tens of variable. More significantly, these methods often lack a solid convergence theory. In contrast, as we have demonstrated over recent years, specialized non-smooth techniques are highly efficient in practice, are based on a sophisticated convergence theory, capable of solving medium size problems in a matter of seconds, and are still operational for large size problems with several hundreds of states.

IV. NON-SMOOTH OPTIMIZATION TECHNIQUES

The benefit of the non-smooth casts ( 3) and ( 4) lies in the possibility to avoid searching for Lyapunov variables, a major advantage as their number (n P + n K ) 2 /2 usually largely dominates n, the number of true decision parameters x. Lyapunov variables do still occur implicitly in the function evaluation procedures, but this has no harmful effect for systems up to several hundred states. In abstract terms, a non-smooth optimization program has the form

minimize f (x) subject to g(x) ≤ 0 x ∈ R n (5)
where f, g : R n → R are locally Lipschitz functions and are easily identified from the cast in (4).

In the realm of convex optimization, non-smooth programs are conveniently addressed by so-called bundle methods, introduced in the late 1970s by Lemaréchal [START_REF]An extension of Davidon methods to nondifferentiable problems[END_REF]. Bundle methods are used to solve difficult problems in integer programming or in stochastic optimization via Lagrangian relaxation. Extensions of the bundling technique to nonconvex problems like (3) or (4) were first developed in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex semi-infinite maximum eigenvalue functions[END_REF], [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF]- [START_REF]Nonsmooth optimization for multidisk H∞ synthesis[END_REF], [START_REF] Apkarian | A trust region spectral bundle method for nonconvex eigenvalue optimization[END_REF], and in more abstract form, in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF].

Fig. 3 shows a schematic view of a non-convex bundle method consisting of a descent-step generating inner loop (yellow block) comparable to a line search in smooth optimization, embedded into the outer loop (blue box), where serious iterates are processed, stopping criteria are applied, and the model tradition is assured. At the core of the interaction between inner and outer loop is the management of the proximity control parameter τ , which governs the stepsize x-y k between trial steps y k at the current serious iterate x. Similar to the management of a trust region radius or of the stepsize in a linesearch, proximity control allows to do shorter trial steps if agreement of the local model with the true objective function is poor, and to allow larger steps if agreement is satisfactory.

Oracle-based bundle methods traditionally assure global convergence in the sense of subsequences under the sole hypothesis that for every trial point x the function value f (x) and a Clarke subgradient φ ∈ ∂f (x) are provided. In automatic control applications it is as a rule possible to provide more specific information, which may be explored to speed up the convergence.

Flow chart of proximity control bundle algorithm

In the of the H 2 -norm, f (x) = T wz (P, K(x)) 2 , computing function value and gradient requires essentially the solution of four Lyapunov equations of size n P +n K , see [START_REF] Apkarian | Mixed H 2 /H∞ control via nonsmooth optimization[END_REF], [START_REF] Rautert | Computational design of optimal output feedback controllers[END_REF]. For the H ∞ -norm, f (x) = T wz (P, K(x)) ∞ , function evaluation is based on the Hamiltonian algorithm of [START_REF] Benner | L∞-norm computation for continuous-time descriptor systems using structured matrix pencils[END_REF], [START_REF] Boyd | A bisection method for computing the H∞ norm of a transfer matrix and related problems[END_REF]. The Hamiltonian matrix is of size n P + n K , so that function evaluations may be costly for very large plant state dimension (n P > 500), even though the number of outer loop iterations of the bundle algorithm is not affected a large n P and generally relates to n, the dimension x. The additional cost for subgradient computation for large n P is relatively cheap as it relies on linear algebra [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF].

V. COMPUTATIONAL TOOLS

The novel non-smooth optimization methods became available to the engineering community since 2010 via the MATLAB Robust Control Toolbox [START_REF]Robust Control Toolbox 4.2[END_REF], [START_REF] Gahinet | Structured H∞ synthesis in MATLAB[END_REF]. Routines HINFSTRUCT , LOOPTUNE and SYSTUNE are versatile enough to define and combine tunable blocks K i (x), to build and aggregate design requirements T (k) wz of different nature, and to provide suitable validation tools. Their implementation was carried out in cooperation with P. Gahinet (MathWorks). These routines further exploit the structure of the problem (4) to enhance efficiency, see [START_REF]Nonsmooth optimization for multiband frequency domain control design[END_REF] and [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF] .

It should be mentioned that design problems with multiple hard constraints are inherently complex. It is well known that even simultaneous stabilization of several plants P (j) with a structured control law K(x) is generally NP-complete, so that exhaustive methods are expected to fail even for small to medium problems. The principled decision made in [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF], and reflected in the MATLAB routines, is to rely on local optimization techniques instead. This leads to weaker certificates, but has the advantage to work successfully in practice. In the same vein, in (4) it is preferable to rely on a mixture of soft and hard requirements, for instance, by the use of exact penalty functions [START_REF] Noll | Spectral bundle methods for nonconvex maximum eigenvalue functions: first-order methods[END_REF]. Key features implemented in the mentioned MATLAB routines are discussed in [START_REF]Nonsmooth optimization for multiband frequency domain control design[END_REF], [START_REF] Gahinet | Structured H∞ synthesis in MATLAB[END_REF], [START_REF] Apkarian | Tuning controllers against multiple design requirements[END_REF].

VI. DESIGN EXAMPLE

Design of a feedback regulator is an interactive process, in which tools like SYSTUNE , LOOPTUNE or HINFSTRUCT support the designer in various ways. In this section we illustrate their enormous potential by solving a multi-model, fixed-structure reliable flight control design problem.

In reliable flight control one has to maintain stability and adequate performance not only in nominal operation, but also in various scenarios where the aircraft undergoes outages in elevator and aileron actuators. In particular, wind gusts must be alleviated in all outage scenarios to maintain safety. Variants of this problem are addressed in [START_REF] Liao | Reliable robust flight tracking control: An LMI approach[END_REF]. In addition to nominal operation, we consider 8 outage scenarios shown in Table I. The design requirements are as follows:

• Good tracking performance in µ, α, and β with adequate decoupling of the three axes.

• Adequate rejection of wind gusts of 5 m/s.

• Maintain stability and acceptable performance in the face of actuator outage.

Tracking is addressed by an LQG-cost [START_REF] Maciejowski | Multivariable Feedback Design[END_REF], which penal-izes integrated tracking error e and control effort u via

J = lim T →∞ E 1 T T 0 W e e 2 + W u u 2 dt . (6)
Diagonal weights W e and W u provide tuning knobs for trade-off between responsiveness, control effort, and balancing of the three channels. We use W e = diag [START_REF]An extension of Davidon methods to nondifferentiable problems[END_REF][START_REF] Liao | Reliable robust flight tracking control: An LMI approach[END_REF][START_REF]An extension of Davidon methods to nondifferentiable problems[END_REF], W u = I 3 for normal operation and W e = diag(8, 12, 8), W u = I 3 for outage conditions. Modeldependent weights allow to express the fact that nominal operation prevails over failure cases. Weights for failure cases are used to achieve limited deterioration of performance or of gust alleviation under deflection surface breakdown.

The second requirement, wind gust alleviation, is treated as a hard constraint limiting the variance of the error signal e in response to white noise w g driving the Dryden wind gust model. In particular, the variance of e is limited to 0.01 for normal operation and to 0.03 for the outage scenarios.

With the notation of section IV, the functions f (x) and g(x) in ( 5 WeightedVariance options in [START_REF]Robust Control Toolbox 4.2[END_REF].

With this setup, we tuned the controller gains K i and K x for the nominal scenario only (nominal design) and for all 9 scenarios (fault-tolerant design). The responses to setpoint changes in µ, α, and β with a gust speed of 5 m/s are shown in Fig. 5 for the nominal design and in Fig. 6 for the faulttolerant design. As expected, nominal responses are good but notably deteriorate when faced with outages. In contrast, the fault-tolerant controller maintains acceptable performance in outage situations. Optimal performance (square root of LQG cost J in ( 6)) for the fault-tolerant design is only slightly worse than for the nominal design (26 vs. 23). The nonsmooth program (5) was solved with SYSTUNE and the faulttolerant design (9 models, 11 states, 27 parameters) took 30 seconds on Mac OS X with 2.66 GHz Intel Core i7 and 8 GB RAM. The reader is referred to [START_REF]Robust Control Toolbox 4.2[END_REF] or higher versions, further examples, and additional details. CROSS REFERENCES

Fig. 1 .

 1 Fig. 1. Black-box full-order controller K on the left, structured 2-DOF control architecture with K = block-diag(K 1 , K 2 ) on the right.

Fig. 2 .

 2 Fig. 2.Synthesis of K = block-diag(K 1 , . . . , K N ) against multiple requirements or models P(1) , . . . , P(M ) . Each K i (x) can be structured.

  wizi k, k 2 SOFT, while enforcing the hard constraints kT (k) wj zj k  1, k 2 HARD.

Fig. 1 .

 1 Fig. 1. Black-box full-order controller K on the left, structured 2-DOF control architecture with K = block-diag(K 1 , K 2 ) on the right.

Fig. 1 .

 1 Fig. 1. Black-box full-order controller K on the left, structured 2-DOF control architecture with K = block-diag(K 1 , K 2 ) on the right.

  wizi , k ∈ SOFT, while enforcing the hard constraints T (k)wj zj ≤ 1, k ∈ HARD.

Fig. 4 .

 4 Fig. 4. Synthesis interconnection for reliable control

,

  ) are f (x) := max k=1,...,9 T (k) rz (x) 2 and g(x) := max k=1,...,9 T (k) wge (x) 2 , where r denotes the setpoint inputs in µ, α and β. The regulated output z isz T := (W with x = (vec(K i ), vec(K x )) ∈ R27 . Soft constraints are the square roots of J in (6) with appropriate weightings W e and W u , hard constraints the RMS values of e, suitably weighted to reflect variance bounds of 0.01 and 0.03. These requirements are covered by the Variance and

Fig. 5 .

 5 Fig. 5. Responses to step changes in µ, α and β for nominal design.

Fig. 6 .

 6 Fig. 6. Responses to step changes in µ, α and β for fault-tolerant design.