domain: spi math info phys

text:
A Quantized Tensor Train Method for High-Frequency Scattering Problems Involving Heterogeneous Dielectric Layers Matthias Baray, David Levadoux, Jean-René Poirier

Abstract-We present a new numerical scheme to solve efficiently scattering problems involving an elongated and flat heterogeneous dielectric material assumed to be invariant along a direction of space. The technique consists in compressing the integral operators of an integrodifferential formulation with a so-called quantized tensor-train algorithm whose use is rather original in this context. We show that it allows to compute and store operators with a notably small memory footprint while having at the same time a fast matrix-vector product leading to a competitive method compared to the more classical H-matrix approach.

Index Terms-QTT decomposition, integrodifferential equations, boundary integral equations, electromagnetic scatteringtransmission problems, FEM-BEM coupling.

I. INTRODUCTION

W E are interested in the interaction between an incident electromagnetic field and a penetrable material as encountered in some applications using high-frequency or optical properties of complex thin layers. For instance, asymptotic ray-tracing codes, used to predict the radar cross-section of some metallic structure coated with a thin layer of material, need backscattering coefficients that can be evaluated with an exact two-dimensional code. In this paper, we present a method to compute the scattered and transmitted fields of an electromagnetic plane wave impinging, in a TM mode, a twodimensional layer made of a long and thin piece of a dielectric inlaid with a stochastic seedling of small metallic square inclusions (see Figure 8). Such a problem is usually cast in an integrodifferential formulation leading, after discretization, to a linear system coupling a Finite Element Method (FEM) with a Boundary Element Method (BEM).

It is well-known that the BEM part makes the solution very accurate compared to a regular FEM using a perfectly matched layer instead of the exact boundary operators. However, a painful drawback is that BEM matrices are dense and need to be compressed. The Fast Multipole Method (FMM) [START_REF] Song | Multilevel fast-multipole algorithm for solving combined field integral equations of electromagnetic scattering[END_REF] and the H-matrix [START_REF] Hackbusch | A sparse matrix arithmetic based on H-Matrices. part I: Introduction to H-Matrices[END_REF] have made significant inroads into this topic and are now widely used. Nevertheless, their developments and uses are touchy, making attractive some newer and more robust compression algorithms. Among these are undoubtedly the tensor methods whose popularity in a wide range of domains does not cease to grow since some ten years ago [START_REF] Kolda | Tensor decompositions and applications[END_REF].

Initially designed to break the "Curse of the dimension" (exponential growth of the storage with the dimension of an array), tensor methods, especially the Quantized Tensor Train (QTT) decomposition [START_REF] Oseledets | Tensor-train decomposition[END_REF], are now able to tackle some difficult problems of computational physics such as the solution of very big linear systems [START_REF]Approximation of 2 d ×2 d matrices using tensor decomposition[END_REF]. Despite this promising surge, it is worth noting that the QTT format is still not widely used in the field of boundary integral equations, for instance, we can cite [START_REF] Corona | A tensor-train accelerated for integral equations in complex geometries[END_REF]- [START_REF] Giannakopoulos | Compression of volume-surface integral equation matrices via Tucker decomposition for magnetic resonance applications[END_REF]. The technique we investigate is a contribution in this field and consists of treating the FEM part of our linear system classically and compressing the BEM parts with a QTT algorithm. At last, we will see that the main advantage of this approach is to drastically reduce the building time and the storage of BEM matrices.

The outline of the paper is as follows. Section II presents the boundary/transmission value problem we want to solve and its translation in a FEM-BEM coupling. Section III is a brief introduction to the QTT theory, essentially aiming to explain how to compress a matrix into a sparse highdimensional tensor and how the associated fast matrix-vector product works. In section IV, we give some comparative results between the QTT and H-matrix techniques related to the BEM matrices of section II. Finally, in section V, we apply our numerical scheme to a big scattering-transmission problem involving a dielectric with stochastic inclusions, as mentioned hereinbefore.

II. SCATTERING-TRANSMISSION PROBLEM

Given a plane wave u i impinging a material Ω -with boundary Γ, we seek two fields u s and u defined resp. on Ω + := R 2 \Ω -and Ω -, satisfying the following transmission problem:

       ∆u + k 2 -u = 0, in Ω - ∆u s + k 2 + u s = 0, in Ω + u = u s + u i , on Γ ∂ n u = ∂ n u s + ∂ n u i , on Γ (1)
where k + , k -are wavenumbers (k + being constant), and ∂ n is the outer normal derivative on Γ. In addition, the scattered field u s has to satisfy the Sommerfeld radiation condition:

lim r→∞ √ r ∂u s ∂r -ik + u s = 0 . (2)
It is well-known [START_REF] Hiptmair | Stabilized Fem-Bem Coupling for Helmholtz Transmission Problems[END_REF] that a weak integrodifferential formulation of (1)-( 2) involving the boundary integral operators [START_REF] Nédélec | Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems[END_REF]: is the Hankel function of the first kind), consists of finding u N = -∂ n u and u |Γ such that for any fields ϕ and ψ defined resp. on Γ and Ω -, one has:

Sp(x) = Γ G(x,
     ⟨Su N , ϕ⟩ + ⟨(D + 1 2 )u |Γ , ϕ⟩ = ⟨u i |Γ , ϕ⟩ ⟨(D t + 1 2 )u N , ψ |Γ ⟩ + H(u, ψ) + ⟨N u |Γ , ψ |Γ ⟩ = ⟨∂ n u i , ψ |Γ ⟩ (3)
where ⟨•, •⟩ are the L 2 scalar products on Γ or Ω -and H(u, ψ) = ⟨∇u, ∇ψ⟩ -⟨k 2 -u, ψ⟩. Volumic and surfacic interactions are discretized with nodal finite elements of order 1 (P 1 approximation). The numbering of unknowns starts with Neumann surfacic basis functions ϕ i (i ∈ S), followed by volumic basis functions in touch with Γ ψ i (i ∈ V Γ ), and ends with pure volumic basis functions ψ i (i ∈ V) (S, V Γ , V being contiguous intervals of N). The numerical approximation of u N and u in (3) is given resp. by i∈S X i ϕ i and i∈V Γ ∪V X i ψ i where X ∈ C N is solution of the linear system M X = B where:

M =   S D + I 2 0 D t + I 2 N + H 1 H 2 0 H 3 H 4   (4)
with

S = (S ij ), S ij = ⟨Sϕ i , ϕ j ⟩ ∀(i, j) ∈ S 2 D = (D ij ), D ij = ⟨Dϕ i , ψ j ⟩ ∀(i, j) ∈ S × V Γ D t = (D t ij ), D t ij = ⟨D t ψ i , ϕ j ⟩ ∀(i, j) ∈ V Γ × S N = (N ij ), N ij = ⟨N ψ i , ψ j ⟩ ∀(i, j) ∈ V 2 Γ H 1 = (H 1ij ), H 1ij = H(ψ i , ψ j ) ∀(i, j) ∈ V 2 Γ H 2 = (H 2ij ), H 2ij = H(ψ i , ψ j ) ∀(i, j) ∈ V Γ × V H 3 = (H 3ij ), H 3ij = H(ψ i , ψ j ) ∀(i, j) ∈ V × V Γ H 4 = (H 4ij ), H 4ij = H(ψ i , ψ j ) ∀(i, j) ∈ V 2 I = (I ij ), I ij = ⟨ψ i , ϕ j ⟩ ∀(i, j) ∈ V Γ × S and B = B 1 B 2 0 with B 1 = (B 1i ), B 1i = ⟨u i |Γ , ϕ i ⟩ ∀i ∈ S B 2 = (B 2i ), B 2i = ⟨∂ n u i , ψ i|Γ ⟩ ∀i ∈ V Γ .
III. QUANTIZED TENSOR TRAIN DECOMPOSITION We present now a brief overview of tensor techniques we will use to compress BEM matrices (S, D, N ) involved in the linear system [START_REF] Oseledets | Tensor-train decomposition[END_REF]. Recall that a simple multi-dimensional array A ∈ C m1ו••×m d is named a tensor of dimension d with mode sizes m k . Therefore, to use tensor machinery, we first have to convert matrices into high-dimensional tensors. This transformation is called quantization [START_REF]Approximation of 2 d ×2 d matrices using tensor decomposition[END_REF], [START_REF] Khoromskij | O(dlog n)-quantics approximation of n -d tensors in high-dimensional numerical modeling[END_REF]. As FMM or H-matrix methods, the starting point is to give a hierarchical structure over the set of unknowns, meaning that we have to give us a bijective application linking any primal index to the corresponding multi-index (d-tuple) in the hierarchy:

T g : [[N ]] → [[n 1 ]] × • • • × [[n d ]] ; i → î with [[N ]] = [[1; d i=1 n i ]], [[n k ]] = [[1; n k ]] and î := (i 1 , . . . , i d ).
This application is equivalent to a group cluster tree (classical structure of H-matrix) that allows dealing with clusters of unknowns of different sizes. For example, a cluster of unknowns of level k can be read as (i 1 , . . . , i k , :, . . . , :). In the sequel, we will only use a binary tree

(n k = 2 ∀ k ∈ [[d]]
) whose inverse application is:

T -1 g ( î) = d k=1 2 d-k (i k -1) + 1, i k ∈ [[2]] . (5)
Given a group cluster tree, we can build canonically another hierarchical structure on the set of the matrix indices, the socalled block cluster tree:

T b : [[2 d ]] 2 → ([[2]] × • • • × [[2]]) 2 ; (i, j) → (T g (i), T g (j)) .
This structure recursively splits the matrix into four equal parts (we have to face a quadtree), and each entry (i, j) of the matrix is identified with a couple of multi-indices ( î, ĵ). To create the tensor related to the matrix, we just have to pair i k and j k and treat them as a single index:

i k j k ∈ [[4]] with i k j k = i k + 2(j k -1). (6)
Therefore, we can associate to any matrix

A ∈ C 2 d ×2 d a d- dimensional tensor A ∈ C 4ו••×4 defined by: A(i 1 j 1 , . . . , i d j d ) = A(T -1 b ( î, ĵ)) . (7)
This last step ends the quantization process, but the question of the matrix storage is not yet answered. To this end, several formats are possible [START_REF] Kolda | Tensor decompositions and applications[END_REF] among these is the Tensor-Train (TT) decomposition [START_REF] Oseledets | Tensor-train decomposition[END_REF], [START_REF]Approximation of 2 d ×2 d matrices using tensor decomposition[END_REF]. Such a format can be viewed as a generalization of the well-known separation of variables in the theory of functions. Indeed, a tensor A ∈ C m1ו••×m d is subject to a TT decomposition if:

A(i 1 , . . . , i d ) = G 1 (i 1 ) . . . G d (i d ) (8)
where the so-called cores i k → G k (i k ) ∈ C r k-1 ×r k are matrix-valued functions and products between cores, matrixmatrix products (see Figure 1 for an illustration). The dimen-sions of the cores, i.e., the r k , are called TT-ranks. Tensor elements being scalars, it follows r 0 = r d = 1.

Applied to a tensor resulting from the quantization of a matrix [START_REF] Poirier | Fast BEM Solution for 2-D Scattering Problems Using Quantized Tensor-Train Format[END_REF] the TT decomposition becomes:

A(i 1 j 1 , . . . , i d j d ) = G 1 (i 1 j 1 ) . . . G d (i d j d ) (9)
The quantization followed by a TT-formatting is called Quantized Tensor Train (QTT) decomposition, and the QTT decomposition of a matrix is called a TT-matrix. An interesting example is the TT-matrix related to the full matrix A ∈ C 2 d ×2 d :

A ij = e ik(xi-xj ) , ∀i, j ∈ [[2 d ]] with k ∈ R, x i = (i -1)h and h = 1 2 d -1 .
With the binary tree (5) we have

x i -x j = h d l=1 2 d-l (i l -j l ) .
Applying the exponential and using ( 6) we obtain the TT-matrix A of A:

A(i 1 j 1 , . . . , i d j d ) = e ik2 d-1 h(i1-j1) • • • e ikh(i d -j d ) = G 1 (i 1 j 1 ) • • • G d (i d j d ) with G l (i l j l ) = e ik2 d-l h(i l -j l ) , l ∈ [[d]].
As in this decomposition, all TT-ranks are equal to 1, the storage of the TT-matrix (i.e. cores G l ) is 4d instead of 2 d × 2 d . Nevertheless, for an arbitrary TT-matrix (structured with a quadtree), the storage is just bounded by 4dr 2 with r = max(r k ), showing that it is crucial to deal with a class of low-rank (r ≪ N ) tensors. The situation is even more critical in practice, where there is no chance to have a finite decomposition due to the complexity of our oscillating kernel functions. Therefore we have to build approximations of them with a small rank r, which is done by some algorithms. Among these, we choose the very attractive so-called TT-cross working only on sparse sampling of the matrix (similar to H-matrix). We refer to [START_REF] Oseledets | TT-cross approximation for multidimensional arrays[END_REF] for details of this technical algorithm.

As we plan to iteratively solve the linear system governed by (4), let us now examine the Matrix-Vector Product (MVP). A QTT-compression applied to the whole matrix (4) is, a priori, not possible due to the unstructured parts (FEM blocks). A simple solution is to use a "mixed" TT-MVP [START_REF]Approximation of 2 d ×2 d matrices using tensor decomposition[END_REF] that involves a TT-matrix with a vector in regular format (decimal numbering). Although there exists a "pure" TT-MVP [START_REF] Oseledets | Tensor-train decomposition[END_REF] (matrix and vector in QTT-format), it would imply to change the format of the vector (regular vs. QTT) at each iteration. Besides, the "mixed" product (simply named TT-MVP hereafter) is very easy both to understand and implement. Indeed, if

A ∈ C 2 d ×2 d and x ∈ C 2 d
, the regular MVP reads as:

y i = 2 d j=1 A ij x j , ∀i ∈ [[2 d ]] . (10)
Replacing in [START_REF] Nédélec | Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems[END_REF] native entries of A and x with their tensorial counterparts A (TT-matrix ( 9)) and X (quantization (5)) we obtain:

y(T -1 g ( î)) = j1,...,j d A(i 1 j 1 , . . . , i d j d )X (j 1 , . . . , j d ) = j1,...,j d G 1 (i 1 j 1 ) • • • G d (i d j d )X (j 1 , . . . , j d ) .

Algorithm 1 TT-MVP

In : A TT-matrix A of A with cores G k , a vector x.

Out : y = Ax 1: Y = x 2: for k = 1 to d do 3: Y = reshape(Y, [2r k-1 , numel(Y) 2r k-1 ]) 4: Y = G k • Y 5: Y = reshape(Y, [2, numel(Y)
G k := reshape(G k , [r k-1 , 2, 2, r k ]) G k := permute(G k , [2, 4, 1, 3]) G k := reshape(G k , [2r k , 2r k-1 ]) .
Therefore, the TT-MVP algorithm simply consists in performing d successive small matrix-matrix products, and its complexity is O(r 2 N log(N )) where r is the higher rank of the TT-matrix. Concerning the storage, the intermediate arrays Y count r k-1 N (line 3) and r k N (lines 4-6) entries leading to a complexity of O(rN ).

IV. PERFORMANCE ANALYSIS OF QTT-FORMAT FOR INTEGRAL OPERATORS

Through some numerical experiments and comparisons with H-matrix methods, this section aims to evaluate the performance rates of the CPU building time, storage, and MVP speed related to BEM matrices in QTT-format. We use two open-source toolboxes, the Python package ttpy [START_REF] Oseledets | ttpy[END_REF] for QTTformat and the Matlab toolbox Gypsilab [START_REF] Aussal | Gypsilab[END_REF] for H-matrix. The precision parameter (error of compression) related to each method is fixed to 10 -4 for QTT and 10 -3 for H-matrix. This choice is purely pragmatic, the H-matrix one is conventional and the QTT one results from experimental tests ensuring approximately the same precision.

The BEM operators under test are defined on the boundary Γ of the layer used in section V: a long and thin rectangular contour of size 1 × 3.10 -4 . The mesh is regular and the unknowns are ordered according to the curvilinear coordinates. CPU building times, CPU MVP times and storage ratings (resp. named CPU build , CPU mvp , MEM build ) are evaluated by summing each score related to S, D, and N matrices. The notation for the slope of a performance plot is p T T in the case of the QTT method and p H for the H-matrix. All computations are performed sequentially on a laptop with a 2.8 GHz Intel Core i5 processor and 8 GB of RAM.

A. Building time and storage

Firstly, we are interested in the building times for a frequency sweeping. The mesh of Γ carries 131 072 unknowns and the wavelength λ ranges from 20h to 20 480h (h being the edge size). n λ = 1/λ is the radioelectric size of Γ. The Figure 2 shows that building a BEM matrix with TTcross is very competitive compared to the H-matrix, around 60 times faster (mind that the building time for the Hmatrix method is measured in minutes). This performance is correlated to a spectacular compression rate as shown in Figure 3, saving 3 orders of magnitude in memory with QTT (the storage for the H-matrix method is measured in gigabytes). Moreover, all these good results are stable from the highfrequency (with an optimal mesh λ/h = 20) to the lowfrequency regime (without remeshing).

After a frequency sweeping with a constant mesh, it is interesting to examine what happens when the frequency and the number of unknowns increase in a correlated manner with a fixed λ/h = 20. We use now a frequency sweeping that corresponds to radioelectric sizes ranging from 819λ to 26 206λ (from 0.1λ to 8λ in the thickness).

The Figure 4 shows that the QTT integral operators are still quicker to build than with H-matrix, almost again sixty times faster. Besides, the QTT method can build the matrices with more than one million unknowns in less than 5 minutes while it can not be done in a reasonable amount of time with H-matrix.

About the memory consumption, Figure 5 confirms the previous results of Figure 3 showing that even exploring the high-frequency regime we always save 3 orders of magnitude compared to H-matrix techniques. Moreover, the memory consumption falls down with QTT to a sub-linear law, contrary to H-matrix which theoretically reaches O(n log n) and in practice here reaches a little over-linear law.

B. Matrix-vector product

We evaluate now the cost of the TT-MVP compared to the H-matrix product. We use a random vector because it does not impact the efficiency of the methods (well-known for Hmatrix and obvious for QTT in view of Algorithm 1).

First, let us check the storage of the intermediate arrays Y (Alg. 1) which can reach rN entries and be too large if r is not small enough. Fortunately, this is not the case in our applications. Indeed, the largest storage is given by the single layer operator at the highest frequency (n λ = 26 206, N = 1 048 576) and amounts to 1.4 GB (with r = 83). This storage is not prohibitive, being around a third of the H-matrix one.

Let us now examine the CPU time of the TT-MVP. Whatever the kind of the frequency sweeping (with a fixed or adaptive mesh, see the former part A), Figures 6 and7 show that the two MVPs are comparable in complexity, but the Hmatrix product is twice as fast.

These results are disappointing regarding the spectacular compression of the QTT format. However, they do not call into question the use of TT-MVP in an iterative solver for our application. Furthermore, it is worth noting that the small memory footprint of the TT-matrices would allow to implement a simple and efficient MPI parallelization of the TT-MVP (copying of the TT-matrix to each processor being not an issue). Finally, let us mention that to our knowledge, generalization of such results to 3D applications is still an issue mainly because of the drastic increase of the TT-ranks. For instance, the first results we made on a thin square plate with an EFIE operator show that despite saving again two orders in magnitude for memory compared to H-matrix, we lose one order of magnitude for the CPU mvp .

V. A QTT-BASED SOLVER FOR A TRANSMISSION PROBLEM WITH A HETEROGENEOUS DIELECTRIC LAYER

In this section, we focus on the resolution of the problem [START_REF] Song | Multilevel fast-multipole algorithm for solving combined field integral equations of electromagnetic scattering[END_REF] where Ω -is a rectangular stochastic material (cf. Figure 8) of permittivity 2+i and radioelectric size 2200λ + × 0.67λ + = 3275.8λ -× 1λ -. We point out that Ω -is an absorbing material as our radiation condition (2) goes with a e -iωt time behavior. To generate the seedling, we use 4 random variables: δ x , δ y ∼ U( -l 20 , l 20 ), γ ∼ U(0, l 5 ) and θ ∼ U(0, π 2 ), U being the uniform distribution and l the material thickness. Each of the 3088 metallic squares are characterized by their center

(x i +δ x , y i +δ y ) (x i = l √ 2 2 +l(i-1) 3 √ 2 4 , y i = l 4 , i ∈ [[ 3088
]]), edge size l 2 -γ and rotation angle θ. We compare the QTT and H-matrix methods using the precision and cluster trees used in section IV. Ω -is meshed with 1 366 895 unknowns and the rectangular boundary Γ with 131 072 unknowns (leading to #S = 131 072, #V Γ = 131 072, and #V = 1 235 823). The material is horizontally elongated and illuminated from below by a plane wave making an incidence angle with the x-axis equal to 26.57 • .

We use the GMRES-like solver described in [START_REF] Soudais | Iterative solution of a 3-d scattering problem from arbitrary shaped multidielectric and multiconducting bodies[END_REF] with a stopping residue of 10 -3 . The native linear system (4) is preconditioned with a Schur complement technique. Precisely, the matrix of the linear system (4) be written as:

M = M 1 M 2 M 3 M 4 with M 1 = S D + I 2 D t + I 2 N + H 1 , M 2 = 0 H 2 , M 3 = 0 H 3 , M 4 = H 4 .
The Schur complement related to the block M 4 is defined as

M 1 -M 2 H -1 4 M 3 .
Before the iterative solution, we factorize H 4 with an LU decomposition. At each iteration, we compute H -1 4 with a forward and backward substitution. Figure 9 is an overview of the entire simulation cost (CPU and memory) to solve the scattering-transmission problem. The left side shows that the QTT method drastically reduces the CPU building times. For both methods, the solver converges in 35 iterations. Nevertheless, the QTT method does not improve the solution time because the TT-MVP is a little less effective than the H-matrix product (actually 1.5 times H-MVP as shown in Figure 7). However, solution times are almost the same because the 1.5 GB required to store H matrices nearly fill the remaining memory of our computer. Moreover, in agreement with the results of the section IV, the memory cost (right side) of the TT-matrices is quite negligible (note the log scale) compared to H-matrix. Tables I and II give details of the CPU building times for BEM matrices and their storage.

The relative error between the QTT and H-matrix solutions (resp. u QTT and u H ) gives: We see that replacing H-matrix compression by QTT method is safe from an engineering point of view. For instance, Figure 10 shows that the two bistatic radar cross-sections (from 0 • to 360 • by turning counterclockwise) are superimposed.

∥u QTT -u H ∥ 2 ∥u H ∥ 2 = 1.239 × 10 -3 .
In light of these rather good results, we can guess what happens if the shape of the tested material is not exactly linear but slightly curved or undulated. The answer is given in Figure 11, where we used sine shape profiles: a half-sine (curved) and a 4-sines train (undulated) whose amplitudes range from -4λ -to +4λ -(in order to see the different shapes, the drawings in the figure are not to scale). We see that the QTT method remains globally efficient even if results are a little disappointing compared to [START_REF] Poirier | Fast BEM Solution for 2-D Scattering Problems Using Quantized Tensor-Train Format[END_REF] or [START_REF] Corona | A tensor-train accelerated for integral equations in complex geometries[END_REF] where rough surfaces or irregular lattices were treated spectacularly. Actually, it appears that applying tensor methods in the context of high-frequency regime is still challenging, and our results probably call to experiment other tensor formats, more suited to these special situations (as perhaps the Tucker format successfully used in [START_REF] Giannakopoulos | Compression of volume-surface integral equation matrices via Tucker decomposition for magnetic resonance applications[END_REF]).

  y)p(y)dy + ∥x -y∥) (where x, y ∈ Γ and H(1) 0

Fig. 1 :

 1 Fig. 1: TT decomposition of a 5-dimensional tensor. A(1, 2, 4, 3, 3) is obtained by the product of the darkened matrices.

  end for 8: y = reshape(Y, [numel(Y), 1]) This last expression gives the Algorithm 1 where the cores G k are reshaped in G k by the well-known Matlab functions:

Fig. 2 :

 2 Fig. 2: CPU build as a function of the radioelectric size.

Fig. 3 :

 3 Fig. 3: MEM build as a function of the radioelectric size.

Fig. 4 :Fig. 5 :

 45 Fig.4: CPU build as a function of the number of unknowns.

Fig. 6 :

 6 Fig. 6: CPU mvp as a function of the radioelectric size.

Fig. 7 :

 7 Fig. 7: CPU mvp as a function of the number of unknowns.

Fig. 8 :

 8 Fig. 8: Test case.

Fig. 9 :

 9 Fig. 9: Overview of the computation times (s) and the storage.

Fig. 10 :

 10 Fig. 10: Bistatic radar cross-sections.

Fig. 11 :

 11 Fig. 11: Overview of the computation times (s) and the storage for different shapes of material.

TABLE I :

 I Integral operators building time (s)

	Format	S	D	N
	QTT	11.902	11.669	10.595
	H-matrix 258.587 684.674 552.472

TABLE II :

 II Integral operators storage (MB)

	Format	S	D	N
	QTT	0.97	0.59	0.30
	H-matrix 548.34 525.55 641.38

ACKNOWLEDGMENT

We would like to thank the Agence Innovation Défense (AID) and the ONERA for funding this work. We are also grateful to Bastiaan Michielsen and Vincent Mouysset for their critical reading and providing valuable feedback.


Matthieu Masson
email: matthieu.masson@laas.fr.

Denis Arzelier
email: arzelier@laas.fr.

Mioara Joldes
email: joldes@laas.fr.

Bruno Revelin
email: bruno.revelin@thalesgroup.com.

Jérôme Thomassin
email: jerome.thomassin@cnes.fr.

Keywords: Collision Avoidance Maneuver, Multiple Short-term Encounters, Impulsive Maneuvers, Quadratically Constrained Nonconvex Optimization

In this article, the problem of modeling and computing the optimal collision avoidance maneuvers for multiple short-term encounters is presented. Several metrics defining the overall collision risk are analyzed. These include imposing lower bounds on the miss-distance or the Mahalanobis distance at each corresponding time of closest approach (TCA) or upper-bounds on the orbital collision probability. The avoidance maneuvers are modeled as impulsive ones in a single direction of the local frame and for a priori fixed dates for operational motivations. Station-keeping constraints are also imposed via linear inequalities on the relative states at each TCA. The specific nature of the imposed constraints is dictated by a practical framework provided by the French Space Agency (CNES). This results in formulating the maneuver design problem as an optimization problem with linear objective and nonconvex quadratic constraints. Different algorithms are presented to solve this problem. Finally, the relative efficiency of the proposed approaches is evaluated and analyzed on two realistic conjunctions built from data extracted from the CNES database and on academic examples of higher complexity.

Multi-maneuvers algorithms for multi-risk collision avoidance via nonconvex quadratic optimization

INTRODUCTION

Nowadays the risk of orbital collision for controlled satellites, especially in low orbits, needs to be assessed and mitigated by space agencies and owners/operators (O/O) of the field. When the predicted risk is too high according to mission requirements, a collision avoidance strategy usually consists in one or several evasive maneuvers. These are planned by taking into account both fuel consumption minimization (for preserving the lifetime of the satellite) and other operational constraints, like station-keeping by the time the threat is gone. More demanding operational constraints like disjunctive thrusts, minimum elapsed time between two thrusts, saturation constraints, etc., are often necessary when low-thrust (or electric) propulsion is used.

In this study, the collision avoidance problem between one active low-thrust satellite involved in short-term encounters with several uncontrolled objects, is considered. The active satellite, also called primary p, is initially set on a reference orbit. The other involved objects, called secondaries s j , with j = 1, . . . , N, are assumed to be uncontrolled space debris. The use of low-thrust propulsion is modeled, like in the work of [START_REF] Hennes | Fast approximators for optimal low-thrust hops between main belt asteroids[END_REF], by a sequence of impulsive thrusts on specifically imposed thrust ranges over which the thrust magnitude is spread. Furthermore, due to the lack of precision in measurements, the quantification of the collision risk is also

The Authors would like to thank CNES for its partial financial support of this work and Thales Alenia Space for giving access to realistic space data. a difficult task and different risk assessment indicators include, for each pair (p, s j ): (1) the miss-distance at the Time of Closest Approach (TCA); (2) the Mahalanobis miss-distance at TCA;

(3) the short-term or (4) maximum collision probability (see for instance [START_REF] Fernández-Mellado | On the use of machine learning and evidence theory to improve collision risk management[END_REF] and Sec. 2 for details).

Hence, in this context, the collision avoidance problem is formulated as a constrained optimization problem: it consists in computing a fixed-time fuel optimal, finite sequence of impulsive maneuvers performed by the active spacecraft such that (one of) the risk reduction criteria above, as well as saturation and operational constraints, are met.

Even in this simplified and restricted theoretical framework, there are relatively few existing works like those of [START_REF] Duncan | Optimal collision avoidance for multiple conjunction events[END_REF]; [START_REF] Kotz | A study on the collision avoidance maneuver optimization with multiple space debris[END_REF], which are based on defining a penalty function as the sum between the objective function and a penalty parameter multiplied by a measure of violation of the constraints. The resulting non-constrained optimization problem is solved either with derivative-free methods or by evolutionary optimization methods.

In this work, we revisit and analyze several reformulations of the Collision Avoidance Maneuver (CAM) design problem, when considering short-term encounters, in a specific practical framework provided by the French Space Agency (CNES). By extending and adapting the methods of [START_REF] Slater | Collision Avoidance for Satellites in Formation Flight[END_REF], [START_REF] Bombardelli | Optimal impulsive collision avoidance in low Earth orbit[END_REF] and [START_REF] Armellin | Collision avoidance maneuver optimization with a multiple-impulse convex formulation[END_REF], the formulation of nonconvex Quadratically Constrained Linear Programming (QCLP) problems is proposed in Section 2.

These are known to be NP-hard [START_REF] Vavasis | Complexity issues in global optimization: A survey[END_REF]) as soon as more than one quadratic constraint is present, with few algorithms existing in the literature ( (Floudas and Visweswaran, 1995, Sec. 8)). We focus in Section 3 on two distinct classes of numerical approaches composed each of two different methods. The first two methods rely on SemiDefinite Programming (SDP) as a common basis, with the first one exploiting recent results in semi-algebraic optimization, while the second one rather uses a randomization step. The second class consists of branch and bound (B&B) algorithms with set-based arithmetic (boxes or ellipsoids). On the one hand, the more sofisticated and mature approach of the IBEX software 1 is considered, which uses contractors (based on a filtering concept borrowed to constraint programming). On the other hand, a prototype experimental algorithm inspired by the ellipsoid method is investigated.

While lacking exhaustiveness, our choice of numerical approaches has a twofold motivation: firstly, we are interested in providing a stand-alone algorithm that can be directly included in existing CNES software tools; secondly, we considered both general well-established routines that converge towards global optimality and more ad-hoc procedures that are efficient on practical CNES examples. Therefore, in Section 4, we compare and evaluate these different approaches on realistic space and academic examples, using IBEX as the reference method because of its greater maturity.

Notations: O p×m and I m denote respectively the null matrix of dimensions p × m, the identity matrix of dimension m.

N n d := {α ∈ N n : |α| ≤ d}, where |α| = ∑ i α i and α = is a multi- index. Indeed, y α = y α 1 1 • • • y α n n with |α| ≤ d and for α ∈ N n d .
For a symmetric real matrix S ∈ R n×n , the notation S 0 (S 0) stands for the negative (positive) semi-definiteness of S while u v means the component-wise inequality for the two vectors u and v of identical dimensions. For the random vector ξ , E(ξ ) is its expectation. Finally, a is the nearest integer greater than or equal to a. Finally, ∂ f is the subdifferential (set of subgradients) of the function f .

A GENERIC CAM DESIGN PROBLEM FORMULATED AS A LINEAR OPTIMIZATION PROBLEM WITH NON-CONVEX QUADRATIC CONSTRAINTS

The objective of this section is to propose a generic framework for CAM design taking into account station-keeping requirements, a variety of operational constraints and a great diversity of possible missions.

However, for the sake of practicality, we closely follow CNES operational procedures and insight, which lead to a specific CAM problem, defined as follows. Besides the low-thrust specific operational constraints mentioned in the introduction, we consider station-keeping requirements and also the necessity to reduce the impact of CAMs on mission service by forbidding attitude changes in the avoidance phase. Specifically, the primary is actuated by one or more impulsive CAMs that must take place over user-predefined time ranges

∆ i , with i = 1, • • • , n.
It is assumed that only the central date t i ∈ ∆ i is allowed per maneuver in order to be able to perform an a posteriori spread of the thrust over the time ranges for an equivalent low-thrust propulsion. It is also assumed that only one identical thrust direction for all ∆ i is possible. More precisely, the thrust vector 1 http://www.ibex-lib.org is defined at t i by ∆V T i = ∆v i β T l , where ∆v i is the magnitude of the impulse and β l is the l-th vector of the canonical basis in the local frame, where the index l is set to 1, 2 or 3 for all i = 1, • • • , n and ∆v i ≥ 0. The magnitude of the thrusts is also bounded as ∆v i ≤ ∆v i . Furthermore, the relative motion of each secondary object s j , with j = 1, . . . , N, is described by a transition matrix obtained by a linearization of the Keplerian dynamics around the mean ballistic orbit of the primary (cf. [START_REF] Yamanaka | New state transition matrix for relative motion on an arbitrary elliptical orbit[END_REF]):

x r j (t|x 0 r j ) = Φ j r (t,t 0 )x 0 r j ,

where x 0 r j denotes the initial relative state of the j-th object wrt. the primary and the transition matrix Φ j r is given by blocks as follows

Φ j r (t,t 0 ) := Φ j r 11 (t,t 0 ) Φ j r 12 (t,t 0 ) Φ j r 21 (t,t 0 ) Φ j r 22 (t,t 0 )
.

Denoting xr j the relative state of the j-th object with respect to the primary (after a maneuver), one has, after n impulsive maneuvers in the β l direction, applied respectively at t i with a magnitude ∆v i :

xr j (t|t 0 ) = x r j (t|t 0 ) + n ∑ i=1 Φ l j r (t,t i )∆v i = Φ j r (t,t 0 )x 0 r j + n ∑ i=1 Φ l j r (t,t i )∆v i , (2)
where Φ l j r (t,t i ) is the 3 + l-th column of the matrix Φ j r (t,t i ). From Equation (2), it is straightforward to derive the relation between the uncontrolled (position r r j and velocity v r j ) and controlled (position rr j and velocity ṽr j ) relative states for each encounter at t j TCA :

rr j (t j TCA ) = r r j (t j TCA ) + n ∑ i=1 Φ l j r 12 (t j TCA ,t i )∆v i , (3)
ṽr j (t j TCA ) = v r j (t j TCA ) + n ∑ i=1 Φ l j r 22 (t j TCA ,t i )∆v i . (4)
Equations ( 3) and (4) allow for the estimation of the miss distance, the Mahalanobis distance as well as the collision probability at each t j TCA to be used as metrics for the CAM triggering.

Miss distance metric formulation

Assuming a rectilinear uniform relative motion around t j TCA of directional unit vector:

e j = ṽr j (t j TCA ) ṽr j (t j TCA ) , (5)
the new miss-distance d new j is given by the norm of the projection ∆r j min of the vector rr j (t j TCA ) on the encounter B-plane orthogonal to e j . Formally, this reads:

∆r j min = (I 3 -e j e T j )r r j (t j TCA ), (6)
where I 3e j e T j is the (symmetric and idempotent) projection matrix on the j-th B-plane. One obtains:

d new 2 j = r r j (t j
TCA ) T (I 3e j e T j )r r j (t j TCA ) +2r r j (t j TCA ) T (I 3e j e T j ) Φl j r 12 ∆v +∆v T Φl j T r 12 (I 3e j e T j ) Φl j r 12 ∆v,

where ∆v T contains the magnitudes of the n impulses:

∆v T = [ ∆v 1 , ∆v 2 , • • • , ∆v n ] ,
and similarly, Φl j r 12 is given by the concatenation of the columns

Φ l r 12 (t j TCA ,t i ) for i = 1, • • • , n.
Hence, d new 2 j is a convex quadratic form:

d new 2 j := ∆v T Q j ∆v + 2q T j ∆v + p j , (8)
where Q j , q j , p j are identified from Equation ( 7). Note that the matrices Q j are semidefinite positive ones and imposing a prescribed minimal miss distance threshold d j after the CAMs, at each t j TCA , defines a concave domain:

∆v T Q j ∆v + 2q T j ∆v + p j ≥ d j , j = 1, • • • , N.
(9) It turns out that the choice of the metric characterizing the collision (miss distance, Mahalanobis distance, 2D probability, maximum probability) has no direct influence on the nature of the obtained problem, since in all cases, the same class of problems is obtained in fine, when considering the same assumptions used in the multiple short-term encounters framework (see also [START_REF] Armellin | Collision avoidance maneuver optimization with a multiple-impulse convex formulation[END_REF] for further details on different metrics) defined by CNES.

Linear constraints: Station-keeping and control bounds requirements

The after-CAMs station-keeping requirements impose the primary to remain in a box defined with respect to the reference bounds on the relative positions for each direction in a local frame at each t j TCA . This generates a set of linear inequalities on the decision variables (the inequality sign is meant componentwise):

u j r -r r j (t j TCA ) Φl j r 12 ∆v u j r -r r j (t j TCA ), j = 1, • • • , N, (10
) for user set bounds u j r and u j r . Similarly, another set of linear inequalities on the decision variables is added due to bounds on the magnitude of the impulsive thrusts: 0 n×1 ∆v ∆v.

(11) Finally, our goal is to minimize the impact of the CAMs on the satellite lifetime and so, the minimum-fuel multi-risk collision avoidance problem is formulated as a non-convex quadratically constrained optimization problem with a linear objective:

min ∆v∈R n c T ∆v s.t. ∆v T Q j ∆v + q T j ∆v + p j ≥ d j , j = 1, • • • , N, A∆v b, ∆v 0 n×1 (12)
where

c T = [ 1, 1, • • • , 1 ] ∈ R n
, and Q j , q j , p j , A and b are identified from Equations ( 9), ( 10), ( 11). The optimization problem (12) consists in minimizing a linear function on a domain given by the intersection between the polyhedron defined by the linear inequalities and the non-convex domain defined by the exterior of a union of (possibly degenerate) ellipsoids.

The feasible domain of this optimization problem is therefore non-convex (and possibly non-connected in some cases).

NUMERICAL ALGORITHMS

As already mentioned, Problem ( 12) is a special case of nonconvex QCLP problems, for which two classes of algorithms (based either on SDP or B&B with set-valued arithmetic) are to be summarized and then compared on numerical examples in the next section. Starting with the first category, note that Problem ( 12) is a polynomial optimization problem of the form:

f * = inf x∈R n f (x), s.t. g (x) ≥ 0, = 1, . . . , r, (13)
where f and g are multivariate polynomials (or affine functions). Let K denote the feasible set of (13).

A moment-SOS hierarchy of SDP relaxations

Lasserre's hierarchy [START_REF] Lasserre | An Introduction to Polynomial and Semi-Algebraic Optimization[END_REF]) is a well-known method to approximate the optimal value of such problems based on semidefinite programming (SDP). From a computational viewpoint, a hierarchy of SDP problems is solved, which provides lower bounds that are convergent to the optimal value of the original problem. However, for large scale problems, the size of the involved SDP may be prohibitive, as well as the numerical extraction of the corresponding optimizers. In (Lasserre, 2015, Section 5.1), it is shown that:

f * = inf ν∈M(K) + K f dν s.t. ν(K) = 1, (14)
where M(K) + is the space of finite nonnegative Borel measures on K. The Problem ( 14) is a linear programming problem in an infinite-dimensional space and as such, it is difficult to get a direct solution in terms of the decision variable ν which is an unknown measure. However, this unknown measure ν may be represented by the infinite sequence of its moments

y α = K x α dν, ∀ α ∈ N n d
and the Problem ( 14) may be recast as:

f * = inf y    ∑ α∈N n d f α y α : ∃ ν ∈ M(K) + s.t. y 0 = 1    .
Thanks to strong results obtained in the framework of the (real) K-moment problem (i.e. constructing the sequences (y α ) that are moment sequences of a nonnegative measure ν with its support contained in set K) for a compact basic semi-algebraic set K (see [START_REF] Lasserre | An Introduction to Polynomial and Semi-Algebraic Optimization[END_REF], Section 2.7.2 and Section 2.7.3) for technical details), a hierarchy of semidefinite relaxations of Problem (13) may be defined as the following:

f * d = inf y L y (c T x) s.t. M d (y) 0 M d-ς l (g y) 0, = 1, • • • , r y 0 = 1, ( 15
)
where

g 0 = 1, ς = (deg(g ))/2 , = 0, • • • , r, d ≥ d 0 = max{ (deg( f ))/2 , max =0,••• ,r ς }. L y is the Riesz linear func- tional, M d (
f * d ≤ f * ; else if d < d then d ← d + 1; Goto 1; end if Global optimality test if rank(M d-ς (y * d )) = rank(M d (y * d )) = k with ς = max ς then f * d = f * and ∃ at least rank(M d (y * d )) = k global minimizers y * that may be extracted. end if if (rank(M d-ς (y * d )) = rank(M d (y * d ))) and (d < d) then d ← d + 1; Goto 1; else Stop. f * d is only a lower bound of the optimum: f * d ≤ f * . end if

Semidefinite relaxation and randomization approach

A straightforward convex relaxation of the Problem (12) may also be derived by adding decision variables allowing to linearize the problem as follows [START_REF] Luo | Semidefinite relaxation of quadratic optimization problems[END_REF]). Let X denote a n × n symmetric matrix, then the Problem (12) can be written as: min

x,X c T x

s.t. trace(Q j X) + q T j x + p j ≥ 0, j = 1, • • • , N, Ax b, x 0, X -xx T = 0. (16)
The last equality constraint X -xx T , in which the non convexity of the problem is aggregated, is now replaced with a convex constraint Xxx T 0 n×n leading to the SDP relaxation:

f * sdp = min x,X c T x s.t. trace(Q j X) + q T j x + p j ≥ 0, j = 1, • • • , N, Ax b, x 0, 1 x T x X 0,
(17) after using a Schur complement argument on the last constraint. This derivation may be considered as a particular instance of a lifting technique. Note also that (17) is the Lagrangian bidual of the nonconvex original problem. In addition, the relaxation (17) is exactly the relaxation (15) of ( 13) for d = 1. As a convex relaxation, the Problem (17) provides in general only a lower bound of the optimal solution of the genuine nonconvex Problem (12). Moreover, an optimal pair (X * , x * ) solution of (17) may not be feasible for the Problem ( 16) and x * cannot be used in general as a solution of (12). However, a simple stochastic interpretation of the pair (X * , x * ) will help in finding a systematic way to generate samples in the feasible set of (12). Then, the best sampled feasible point provides a good approximation x * srr of the optimal solution of (12). Indeed, for an optimal pair (X * , x * ) solution of ( 17), the matrix Xx xT is a covariance matrix. Then, by defining the random vector ξ distributed as the normal vector ξ ∼ N (µ = x * , Σ = X *x * x * T ), it is simple to show that (µ, Σ) is a pair solution of the following stochastic programming problem.

min µ,Σ E(ξ ) = c T µ s.t. E(g j (ξ )) ≥ 0, j = 1, • • • , N, E(Aξ ) = Aµ b, E(ξ ) = µ 0, Σ 0. ( 18
)
where E(•) is the expectation operator. Loosely speaking, ξ is an "average" solution of ( 12) over the Gaussian distribution, meaning that the expected value of the objective is minimized and the constraints are satisfied in expectation. This interpretation provides a procedure to build a representative sampling of the feasible set of (13). It consists in drawing L samples ξ (i) from the truncated (0 ξ (i) ∆ v) Gaussian distribution ξ ∼ N (µ, Σ). If the sample ξ (i) does not satisfy the constraints in ( 12), project it onto the feasible set K as ξ (i). This procedure is now detailed in the Algorithm 2 given below.

Algorithm 2 SRR(c, Q j , q j , p j , A, B)

Require: Cost vector c, matrices Q j , j = 1, • • • , N, A, vectors q j , j = 1, • • • , N, b, points p j , j = 1, • • • , N. Ensure: Approximation x *
srr of the optimal solution of Problem (13).

Solve of SDP relaxation (17) to get x * et X * ; Feasibility test of the solution of the relaxation if x * ∈ K then

x * ← x * ; Stop ; else Draw L samples ξ (i), i = 1, • • • , L from the multivariate Gaussian N ( x * , X *x * x * T ) tuncated in the set {0 x ∆ v} ; Feasibility test of the samples

for i ← 1 to L do if ξ (i) ∈ K then x(i) ← ξ (i) ; else for j ← 1 to N do Compute val( j) = ξ (i) + 0.5 * Q i * q i T * Q -1 j * ( q T j Q -1 j q j 4 - p j ) -1 * ξ (i) + 0.5 * Q i * q i ; end for val min ← min j=1,••• ,N val ; jmin ← arg min j=1,••• ,N val ; if val min < 1 then ξ (i) = (( ξ (i)+0.5 * Q jmin * q jmin )/ √ val min )-0.5 * Q jmin * q jmin ; end if Feasibility test of the scaled sample ξ (i) in K if ξ (i) ∈ K then x(i) ← ξ (i) ; end if end if end for x * srr ← arg min i c T x(i) ; end if

B&B algorithms with set-based arithmetic

An alternative classical approach to problem ( 12) is based on subsequent subdivision of an initial search domain (defined in the space of decision variables and including the initial feasible domain), via B&B algorithms. At each iteration, one discards subdomains found to be infeasible or subdomains whose minimum attained objective value (by any of its feasible points) is greater than the current best objective value found so far by the algorithm. The subdivision process stops for a subdomain when the difference between a lower bound and an upper bound of the objective function over the subregion is less than a given tolerance. The list of remaining subdomains contains all the optimal solutions. The main challenge lies in obtaining good lower and upper bounds of the objective function in a given subregion.

A solution consists in representing the domains by boxes (axesaligned interval tensors in R n ) and evaluating the involved functions by interval analysis (see [START_REF] Fernández | Interval tools in branch-and-bound methods for global optimization[END_REF] for a recent review of interval branch-and-bound based methods for global optimization).

IBEX Software. For instance, [START_REF] Ninin | Global optimization based on contractor programming: An overview of the ibex library[END_REF] describes Ibex implementation with constraint propagation techniques suitable for constrained global optimization problems. To accelerate convergence, contractors are used at each iteration to prune the width of boxes. Inspired by the filtering concept in constraint programming, the goal of a contractor is to eliminate unfeasible parts of a domain. For instance, some contractors of IBEX (which are used by default) are based on linear relaxation: given a system of linear inequalities, a box is contracted to the hull of the polytope (the set of feasible points); non-linear constraints are automatically linearized by Taylor or affine arithmetic. Also, constraint propagation algorithm HC4 and shaving contractors (which aim at elliminating a slice of a box by calling a contractor on it) are used by default. [START_REF] Ninin | Global optimization based on contractor programming: An overview of the ibex library[END_REF] mentions that the performance of the default optimizer is comparable to the global optimizer BARON, while also performing all the operations in a set-based arithmetic, which ensures that no solution is lost. As it will be shown in Section 4, this performance claim is confirmed on the tests specific to our problem.

Ellipsoid-based B&B. Since our specific problem is the QCLP-type, a natural set-based extension of interval algorithms is to represent the underlying subdomains by ellipsoids. Starting from an ellipsoid containing the initial feasible domain, the idea is to generate a sequence of ellipsoids of smaller volumes by employing cutting planes, like in the classical ellipsoid algorithm of Yudin and Nemirovski (see [START_REF] Bland | The ellipsoid method: A survey[END_REF] for a detailed historical and technical description).

More precisely, the minimization of a convex function f : R n → R, on an elliptic domain ε 0 = {x ∈ R n : (x-x 0 ) P 0 (x-x 0 ) ≤ 1}, given a semidefinite positive matrix P 0 ∈ R n×n consists in the convergent iteration (starting with k = 0) of the following two steps:

-Calculate a sub-gradient h k ∈ ∂ f (x k ) at the center x k of the ellipsoid, which allows for identifying a cutting plane and thus a half-space containing x k but no minimizer of f ; -Calculate and return the ellipsoid ε k+1 of minimal volume containing the half-ellipsoid

ε k ∩ {x ∈ R n : h T k (x -x k ) ≤ 0} which contains a minimizer of f .
This procedure, called in what follows ELLIPSOIDCUT can be easily modified to reduce the volume of an ellipsoid wrt. convex-constraints: when no reduction is possible, either the center is feasible, or the ellipsoid is completely empty.

As expected, since our problem is not convex, a subdivision strategy is necessary: an ellipsoid is bisected by a cutting plane (as above) function of its longest axis. The resulting Algorithm 3 keeps track of a list of ellipsoids which are guaranteed to contain all the optimal solutions.

Algorithm 3 BRANCHANDBOUNDELLIPSOID Require: Initial quadratic domain (x -o 0 ) P 0 (x -o 0 ) ≤ 1,

NUMERICAL EXAMPLES

Prototype implementations of the three algorithms were initially2 coded and tested in Matlab c R2020b. For completeness' sake, we have also included results obtained with the IBEX Optimization tool (with default parameters). Timings are provided, but since this implementation is in C++, a performance comparison is difficult to make.

First, in Section 4, two numerical examples with a limited number of conjunctions (from 1 to 3) and two maneuvers are built by using real space data. Even though the proposed numerical examples are not based on real conjunction data, two realistic operational situations are proposed by considering an actual telecom satellite as a model for the definition of the primary object.

Secondly, the efficiency and scalability of the algorithms are tested on some random academic examples of nonconvex quadratic problems (13). These tests still focus only on some expected maximum realistic operational situations: up to 8 maneuvers and less than 20 conjunctions. The details of the design of this set of random examples are given in Subsection 4.2.

Two realistic examples with 2 maneuvers

This section describes the common framework used for building two realistic encounter geometries presented in the following subsections. Specific details pertaining to a given test case, e.g. TCAs, maneuvers times or station-keeping constraints, are given in later subsections. The full orbital elements set of the model for the primary object is given in Table 1. Table 1. Orbital parameters of the primary object.

In order to set up the orbits of the secondaries, the orbit of the primary is propagated from the epoch to the TCA corresponding to the associated conjunction. The relative positions and velocities at this TCA are added to the position and velocity of this propagated orbit, creating the cartesian orbit of the associated secondary. The relative positions and velocities in the NTW frame (as it is defined in (Vallado, 2001, Section 3.3.3)) are given in Table 2 for the three secondaries considered in the following. This table contains additional data consisting in the miss-distance for each possible conjunction, the semi-major axis and inclination of each secondary. Even if the algorithms Table 2. Data for the secondaries orbits.

are not limited to two maneuvers, we chose to stick to two maneuvers, allowing a simple graphical illustration in the ∆v 1 , ∆v 2 plane of both the feasible domain and the optimal solution. For each numerical example, the two impulsive maneuvers are given along one single direction of the local NTW frame as detailed in Section 2. For each example, the parameters of the three algorithms are the same: ε vol. = 0.0001 for the ellipsoidal B. and B. algorithm, 5000 draws for the SDR method and relaxation of order 2 for the moment relaxation approach.

Test case 1: two impulsive maneuvers in the T direction and 1 TCA. The first example is essentially designed to be a proofof-concept example for the different algorithms proposed. Therefore, only one risk (a single secondary) corresponding to CDM 2 in Table 2 is considered with a TCA given at t 0 + 24 h.. The two maneuvers are respectively set to TCA minus 1 orbit (t 0 +22.33 h.) and TCA minus 0.5 orbit (t 0 +23.16 h.) and are in the Tangential direction (T) with a limit in magnitude of 15 m/s for each thrust. The station-keeping box on the positions and velocities of the primary is defined in the RTN (Radial, Transverse and Normal) or Gaussian coordinate system (see (Vallado, 2001, Section 3.3.3)) by the vector [ 5 6 5 ] km for the positions and [ 10 10 10 ] m/s for the velocities. The obtained results are summarized in Table 3 while Figure 1 exhibits the geometry of the nonconvex quadratic optimization problem in the decision variables plane as well as the behavior of the SRR algorithm. It is easily seen in Table 3 that all algorithms perform well in finding the global optimal solution in a reasonable timing except for the ellispoidal B&B algorithm which gives, though very quickly, a worst quality approximation of the optimal solution. This could be fixed by choosing a more accurate precision parameter ε vol. (0.000001 would be sufficient here) at the expense of a slight rise of the CPU time.

It is worth noticing that the relative trajectories presented in Figure 4 comply with the restrictions in all cases even if the deviation from the reference is more important in the case of the approximate optimal solution given by the ellipsoidal B&B algorithm.

Test case 2: Two impulsive maneuvers in the N direction and 3 TCAs. For this second example, a sequence of three shortterm close approaches between the primary and three secondaries (characterized by the 3 CDMs provided by CNES and given in Table 2) are defined. The three TCAs are respectively set at t 0 + 24 h., t 0 +25.5 h. and t 0 +34 h.. The two maneuvers are respectively set to TCA1 minus 1.5 orbit (t 0 +21.58 h.) and TCA1 minus 1 orbit (t 0 +22.39 h.) and are in the Normal direction (N) with a limit in magnitude of 15 m/s for each thrust.

The station-keeping box on the positions and velocities of the primary is defined in the RTN by the vector [ 1 5 6 ] km for the positions and [ 10 10 10 ] m/s for the velocities. The obtained results are summarized in ing of Table 5 with an additional difficulty for the ellipsoidal B&B algorithm when the number of nonconvex quadratic constraints increases. This feature will be more prominent when dealing with academic examples with a more important number of decision variables and nonconvex constraints. In this particular case, the worst quality approximaton of B&B algorithm does not have a significant impact on the relative trajectories. n corresponds to the number of maneuvers and respectively N quantifies the number of conjunctions. Their maximum values are set to 8 and respectively 17, with the goal of providing an insight about the expected computational efficiency and accuracy for some more involved operational examples.

Table 7 provides a synthetic overview of timings and accuracy obtained with IBEX (C++ implementation, IbexOpt version 2.8, with default parameters). It follows that for these cases, this approach is very efficient, while providing set enclosures of the optimal value.

The results obtained with the other previously presented methods (prototyped and tested with Matlab c R2020b and the Mosek optimization toolbox (version 10) software for SDP solving) are summarized in Tables 8, 9 and 10. While for lowdimensional examples, the results are very similar, when the dimension increases, the reference software IBEX is clearly more efficient. Among the prototype algorithms, the SRR method provides the fastest results, at the expense of loosing theoretical guarantees of optimality.

It is interesting to note that the B&B methods seem to perform better on these particular instances. While being preliminary, this study of computational feasibility for a practical CNES problem reveals a potential lack of numerical efficiency for the first class of SDP-based implementations. Similarly, it goes to show that the interesting theoretical properties of the ellipsoid method need to be further tuned in order to be as efficient as the more mature interval boxes representation of set-based numerics.

CONCLUSION

While several works already dealt with the optimal collision avoidance for one short-term encounter [START_REF] Bombardelli | Optimal impulsive collision avoidance in low Earth orbit[END_REF]), which can be solved by efficient convex optimization methods, the case of multiple short-term encounters boils down to a difficult non-convex NP-hard optimization problem. We formulated these problems by exploring the practical simplifications provided by CNES and observed that the indicator characterizing the collision (miss distance, Mahalanobis distance, 2D probability, maximum probability) has no direct influence on the QCLP nature of the problem obtained. As soon as at least two encounters are considered, the maneuver design is formulated as an optimization problem with linear objective and nonconvex quadratic constraints. The QCLP class of problems is then solved with several existing numerical methods whose efficiency is analyzed on realistic space examples and more academic and complex examples. These preliminary numerical analysis clearly show that the IBEX implementation is presently the most suitable and efficient tool to tackle the specific academic and realistic space examples dealt with in this article. Nevertheless, it is important to notice that the maturity of the implementation of the alternative tools based on SDP is not in any way comparable to the one of IBEX. A precise analysis of the impact of this implementation would definitely be an interesting future direction of research.

  ∆v 1 (m/s) ∆v 2 (m/s) CPU time (s) Obj. (mResults of the different algorithms for Example 1.

Figure 1 .

 1 Figure 1. Polytopic feasible region (light green with black edges), one quadratic constraint (yellow ellipse), optimal solution (red point) and feasible solutions (blue points).

Figure 2 .

 2 Figure 2. Polytopic feasible region (light green), three quadratic constraints (yellow, green and purple), optimal solution (red) and feasible solutions (blue).

4. 2

 2 Other examples with n > 2 For a fixed number n of decision variables, N quadratic constraints are randomly generated in this n-dimensional space. The N -1 ellipsoid centers are drawn from a uniform distribution on the interval [-2, 2] while the last quadratic constraint is an n-dimensional ball centered in 0 and of radius equal to 0.5. Finally, linear constraints 0 ≤ x i ≤ xi , for i = 1, • • • , n with xi uniformly distributed in the interval [1, 3] are added to the quadratic constraints. Similar to the previously presented tests, (a) and (b) Miss distance (m.) around TCA for SRR algorithm and BB algorithm (c) and (d) Relative positions and velocities trajectories from MAP reference in QSW frame for SRR and BB maneuvers.

  y) is the moment matrix and M d-ς (g y) is the localizing matrix with respect to y and g (Lasserre, 2015, Section 2.7). Obviously, f * d ≤ f * and f * d ≤ f * d+1 for all d since the set of constraints for order d + 1 is included in the one for d. One of the main features of this hierarchy of relaxations lies in the finite convergence (a finite order d exists s.t. f * d = f * ) of the hierarchy of semidefinite relaxations (15) under some classical conditions in nonlinear programming (strict complementarity condition, constraint qualification sufficient conditions, secondorder sufficient condition of optimality). Under these conditions, this finite convergence result and a certificate of global optimality are mainly entailed by Putinar's Positivstellensatz and a technical assumption on the polynomials defining the set K, which is fulfilled in our case due to the compactness of K. The algorithm based on the moments hierarchy is now summarized below. Algorithm 1 A HIERARCHY OF SDP RELAXATIONS Require: Cost function c, set K defined by quadratic and linear constraints g , = 1, • • • , r, maximum order of the relaxation d. Ensure: Optimal cost f * and global minimizer x * of Problem (13) or lower bound f * d of f * .

	Solution of SDP problem associated to relaxation of order d
	Solve the SDP Problem (15) to get f * d and y * d if it exists; if y * d does not exist then f * d is only a lower bound of the optimum:

  Nonconvex Quadratic constraints ELLIPSOIDSLIST, Linear constraints LINEARCON-STRLIST, linear objective c, subdivision accuracy ε VMin . Ensure: Minimal cost c m , value o m and list of ellipsoids containing all the global optimizers L keep .

	Start with initial ellipsoid
	L ← ∅; Push (L, (P 0 , o 0 )), c m ← +∞
	while L = ∅ do
	(P, o) ←Pop(L)
	if VOLUMEELLIPSOID(P) ≤ ε VMin then
	Push(L keep , (P, o))
	else
	Volume contractor wrt. linear constraints
	(P, o) ←ELLIPSOIDCUT((P, o), LINEARCONSTRLIST)
	if (P, o) = ∅ then
	Test center feasibility wrt. non-convex constraints
	if ISFEASIBLE(o,ELLIPSOIDLIST) then
	{Update current minimum}
	if c m > c T o then
	c m ← c T o; o m ← o
	end if
	Volume contractor wrt. objective
	(P, o) ←ELLIPSOIDCUT((P, o), c)
	Push((P, o), L)
	else
	Subdivide ellipsoid
	(P 1 , o 1 , P 2 , o 2 ) ← BISECTELLIPSOID(P, o)
	Push((P 1 , o 1 ), (P 2 , o 2 ), L)
	end if
	end if
	end if
	end while
	return (c m , o m )

Table 5 .

 5 Table 5, while Figure2exhibits the geometry of the nonconvex quadratic optimization problem in the decision variables plane as well as the behavior of the SRR algorithm. Identical conclusions may be drawn from the read-Results of the different algorithms for Example 2.

	s)

∆v 1 (m/s) ∆v 2 (m/s) CPU time (s) Obj. (m/

Table 7 .

 7 IBEX results for higher-dimensional examples.

	Case				Solution features
	#	n	N		CPU (s)	f * d	f * IBEX
	3D	3	17		0.05	0.5402	0.5414
	4D	4	10		0.12	0.5	0.5
	5D	5	10		0.23	1.42	1.4155
	8D	8	2		0.22	3.1632	3.1697
	Case				Solution features
	# 3D	n 3 17 N	L 10 4	CPU (s) 5.15	f * d 0.5402	f * SRR 0.5667
	4D	4 10	10 5	4.03	0.5	0.5732
	5D	5 10	10 5	4.8	1.142	1.7956
	8D	8	2	10 6	11173	3.1623	5.77

Table 8 .

 8 SRR results for higher-dimensional examples (L corresponds to the number of random samples).

	Case			Solution features
	#	n	N	d	CPU (s)	f * d
	3D	1 17	4	2.12	0.5402
	4D	1 10	4	4.76	0.5
	5D	1 10	4	17.64	1.142
	8D	1	2	4	1633.5	3.1623

Table 9 .

 9 Lasserre's hierarchy results for higher-dimensional examples (d corresponds to the relaxation degree) .

	Case			Solution features
	#	n	N	ε VMin	CPU (s)	f * d	f * ell
	3D	3 17	0.004	10.52	0.5402	0.74
	3D	3 17	0.0004	196.33	0.5402	0.59
	4D	4 10	0.0004	5.45	0.5	0.61
	4D	4 10 0.00004	171.6	0.5	0.54
	5D	5	5	0.004	0.4	1.142	2.55
	5D	5	5	0.0004	2.08	1.142	1.91
	5D	5	5	0.00004	55.15	1.142	1.72
	8D	8	2	0.0004	6181.3	3.16323	4.19

Table 10 .

 10 Ellipsoid B&B results for higher-dimensional examples (ε VMin is the volume threshold for processed elliposids).

A more efficient implementation is available in CNES flight dynamics Java library PATRIUS.


Fabien Halter
email: fabien.halter@cnrs-orleans.fr

Alka Karan

Guillaume Dayma

Christian Chauveau

Investigation of the Chemiluminescence Signature of Ammonia Flames

The interest in studying the chemiluminescence signature of flames lies in that it often provides information correlating the excited species to different parameters like the position of maximum HRR, blow-off characteristics, etc. For ammoniaair flames, NH2* has been determined as one of the indicators of HRR. Excited species have not been accommodated in ammonia kinetic mechanisms due to a lack of required data. The objective of this work is to provide concentration for these excited species (NH2* and OH*) to ultimately include them into existing kinetic mechanisms (not done in the current study).

Introduction

Ammonia combustion has been receiving widespread attention in the past few years. Various properties of ammonia flames have been characterized especially for blended fuels. It is important to know the position of maximum heat release rate (HRR) in gas combustors as it provides details on flame stabilization. The interest of using ammonia as a fuel in gas turbines and spark ignition engines has been well established [START_REF] Frigo | Analysis of the behaviour of a 4stroke SI engine fueled with ammonia and hydrogen[END_REF]. In combustion, chemiluminescence is often used to correlate excited species to the position of maximum HRR. Chemiluminescence is also used to predict the flame topology before the blowout [START_REF] Wiseman | A comparison of the blow-out behaviour of turbulent premixed ammonia/hydrogen/nitrogen-air and methane-air flames[END_REF], to quantify and study NO emissions [START_REF] Pugh | An investigation of ammonia primary flame combustor concepts for emissions reduction with OH*, NH2* and NH* chemiluminescence at elevated conditions[END_REF][START_REF] Zhu | NO and OH* emission characteristics of very-lean to stoichiometric ammonia-hydrogen-air swirl flames[END_REF] and to study radical formation in high-temperature reactions [START_REF] Vigueras-Zúñiga | Methane/Ammonia radical formation during high temperature reactions in swirl burners[END_REF]. Studies on excited species are widely carried out for blended mixtures like ammonia-methane or ammonia-hydrogen [START_REF] Zhu | Chemiluminescence signature of premixed ammonia-methane-air flames[END_REF][START_REF] Valera-Medina | Premixed ammonia/hydrogen swirl combustion under rich fuel conditions for gas turbines operation[END_REF]. Performing a qualitative analysis is the most common way of assessing excited species. From [START_REF] Karan | Insight into the inner structure of stretched premixed ammonia-air flames[END_REF], it is seen that the product of the concentrations of NH2* and NH* is a good indicator of the HRR for ammonia-air flames. It was also seen that the signal of NH2* is far superior to NH* and so monitoring NH2* alone is sufficient to estimate the position of HRR. Literature in ammonia combustion provides only a qualitative analysis of species and so these species have not been accommodated into the kinetic schemes due to a lack of required data. The objective of this work is to give a meticulous description of the process used to determine the concentration of NH2* and OH*. The data obtained from this work can be used to add these excited species into any of the existing kinetic mechanisms.

Method 2.1 Chemiluminescence set up

The experimental set-up used to study the chemiluminescence of the excited species -NH2* and OH* consists of a low-speed camera with high resolution (Back-illuminated sCMOS Kuro with a frame size of 2048 × 2048 pix 2 ). A UV objective lens is used to measure OH* as its emission peaks in the UV range whereas a VR (visible range) objective lens is used for NH2* as it is found in the VR of the emission spectra. A filter of interest is mounted upon the lens to record the emissions from these excited species. The wavelengths of the filters for NH2* and OH* are centered at 632 nm and 310 nm respectively each with a bandwidth of 10 nm. The filter diameter of NH2* and OH*is 50 mm and 25 mm respectively. These experiments are performed on a Bunsen burner at atmospheric conditions for premixed ammonia-air mixtures with an equivalence ratio ranging between 0.9 and 1.4. Since ammonia flames are weak in nature, a coflow of stoichiometric methane-air flame is used to stabilize the main flame. The effect of the pilot flame on the main flame has been studied by changing the power emitted by the pilot flame and it was observed that there was no direct influence of this pilot flame on the parameters measured in the main flame. The mean velocity of the fresh gases at the exit of the burner ranged from 0.18 to 0.2 m s -1 . The exposure time of the camera was set to 200 ms. An instantaneous shot is used for the study as it was seen that an average flame perturbation of over 200 images is small i.e. the induced error % in thickness is 13% for an exposure time of 200 ms. The conversion factor is 1.75 × 10 -2 mm pix -1 for NH2* and 2.45 × 10 -2 mm pix -1 for OH*.

Procedure to convert the intensity of the species to the concentration

Under this section, the procedure to evaluate the concentration from the chemiluminescence image of the species is described. This is a 4-step process that involves considering the response of the optics, followed by the response of the calibration lamp, then the procedure to correlate the intensity of pixels as seen by the camera to the concentration of the species and finally to adapt the recommended procedure to the given conditions. A region of interest (ROI) was fixed by using the longest flame possible in these experiments (at φ = 1.4).

Step 1: Optics calibration Both the camera and the filter responses are studied in this section. An LED lamp of 95% homogeneity with dimensions 172 × 132 (mm 2 ) is placed at the position of the burner. Two images of the LED lamp are captured, one with the filter mounted on and the other without the filter. On choosing a line across the images, the normalized intensities of these images are plotted as a function of the width of the image. These plots resembled gaussian-like behaviour and so these curves were fitted with the appropriate Gaussian curves and it was seen that the intensity was reduced by 2% upon the use of the filter.

To study the camera response, two methods are used. The first method is to maintain the exposure time of the camera and increase the power of the LED lamp and the other method is to fix the power of the lamp and increase the exposure time of the camera. Both the methods were deployed and it was seen that the intensity of these images (i.e. the camera response) was linear.

Step 2: Tungsten lamp calibration A quartz-tungsten-halogen lamp of power 45 W and current 6.5 A is used. The irradiation of this lamp has been provided as a function of wavelength. Since this lamp was too bright, two neutral density (ND) filters 0.8 and 1 were used. The transmission of these ND filters varies as a function of wavelength and is provided in the datasheet of these filters. As the optical configuration changed, Step 1 was repeated using the tungsten lamp as the source of light. Once again, a linear response was observed for the camera.

Step 3: Correlation between the intensity of the species and its concentration

From [START_REF] Zhao | Experimental and numerical study of OH* chemiluminescence in hydrogen diffusion flames[END_REF] and [START_REF] Liu | Quantitative Measurement of OH* and CH* Chemiluminescence in Jet Diffusion Flames[END_REF], equations correlating the intensities of the inversed Abel transformed image (𝐼 𝑂𝐻 * ,𝐴𝑏𝑒𝑙 ) and the calibration lamp (𝐼 𝑐 ) to the respective photons emitted are expressed as:

(1)

The ratio of the number of photons emitted can be directly related to the concentration of the species and is given as:

(2)

Where 𝑁 𝐴 is the Avogadro's constant, λ is the wavelength of interest, t is the exposure time, 𝜏 is the efficiency of the optics, h is the Planck's constant, ν is the frequency of light, 𝜂 is the efficiency of the filter and S is the irradiation of the lamp reported in 𝑊 𝑚 2 . µ𝑚 . 𝑠𝑟 On equating these two equations 1 and 2, the intensity of the species can be directly related to its concentration and is given as:

(3)

In equation 3, subscript 1 represents the chemiluminescence case and subscript 2 represents the calibration lamp case. 𝜆1 and 𝜆2 represent the bandwidth of the filter. Since the efficiency of the filter does not change within the bandwidth, it can be taken out of the integral and can be cancelled without the numerator. The optical efficiency for the two cases is not the same as NDs are used for the calibration lamp.

Step 4: Adapting the irradiation term to the experiment system

The irradiation of the tungsten lamp is provided as a function of wavelengths starting from 300 nm as a step of 5 nm. These values are given at a distance of 50 cm and are reported in 𝑊 𝑚 2 . 𝑛𝑚 . The distance between the burner and the filter is 32 cm and so the irradiation of interest is 32 cm. To determine this value, the reported irradiation is multiplied by the area of the sphere at 50 cm and is then divided by the area of the sphere at 32 cm. In equation 3, the required irradiation in terms of 𝑊 𝑚 2 . µ𝑚 . 𝑠𝑟 and so, this new irradiation term is divided by the solid angle. The solid angle is defined as the area of the filter divided by the square of the distance of interest. On inputting all the values into equation 3, the concentration of the species can be determined.

Results and Discussion

The concentrations of chemiluminescence of NH2* and OH* for φ=0.9-1.4 are computed. Fig. 2 and Fig. 3 represent the concentrations of NH2* and OH* at φ = 1.4. The quenching effects have not been accounted for. These flames undergo stretch effects as well and so the concentration varies as a function of stretch. However, [START_REF] Karan | Insight into the inner structure of stretched premixed ammonia-air flames[END_REF] pointed out that the stretch values are quite low under these conditions.

Concluding Remarks

The concentrations for two excited species: NH2* and OH* for equivalence ratios ranging between 0.9 and 1.4 are provided from the chemiluminescence of ammoniaair flames. The procedure used to determine the species has been described and was validated using the OH* and CH* concentrations for methane-air flames provided in the literature [START_REF] Panoutsos | Numerical evaluation of equivalence ratio measurement using OH * and CH * chemiluminescence in premixed and non-premixed methane-air flames[END_REF] and by running a Chemkin simulation using the GRI 3.0 Mech. The quenching effects and the stretch effects have not been accounted for. These concentrations can be used to add these excited species to the existing kinetic mechanisms for ammonia combustion so that different parameters like HRR can be monitored in numerical simulations by following the profiles of NH2*.

Fig. 1 :

 1 Fig. 1: Schematic diagram of the chemiluminescence setup.These experiments are performed on a Bunsen burner at atmospheric conditions for premixed ammonia-air mixtures with an equivalence ratio ranging between 0.9 and 1.4. Since ammonia flames are weak in nature, a coflow of stoichiometric methane-air flame is used to stabilize the main flame. The effect of the pilot flame on the main flame has been studied by changing the power emitted by the pilot flame and it was observed that there was no direct influence of this pilot flame on the parameters measured in the main flame. The mean velocity of the fresh gases at the exit of the burner ranged from 0.18 to 0.2 m s -1 . The exposure time of the camera was set to 200 ms. An instantaneous shot is used for the study as it was seen that an average flame perturbation of over 200 images is small i.e. the induced error % in thickness is 13% for an exposure time of 200 ms. The conversion factor is 1.75 × 10 -2 mm pix -1 for NH2* and 2.45 × 10 -2 mm pix -1 for OH*.

Fig. 2 :

 2 Fig. 2: Concentration of NH2* at ϕ= 1.4. The colour bar indicates the concentration in mol m -3 .

Fig. 3 :

 3 Fig. 3: Concentration of OH* at ϕ= 1.4. The colour bar indicates the concentration in mol m -3 .