domain: spi phys

text:
C Musso

F Dambreville
email: frederic.dambreville@onera.fr

C Chahbazian
email: clement.chahbazian@onera.fr

Filtering and sensor optimization applied to angle-only navigation

Passive target estimation is a widely investigated problem of practical interest for which particle filters represent a popular class of methods. We propose an adaptation of the Laplace Particle Filter applied to angle-only navigation using landmarks. In this specific context, a high number of aiding landmarks or features could be hard to handle in terms of computational cost. Hence, this paper introduces a Cross-entropy algorithm that selects landmarks having a high contribution to the state estimation. This parsimonious approach reduces the resources required for navigation systems while holding a good accuracy. These methods are discussed through numerical results on an Angle-only navigation scenario. • Sampling period: ∆T = 1 s Number of angle measurements (2): 200 Standard deviation of the measurements: σ = 1/10 deg • Number of particles: N = 1000 Coefficient of the resampling threshold: c = 0.2 where N th = cN • Initial state X 0 : x 0 = 0, y 0 = 0, z 0 = 5000m, ẋ = 200m/s, ẏ = 200m/s, ż = 0m/s

I. INTRODUCTION

Angles-only navigation is a popular solution in robotics and transportation applications as it copes with widespread visual sensors or radio frequency receivers. Such navigation problems involve highly non-linear measurement equations and require a robust estimation process. Hence, Particle filtres appear as natural candidates. Recent works on the Laplace Particle Filter (LPF) which introduced accurate importance functions for resampling steps [START_REF] Musso | A laplacebased particle filter for track-before-detect[END_REF], [START_REF] Bui | Approximation particulaire et méthode de laplace pour le filtrage bayésien[END_REF] showed improved accuracy and precision with a reduced number of particles, making it a good choice for highly nonlinear navigation scenarios. However, angle-only navigation scenarios often provide a large number of landmarks, which can be very costly to process during the flight. For instance, Micro Aerial Vehicle navigation scenarios with stereo cameras [START_REF] Burri | The euroc micro aerial vehicle datasets[END_REF] give many features which require a large amount of signal processing resources. To that extent, a solution would be to select landmarks that are useful to the state estimation and discard those which are not. The idea introduced in this paper is to use a Cross-entropy optimization idea [START_REF] Reuven | The cross-entropy method: a unified approach to combinatorial optimization[END_REF], [START_REF] Dambreville | Cross-entropic learning of a machine for the decision in a partially observable universe[END_REF] to reduce the number of landmarks taken into account in the estimation process while maintaining good accuracy. The optimization is performed before the navigation and gives the subset of landmarks which is optimal with respect to a criterion on the Cramer-Rao Lower Bound (CRLB) detailed in the latter. After introducing the estimation problem section (II), the paper describes the principle of LPF and its adaptation for the angle-only navigation in section (III). Eventually, a specific modified Cross-entropy algorithm is proposed for landmarks selection optimization in section (IV).

II. PROBLEM STATEMENT

We consider an aircraft following a constant velocity motion. At time k, the state vector we aim to estimate is composed of the positions and the velocities:

X k = (x k , y k , z k , ẋk , ẏk , żk ) T . (1)
The aircraft measures the azimuth and the elevation of the aircraft to beacon line of sight (Fig. 1):

                   Y k (1) = arctan y k -y L x k -x L + ε k , = θ k + ε k Y k (2) = arctan z k -z L (x k -x L ) 2 + (y k -y L ) 2 + ν k = ϕ k + ν k , (2)
where ε k and ν k are zero-mean Gaussian noises with variance σ 2 k and are assumed to be mutually independent. The landmark coordinates (x L , y L , z L ) are assumed to be known.

X k = f k (X k-1 ) + η k , Y k = h k (X k ) + k , (3)
where f k and h k being possibly non-linear functions and where k , η k are-zero mean Gaussian noises. Nonlinear filtering aims to estimate the posterior density:

p k (x) P(X k = x|Y 1:k ), (4)
where

Y 1:k = (Y 1 , ..., Y k ).
A large variety of particle filters (PF) has been designed for this purpose (see tutorial [START_REF] Sanjeev | A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking[END_REF]). We recall the basics of the PF and Laplace Particle Filter (LPF) presented in [START_REF] Bui | Particle filtering and the laplace method for target tracking[END_REF]. As will be seen in the simulation (section III-D), the LPF is well adapted in such a context of angleonly filtering [START_REF] Bui | Multidimensional laplace formulas for nonlinear bayesian estimation[END_REF].

A. Generic particle filter

The PF estimates the state given all the measurements until time (k) through a Dirac mixture of N weighted particles X i k , w i k :

p k (x) ≈ N i=1 w i k δ x=X i k pk (x). (5)
The particles X i k evolve according to the dynamic model (3) and are corrected through the likelihood:

g k (x) P(Y k |X k = x). (6)
The following steps summarize the Sampling Importance Resampling algorithm (SIR) with prior proposal [START_REF] Neil J Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF]:

1) Initialisation (k=1) For i = 1, . . . , N . Generate the particles X i k-1 (positions and velocities) according to the prior with w i k-1 ≡ 1/N .

2) Prediction For i = 1, . . . , N . Propagate the particles by applying the dynamic model (3)

X i k|k-1 = f k (X i k-1 ) + η i k .
3) Correction For i = 1, . . . , N . Compute the likelihood (6) g k (X i k|k-1 ) = P(Y k |X i k|k-1 ) and the weights The corrected particles are:

w i k ∝ g k (X i k|k-1 ) w i k-1 such that N i=1 w i k = 1. Compute N ef f = 1 N i=1 [w i k ] 2 and the threshold N th = cN . * If N ef f ≥ N th . The corrected particles are X i k , w i k = X i k|k-1 , w i k * If N ef f < N th discard/
X i k , w i k = Xi k , 1 N .

4) State estimation

The state is estimated by

Xk = N i=1 w i k X i k .
Go to the prediction step k → k + 1

This PF algorithm can be significantly improved when considering importance functions (IF) which take into account the current measurement [START_REF] Godsill | Particle filtering: the first 25 years and beyond[END_REF] such as the Rao-Blackwellized particle filter (RBPF) [START_REF] Doucet | Rao-blackwellised particle filters for dynamic bayes nets[END_REF] or such as Sequential Markov Chain Monte Carlo (SMCMC) [START_REF] Brockwell | Sequentially interacting markov chain monte carlo methods[END_REF]. We propose here an adaptation of the LPF to the angle-only navigation. We recall briefly this algorithm, details are provided in [START_REF] Bui | Approximation particulaire et méthode de laplace pour le filtrage bayésien[END_REF], [START_REF] Bui | Multidimensional laplace formulas for nonlinear bayesian estimation[END_REF].

B. Laplace-based importance function

For the sake of clarity, in this section, the time index (k) is omitted. Consider an unknown d-dimensional state x distributed according to a prior q and observed through the measurement:

y = h(x) + , (7)
a non-linear function and a zero-mean Gaussian noise. In the particle filter framework, the prior q (which plays the role of the predicted density) is not expressed in closed form, only a sample from q is available. The likelihood p(y|x) is denoted by g(x), the posterior is written as:

P(X = x|Y = y) p(x|y) ∝ g(x) q(x). (8)
We aim to compute an estimate of the posterior. The importance sampling (IS) estimator p of p(x|y) is obtained by drawing N samples X i from a proposal distribution q so that:

p(x|y) ≈ N i=1 w i δ x=X i p(x), (9)
where

w i = wi N i=1
wi with the following importance weights:

wi = g(X i )q(X i ) q(X i ) . ( 10
)
It is well known that the optimal importance function is the posterior qopt (x) = p(x|y) which is unavailable. We choose a proposal, for instance a Gaussian, which has moments nearly equal to those of the posterior. The posterior expectation E[X|Y ] and the posterior covariance matrix V[X|Y ] are well approximated by the Laplace formula if the posterior has a predominant mode [START_REF] Bui | Multidimensional laplace formulas for nonlinear bayesian estimation[END_REF]. The following approximations are in general sufficient:

E[X|Y ] ≈ x, V[X|Y ] ≈ Ĵ-1 , (11)
where x is the maximum a posteriori (MAP), Ĵ = J(x) and where J is the observed information matrix:

   x = arg max x∈R d {g(x) q(x)}, J(x) = -(log g) (x) -(log q) (x). (12)
Assuming that the prior q(x) is Gaussian, J(x) becomes:

J(x) = -(log g) (x) + P -1 , ( 13
)
where P is the covariance matrix of the prior q. For example, we choose q as a Gaussian having x as mean and Ĵ-1 as covariance matrix:

q ∼ N (x, Ĵ-1 ). (14)

C. LPF adapted to angle-only navigation

We aim to approximate the MAP x defined as:

x = arg max x∈R 6 {g(x) q(x)}. (15)
Denote by x = (x 1 , x 2 ) T the state vector (1) where x 1 and x 2 are composed of, respectively, the positions and the velocities (1). The measurement function consists of the set of azimuth and elevation angles and depends solely on x 1 (2) (time index is still omitted). By noting that g(x) = p(y|x) = p(y|x 1 ), the maximization can be split as follows:

max (x1,x2) p(x 1 , x 2 |y) = max x1 q(x 1 ) p(y|x 1 ) max x2 q(x 2 |x 1 ) , (16)
where y is the measurement vector with covariance matrix R . Now, we assume that the prior q(x 1 , x 2 ) has a Gaussian distribution with mean (E 1 , E 2 ) T and with covariance matrix:

P = P 11 P 12 P 21 P 22 . (17)
The conditional pdf q(x 2 |x 1 ) being a Gaussian with conditional mean E 2|1 (x 1 ), the maximization of the second part of (16) with respect to x 2 leads to:

x2 (x 1 ) = E 2|1 (x 1 ) = E 2 + P 21 P -1 11 (x 1 -E 1 ) . (18)
Eventually, the initial maximization of dimension 6 (15) is reduced to a minimization on x 1 of dimension 3:

H(x 1 ) [x 1 -E 1 ] T P -1 11 [x 1 -E 1 ] (19)
+ [y -h(x 1 )] T R -1 [y -h(x 1 )].
One way of obtaining simply an approximation of the MAP is therefore as follows:

1) Minimize with respect to x 1 , the positions (19)

x1 = arg min x1∈R 3 H(x 1 )
2) Estimate x2 , the velocities (18)

x2 = E 2|1 (x 1 )
3) Evaluate x, the approximation of the MAP (15)

x = (x 1 , x2 ) T (20)
The LPF algorithm applied to the angles-only navigation is described below. It differs from the SIR particle filter only on the resampling step.

1) Initialisation (k=1) For i = 1, . . . , N . Generate the particles

X i k-1 = (x i 1,k-1 , x i 2,k-1
) (positions and velocities) according to the prior with w i k-1 ≡ 1/N . 2) Prediction For i = 1, . . . , N . Propagate the particles by applying the dynamic model ( 3)

X i k|k-1 = f k (X i k-1 ) + η i k .
Compute the mean and covariance of the predicted state Xk|k-1 and P k|k-1 .

3) Correction For i = 1, . . . , N . Compute the likelihood (6) g k (X i k|k-1 ) = P(Y k |X i k|k-1 ) and the weights

w i k ∝ g k (X i k|k-1 ) w i k-1 such that N i=1 w i k = 1. Compute N ef f = 1 N i=1 [w i k ] 2 and the threshold N th = cN . * If N ef f ≥ N th . The corrected particles are X i k , w i k = X i k|k-1 , w i k * If N ef f < N th then
( Xi k ) = φ( Xi k , Xk|k-1 , P k|k-1 )
where φ is a Gaussian pdf -Compute the importance weights:

wi k ∝ g k ( Xi k ) q k ( Xi k ) qk ( Xi k ) such that N i=1 wi k = 1. The corrected particles are: X i k , w i k = Xi k , wi k .

4) State estimation

The state is estimated by

Xk = N i=1 w i k X i k
Go to the prediction step k → k + 1. Remark: using the predicted covariance matrix P k|k-1 instead of Ĵ-1 k (14) could provide a better robustness.

D. Simulations

We consider an aircraft following a constant velocity noisefree motion. The simulation parameters are the following: Initial uncertainty state covariance : P 0 = diag(3000m, 3000m, 100m, 20m/s, 20m/s, 20m/s) 2 . Position particles initialization are generated around the MAP obtained by the first triangulation (using 2 landmarks) with the Cramer-Rao Lower Bound (CRLB) evaluated at this MAP as covariance matrix.

• Angle measurements are provided using successively 3 pairs of landmarks at altitude 0. The x -y landmark positions, displayed in km, are as follows:

-for k ∈ [1s, 60s]:

L 1 = (1, 15), L 2 = (-1, 15) -for k ∈ [61s, 149s]: L 1 = (5, 20), L 2 = (-5, 25) -for k ∈ [150s, 200s]: L 1 = (50, 30), L 2 = (
20, 50) 100 Monte Carlo trials have been performed. The root mean square error (RMSE) is computed for the 6 components of the state vector X k . The error averages are compared to the CRLB, the inverse of the information matrix J k . The information matrix of the azimuth and elevation angles relative to X k up to time k using L(k) landmarks at time k (here L(k) = 2) can be expressed recursively as follows:

J k = L(k) j=1 1 σ 2 k ∂θ j k ∂X k ∂θ j k ∂X k T (21) + L(k) j=1 1 σ 2 k ∂ϕ j k ∂X k ∂ϕ j k ∂X k T + Φ -1 (k -1, k) T J k-1 Φ -1 (k -1, k),
with J 0 = P -1 0 being the initial information matrix, and

Φ(k -1, k) = ∂X k
∂X k-1 the transition matrix. The filtering errors are close to those of the CRLB (Fig. 3,4). The LPF resampling rate is about 10 %, the computational time is low, about 2s for 200 iterations on a 2,6 GHz Intel Core i5.

A. Introduction

The signal process to provide the angles from the landmarks is energy and time-consumming. Therefore, we aim to select a subset of the landmarks without degrading the navigation performance too much. To evaluate the impact of the location of the landmarks on the estimated position, we consider the following simple example (Fig. 5). Let L 1 and L 2 be the landmarks. The triangulation is done using the azimuths θ 1 , θ 2 with variance σ 2 . P r (Aircraft) is the orthogonal projection of the aircraft on the plane. The following approximation, valid for relatively small values of |B/r|, gives the order of magnitude of the error of estimation of |r| and therefore of the error of the position estimation:

σ( | r|) ≈ √ 2 σ |r| 2 |B| | sin α| . (22)
Indeed, we have r 1 = r + B/2 and r 2 = r -B/2.

Using the complex number notation,

r 1 = |r 1 | e i( π 2 -θ1) , r 2 = |r 2 | e i( π 2 -θ2
) , a taylor expansion gives:

                     θ 1 = π/2 -Im [log(r 1 )] ≈ π/2 -Im log(r) + B 2r - 1 2 B 2r 2 , θ 2 = π/2 -Im [log(r 2 )] ≈ π/2 -Im log(r) - B 2r - 1 2 B 2r 2 ,
where "Im" denotes the imaginary part. We obtain now:

θ 2 -θ 1 = Im B r + o |B/r| 3 = |B/r| sin(α) + o |B/r| 3 .
Thus, the variance of the inverse of the estimation of |r| is given by:

V 1 | r| ≈ 2σ 2 |B| 2 sin 2 (α)
, leading to the desired standard deviation of the position in the x-y plane (22) observing that:

V f ( | r|) ≈ f 2 (r) V ( | r|)
where f (r) = 1 r . We note (22) that the Std of | r| is proportional to the square of the distance from the middle of landmark positions and P r (Aircraft). The same is true for the x -y position errors. The approximation (22) suggests that the landmarks should be chosen in such a way to ensure (22):

• a small distance |r| • an angle α close to ±π/2 • a baseline B large enough (Fig. 5) In the next section, the simulation results of the landmarks selection optimization will reflect the trade-off between these 3 criteria.

B. Cross entropy

Subsequently, we address the optimization of landmarks use in order to maximize the localization accuracy. Our idea is to restrict the number of landmarks in use at each step, and then optimize the landmarks choice trajectory. More precisely, if Q is the total number of landmarks, M is the number of selected landmarks at each measure stage and K is the number of measure stage, then we have to optimize a trajectory in ([[1, Q]] M ) K which describe the landmarks selection for our localization process. Here, we have to deal with a discrete combinatorial optimization problem, which is solved using a metaheuristic. Our approach is inspired from the cross-entropy (CE) method [START_REF] Reuven | The cross-entropy method: a unified approach to combinatorial optimization[END_REF], [START_REF] Dambreville | Cross-entropic learning of a machine for the decision in a partially observable universe[END_REF]. Consider a function f : x ∈ E → f (x) ∈ R to be maximized. Assume a law family (π λ ) λ∈Λ defined over E. In particular, it is assumed that this family is sufficiently expressive to approximate a Dirac around any element of E. Let ρ ∈ [0, 1[ be a smoothing parameter. The cross-entropy method is an optimization algorithm that will tune the family parameter to approximate an optimum of function f . It is broadly based on the following principle:

• Initialize λ 0 ∈ Λ and set t = 0.

• Repeat the steps until optimal criterion is achieved:

1) Generate N samples ( x 1:N ) from law π λt , 2) Compute the evaluations f ( x n ) for 1 ≤ n ≤ N , 3) Assign weights to the samples in accordance with the evaluation. It is obtained a particles cloud ( x n , w n )| 1≤n≤N where w ≥ 0 and While general formulation (23) makes sense intuitively, its mathematical formalization is not clear when dealing with Dirac distribution. But, by splitting the logarithm within D( | ) and discarding the constant part of the expression, we are able to reduce (23) to mathematically sound maximizations:

• If π λ = I i=1 λ i δ xi with λ = (λ i ) ∈ R I
+ and i λ i = 1 (discrete case), then there is i n such that x n = x in by construction and minimization ( 23) is equivalent to:

max λ ρ I i=1 λ t,i ln(λ i ) + (1 -ρ) N n=1
w n ln(λ in ) , (24)

• If π λ is a density, minimization (23) is equivalent to:

max λ ρ π λt (x) ln(π λ (x))dx + (1 -ρ) N n=1 w n ln(π λ ( x n )) .
(25) This algorithm is somehow related to evolutionary algorithms, but the population distribution is characterized by the law parameter. The samples weights w n are chosen so as to select the best samples in regards to f . In [START_REF] Reuven | The cross-entropy method: a unified approach to combinatorial optimization[END_REF], only a quantile of the samples is selected, so that w n is uniform on the best quantile and is zero otherwise. However, in our implementation, the samples are sorted in regards to their evaluations by f and the weight of a sample is defined in proportion to its sorting rank.

C. Optimization

The criterion f to be maximized is based on the information matrix. Subsequently, f is defined as the 12th root of the determinant of the last Fisher information matrix computed recursively by (21):

• 12
2 is actually the size of the information matrix, so that f is comparable to a distance. • The information matrix is related to the last step of the time interval of the whole filtering process. Thus, this is a prior optimization, and not a dynamic optimization. The landmarks are available all the time; this paper does not deal with randomly intermittent landmarks, although this difficulty may be considered for future works. The variable to be optimized is the choice of landmarks subsets to be observed at each iteration of the localization process. Duplicate landmarks for the same iteration stage is simply omitted in the recursive Fisher computation. We intend then to evaluate the overall performance of localization for different cardinality bounds, M , on the number of landmarks in use.

In practice, we used N = 1000 particles in our optimization and ρ = 0.95 as a smoothing parameter. Each particle is a landmark trajectory (q m,s ) ∈ ([[1, Q]] M ) S where Q = 35 is the total number of landmarks, M is the number of selected landmarks at each measure stage and K = 10S is the number of iterations of our localization process: it is assumed that the same set of landmarks is used during a 10-length period. As a consequence, measure stage k and landmark stage s are related by relation s = . We used a family based on discrete Markov laws in order to generate the landmark trajectories:

1) Landmarks for s = 1 are generated using a sequence of conditional laws, 2) i-th landmark of stage s + 1 ∈ [[2, S]] is generated from i-th landmark of stage s using a conditional law. Family (π λ ) λ∈Λ is thus discrete and expressed by:

π λ (q) = λ(q 1,1 ) M m=2 λ(q m,1 |q m-1,1 ) S s=2 λ(q m,s |q m,s-1 ) , (26)
where parameters q 1,1 → λ(q 1,1 ), q m,1 → λ(q m,1 |q m-1,1 ) and q m,s → λ(q m,s |q m,s-1 ) are discrete probabilities. Update of the law by means of (24) is rather easy, since the optimization is entirely separable owing to the Markovian nature of (26):

λ o (q 1,1 ) = ρλ t (q 1,1 ) + (1 -ρ) N n=1 w n δ q1,1= qn,1,1 N n=1 w n , (27) λ o (q m,1 |q m-1,1 ) ∝ ρλ t (q m-1,1 )λ t (q m,1 |q m-1,1 ) + (1 -ρ) N n=1 I qn,m-1,1=qm-1,1 w n δ qm,1= qn,m,1 N n=1 w n , (28) λ o (q m,s |q m,s-1 ) ∝ ρλ t (q m,s-1 )λ t (q m,s |q m,s-1 ) + (1 -ρ) N n=1 I qn,m,s-1=qm,s-1 w n δ qm,s= qn,m,s N n=1 w n , (29)
where I true = 1 -I false = 1 and where λ t (q m-1,1 ) and λ t (q m,s-1 ) are computed from Markovian definition (26). In this implementation, the smoothing parameter ρ is included within divergence optimization (23). In the reference implementation of CE [START_REF] Reuven | The cross-entropy method: a unified approach to combinatorial optimization[END_REF], this smoothing parameter is applied afterward, that is on the parameter vector λ. To the benefit of our approach, the update is able to adapt optimally to the number of samples that actually match a condition. For example, the update process will produce λ o (q m,s |q m,s-1 ) = λ t (q m,s |q m,s-1 ) when there is no sample such that q n,m,s-1 = q m,s-1 . The CE process is stopped after 5000 unsuccessful steps.

D. Simulations

In our experiment, landmarks have been generated randomly at some distance along the trajectory and K = 100 observation stages (S = 10) are processed. Figure 6 presents the mean trajectory of the aircraft, from (0km, 0km, 10km) to (10km, 10km, 10km), together with 35 landmarks (dots in black). The initial velocities are (200m/s, 200m/s, 0m/s). The initial uncertainty state covariance is still:

P 0 = diag(3000m, 3000m, 100m, 20m/s, 20m/s, 20m/s) 2 .
Still in figure 6, the optimal landmarks selection is presented when M = 5 and for stages [ [START_REF] Musso | A laplacebased particle filter for track-before-detect[END_REF][START_REF] Godsill | Particle filtering: the first 25 years and beyond[END_REF] best is the optimal found value, and iter is the number of needed iterations. In this landmark configuration with a lot of redundancy, a rather good level of performance is already obtained after the first 1000 samples generation. The optimal performances need however some substantial amount of iterations within the CE algorithm. We have in mind some approaches to speed up this convergence based on geometric criteria (22).

Figure 7 presents the position standard deviations obtained after optimization. Figure 8 presents the velocity standard deviations obtained after optimization. The curves are presented with regards to M , the number of used landmarks at each step. These figures together with the previous table confirm that the performance improvement rate will decrease with the number of used landmarks.

E. Future works

Navigation systems often estimate the attitude of the aircraft in addition to the position and the velocity. However, adding the vehicle's attitude to the estimation process is expected to complexify the optimization and degrade the triangulation's performance. The impact of attitude uncertainties on triangulation can be described via the study of the position Cramer Rao Bound. Let the state vector X k composed of the position and the attitude:

X k = (x k , y k , z k , ψ 1 , ψ 2 , ψ 3 ) T , (30)
where (ψ 1 , ψ 2 , ψ 3 ) are respectively the precession, nutation and spin Euler angles. The velocity is overlooked here as the vehicle is assumed stationary. This augmented state vector follows the model (3). Focusing on the update step, the measurement model becomes for the i th observed landmark:

                     Y i k (1) = arctan ∆ i y,k ∆ i x,k + k = θ i k + k Y i k (2) = arctan   -∆ i k,z (∆ i x,k ) 2 + (∆ i y,k ) 2   + ν k = ϕ i k + ν k , (31)
where ε k and ν k are zero-mean Gaussian noises with variance σ 2 k . Also, ∆ i k denotes the relative distance between the i th landmark and the aircraft resolved in the aircraft frame:

∆ i k =   ∆ i x,k ∆ y,k ∆ i z,k   = C b e,k     x i L y i L z i L   -   x k y k z k     , (32)
and C b e,k denotes the rotation matrix from the ground frame towards the vehicle frame. It is obtained from the Euler angles (ψ 1,k , ψ 2,k , ψ 3,k ). From now on, the measurements for 2 landmarks are such that:

Y k ∼ N (h(X k ), R k ) h(X k ) = [θ 1 k , ϕ 1 k , θ 2 k , ϕ 2 k ] T , (33)
where

Y k = Y 1 k (1), Y 1 k (2), Y 2 k (1), Y 2 k (2)
T (31) and where

R k = σ 2 k I 4 .
The prior law of the state vector is the following:

X k ∼ N ( Xk , P k ), (34)
where P k is a diagonal matrix such that:

P k = I 3 σ 2 attitude 0 0 I 3 σ 2 position . ( 35
)
The impact of σ attitude in the prior density is studied thanks to the CRLB of the norm of the position error:

ρ(x k , ŷk , ẑk ) = (x k -x k ) 2 + (ŷ k -y k ) 2 + (ẑ k -z k ) 2 , ( 36
)
The CRLB is defined as the inverse of the information matrix:

CRLB = J k (X k ) -1 . (37)
According to (21) with L(k) = 2 and Φ = I 6 , J k is given by:

J k (X k ) = 1 σ 2 k ∂h(X k ) ∂X k ∂h(X k ) ∂X k T + P -1 k . (38)
Then, the variance of ρ from the CRLB is computed as follows:

V(ρ(X k )) = ∂ρ(X k ) T ∂X k J k (X k ) -1 ∂ρ(X k ) ∂X k . ( 39
)
Simulations are run on a static triangulation problem from two landmarks. For the sake of simplicity, the landmarks are symmetric with respect to the (x, z) plane and the prior error on position is σ position = 100m. As displayed (Fig. 9), the prior error on attitude sharply increases the CRLB before reaching an asymptote, showing that the performance of the estimation problem strongly depends on a correct estimation of the vehicle attitude. In order to broaden the results, similar simulations were run with a variable distance between the fixed points. The second plot shows the impact of σ attitude on the position estimation error as a function of the landmarks spacing. Hence, taking into account the attitude error in the landmarks optimization process would complete the approach described in this paper and will be the object of future work.

V. CONCLUSION

The paper proposes an adapted version of the Laplace Particle Filter for the angle-only navigation using a set of landmarks. Some simulation results using 3 pairs of landmarks show the accuracy of this algorithm by comparing its performance with the Cramer Rao Lower Bound. We propose also an optimization algorithm for the selection of the landmarks. This algorithm is based on the idea of the Cross-entropy method with more complex Markov simulation laws. Simulation results show that one can reduce drastically the number of landmarks without degrading significantly the navigation performance. Lastly, the impact of attitude uncertainties on triangulation is analyzed. Further works will take into account the attitude angles uncertainties for the landmarks selection optimization.

Fig. 1 :

 1 Fig. 1: Illustration of the angle of arrival measurements from a fixed landmark

Fig. 2 :

 2 Fig. 2: Illustration of simultaneous measurements from a set of landmarks

Fig. 3 :

 3 Fig. 3: Position errors (m)

Fig. 4 :

 4 Fig. 4: Velocity errors (m/s)

1 Fig. 5 :

 15 Fig. 5: Triangulation in the x-y plane

N n=1 w n = 1 , 4 )

 14 Let δ x= xn be Dirac distribution on x n . Then, compute: λ o ∈ arg min λ∈Λ D ρπ λt + (1 -ρ) N n=1 w n δ x= xn π λ , (23) so as to minimize the Kullback-Leiber divergence D( | ) of π λ with the mixture ρπ λt + (1 -ρ) N n=1 w n δ x= xn ; the smoothing parameter ρ is here to tune down the innovation implied by the selected samples, 5) Set λ t+1 = λ o and t ← t + 1.

  k+9

10

 10

  ] (red), [[31, 40]] (blue), [[61, 70]] (green) and [[91, 100]] (cyan). It is noteworthy that the optimizer favors the choice of the closest landmarks while keeping a large baseline.

Fig. 6 :

 6 Fig. 6: Optimal landmarks selection at different stages Subsequent table presents some convergence statistics for M ∈ [[1, 10]] and should be compared to the value f [35] = 8.13 obtained for M = 35: M 1 2 3 4 5 6 7 8 9 10 f [M ] 1000 1.49 2.07 2.51 2.86 3.22 3.45 3.70 3.95 4.14 4.31 f [M ] best 1.60 2.25 2.76 3.18 3.53 3.85 4.13 4.39 4.63 4.72 iter 1000 7.9 12 18 22 27 32 38 41 44 26

Fig. 7 :

 7 Fig. 7: Std of the positions: σ x , σ y , σ z (m). Curves for M = 1 : 10.

Fig. 8 :

 8 Fig. 8: Std of the velocities: σ ẋ , σ ẏ , σ ż (m/s). Curves for M = 1 : 10.

Fig. 9 :

 9 Fig.9: The first plot shows the CRLB for the position norm as a function of the prior standard deviation on the attitude. The second plot shows CRLB for the position norm as a function of the distance between the two landmarks, for various prior standard deviations on the attitude.


Adel Djoudi

Martin Hána

Nikolai Kosmatov

Milan Kříženecký

Franck Ohayon

Patricia Mouy

Arnaud Fontaine

David Féliot

A Bottom-Up Formal Verification Approach for Common Criteria Certification: Application to JavaCard Virtual Machine

I. INTRODUCTION

The quality and security of critical software have become nowadays a major concern. The Common Criteria (CC) for Information Technology Security Evaluation [1] provide an international standard for computer security certification. Its highest assurance levels EAL6-EAL7 require a formal Security Policy Model (SPM) and an associated mathematical proof of security properties (i.e. confidentiality, integrity). Thales recently conducted a formal verification of a JavaCard platform module [START_REF] Oracle | Java Card Platform: Virtual Machine Specification, Classic Edition[END_REF] in a novel EAL6 certification project of a smart card product. This certification project was evaluated by CEA-Leti (an evaluation center, or ITSEF) with the supervision of ANSSI (the French national cybersecurity agency, and the French certification body).

Historically, since the verification of real-life code was not feasible for large industrial projects, the certification usually followed a top-down approach, where a separate abstract model was used to verify the specified security properties, and then refined to the code. A classical approach of applying formal verification on a JavaCard platform consists in building a high-level formal model of target sub-modules. The need to bridge the gap between the formal model and the implementation and to provide stronger guarantees for the real-life code was reported by experts [START_REF] Beckert | Mind the gap: Formal verification and the common criteria (discussion paper)[END_REF].

In our work, we adopt a novel bottom-up methodology relying on verification of the real-life code of a JavaCard Virtual Machine using the Frama-C verification platform [START_REF] Kirchner | Frama-C: A software analysis perspective[END_REF]. We expressed all specified features and properties as annotations in a formal specification language, called ACSL (ANSI C Specification Language), inserted directly in the source code. We annotated over 7,000 lines of C code in ACSL, and over 50,000 proof goals were generated and formally proved by the tool. An earlier paper [START_REF] Djoudi | Formal verification of a JavaCard virtual machine with Frama-C[END_REF] focused on technical and scalability issues of the proof without addressing the certification methodology. In this paper, we focus on methodology aspects: we describe this bottom-up approach, discuss its benefits and challenges and compare it to previous top-down approaches.

II. COMMON CRITERIA CERTIFICATION PROCESS

A. Overview of Common Criteria Evaluation

The international standard ISO/IEC 15408-which defines the Common Criteria (CC)-is an international agreement on security evaluation of IT products. It contains a Common Evaluation Methodology (CEM) describing the general evaluation process from EAL1 up to EAL5. SOG-IS 1 is the European mutual recognition agreement that was concluded in 2010 and involves ten countries. For EAL6 and EAL7, formal methods have to be used, but the CEM does not detail how to use them to demonstrate the respect of the security properties. For this reason, the CEM is completed by an additional interpretation by the National Certification Body. Within the SOG-IS, only three countries mutually emit and recognize certificates up to EAL7: France, Germany and Netherlands. Each of them defines its own interpretation for EAL6 and EAL7 (AIS-34 in Germany and Note 12 [START_REF] Sgdn | Application note: Target evaluation's security policies formal modelling[END_REF] in France). Recognition agreements beyond the EU (between more than 30 countries) are defined in the Common Criteria Recognition Arrangement (CCRA).

A Common Criteria evaluation is initiated by the owner of the product to evaluate. The product owner (generally, an industrial) establishes a contract with an approved evaluation center (ITSEF) and registers the evaluation with the Certification Body. At the end of the evaluation, an Evaluation Technical Report (ETR) is produced by the ITSEF to be reviewed by the Certification Body.

The national certification bodies recognized to deliver EAL6-EAL7 certificates do not prescribe the use of any particular formal method or tool [START_REF] Lavatelli | Formal models for high assurance: why and how[END_REF]. Up to now, only methods relying on high-level models in B and Coq were recognized by ANSSI in the Common Criteria context. A guidance document about the use of formal methods was published by ANSSI [START_REF] Sgdn | Remarques relatives à l'emploi des méthodes formelles (déductives) en sécurité des systèmes d'information[END_REF]. A more recent guidance written by ANSSI and INRIA research teams was published but with a focus on Coq [START_REF] Anssi-Inria | Requirements on the use of Coq in the context of common criteria evaluations[END_REF] with a dedicated paper [START_REF] Bertot | The use of Coq for common criteria evaluations[END_REF], which details the rationale of the guidelines and requirements from ANSSI. The respect of these guidelines has to be verified during the evaluation process.

The adoption of a new formal method or tool requires a pilot evaluation. It implies a tripartite work between the developer, the ITSEF and the certification body, and additional effort to inspect formal assurance.

For an evaluation, the developer has to supply the following evidence:

• the source code of the Formal Security Model (SPM),

• an explanatory document presenting a description of the model, a complete list of the associated hypotheses (explicit ones but also implicit ones due to modeling choices), their justification and their consistency, • explicit links between the model, the security target and the security properties (further discussed below),

• a clear justification of the level of confidence for the method and tools used, and • all the necessary information to allow the evaluator to reproduce the proof(s).

B. Security Specification

A Common Criteria certification of a product helps the customer to determine whether the security of a product is sufficient to meet their needs and to ensure that the security properties are satisfied. For the product owner, the Common Criteria help to identify security issues, define security objectives, establish security requirements relying on a standardized catalog and then define a precise Target of Evaluation (TOE) that is usually only part of the entire IT product. The security specification plays an important role in the certification process.

The Common Criteria offer a catalog of Security Functional Requirements (SFRs) and Security Assurance Requirements (SARs) [1]. The SFRs define the TOE security characteristics. The SARs define confidence degree in the enforcement of the security objectives of the TOE. The assurance level is increased by increasing the scope, depth and rigor of the evaluation effort in six assurance classes: Security Target Evaluation (ASE), Development (ADV), Guidance (AGD), Life Cycle Support (ALC), Tests (ATE) and Vulnerability Assessment (AVA). At EAL6-EAL7 levels the ADV assurance requirement class mandates to build a formal security policy model (SPM) as the most rigorous way to identify and eliminate ambiguous, inconsistent, unenforceable or contradictory security policy elements [1]. For instance, the Common Criteria Action Element ADV SPM.1.1D mandates the developer to identify the security policies that should be formally modeled. Note that the National Certification Body provides guidance for the interpretation of such statements [START_REF] Sgdn | Note d'application : Modélisation formelle des politiques de sécurité d'une cible d'évaluation[END_REF] in order to satisfy the target security objectives. The identification of the security policy implies a list of SFRs to be formally modeled.

C. Application to JavaCard

JavaCard system is a well-known security-critical product. Many JavaCard products have been subject to Common Criteria evaluation. The Security Target Evaluation class (ASE) enforces the definition of a Security Target (ST) for an identified product (e.g. a particular JavaCard product). A Security Target is an implementation-dependent statement of security needs. It states what is to be evaluated before the evaluation is performed (and thus helps to understand after the evaluation what was actually evaluated). The Security Target may claim conformance (strict or demonstrable) to a Protection Profile for a generic TOE type (such as a JavaCard system) [START_REF] Oracle | Java Card system -open configuration protection profile, version 3.1[END_REF]. A protection profile provides a standardized statement of Security Policies to be tailored according to the defined scope of evaluation in Security Targets. Over the years, the JavaCard System Protection Profile has established as one of the most recognized smartcard industry reference and is typically mandated in tenders or requested explicitly by customers. For instance, the Open Configuration Protection profile of JavaCard systems defines the Firewall Security policy/aspect [START_REF] Oracle | Java Card system -open configuration protection profile, version 3.1[END_REF] as follows:

"#.FIREWALL: The Firewall shall ensure controlled sharing of class instances, and isolation of their data and code between CAP files (that is, controlled execution contexts) as well as between CAP files and the JCRE context. An applet shall not read, write, compare a piece of data belonging to an applet that is not in the same context, or execute one of the methods of an applet in another context without its authorization." The protection profile also instantiates this security aspect with a security objective (O.Firewall) and provides a rationale for the list of SFRs to be satisfied in order to meet this objective. The Security Target of a product may then instantiate the Common Criteria requirement ADV SPM.1.1D as follows:

"ADV SPM.1.1D: The developer shall provide a formal security policy model for the Firewall Security Policy." and provide a rationale (according to Common Criteria security components ASE OBJ and ASE REQ) for the (sub)set of SFRs that are formally modeled in order to meet this requirement. In the sequel of this document, we consider the two following (simplified) SFRs for illustration purposes:

SFR1 The Target of Evaluation Security Functions shall enforce the Firewall access control policy to provide restrictive default values for security attributes that are used to enforce the security policy. The objects' security attributes of the access control policy are created and initialized at the creation of the objects. Afterwards, these attributes are no longer mutable. SFR2 The Target of Evaluation Security Functions shall enforce the following rule to determine if an operation among controlled objects is allowed: the Currently Active Context may freely perform any memory access operation upon any object whose Lifetime attribute has value "Persistent" only if the object's owning Context attribute has the same value as the Currently Active Context. Further details about the interpretation of these SFRs are provided in the following sections. In particular, we illustrate the mapping of these SFRs to our formal model in Section V. JavaCard applets are compiled to bytecode, which is executable in the JavaCard Virtual Machine (JCVM). A binary file (called CAP file), loadable to the platform, encapsulates mainly the bytecode together with class definitions. It may contain several Java packages and applets. A unique context is associated to each CAP file during loading to the card. Prior to loading a CAP file, a ByteCode Verifier (BCV) is run off-card to perform a static analysis (typelevel abstract interpretation) of the applets [START_REF] Leroy | Bytecode verification on Java smart cards[END_REF]. This ensures that the code does not attempt to perform ill-typed operations that may bypass security protections ensured by the JCVM. Indeed, the virtual machine ensures bytecode interpretation and offers higherlevel, more secure abstractions of data than the hardware processor, such as object references instead of memory addresses. Memory access operations on the stack, heap or code area are performed by the JCVM for each interpreted bytecode.

III. BOTTOM-UP APPROACH BASED ON DEDUCTIVE VERIFICATION

A. JavaCard Virtual Machine

The JCVM firewall enforces runtime-protection of applet security properties against major security concerns: developer mistakes and design oversights allowing sensitive data "leakage" from the applet owning the data to another applet without explicit permissions [START_REF] Oracle | Java Card Platform: Runtime Environment Specification, Classic Edition, Version 3.1[END_REF]. In general, the firewall blocks access to the data in one CAP file to an applet in another CAP file (having different execution contexts), except some well-defined exceptions (such as global arrays, ArrayViews or class variables).

B. Target of Evaluation

The target of evaluation (TOE) of our project includes a set of Security Functions (TSF) that support the O.Firewall Security Objective enforcement. The implementation of these Security Functions in our project consists in a large subset of C functions of the JCVM. This subset features all possible functionalities of a single JCVM run that ensures the execution of any applet being selected. The JCVM specification [START_REF] Oracle | Java Card Platform: Virtual Machine Specification, Classic Edition[END_REF] describes the data life cycle of an applet: transient deselect data is erased during owning application deselection, transient reset data is erased only when the smart card is reset, and persistent data is preserved anytime. Data and the associated life cycle of an undefined number of applets are carefully considered in our Security Policy Model (SPM) in order to ensure isolation properties. The SPM is deemed to comply with security functional requirements (SFRs) of the O.Firewall security objective defined in the JavaCard System Protection Profile [START_REF] Oracle | Java Card system -open configuration protection profile, version 3.1[END_REF], according to which "an applet shall not read, write, compare a piece of data belonging to an applet that is not in the same context, or execute one of the methods of an applet in another context without its authorization".

C. Formal Security Policy Model (SPM)

Before this work, ANSSI's requirements for formal methods in the context of Common Criteria evaluations were only successfully fulfilled using B and Coq methods. This project was part of a pilot evaluation due to the adoption of a new formal method based on Frama-C [START_REF] Kirchner | Frama-C: A software analysis perspective[END_REF], a verification platform for C code. Frama-C uses ACSL (ANSI C Specification Language) [START_REF] Baudin | ACSL: ANSI/ISO C Specification Language[END_REF], a formal specification language for C programs. It allows the user to specify annotations that express the expected program properties. The Frama-C/WP plugin can then be used to prove them for each function of the code: this technique is called (modular) deductive verification. In order to build the formal security policy model (SPM), we follow a bottom-up approach, in which the C code implementation is enriched with annotations instead of merely making a separate model and a formal representation of a set of SFRs being claimed (see Fig. 1). The security properties are proved to be enforced by the actual implementation design. Examples of ACSL annotations and security properties are given in the next section.

IV. FORMAL SPECIFICATION OF CONTRACTS AND SECURITY PROPERTIES

A. Function Contracts

Deductive verification with Frama-C/WP requires that each C function be annotated with a (function) contract expressed in ACSL language. Such a contract contains preconditions (in requires clauses), which express properties of program variables that must be respected before the function is called, and postconditions (in ensures clauses), which express properties that must be ensured after the function terminates. A special kind of postconditions is expressed by assigns clauses, which give the list of variables that the function is allowed to modify. All other variables cannot be modified by the function. Each function should then be proved by Frama-C/WP to respect its contract. All functions in our Formal Security Policy Model are annotated by ACSL contracts (see [START_REF] Djoudi | Formal verification of a JavaCard virtual machine with Frama-C[END_REF] for detailed examples). For lack of space, we illustrate an ACSL contract in Fig. 2 for a very simple C function swap that swaps the values * p1 and * p2 referred by the two pointers given as function arguments. Line 1 expresses a precondition stating that input pointers p1 and p2 must refer to valid memory locations, that is, locations that can be safely accessed, and that these locations are separared, that is, the underlying bytes are disjoint. The validity is necessary here to ensure the absence of runtime errors, also known as undefined behaviors. The separation assumption is necessary to avoid overwriting some bytes of * p2 when modifying * p1, and vice versa. one to refer to the value of expression e before the call. Finally, line 3 states that * p1 and * p2 are the only memory locations the function is allowed to modify.

B. Security Properties Expressed as Predicates

Some security properties (typically, for integrity) can be specified in ACSL as invariant properties maintained by relevant functions and directly proved by Frama-C/WP. We illustrate one of such properties, expressed by predicate object_headers_intact shown in Fig. 3, stating that object headers are not modified between program points L1 and L2. For lack of space, we give here only the main ideas of this definition and refer the reader to [START_REF] Djoudi | Formal verification of a JavaCard virtual machine with Frama-C[END_REF] for detailed explanations. This predicate states that for any allocated object (represented by its index i) and for any offset off within the object header, the byte at offset off in the object header of object i has the same value in program points L1 and L2. A typical usage of this predicate is to include the postcondition ensures object_headers_intact{Pre, Post}; in the contract of every function. This postcondition ensures that the object headers are not modified by the function between program points Pre and Post, which refer to the states before and after the function call.

C. Security Properties Expressed as Metaproperties

Other properties (for confidentiality or some cases of integrity) are stated as metaproperties [START_REF] Robles | Methodology for specification and verification of high-level properties with MetAcsl[END_REF]. The main principle of a metaproperty is to state a global property for a specified set of target functions and a specified context. Two examples of key metaproperties are shown in Fig. 4. Again, we give here the main ideas of their definition, a more detailed presentation being available in [START_REF] Djoudi | Formal verification of a JavaCard virtual machine with Frama-C[END_REF]. Metaproperty persi_objects_integrity states that the data of a persistent object (represented by its index i) cannot be modified unless the current context (JCC) is the object owner (which is stored at offset 0 in the object header). Similarly, metaproperty persi_objects_confidentiality states that the data of a persistent object (represented by its index i) cannot be read unless the current context is the object owner.

Each metaproperty is instantiated by the Frama-C/MetAcsl plugin into assertions in relevant program points in all target functions. For example, the first metaproperty has the writing context, therefore the corresponding property must be checked each time when a memory location is modified. So an assertion is automatically added by the Frama-C/MetAcsl plugin at all those memory locations, where the metavariable \written is replaced by the written memory location. The proof of the resulting assertions ensures that the metaproperty is globally respected by the code.

The proof of real-life code in our project requires a careful combination of several ingredients (see Fig. 6): macros, companion ghost code, global preservation properties in addition to lemmas and proof scripts. Macros reduce redundancy in specifications and facilitate updates and maintenance. Ghost code is mainly used to describe the memory model and to offer an alternative encoding of low-level operations, amenable to automatic provers. This combination made it possible to efficiently reason about non-trivial code fragments involving bitwise operations without the use of external interactive tools (e.g. Coq) with a high level of automatic proof.

V. TRACEABILITY BETWEEN CC REQUIREMENTS AND IMPLEMENTATION a) General Approach: As described in Section II, Security Objectives and related Security Functional Requirements are summarized in a Protection Profile describing a particular product type. In this project, the Security Target is related to the Protection Profile for JavaCard System [START_REF] Oracle | Java Card system -open configuration protection profile, version 3.1[END_REF]. To establish a correspondence between formal (SPM) and informal concepts (ST), the developer must establish and describe the links between them, as mandated by [START_REF] Sgdn | Note d'application : Modélisation formelle des politiques de sécurité d'une cible d'évaluation[END_REF]. In fact, ANSSI and BSI (the German national certification authority) have driven the use of formal methods in Common Criteria evaluations with the publication of guidance for developers and evaluators (Note-12 [START_REF] Sgdn | Note d'application : Modélisation formelle des politiques de sécurité d'une cible d'évaluation[END_REF] and AIS34 [START_REF] Bsi | Application notes and interpretation of the scheme (AIS)[END_REF]). Formal analyses in CC context consist in giving a proof that the TOE Security Functions correctly implement the expected security objectives. Several "representations" of the TSFs are provided as shown in Fig. 5: The CC Assurance Development (ADV) components : Security Policy Model (ADV SPM), Functional Specification (ADV FSP), Target of evaluation Design (ADV TDS), Implementation representation (ADV IMP) provide a list of requirements to be fulfilled by each of the Security Function representations. Developers are also expected to establish a "formal equivalence" of these various representations of the TSF. The primary objective is to formally establish the correctness of the SPM w.r.t. security objectives. The rationale is to authorize reasoning at an abstract level (SPM) and to propagate the result toward the implementation (IMP).

The lack of proof of refinement until the implementation is the rub with the top-down strategy [START_REF] Chetali | About the world-first smart card certificate with EAL7 formal assurances[END_REF], [START_REF]Industrial use of formal methods for a high-level security evaluation[END_REF], [START_REF] Bolignano | Formal proof of a secure OS full trusted computing base[END_REF]. This is all the more unfortunate since most of the time the formal model is provided a posteriori, only for CC evaluation, and is not used to guide the design of the implementation.

Since our verification approach is based on formal properties written in ACSL annotations directly on the implementation, we link all informal concepts to formal ACSL annotations. Design (ADV TDS), Security Policy Model (ADV SPM), Functional Specification (ADV FSP)) are structured with respect to the actual product source code.

The proposed bottom-up approach intrinsically encompasses the refinement from the functional specification through the design to the implementation that is usually required in top-down methodologies.

b) Mapping Security Objectives and Requirements: The Security Objective defined in the Security Target (in our case, O.FIREWALL) is mapped to global ACSL annotations expressing the required isolation properties (either as ACSL predicates used as invariants, or as metaproperties for Frama-C/MetAcsl). A second mapping is then done for Security Functional Requirements (SFRs), describing the characteristics of the model (i.e. the contained functionality) for which the Security Objective is proved. Both mappings (Security Objectives and SFRs) are helpful to clearly articulate the security behavior chosen to be modeled, in other words, to clearly define the scope of the Security Policy Model (SPM). Security Functional Requirements are linked against some ACSL annotations in the SPM, it can be for example the contract of the firewall function or assignment specification of important variables decisive for firewall result (current context or context of objects). To make the mapping easier to evaluate, separate tables are provided to "translate" the Common Criteria terms of low granularity (e.g. Sharing, Currently Active Context) into their counterparts inside the Security Policy Model. It helps the evaluator to check the correctness of the mappings for particular SFRs based on these terms. For instance, SFRs introduced in II-C are mapped to their formalization in ACSL (depicted in Fig. 3 and4) as follows:

SFR1 As sensitive security attributes are stored in object headers, global invariants are used to prove that the headers of objects are intact during the entire VM run. The predicate object_headers_intact{Pre, Post} used in a function contract, ensures that the content of object headers is the same at the end of the function as it was before executing the function. Thus, preserving this predicate throughout all relevant function calls during the VM run ensures the integrity of sensitive security attributes. SFR2 Metaproperties persi_objects_integrity and persi_objects_confidentiality target integrity and confidentiality of object data respectively. In case the owner of a persistent object is different from the active context (JCC), any accessed memory location must be separated from the body of this object (see also Section IV-C).

VI. SPECIFICATION EFFORT AND PROOF RESULTS

Fig. 7. Specification effort for real-life code.

Applying deductive verification on large real-life industrial code requires a lot of care in order to avoid or at least mitigate scalability issues. In this project, we managed to ensure the scalability of our approach despite of the big size of the analyzed C code (7014 lines of C code split into 381 functions). As low-level operations are difficult to handle by automated provers, it was necessary to make abstraction of such low-level operations in a sound way. For that purpose, we introduced 162 lines of ghost code, carefully chosen to help automatic provers. As shown in Fig. 7, 35,480 lines of user-provided ACSL annotations were required in order to formally specify the security policy of the targeted virtual machine. However, we succeeded to reduce this effort up to 13,432 lines of ACSL annotations thanks to the usage of parameterized macros that gather redundant annotations. The effort is still considerable. 396,603 lines of ACSL annotations were automatically generated from 36 metaproperties (written in 350 lines of ACSL) only by Frama-C/MetAcsl. 2,290 lines of ACSL annotations were automatically generated by Frama-C/RTE.

Our proof scales reasonably well with an increasing number of proof goals. In particular, thanks to the translation of metaproperties into annotations that do not overload proof contexts, the metaproperty-based approach scales very well, despite a great number of generated annotations. Overall, it takes almost 3h30m to prove 52,198 proof goals. 99% of proof goals are automatically discharged by automatic provers. However, an important manual effort is required to maintain the proof of remaining proof goals when the code or the formal model are updated.

VII. RELATED WORK

A classical approach of applying formal verification on a JavaCard platform relies on a high-level formal model [START_REF] Chetali | About the world-first smart card certificate with EAL7 formal assurances[END_REF], [START_REF]Industrial use of formal methods for a high-level security evaluation[END_REF], [START_REF] Bolignano | Formal proof of a secure OS full trusted computing base[END_REF]. Several case studies have adopted this approach. An executable formal semantics of the JavaCard Virtual Machine (JCVM) and the ByteCode Verifier (BCV) is proposed in [START_REF] Barthe | A formal executable semantics of the JavaCard platform[END_REF] with 15,000 lines of specification in the Coq proof assistant. An operational semantics of a language modeling the JCVM behavior is proposed in [START_REF] Siveroni | Operational semantics of the Java Card Virtual Machine[END_REF], [START_REF] Éluard | An operational semantics of the Java Card firewall[END_REF]. Authors of [START_REF] Nguyen | Certifying native java card API by formal refinement[END_REF] describe a refinement-based approach, relying on the Coq proof assistant, to show that a native JavaCard API function fulfills its specification. In general, in such approaches, the traceability of formally proven properties may require an important effort to be justified because of the gap between the formal model and the source code [START_REF] Beckert | Mind the gap: Formal verification and the common criteria (discussion paper)[END_REF]. Among tools designed and/or used for the purpose of providing formal guarantees about JavaCard platform security properties (but not in a Common Criteria context) we can list: Key [START_REF] Mostowski | Fully verified Java Card API reference implementation[END_REF], KRAKATOA [START_REF] Marché | The KRAKATOA tool for certification of Java/JavaCard programs annotated in JML[END_REF] and Caduceus [START_REF] Andronick | Formal verification of security properties of smart card embedded source code[END_REF].

This work is also broadly related to other projects in which real-world software is verified. For instance, formal verification of the seL4 microkernel (comprising 8700 lines of C and 600 lines of assembly) was performed in a certification context [START_REF] Klein | seL4[END_REF]. Heitmeyer et al. [START_REF] Heitmeyer | Formal specification and verification of data separation in a separation kernel for an embedded system[END_REF] report on evidence for a Common Criteria evaluation of an embedded software system, which uses a separation kernel (of over 3,000 lines of C and assembly code). Although the separation kernel enforcing data separation was annotated with pre-and postconditions in the style of Hoare and Floyd, the machine-checked proof is not directly applied on the C code. A mechanized formal proof is performed on a Top Level Specification (TLS) of the separation-relevant behavior of the kernel. A correspondence between the annotated code and the TLS was established separately.

In general, while the usage of formal methods is, indeed, a costly validation technique, it is often seen as counterproductive to the current industrial development processes, even for having an advantage over competitors. Especially when the obtained certificates need to be renewed regularly, in order to try to cope with the dynamic landscape of a product's life-cycle [START_REF] Kaluvuri | A quantitative analysis of Common Criteria certification practice[END_REF]. We believe that our approach is well-suited to optimize, albeit still high, the investments to reach EAL7 certifications using formal methods. Further enhancements of our approach for automatic generation of Common Criteria documentation like in [START_REF] Heck | A software product certification model[END_REF], [START_REF] Tuong | Deeply integrating C11 code support into Isabelle/PIDE[END_REF], [START_REF] Bezzecchi | Making agile development processes fit for V-style certification procedures[END_REF], may be a step further for a better integration in current industrial development processes.

VIII. DISCUSSION AND LESSONS LEARNED

A. From the Developers' Point of View

An important drawback of top-down certification approaches, based on a high-level model (e.g. in Coq) is the traceability issue: the difficulty to relate the model with the real-life code. It can be complex to ensure that the model faithfully represents the behavior of the code. Another issue is the maintainability of the model for the developers: changes can be difficult to integrate and can require a very significant review of the whole proof.

a) Benefits: The presented bottom-up approach facilitates traceability since the SPM in that case is (a subset of) the real-life code, with the same structure (same functions, variables, data structures, etc.). Another benefit of the proposed approach is a better maintainability (in particular, in case of minor code updates or scope extensions) and a more straightforward extension from EAL6 to EAL7. Compared to a top-down approach, where significant model and proof changes are often required for more complex properties or a larger scope, integration of new properties or functions for an EAL7 certification in the presented approach can more significantly rely on the proof performed for the EAL6 level. Small design changes can easily be integrated in order to check if the security properties are affected.

The proposed bottom-up verification approach strongly benefits from automation, which is particularly important for a large industrial product. The link between the SPM and the real-life code in our project is explicit and can be automatically exploited by various tools. For example, they include syntactic code comparison and identification of possible differences-code transformations-between the SPM and the real-life code. It is important to efficiently identify and review such transformations (used in the verified code of the SPM e.g. to avoid some tool limitations or to realize some scope restrictions). Construction of a control-flow graph can help to identify the function hierarchy, in particular, to automatically distinguish fully verified functions (with a contract and a body), functions that are included as stubs (with a contract and a declaration but without a body) and excluded functions.

Furthermore, most security-related ACSL annotations in our approach are generated automatically from a few high-level security properties (stated as metaproperties [START_REF] Robles | Methodology for specification and verification of high-level properties with MetAcsl[END_REF]) by the Frama-C/MetAcsl plugin and then verified by Frama-C/WP. In this project, over 22,000 assertions are automatically generated by Frama-C/MetAcsl from (only!) 36 manually written metaproperties. Similarly, properties on the absence of runtime errors (RTE, also known as undefined behaviors) are generated automatically by another plugin, Frama-C/RTE, before being verified by Frama-C/WP.

Another major advantage of the approach is that it strongly relies on automatic proof. In the presented project, about 99% of proof goals are proved automatically by the Alt-Ergo solver or by the Frama-C/WP plugin (and its internal simplification engine Qed). A huge effort would be required to prove them interactively (that is, basically, manually) in a proof assistant. Finally, the proposed approach is suitable for a continuous integration process and it is planned to use it in the future in a continuous integration environment.

b) Challenges and Points of Attention: Those benefits also come with challenges for developers. An EAL6 certification with a bottom-up approach takes a more significant effort, already going closer to the implementation than actually required by the Common Criteria. Source-code level formal verification can be sensitive to tool scalability issues. Indeed, the tool has to deal simultaneously with high-level program properties and low-level properties (such as absence of runtime errors, presence of casts and bit-level operations), that can lead to a large number of relatively complex properties. A specific expertise and a good understanding of the capacities of the chosen verification tools and automatic solvers are required for the developers. Advanced tools like Frama-C/WP offer interactive proof features (proof scripts) that help the developer to finish most complex proofs.

A significant effort of manual annotation of the code is another challenge. In the presented project, ∼12,000 lines of annotations were written manually for ∼7,000 lines of C code (i.e. a factor of 1.7). Some of them rely on carefully chosen macros to avoid annotation redundancy, without using it the developers would have to write ∼35,000 lines of annotations (i.e. a factor of 5 compared to the C code). Another challenge is an efficient and meaningful organization of annotations and global properties-sometimes not obvious in modular verification-that can have an impact both on the capacity to prove and to define the mapping (see Section V).

B. From the Evaluators' Point of View a) Benefits: 1) the main benefit of the bottom-up approach for an evaluator is the immediate understanding of the formal entities, such as the modeling of the heap and the stack in a JCVM, because the program 2 has the same representation as the JCVM implementation, which has already been evaluated (ADV IMP); 2) as SFRs are directly represented in the program as formal annotations, the correspondence from the SFRs to the model can be easily understood by the evaluator, in particular to check that SFRs are modeled precisely enough to allow the verification of the security objectives; 3) only a single model needs to be reviewed because all the security and functional properties can be verified by the same model built on top of the implementation; 4) there is no refinement, no abstract level, no relation between multiple models to be evaluated; 5) the formal modeling of the subsystems is implicitly provided by the program; 6) apart from code transformations, the verified program has been written by a separate team (i.e. not the modeling team) which makes the relevance of the model easier to evaluate than when properties are verified on a dedicated model (in particular, a purely logical model); 7) the evaluator can directly check that well-known attack paths (e.g. type confusions) are not ignored by the model; 8) the bottom-up approach perfectly fits into the continuation of the evaluation process, unlike the top-down approach that involves a new (formal) design, which is more difficult to evaluate because of the traceability issue (as mentioned in Section VIII-A).

b) Limitations and Points of Attention: The limitations of the bottom-up approach are caused by the code complexity which is directly transferred to the model making the whole proof too complex to be fully reviewed by the evaluator. However, deductive verification ensures that properties are correctly propagated between functions in the callgraph. If a property is verified in the contract of a caller, deductive verification ensures that the contracts of the callees are complete enough to verify the property at the caller level. Therefore, the evaluator can have confidence in the proof, as long as the properties are correctly defined, and the program is correctly modeled in terms of code transformations, ghost code, and hypotheses. Hence, the main points of attention are the following.

• Properties: The bottom-up approach may make the properties more complex to evaluate because they directly rely on the implementation. Metaproperties, which are defined globally and whose number is limited, can be reviewed more easily but their statement remains as complex as the low-level annotations. • Code transformations: In some cases, the real-life code complexity cannot be fully supported by the existing tools, and requires manual transformations of the code. While some bugs in code transformations can be detected by the logical part of the model, some other bugs can lead to missing states and be hard to detect. An exhaustive manual review of the code remains difficult. Therefore, code transformations should rigorously follow a precise methodology even for the (apparently) simple cases. • Ghost code: Companion ghost code increases the complexity of the model, but also helps the evaluator to understand how the proof is conducted. The evaluator should detect non-companion ghost code that is used to create new concepts in the model either to express properties that cannot be directly expressed with the concepts of the implementation, or to simplify the model. The evaluator needs to check that these new concepts are valid and consistent with the implementation. • Hypotheses: Hypotheses can be introduced as preconditions of the entry point function, or contracts of stubs (whose code is not provided so their contracts are admitted). Stubs can include functions from the implementation removed for simplification, or new functions, specifically declared to introduce some local hypotheses at some points in the program. Local hypotheses make the understanding of the chain of reasoning and the detection of contradictory hypotheses more difficult. Additional tool-related hypotheses (like memory model assumptions in Frama-C/WP) also require specific attention of the evaluators.

C. Conclusion

The bottom-up approach brings many benefits to the certification process in terms of model understanding and confidence in verified properties. It helps to reduce the gap between the formally proved properties and the implementation, and should facilitate the step from EAL6 to EAL7 for the developers. Modern verification tools like Frama-C/WP and Frama-C/MetAcsl are capable to deal with real-life code after only a limited number of code modifications. The application of bottom-up approaches can require some adjustments in the existing terminology and certification guidelines (e.g. [START_REF] Sgdn | Note d'application : Modélisation formelle des politiques de sécurité d'une cible d'évaluation[END_REF]) released by certification bodies. Indeed, historically, they were designed with top-down approaches in mind, and their application on the bottom-up approach requires some clarifications. The existing evaluation methodology has to be extended with additional tasks for a careful analysis of properties, code transformations, ghost code and hypotheses.

Fig. 1 .

 1 Fig. 1. Top-down/bottom-up approaches.

  Fig. 2. A C function swap which permutes the values referred to by two given pointers p1 and p2, and its ACSL contract.

Fig. 3 .Fig. 4 .

 34 Fig.3. Example of a security-related predicate expressing that object headers are not modified between program points L1 and L2.

•Fig. 5 .

 5 Fig. 5. Representations of the TSFs

  Fig. 6. Structure of the SPM with respect to CC requirements.

See https://www.sogis.eu/