Minh Ngoc Dao

Dominikus Noll

Pierre Apkarian

ROBUST EIGENSTRUCTURE CLUSTERING BY NONSMOOTH OPTIMIZATION

Keywords: Structured feedback control, eigenstructure assignment, modal shaping, nonsmooth optimization, frequency-domain, robust design

We extend classical eigenstructure assignment to more realistic problems where additional performance and robustness specifications arise. Our aim is to combine time-domain constraints, as reflected by pole location and eigenvector structure, with frequency-domain objectives such as the H 2 , H ∞ or Hankel norms. Using pole clustering, we allow poles to move in polydisks of prescribed size around their nominal values, driven by optimization. Eigenelements, that is poles and eigenvectors, are allowed to move simultaneously and serve as decision variables in a specialized nonsmooth optimization technique. Two aerospace applications illustrate the power of the new method.

Introduction

Since its introduction by Wonham [START_REF] Wonham | On pole assignment in multi-input controllable linear systems[END_REF] and Moore [START_REF] Moore | On the flexibility offered by state feedback in multivariable systems beyond closed loop eigenvalue assignment[END_REF], eigenstructure assignment has developed into a powerful controller design tool in the aerospace sector and in other high technology fields. Eigenstructure assignment aims at shaping the responses of the closedloop system to certain input signals by way of two mechanisms. The placement of closedloop modes to stabilize and achieve satisfactory transients, and eigenvector structure to decouple responses to specific initial conditions. In this paper we are concerned with the design of output feedback control laws, where only partial eigenstructure assignment or pole placement can be expected. In that case the standard approach to first selecting a partial set of closed-loop modes λ 1 , . . . , λ p , and then using the remaining degrees of freedom to shape the corresponding closed-loop eigenvectors (v i , w i ), is prone to failure to stabilize the system, as the remaining closed-loop modes cannot be influenced directly.

As a remedy we propose to assign the eigenelements (λ i , v i , w i ) simultaneously. We allow eigenelements (λ i , v i , w i ) to move in the neighborhood of their nominal values (λ 0 i , v 0 i , w 0 i ) in such a way that closed-loop stability and performance can be further improved. The price for this gain of flexibility is that eigenelement assignment can no longer be achieved by linear algebra methods alone. Instead, a combination of nonlinear optimization and linear algebra is required.

Over the years there have already been attempts to enhance eigenspace control using off-the-shelf optimization. An early approach is Sobel and Shapiro [START_REF] Sobel | A robust pitch pointing control law[END_REF], where hand-tuning of eigenvalues was shown to improve stability margins of the controlled system. In [START_REF]Flight control examples of robust eigenstructure assignment[END_REF] the same authors elaborate on this idea and suggest a first-order gradient method. In [START_REF] Apkarian | Structured stability robustness improvement by eigenspace techniques: a hybrid methodology[END_REF][START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF], a sequential quadratic programming (SQP) technique with finite-difference gradients was used to improve µ robustness indicators, with eigenvalues and some eigenvectors as decision variables. In [START_REF] Patton | Robust control design via eigenstructure assignment, genetic algorithms and gradient-based optimisation[END_REF], Patton and Liu make full use of the freedom offered by eigenstructure assignment to improve the frequency-domain sensitivities functions S and KS. They use a genetic algorithm in tandem with gradient-based techniques. The same idea is applied to a variety of problems in their monograph [START_REF] Liu | Eigenstructure Assignment for Control System Design[END_REF]. In the same vein, reference [START_REF] Kshatriya | Optimized partial eigenstructure assignment-based design of a combined PSS and active damping controller for a DFIG[END_REF] exploits the Nelder-Mead direct search method to optimize assignable eigenvalues and eigenvectors, while safeguarding stability of unassigned eigenvalues via constraints. In [START_REF] Magni | A toolbox for robust modal control design (RMCT)[END_REF], eigenstructure assignment with dynamic compensators and linear programming (LP) or quadratic programming (QP) are used to achieve stability and performance for an entire family of plants. Merits of these approaches have been demonstrated in numerous applications. See [START_REF] Liu | Eigenstructure Assignment for Control System Design[END_REF] and references therein.

In this work, we suggest a novel approach to eigenstructure assignment based on a nonsmooth optimization technique, which has the following features:

• Unassigned poles are constrained to be stable, which secures stability of the closedloop system. • Additional performance or robustness requirements such as H 2 or H ∞ are handled rigorously by accounting for their nonsmoothness. Nonsmoothness arises due to the spectral abscissa, and via H ∞ -norm or Hankel norm based requirements, but also when max-function of differentiable functions such as the H 2 -norm are built. The key observation is that disregarding nonsmoothness is a serious source of numerical trouble. Avoiding this pitfall is a central motivation of this work. Our investigation leads to a theoretically justified nonsmooth method with local convergence certificate, which has good performance in practical applications. The focus of this paper is on control aspects. A thorough convergence analysis of the proposed algorithm is given in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF][START_REF] Noll | Cutting plane oracles to minimize non-smooth non-convex functions[END_REF][START_REF] Gabarrou | Design of a flight control architecture using a non-convex bundle method[END_REF][START_REF] Dao | Minimizing memory effects of a system[END_REF] for the interested readers.

The structure of the paper is as follows. Section 2 recalls the basics of eigenstructure assignment using static output feedback and its variation as pole clustering, where poles are allowed to move in small polydisks around their nominal values. Section 3 extends the pole clustering problem to a variety of performance or robustness criteria and gives a pseudo-code of our algorithmic approach to those problems. Overdetermined and underdetermined eigenproblems are discussed in Section 4. Section 5 shows how subgradients are computed for typical design requirements. Our nonsmooth solver, along with its convergence properties, is presented in Section 6. Sections 7 and 8 illustrate our approach. We design a launcher and an aircraft control system, two cases where poles and eigenvector structure play an important role.

Partial eigenstructure assignment

Consider a linear time-invariant system described by the equations

(1) ẋ = Ax + Bu y = Cx with x ∈ R n , u ∈ R m and y ∈ R p . Given a self-conjugate set Λ = {λ 1 , . . . , λ p } ⊂ C -,
partial pole placement consists in computing a static output feedback control law u = Ky for (1) such that λ 1 , . . . , λ p become eigenvalues of the closed-loop system

ẋ = (A + BKC)x.
As is well-known [START_REF] Moore | On the flexibility offered by state feedback in multivariable systems beyond closed loop eigenvalue assignment[END_REF], solving the set of linear equations

A -λ i I n B v i w i = 0, with v i ∈ C n , w i ∈ C m , i = 1, . . . , p leads to a (static) control law (2) K = [w 1 , . . . , w p ] (C [v 1 , . . . , v p ]) -1 ∈ R m×p
with the desired closed-loop modes, provided the v i are chosen in such a way that the p×p matrix C [v 1 , . . . , v p ] is invertible, i.e., if span{v 1 , . . . , v p } ∩ ker(C) = {0}. Note that the outlined technique is readily specialized to state-feedback C = I and extended to nonzero feedthrough D = 0 and to dynamic compensators through a preliminary augmentation of the plant [START_REF] Sobel | Eigenstructure assignment[END_REF].

In the case m > 1, it is possible to achieve more. One may then additionally shape the v i , or w i , e.g. by arranging v ij = 0 or w ik = 0 for certain j, k. Formally this can be expressed by linear equations

(3) A -λ i I n B M i N i v i w i = 0 r i , with suitable M i ∈ C m i ×n , N i ∈ C m i ×m , r i ∈ C m i , m i 0, i = 1, . . . , p, leaving at least one degree of freedom in each triplet (λ i , v i , w i ) ∈ C 1+n+m
. This is usually referred to as partial eigenstructure assignment. Typical choices of M i , N i , r i can be found in our experimental Sections 7 and 8.

The traditional approach to eigenstructure assignment consists in first choosing the set Λ ⊂ C -, then introducing the desired structural constraints on the eigenvectors v i , w i via the matrices M i , N i and the vector r i , using the remaining degrees of freedom, and then computing v i , w i accordingly. Unfortunately, fixing the λ i may be too restrictive, because partial eigenvalue placement does not guarantee stability in closed-loop, so that some post-processing based on trial-and-error is often required. Greater flexibility in the design is achieved by moving (λ i , v i , w i ) simultaneously.

What we have in mind is to interpret the eigenstructure equations (3) as mathematical programming constraints and then optimize closed-loop stability subject to these constraints. With the definition α(A) := max{Re λ : λ eigenvalue of A} of the spectral abscissa, this leads us to an optimization program of the form (4) minimize α(A + BKC)

subject to A -λ i I n B M i N i v i w i = 0 r i for i = 1, . . . , p | Re λ i -Re λ 0 i | δ i , | Im λ i -Im λ 0 i | δ i , i = 1, . . . , p K = W (CV ) -1 as in (2).
Here the λ 0 i ∈ C -are nominal closed-loop poles, and the δ i are tolerances which allow the poles to move around their nominal values. As soon as K with α(A + BKC) < 0 is reached, the optimization of (4) can be stopped with an internally stabilizing solution of the partial eigenstructure assignment procedure.

Including performance criteria

While ( 4) is a natural approach to optimize closed-loop stability in partial eigenstructure assignment, it seems even more attractive to include also closed-loop performance or robustness criteria into the set-up. Given a linear time-invariant plant P in standard form [START_REF] Apkarian | Structured stability robustness improvement by eigenspace techniques: a hybrid methodology[END_REF] P :

   ẋ = Ax + B 1 w + Bu z = C 1 x + D 11 w + D 12 u y = Cx + D 21 w
where x ∈ R n is the state vector, u ∈ R m the vector of control inputs, w ∈ R m 1 the vector of exogenous inputs, y ∈ R p the vector of measurements and z ∈ R p 1 the controlled or performance vector, let u = Ky be a static output feedback control law for [START_REF] Apkarian | Structured stability robustness improvement by eigenspace techniques: a hybrid methodology[END_REF]. Then the closed-loop performance channel w → z has the state-space representation

T w→z (K) : ẋ = (A + BKC)x + (B 1 + BKD 21 )w z = (C 1 + D 12 KC)x + (D 11 + D 12 KD 21 )w.
Note the slight abuse of notation in [START_REF] Apkarian | Structured stability robustness improvement by eigenspace techniques: a hybrid methodology[END_REF] because the state-space data of P may include filters, weightings or other dynamic elements that are not present in [START_REF] Wonham | On pole assignment in multi-input controllable linear systems[END_REF]. We assume the distinction will be clear from the context.

Given a self-conjugate eigenvalue set Λ 0 = {λ 0 1 , . . . , λ 0 p } ⊂ C -and tolerances δ i , we now consider the following extension of (4):

(6) minimize T w→z (K) subject to A -λ i I n B M i N i v i w i = 0 r i for i = 1, . . . , p | Re λ i -Re λ 0 i | δ i , | Im λ i -Im λ 0 i | δ i , i = 1, . . . , p K = K(λ, v, w) as in (2)
where λ 0 i are nominal closed-loop pole positions, and (3) again conveys additional structural constraints on v, w. As compared to (4), the cost function T w→z (K) in ( 6) may now be used to enhance stability and to achieve additional performance or robustness specifications of the design.

Standard choices of • include the H ∞ -norm • ∞ , the H 2 -norm • 2 , or the Hankel norm • H . One generally expects that T w→z (K) < ∞ implies closed-loop stability, but should this fail, it is possible to add a stability constraint c(λ, v, w) = α(A+BKC)+ε 0 to the cast [START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF], where ε > 0 is some small threshold. Altogether we propose the following Algorithm 1. Optimized partial eigenstructure assignment

Input: Nominal modal set Λ 0 = {λ 0 1 , . . . , λ 0 p } with distinct λ 0 i . Output: Optimal modal set Λ = {λ 1 , . . . , λ p }, v i , w i , K * .
Step 1 (Nominal assignment). Perform standard eigenstructure assignment based on Λ 0 and structural constraints

M i , N i , r i . Obtain nominal eigenvectors v 0 i , w 0 i , i = 1, . . . , p. Assure that C[v 0 1 , . . . , v 0 p ] is invertible and obtain nominal K 0 = W 0 (CV 0 ) -1 .
Step 2 (Stability and performance). If K 0 assures closed-loop stability and good performance T w→z (K 0 ) , stop the algorithm. Otherwise, goto step 3.

Step 3 (Tolerances). Allow tolerances

| Re λ i -Re λ 0 i | δ i , | Im λ i -Im λ 0 i | δ i , i = 1, . . . , p.
Step 4 (Parametric clustering). Solve the optimization program (6) using a nonsmooth descent algorithm with (λ 0 , v 0 , w 0 ) as initial seed.

Step 5 (Synthesis). Return optimal Λ = {λ 1 , . . . , λ p }, v, w, and K * .

Structure of eigenproblems

In this section we discuss practical ways to deal with the general nonlinear constraint (3) in [START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF]. We assume that (A, B) is controllable, which is equivalent to [A -λI n B] having full row rank n for all λ in C (see, e.g., [START_REF] Zhou | Robust and Optimal Control[END_REF]Theorem 3.1]). To deal with (3), we observe that the m i 's can be distinct and the possibility m i = 0 is not excluded. We now distinguish two cases.

The first case is when m i m. Here pole assignment is ensured by pre-solving for v i in (3). We get

v i = (λ i I -A) -1 Bw i .
In this case eigenvector decoupling is only possible in the least-square sense by minimizing the Euclidean norm of

M i v i + N i w i -r i . Upon defining the transfer function F i (λ) := M i (λI -A) -1 B + N i ,
and assuming for simplicity that F i (λ) has full-column rank for λ in the neighborhood of the nominal λ 0 , we have

w i = F i (λ i ) † r i ,
where F i (λ i ) † denotes the Moore-Penrose inverse or left-inverse of F i at λ i . Altogether we have derived the expression ( 7)

v i w i = (λ i I -A) -1 B I F i (λ i ) † r i .
Vectors v i and w i are now defined explicitly as functions of λ i . It follows that a parametrization of the control law (2) in the sense of structured synthesis introduced in [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF] has been obtained. Tunable variables in this parametrization are the desired assignable eigenvalues Λ = {λ 1 , . . . , λ p }.

The rationale in this first case is as follows. We want to gain some flexibility in the assignment by allowing λ i to move in a neighborhood of the nominal λ 0 i . Now if the (v 0 i , w 0 i ) are computed from [START_REF] Patton | Robust control design via eigenstructure assignment, genetic algorithms and gradient-based optimisation[END_REF] for the nominal value λ 0 i , the (v i , w i ), depending continuously on λ i via [START_REF] Patton | Robust control design via eigenstructure assignment, genetic algorithms and gradient-based optimisation[END_REF], will move in a neighborhood of the nominal (v 0 i , w 0 i ), so that optimization may decrease the cost function and thereby enhance stability and performance. The outlined approach therefore generalizes eigenstructure assignment with approximate decoupling as discussed in [START_REF] Sobel | Eigenstructure assignment[END_REF].

If F i (λ) is not guaranteed to have full-column rank in the neighborhood λ 0 , the cast in ( 6) could be modified as follows:

(8) minimize max T w→z (K) , µ max i=1,...,p F i (λ i )w i -r i 2 subject to | Re λ i -Re λ 0 i | δ i , | Im λ i -Im λ 0 i | δ i , i = 1, . . . , p. K = [w 1 , . . . , w p ] (C [v 1 , . . . , v p ]) -1 K closed-loop stabilizing
where µ is a penalty parameter used to weigh the relative importance of robustness or performance as expressed through T w→z (K) against eigenvector shaping. Here the objective becomes a max-function which is truly nonsmooth and thus requires special handling.

The second case is when m i < m. Here we partition

B = [B i Q i ], N i = [P i R i ], w i = u i t i , such that B i , P i have m i columns and u i ∈ C m i . Then (3) becomes A -λ i I n B i M i P i v i u i = 0 r i - Q i R i t i .
Assuming that the matrix

A i (λ i ) = A -λ i I n B i M i P i ∈ C (n+m i )×(n+m i )
is invertible in a neighborhood of the nominal λ 0 i , we get the parametrization

v i = v i (λ i , t i ), u i = u i (λ i , t i ), which in explicit form is (9) v i u i = A i (λ i ) -1 -Q i t i r i -R i t i .
The idea is now the same as in the first case. Allow λ i to move around their nominal values λ 0 i , and t i ∈ C m-m i around the nominal t 0 i . That also allows the dependent variables v i , u i to move in a neighborhood of their nominal values v 0 i , u 0 i , and optimization uses this to enhance stability and robustness. In this second case we have enough degrees of freedom to achieve true decoupling of some of the channels by satisfying M i v i + N i w i = r i exactly.

In order to apply nonlinear and nonsmooth optimization techniques to programs of the form [START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF] it is necessary to provide derivative information at acceptable cost. As we shall see, this may be implemented by simple linear algebra techniques. We have the following propositions with proofs given in the Appendix.

Proposition 1 (Over-specified eigenstructure).

Let K = W (CV ) -1 with W = [w 1 . . . w p ] and V = [v 1 . . . v p ]. If m i m then (10) ∂K ∂λ i = 0 • • • ∂w i ∂λ i -KC ∂v i ∂λ i • • • 0 (CV ) -1 ,
where v i , w i are given in [START_REF] Patton | Robust control design via eigenstructure assignment, genetic algorithms and gradient-based optimisation[END_REF] and

(11)    ∂v i ∂λ i ∂w i ∂λ i    = (λ i I -A) -1 BF i (λ i ) † M i -I F i (λ i ) † M i (λ i I -A) -2 BF i (λ i ) † r i .
Proof. See Appendix.

Proposition 2 (Under-specified eigenstructure).

Let K = W (CV ) -1 with W = [w 1 . . . w p ] and V = [v 1 . . . v p ]. Suppose m i < m, partitioning w i = u i t i with u i ∈ C m i and t i = [t 1i , . . . , t (m-m i )i ] ∈ C m-m i , then ∂K ∂λ i = 0 • • • ∂u i ∂λ i 0 -KC ∂v i ∂λ i • • • 0 (CV ) -1 , ∂K ∂t ki = 0 • • • ∂u i ∂t ki e ki -KC ∂v i ∂t ki • • • 0 (CV ) -1 , (12)
where e ki ∈ R m-m i is the vector all of whose components are zero, except the kth component which is one, and

   ∂v i ∂λ i ∂v i ∂t ki ∂u i ∂λ i ∂u i ∂t ki    = I n 0 n×m i 0 m i ×n I m i A i (λ i ) -1 v i 0 -s ik , with s ik the kth column of Q i R i . (13)
Proof. See Appendix.

Remark 1. As derivatives have to be evaluated repeatedly in minimization programs, it is desirable to pre-calculate as many elements as possible in [START_REF] Magni | A toolbox for robust modal control design (RMCT)[END_REF] and [START_REF] Noll | Cutting plane oracles to minimize non-smooth non-convex functions[END_REF]. This is what we discuss next. Substantial speed-up can be achieved in the under-specified case m i < m since A i (λ i ) is a reduced rank modification of a constant matrix, that is, not depending on λ i . We therefore pre-compute

A B i M i P i -1 = P i 11 P i 12 P i 21 P i 22 ,
where P i 11 and P i 22 are of size n × n and m i × m i , respectively. Using the Sherman-Woodbury-Morrison formula [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] for

A i (λ i ) = A B i M i P i + -I 0 (λ i I) I 0 gives A i (λ i ) -1 = (I n -λ i P i 11 ) -1 P i 11 (I n -λ i P i 11 ) -1 P i 12 P i 21 (I n -λ i P i 11 ) -1 P i 22 + P i 21 λ i (I n -λ i P i 11 ) -1 P i 12 .
As a consequence, there is only need to compute the inverse of the smaller matrix (I nλ i P i 11 ) to get the entries in (13). Remark 2. Our algorithm can be extended to include nonlinear constraints on v i . We just add those to program [START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF]. Note also that the algorithm will return the standard nominal modal set λ 0 = {λ 0 1 , . . . , λ 0 p } if we choose δ i = 0, i = 1, . . . , p, so we present a genuine extension of the traditional assignment procedure.

System norms and their subdifferential in closed-loop

To solve program ( 6) algorithmically, we have to compute function values and subgradients of the cost function f (x) := T w→z (K(x)) 2 , where • is the H ∞ -norm • ∞ , the H 2 -norm • 2 or the Hankel norm • H , and where x represents the decision variables. Here x regroups λ i if m i m, and (λ i , t i ) if m i < m, i = 1, . . . , p. The gradients given in [START_REF] Magni | A toolbox for robust modal control design (RMCT)[END_REF], respectively [START_REF] Noll | Cutting plane oracles to minimize non-smooth non-convex functions[END_REF], are generally complex gradients. Algorithmic implementation requires passing from complex to real gradients. This is done using Wirtinger formulas [19, Section 2.3]. For a complex variable z, we have that

∂K/∂ Re z = ∂K/∂z + ∂K/∂ z = 2 Re(∂K/∂z), ∂K/∂ Im z = (∂K/∂z -∂K/∂ z) = -2 Im(∂K/∂z).
For simplicity of the notation, it is assumed from now on that x is a real q-dimensional vector regrouping real and imaginary parts of all free parameters (λ i , t i ). Partial derivatives with respect to x will be denoted K i (x) := ∂K(x)/∂x i in the sequel of the paper. In consequence it now remains to compute Clarke subgradients of T w→z (K) 2 with respect to K. By the generalized chain rule [START_REF] Clarke | Generalized gradients of Lipschitz functionals[END_REF], this requires subgradients of the norm in question, and the derivative of the transfer function T w→z (K) with respect to K.

Concerning the closed-loop, and to prepare the following, by setting

A cl = A + BKC, B cl = B 1 + BKD 21 , C cl = C 1 + D 12 KC, D cl = D 11 + D 12 KD 21 ,
the controllability Gramian X and the observability Gramian Y can be obtained from the Lyapunov equations [START_REF] Zhou | Robust and Optimal Control[END_REF]

A cl X + XA cl + B cl B cl = 0, (14)
A cl Y + Y A cl + C cl C cl = 0. (15) 5.1. The H ∞ -norm. Consider a stable LTI system G : ẋ = Ax + Bw z = Cx + Dw with state x ∈ R n , input w ∈ R m , and output z ∈ R p . It is well-known that the H ∞ -norm of G is defined as G ∞ = sup ω∈R σ max (G(jω)) = sup ω∈R λ max (G(jω) H G(jω)),
where σ max denotes the maximum singular value of a matrix, and λ max denotes the maximum eigenvalue of a matrix. We now replace G by T w→z (K) and rewrite

f (K) = T w→z (K) 2 ∞ = sup ω∈R f (K, ω),
with f (K, ω) := λ max T w→z (K, ω) H T w→z (K, ω) . Using the notation

T w→z (K, s) G 12 (K, s) G 21 (K, s) * = C cl C (sI -A cl ) -1 B cl B + D cl D 12 D 21 * ,
and following [21, Lemma 1], closed-loop stability implies that either f (K) = f (K, ω) for all ω or f (K) = f (K, ω) for a finite number of active frequencies ω 1 , . . . , ω q . From [17,

Section IV] we now obtain the Clarke subgradients of f at K as

Φ U = 2 q k=1 Re G 21 (K, ω k )T w→z (K, ω k ) H R k U k R H k G 12 (K, ω k ) ,
where R k is a matrix whose columns form an orthonormal basis of the eigenspace of dimension r k ∈ N associated with λ max T w→z (K, ω k ) H T w→z (K, ω k ) , and where

U k ∈ S r k , U k 0, q k=1 Tr(U k ) = 1.
The symbol S m stands for the space of m × m symmetric or Hermitian matrices, and Tr(M ) denotes the trace of M . By the application of the chain rule in [START_REF] Clarke | Generalized gradients of Lipschitz functionals[END_REF], we deduce that the Clarke subdifferential of f at x is the set

∂f (x) = Tr(K 1 (x) Φ U ), . . . , Tr(K q (x) Φ U ) : Φ U ∈ ∂f (K) . 5.2. The H 2 -norm. The H 2 -norm of a system G of the form (16) G : ẋ = Ax + Bw z = Cx
is defined as

G 2 = 1 2π +∞ -∞ Tr(G(jω) H G(jω))dω 1/2
. Suppose D cl does not explicitly depend on K, which is e.g. the case for D 12 = 0 or D 21 = 0. Then it is reasonable to assess the closed-loop system via the H 2 -norm of (A cl , B cl , C cl , 0). We have 14) and [START_REF] Sobel | Eigenstructure assignment[END_REF], it follows from [22, Theorem 3.2] that f is differentiable at each closed-loop stabilizing K, and

f (K) = T w→z (K) 2 2 = Tr(B cl Y B cl ) = Tr(C cl XC cl ). Using (
∇f (K) = 2 B Y + D 12 C cl XC + 2B Y B cl D 21 .
Therefore, ∇f (x) = Tr(K 1 (x) ∇f (K)), . . . , Tr(K q (x) ∇f (K)) for all x for which K(x) is closed-loop stabilizing.

5.3. The Hankel norm. For a stable system G of the form (16), we think of w(t) as an excitation at the input which acts over the time period 0 t T . Then the ring of the system G after the excitation has stopped at time T is z(t) for t > T . If signals are measured in the energy norm, this leads to the Hankel norm of G defined as

G H = sup T >0 ∞ T z zdt 1/2 : z = Gw, T 0 w wdt 1, w(t) = 0 for t > T .
The Hankel norm [START_REF] Glover | All optimal Hankel-norm approximations of linear multivariable systems and their L ∞error bounds[END_REF][START_REF] Dao | Minimizing memory effects of a system[END_REF] can be understood as measuring the tendency of a system to store energy, which is later retrieved to produce undesired noise effects known as system ring. Minimizing the Hankel norm T w→z (K) H therefore reduces ringing in the closedloop channel w → z.

If we assume as above that D cl does not explicitly depend on K, it is reasonable to assess the channel w → z via the objective

f (K) = T w→z (K) 2 H = λ max (XY )
, where X and Y are the closed-loop Gramians ( 14) and [START_REF] Sobel | Eigenstructure assignment[END_REF]; see also [START_REF] Dao | Minimizing memory effects of a system[END_REF]Lemma 1]. Due to positive semidefiniteness of B cl B cl and C cl C cl , closed-loop stability assures positive semidefiniteness of X and Y in ( 14) and [START_REF] Sobel | Eigenstructure assignment[END_REF]. Therefore, although the product XY need not be symmetric, we have

λ max (XY ) = λ max (X 1 2 Y X 1 2 ) = λ max (Y 1 2 XY 1 2
), which brings us back to the realm of eigenvalue theory for symmetric matrices. Let

Z := X 1 2 Y X 1 2
and take R to be a matrix whose columns form an orthonormal basis of the eigenspace of Z of dimension r ∈ N associated with λ max (Z). We write M i (x) := ∂M (x)/∂x i as before, and M

1 2 i short for (M 1
2 ) i , i = 1, . . . , q. Then according to [14, Proposition 1], the Clarke subdifferential of f at x is

∂f (x) = {(Tr(RU R Z 1 (x)), . . . , Tr(RU R Z q (x))) : U ∈ S r , U 0, Tr(U ) = 1}, with (17) Z i (x) = X 1 2 i (x)Y X 1 2 + X 1 2 Y i (x)X 1 2 + X 1 2 Y X 1 2 i (x).
Here X i (x), Y i (x) and X 1 2

i (x) are the solutions of the following Lyapunov equations

A cl X i (x) + X i (x)A cl = -BK i (x)CX -X(BK i (x)C) (18)
-BK i (x)D 21 B cl -B cl (BK i (x)D 21 ) , A cl Y i (x) + Y i (x)A cl = -(BK i (x)C) Y -Y BK i (x)C (19) -(D 12 K i (x)C) C cl -C cl D 12 K i (x)C, X 1 2 X 1 2 i (x) + X 1 2 i (x)X 1 2 = X i (x). ( 20
)

Nonsmooth solver

Step 4 of our main Algorithm 1 requires a subroutine to solve [START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF]. Here we use a nonsmooth descent algorithm, presented as Algorithm 2, which we now discuss briefly. To extend the scope, we consider a constrained optimization programs of the more abstract form [START_REF] Bompart | Second-order nonsmooth optimization for H ∞ synthesis[END_REF] minimize f (x) subject to h(x) 0 Ax b where x ∈ R q is the decision variable, and f and h are potentially nonsmooth and nonconvex. This covers program [START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF], where f (x) = T w→z (K(x)) 2 for one of the norms discussed in Section 5, while h(x) 0 could represent the stability constraint [START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF] can easily be converted to the form Ax b. According to the cases discussed in Section 4, the decision variable x regroups either the λ i with [START_REF] Patton | Robust control design via eigenstructure assignment, genetic algorithms and gradient-based optimisation[END_REF], or (λ i , t i ) as in [START_REF] Kshatriya | Optimized partial eigenstructure assignment-based design of a combined PSS and active damping controller for a DFIG[END_REF]. The cast ( 8) is also covered by [START_REF] Bompart | Second-order nonsmooth optimization for H ∞ synthesis[END_REF].

α(A + BKC) + ≤ 0. The polydisk constraints | Re λ i -Re λ 0 i | δ i , | Im λ i -Im λ 0 i | δ i in
To solve [START_REF] Bompart | Second-order nonsmooth optimization for H ∞ synthesis[END_REF] we use a progress function at the current iterate x,

F (•, x) = max{f (•) -f (x) -νh(x) + , h(•) -h(x) + },
for some fixed parameter ν > 0, which is successively minimized subject to the linear constraints. Antecedents of this idea can for instance be found in Polak [24, Section 2.2.2] in the smooth case, or Polak and Wardi [START_REF] Polak | Nondifferentiable optimization algorithm for designing control systems having singular value inequalities[END_REF] in a nonsmooth setting, and in our own contributions [START_REF] Apkarian | Mixed H 2 /H ∞ control via nonsmooth optimization[END_REF][START_REF] Gabarrou | Design of a flight control architecture using a non-convex bundle method[END_REF][START_REF] Dao | Minimizing memory effects of a system[END_REF], where more details and convergence proofs can be found.

Algorithm 2. Nonsmooth optimization subroutine

Parameters: 0 < γ < γ < 1, 0 < γ < Γ < 1, 0 < q < ∞, q < T < ∞.
Step 1 (Initialize outer loop). Choose initial iterate x 1 with Ax 1 b and matrix

Q 1 = Q 1 with -qI Q 1 qI. Initialize memory control parameter τ 1 > 0 such that Q 1 + τ 1 I 0. Put outer loop counter j = 1.
Step 2 (Stopping test). At outer loop counter j, stop if x j is a KKT-point or a critical point of constraint violation. Otherwise, goto inner loop.

Step 3 (Initialize inner loop). Put inner loop counter k = 1 and initialize τ 1 = τ j . Build initial working model

Φ 1 (•, x j ) = g 0j (• -x j ) + 1 2 (• -x j ) Q j (• -x j ), where g 0j ∈ ∂ 1 F (x j , x j ).
Step 4 (Trial step generation). At inner loop counter k find solution y k of the tangent program minimize Φ k (y, x j ) + τ k 2 yx j 2 subject to Ay b, y ∈ R n .

Step 5 (Acceptance test). If

ρ k = F (y k , x j ) Φ k (y k , x j ) γ,
put x j+1 = y k (serious step), quit inner loop and goto step 8. Otherwise (null step), continue inner loop with step 6.

Step 6 (Update working model). Generate a cutting plane m k (•, x j ) = a k + g k (•x j ) at null step y k and counter k using downshifted tangents. Compute aggregate plane m * k (•, x j ) = a * k +g * k (•-x j ) at y k , and then build new working model Φ k+1 (•, x j ) by including cutting plane and aggregate plane.

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

ρ k = Φ k+1 (y k , x j ) Φ k (y k , x j ) and put τ k+1 = τ k if ρ k < γ, 2τ k if ρ k γ.
Increase inner loop counter k and loop back to step 4.

Step 8 (Update Q j and memory element). Update matrix Q j → Q j+1 respecting Q j+1 = Q j+1 and -qI Q j+1 qI. Then store new memory element

τ j+1 = τ k if ρ k < Γ, 1 2 τ k if ρ k Γ. Increase τ j+1 if necessary to ensure Q j+1 +τ j+1 I 0. If τ j+1 > T then re-set τ j+1 = T .
Increase outer loop counter j and loop back to step 2.

Convergence theory of Algorithm 2 is discussed in [START_REF] Gabarrou | Design of a flight control architecture using a non-convex bundle method[END_REF][START_REF] Dao | Minimizing memory effects of a system[END_REF]. The following result is slightly more general than the main convergence theorem in [START_REF] Gabarrou | Design of a flight control architecture using a non-convex bundle method[END_REF] or [START_REF] Dao | Minimizing memory effects of a system[END_REF], but can be obtained based essentially on the same convergence analysis: Theorem 1. Suppose f and h in program (21) are lower-C 1 functions in the sense of [START_REF] Spingarn | Submonotone subdifferentials of Lipschitz functions[END_REF] such that the following conditions hold:

(i) f is weakly coercive on the constraint set Ω = {x ∈ R q : h(x) 0, Ax b}, i.e., if x j ∈ Ω and x j → ∞, then f (x j ) is not monotonically decreasing. (ii) h is weakly coercive on P = {x ∈ R q : Ax b}, i.e., if x j ∈ P and x j → ∞, then h(x j ) is not monotonically decreasing.

Then the sequence of serious iterates x j ∈ P generated by Algorithm 2 is bounded, and every accumulation point x * of the x j satisfies x * ∈ P and 0 ∈ ∂ 1 F (x * , x * ) + A η * for some multiplier η * 0 with η * (Ax *b) = 0. In other words, x * is either a critical point of constraint violation, or a Karush-Kuhn-Tucker point of program [START_REF] Bompart | Second-order nonsmooth optimization for H ∞ synthesis[END_REF].

Note that the functions f, h used in ( 6) are indeed lower-C 1 functions, see [START_REF] Gabarrou | Design of a flight control architecture using a non-convex bundle method[END_REF][START_REF] Dao | Minimizing memory effects of a system[END_REF], so our convergence theory applies. Convergence for even larger classes of nonsmooth functions is discussed in [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF][START_REF] Noll | Cutting plane oracles to minimize non-smooth non-convex functions[END_REF]. For additional insight into this type of nonconvex bundle method see [START_REF] Apkarian | Mixed H 2 /H ∞ control via nonsmooth optimization[END_REF][START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF][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].

Control of a launcher in atmospheric flight

We consider attitude control of a satellite launcher in atmospheric flight. The linear model

ẋ = Ax + Bu y = Cx is specified as A =            Z w Z q + U 0 Z θ Z v 0 Z ψ Z p Z φ M w M q 0 0 M r 0 M p 0 0 T q 0 0 T r 0 0 0 Y w 0 Y θ Y v Y r Y ψ Y p Y φ 0 N q 0 N v N r 0 N p 0 0 P q 0 0 P r 0 0 0 0 L q 0 0 L r 0 L p 0 0 F q 0 0 F r 0 1 0            , B =   Z βz M βz 0 0 0 0 0 0 0 0 0 Y βy N βy 0 0 0 0 0 0 0 0 0 L βr 0   .
The states and controls are defined in Tables 1 and2, while the vector of measurements is y = [q θ r ψ p φ] ∈ R 6 . The model has been obtained from linearization of the nonlinear equations [START_REF] Mclean | Automatic Flight Control Systems[END_REF] about a steady state flight point U 0 = 88.11 m/s, v 0 = 0.678 m/s, w 0 = -1.965 m/s, p 0 = -0.0006 rad/s, q 0 = 0.0026 rad/s, r 0 = 0.0046 rad/s, θ 0 = 8.38 • , ψ 0 = 3.48 • , φ 0 = 11.99 • , the procedure being explained in [START_REF] Greensite | Elements of Modern Control Theory[END_REF]. Numerical data in A, B are gathered in Table 3. • Decoupling of the 3 axes (θ, q), (ψ, r), and (φ, p).

Z w -0.0162 M w 0.0022 Y w -6e-4 N q 5e-4 Z q 87.9 -88.11 M q 0.0148 Y θ -2.11 N v -M w Z θ -9.48 M r -0.0005 Y v Z w N r 0.0151 Z v 0.0006 M p 0.0042 Y r -87.9 N p -0.0024 Z ψ -2.013 T q 0.
• Well-damped responses to set-points in θ, ψ, and φ, the selector outputs.

• Settling times around 2.5 seconds.

We use a set-point tracking control architecture with MIMO PI feedback as shown in Figure 1. Tunable matrix gains are therefore K P and K I .

Tracking performance is incorporated into program (6) by minimizing the tracking error transfer function T w ref →e (K). Pole placement with integral action is easily formulated using the augmented state-space matrices The control law is structured conformably upon defining

A a = A 0 -HC 0 , B a = B 0 , C a = C 0
W = [w 1 . . . w 9 ] , V = [v 1 . . . v 9 ] , [A a -λ i I 11 |B a ] v i w i = 0. K a = W (C a V ) -1 , K a = [-K P K I ] .

Study 1.

In a first study we compare traditional and optimized partial pole placement without shaping of eigenvectors. We start by choosing reference values ξ, ω to achieve appropriate second-order system responses. We have chosen the desired damping

ξ = √ 2
2 , and natural frequencies ω 1 = 2.1, ω 2 = 2.2, ω 3 = 1.8, which leads to the nominal modal set Λ = {λ 0 1 , . . . , λ 0 9 }, with

λ 0 1 = -ω 1 ξ +  1 -ξ 2 , λ 0 2 = -ω 1 ξ - 1 -ξ 2 , λ 0 3 = -ω 2 ξ +  1 -ξ 2 , λ 0 4 = -ω 2 ξ - 1 -ξ 2 , λ 0 5 = -ω 3 ξ +  1 -ξ 2 , λ 0 6 = -ω 3 ξ - 1 -ξ 2 , λ 0
7 = -3.5, λ 0 8 = -4, λ 0 9 = -4.5. Classical pole placement now leads to the initial controller K 0 in Algorithm 1. To find the optimal controller K * , we follow Algorithm 1 and minimize the tracking error w ref → e subject to the pole placement constraint in (6) via Algorithm 2, which returns the optimal controller K * . We have run program [START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF] with three different norms, the Hankel norm, the H ∞ -, and the H 2 -norm. The improvements in the cost function can be seen in Table 4. The wandering of the poles during optimization shown in Figure 3 corresponds to the case of the Hankel norm. Figure 2 shows that decoupling is substantially improved in all three cases. Note the sluggish responses for the initial controller are due to unassigned modes of classical eigenstructure assignment. This is in contrast with the proposed approach in which modes that are left unspecified are indirectly assigned to achieve additional performance requirements.

Study 2.

In our second study we compare standard and optimized eigenstructure assignment. We achieve preliminary decoupling of the modes by choosing structural constraints on eigenvectors v i . These constraints comply with decoupling requirements of the launcher motion. The eigenvectors v 1 and its complex conjugate v 2 are chosen to have zero entries in the rows corresponding to ψ and φ. The eigenvector v 3 and its conjugate These structural constraints define the matrices M i , N i of (3) in each case. We have again tested the Hankel, H ∞ and H 2 -norms in the objective f of (6).

The optimal controller K * computed by Algorithm 1 for the Hankel norm gives the value T (P perf , K * ) H = 0.7360, while the initial controller K 0 leads to T (P perf , K 0 ) H = 0.7787. Similar improvements are obtained for the other norms. The step responses are shown in Figure 4.

In conclusion, the launcher application shows that decoupling can be significantly enhanced through optimization even without shaping of the v i (study 1) if the performance channel T w ref →e is used within optimization program [START_REF] Morris | Synthesizing robust mode shapes with µ and implicit model following[END_REF]. The second study shows that even when 0's are assigned to specific v ik 's, the use of optimization is still useful, as it significantly enhances decoupling as demonstrated by simulation.

Application to autopilot design for a civil aircraft

In this section, we consider the longitudinal dynamics for the robust civil aircraft model (RCAM) at a nominal condition with the aircraft in its standard configuration: aircraft air speed of 80 m/s, aircraft altitude of 305 m (1000 ft), aircraft mass of 120 tons, aircraft centre of gravity at 23% horizontal MAC and 0% vertical MAC, flight path angle of 0 • Step Response Time (seconds) where states are described in Table 5, the input vector is u = [δ t δ th ] with δ t the tailplane deflection and δ th the throttle position. The vector of measurements is y = [q n z w V z V c ] , where n z is vertical acceleration, w V vertical velocity, and V c the air the state corresponding to the first order engine model speed. Data borrowed from [START_REF]Robust flight. Control design challenge problem formulation and manual: the research civil aircraft model (RCAM)[END_REF] are given as The controller structure of the longitudinal autopilot with tunable gains K I , K P is similar to the launcher structure given in Figure 1. The output is now y = [q n z w V z] T , and the selector produces y H = [z V c ] T . We next design a closed-loop controller such that altitude is decoupled from air speed command and conversely. This leads to decoupling altitude and altitude-tracking modes from forward speed u B , and decoupling of the air speed track mode from the upwards velocity w B . Other modes are also decoupled from some states to reduce the mutual influence of the aircraft variables. Accordingly, we take the nominal modes as follows: which defines the data M i , N i and r i in (3). The optimal controller K * computed by Algorithm 1 gives T (P perf , K * ) H = 0.6270, while the initial controller K 0 obtained by standard assignment had T (P perf , K 0 ) H = 1.5041. The closed-loop eigenvalues returned by the algorithm are λ 1,2 = -0.8 ± 0.8, λ 3,4 = -0.35 ± 0.05, λ 5 = -0.3, λ 6 = -0.05, λ 7 = -0.37, which shows that some of the poles took indeed the opportunity to wander away from their nominal values once they were allowed to do so. Step responses are compared in Figure 5. The interpretation of the results is that optimization is useful to further enhance decoupling even when eigenvectors are already shaped.

A =          -0.9825 0

Conclusion

We have presented a new approach to partial eigenstructure assignment in output feedback control in which the eigenelements (λ, v, w) are allowed to move simultaneously in a neighborhood of their nominal values (λ 0 , v 0 , w 0 ) obtained by standard partial assignment. The flexibility gained in allowing this is apparent on two fronts. First, stability of unassigned modes is guaranteed, leading to an internally stable closed-loop system. Secondly, criteria such as H ∞ , H 2 and Hankel norms can be incorporated into our formulation to improve performance and/or robustness of the controlled system. The efficiency of the new approach was demonstrated on two aerospace applications, control of a launcher in atmospheric flight, and attitude control of a civil aircraft.

Figure 1 .

 1 Figure 1. Launcher control architecture with MIMO PI-controller

Figure 2 .

 2 Figure 2. Control of a launcher, study 1. Initial and final controllers obtained respectively by standard and optimized eigenstructure assignment in the case where eigenvectors are not structured (m i = 0). Decoupling is improved for each norm

Figure 3 .

 3 Figure 3. Control of launcher, study 1. Itineraries of closed-loop poles in optimized eigenstructure assignment based on the Hankel program (6)

Figure 4 . 1 (

 41 Figure 4. Control of launcher, study 2. Initial and final controller obtained respectively by standard and optimized eigenstructure assignment based on Hankel program with m i = m or m i = m -1

  1927 -9.7758 -0.0325 0.0743 0 0.1836 19.6200 77.3571 -0.7674 -0.2265 -0.6683 0 -6

  The corresponding desired eigenvectors are shaped as v 1,2 = [ * * 0 * * * * * * ] , v 3,4 = [ * * * 0 * * * * * ] , v 5 = [ * * 0 * * * * * * ] , v 6 = [ * * * 0 * * * * * ] , v 7 = [ * * 0 * * * * * * ] ,

Figure 5 .

 5 Figure 5. Aircraft attitude control. Responses to a step command in altitude and in air speed. Optimal controller computed by optimized eigenstructure assignment (m i = 1, m = 2) reduces coupling

Table 1 .

 1 States definitions

		name meaning
		w	vertical velocity (m/s)
		q	pitch rate (deg/s)
		θ	pitch angle (deg)
		v	lateral velocity (m/s)
		r	yaw rate (deg/s)
		ψ	yaw angle (deg)
		p	roll rate (deg/s)
		φ	roll angle (deg)
		Table 2. Controls definitions
	name meaning
	β z	deflection of pitch nozzle actuator (deg)
	β y	deflection of yaw nozzle actuator (deg)
	β r	deflection of roll nozzle actuator (deg)

Table 3 .

 3 Numerical coefficients at steady state flight point

Table 4 .

 4 Launcher study 1. Cost for initial K 0 and optimal K * controllers

	Hankel	H ∞	H 2
	K 0 66.7208 2.3714 45.3537
	K * 0.7135 1.4058 3.0845

Table 5 .

 5 States of the longitudinal model

	name meaning
	q	pitch rate
	θ	pitch angle
	u B	forward speed
	w B	upwards velocity
	z	altitude
	x t	the state corresponding to the first order tailplane model
	x th

Appendix

Proof of Propositions 1 and 2. Let us start by discussing the case m i m. Derivatives of w i with respect to λ i can be derived from the normal equations F i (λ i ) H F i (λ i )w i = F i (λ i ) H r i or directly from the expression of w i in [START_REF] Patton | Robust control design via eigenstructure assignment, genetic algorithms and gradient-based optimisation[END_REF]. Assuming F i (λ i ) is full-column rank, we rewrite

. The partial derivative of w i with respect to λ i is then readily derived by exploiting that for an invertible matrix M depending smoothly on a parameter t, the derivative of its inverse is obtained as

The derivative of v i is obtained in much the same way using the upper part of [START_REF] Patton | Robust control design via eigenstructure assignment, genetic algorithms and gradient-based optimisation[END_REF]. Finally, collecting the results for w i and v i leads to expression [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF]. This allows us now to express the terms ∂K/∂λ i where K = W (CV ) -1 . Using again the derivative of a matrix inverse, we have

Combining with

yields [START_REF] Magni | A toolbox for robust modal control design (RMCT)[END_REF].

Next consider the under-specified case m i < m. We have that (9) yields (13) analogously to the over-specified case. Finally, formulas for ∂K/∂λ i and ∂K/∂t ki in [START_REF] Noll | Cutting plane oracles to minimize non-smooth non-convex functions[END_REF] are obtained from the fact that K = W (CV ) -1 with V = [v 1 . . . v p ] and


Adrien Kaladji
email: kaladrien@hotmail.fr

Alexandre Villena

Remy Pascot

Florent Lalys

Anne Daoudal

Elodie Clochard

Antoine Lucas

Alain Cardon

Fusion imaging for EVAR with mobile c-arm

Keywords: Fusion imaging, mobile c-arm, flat panel, EVAR, hybrid room, registration

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Fusion imaging is a technique that facilitates endovascular navigation but is only available in hybrid rooms. The goal of this study was to evaluate the feasibility of fusion imaging with a mobile C-arm in a conventional operating room through the use of an angio-navigation station.

Methods

From May 2016 to June 2017, the study included all patients who underwent an aortic stent graft procedure in a conventional operating room with a mobile flat-panel detector (Cios Alpha, Siemens) connected to an angio-navigation station (EndoNaut, Therenva). The intention was to perform preoperative 3D CT/perioperative 2D fluoroscopy fusion imaging using an automatic registration process. Registration was considered successful when the software was able to correctly overlay preoperative 3D vascular structures onto the fluoroscopy image. For EVAR, contrast dose, operation and fluoroscopy time were compared to those of a control group drawn from the department's database who underwent a procedure with a C-arm image intensifier.

Results

The study included 54 patients and the procedures performed were: 49 EVAR, 2 TEVAR, 2 IBD, 1 FEVAR. Of the 178 registrations that were initialised, it was possible to use the fusion imaging in 170 cases, i.e. a 95.5% success rate. In the EVAR comparison, there were no difference with the control group (n=103) for fluoroscopy time (21.9 ± 12 vs. 19.5 ± 13 min, p=0.27), but less contrast agent was used in the group undergoing a procedure with the angionavigation station (42.3 ± 22 ml vs. 81.2 ± 48 ml, p<0.001) and operation time was shorter (114 ± 44 vs. 140.8 ± 38 min, p<0.0001).

Conclusion

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Introduction

Fusion imaging is a technique used in vascular surgery in which a preoperative 3D CT scan is projected onto 2D fluoroscopy. This makes it possible to navigate the vascular tree without necessarily having to use iodinated contrast agent. An increasing number of studies are reporting on the use of fusion imaging in aortic endovascular procedures and its benefits [START_REF] Dias | The effects of combining fusion imaging, low-frequency pulsed fluoroscopy, and lowconcentration contrast agent during endovascular aneurysm repair[END_REF][START_REF] Tacher | Image guidance for endovascular repair of complex aortic aneurysms: comparison of twodimensional and three-dimensional angiography and image fusion[END_REF][START_REF] Hertault | Impact of hybrid rooms with image fusion on radiation exposure during endovascular aortic repair[END_REF][START_REF] De Ruiter | Meta-analysis of Cumulative Radiation Duration and Dose During EVAR Using Mobile, Fixed, or Fixed/3D[END_REF][START_REF] Maurel | Techniques to reduce radiation and contrast volume during EVAR[END_REF]. This imaging modality is however based on a certain number of prerequisites that limit the number of centres using it. These prerequisites are currently only fulfilled by floor-mounted fixed imaging systems that make automatic adjustments for movements of the operating table.

Fusion imaging involves alignment (registration) of the preoperative CT scan with the fluoroscopy image. Once registration has been performed, whenever the table or C-arm moves the system is capable of realigning the preoperative CT scan with the 2D fluoroscopy image. In practice, only surgeons performing procedures in hybrid rooms have access to fusion imaging and other modern navigation tools. The cost of hybrid rooms is a major barrier to adoption of this technique.

The FUTUR study (whose French acronym stands for 'feasibility of computer-assisted aortic and iliac endovascular procedures with mobile C-arm') set out to overcome the challenge of performing fusion imaging with a mobile C-arm in a conventional operating room and making it compatible with the clinical workflow.

Patients and Methods

The protocol and informed consent form were approved by the local institutional review board, and all subjects gave informed consent.

Study design

The study was a single-centre, prospective, consecutive feasibility pilot study. The primary

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 5
objective was to evaluate the feasibility of fusion imaging during aortic endovascular procedures with a mobile C-arm through the use of the EndoNaut® angio-navigation station (Therenva; Rennes, France). Secondary objectives were to evaluate the efficiency of the system when deploying infrarenal aortic stent grafts to treat unruptured atheromatous aneurysms (EVAR).

The primary endpoint was the feasibility rate of fusion defined as the convergence of bone registration and the projection of the vascular 3D structure onto the fluoroscopy image (number of initialised registrations/number of convergent registrations). If the accuracy and robustness of the registration algorithm was already quantified in Duménil et al. [START_REF] Dumenil | A versatile intensitybased 3D/2D rigid registration compatible with mobile C-arm for endovascular treatment of abdominal aortic aneurysm[END_REF], the convergence of the algorithm was qualitatively performed by visual inspection of the main angio-navigation system user solely. Registration failures were easily detectable even from a non-expert eye as the algorithm gave an aberrant solution to the registration problem.

Secondary endpoints concerned radiation dose as measured by fluoroscopy time (FT in min), dose-area product (DAP in Gy.cm 2 ) and air kerma (AK in mGy).

Inclusion criteria

Patients eligible for endovascular treatment of aneurysm disease of the aorta.

Procedure performed in a conventional operating room equipped with a mobile flat-panel detector (30×30 cm) (Cios-Alpha, Siemens®, Munich, Germany) and a floating table.

Patients who received written and verbal information about the protocol and did not object to participating in the trial.

Non-inclusion criteria

-Patients who also required a conventional surgical revascularisation procedure or who required an endovascular revascularisation procedure in another site. -Non-analysable CT angiogram (no or poor injection).

-Procedure performed in a hybrid room or in an operating room not equipped with a mobile flat-panel detector.

Fusion imaging principle and operation of EndoNaut® station

The EndoNaut station is connected to the C-arm and retrieves the video signal generated by the C-arm (Fig. 1). The station is positioned in front of the surgeon and becomes his or her primary navigation interface. Sizing and planning data are derived from the preoperative analysis carried out with EndoSize® software (Therenva; Rennes, France) and transmitted to the station by importing a dedicated file from a USB storage device. The following data are exported from the sizing software: aortoiliac mesh with key points used for sizing represented as rings (below the lowest renal artery, iliac bifurcations -Fig. 1). All measurements and aortic 3D screenshots (reporting c-arm angulation) are also exported. Human-computer interaction is via a touch-sensitive tablet in one touch mode. Preoperative CT-scan is searchable during all the procedure thanks to the tablet which controls every "image action" (registration, navigation, measurements, numerical zoom...). To use fusion imaging functionality, the registration must be initialised after setting C-arm/operating table angulation as determined during preprocedural planning. The user then has to align, in an approximate manner and in only one view, the bone 3D volume of the preoperative CT scan with the 2D bony structures of the fluoroscopy image. Next, the software performs perfect, precise alignment of bony structures using rigid 3D/2D registration. The duration of this geometric transformation depends on the size of the image matrix and the C-arm used. For the C-arm used in the present study, the average duration was 15 +/-3s and never exceeded 22 seconds.

Each registration is valid until either the table or the C-arm position/angulation (i.e. C-arm

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

pose) changes. A new registration becomes mandatory in that cases, and on the contrary to fixed C-arm in hybrid room the fusion mask cannot be used to automatically position the table or the gantry. After completion of registration, the vascular tree is projected with the planning data defined using EndoSize®. Rings are visible in the planned deployment areas in the upper and lower landing zones (Fig. 2A). For tortuous anatomy, in order to anticipate anatomical deformations [START_REF] Maurel | Evaluation of visceral artery displacement by endograft delivery system insertion[END_REF][START_REF] Kaladji | Prediction of deformations during endovascular aortic aneurysm repair using finite element simulation[END_REF][START_REF] Kauffmann | Source of errors and accuracy of a two-dimensional/three-dimensional fusion road map for endovascular aneurysm repair of abdominal aortic aneurysm[END_REF] caused by the extrastiff guidewire and the delivery system, a previously simulated deformed 3D model is projected (Fig. 2B). Several studies have specifically addressed simulation [START_REF] Gindre | Estimation of clinically relevant indicators for EVAR using patient-specific finite element simulation[END_REF][START_REF] Gindre | Finite element simulation of the insertion of guidewires during an EVAR procedure: example of a complex patient case, a first step toward patient-specific parameterized models[END_REF][START_REF] Gindre | Patient-Specific Finite-Element Simulation of the Insertion of Guidewire During an EVAR Procedure: Guidewire Position Prediction Validation on 28 Cases[END_REF]. These simulations are not carried out by the station but on a workstation dedicated to sizing and simulation. For TEVAR (thoracic stentgrafts), FEVAR (fenestrated stentgrafts) and IBD procedures (iliac branched device), rings are projected onto the planned landing zones and the ostia of target vessels (Fig. 3 and4). The precision of the fusion in terms of location of target vessels (renal, internal iliac arteries) is always verified by injection of a small volume of iodinated contrast agent. If there is a misalignment of the projection of the 3D fusion mask and angiography due to the insertion of stiff guidewire, adjustments are made.

Control group

For the EVAR procedures, contrast dose, fluoroscopy an operation time were compared to those of a control group of patients who received an aorto-bi-iliac stent graft for non-ruptured AAA using a mobile image intensifier system (OEC 9800, General Electric; GE, USA) without an angionavigation station. These patients were drawn from the local EVAR database and were operated upon between January 2013 and December 2014, the date on which the mobile flat-panel detector in the study came into use.

Statistical analysis

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 8
Quantitative data are expressed as mean ± standard deviation and qualitative data as a number and corresponding percentage. Data were compared using Student's t-test. Significance was set at p<0.05. Every percutaneous access was ultrasound-guided. Patients were deemed unfit for percutaneous access (n=4) when the majority of the anterior wall of the common femoral arteries were calcified.

M A N U S C R I P T A C C E P T E D

Results

Patients and procedures

Primary endpoint

The mean number of registrations initialised per patient was 3.3 ± 1.1. Of the 178 registrations initialised across all patients, it was possible to project the vascular tree onto fluoroscopy in 170 cases, i.e. a feasibility rate of 95.5%. The feasibility rate seems to be similar for all procedure types (Table 1).

Secondary endpoints and comparison with control group

Table 1 shows the radiation and contrast dose for each procedure type. The control group consisted of 103 patients who underwent EVAR (Fig. 4). Demographic data are compared

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 10
between the EVAR group of the FUTUR study and the control group in the table 2. For radiation data, no significant difference between the groups was found. On the basis of comparable body mass indices, fluoroscopy time (21.9 ± 12 vs. 19.5 ± 13 min, p=0.27) and dose-area product (70.6 ± 48 vs. 67.3 ± 74 Gy.cm 2 , p=0.77) did not vary between the two groups. However, less contrast agent was used in the "FUTUR" group (42.3 ± 22 ml vs. 81.2 ± 48 ml, p<0.0001).

Discussion

Herein we report what is to our knowledge the first use of software capable of performing 3D monomodal (CT)/2D fusion with a mobile C-arm for the placement of an aortic stent graft. In all series reporting the use of fusion imaging a fixed system has been necessary [START_REF] Dias | The effects of combining fusion imaging, low-frequency pulsed fluoroscopy, and lowconcentration contrast agent during endovascular aneurysm repair[END_REF][START_REF] Hertault | Impact of hybrid rooms with image fusion on radiation exposure during endovascular aortic repair[END_REF][START_REF] Sailer | CTA with fluoroscopy image fusion guidance in endovascular complex aortic aneurysm repair[END_REF][START_REF] Kaladji | Endovascular aortic repair of a postdissecting thoracoabdominal aneurysm using intraoperative fusion imaging[END_REF][START_REF] Koutouzi | EVAR Guided by 3D Image Fusion and CO2 DSA: A New Imaging Combination for Patients With Renal Insufficiency[END_REF][START_REF] Tacher | Feasibility of Three-Dimensional MR Angiography Image Fusion Guidance for Endovascular Abdominal Aortic Aneurysm Repair[END_REF][START_REF] Alomran | Image fusion for hybrid repair of dislocated superior mesenteric branch of a branched endovascular aortic graft[END_REF][START_REF] Kaladji | Interest of fusion imaging and modern navigation tools with hybrid rooms in endovascular aortic procedures[END_REF][START_REF] Koutouzi | Orthogonal Rings, Fiducial Markers, and Overlay Accuracy When Image Fusion is Used for EVAR Guidance[END_REF][START_REF] Kobeiter | Zero-contrast thoracic endovascular aortic repair using image fusion[END_REF].

The main reason is that the operating table has to be connected to the imaging system so that the latter can realign the preoperative CT scan with the fluoroscopy image whenever the table moves. In this regard, the system we report is not capable of realignment and this may constitute in itself a limitation. However, in routine practice, fusion imaging is only useful at very specific moments during placement of, for example, a bifurcated stent graft. During deployment of the main body, the angle of the C-arm has already been determined in advance, and the practitioner can focus on gradual deployment of the stents and needs to see the position of the renal arteries to perfectly position the stent graft. At this stage, there are no table or C-arm movements; on the contrary, good visualisation of the stent graft and aorta is required. This is also true during deployment of the iliac limbs, in order to ensure precise placement in relation to the internal iliac artery. During insertion of the guidewire, catheters or stent graft, and even during catheterisation of the contralateral stump, table movements are frequent, but projection of the vascular tree is not essential to the practitioner. So although the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

11 software reported here requires that registration be repeated after every table movement, in our opinion this does not constitute a limitation in routine practice.

The technological challenge presented by fusion imaging with a mobile C-arm is to offer registration that is as fast as possible and above all only requires one 2D fluoroscopy view so that the system is compatible with the clinical workflow. The EndoNaut® station is the only software that uses only one 2D view to perform registration with the 3D CT scan, as has been previously reported [START_REF] Dumenil | A versatile intensitybased 3D/2D rigid registration compatible with mobile C-arm for endovascular treatment of abdominal aortic aneurysm[END_REF]. If the current system only uses information extracted from the preoperative 3D CT scan, other types of imaging modalities (MRI, non-enhanced CT, ultrasound) could also be considered and is a subject of further improvement. The first publications for fusion imaging described systems that could only perform 3D/3D registration, which required CBCT at the start of the procedure [START_REF] Kobeiter | Zero-contrast thoracic endovascular aortic repair using image fusion[END_REF][START_REF] Dijkstra | Intraoperative C-arm cone-beam computed tomography in fenestrated/branched aortic endografting[END_REF][START_REF] Kaladji | Safety and accuracy of endovascular aneurysm repair without pre-operative and intra-operative contrast agent[END_REF]. Given the radiation emitted by CBCT, systems evolved towards 3D/2D registration(3) that restricted CBCT to post-procedure assessment [START_REF] Biasi | Intraoperative DynaCT detection and immediate correction of a type Ia endoleak following endovascular repair of abdominal aortic aneurysm[END_REF][START_REF] Biasi | Intra-operative DynaCT improves technical success of endovascular repair of abdominal aortic aneurysms[END_REF][START_REF] Hertault | Benefits of Completion 3D Angiography Associated with Contrast Enhanced Ultrasound to Assess Technical Success after EVAR[END_REF][START_REF] Tornqvist | Intra-operative cone beam computed tomography can help avoid reinterventions and reduce CT follow up after infrarenal EVAR[END_REF]. Hence 2D fluoroscopy acquisition was necessary but 2 views were needed to achieve sufficient precision and this is no longer the case in the present study. With regard to the duration of the registration calculation, we did not report it in a quantitative manner because short duration is a prerequisite for studies such as the present one designed to demonstrate feasibility before envisaging studies to demonstrate clinical benefits. The reported results demonstrate that fusion imaging is feasible in the vast majority of EVAR cases. It was not feasible in cases with C-arm/operating table procedural angulations that were extreme for a conventional operating room. Beyond 35-40° inclination (regardless of orientation), the fluoroscopy image is contaminated by bony structures such as the upper limbs which are not visible on the CT scan and which jeopardize alignment of the two reference images. Nevertheless, this problem can be overcome by abducting and externally rotating the arms during the procedure, for example.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 12

The majority of patients in our study underwent EVAR, as this is the predominant activity and it is also the procedure with the highest reproducibility. The other stent graft cases we reported are anecdotal and were mentioned in order to demonstrate that the station also works with other procedures. Similarly, we did not report cases of iliac recanalization or renal and mesenteric angioplasty, which are now systematically performed with the station.

Unsurprisingly, the radiation dose was not shown to be lower than in the control group. This finding must be interpreted with the utmost caution because the comparison involves completely different systems. The flat-panel detector does not use the same technology as an image intensifier, nor does it have the same physical characteristics, both factors that affect the delivered dose. Hence, drawing conclusions from this finding is fraught with bias.

Moreover, we are reporting the experience of one teaching hospital, one in which experienced senior practitioners, senior practitioners who have become independent more recently and junior practitioners under supervision participate in procedures. There is heterogeneity in the use of X-rays beyond adherence to ALARA principles. For example, a junior practitioner will attempt to track progress of the device through the aorta fluoroscopically, whereas a senior practitioner will insert the device "blindly". The real comparison for radiation dose and X-ray use would be a single-operator randomised series to overcome the anatomical differences between patients that give rise to procedural difficulties. In comparison to image intensifiers, flat-panel detectors offer higher image quality (more pixels per mm 2 ), and the image matrix is larger, so a higher radiation dose could have been expected, but this was not the case.

On the other hand, when we consider use of contrast agent, there is a clear difference. And it is more logical to consider that the benefit is related to the station and not to the different characteristics of the C-arms used. As a reminder, in the series of Hertault et al.(3), use of a hybrid operating room did not lead to a significant reduction in the volume of contrast agent used for bifurcated stent grafts (only for fenestrated/branched stent grafts). In this series, the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 13
X-ray dose reduction was very clear but it was comparing two different systems, hence it was not the fusion imaging that led to this reduction.

Finally, the recent meta-analysis by De Ruiter et al.(4) concluded: "For equivalent fluoroscopy times, the use of a fixed C-arm in noncomplex procedures leads to higher patient radiation doses compared to a mobile C-arm". We therefore believe that the combination of a mobile flat-panel detector with fusion imaging is a completely acceptable alternative to a hybrid operating room.

Limitations

In the present study, we did not report any quantitative criteria for fusion precision. Several publications [START_REF] Tacher | Image guidance for endovascular repair of complex aortic aneurysms: comparison of twodimensional and three-dimensional angiography and image fusion[END_REF][START_REF] Maurel | Evaluation of visceral artery displacement by endograft delivery system insertion[END_REF][START_REF] Kauffmann | Source of errors and accuracy of a two-dimensional/three-dimensional fusion road map for endovascular aneurysm repair of abdominal aortic aneurysm[END_REF][START_REF] Sailer | CTA with fluoroscopy image fusion guidance in endovascular complex aortic aneurysm repair[END_REF] have shown that, regardless of the registration method used, fusion precision is found to be lacking when compared to subtraction angiography. Maurel et al. [START_REF] Maurel | Evaluation of visceral artery displacement by endograft delivery system insertion[END_REF] clearly showed, with the aid of perioperative CBCT, displacement of the ostia of renal and visceral arteries in different planes due to deformations caused by insertion of a rigid material (extrastiff guide, introducer) in the aorta. Several solutions have been considered to resolve this phenomenon, such as "Image-based tracking fusion system" approaches. We decided not to measure fusion precision here because currently no system is able to predict deformations in a very precise manner. Even in cases with non-tortuous anatomy, there may be a mismatch between the fused image and angiography, which in our opinion remains essential for validating the fusion (7 ml to 30 ml/s suffice at the proximal landing zone). As long as registration is based on bony structures (hence with rigid transformations), there will be a mismatch with arterial structures which by nature are soft and deformable, which therefore calls for approaches based on elastic registration [START_REF] Nasr | Characterization of the Physiological Displacement of the Aortic Arch Using Non-Rigid Registration and MR Imaging[END_REF], for example, or digital simulation approaches that predict deformations using a biomechanical model. We have already published several studies on this simulation approach [START_REF] Gindre | Finite element simulation of the insertion of guidewires during an EVAR procedure: example of a complex patient case, a first step toward patient-specific parameterized models[END_REF][START_REF] Gindre | Patient-Specific Finite-Element Simulation of the Insertion of Guidewire During an EVAR Procedure: Guidewire Position Prediction Validation on 28 Cases[END_REF], which is integrated into the planning and fusion software. The objective of the present study was not to test the precision of registration with a simulated and deformed model. A specific methodology, which was not possible nor envisaged in the study design, needs to be developed for this purpose (using, among others, perioperative CBCT). Moreover, we have already quantified the accuracy and robustness of the registration algorithm in the dedicated methodological article presenting the detailed principle of the 3d/2d registration [START_REF] Dumenil | A versatile intensitybased 3D/2D rigid registration compatible with mobile C-arm for endovascular treatment of abdominal aortic aneurysm[END_REF]. In this already published paper, a thorough and precise validation scheme was proposed, and the mean registration error on the bony landmarks was found to be < 0.5mm. Finally, the learning curve of the team and the fact that we are a training center can lead to a bias in the reproducibility of the results and the efficiency in the EVAR procedures. In this study, there was investigators with different experiences (range from 10 to 300 EVAR procedures) but subgroups would not give enough statistical power to reach statistical differences in radiation parameters.

Conclusion

Fusion of the 3D preoperative CT scan with 2D fluoroscopy is possible with a mobile C-arm and compatible with the clinical workflow. In addition to the unquantifiable visual comfort and the possibility of using all modern navigation tools, this technique also appears to reduce the volume of contrast agent for EVAR procedures. The combination of a mobile flat-panel detector with a computer-assisted surgery station is an acceptable alternative to hybrid operating rooms for complex aortic procedures with high-performance imaging. This avoids a procedure with the C-arm sideways, which would increase the radiation dose to the operator.

6 -

 6 Patients who underwent MR angiography during preoperative evaluation.

From 1

 1 May 2016 to 30 June 2017, 54 patients (49 men, 92.6%) of mean age 73.4 ± 9.3 years were included in the study. Mean BMI was 28.3 ± 7.7 kg/m 2 . The procedures performed were: 49 bifurcated stent grafts including 2 with unilateral internal iliac artery embolisation and one with an anchoring system (Aptus, Medtronic); 2 branched iliac; 2 thoracic; and 1 fenestrated stent graft. The mean duration of the EVAR procedures was 103.6 ± 25.4 min. In this cohort, proximal and distal seals were achieved with the successful introduction and deployment of the device in the absence of surgical conversion or mortality, type I or III endoleaks, or graft limb obstruction. Considering this definition, the technical success was 100%. Percutaneous access was performed for 50 (92.6%) patients (bilateral) without open conversion by using 6 Fr Perclose Proglide device (Abbott Vascular, Redwood City, CA).

Fig. 1 :

 1 Fig. 1: The angio-navigation station is placed in front of the surgeon (black arrow), who

Fig. 2 :

 2 Fig. 2: Fusion of preoperative CT scan for EVAR using a non-deformed model (A), which is

Fig. 3 :

 3 Fig. 3: Regardless of the procedure, low-volume angiography is systematically performed to

Fig. 4 :

 4 Fig.4: For TEVAR procedures in the descending thoracic artery, the prosthesis is deployed

Fig. 5 :

 5 Fig. 5: Box plots comparing variables of patients operated on with the angio-navigation

  Glomerular filtration Rate in ml/min/1.73m -2 M A N U S C R I P T A C C E P T E D

Table 1 :

 1 Fusion feasibility rate and contrast and radiation dose by procedure type DAP = Dose-area product, AK = Air kerma

		Registration success	Contrast agent (mL)	Fluoroscopy time (min)	DAP (Gy.cm 2 )	AK (mGy)
	EVAR (n=49)	95.1% (155/163)	41.9 ± 23.1	21.6 ± 12.3	70.9 ± 48.2 254.2 ± 161.5
	TEVAR (n=2)	100% (4/4)	23.8	6	20	79
	IBD (n=2)	100% (6/6)	39.8	39.8	67.7	363.5
	FEVAR (n=1)	100% (5/5)	65	38	163.1	603

*

Table 2 EVAR Fusion

 2Fusion

			EVAR No
			fusion	P value
		(n=49)	(n=103)
	Age (years; mean ± SD)	73.7 ± 9.3	73.9 ± 9	0.89
	Gender (Male)	44 (89.8%) 97 (94.2%)	0.26
	BMI (kg/m 2 )	28.3 ± 5	26.9 ± 4	0.19
	Symptomatic PAD*	0 (0%)	4 (3.8%)	0.21
	Coronary artery lesions	12 (24.5%)	28 (27.2%)	0.84
	Severe respiratory failure	1 (2%)	4 (3.9%)	0.48
	Renal failure (eGFR<30)**	2 (4.1%)	6 (5.8%)	0.49
	Treated hypertension	30 (61.2%)	82 (79.2%)	0.06
	History of tobacco	36 (73.5%)	75 (72.8%)	0.9
	Diabetes	6 (12.2%)	7 (6.8%)	0.53
	Treated dyslipidemia	35 (71.4%)	64 (62.1%)	0.53
	Anticoagulant therapy	3 (6.1%)	4 (3.9%)	0.41

Acknowledgements

Funding

None

Conflict of interest