Jérôme Dedecker

Florence Merlevède

Almost sure invariance principle for the Kantorovich distance between the empirical and the marginal distributions of strong mixing sequences

Keywords: . Empirical process, Wasserstein distance, Almost sure invariance principle, Compact law of the iterated logarithm, Bounded law of the iterated logarithm, Conditional Value at Risk Mathematics Subject Classification (2010). 60F15, 60G10, 60B12

We prove a strong invariance principle for the Kantorovich distance between the empirical distribution and the marginal distribution of stationary α-mixing sequences.

Running head. ASIP for the empirical W1 distance.

Introduction and notations

Let (Xi) i∈Z be a strictly stationary sequence of real-valued random variables. Define the two σ-algebras F0 = σ(Xi, i ≤ 0) and G k = σ(Xi, i ≥ k), and recall that the strong mixing coefficients (α(k)) k≥0 of Rosenblatt [START_REF] Rosenblatt | A central limit theorem and a strong mixing condition[END_REF] be the empirical measure based on X1, . . . , Xn. In this paper, we prove a strong invariance principle for the Kantorovich distance W1(µn, µ) between µn and µ under a condition on the mixing coefficients α(k). Recall that the Kantorovich distance (also called Wasserstein distance of order 1) between two probability measures µ and ν is defined by

W1(µ, ν) = inf π∈M (µ,ν)
|x -y|π(dx, dy) ,

where M (µ, ν) is the set of probability measures on R 2 with marginals µ and ν. We shall use the following well known representation for probabilities on the real line:

W1(µ, ν) = |Fµ(x) -Fν (x)|dx , (1.2)
where Fµ is the cumulative distribution function of µ. Let H : t → P([X0| > t) be the tail function of |X0|. In the case where (Xi) i∈Z is a sequence of independent and identically distributed (i.i.d.) random variables, del Barrio et al. [START_REF] Del Barrio | Central limit theorems for the Wasserstein distance between the empirical and the true distributions[END_REF] used the representation (1.2) and a general result of Jain [START_REF] Jain | Central limit theorems and related questions in Banach space[END_REF] for Banach-valued random variables to prove a central limit theorem for √ nW1(µn, µ). More precisely, they showed that √ nW1(µn, µ) converges in distribution to the L1(dt) norm of an L1(dt)-valued Gaussian random variable, provided that

∞ 0 H(t) dt < ∞ . (1.3)
They also proved that √ nW1(µn, µ) is stochastically bounded iff (1.3) holds, proving that this condition is necessary and sufficient for the weak convergence of √ nW1(µn, µ). Still in the i.i.d. case, we easily deduce from Chapters 8 and 10 in Ledoux and Talagrand [START_REF] Ledoux | Probability in Banach spaces. Isoperimetry and processes[END_REF] that: if (1.3) holds, then the sequence

√ n √ 2 log log n W1(µn, µ) (1.4)
satisfies a compact law of the iterated logarithm.

For strongly mixing sequences in the sense of Rosenblatt [START_REF] Rosenblatt | A central limit theorem and a strong mixing condition[END_REF], we proved in [START_REF] Dedecker | Behavior of the Wasserstein distance between the empirical and the marginal distributions of stationary α-dependent sequences[END_REF] the central limit theorem for √ nW1(µn, µ) under the condition

∞ 0 ∞ k=0 (α(k) ∧ H(t)) dt < ∞ (1.5)
(where a ∧ b means the minimum between two reals a and b), and we give sufficient conditions for (1.5) to hold. Note that, in [START_REF] Dedecker | Behavior of the Wasserstein distance between the empirical and the marginal distributions of stationary α-dependent sequences[END_REF], we used a weaker version of the α-mixing coefficients, that enables to deal with a large class of non-mixing processes in the sense of Rosenblatt [START_REF] Rosenblatt | A central limit theorem and a strong mixing condition[END_REF].

In Section 2 of this paper, we prove a strong invariance principle for W1(µn, µ) under the condition (1.5). The compact law of the iterated logarithm for (1.4) easily follows from this strong invariance principle. In Section 3, we apply our general result to derive the almost sure rate of convergence of the empirical estimator of the Conditional Value at Risk (CV aR) for stationary α-mixing sequences.

In the rest of the paper, we shall use the following notation: for two sequences (an) n≥1 and (bn) n≥1 of positive reals, an bn means there exists a positive constant C not depending on n such that an ≤ Cbn for any n ≥ 1.

Main result

Our main result is the following strong invariance principle for W1(µn, µ). Theorem 2.1. Assume that (1.5) is satisfied. Then, enlarging the probability space if necessary, there exists a sequence of i.i.d. L1(dt)-valued centered Gaussian random variables (Zi) i≥1 with covariance function defined as follows: for any f, g ∈ L∞(dt),

Γ(f, g) = Cov f (t)Z1(t) dt, g(t)Z1(t) dt = k∈Z f (t)g(s)Cov(1 X 0 ≤t , 1 X k ≤s ) ds dt , (2.1
) and such that nW1(µn, µ) - Remark 2.2. In [START_REF] Cuny | Invariance principles under the Maxwell-Woodroofe condition in Banach spaces[END_REF], Cuny proved a strong invariance principle for W1(µn, µ). under the condition

∞ k=0 1 √ k + 1 ∞ 0 α(k) ∧ H(t) dt < ∞ (2.2)
(in fact, he proved the result for a weaker version of the α-mixing coefficient, the same as that used in [START_REF] Dedecker | Behavior of the Wasserstein distance between the empirical and the marginal distributions of stationary α-dependent sequences[END_REF] for the central limit theorem). It follows from Section 5 of [START_REF] Dedecker | Behavior of the Wasserstein distance between the empirical and the marginal distributions of stationary α-dependent sequences[END_REF], that the condition (1.5) is always less restrictive than (2.2). As a consequence of Theorem 2.1, we get the compact law of the iterated logarithm. Let K be the unit ball of the reproducing kernel Hilbert space (RKHS) associated with Γ, and C be the image of K by the L1(dt) norm. The following corollary holds: Corollary 2.1. Assume that (1.5) is satisfied. Then the sequence

√ n √ 2 log log n W1(µn, µ)
is almost surely relatively compact, with limit set C. The proof of Theorem 2.1 is based on two ingredients: a martingale approximation in L1(dt), as in [START_REF] Dedecker | Behavior of the Wasserstein distance between the empirical and the marginal distributions of stationary α-dependent sequences[END_REF], and the following version of the bounded law of the iterated logarithm, which has an interest in itself. Proposition 2.1. Assume that (1.5) holds, and let

V = ∞ 0 ∞ k=0 (α(k) ∧ H(t)) dt .
(2.3)

Then, there exists a universal constant η such that for any ε > 0, Combining the arguments in [START_REF] Berthet | Central limit theorem and almost sure results for bivariate empirical W1 distances[END_REF] and the proof of Theorem 2.1, one can prove the following strong invariance principle for n (W1(µn,X , µn,Y ) -W1(µX , µY )).

n≥2 1 n P max 1≤k≤n kW1(µ k , µ) > (ηV + ε) n log log n < ∞ . ( 2
Let ϕ be the continuous function from L1(dt) to R defined by

ϕ(x) = sign{FX (t) -FY (t)} x(t)1 F X (t) =F Y (t) + |x(t)|1 F X (t)=F Y (t) dt ,
where FX (resp. FY ) is the cumulative distribution function of µX (resp. µY ). Assume that

∞ 0 ∞ k=0 (α(k) ∧ HX (t)) dt < ∞ and ∞ 0 ∞ k=0 (α(k) ∧ HY (t)) dt < ∞ .
Then, enlarging the probability space if necessary, there exists a sequence of i.i.d. L1(dt)-valued centered Gaussian random variables (Zi) i≥1 with covariance function given by: for any f, g ∈ L∞(dt),

Γ(f, g) = Cov f (t)Z1(t) dt, g(t)Z1(t) dt = k∈Z f (t)g(s)Cov(1 X 0 ≤t -1 Y 0 ≤t , 1 X k ≤s -1 Y k ≤s ) ds dt ,
and such that

n (W1(µn,X , µn,Y ) -W1(µX , µY )) -ϕ n k=1 Z k = o( n log log n) almost surely.
3 Rates of convergence of the empirical estimator of the Conditional Value at Risk

The Conditional Value at Risk at level u ∈ (0, 1] of a real-valued integrable random variable X (CV aRu(X)) is a "risk measure" (according to the definition of Acerbi and Tasche [START_REF] Acerbi | On the coherence of Expected Shortfall[END_REF]), which is widely used in mathematical finance. It is sometimes called Expected Shortfall of Average Value at Risk. We refer to the paper [START_REF] Acerbi | On the coherence of Expected Shortfall[END_REF] for a clear definition of that indicator, and for its relation with other well known measures, such as the Value at Risk, the Worst Conditional Expectation, the Tail Conditional Expectation... According to Acerbi and Tasche [START_REF] Acerbi | On the coherence of Expected Shortfall[END_REF], CV aRu(X) can be expressed as

CV aRu(X) = - 1 u u 0 F -1 X (x)dx ,
where FX is the cumulative distribution function of the variable X, and F -1 X is its usual cadlag inverse:

F -1 X (u) = inf{x ∈ R : FX (x) ≥ u}.
Concerning the difference between the Conditional Value at Risk of two random variables X and Y , the following elementary inequality holds (see for instance [START_REF] Rio | About the conditional value at risk of partial sums[END_REF]):

|CV aRu(X) -CV aRu(Y )| ≤ 1 u 1 0 |F -1 X (x) -F -1 Y (x)|dx = 1 u W1(µX , µY ) , (3.1)
where µX (resp. µY ) is the distribution of X (resp. Y ).

Consider now the problem of estimating CV aRu(X) from the random variables X1, ..., Xn, where (Xi) i∈Z is a stationary sequence of α-mixing random variables with common distribution µ = µX . A natural estimator is then

CV aRu,n = - 1 u u 0 F -1 n (x)dx ,
where Fn is the empirical distribution function based on X1, . . . , Xn. From (3.1), we get the upper bound [START_REF] Ledoux | Probability in Banach spaces. Isoperimetry and processes[END_REF]) that the constant κ(Γ) can be expressed as Let (Ω, A, P) be the underlying probability space. By a standard argument, one may assume that Xi = X0 • T , where T : Ω → Ω is a bijective, bi-measurable transformation, preserving the probability P.

CV aRu(X) -CV aRu,n ≤ 1 u 1 0 |F -1 X (x) -F -1 n (x)|dx = 1 u W1(
κ(Γ) = sup f : f ∞≤1 Var f (t)Z(t)dt 1/2 ≤ |Z(t)|dt 4
Let also Fi = σ(X k , k ≤ i). Let Y0(t) = 1 X 0 ≤t -F (t), and Y k (t) = Y0(t) • T k = 1 X k ≤t -F (t)
. With these notations and the representation (1.2) one has that

nW1(µn, µ) = n k=1 Y k (t) dt . (4.1)
From Section 4 in [START_REF] Dedecker | Behavior of the Wasserstein distance between the empirical and the marginal distributions of stationary α-dependent sequences[END_REF], we know that, if (1.5) holds, then

Y0(t) = D0(t) + A(t) -A(t) • T, (4.2)
where D0 is such that E(D1(t)|F-1) = 0 almost surely and D0(t) 2 dt < ∞, and A is such that A(t) 1 dt < ∞. Moreover, the covariance operator of D0 is exactly Γ: for any f, g ∈ L∞(dt),

Γ(f, g) = Cov f (t)D0(t) dt, g(t)D0(t) dt . (4.3) Let D k (t) = D0 • T k . From (4.2), it follows that n k=1 Y k = n k=1 D k + A • T -A • T n . (4.4)
From [4, Proposition 3.3], we know that, enlarging the probability space if necessary, there exists a sequence of i.i.d. L1(dt)-valued centered Gaussian random variables (Zi) i≥1 with covariance function Γ such that Now, as will be clear from the proof, Proposition 2.1 also holds on the space L 1 ([-M, M ] c , dt), and implies that there exists a universal constant η such that, for any positive ε,

lim sup n→∞ 1 √ n log log n [-M,M ] c n k=1 Y k (t) dt ≤ ε + η ∞ M ∞ k=0
min {α(k), H(t)} dt almost surely. Hence the proof of (4.6) will be complete if we prove that, for any M > 0, lim sup

n→∞ 1 √ n log log n M -M |A(t) • T n | dt = 0 almost surely. (4.9)
To prove (4.9), we work in the space H = L2([-M, M ], dt), and we denote by • H and •, • the usual norm and scalar product on H. Since E( D0 2 H ) < ∞, we know from [START_REF] Cuny | Invariance principles under the Maxwell-Woodroofe condition in Banach spaces[END_REF] that n k=1 D k satisfies the compact law of the iterated logarithm in H. Since k≥0 α(k) < ∞ and Y0 is bounded in H, we infer from [START_REF] Dedecker | On the almost sure invariance principle for stationary sequences of Hilbert-valued random variables[END_REF] that n k=1 Y k satisfies also the compact law of the iterated logarithm in H. Now, arguing exactly as in the end of the proof of [START_REF] Dedecker | On the almost sure invariance principle for stationary sequences of Hilbert-valued random variables[END_REF]Theorem 4], one has: for any f in H

Proof of Proposition 2.1

For any n ∈ N, let us introduce the following notations:

R(u) = min{q ∈ N * : α(q) ≤ u}Q(u) and R -1 (x) = inf{u ∈ [0, 1] : R(u) ≤ x} .
For a positive real a that will be specified later, let

mn = a n log log n , vn = R -1 (mn) , Mn = Q(vn) . (4.14)
For any M > 0, let gM (y) = (y ∧ M ) ∨ (-M ). For any integer i, define

X i = gM n (Xi) and X i = Xi -X i . (4.15)
We first recall that, by the dual expression of W1(µn, µ), nW1(µn, µ) = sup

f ∈Λ 1 n i=1 (f (Xi) -E(f (Xi))) .
where Λ1 is the set of Lipschitz functions such that |f (x) -f (y)| ≤ |x -y|. Hence,

nW1(µn, µ) ≤ sup f ∈Λ 1 n i=1 f (X i ) -E(f (X i )) + sup f ∈Λ 1 n i=1 f (Xi) -f (X i ) -E(f (Xi) -f (X i )) .
Therefore, setting,

F n (t) = 1 n n k=1 1 {X k ≤t} and F (t) = P(X 1 ≤ t) ,
and noticing that

k F k -F 1 = sup f ∈Λ 1 k i=1 f (X i ) -E(f (X i ) , we get max 1≤k≤n kW1(µ k , µ) ≤ max 1≤k≤n k F k -F 1 + n i=1 (|X i | + E(|X i |) . (4.16)
Now, note that

n≥2 1 √ n log log n E(|X n |) ≤ n≥2 1 √ n log log n +∞ 0 P |X0|1 |X 0 |>Q(vn) > t dt ≤ n≥2 1 √ n log log n +∞ Q(vn) H(t)dt ≤ n≥2 1 √ n log log n vn 0 Q(u)du ≤ n≥2 1 √ n log log n 1 0 Q(u)1 mn≤R(u) du 1 0 R(u)Q(u)du .
But, according to Propositions 5.1 and 5.2 in [START_REF] Dedecker | Behavior of the Wasserstein distance between the empirical and the marginal distributions of stationary α-dependent sequences[END_REF], condition (1.5) implies that

1 0 R(u)Q(u)du < ∞ . (4.17)
Hence, to prove (2.4) it suffices to show that there exists an universal constant η such that for any ε > 0, Taking into account the definition of mn, it follows that there exists n0 depending on a, V and η, such that for any n ≥ n0, 8mn ≤ κV √ n log log n. This proves the proposition in the case where qn = n.

From now on, we assume that qn < n. Therefore qn = min{k ∈ N * : α(k) ≤ vn} and then α(qn) ≤ vn. For any integer i, define

Ui(t) = iqn k=(i-1)qn+1 1 X k ≤t -E 1 X k ≤t .
k F k -F 1 ≥ 4V κ n log log n ≤ I1(n) + I2(n) + I3(n) , (4.24)
where

I1(n) = P   Mn -Mn [n/qn] i=1 |Ui(t) -U * i (t)| dt ≥ V κ n log log n   I2(n) = P Mn -Mn max 2≤2j≤[n/qn] j i=1 U * 2i (t) dt ≥ V κ n log log n I3(n) = P Mn -Mn max 1≤2j-1≤[n/qn] j i=1 U * 2i-1 (t) dt ≥ V κ n log log n .
Using Markov's inequality and (4.22), we get

I1(n) n √ n log log n Mnα(qn) n √ n log log n vnQ(vn) n √ n log log n R -1 (mn) 0 Q(u)du .
Hence, by (4.17),

n≥2 1 n I1(n) n≥2 1 √ n log log n R -1 (mn) 0 Q(u)du 1 0 R(u)Q(u)du < ∞ .
To handle now the term I2(n) (as well as I3(n)) in the decomposition (4.24), we shall use again Markov's inequality but this time at the order p ≥ 2. Hence for p ≥ 2, taking into account the stationarity, we get Note that u < H(t) ⇐⇒ t < Q(u). Consequently u < H(t) implies that Q -2 (u) < t -2 . Hence

I2(n) ≤ 1 (V κ) p (n log log n) p/2   Q(vn) -Q(vn) max 2≤2j≤[n/qn] j i=1 Ũ2i (

  are defined by α(k) = sup A∈F 0 ,B∈G k |P(A ∩ B) -P(A)P(B)| . (1.1)Let µ be the common distribution of the Xi's, and let

Z

  k (t) dt = o( n log log n) almost surely.

. 4 )

 4 Remark 2.3. (The bivariate case). Let (Xi, Yi) i∈Z be a stationary sequence of R 2 -valued random variables, and define the coefficients α(k) as in (1.1), with the two σ-algebras F0 = σ(Xi, Yi, i ≤ 0) and G k = σ(Xi, Yi, i ≥ k). Let µX (resp. µY ) be the common distribution of the Xi's (resp. the Yi's), and let

Proofs 4 . 1

 41 Proof of Theorem 2.1

1 √ 1 √

 11 result will follow from (4.1), (4.4) and (4.5) if we can prove that lim n→∞ n log log n |A(t) • T n | dt = 0 almost surely. (4.6) To prove (4.6), we start by considering the integral over [-M, M ] c , for M > 0. Applying again [4, Proposition 3.3], we infer that lim sup n→∞ 2n log log n [-M,M ] c n k=1 D k (t) dt ≤ [-M,M ] c D0(t) 2 dt almost surely. (4.7)

  n [-M,M ] c |A(t) • T n | dt = 0 almost surely.

  ) i≥1 be a complete orthonormal basis of H and PN (f ) = N k=1 f, e k e k be the projection of f on the space spanned by the first N elements of the basis. From (4.10), we get thatlim n→∞ PN (A • T n ) √ n log log n = 0 almost surely. (4.11) On another hand, applying again [4, Proposition 3.3] (as done in (4.7)), we get lim PN )A • T n H √ n log log n = 0 almost surely, which, together with (4.11), implies (4.9). The proof of Theorem 2.1 is complete.

  k -F 1 > ηV n log log n < ∞ . (4.18) For this purpose, let qn = min{k ∈ N * : α(k) ≤ vn} ∧ n . (4.19)Since R is right continuous, we have R(R -1 (w)) ≤ w for any w, henceqnMn = R(vn) = R(R -1 (mn)) ≤ mn . (4.20) Assume first that qn = n. Bounding f (X i ) -E(f (X i )) by 2Mn, we obtain max 1≤k≤n k F k -F 1 ≤ 2nMn = 2qnMn ≤ 2mn . (4.21)

( 4 . 2 Qc11

 42 25)Using similar arguments as to handle the quantity I2(n) in the proof of[START_REF] Dedecker | Behavior of the Wasserstein distance between the empirical and the marginal distributions of stationary α-dependent sequences[END_REF] Proposition 3.4], we have κ) p (n log log n) p/nκ p (log log n) p/2 .Let nowp = pn = max{c log log n, 2}, where c will be specified later. Set n1 = min{n ≥ 2 : c log log n ≥ 2}. It follows thatn≥n 1 1 n(V κ) p (n log log n) p/2 Q(vn) -Q(vn) J1(t) 1/p dt p log log n ,which is finite provided we take κ such that 2c 1 √ 2c κ = α -1 with α > 1 and c > (log α) -1 . On another hand, proceeding as in (4.26), we deduce that, for any t > 0, {X i ≤t} -P(X i ≤ t) {X i ≤t} -P(X i ≤ t)

- 2 )

 2 -1 Q(vn) 1-2/p .

  Let kn = [n/qn]. For any t, applying Rio's coupling lemma (see[START_REF] Rio | Théorie asymptotique des processus aléatoires faiblement dépendants[END_REF] Lemma 5.2]) recursively, we can construct random variables (U * i (t)) 1≤i≤kn such that • U * i (t) has the same distribution as U i for all 1 ≤ i ≤ kn,

	and notice that
							Mn	j
		max 1≤k≤n	k F k -F 1 ≤ 2qnMn +	-Mn	1≤j≤[n/qn] max	i=1	Ui(t) dt .
	• the random variables (U * 2i (t)) 2≤2i≤kn are independent, as well as the random variables
	(U * 2i-1 (t)) 1≤2i-1≤kn ,
	• we can suitably control Ui(t) -U * i (t) 1 as follows: for any i ≥ 1,
				Ui(t) -U * i (t) 1 ≤ 4qnα(qn) .	(4.22)
	Substituting U * i (t) to Ui(t), we obtain
						j
	max 1≤k≤n	k F k -F 1 ≤ 2qnMn +	max 2≤2j≤[n/qn]	i=1	U * 2i (t)
							j	[n/qn]
				+	max 1≤2j-1≤[n/qn]	i=1	U * 2i-1 (t) +	i=1	|Ui(t) -U *

i (t)| . (4.23) Therefore, setting κ = η/4, for n ≥ n0, P max 1≤k≤n

  t) Applying Rosenthal's inequality (see for instance[START_REF] Pinelis | Optimum bounds for the distributions of martingales in Banach spaces[END_REF] Theorem 4.1]) and taking into account the stationarity, there exist two positive universal constants c1 and c2 not depending on p such that

		j		p
	max 2≤2j≤[n/qn]	i=1	U * 2i (t)	p	≤ c p 1 p p/2 (n/qn) p/2 U2(t) p 2 + c p 2 p p (n/qn) U2(t) p p := J1(t) + J2(t) .
						 p
					dt		.
					p

, where Z is an L1(dt)-valued centered random variable with covariance function Γ.

Set n2 = min{n ≥ 2 : c log log n ≥ 4}. It follows that

Hence, since qnMn ≤ mn, we get

which is finite by taking into account (4.17), and if we choose a = (c1κV )/(2c2

This ends the proof of the proposition.


Félix Gontier
email: felix.gontier@inria.fr

Romain Serizel

Christophe Cerisara
email: christophe.cerisara@loria.fr

SPICE+: EVALUATION OF AUTOMATIC AUDIO CAPTIONING SYSTEMS WITH PRE-TRAINED LANGUAGE MODELS

Keywords: Audio captioning, Evaluation, DCASE

Audio captioning aims at describing acoustic scenes with natural language. Systems are currently evaluated by image captioning metrics CIDEr and SPICE. However, recent studies have highlighted a poor correlation of these metrics with human assessments. In this paper, we propose SPICE+, a modification of SPICE that improves caption annotation and comparison with pre-trained language models. The metric parses captions to semantic graphs with a deep dependency annotation model and a refined set of linguistic rules, then compares sentence embeddings of candidate and reference semantic elements. We formulate a score for general-purpose captioning evaluation, that can be tailored to more specific applications. Combined with fluency error detection, the metric achieves competitive performance on the FENSE benchmark, with 84.0% accuracy on AudioCaps and 74.1% on Clotho. Further experiments show that the metric behaves similarly to the full sentence embedding similarity, while the decomposition into semantic elements allows better interpretability of scores and can provide additional information on the properties of captioning systems.

INTRODUCTION

Automatic audio captioning (AAC) is the process of describing acoustic scenes with fluent sentences [START_REF] Drossos | Automated audio captioning with recurrent neural networks[END_REF]. The task has received increasing attention in recent years, in particular through a dedicated DCASE challenge task 1 . However, the adequacy of current metrics for evaluating AAC systems is an prominent subject of discussion within the community.

Current evaluation metrics are borrowed from the image captioning field. The first metrics, BLEU [START_REF] Papineni | Bleu: a method for automatic evaluation of machine translation[END_REF], ROUGE [START_REF] Lin | Automatic evaluation of machine translation quality using longest common subsequence and skip-bigram statistics[END_REF], and METEOR [START_REF] Denkowski | Meteor universal: Language specific translation evaluation for any target language[END_REF], originate in machine translation and summarization tasks, and measure the n-gram overlap between candidate and reference sentences. Several studies have subsequently revealed poor correlation of these metrics with human assessments of image caption quality [START_REF] Hodosh | Framing image description as a ranking task: Data, models and evaluation metrics[END_REF][START_REF] Elliott | Comparing automatic evaluation measures for image description[END_REF][START_REF] Kulkarni | Baby talk: Understanding and generating simple image descriptions[END_REF]. As a re-This work was funded under the ANR project LEAUDS (Grant No. ANR-18-CE23-0020). 1 dcase.community/ sult, captioning-specific alternatives CIDEr [START_REF] Vedantam | Cider: Consensus-based image description evaluation[END_REF] and SPICE [START_REF] Anderson | Spice: Semantic propositional image caption evaluation[END_REF] were proposed. CIDEr retains the approach of n-gram comparison, but integrates reference consensus information and Term Frequency Inverse Document Frequency (TF-IDF) to focus on relevant terms in the description. SPICE instead parses reference captions to a graph containing semantic elements, their attributes, and relations to one another, and evaluates the candidate graph via synonym lemma matching. SPI-DEr [START_REF] Liu | Improved image captioning via policy gradient optimization of spider[END_REF], the combination of SPICE and CIDEr, was found to outperform both in correlation to human assessments, and is now the main metric for evaluating AAC systems. The properties indicative of acceptable captions are not well defined in the current formulation of the AAC task. Examples from Clotho [START_REF] Drossos | Clotho: An audio captioning dataset[END_REF] and AudioCaps [START_REF] Kim | AudioCaps: Generating captions for audios in the wild[END_REF] associate several human reference captions to examples, each describing events relevant to the annotator instead of the entire scene. Likewise, the same acoustic content can be described with different sentence structure, which reference captions do not cover exhaustively. Evaluation through n-gram overlap, suited for language generation with specific wording, thus seems less relevant for audio captioning. Labbé et al. [START_REF] Labbé | Is my automatic captioning system so bad? spider-max: a metric to consider several caption candidates[END_REF] further note that SPIDEr shows high variability to small differences in caption formulation. As a workaround, they propose SPIDEr-max, where the maximum SPIDEr score out of several candidates produced by the system is retained.

In [START_REF] Zhou | Can audio captions be evaluated with image caption metrics?[END_REF], the authors propose a benchmark to evaluate the accuracy of metrics in choosing the best caption in candidate pairs, against human preference. Most n-gram metrics as well as SPICE perform poorly, cementing the need for further research in caption evaluation. The authors propose FENSE, a metric based on Sentence-BERT [START_REF] Reimers | Sentence-bert: Sentence embeddings using siamese bert-networks[END_REF] embeddings similarity and penalized by fluency error detection. Although this method achieves high correlation with human judgements, Martín-Morató et al. [START_REF] Martín-Morató | A summarization approach to evaluating audio captioning[END_REF] argue that sentence-wise semantic similarity does not necessarily reflect the correspondence of described sounds. Instead, they present CBscore, a summarization-oriented metric that extracts AudioSet tags [START_REF] Gemmeke | Audio set: An ontology and human-labeled dataset for audio events[END_REF] from captions, then derives a relevance-weighted precision score from the consensus in the reference set. However, the approach is limited to sound classes, and discards their qualifiers and relationships in the description.

In this paper, we propose to combine the semantic graph parsing paradigm of SPICE with the strengths of deep language models. The metric uses a pre-trained language model to produce caption annotations, which are then parsed to semantic graphs with a broad set of linguistic rules. Candidate and reference graphs are compared through embedding similarity. The contributions of this paper are: (i) we identify limitations in the current implementation of SPICE, namely that semantic graphs are incomplete, and candidate evaluation through synonym matching is too restrictive, (ii) we replace components of SPICE with pre-trained models to produce comprehensive semantic representations, and (iii) we propose a score based on the embedding similarity between candidate and reference sub-graphs, that is suitable in the current general task setting but can be adapted to specific sub-tasks. 2

PROPOSED METHOD

The proposed metric evaluates a candidate caption against an arbitrarily large set of references, following the same process as SPICE [START_REF] Anderson | Spice: Semantic propositional image caption evaluation[END_REF]. First, the caption is annotated with partof-speech tagging, lemmatization, and dependency parsing, which produces a tree-like representation of grammatical relations in the sentence. The annotation is then parsed into semantic graphs with a set of linguistic rules. Graphs are composed of nodes (nouns), attributes (adjective, verbs) and relations. SPICE views graphs as sets of node, node-attribute, and node-relation-node tuples. The candidate graph is evaluated against a single graph comprising all reference captions by comparing each tuple pair independently. SPICE+ alters the three processing steps, by performing dependency annotation with a pre-trained language model, adding new linguistic rules and post-processing to semantic graph construction, and evaluating candidates on deep embedding comparison 3 . [START_REF] Nivre | Universal Dependencies v2: An evergrowing multilingual treebank collection[END_REF] using the Stanford probabilistic context-free grammar (PCFG) parser [START_REF] Klein | Accurate unlexicalized parsing[END_REF]. From empirical observations on Clotho examples, this parser has a high error rate, which propagates to poor rule-based parsing and leads to incomplete or incorrect semantic graphs. Since then, better parsers have been developed in the natural language processing community [START_REF] Zeman | CoNLL 2018 shared task: Multilingual parsing from raw text to Universal Dependencies[END_REF]. Thus, SPICE+ replaces the original PCFG parser with a language model trained on Universal Dependencies parsing.

Semantic graph parsing

SPICE annotates captions on Universal Dependencies

A set of linguistic rules is then matched to the dependency annotation to obtain a semantic graph. Five rules are added compared to SPICE to handle formulations often found in audio captions, such as open clausal complements (e.g. "rain continues to pour"), passive nominal subjects (e.g. "a radio is being played"), and clausal noun modifiers (e.g. "peo- 2 Code is available at: github.com/felixgontier/spice-audio. 3 Examples on the companion page: felixgontier.github.io/spice-audio ple talking"). In addition, two post-processing steps are applied on semantic graphs, that led to better results in empirical experiments. First, coordinating conjunctions systematically duplicate attributes and relations. For instance, "men and women are talking and laughing" is parsed to two nodes man and woman, each with both attributes talk and laugh. Lastly, for expressions in the form "A of B" (e.g. "group of people", "body of water") attributes and relations are associated to node B instead of node A.

Candidate evaluation

The original implementation of SPICE compares candidate and reference graphs as sets of tuples through synonym comparison. The process outputs a binary matrix, from which true positives, false positives, and false negatives are derive to compute an F-score. This approach discards related but nonsynonymous words, which results in a prominent spike at 0 in score distributions [START_REF] Martín-Morató | A summarization approach to evaluating audio captioning[END_REF]. At the same time, in audio captioning there may be several valid ways to describe a sound, e.g. a metal tube clanging, something hitting a pipe, that are not exhaustively represented in the reference set.

SPICE+ retains a tuple comparison approach, but replaces binary synonym matching with a continuous similarity measure. First, a pre-trained language model extracts sentence embeddings for each tuple in the candidate and reference graphs. For attribute or relation tuples, the element lemmas are concatenated and fed directly to the model. Then, each candidate embedding X i is compared to all reference embeddings Y j through cosine similarity, yielding a similarity matrix S i,j . Pseudo-precision and pseudo-recall scores are obtained by averaging maximum similarities along i and j:

P r = 1 N i Ni i=1 max j S i,j , Re = 1 N j Nj j=1 max i S i,j (1)
The maximum operation is applied because an element of the candidate is not expected to match all elements of the reference set. The final metric is the pseudo F-score, which is referred as SPICE+emb in experiments.

Hyper-parameters

Several parameters of the metric can be set to match specific application needs. First, tuples can be weighted in eq. ( 1), either to account for consensus in the references (e.g. with a relevance measure [START_REF] Martín-Morató | A summarization approach to evaluating audio captioning[END_REF]), or with a pre-defined term-importance dictionary depending on the task focus. The balance between precision and recall in the final score can also be changed to penalize false positives or false negatives respectively. Lastly, a decision on matches can be made by applying a threshold on S. This sets a tolerance on lexical specificity: a high threshold only matches synonymous or identical terms, whereas a lower threshold accepts different but related concepts. Table 1. Metric performance on the FENSE benchmark, with 95% confidence intervals and highest accuracy in bold.

EXPERIMENTS

Experimental setup

The FENSE benchmark [START_REF] Zhou | Can audio captions be evaluated with image caption metrics?[END_REF] is used as the main evaluation tool for the proposed metric. The benchmark evaluates the ability of metrics to assess the quality of one candidate caption against another in four cases: human-correct (HC) where both captions are references for the correct example, humanincorrect (HI) where one caption is a reference for a different example, human-machine (HM) where one caption is from a captioning system, and machine-machine (MM) where both captions are produced by a trained system. 250 pairs (921 and 1000 for MM resp.) are formed on AudioCaps and Clotho. Four human annotations of relative quality are collected per pair, and metric performance is given by its accuracy in attributing higher scores to the preferred caption. Note that the benchmark removes all punctuation from Clotho captions, which increases the error rate of dependency annotation. Experiments consider the general captioning task setting of current datasets. As such, the evaluated metric is the pseudo F-score with default precision-recall balance. Furthermore, no decision threshold or tuple weighting is applied to the similarity matrix S. The dependency and part-of-speech annotations are done with UDify [START_REF] Kondratyuk | 75 languages, 1 model: Parsing universal dependencies universally[END_REF]. The Stanford parser [START_REF] Klein | Accurate unlexicalized parsing[END_REF] is retained for lemma annotations following empirical experimentation. Tuple embeddings are computed with Sentence-BERT [START_REF] Reimers | Sentence-bert: Sentence embeddings using siamese bert-networks[END_REF] (checkpoint paraphrase-TinyBERT-L6-v2) to follow FENSE [START_REF] Zhou | Can audio captions be evaluated with image caption metrics?[END_REF]. However, similar results are obtained using MPNet [START_REF] Song | Mpnet: Masked and permuted pre-training for language understanding[END_REF].

Ranking capabilities

The ranking performance of SPICE+ is first evaluated on the FENSE benchmark against current metrics CIDEr, SPICE, and SPIDEr, as well as the FENSE metric based on S-BERT. CB-score [START_REF] Martín-Morató | A summarization approach to evaluating audio captioning[END_REF] is omitted as it is precision-oriented, and is not intended for ranking in the FENSE setting. results with or without adding a penalty on fluency using the FENSE error detection model [START_REF] Zhou | Can audio captions be evaluated with image caption metrics?[END_REF]. SPICE+ slightly improves on SPICE, indicating that the pre-trained dependency annotator and additional graph parsing rules are beneficial. However, its absolute performance is still poor with less than 50% MM ranking accuracy on both datasets without accounting for fluency. Adding embedding similarity significantly increases performance, confirming the limitations of synonym tuple matching. SPICE+emb also significantly outperforms CIDEr and SPIDEr on MM and total accuracy, the current accepted metrics for audio captioning. Combining SPICE+emb with CIDEr to obtain a SPIDEr equivalent did not improve performance in our experiments.

Overall, SPICE+emb is worse than Sentence-BERT on total accuracy, with a smaller difference when penalizing fluency through error detection. The metric evaluates semantic elements rather than the entire caption, and thus cannot account for fluency. However, this approach provides more information on system performance by analyzing the contributions of elements to the final score, which is not possible with full sentence embedding comparison. The currently used SPIDEr score, as well as Sentence-BERT and SPICE+emb, are shown in Table 2 for systems of the DCASE2022 challenge on the Clotho-testing dataset 4 . SPICE+emb appears correlated to Sentence-BERT, and significantly alters system ranking compared to SPIDEr.

Comparison to SPICE behavior

Next, the behavior of SPICE+ on absolute scoring of captions is evaluated. Figure 1 presents the distributions of scores for HI pairs on Clotho, i.e. correct and incorrect captions. For SPICE, most incorrect captions are concentrated at 0, and give little information on the degree of similarity between candidates and references. Interestingly, the distribution is similar for SPICE+ with synonym matching, with only a slight shift of correct caption scores to higher values. In contrast, both embedding-based metrics SPICE+emb and Sentence-BERT provide a bimodal distribution, with some overlap. These metrics reliably assign clear correct and incorrect candidates with high and low scores respectively, which is consistent with the high performance (> 90%) on the HI case in Table 1.

Captioning system understanding

Systems are conditioned by training datasets to produce sentences with certain properties (e.g. length), which limits their ability to fully describe the acoustic scene in a single caption. When tasked with producing multiple captions, some systems may output variations of the same sentence while others may describe the scene more comprehensively at the cost of lower precision. Because SPICE+ compares semantic elements extracted from captions, it can also enable better understanding the behavior of captioning systems in such conditions. 4 Full results and technical reports at dcase.community/challenge2022 As a first experiment in this direction, SPICE+emb is computed for the DCASE2022 challenge baseline with varying number of generated captions from 1 to 20, obtained with a beam size of 25. For multiple candidates, semantic graph parsing is symmetrical, i.e. a single graph is constructed from all captions. The pseudo-precision and pseudo-recall tradeoff is shown in Figure 2. Both scores are stable and the recall does not increase to 1, indicating that the system outputs captions with similar content and is unable to match the diversity found in the reference set. SPICE+ is thus able to identify a limitation of the baseline system, and could provide a more comprehensive analysis of captioning systems in this manner.

CONCLUSION

We introduced SPICE+, a modification of SPICE including pre-trained language models for dependency annotation and graph comparison. The metric achieves competitive accuracy on a candidate ranking benchmark, but should be used in combination with a fluency measure due to its sensitivity to grammatical errors. Compared to full sentence embedding comparison, SPICE+ is more interpretable as it evaluates individual semantic elements. It is also adaptable through its hyper-parameters to the needs of specific applications that may arise in the community.

Still, some improvements should be investigated in future work. Parts of the captions are currently omitted from the graph, such as adverbs. The similarity between same-graph tuples is not accounted for in the score computation, which gives more importance to nodes with multiple attributes or relations. Lastly, the metric could be specialized to audio by adding a focus on salience and temporal relationships, which are not prevalent in current datasets of short acoustic scenes.

Fig. 1 .

 1 Fig. 1. Distributions of scores for correct and incorrect captions.

Fig. 2 .

 2 Fig. 2. SPICE+ precision and recall for varying number of candidates generated by the DCASE2022 challenge baseline.

  ± 6.2 96.0 ± 2.4 90.4 ± 3.7 61.2 ± 3.1 71.0 ± 2.2 51.4 ± 6.2 91.8 ± 3.4 70.3 ± 5.7 56.0 ± 3.1 63.2 ± 2.3 SPICE 50.2 ± 6.2 83.8 ± 4.6 77.8 ± 5.2 49.1 ± 3.2 59.7 ± 2.4 44.3 ± 6.2 84.4 ± 4.5 65.5 ± 5.9 48.9 ± 3.1 56.3 ± 2.3 SPIDEr 56.7 ± 6.1 96.0 ± 2.4 90.4 ± 3.7 63.4 ± 3.1 72.2 ± 2.1 53.3 ± 6.2 93.4 ± 3.1 70.3 ± 5.7 57.0 ± 3.1 64.2 ± 2.2 S-BERT (FENSE) 64.0 ± 6.0 99.2 ± 1.1 92.5 ± 3.3 73.6 ± 2.8 79.6 ± 1.9 60.0 ± 6.1 95.5 ± 2.6 75.9 ± 5.3 66.9 ± 2.9 71.8 ± 2.1 SPICE+ 59.1 ± 6.1 85.4 ± 4.4 83.7 ± 4.6 49.0 ± 3.2 62.0 ± 2.3 46.7 ± 6.2 88.1 ± 4.0 70.3 ± 5.7 48.7 ± 3.1 57.8 ± 2.3

	Standalone			AudioCaps					Clotho
	Method	HC	HI	HM	MM	Total	HC	HI	HM	MM	Total
	CIDEr 56.2 SPICE+emb 63.5 ± 6.0 96.4 ± 2.3 91.6 ± 3.4 70.0 ± 3.0 77.0 ± 2.0 61.0 ± 6.0 94.7 ± 2.8 76.3 ± 5.3 61.6 ± 3.0 68.9 ± 2.2
	With fluency			AudioCaps					Clotho
	Method	HC	HI	HM	MM	Total	HC	HI	HM	MM	Total
	CIDEr	56.2 ± 6.2 95.5 ± 2.6 90.4 ± 3.7 75.4 ± 2.8 78.6 ± 2.0 51.9 ± 6.2 91.8 ± 3.4 76.3 ± 5.3 68.8 ± 2.9 71.3 ± 2.1
	SPICE	50.7 ± 6.2 83.4 ± 4.6 78.7 ± 5.1 54.8 ± 3.2 62.8 ± 2.3 44.3 ± 6.2 84.4 ± 4.5 71.1 ± 5.6 57.7 ± 3.1 62.1 ± 2.3
	SPIDEr	57.1 ± 6.1 95.5 ± 2.6 90.0 ± 3.7 76.3 ± 2.7 79.1 ± 1.9 53.8 ± 6.2 93.4 ± 3.1 76.7 ± 5.2 68.9 ± 2.9 71.9 ± 2.1
	S-BERT (FENSE) 64.5 ± 5.9 98.4 ± 1.6 91.6 ± 3.4 84.6 ± 2.3 85.3 ± 1.7 60.5 ± 6.1 94.7 ± 2.8 80.2 ± 4.9 72.8 ± 2.8 75.7 ± 2.0
	SPICE+	59.1 ± 6.1 85.0 ± 4.4 84.5 ± 4.5 57.7 ± 3.2 66.8 ± 2.3 46.7 ± 6.2 88.1 ± 4.0 75.4 ± 5.3 56.8 ± 3.1 63.2 ± 2.3
	SPICE+emb	63.5 ± 6.0 95.5 ± 2.6 91.6 ± 3.4 83.4 ± 2.4 84.0 ± 1.8 61.0 ± 6.0 93.9 ± 3.0 81.0 ± 4.9 70.0 ± 2.8 74.1 ± 2.1

Table 1

 1

	shows the

Table 2 .

 2 DCASE2022 captioning systems evaluation on Clotho-testing. SPIDEr is the challenge ranking metric.


Hubert Ménou
email: hubert.menou@minesparis.psl.eu

Éric Bourgeois
email: eric.bourgeois@cnes.fr

Nicolas Petit
email: nicolas.petit@minesparis.psl.eu

MANAGING CONSTRAINTS FOR INFEASIBLE ROCKET LANDING PROBLEMS

For reusable launchers, in-flight landing trajectory planning is a challenging task. The naturally disturbed atmospheric environment may render the landing objectivessuch as the landing site location or the incidence safety boundsmutually incompatible. This article discusses a methodology to recover feasibility in over-constrained optimization problems used for Powered Descent Guidance. Nominal and emergency guidance problems are first discussed from a high-level viewpoint, and then developed to handle a Quadratic Programming based guidance method, which is illustrated on numerical examples for a planar rocket model.

INTRODUCTION

Powered landing for reusable launchers occurs in a naturally disturbed and uncertain environment. For instance, unexpected wind gusts can push a rocket far from the nominal landing site. The online Powered Descent Guidance (PDG) problem, aiming at providing a landing trajectory, must answer the following question at many instants during the flight: provided the current rocket position, orientation, and velocities, is landing still feasible at the landing site, and if not, what is a satisfactory solution?

Though the PDG problem is formulated in several different ways in the literature [START_REF] Song | Survey of autonomous guidance methods for powered planetary landing[END_REF]- [START_REF] Malyuta | Convex Optimization for Trajectory Generation[END_REF], a general common feature of the formulation is that they start by solving a feasibility problem for a constrained dynamical system. We are interested in the landing scenarios when this problem is not feasible anymoreas shown in Figure 1. Our work focuses on designing a method to keep providing a landing trajectory when the ideal target is not reachable anymore. Such a trajectory is called an emergency trajectory, and is computed by relaxing a constraint subset, using negotiable parameters.

This article is structured as follows. First, a highlevel description of the nominal and emergency guidance scenario is introduced in Section 2. More precisely, a systematic methodsee Algorithm 1 belowchooses the constraints that can or cannot be loosened. In Section 3, we show that if the initial PDG problem is formulated as a Quadratic Program (QP)using previous work of the authors [START_REF] Ménou | Sensitivity Analysis for Powered Descent Guidance: Overcoming degeneracy[END_REF] then the negotiability problem can be formulated as a Linear Program (LP). The latter makes real-time resolution possible by relying on off-the-shelf solvers. Numerical simulations are provided in Section 4, to demonstrate the above-mentioned concepts for various choices of negotiable parameters.

Figure 1:

The need for emergency guidance. Here, it is impossible to design a landing trajectory satisfying all the objectives for the left-most rocket initial conditions.

THE EMERGENCY PROBLEM

First, the variables describing the PDG problem are defined, and a generic description of nominal guidance method is introduced. Then, a high-level description of a procedure providing emergency guidance is presented.

Powered Descent Guidance variables

As will be detailed below, the PDG problem is described by its decision variable z, its input 𝜉, and several objective parameters 𝑝.

Consider a rocket described by its state 𝑥, conveying its positions, its speeds, and its mass. It is assumed that it is controlled by a variable 𝑢, conveying its engine flow 𝑞, and its incidence 𝛼.

The inputs of the PDG problem are the initial conditions of the rocketdenoted 𝑥 0and parametersdenoted 𝜂that influence the Equations of Motion (EoM), such as the wind magnitude (𝑤) or the Engine Specific Impulse (𝐼𝑠𝑝). The dynamic equation writes

𝑥̇= 𝑓(𝑥, 𝑢, 𝜂).

To alleviate the writing, these inputs are conveyed by a variable 𝜉 = (𝑥 0 , 𝜂).

We consider that 𝑢, which is an infinite dimensional variable in general, is approximated by a finite dimensional variable 𝜇. Basically, 𝜇 conveys the values of 𝑢 at some prescribed collocation points [START_REF] Ménou | Sensitivity Analysis for Powered Descent Guidance: Overcoming degeneracy[END_REF], [START_REF] Kraft | On Converting Optimal Control Problems into Nonlinear Programming Problems[END_REF].

Another important variable for the PDG problem is the time-of-flight 𝑡 𝑓 which is implicitly defined by the landing condition constraint "Altitude(𝑡 𝑓 ) = 0." Therefore, the finite-dimensional decision variable of our generic PDG problem is

𝑧 = ( 𝜇 𝑡 𝑓 ).
Finally, the PDG problem has several objectiverelated parameters, conveyed by a variable 𝑝. For instance, the position of the landing site or the incidence bound (respectively denoted 𝑧 𝑓 and 𝛼 max ) are problemspecific choices that can be present in this parameter 𝑝. Two examples will be presented in Section 4. The value of these parameters is 𝑝 nom for a nominal landing scenario.

The nominal PDG problem

Let us denote by ℱ(𝜉, 𝑝) the set of decision variables 𝑧 that satisfy the constraints above. Basically, ℱ conveys the dynamic equation 𝑓, the various bounds on the engine and the incidence, the landing site location constraint, and other bounds such as normal load limits. Also, denote by 𝐽(𝑧) a cost on 𝑧 that defines how desirable is the trajectory defined by 𝑧. Then, the problem of finding the best landing trajectory boils down to solving

PDG(𝜉) ∶= min 𝑧∈ℱ(𝜉,𝑝 nom ) 𝐽(𝑧) .
Most of the current research efforts in PDG focus on finding the best way to express this problem and using the fastest and most reliable computational methods to solve it [START_REF] Song | Survey of autonomous guidance methods for powered planetary landing[END_REF], [START_REF] Sagliano | Onboard Guidance for Reusable Rockets: Aerodynamic Descent and Powered Landing[END_REF], [START_REF] Malyuta | Convex Optimization for Trajectory Generation[END_REF].

PDG(𝜉

) is what we call the nominal guidance problem, when ℱ(𝜉, 𝑝 nom ) is a non-empty set. However, it is possible that for some values of 𝜉, the latter becomes infeasible, i.e., ℱ(𝜉, 𝑝 nom ) = ∅. Altering the value of 𝑝 to recover feasibility is what we call emergency PDG. Such a procedure is presented below.

A generic methodology for emergency guidance

Consider an infeasible case, and let us formulate a question worth answering, e.g., what should be done if the rocket starts its descent too far from the landing site?

To answer it, we propose a general approach, detailed below, that we apply on a specific landing procedure in the next section.

Emergency landing is required when the given input 𝜉 makes ℱ(𝜉, 𝑝 nom ) empty. This is due to the impossibility to achieve the target values 𝑝 nom . We propose to first compute the smallest perturbation of 𝑝 that recovers feasibility, and then to re-optimize it in view of the original PDG(𝜉) problem.

The perturbation in the variable 𝑝 is denoted Δ𝑝 such that 𝑝 = 𝑝 nom + Δ𝑝, where Δ𝑝 is bounded between Δ𝑝 -and Δ𝑝 + . This perturbation Δ𝑝 represent our negotiable parameters. Thus, we can introduce Negotiate(𝜉), also called the negotiation problem, which is min 𝑧, Δ𝑝 γ(Δ𝑝)

s. t. : 𝑧 ∈ ℱ(𝜂, 𝑝 nom + Δ𝑝), Δ𝑝 -≤ Δ𝑝 ≤ Δ𝑝 + .
for some arbitrary non-negative penalty function 𝛾 satisfying

𝛾(Δ𝑝) = 0 ⇔ Δ𝑝 = 0.
Typically, 𝛾 is a norm such as the 1-norm or the infinity norm. We are interested in the optimal value of Negotiate(𝜉), denoted 𝒫 * , which allows us to define the problem Refine(𝜉, 𝒫 * ) such that min 𝑧, Δ𝑝 𝐽(𝑧)

s. t.: 𝑧 ∈ ℱ(𝜉, 𝑝 nom + Δ𝑝), Δ𝑝 -≤ Δ𝑝 ≤ Δ𝑝 + , 𝛾(Δ𝑝) = 𝒫 * .
The fact that Δ𝑝 is free in the constraints of Refine(𝜉, 𝒫 * ) means that we seek the optimal value of z such that the objective parameters are as close as possible in the sense imposed by 𝛾from their nominal value 𝑝 nom . The high-level procedure for Emergency Powered Descent Guidance (E-PDG) is described by Algorithm 1.

Input: 𝜉

Step 1: Solve Negotiate(𝜉), which gives an optimal penalty 𝒫 * on the negotiable parameters.

Step 2: Solve Refine(𝜉, 𝒫 * ), which gives an optimal landing strategy 𝑧 * . Output. 𝑧 * Algorithm 1: Emergency Powered Descent Guidance (E-PDG)

Note that this algorithm does not expect Negotiate(𝜉) to have a unique solution (otherwise Step 2 is useless), since only the value of 𝒫 * matters.

Moreover, E-PDG is compatible with the nominal guidance. In other words, for all the values of 𝜉 for which PDG(𝜉) is feasible, E-PDG and PDG(𝜉) will return the same value. Indeed, if ℱ(𝜉, 𝑝 nom ) is feasible for some input 𝜉, then Negotiate(𝜉) will return 𝒫 * = 0, which implies Δ𝑝 = 0. In this case, Refine(𝜉, 0) is directly equivalent to PDG(𝜉). This is one salient advantage of E-PDG: it can be plugged in place of any arbitrary nominal guidance method.

QP/LP EMERGENCY GUIDANCE

In practice, we applied the E-PDG algorithm to a Quadratic Programming based guidance method developed by the authors [START_REF] Ménou | Sensitivity Analysis for Powered Descent Guidance: Overcoming degeneracy[END_REF].

PDG using Quadratic Programming

Extending the notations from [START_REF] Ménou | Sensitivity Analysis for Powered Descent Guidance: Overcoming degeneracy[END_REF], the main problem PDG(𝜉) is described using a quadratic cost and linear constraints, which we will denote QP(𝜉) min 𝑧 where 𝐺, ℎ, 𝐻, 𝐴, 𝑏, 𝐵 are matrices that depend on a mission-specific reference trajectory. Here, ℎ and 𝑏 are affine functions of the input 𝜉. Here, 𝑊simply denotes a weighting matrix, assumed positive-definite.

Application of E-PDG

The penalty 𝛾 is here the 1-norm. The negotiation problem is denoted Negotiate LP (ξ) and consist in Note, as suggested by its name, that Negotiate LP (ξ) can be solved using Linear Programming, by decomposing Δ𝑝 into its positive and negative parts. See Example 1.13 in Betts' book [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF] for further details. Likewise, Refine QP (𝜉, 𝒫 * ) is solved using Quadratic Programming [START_REF] Andersen | CVXOPT: Convex Optimization[END_REF].

It is noteworthy that the problem Negotiate LP (ξ), aiming at recovering feasibility, builds upon mathematical programming techniques known as righthand side alteration methods [START_REF] Chinneck | Feasibility and Infeasibility in Optimization: Algorithms and Computational Methods[END_REF]. Here, as conveyed by the matrices 𝐻 𝑝 and 𝐵 𝑝 and due to the limited number of levers, only a sub-space of the constraint right-hand side can be altered.

A geometric interpretation

The connections between problems QP(𝜉), Negotiate LP (ξ) and Refine QP (𝜉, 𝒫 * ) can be easily understood from a geometric point of view. On a toyproblem, let us see with a simplifying perspective how these work by projecting the constraint on a planar representation.

The constraints of QP(𝜉) define a polytope: the inequalities define half-spaces, and the equalities define hyperplanes. Thus, the feasible set for the variable 𝑧 in 𝑄𝑃(0) is the set represented in blue in Figure 2. When the inputs are non-zero, for example when Δ𝑥 0 ≠ 0, the right-hand side vectors ℎ and 𝑏 of the constraints of QP(𝜉) are modified. Geometrically, it corresponds to a pure translation of the lines representing the constraints, as shown in Figure 3. However, if Δ𝑥 0 is too large (for example when the rocket starts its descent too far from the landing site, as pictured in Figure 1), then satisfying all the constraints simultaneously becomes impossible, as Figure 4 suggests. Hence the need for the E-PDG algorithm. Let us assume that the final horizontal position Δ𝑧 𝑓 can be negotiated, i.e., we take 𝑝 = Δ𝑧 𝑓 . Modifying its value will have a direct impact on the right-hand side of the equality constraints from QP(𝜉). By solving Negotiate LP (ξ), we find the smallest value of Δ𝑧 𝑓 that recovers feasibility, as illustrated in Figure 5.

NUMERICAL SIMULATIONS

Let us apply the above-mentioned guidance algorithm, to two scenarios conveying different negotiation parameters.

Both scenarios presented below consider inputs defined with respect to a reference trajectory, shown in plain black in Figures 6 and7. Several input values 𝜉 that differ in their initial horizontal position Δ𝑧 0 are presented. Note that to ease the comparison between them, Figures 6 and7 are plotted for the same input list.

FINAL HORIZONTAL POSITION

Let us assume that the landing site is very large (for instance if the rocket landing site is in a flat desert). Thus, if necessary, the final horizontal position Δ𝑧 𝑓 can be negotiated. As shown in Figure 6, if Δ𝑥 0 becomes large enough, negotiation is needed, and it yields the farthest trajectories (in dashed blue).

INCIDENCE BOUND

Contrary to the previous example, let us assume that the final horizontal position is not negotiable, which happens if one whishes to land on a platform at see for instance. Instead, we assume that the safety bound 𝛼 𝑚𝑎𝑥 imposed on the rocket incidence 𝛼, can be negotiated by a parameter 𝑝 = Δ𝛼 𝑚𝑎𝑥 such that |𝛼| ≤ 𝛼 𝑚𝑎𝑥 + Δ𝛼 𝑚𝑎𝑥 .

Considering the same input values as in the previous example, we see that the farthest trajectories from the reference need to negotiate the value of Δ𝛼 𝑚𝑎𝑥 to recover feasibility, leading to a sharped turn for the dashed pink trajectories shown in Figure 7. the incidence bound is allowed to be negotiated.

FUTURE WORK

This paper presented a generic methodology to provide emergency guidance for infeasible Powered Descent Guidance problems and has been illustrated on a 2D rocket model. Future work will focus on extending the E-PDG algorithm to a richer framework that enforces a hierarchy of importance between the components of 𝑝 (see for example [START_REF] Ménou | Nominal And Emergency Rocket Landing Guidance Using Quadratic Programming[END_REF], which has been accepted after presenting our work at the 25 th ESA Symposium, but before the final editing of this paper). Moreover, extension to more complex rocket models, with 3D descriptions and nontrivial engine dynamics are under study.

  + 𝐻 𝑝 (𝑝 nom + Δ𝑝) 𝐴𝑧 = 𝑏 + 𝐵 𝑝 (𝑝 nom + Δ𝑝) Δ𝑝 -≤ Δ𝑝 ≤ Δ𝑝 + and returns its optimal value 𝒫 * . The refine problem is denoted Refine QP (𝜉, 𝒫 * ) + 𝐻 𝑝 (𝑝 nom + Δ𝑝) 𝐴𝑧 = 𝑏 + 𝐵 𝑝 (𝑝 nom + Δ𝑝) Δ𝑝 -≤ Δ𝑝 ≤ Δ𝑝 + ‖Δ𝑝‖ 1 = 𝒫 *

Figure 2 :

 2 Figure 2: The constraints of QP(0).

Figure 3 :

 3 Figure 3: The constraints of QP(𝜉),when Δ𝑥 0 ≠ 0.

Figure 4 :

 4 Figure 4: The constraints of QP(𝜉),when Δ𝑥 0 makes the constraints infeasible.

Figure 5 :

 5 Figure 5: The constraints of QP(𝜉), when Δ𝑝 = (Δ𝑧 𝑓 ) helps recover feasibility.

Figure 6 :

 6 Figure 6: Nominal and emergency trajectories, when landing site location is allowed to be negotiated.

Figure 7 :

 7 Figure 7: Nominal and emergency trajectories, whenthe incidence bound is allowed to be negotiated.


Hubert Ménou

Eric Bourgeois

Nicolas Petit

Sensitivity Analysis for Powered Descent Guidance: Overcoming degeneracy

This paper exposes a method to handle inaccurate modeling and/or initial state errors during the Powered Descent Guidance (PDG), a critical phase of atmospheric rocket landing. For this, we develop a replanification method having reliable online computational capabilities. From a reference descent scenario, an optimal correction problem is formulated. After revisiting results on Non Linear Programming sensitivity for degenerate optimization problems, we conclude that Quadratic Programming (QP) provides a local solution to the replanification problem. Using three illustrative PDG scenarios, we stress degeneracy and show how QP is used to evaluate the upper Dini derivatives at stake. Further, we discuss to what extent QP also provides a quantitatively reasonable solution outside a small neighboorhood of the reference scenarios.

I. Introduction

Due to its critical nature, Powered Descent Guidance (PDG), the final phase of atmospheric rocket landing, is often addressed as an Optimal Control Problem (OCP) having many parameters such as initial states, aerodynamic coefficients and engine efficiency [START_REF] Steinfeldt | Guidance, Navigation, and Control System Performance Trades for Mars Pinpoint Landing[END_REF], [START_REF] Blackmore | Autonomous Precision Landing of Space Rockets[END_REF]. No known analytic solution is available for the general optimal atmospheric PDG problem. Generating a numerical solution in real-time is challenging, especially with nonlinear dynamics and multiple constraints [START_REF] Szmuk | Successive Convexification for 6-DoF Mars Rocket Powered Landing with Free-Final-Time[END_REF]. To this purpose, it is of interest to explore a parametric problem description and perform first-order extrapolations in the vicinity of reference solutions.

Sensitivity of parametric OCPs have been studied for a long time [START_REF] Bensoussan | Perturbation methods in optimal control[END_REF]. Assuming that one knows a solution for some prescribed parameter value, the idea is to deduce the neighboring solutions for small parameter variations [START_REF] Büskens | Sensitivity Analysis and Real-Time Control of Parametric Optimal Control Problems Using Nonlinear Programming Methods[END_REF]. Due to the practical drawbacks of indirect methods [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF]Sec. 4.3], we will focus on solving OCPs via direct methods: OCPs are described, or at least approximated, using finite-dimensional variables, and then solved using Nonlinear Programming (NLP) [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF]- [START_REF] Vlassenbroeck | A chebyshev polynomial method for optimal control with state constraints[END_REF]. Therefore, sensitivity analysis boils down to studying parametric NLP.

Classical results state under mild conditions that if, for the reference parameter, the multipliers are (strictly) H. Ménou (hubert.menou@mines-paristech.fr) and N. Petit (nicolas.petit@mines-paristech.fr) are with MINES Paris-Tech, Centre Automatique et Systèmes, PSL University, 60 bd. St Michel, 75272 Paris Cedex, France, and E. Bourgeois (eric.bourgeois@cnes.fr) is with CNES, Direction des lanceurs, 52, rue Jacques Hillairet, 75612 Paris Cedex, FRANCE.

This work was supported by CNES and Mines ParisTech.

positive, then the optimal solution is continuously differentiable [START_REF] Fiacco | Introduction to Sensitivity and Stability Analysis in Nonlinear Programming[END_REF]. Such problems are called non-degenerate.

Computing the first order expansion of their solution only requires to invert a linear system, resulting from the application of the Implicit Function Theorem (IFT) on the Karush-Kuhn-Tucker (KKT) conditions [9, Thm. 3.2.2], [START_REF] Büskens | Sensitivity Analysis and Real-Time Control of Parametric Optimal Control Problems Using Boundary Value Methods[END_REF]. This is a powerful result, with many applications [START_REF] Braun | Post-optimality analysis in aerospace vehicle design[END_REF], [START_REF] Reiter | Real-Time Computation of Inexact Minimum-Energy Trajectories Using Parametric Sensitivities[END_REF]. However, this positive multiplier assumption is not always satisfied, and degenerate problems with non-differentiable solutions are often encountered1 . This paper focuses on two elements. First, it revisits a method to handle sensitivity analysis at degenerate points, for which weaker kind of differentiability is available -namely, existence of the upper Dini derivativeand details how to compute it. Second, it reports the key observation that PDG scenarios can be degenerate.

The paper is organized as follows. First, in Section II, a general class of OCP is formulated and then discretized. Secondly, in Section III, NLP sensitivity results are recalled and several computational aspects are discussed. Then, this framework is applied to rocket landing in Section IV. A rocket model is introduced and the necessity to be able to handle degenerate scenarios is proved. Also, a way to assess the proper reference parameter vicinity -ie. where the corrections remain valid -is presented. Finally, comments and extensions are provided in Section V.

Notations

Vector inequalities must be understood element-wise. For any map h : R p → R q , denote its jacobian h x := ∂h ∂x ∈ R q×p . When it exists, the upper Dini derivative of h at x in the direction d is denoted by

D + d h(x) := lim t↓0 h(x + td) -h(x) t ∈ R q .
II. Optimal Control Problem First, a class of optimal control problems is introduced, and then discretized via a direct approach.

A. General correction problem

Consider a controlled system defined by its dynamic function f . For a control u : [0, 1] → R m , a parameter p ∈ R Np , and for t ∈ [0, 1], the flow Φ f of f is defined such that Φ f (t, x 0 , p; u) = x(t) where the state x of dimension n satisfies the Initial Value Problem (IVP)

x(0) = x 0 and ẋ = f (x, u, p). ( 1
)
If f is twice differentiable, so is Φ f [START_REF] Bonnans | Course on Optimal Control, Part I: the Pontryagin approach[END_REF]Prop. 2.46]. Define a parametric control u(z 1 , t) where z 1 ∈ R N1 is the parameter, such that u(0, .) ≡ 0 over [0, 1]. Typically, u can denote piecewise polynomial functions in t with linear dependency on z 1 [START_REF] Kraft | On Converting Optimal Control Problems into Nonlinear Programming Problems[END_REF]Sec. 3]. Also, denote the vector z 2 ∈ R N2 s.t. z = (z 1 , z 2 ) ∈ R N . Provided a reference control ū, we are interested in the following optimal control problem min

z∈R N J(z, p) (2a) s.t: ẋ(t) = f (x(t), ū(t) + u(z 1 , t), p), ∀t ∈ [0, 1] (2b) x(0) = h(p, z 2 ) (2c) g(x(t), ū(t) + u(z 1 , t)) ≤ 0, ∀t ∈ [0, 1] (2d) Ψ(x(1)) = 0 (2e)
where J, h, g, Ψ are twice-differentiable. Remark 1: To handle problems with a final time t f other than 1, it is sufficient to consider the scaled time t = t/t f and multiply the dynamic equation by t f . For free-final time problems, also add t f as an extra state of trivial dynamic in f , as shown later in Section IV-B.

We assume that ū, a reference control, is such that: Assumption 1: If p = 0, then z = 0 satisfies the mixed state-control constraints (2d) and the terminal constraints (2e).

It is important to note that (2) is a correction problem. Indeed, ū does not need to be the solution of any optimization problem. It can be whatever control that steers the system from x 0 to the condition Ψ(x(1)) = 0, while enforcing g(x, ū) ≤ 0. Naturally, the choice of ū has a direct impact on the sensitivity analysis that follows.

Also, note that J is general enough to convey a Mayer cost [START_REF] Kraft | On Converting Optimal Control Problems into Nonlinear Programming Problems[END_REF] by considering a map ϕ : R n → R and

J(z, p) = ϕ (Φ f (1, h(p, z 2 ), p; u(z 1 , .))) .
(

From a practical perspective, computing the double derivative of (3) requires to evaluate the second derivative of f (see eg. [START_REF] Bonnans | Course on Optimal Control, Part I: the Pontryagin approach[END_REF]Sec. 3.2]) which may require a significant analytical and numerical efforts in nontrivial contexts. To simplify this matter, J may be taken directly quadratic in z and mixed linear in p and z s.t.

J(z, p) = 1 2 z ⊤ W z + p ⊤ Qz (4)
for some weighting matrices W (symmetric) and Q.

B. Discretization and NLP

To obtain a finite dimensional optimization problem with a finite constraint number, (2d) can be discretized following the direct collocation approach [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF].

Denote by 0 = t 0 < t 1 < . . . < t Nt = 1 an arbitrary subdivision of [0, 1] in N t intervals. We choose to enforce (2d) only at those prescribed time instances.

Consequently, Problem 2 writes as a nonlinear program

NLP(p) := min z∈R N J(z, p) (5a) s.t: G(z, p) ≤ 0 (5b) H(z, p) = 0 (5c) where x k (z, p) := Φ f (t k , h(p, z 2 ), p; u(z 1 , .
)) fully conveys the dynamics equation and

G(z, p) :=    g(x 0 (z, p), ū(t 0 ) + u(z 1 , t 0 )) . . . g(x Nt (z, p), ū(t Nt ) + u(z 1 , t Nt ))    , (6a) H(z, p) := Ψ(x Nt (z, p)). ( 6b
)
We will now focus on NLP(p).

III. Sensitivity Analysis Results

The shape of NLP(p) is classic [START_REF] Fiacco | Nonlinear Programming[END_REF]. Provided a solution z 0 corresponding to a parameter p 0 , the goal is to approximate the optimal value z * (p) near p = p 0 .

In this section, results defining the first order expansion of z * at p 0 are recalled, and a method to effectively compute it is discussed. Without loss of generality, assume p 0 = 0 and z 0 = 0.

A. Directional differentiability

For multipliers η and λ, introduce the Lagrangian

L(z, η, λ, p) := J(z, p) + η ⊤ G(z, p) + λ ⊤ H(z, p). (7)
Recall that a tuple (z, η, λ, p) is said to satisfy the KKT conditions if η ≥ 0 and

L z (z, η, λ, p) = 0, H(z, p) = 0, η ⊤ G(z, p) = 0. (8)
At p 0 , the multipliers satisfying the KKT conditions are denoted η 0 and λ 0 . The goal is to highlight general conditions under which the following (directional) expansion is defined

z * (p) = z 0 + D + p z * (0) + o(∥p∥). ( 9
)
Strict Complementarity Slackness (SCS), ie. η 0 > 0, is often assumed to prove local uniqueness and differentiability of z * , by applying the IFT on the KKT conditions [START_REF] Fiacco | Introduction to Sensitivity and Stability Analysis in Nonlinear Programming[END_REF]. However, this assumption does not necessarily hold, even on simple examples, as shown in the introduction. To circumvent this limitation, one can use Strong Second Order Sufficiency Conditions (SOSC)2 , as expressed in [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF].

Definition 1 (Strong SOSC):

There exist a > 0 s.t.

ν ⊤ L zz (z 0 , η 0 , λ 0 , p 0 )ν ≥ a∥ν∥ 2 (10)
for all ν ∈ R Nz s.t. H z (z 0 , p 0 )ν = 0 and G j z (z 0 , p 0 )ν = 0 for every component j s.t. G j (z 0 , p 0 ) = 0 and (η 0 ) j > 0.

As shown by Jittorntrum [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF], relaxing the SCS using Strong SOSC instead eventually leads to directional differentiability properties, in the sense of the upper Dini derivative. For any direction p, consider problem NLP(tp), where t ≥ 0 is a scalar.

Theorem 1 (From [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF]): Assume that at z 0 (i) the linear independence condition holds (ie. the columns of H z (z 0 , p 0 ) and G j z (z 0 , p 0 ) for the active constraints j are linearly independent), (ii) Strong SOSC holds, (iii) J, G, and H are twice continuously differentiable in a neighborhood of z 0 . Then, there exists a unique continuous function t → (z * (tp), η * (tp), λ * (tp)) as the (local) solution of NLP(tp), for t ≥ 0, such that z * (0) = z 0 , η * (0) = η 0 and λ * (0) = λ 0 . Furthermore, its right-hand derivative exists.

Thus, for any p and since z 0 = 0, the following holds

z * (tp) = D + p z * (0)t + o(t). ( 11
)
Remark 2: Directional differentiability results also exist via the indirect approach for parametric OCP [START_REF] Deshpande | Directional Input Adaptation in Parametric Optimal Control Problems[END_REF].

B. Computation method

When SCS is assumed, computing the expansion of z * only requires to invert a linear system [START_REF] Fiacco | Introduction to Sensitivity and Stability Analysis in Nonlinear Programming[END_REF]. However, this system becomes singular when SCS is not satisfied and the correction cannot be computed, hence the need for the above-mentioned tools, and the computational method detailed below. Jittorntrum's proof for Theorem 1 relies on a side Quadratic Program (QP) to serve intermediate theoretical purposes, which actually has powerful practical use.

Proposition 1 (Adapted from [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF]Eq. 24]):

Solving 3 QP(p) := min ∆z∈R N 1 2 ∆z ⊤ L zz [0]∆z + p ⊤ L pz [0]∆z (12a) s.t: G[0] + G z [0]∆z + G p [0]p ≤ 0 (12b) H z [0]∆z + H p [0]p = 0 (12c)
provides the value of ∆z := D + p z * (0). Similar methods can be found in more recent work as well, see e.g. Bonnans & Shapiro [17,Sec. 5.2].

The primary goal of QP(p) is to compute an approximation of z * (p) in the direction p, which works locally according to Theorem 1. However, inequality constraints that are inactive at z 0 can become active when p increases (in norm). Provided that the approximation (12b) represent well the constraints, it will help QP(p) to provide reasonable approximations of NLP(p). In other words, QP(p) manages non-local constraints using first-order constraint approximations [START_REF] Büskens | Sensitivity Analysis and Real-Time Control of Parametric Optimal Control Problems Using Boundary Value Methods[END_REF].

The main advantage of solving QP(p) instead of NLP(p) is that one can rely on fast, mature and dependable 4 QP solvers [START_REF] Andersen | CVXOPT: Convex Optimization[END_REF], [START_REF] Grant | CVX: Matlab Software for Disciplined Convex Programming[END_REF], and thus can use it for realtime applications. However, note that QP(p) cannot be used as a brick for Successive QP, as the intermediate matrices would take too much time to compute.

As a summary of the preceding discussion, provided a reference control ū for Problem 2, computing in advance J, G, H, their first and second derivatives at (z 0 , p 0 ), and checking the assumptions of Theorem 1, one is able to reliably compute the expansion [START_REF] Braun | Post-optimality analysis in aerospace vehicle design[END_REF] with low computational effort using QP(p).

A useful simplifying aspect should be noted. If J is of the form (4), Assumption 1 implies that η 0 = 0 and λ 0 = 0 necessarily hold. Indeed, z = 0 globally minimizes z → J(z, 0) and z = 0 satisfies the constraints. Therefore, considering that the second order derivatives of G and H only appear in the cost of QP(p), computing QP(p) requires only the first order derivatives of G and H. Thus, (12a) boils down to (4). This remark actually remains valid for general J provided that J z (z 0 , p 0 ) = 0.

IV. Application to Rocket Landing

The preceding results are now applied to the PDG problem. First, a rocket model and the choices corresponding to Problem (2) are described. Then, three examples showing the importance of the above-mentioned results are discussed. To alleviate the writing, variable dependencies are only written when truly necessary to avoid ambiguity.

A. Rocket model

The approach is applied to PDG with a planar rocket model [START_REF] Vinh | Hypersonic and planetary entry flight mechanics[END_REF], [START_REF] Souza | Survey of planetary entry guidance algorithms[END_REF], pictured in Figure 1. For h ≥ 0 the altitude, v h < 0 the vertical speed, z the horizontal position, v z the horizontal speed and m ≥ m dry the total mass of the rocket, the dynamic equations are ḣ = v h , ż = v z and vh

= -g + 1 m (F L sin θ + (T -F D ) cos θ) (13a) vz = 1 m (F L cos θ -(T -F D ) sin θ) (13b) ṁ = -q (13c)
where q is the engine flow, θ the attitude, and

(Thrust) T := gI SP q -P atm (h)S E , (14a)

(Speed) V r := v 2 h + (v z -w(h)) 2 , ( 14b
) (Drag) F D := 1 2 ρ(h)V 2 r S R C D (M a , α), ( 14c
) (Lift) F L := 1 2 ρ(h)V 2 r S R C L (M a , α), ( 14d
) (Mach) M a := V r /S SP (h), ( 14e
) (Incidence) α := arctan v z -w(h) |v h | -θ. (14f)
Here, g denotes the gravity, I SP the engine specific impulse, P atm and ρ the atmospheric pressure and density, C D and C L the drag and lift coefficients, S SP the sound speed, w the horizontal wind and S R and S E characteristic rocket surfaces. The pressure bias in the thrust expression conveys the air flow envelop effect on the rocket when landing. When the rocket goes down vertically with no wind, α = θ = 0. The wind profile is assumed to be piece-wise affine, null at h = 0 and defined by its values w 1 at 5 km and w 2 at 10 km.

B. OCP and NLP

The elements needed to build Problem 2 and NLP(p) are presented here. The state x is (h, v h , z, v z , m) ⊤ and the control u is (q, α) ⊤ . Depending on the examples presented below, p may convey any set of parameters describing the rocket model or its initial condition. It yields a dynamic equation ẋ = f (x, u, p).

The terminal condition (2e) here corresponds to

h(t f ) = z(t f ) = v z (t f ) = 0 and v h (t f ) = -ε f
where ε f > 0 denotes a non-zero small landing speed. For any time t, the constraints (2d) correspond to

q min ≤ q(t) ≤ q max , α min ≤ α(t) ≤ α max . ( 15
)
Since the terminal condition implicitly determines the free final-time t f , we add it to the dynamics as an extra state x := (x ⊤ , t f ) ⊤ . Thus, for all s in [0, 1], it yields

ẋ(s) = f (x(s), u(s), p) := t f (s)f (x(s), u(s), p) 0 . ( 16
)
For the next examples, consider N t = 4 and a parametric correction conveyed by a Cubic Spline [START_REF] Kraft | On Converting Optimal Control Problems into Nonlinear Programming Problems[END_REF], which requires the knowledge of u at each time t 0 , . . . , t Nt (taken uniformly distributed) and its slopes at t 0 and t Nt . Thus, z equals z := (u(t 0 ) ⊤ , . . . , u(t Nt ) ⊤ , u(t 0 ) ⊤ , u(t Nt ) ⊤ , τ ) ⊤ [START_REF] Bonnans | Perturbation Analysis of Optimization Problems[END_REF] where u(t k ) := (q(t f ), α(t k )) ⊤ , τ denotes the free-final time change and N = 2(N t +3)+1 = 15. The cost (2a) is then taken quadratic in z, where W from ( 4) is diagonal, with only positive scaling weights w q (resp. w α and w t f ) on components associated to the flow (resp. the incidence and the final time), ie. J(z, p) = 1 2 z ⊤ W z and W = Diag w q , w α , w q , w α , . . . , w t f . (

) 18
After having specified the framework for Problem 2 and NLP(p), let us focus on the reference control ū that will be used below. Consider a reference scenario of trajectory x, shown in black in Figure 3-(e), starting at x 0 , flying for t f , with a constant engine flow q ≡ q cst . The reference incidence control ᾱ is such that it switches once between negative and positive values, as shown in Figure 3-(a) in black. ᾱ is build such that it reaches its maximum value α max at some time t 2 , and remains far from its minimum value α min = -α max on the other side.

The assumptions of Theorem 1 are satisfied. First, the linear independence assumption is checked numerically, on the terminal condition and the only active inequality constraint, namely α(t 2 ) = α max . Second, as noted at the end of Section III, such a quadratic cost J leads to null multipliers for p = 0 (even though one control constraint is active). Thus, L zz [0] = W is positive definite, showing Fig. 1. Forces and angles at stake. Fig. 2. Sensitivity of the terminal horizontal position with respect to variations of the I SP . In red is z(t f ) when no correction is applied to I SP variations. In dashed orange is the first order estimation of the red curve, extracted from the problem linearization. In blue is z(t f ) when applying the corrections given by QP(∆I SP ). Terminal condition (2e) is thus locally well satisfied after correction.

that Strong SOSC holds. It also shows that SCS is not satisfied since η 0 = 0.

In practice, QP(p) has been implemented using cvxopt [START_REF] Andersen | CVXOPT: Convex Optimization[END_REF], and the flow derivatives required in (6) using material from [START_REF] Bonnans | Course on Optimal Control, Part I: the Pontryagin approach[END_REF]. We can now dive into the examples.

C. Direct application to I SP error correction

The first example aims at showing that in some neighborhood of the reference parameters, the expected expansion [START_REF] Braun | Post-optimality analysis in aerospace vehicle design[END_REF] provides accurate corrections enforcing the terminal constraints.

In this example, let us consider a single parameter : the error on the engine specific impulse p = ∆I SP := I SP -I 0 SP , where I 0 SP is the reference value. Due to the terminal constraints, the four first states, namely h, v h , z, v z are expected to be null at the final time (except v h (t f ) = -ε f ). If the landing occurs with p ̸ = 0, but without correcting it, non-zero values are Fig. 3. Illustration of the non-differentiability of the optimal correction on a standard scenario. The landing trajectory is shown in (e). The landing trajectories have been re-played on the nonlinear system, after using the correction from QP(p). Blue features correspond to negative parameter values ∆z 0 < 0 and red features to ∆z 0 > 0. The focus is mainly on the incidence correction coefficient ∆α(t 2 ). Its greatest magnitude value, defined in (b), is denoted α. The non-smooth behavior of ∆α(t 2 ) w.r.t. parameter changes is clearly seen in (d).

expected to appear. However, if the change in parameter is corrected using Proposition 1, the terminal constraints are assumed to be locally satisfied.

These two behaviors are well seen in Figure 2, representing the terminal horizontal position z(t f ). The other terminal constraints on h, v h and v z have voluntarily been omitted, as their correction curve is much flatter than for z. In other words, it means that even though this terminal constraint component has the "worst" correction curve, it still demonstrates that the first order corrections brought by QP(p) work well in a non-trivial neighborhood of p = 0 in practice.

D. Non-smooth sensitivity to initial horizontal position

The second example aims at showing that the optimal solutions of NLP(p) are indeed only Dini-differentiable and not smooth, even in a standard landing scenario.

Let us consider a single parameter, the error in initial horizontal position p = ∆z 0 := z 0 -z0 .

The directional-derivatives of α at the middle point t 2 in the directions ∆z 0 = 1 and ∆z 0 = -1 differ, as shown in Figure 3-(d). This behavior has a real-world interpretation. When ∆z 0 > 0, using more incidence on α(t 2 ) is not an option as the constraint is already active and becomes strictly active (the associated multiplier becomes positive when ∆z 0 > 0). However, when ∆z 0 < 0, lowering α(t 2 ) is possible, allowing the corrections presented.

E. Non-local I SP and initial horizontal position error correction

This last example focuses on assessing whether the corrections from QP(p) can be used non-locally.

Let us first define a metric that assesses a correction correctness. When a parameter p ̸ = 0 is estimated, a correction is computed using QP(p), and then used on the non-linear system in an open-loop fashion. The firstorder changes in the terminal constraint is supposed to be null locally near p ̸ = 0. However, for significant values of p, nonlinear behavior will eventually make the corrections fail. Thus, we seek a performance function Γ s.t. Γ(p) ≤ 1 if and only if the terminal constraints are acceptable. Given some scaling factors w h , w v h , w z , w vz , consider the map

µ(x) := max |h| w h , |v h + ε f | w v h , |z| w z , |v z | w vz . ( 19
)
When µ is applied on the final state x(t f ), the scaling factors define what terminal constraint errors are acceptable. Thus, using the notations from ( 6), a suitable performance function is given by evaluating the terminal point

Γ(p) := µ x Nt D + p z * (0), p . (20)
Then, building upon the two previous example, consider the two-dimension parameter p = (∆z 0 , ∆I SP ) ⊤ . Evaluating Γ on a square grid of 51 × 51 elements, it yields the level-set map shown in Figure 4. It shows a non-trivial area where first-order corrections from QP(p) implies acceptable terminal condition satisfaction on the the true nonlinear system.

Note that this kind of green valleys that extend on the figure corners, or "cross pattern", conveys the fact that some errors compensate each other in a constructive way. Also, the non-smooth behavior of the contours at the bottom of Figure 4 correspond to the activation of one engine flow control constraint, hence this sharp rupture.

V. Extensions and discussions

First, note that the choice of parameters in Section IV is arbitrary, and the number of parameters was chosen small for visualization purposes. It does not change the approach to add a great number of parameters (e.g. detailed atmospheric profiles, non-linear gravity model).

Second, linear parametric corrections u are advantageous. Indeed, if u(z, t) is taken linear in z (e.g. Cubic Splines or Hermite polynomials [START_REF] Kraft | On Converting Optimal Control Problems into Nonlinear Programming Problems[END_REF]) and if ū is described the same way, then the linearization (12b) of bounded control constraints, such as [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF], is actually exact.

Moreover, more intricate mixed state-control constraints are under investigation (e.g. acceleration bounds), and will be the matter of future publications.

Finally, consider an offline/online procedure consisting in first computing a library of reference trajectories and its sensitivities before the flight, and then using QP(p) on the closest reference trajectory during the flight. Such an approach has already been described in previous works [START_REF] Büskens | Sensitivity Analysis and Real-Time Control of Parametric Optimal Control Problems Using Boundary Value Methods[END_REF] and used in robotics [START_REF] Reiter | Real-Time Computation of Inexact Minimum-Energy Trajectories Using Parametric Sensitivities[END_REF] for example. The need to be able to guarantee that the pre-computed trajectory library is rich enough to cover all needs is crucial for critical applications [START_REF] Clark | A Study on Run Time Assurance for Complex Cyber Physical Systems[END_REF]. Thus, checking that the union of acceptable level sets for every reference trajectory -as in Figure 4 -represents a finite-covering of the possible inputs would demonstrate the completeness of the approach.

VI. Conclusion

This paper has revisited useful results from the literature, regarding sensitivity analysis for optimal control problem. Its application to Powered Descent Guidance proved the need and how to circumvent the strict complementarity slackness assumption. Practical methods have been illustrated to assess the quality of the approximated corrections on the true nonlinear system.

Fig. 4 .

 4 Fig. 4. Level sets for the performance function defined in (20). The contours in green denotes values lower than 1, where the parameter changes are corrected by QP(p) within an acceptable tolerance. In orange shades are the values greater than 1.

A simple example being minimizing z

such that z ≥ p, whose solution is z = 0 for p < 0 and z = p for p ≥ 0.

For a useful interpretation of this condition, see [10, Sec. 2.3].

The compact notation L[0] = L(z 0 , η 0 , λ 0 , p 0 ) is used throughout the paper. Similar notation is used for H and G.

[START_REF] Bensoussan | Perturbation methods in optimal control[END_REF] Not suffering from initial guess sensitivity to converge.