Michael Saint-Guillain

Christine Solnon

Yves Deville

The Static and Stochastic VRP with Time Windows and both random Customers and Reveal Times

Static and stochastic vehicle routing problems (SS-VRP) aim at modeling and solving real life problems by considering uncertainty on the data. In particular, customer data may not be known with certainty. Before the beginning of the day, probability distributions on customer data are used to compute a firststage solution that optimizes an expected cost. Customer data are revealed online, while the solution is executed, and a recourse strategy is applied on the first-stage solution to quickly adapt it. Existing SS-VRP variants usually make a strong assumption on the time at which a stochastic customer reveals its data (e.g., when a vehicle arrives at the corresponding location). We introduce a new SS-VRP where customer reveal times are stochastic. We define first-stage solutions and a recourse strategy for this new problem. A key point is to introduce waiting locations that are used in the first stage-solution to wait for the realization of customer stochastic data. We show how to compute the expected cost of a first-stage solution in pseudo polynomial time, in the particular case where the vehicles are not constrained by a maximal capacity. We also introduce a local search-based approach for optimizing the first-stage solution, and introduce a scale parameter to tune the precision and cost of the expected cost computation. Experimental results on small to large instances demonstrate its efficiency and flexibility.

Introduction

The Vehicle Routing Problem (VRP) aims at modeling and solving a real life common operational problem, in which a set of customers must be visited using a fleet of vehicles. Each customer comes with a certain demand. In the VRP with Time Windows (VRPTW), each customer must be visited within a given time window. A feasible solution of the VRPTW is a set of vehicle routes, such that every customer is visited exactly once during its time window and that sum of the demands along each route does not exceed the corresponding vehicle's capacity. The objective is then to find an optimal feasible solution, where optimality is usually defined in terms of travel distances.

The classical deterministic VRP(TW) assumes that customer data are known with certainty before the computation of the solution. Contrary to standard academic formulations, real world applications usually have missing part of the problem data when computing a solution. For instance, only a subset of the customer demands may be known before online execution. Missing demands hence arrive in a dynamic fashion, while vehicles are on their route. In such a context, a solution should contain operational decisions that deal with current known demands, but should also anticipate potential unknown demands. Albeit uncertainty may be considered for various attributes of the VRP (e.g., travel times), we focus on situations where the customer data are unknown a priori, and we assume that we have some probabilistic knowledge on missing data (e.g., probability distributions computed from historical data). This probabilistic knowledge is used to compute a first-stage solution which is adapted online when random variables are realized. Two different kinds of adaptations may be considered: Dynamic and Stochastic VRPTW (DS-VRPTW) and Static and Stochastic VRPTW (SS-VRPTW).

In the DS-VRPTW, the solution is re-optimized at each time-step, and this reoptimization involves solving an N P-hard problem so that it is usually approximated with meta-heuristics as proposed, for example, in [START_REF] Bent | Waiting and Relocation Strategies in Online Stochastic Vehicle Routing[END_REF][START_REF] Ichoua | Exploiting knowledge about future demands for real-time vehicle dispatching[END_REF][START_REF] Saint-Guillain | A Multistage Stochastic Programming Approach to the Dynamic and Stochastic VRPTW[END_REF]. Note that the DS-VRPTW assumes a probabilistic knowledge on the potential requests. In contrary, in [START_REF] Branke | Waiting Strategies for Dynamic Vehicle Routing[END_REF] for instance no prior knowledge is provided on the potential requests, which are then assumed to be uniformly distributed in the Euclidean plan.

In the SS-VRP(TW), no expensive reoptimization is allowed during online execution. When unknown information is revealed, the first stage solution is adapted online by applying a predefined recourse strategy whose time complexity is polynomial. In this case, the goal is to find a first stage solution that minimizes its total cost plus the expected extra cost caused by the recourse strategy. For example, in [START_REF] Bertsimas | A vehicle routing problem with stochastic demand[END_REF], the first stage solution is a set of vehicle tours which is computed offline with respect to probability distributions of customer demands. Real customer demands are revealed online, and two different recourse strategies are proposed: in the first one, each demand is assumed to be known when the vehicle arrives at the customer place, and if it is larger than or equal to the remaining capacity of the vehicle, then the first stage solution is adapted by adding a round trip to the depot to unload the vehicle; in the second recourse strategy, each demand is assumed to be known when leaving the previous customer and the recourse strategy is refined so that customers with null demands are skipped.

In this paper, we focus on the SS-VRPTW, and introduce a new variant where no strong assumption is made on the moment at which customer requests are revealed during the operations (contrary to most existing work that assume that customer requests are known either when arriving at the customer place, or when leaving the previous customer). In this new variant, called the SS-VRPTW with random Customers and Reveal time (SS-VRPTW-CR), the reveal times of customer requests are random variables. To handle uncertainty on reveal times, we introduce waiting locations when computing first-stage solutions: the routes computed offline visit waiting locations and a waiting time is associated with each waiting location. When a customer request is revealed, it is either accepted (if it is possible to serve it) or rejected. The recourse strategy then adapts routes so that all accepted requests are guaranteed to eventually be served. The goal is to compute the first-stage solution that minimizes the expected number of rejected requests.

Our motivating application is an on-demand health care service for elderly or disabled people. Health care services are provided directly at home by mobile medical units. Every person who's registered to the service can request a health care support at any moment of the day with the guarantee to be satisfied within a given time window. From historical data, we know, for each customer region and each time unit, the probability that a request appears. Given this stochastic knowledge, we compute a firststage solution. When a request appears (online), the recourse strategy is used to decide whether the request is accepted or rejected and to adapt medical unit routes. When a request is rejected, the system must rely on an external service provider in order to satisfy it. Therefore, the goal is to minimize the expected number of rejected requests.

Organization. In section 2, we review the existing studies on VRPs that imply stochastic customers. Section 3 formally defines the general SS-VRPTW-CR. Section 4 describes a recourse strategy for this problem. Section 5 shows how the expected number of rejected requests can be efficiently computed from a first stage solution and for a specific recourse strategy. Section 6 describes a local search-based approach for approximating an optimal first stage solution. Experimental results are analysed in section 7. Further research directions are finally discussed in section 8.

Related work

The most studied cases in SS-VRPs are stochastic customers (presence of customers are random variables), stochastic demands (quantities required by customers are random variables), and stochastic times (travel and/or service times are random variables). Since the SS-VRPTW-CR belongs to the first case, we focus this review on customers uncertainty.

The Traveling Salesman Problem (TSP) is a special case of the VRP with only one uncapacitated vehicle. [START_REF] Jaillet | Probabilistic traveling salesman problems[END_REF] introduced the TSP with stochastic Customers (SS-TSP-C), and provided mathematical formulations and a number of properties and bounds. In particular, he showed that an optimal solution for the deterministic problem can be arbitrarily bad in case of uncertainty. [START_REF] Laporte | A priori optimization of the probabilistic traveling salesman problem[END_REF] developed the first exact solution method for the SS-TSP-C, using the integer L-shaped method [START_REF] Laporte | The integer L-shaped method for stochastic integer programs with complete recourse[END_REF] to solve instances up to 50 customers. Heuristics for the SS-TSP-C have then been proposed (e.g. [START_REF] Jezequel | Probabilistic vehicle routing problems[END_REF][START_REF] Bertsimas | Computational approaches to stochastic vehicle routing problems[END_REF][START_REF] Bianchi | Extension of the 2-p-opt and 1-shift algorithms to the heterogeneous probabilistic traveling salesman problem[END_REF]) as well as meta-heuristics such as simulated annealing [START_REF] Bowler | Characterisation of the probabilistic travelling salesman problem[END_REF] or ant colony optimization [START_REF] Bianchi | An ant colony optimization approach to the probabilistic traveling salesman problem[END_REF].

The first SS-VRP with stochastic Customers (SS-VRP-C) has been studied by [START_REF] Jezequel | Probabilistic vehicle routing problems[END_REF] as a generalization of the SS-TSP-C. [START_REF] Waters | Vehicle-scheduling problems with uncertainty and omitted customers[END_REF] compared different heuristics. [START_REF] Bertsimas | A vehicle routing problem with stochastic demand[END_REF] considered a VRP with stochastic Customers and Demands (SS-VRP-CD). A customer demand is assumed to be revealed either when the vehicle leaves the previous customer or when it arrives at the customer's own location. Two different recourse strategies are proposed. For both strategies, closed-form mathematical expressions are provided to compute the expected total distance, provided a first stage solution. [START_REF] Gendreau | An exact algorithm for the vehicle routing problem with stochastic demands and customers[END_REF] and [START_REF] Séguin | Problemes stochastiques de tournees de vehicules[END_REF] developed the first exact algorithm for solving the SS-VRP-CD for instances up to 70 customers, by means of an integer L-shaped method. [START_REF] Gendreau | A Tabu Search Heuristic for the Vehicle Routing Problem with Stochastic Demands and Customers[END_REF] later proposed a tabu search to efficiently approximate the solution. Experimentations are reported on instances with up to 46 customers. [START_REF] Gounaris | An adaptive memory programming framework for the robust capacitated vehicle routing problem[END_REF] later developed an adaptive memory programming metaheuristic for the SS-VRP-C and assess it on benchmarks with up to 483 customers and 38 vehicles.

Particularly close to the SS-VRPTW-CR is the SS-TSP-C with Deadlines [START_REF] Campbell | Probabilistic traveling salesman problem with deadlines[END_REF]. Unlike the SS-VRPTW-CR, the set of customers is revealed at the beginning of the operations. A recent literature review on SS-TSP-C may be found in [START_REF] Henchiri | A probabilistic traveling salesman problem: a survey[END_REF].

[21] considered a variant of the SS-VRPTW-C, the Courier Delivery Problem with Uncertainty. Unlike the SS-VRPTW-CR, authors assume that customer presences are not revealed at some random moment during the operations, but all at once at the beginning of the day (that is, after computing the first stage solution).

Input data. We consider a complete directed graph G = (V, A) and a discrete time horizon H = [1, h], where interval [a, b] denotes the set of all integer values i such that a ≤ i ≤ b. To every arc (i, j) ∈ A is associated a travel time (or distance) d i,j ∈ N. The set of vertices V = {0} ∪ W ∪ C is composed of a depot 0, a set of m waiting locations W = [1, m] and a set of n customer regions C = [m + 1, m + n]. We note

W 0 = W ∪ {0} and C 0 = C ∪ {0}. The fleet is composed of K uncapacitated vehicles.
We consider the set R = C × H of potential customer requests such that an element (c, Γ ) ∈ R represents a potential request revealed at time Γ ∈ H for customer region c. To each potential request r = (c, Γ ) ∈ R is associated a deterministic demand q r ∈ [1, Q], a deterministic service duration s r ∈ H and deterministic time window [e r , l r ] with Γ ≤ e r ≤ l r ≤ h. We note p r the probability that r appears on vertex c at time Γ , and assume independence between request probabilities. Although our formalism imposes Γ ≥ 1 for all potential requests, in practice a request may be known with certainty that is, with probability 1.

To simplify notations, a request r = (c, Γ ) can be written in place of its own region c. For instance, the distance d v,c can also be intuitively written d v,r . Furthermore, we use Γ r to denote the reveal time of a request r ∈ R and c r for its customer region.

First stage solution. The first-stage solution is computed offline, before the beginning of the time horizon. It consists in a set of K vehicle routes visiting a subset of the m waiting vertices, together with time variables denoted τ indicating how long a vehicle should wait on each vertex. More specifically, we denote (x, τ ) a first stage solution to the SS-VRPTW-CR, where:

x = {x 1 , ..., x K } defines a set of K sequences of waiting vertices of W , such that each sequence x k starts and ends with 0 and each vertex of W occurs at most once in x. We note W x ⊆ W the set of waiting vertices visited in x. τ : W x → H associates a waiting time τ w ≥ 1 to every waiting vertex w ∈ W x ; Each sequence 0, w 1 , ..., w m , 0 in x is such that the vehicle is back to the depot before the end of the day.

In other words, x defines a Tour Orienteering Problem (TOP, see [START_REF] Chao | The team orienteering problem[END_REF]) to which each visited location is assigned a waiting time by τ . Given a first stage solution (x, τ ), we define on(w) = [on(w), on(w)] for each vertex w ∈ W x such that on(w) (resp. on(w)) is the arrival (resp. departure) time on w. In a sequence 0, w 1 , ..., w m , 0 in x, we then have on(w i ) = on(w i-1 )+d wi-1,wi and on(w i ) = on(w i )+τ wi for i ≥ 1 and assume both on(0) = 1 and on(0 Recourse strategy and second stage solution. A recourse strategy R states how a second stage solution is gradually constructed as requests are dynamically revealed. In this paragraph, we define the properties of a recourse strategy. An example of recourse strategy is given in Section 4.

) = on(w m ) + d w m ,0 ≤ h.
Let ξ ⊆ R be the set of requests that reveal to appear by the end of the horizon H. The set ξ is also called a scenario. We note ξ t ⊆ ξ the set of requests appearing at time t ∈ H, i.e., ξ t = {r ∈ ξ : Γ r = t}. We note ξ 1..t = ξ 1 ∪ ... ∪ ξ t the set of requests appeared up to time t.

A second stage solution is incrementally constructed at each time unit by following the skeleton provided by the first stage solution (x, τ ). At a given time t of the horizon, we note (x t , A t ) the current state of the second stage solution:

x t defines a set of vertex sequences describing the route operations performed up to time t. Unlike x, we define x t on a graph that also includes the customer regions. Operations described in x t must satisfy the time window and vehicle capacity constraints imposed by the VRPTW. A t ⊆ ξ 1..t is the set of accepted requests up to time t. Requests of ξ 1..t that do not belong to A t are said to be rejected.

We distinguish between requests that are accepted and those that are both accepted and satisfied. Up to a time t, an accepted request is said to be satisfied if it is visited in x t by a vehicle. Accepted requests that are no yet satisfied must be guaranteed to be eventually satisfied according to their time window.

Figure 1(right) illustrates an example of second stage solution being partially constructed at some moment of the time horizon. Before starting the operations (time 0), x 0 is a set of K sequences that only contain vertex 0, and A 0 = ∅. At each time unit t ∈ H, given a first stage solution (x, τ ), a previous state (x t-1 , A t-1 ) of the second stage solution and a set ξ t of requests appearing at time t, the new state (x t , A t ) is obtained by applying a specific recourse strategy R:

(x t , A t ) = R (x, τ ), (x t-1 , A t-1 ), ξ t . (1)
A necessary property of a recourse strategy is to avoid reoptimization. We consider that R avoids reoptimization if the computation of (x t , A t ) is achieved in polynomial time.

We note cost(R, x, τ, ξ) = |ξ \ A h | the final cost of a second stage solution with respect to a scenario ξ, given a first stage solution (x, τ ) and under a recourse strategy R. This cost is the number of requests that are rejected at the end h of the time horizon.

Optimal first stage solution. An optimal first stage solution (x, τ ) to the SS-VRPTW-CR minimizes the expected cost of the second stage solution under a given strategy R, satisfying statements (2)-( 3):

(SS-VRPTW-CR) Minimize x,τ Q R (x, τ ) (2)
s.t. (x, τ ) is a first stage solution. ( 3
)
The objective function Q R (x, τ ), which is nonlinear in general, determines the expected number of rejected requests, i.e. requests that fail to be visited under recourse strategy R and first stage solution (x, τ ):

Q R (x, τ ) = ξ⊆R Pr(ξ) cost(R, x, τ, ξ) (4)
where Pr(ξ) defines the probability of scenario ξ. Since we assume independence between requests, we have Pr(ξ)

= r∈ξ p r • r∈R\ξ (1 -p r ).

Description of a recourse strategy

In order to avoid reoptimization, the set R of potential requests is ordered. Furthermore, given a first stage solution (x, τ ) that visits the set W x of waiting locations, each potential request r = (c, Γ ) ∈ R is assigned to exactly one waiting vertex (and hence, a vehicle) in W x . Informally, the recourse strategy accepts a new request if it is possible for the vehicle associated to its corresponding waiting vertex location to adapt its first stage tour to visit the customer. The vehicle will then travel from the waiting location to the customer and return to the waiting location. Time window constraints should be respected, and the already accepted requests should not be perturbed. In the recourse strategy we propose here, we assume the vehicles not to be constrained by a maximal capacity.

Request ordering. Before computing first-stage solutions, we order R by increasing reveal time Γ r first, end of time window l r second and request number r to break ties. Let < R denote this total strict order on R. Whereas the remaining of the paper is based on the assumption of total order on Γ r , the ordering criteria may be modified without loss of generality (e.g., replacing l r by e r ), as long as the total order remains strict and primarily based on Γ r , i.e.

∀r 1 , r 2 ∈ R, Γ r1 < Γ r2 ⇒ r 1 < R r 2 .
Request assignment according to a first stage solution. Given a first-stage solution (x, τ ), we assign each request of R to a waiting vertex visited in (x, τ ). This assignment is computed for each first stage solution (x, τ ) before the application of the recourse strategy. As an optimally fair distribution of the potential requests might be excessively expensive to compute, we propose the following heuristic.

Let t min

r,w = max{on(w), Γ r , e r -d w,r } be the minimum time for leaving waiting location w to satisfy request r. Indeed, a vehicle cannot handle r before (1) the vehicle is on w, (2) r is revealed, and (3) the beginning e r of the time window minus the time d w,r needed to go from w to r.

Let t max r,w = min{l r -d w,r , on(w) -d w,r -s r -d r,w } be the latest time at which a vehicle can handle r (which also involves a service time s r ) from waiting location w and still leave it in time t ≤ on(w).

Given a first stage solution (x, τ ), we assign each request r ∈ R either to a waiting vertex of W x or to ⊥ (to denote that r is not assigned). We note w(r) this assignment which is computed as follows:

Let W x r = {w ∈ W x : t min r,w ≤ t max r,w } be the set of feasible waiting locations for r If W x r = ∅ then set w(r) to ⊥ (r is always rejected) Else set w(r) to the feasible vertex of W x r that has the least number of requests already assigned to it (break further ties w.r.t. vertex number) Once finished, the request assignment ends up with a partition {π ⊥ , π 1 , ..., π K } of R, where π k is the set of requests assigned to the waiting vertices visited by vehicle k and π ⊥ is the set of unassigned requests (such that w(r) = ⊥). We note π w k ⊆ π k the set of requests assigned to w ∈ W x in route k. We note fst(π w k ) the first request of π w k according to order < R , and for each request r ∈ π w k such that r = fst(π w k ) we note prv(r) the request of π w k that immediately precedes r according to order < R .

Using the recourse strategy to adapt a first stage solution at a current time t. At each time step t, the recourse strategy is applied to compute the second stage solution (x t , A t ), given the first stage solution (x, τ ), the second stage solution (x t-1 , A t-1 ) at the end of time t -1, and the incoming requests ξ t .

A t is the set of accepted requests. It is initialized with A t-1 . Then, each incoming request of ξ t is considered (taken by increasing order of < R ) and either accepted (added to A t ) or rejected (not added to A t ) by applying the following decision rule:

Let k be the vehicle associated with r (i.e., r ∈ π k ) Let y : R → H be the function returning the time at which k finishes to satisfy all accepted requests that precede r (according to < R ) and reaches waiting vertex w(r). Namely, y(r) is the time at which k is available for r and is defined by: If r = fst(π w k ), then y(r) = on(w) else if prv(r) / ∈ A t then y(r) = y(prv(r)) else y(r) = max(y(prv(r)) + d w,prv(r) , e prv(r) ) + s prv(r) + d prv(r),w If y(r) allows k to reach r during its time window and to arrive in time to its next waiting location (i.e., y(r) ≤ t max r,w(r) ) then r is accepted and added to A t ; otherwise it is rejected.

Once A t has been computed, vehicle operations for time unit t must be decided. Vehicles operate independently of each other. If vehicle k is traveling between a waiting location and a customer region, or if it is serving a request, then its operation remains unchanged; Otherwise, let w be the current waiting location (or the depot) of vehicle k:

If t = on(w), the operation for k is "travel from w to the next waiting vertex (or the depot), as defined in the first stage solution" Otherwise, let P = {r ∈ π w k |c r / ∈ x t ∧ (r ∈ A t ∨ t < Γ r )} be the set of requests of π w k that are not yet satisfied and that are either accepted or with unknown revelation If P = ∅, then the operation for k is "travel back to the depot" Otherwise, let r next be the smallest element of P according to < R → If t < t min r next ,w , then the operation for k is "wait until t + 1" → Otherwise, the operation is "travel to r next , serve it and come back to w" Figure 2 shows an example of second stage solution at a current time t = 17, from an operational point of view. A filled cross represents a request that appeared, an empty one a request that is either still unknown (e.g., r8) or revealed as being absent (that is didn't appear, e.g., r5). Here π k = ra, r1, . . . , r9 is the sequence of requests assigned to the vehicle, according to (x, τ ). We assume qr = sr = 0, ∀r ∈ R. sat(r) represents, for a request r, the time at which r gets satisfied.

Expected cost of second stage solutions

Provided a recourse strategy R and a first stage solution (x, τ ) to the SS-VRPTW-CR, a naive approach for computing Q R (x, τ ) would be to literally follow equation (4), therefore using the strategy described by R in order to confront (x, τ ) to each and every possible scenario ξ ⊆ R. Because there is an exponential number of scenarios with respect to |R|, this naive approach is not affordable in practice. In this section, we show how the expected number of rejected requests Q R (x, τ ) under the recourse strategy described in Section 4 may be computed in O(nh 2 ) using closed form expressions, in the special case where vehicles are of infinite capacity.

We assume that the potential request probabilities are independent from each other such that, for any couple of requests r, r ∈ R, the probability p r∩r that both requests appear is given by p r∩r = p r • p r .

Expected cost. Q R (x, τ ) is equal to the expected number of rejected requests, which in turn is equal to the expected number of requests that reveal to appear minus the expected number of accepted requests. The expected number of revealed requests is given by the sum of all request probabilities, whereas the expected number of accepted requests is equal to the sum, for every request r, of the probability that it belongs to A h :

Q R (x, τ ) = r∈R p r - r∈R Pr{r ∈ A h } = r∈R p r -Pr{r ∈ A h } (5)
where the right-hand side of the equation comes from the independence hypothesis.

If we consider a request r ∈ π k , the probability Pr{r ∈ A h } only depends on the time at which vehicle k is available for r, which itself depends on previous operations. Recall the y : R → H function described in section 4: y(r) is that time. Whereas y(r) is deterministic for a specific scenario, it is not anymore in the context of the computation of Pr{r ∈ A h } and we are thus interested in its probability distribution. More specifically, we compute the probability that, at a time t ∈ H, a request r already appeared and the vehicle leaves w(r) to satisfy it. Let's call this probability g 1 (r, t): g 1 (r, t) ≡ Pr{request r appeared at time t ≤ t and departureTime(r) = t} where departureTime(r) is the time at which the vehicle leaves vertex w(r) in order to serve r, if r has been accepted. According to the recourse strategy, for a specific scenario we see that departureTime(r) = max(y(r), t min r,w(r) ).

The probability Pr{r ∈ A h } that a request r gets satisfied is the probability that both r appears and that departureTime(r) ≤ t max r,w , that is:

Pr{r ∈ A h } = t max r,w t=1 g 1 (r, t) = t max r,w t=t min r,w g 1 (r, t). ( 6
)
The calculus of g 1 (r, t) is less obvious. Since departureTime(r) depends on previous operations on the same waiting location w = w(r), we calculate g 1 (r, t) recursively starting from the first request r 1 = fst(π w k ) assigned to the waiting location, up to the current request r. The second stage solution strictly respects the first stage schedule when visiting the waiting vertices, that is, these are guaranteed to be visited according to their arrival (on) and departure (on) times. The base case is then:

g 1 (r 1 , t) = p r1 , if t = max(on(w), t min r1,w ) 0 otherwise. ( 7
)
Indeed, if r 1 appeared then the vehicle leaves w at time t min r1,w , unless it has not yet reached w at that time. The general case of a request r > R r 1 , r ∈ π w k , depends on the time at which the vehicle gets rid of the preceding request prv(r). Let f (r, t) be the probability that, at time t, the vehicle either reaches back w after having served r, or discard r because it is not satisfiable or because it has revealed not to appear. It represents the time at which the vehicle becomes available for the next request after r in π w k , if any (computation of f is detailed below). We define g 1 (r, t) based on fprobabilities of the previous request prv(r):

g 1 (r, t) =        p r • f (prv(r), t) if t > t min r,w p r • t min r,w t =on(w) f (prv(r), t ) if t = t min r,w 0 otherwise. (8)
Indeed, if t > t min r,w the vehicle leaves w to serve r as soon as it gets rid of the previous one prv(r). In such case, g 1 (r, t) is the probability that both r has already appeared and the vehicle is available for it at time t, that is, finished with request prv(r) at time t. At any time below t min r,w , the probability that the vehicle leaves w must obviously be zero, since t min r,w is the minimum time for serving r from location w = w(r). At time t = t min r,w , we must consider the possibility that the vehicle was waiting for being able to serve r, but from an earlier time t < t min r . The overall probability that the vehicle leaves w for request r at time t = t min r,w is then p r times all the f -probabilities that the vehicle was actually available from a time on(w) ≤ t ≤ t min r,w . The f -probabilities of a request r depend on what exactly happened to r. Namely, from a time t there are two cases: either r consumed operational time, or it didn't at all:

f (r, t) = g 1 (r, t -S r ) • δ w (r, t -S r ) + g 1 (r, t) • 1 -δ w (r, t) + g 2 (r, t). ( 9
)
where S r = d w,r +s r +d r,w and the function δ w (r, t) returns 1 iff request r is satisfiable from time t and vertex w, i.e., δ w (r, t) = 1 if t ≤ t max r,w , and δ w (r, t) = 0 otherwise. The first term in the summation of the right hand side of equation ( 9) gives the probability that request r actually appeared and got satisfied. In such a case, departureTime(r) must be the current time t, minus delay S r needed for serving r.

The second and third terms of equation ( 9) add the probability that the vehicle was available time t, but that request r did not consume any operational time. There are only two possible reasons for that: either r actually appeared but was not satisfiable (second term), or r did not appear at all (third term), where g 2 (r, t) is the probability that r did not appear and is discarded at time t, and is computed as follows. For the base case of first potential request r 1 = fst(π w k ), we have:

g 2 (r 1 , t) = 1 -p r1 if t = max(on(w), Γ r1 ) 0 otherwise (10)
The general case for r ≥ r 1 , r ∈ π w k , is quite similar to the one of function g 1 . We just consider the probability 1 -p r that r doesn't reveal and replace t min r,w by Γ r :

g 2 (r, t) =      (1 -p r ) • f (prv(r), t) if t > max(on(w), Γ r ) (1 -p r ) • max(on(w),Γr) t =on(w) f (prv(r), t ) if t = max(on(w), Γ r ) 0 otherwise. (11)
A note on implementation. Since we are interested in computing Pr{r ∈ A h } for each request r separately, by following the definition of g 1 , we only require the f -probability associated to prv(r) to be already computed. This suggests a dynamic programming Algorithm 1: Local search to compute a first stage solution of SS-VRPTW-CR 1 Let (x, τ ) be an initial feasible first stage solution.

2 Initialize the neighborhood operator op to 1 3 while some stopping criterion is not met do 4 Select a solution (x , τ ) at random in Nop(x, τ ) 5 if some acceptance criterion is met on (x , τ ) then set (x, τ ) to (x , τ ) and op to 1 ; 6 else change the neighborhood operator op to op % nop + 1 ;

7 return the best first stage solution computed during the search approach. Computing all the f -probabilities can then be incrementally achieved while filling up a 2-dimensional matrix containing all the f -probabilities.

Computational complexity. Complexity of computing the expected cost is equivalent to the one of filling up a |π w k | × h matrix for each visited waiting location w ∈ W x , in order to store all the f (r, t) probabilities. By processing incrementally on each waiting location separately, each matrix cell can be computed in constant time using equation [START_REF] Jezequel | Probabilistic vehicle routing problems[END_REF]. In particular, once the probabilities in cells (prv(r), 1 • • • t) are known, the cell (r, t) such that r = fst(π w k ) can be computed in O(1) according to equations ( 8) - [START_REF] Bianchi | Extension of the 2-p-opt and 1-shift algorithms to the heterogeneous probabilistic traveling salesman problem[END_REF]. Given n customer regions and a time horizon of length h, we have at most

|R| = nh ≥ w∈W x |π w k | potential requests. It then requires at most O(|R|h) = O(nh 2 ) constant time operations to compute Q R (x, τ ).

Local Search for the SS-VRPTW-CR

Algorithm 1 describes a Simulated Annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] local search approach for approximating the optimal first stage solution (x, τ ), minimizing Q R (x, τ ). The computation of Q R (x, τ ) is performed according to equations of section 5 and is considered from now as a black box. Starting from an initial feasible first stage solution (x, τ ), Algorithm 1 iteratively modifies it by using a set of n op = 9 neighborhood operators. At each iteration, it randomly chooses a solution (x , τ ) in the current neighborhood (line 4), and either accepts it and resets the neighborhood operator op to the first one (line 5), or rejects it and changes the neighborhood operator op to the next one (line 6). At the end, the algorithm simply returns the best solution (x * , τ * ) encountered so far.

Initial solution and stopping criterion. The initial first stage solution is constructed by randomly adding each waiting vertex in a route k ∈ [1, K]. All waiting vertices are thus initially part of the solution. The stopping criterion depends on the computational time dedicated to the algorithm. Neighborhood operators. We consider 4 wellknown operators for the VRP: relocate, swap, inverted 2-opt, and cross-exchange (see [START_REF] Kindervater | Vehicle routing: handling edge exchanges[END_REF][START_REF] Taillard | A tabu search heuristic for the vehicle routing problem with soft time windows[END_REF] for detailed description). In addition, 5 new operators are dedicated to waiting vertices: 2 for either inserting or removing from W x a waiting vertex w picked at random, 2 for increasing or decreasing the waiting time τ w at random vertex w ∈ W x , and 1 that transfers a random amount of waiting time units from one waiting vertex to another.

Acceptance criterion. We use a Simulated Annealing acceptance criterion. Improving solutions are always accepted, while degrading solutions are accepted with a probability that depends on the degradation and on a temperature parameter, i.e., the probability of

accepting (x , τ ) is e -1-Q R (x,τ )/Q R (x ,τ ) T
. The temperature T is updated by a cooling factor 0 < α < 1 at each iteration of Algorithm 1: T ← α • T . During the search process, T gradually evolves from an initial temperature T init to nearly zero. A restart strategy is implemented by resetting the temperature to T ← T init each time T decreases below a fixed limit T min .

Experimentations

Test instances. We have randomly generated instances for the SS-VRPTW-CR. Each test instance is drawn in a square of [100, 100] distance units, and is characterized by:

The number |C| ∈ {30, 50, 80} of customer regions, randomly distributed in the square. Each customer c ∈ C region hosts n TS potential requests. The number |W | ∈ {20, 30, 50} of waiting vertices, randomly distributed in the square. The size h = 480 of the horizon (corresponding to the number of minutes in an 8 hour day). The number n TS = 24 of time slots. Time slots are introduced because it is not realistic to detail request probabilities for each time unit of the horizon (i.e., every minute). We set n TS to 24 so that the probability that a request appears at a customer region is specified for 20 minute time slots. The number K ∈ {1, 3, 5, 10} of available vehicles.

Travel times between vertices correspond to Euclidean distances, divided by a velocity parameter specified in the instance. Figure 3 shows an example of test instance. As a convention, the first time slot is associated to time unit 1 whereas time slot i is associated to time unit 1 + (i -1) • h/n TS . In our instance a.1, a potential request r associated to time slot 2 has a reveal time Γ r = 21 and no potential request is associated to time units [START_REF] Ichoua | Exploiting knowledge about future demands for real-time vehicle dispatching[END_REF][START_REF] Henchiri | A probabilistic traveling salesman problem: a survey[END_REF], [START_REF] Chao | The team orienteering problem[END_REF]40], etc. All the test instances are available at http://becool.info. ucl.ac.be/resources/benchmarks-ss-vrptw-cr. Potential requests. Each potential request r, associated with a customer region and a time slot, comes with a deterministic service time s r = 10. The time window [e r , l r ] is such that e r is chosen uniformly in [Γ r , min(Γ r + h nTS ], h), and l r is chosen uniformly in [max(e l , d 0,r ), max(e l + 10, d 0,r )].

Scale parameter. A scale parameter is introduced in order to optimize expectations on coarser data, and therefore to speed-up computations. When equal to 1, expectations are computed while considering the original horizon. When scale > 1, expectations are computed from a coarse version of the initial horizon, scaled down by the factor scale. If scale =10 for instance, then the horizon is scaled to h = 48. All the time data, such as travel and service times, but also time windows and reveal times, are then scaled as well (rounding up to nearest integer). When working on a scaled horizon (i.e. scale > 1), Algorithm 1 deals with an approximate but easier objective function Q R (x, τ ), in O(n( h scale ) 2 ), and a reduced search space due to a coarse time horizon. Algorithm 1 is then modified by simply adapting line 7 to return the best solution encountered so far, according to the initial horizon. Each time a new best solution is found during the search, its true expected cost is computed after scaling it up back to the original horizon, multiplying arrival, departure and waiting times by a factor scale. Experimental plan. All experiments are done under a cluster composed of 32 64bits AMD Opteron 1.4GHz cores. The code is developed in C++11 and compiled with LLDB using -O3 optimization flag. The Simulated Annealing parameters are set to T init = 5, T min = 10 -6 and α = 0.995. Results. Table 1 shows average experimental results over 10 runs on 12 instances: Instances a.x (resp. b.x, c.x and d.x) are such that |C| = 30 (resp. 30, 50 and 80), |W | = 5 (resp. 20, 30 and 50), and K = 1 (resp. 3, 5 and 10). Results are reported with scale ∈ {1, 2, 5, 10} and with a CPU time limit ∈ {30, 180, 1800} seconds. Provided a limited computational time of 30 seconds, using a scaled horizon leads to better results. This is easily explained by the limited number of local search iterations performed when scale=1. As the computational time increases to 3 minutes, working on the original horizon size tends to provide better results. This trend is confirmed by moving to 30 minutes. As the available computational time increases, the accuracy in the objective function eventually overtakes the computational efficiency provided by scaled horizons, especially for large instances such as c.x and d.x.

With their unique vehicle and because requests are uniformly distributed, instances a.x may suffer from an evaluation function being roughly uniform. Sending the vehicle at some location to wait there is, most of the time, more or less equivalent to another location. Consequently, local optima are numerous and little diversified, the more promising ones being hard to detect when using scale 1 for only 30 minutes. In the contrary, using more vehicles (e.g. instances b.x) leads to a less uniform evaluation function. For example, concentrating all the vehicles in the same region would surely leads to poor results. On instances a.x, the diversification brought by scaled horizons then still prevails after 30 minutes. Given larger computation times (5 hours), results on scales 5 and 10 do not show a significant improvement: Scales 1 and 2 in contrary tend to take promising benefits of a larger computation time.

Figure 4 shows how, for instance c.1, the real objective function (i.e. according to the original horizon) evolves in average (over 10 runs) during an execution of Algorithm 1. By reducing both the granularity of the search space and the complexity of the objective function, the parameter scale can therefore be used as a tradeoff between responsiveness and good quality solutions on the long term. Figure 4 also shows that the parameter scale can dynamically be reduced during the search. We introduced a new stochastic VRP, the SS-VRPTW-CR. Unlike existing SS-VRPs with random customers, we don't make any assumption on the moment at which a customer reveals its presence or absence. Instead, this is treated as a random variable as well. We proposed a recourse strategy for a special case of the SS-VRPTW-CR, when there is no maximal vehicle capacities. We showed how the exact expected cost can be computed in pseudo-polynomial time under this recourse strategy, and how to integrate in an efficient meta-heuristic method. Experiments are driven on generated test instances of various sizes. The average results show how a scale parameter, controlling the of the time horizon, can be used to tune the optimization process in the case of limited computational times.

Maximal vehicle capacity constraints. The recourse strategy and equations we give can be extended to take care of vehicle capacities. We are currently working on a generalized version of these equations.

Contribution to online optimization. Another potential application of the SS-VRPTW-CR goes to online optimization problems such as the DS-VRPTW. Because of the huge complexity of reoptimization, heuristic methods are often preferred, including the so called Sample Average Approximation (SAA, see [START_REF] Ahmed | The sample average approximation method for stochastic programs with integer recourse[END_REF]). SAA relies on Monte Carlo sampling, making decisions based on a subset of the scenarios. Thanks to recourse strategies, the SS-VRPTW-CR provides an upper bound on the expected cost of a first stage solution under optimal reoptimization. The SS-VRPTW-CR could therefore be used as a subroutine in order to heuristically solve the DS-VRTPW, whilst considering the whole set of scenarios instead of only a subset of sampled ones. In such a context, the scale parameter we introduce in the experiments can be of great contribution.

Figure 1 (

 1 left) illustrates an example of first stage solution on a basic SS-VRPTW-CR instance.

Fig. 1 .

 1 Fig. 1. On the left: first stage solution with K = 3 vehicles. The depot, waiting vertices and customer regions are represented by a square, circles and crosses respectively. Arrows represent vehicle routes and integers indicate waiting times at waiting locations. Values preceded by 'd' indicate travel times. Waiting vertices h, i and m are not part of the first stage solution. Here on(D) = 1, on(a) = 3, on(a) = 9, on(b) = 12, on(b) = 16, etc. On the right: partial second stage solution (plain arrows). Filled crosses are accepted requests. Some accepted requests, such as r1, have been satisfied (or the vehicle is currently traveling towards the location, e.g., r2), while some others are not yet satisfied (e.g., r3).

Fig. 2 .

 2 Fig.2. Example of second stage solution at time t = 17, under strategy R1. A filled cross represents a request that appeared, an empty one a request that is either still unknown (e.g., r8) or revealed as being absent (that is didn't appear, e.g., r5). Here π k = ra, r1, . . . , r9 is the sequence of requests assigned to the vehicle, according to (x, τ ). We assume qr = sr = 0, ∀r ∈ R. sat(r) represents, for a request r, the time at which r gets satisfied.

Fig. 3 .

 3 Fig. 3. Map representation of instance a.1. The depot (square) is located at the center. The instance counts 30 customer regions (stars) and 20 waiting vertices (circles).Although it is not visible here, the instance has a time horizon of 480 units and counts 24 time slots. If the operational day lasts 8 hours, a time unit represents a 1 minute in real time and each time slot lasts 20 minutes.

  , 2, 5, 10}: 15.3 15.1 14.9 16.5 16.2 15.8 16.1 17.3 14.8 14.6 14.0 16.1

Fig. 4 .

 4 Fig.[START_REF] Branke | Waiting Strategies for Dynamic Vehicle Routing[END_REF]. Average evolution of the best real objective value in Algorithm 1, during 3600 seconds on instance c.1. During the first second, objective values rapidly decrease when optimizing on scaled horizon. Thereafter, depending on the available computation time, some scale factors reveal to be more efficient than others. For less than 1 minute, scale=5 leads to better results. With at least 10 minutes, using the original horizon is definitely better.

Table 1 .

 1 Experimental results while varying horizon scale and computational time.

			30 seconds			3 minutes			30 minutes
	scale:	1	2	5	10	1	2	5	10	1	2	5	10
	a.1	19.7 18.2 15.0 16.2	16.7 16.4 14.9 16.2	16.1 15.3 14.9 16.5
	a.2	22.2 20.1 16.4 17.5	17.7 16.8 16.2 17.4	16.7 16.3 16.2 17.5
	a.3	20.3 20.0 14.4 16.1	16.1 15.9 14.0 16.0	15.6 15.2 14.0 16.2
	b.1	21.1 16.9 11.1 11.0	8.2	6.3	7.1	9.9	5.6	5.7	6.9	9.6
	b.2	22.1 17.5	9.6 11.4	5.6	7.7	6.8	9.9	5.1	7.4	6.4	9.3
	b.3	22.2 17.0 10.8 11.9	8.7	7.2	7.8 11.1	6.2	8.3	7.2 10.8
	c.1	42.1 38.7 23.6 25.9	15.3 13.9 12.2 19.3	8.3	9.3 11.0 16.7
	c.2	43.9 37.2 25.8 27.8	14.0 14.9 13.5 19.8	9.9 10.4 11.6 17.9
	c.3	42.4 39.2 24.3 24.8	17.5 14.9 11.7 17.5	15.2	8.8 10.3 15.8
	d.1	71.6 67.9 54.8 54.3	46.5 30.4 26.1 39.6	11.8 13.0 19.2 32.7
	d.2	72.2 67.9 52.8 56.2	40.9 34.7 25.7 40.1	12.6 19.0 20.3 31.4
	d.3	73.0 67.4 53.8 51.5	40.8 37.2 23.0 35.9	17.4 11.9 18.4 28.4

Acknowledgments Christine Solnon is supported by the LABEX IMU (ANR-10-LABX-0088) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).


Martin Rayrole
email: martin.rayrole@fr.thalesgroup.com

David Faura
email: david.faura@fr.thalesgroup.com

Marc Gatti
email: marc-j.gatti@fr.thalesgroup.com

REDUCING CERTIFICATION GRANULARITY TO INCREASE ADAPTABILITY OF AVIONICS SOFTWARE

A strong certification process is required to insure the safety of airplanes, and more specifically the robustness of avionics applications.

To implement this process, the development of avionics software must follow long and costly procedures. Most of these procedures have to be reexecuted each time the software is modified.

In this paper, we propose a framework to reduce the cost and time impact of a software modification. With this new approach, the piece of software likely to change is isolated from the rest of the application, so it can be certified independently. This helps the system integrator to adapt an avionics application to the specificities of the target airplane, without the need for a new certification of the application.

Introduction

The Integrated Modular Avionics (IMA) principles offers the possibility to host several avionics software applications in a single avionics module, by using partitioning properties specified in avionics standards such as ARINC 653 [START_REF]Design Guidance for Integrated Modular Avionics[END_REF]. These principles reduce the place, the weight and the consumption of avionics equipments, but they increase the interference between applications. This interference makes the certification more complex, because one has to prove that an application failure cannot be propagated to another one on the same module. This additional certification work is significant because, as presented in [START_REF] John | Certification Concerns with Integrated Modular Avionics (IMA) Projects[END_REF], the impact of IMA in the certification process concerns most of the functional domains of an avionics platform.

In this paper, we propose a framework to reduce this additional certification work. The proposed framework allows isolating a piece of software from the rest of the application. At run time, a mechanism guaranties a strong (space and time) partitioning between these two parts of the application. These properties offer the possibility to modify and re-certify the isolated piece of software, without recertifying the rest of the application.

A New Approach for Building Avionics Applications

With the "state of the art" technologies, the certification of software obliges to freeze the entire source code before delivering the software to a system integrator. No additional source code can be added to certified software without re-launching the certification process for this software. To deal with this constraint, two solutions are used to adapt a software behavior to the integration context: defining configuration data, and using incremental certification.

Configuration Data

Some configuration data can be defined by the avionics application supplier to select the software behavior among a set of pre-defined behaviors.

Configuration data are easy to use by the system integrator, but they are limited to tuning capabilities that have been foreseen in detail during the software development phase.

Incremental Certification

The incremental certification is an efficient solution to deal with this certification complexity of IMA. The incremental certification mechanisms allow to independently certify the avionics platform and each of the partitions hosted in this platform. Thanks to these mechanisms, one partition can be modified in order to be adapted to a particular context, without affecting the certification of other partition housed in the same calculator.

The properties of incremental certification can be used to isolate two sub-parts of an avionics application: the piece of software likely to change can be run in an additional partition. The application is then composed of two partitions that can communicate through inter-partition communication mechanisms.

The creation of an additional partition offers good capabilities to adapt the software behavior, but it is heavy to implement and not adapted to small pieces of code.

Proposed Approach

The proposed approach offers wide capabilities to adapt an avionics application, with a minimum impact on the certification work.

For that, a new framework is defined. This framework offers the capability to modify a software behavior by adding a piece of source code, which can be certified independently from the software. This source code is compiled with a certified macro-compiler, and dataloaded on the embedded avionics calculator (called the "target platform"). During the target platform startup procedure, the compiled file is checked to verify its compatibility with the platform configuration. At run time, the compiled file is read, interpreted and executed inside a container, which guaranties the strong memory and time segregation between the compiled file execution and the other platform activities.

In the next chapters, we will present the proposed framework. Then we will describe the two main tools that support the framework: a "macrocompiler" and a "container". Finally, we will show how the framework can be used in a use case example.

Framework Description

Framework actors

Four actors are involved in the proposed framework:

• The "Module Supplier" provides the target platform. It also provides the couple of tools (a "macro-compiler" and a "container") that supports the framework. • The "Function Supplier" provides the avionics application. This application must be designed to be integrated within the framework.

• The "Adaptation Supplier" provides the piece of code that implements an adaptation of the application. This is a new actor with regard to the existing avionics integration processes. • The "System Integrator" provides the platform configuration. This configuration defines the resources allocated to the avionics application and to the container. The figure 1 gives an overview of the proposed framework, with the components provided by each actor.

Figure 1. Framework Overview and Involved Actors

Depending on the context, the "Adaptation Supplier" activities are supported by one of the three other actors.

The Adaptation Supplier activities are typically supported by the System Integrator, to adapt the application behavior to the specificities of the airplane architecture. However, the Adaptation Supplier activities may also be supported by the Function Supplier to adapt a generic application to its clients' needs. It may also be supported by the Module Supplier to adapt a generic platform to its clients' needs or to hardware specificities.

Framework Usage Scenario

In this paragraph, the numbers in parentheses refers to the numbers in figure 1.

The Function Supplier first provides an application (1) that realizes an avionics function. This application has been designed to be integrated in the framework, which means that a part of its functionalities is supposed to be developed by an Adaptation Supplier.

The Adaptation Supplier writes its own piece of source code (3) that implements an adaptation of the application behavior. He then compiles his source code by using the macro-compiler (6) provided by the Module Supplier. A set of "Elementary performance data" files ( 5) is also provided by the Module Supplier and used as an input to the macrocompiler.

The output of the macro-compilation is a "compiled file" (4) that is sent to the System Integrator.

The System Integrator uploads the compiled file on the target platform, as well as the application code, the platform configuration data (2), and the container (7). The container is a piece of software provided by the Module Supplier that can execute a compiled file without interfering with other platform activities.

At runtime, the application reads configuration data to know the amount of resources allocated to it. Then, each time it is needed, the application requests the container to execute the compiled file.

Supporting Tools Certification

The macro-compiler and the container must be developed according to DO-178 [START_REF]Software Considerations in Airborne Systems and Equipment Certification[END_REF] standard, and certified at least to the same Design Assurance Level (DAL) as the target platform.

Macro-Compiler Functionalities

The macro-compiler generates a compiled file from an adaptation source code and elementary performance data files, as shown in figure 2.

Figure 2. Macro-compiler Inputs and Output

The compiled file contains, among other information, a "macro-code" that is the result of the compilation of the adaptation source code. The instructions composing the macro-code are taken from a set of instructions independent from the target platform.

Adaptation Source Code (Macro-Compiler Input)

The adaptation source code respects a programming language syntax defined by the Module Supplier. This programming language is independent from the target platform and from the application.

This programming language has some specific properties that are detailed in the following paragraphs.

WCET calculation

The programming language must allow an automatic calculation of the WCET of the macrocode. Several techniques exist for calculating a WCET. A classification of these techniques has been proposed in [START_REF] Raimund | Classification of WCET Analysis Techniques[END_REF].

The capability to automatically calculate the WCET may have some impacts on the programming language syntax. For example, the maximum number of loop iterations may be explicitly specified in the source code. The macro-compiler could also implement a method, like in [START_REF] Andreas | Loop Bound Analysis based on a Combination of Program Slicing, Abstract Interpretation, and Invariant Analysis[END_REF], to bound the maximum number of loop iterations without an explicit declaration.

Variables Visibility

The programming language must distinguish between external and local variables.

External variables are visible inside the adaptation code, and inside the calling avionics application. External variables are typically used as a communication mean between the adaptation code and the calling application. They can be used by the calling application to give input data to the adaptation code, and by the adaptation code to return data to the calling application.

Local variables are visible only inside the adaptation code.

Memory Access Control

The programming language must allow a strict control of the memory access. This control guaranties that the adaptation code will never access memory outside the defined variables.

A simple solution to offer this guaranty is to forbid the usage of pointers, and to add the detection of out of bound array access.

Elementary Performance Data (Macro-Compiler Input)

The macro-compiler needs some performance data to calculate the WCET of the macro-code. These data are stored in "Elementary Performance Data" files that are given as input to the macro-compiler.

Each

file contains performance data corresponding to one platform type. A "platform type" represents all target platforms on which the execution of a compiled file takes exactly the same time. In the proposed framework, a "platform type" is characterized by

• The type and the version of the platform hardware.

• The type and the version of the Operating

System. • The version of the container. These files contain, for each platform type, the WCET of each instruction of the macro-code instruction set. Each WCET value represents the worst case time for reading, interpreting and executing the instruction on the platform type.

Compiled File (Macro-Compiler Output)

The compiled file generated by the macrocompiler is usable for a pre-defined set of platform types. This set of platform types is defined by the set of elementary performance data files given as input to the macro-compiler.

The generated compiled file contains • A header that gives some characteristics of the compiled file, including the memory size needed to execute the macro-code.

The exact content of this header is described below. • The macro-code, which is the translation of the adaptation source code into a list of instructions. These instructions are taken from a set of elementary instructions independent from the platform type. • The list of platform types able to execute the macro-code, and the macro-code WCET for each platform type. This WCET represents the worst case time for interpreting and executing the macro-code inside the container. This WCET also includes the time needed to analyze the request sent by the calling application.

Compiled File Header

The compiled file header contains all the information required by the container to parse the file, and all the information needed to guaranty the strong time and memory partitioning of the macrocode execution.

For that, the header contains at least the following information:

• The type and version of the macrocompiler, which determines the exact format of the compiled file.

Container Functionalities

The container is in charge of reading the compiled file, and of interpreting and executing the macro-code, when requested by the calling application. These operations are done with two guaranties:

• The execution time of the macro-code execution cannot exceed the WCET stored in the compiled file. • The accesses done by the container are limited to the container's data and external variables.

Container Memory Mapping

At runtime, the container is an embedded piece of software, running in the same partition as the avionics application that will request the execution of the macro-code.

The container needs to access the following memory spaces:

• The macro-code to execute,

• The local variables defined in the macrocode, • The external variables defined in the macro-code, • The container's code that receives and executes requests from the calling application, • Some additional data internally used by the container's code All these memory spaces are contained inside the memory allocated to the container, except the external variables. These external variables are used as a communication mean between the calling application and the container and, for that purpose, are shared between these two pieces of software.

In the proposed framework, the external variables are stored in the memory allocated to the calling application. So the container must have access to this memory, but the container behavior guaranties that it will never access to the calling application's data outside of these external variables. This guaranty is provided by the programming language properties and the macro-compiler (see "Memory Access Control" paragraph above) and by the container itself during the application start-up procedure (see "Initialization Phase" paragraph below).

Figure 3. Container Memory Mapping

Sequence Diagrams

The calling application sends a request to the container during the application start-up phase, to initialize the macro-code execution. Then, the calling application sends a request to the container, each time it is needed, to execute the macro-code.

These two phases are presented in the next paragraphs.

Initialization Phase

The figure 4 shows the sequence diagram of the initialization phase.

Figure 4. Initialization Phase Sequence Diagram

During its start-up procedure, the calling application requests a macro-code initialization to the container. The main objective of this initialization is to detect any configuration error that could prevent the correct execution of the macro-code. The initialization phase is also used to make the macrocode ready to be executed.

From a functional point of view, the calling application could request the initialization of the macro-code just before its first execution. However, for safety considerations, it is most recommended to initialize the macro-code during the application startup procedure. This recommendation avoids detecting a macro-code incompatibility during operational mode, where failures are much more critical.

The macro-code initialization request contains the following parameters:

• The compiled file address in non-volatile memory.

• The starting and ending memory addresses of the external variables.

• The allocated time for the macro-code execution. All these parameters are calculated by the calling application from configuration data, which include the amount of resources allocated to the application's partition.

When receiving this request, the container executes the following tasks:

• Header parsing.

The header of the compiled file is read from the non-volatile memory, and its content is parsed. • Memory check.

The memory sizes for macro-code and local variables (in the compiled file header) are compatible with the memory allocated to the container. And the memory size for external variables (also in the compiled file header) is exactly equal to the memory size received in the request.

• Macro-code load.

The macro-code is copied from the nonvolatile memory to the Random Access Memory (RAM).

• Initialization WCET check.

The number of platform types contained in the compiled file header is checked to be less or equal to the maximum number of platform types configured for the container. This verification is needed to bound the execution time of the initialization phase.

• Compatibility check.

The platform on which the container is currently running appears in the list of platform types stored in the compiled file.

• Execution WCET check.

For the current platform type, the WCET stored in the compiled file is less or equal to the allocated time received in the macro-code initialization request.

• Macro-code context creation.

A macro-code context is created to memorize all data needed to launch the macro-code execution: the memory addresses for the macro-code, for the external variables, and for the local variables. Finally, the container sends a response to the calling application with the initialization status and, if it is successful, with a macro-code context identifier.

Execution Phase

The figure 5 shows the sequence diagram of the execution phase.

Figure 5. Execution Phase Sequence Diagram

Once the macro-code has been initialized, the calling application can request its execution. Only the macro-code context identifier needs to be sent in the execution request.

When receiving an execution request, the container interprets and executes the macro-code. At the end of the execution, a response message is returned to the calling application, with the macrocode execution status. The calling application can then access the external variables to read the result values of the macro-code execution.

Multi Macro-Code Support

The presented framework is compatible with the use of several macro-codes in one calling application.

In such a case, a compiled file is generated for each adaptation source code, and all these compiled files are dataloaded on the target platform.

During its start-up procedure, the calling application requests the initialization of each macrocode. The container creates a macro-code context for each.

Once in operational state, the calling application can request the execution of any of the macro-codes, independently one from each other.

It should be noticed that the Adaptation Supplier activities, for each source code, can be supported by different actors. For example, an avionics application can use two macro-codes: a macro-code to adapt the behavior of the application to the client's need, and another macro-code to adapt the behavior of the application to the airplane. In this example, the first macro-code may be provided by the Module Supplier, and the second one by the System Integrator.

Use Case Example

In this example, an application named "Rack Manager" is used to manage an avionics rack (i.e. a set of avionics calculators integrated in a common structure). In case of power supply partial failure, the Rack Manager must stop some calculators installed in the rack. A decision rule is used to choose which calculators have to be stopped.

The Function Supplier provides the Rack Manager application, and he wants to offer to the System Integrator the possibility to define his own decision rule. So the System Integrator writes an adaptation source code to implement his decision rule, and he uses the macro-compiler to generate a compiled file.

The following figure gives an example of such an adaptation source code, using a C-like syntax.

Conclusion

The incremental certification has reduced the granularity of the certification from the module to the partition. In this paper, we have proposed a framework to reduce the certification granularity to a sub-part of software.

With this framework, an avionics application can be adapted to a client or airplane specificities, without re-certifying the entire application. This reduces the cost and the reactivity of software evolutions during the system integration phase.

Figure 6 .

 6 Figure 6. Adaptation Source Code Example

Figure 7 .

 7 Figure 7. Use Case Sequence Diagram Example

  The number of platform types considered by the macro-compiler for this adaptation code. This number matches the number of platform type and WCET stored at the end of the compiled file.

	external and local variables are stored in
	two different memory spaces.
	•
	• The macro-code length.
	• The memory size needed for external
	variables (with no padding bytes).
	• The memory size needed for local
	variables (with no padding bytes). The
	distinction between external and local
	variables is important because, as
	explained later (see "Erreur ! Source du
	renvoi introuvable." paragraph below),


Anni Määttänen
email: anni.maattanen@latmos.ipsl.fr

Franck Lefèvre

Sabrina Guilbon

Constantino Listowski

LATMOS/ Franck Montmessin

Vertical distribution of ozone at the terminator on Mars

Introduction

The SPICAM/Mars Express UV solar occultation dataset gives access to the ozone vertical distribution via the ozone absorption in the Hartley band (220-280 nm). The first results on the ozone vertical distribution from SPICAM were published by [START_REF] Lebonnois | Vertical distribution of ozone on Mars as measured by SPICAM/Mars Express using stellar occultation[END_REF]] using the stellar occultation data. We will present observations of ozone profiles derived from SPICAM solar occultations in the ultraviolet [START_REF] Määttänen | A complete climatology of the aerosol vertical distribution on Mars from MEx/SPICAM UV solar occultations[END_REF]]. The occultation datasets allow the retrieval of vertical profiles of the abundance of gaseous species, aerosol optical depth and the Ångström coefficient, which can be, through some hypotheses, related to particle size [START_REF] Montmessin | Stellar occultations at UV wavelengths by the SPICAM instrument: Retrieval and analysis of Martian haze profiles[END_REF]]. A data analysis code including a correction for mechanical oscillations induced by the MARSIS radar antennae, has been developed [START_REF] Määttänen | A complete climatology of the aerosol vertical distribution on Mars from MEx/SPICAM UV solar occultations[END_REF]] and was used in the analysis. All solar occultations of adequate quality between orbit numbers 0-10000 (Mars Years 27-30) have been analysed. The dataset spans 4 MY with a good seasonal and spatial coverage. We will present a global overview of the results on the vertical ozone distribution and describe a new method we have developed for comparing the observations to the model.

Dataset

SPICAM (SPectroscopie pour l'Investigation des Caractéristiques Atmosphériques de Mars) has been observing Mars in the ultraviolet (118-320 nm) and in the near infrared (1-1.7 µm) ranges since 2004. Details on the instrument and the observation modes can be found in [START_REF] Bertaux | SPICAM on Mars Express: Observing modes and overview of UV spectrometer data and scientific results[END_REF]]. In the occultation mode, the UV channel gives access to the vertical distribution of CO2, O3, and aerosols. Here we focus on the UV solar occultation dataset of about 650 profiles. The occultation technique is self-calibrated, since the spectra are normalised with the observed solar spectrum to acquire atmospheric transmissions. The transmission spectra are fitted with the Beer-Lambert law, which takes into account extinction by gaseous species (CO2 and O3 in the UV) and aerosols. Aerosol extinction is modeled with the so-called alpha-model, providing access to the Ångström coefficient, which depends on the size of the aerosols. For results on the aerosol vertical distribution see [START_REF] Määttänen | A complete climatology of the aerosol vertical distribution on Mars from MEx/SPICAM UV solar occultations[END_REF]].

Seasonal behavior of the ozone profiles

Figure 1 displays the climatology of the O3 profiles derived from the SPICAM solar occultations and the ozone profiles of the Mars General Circulation Model of the Laboratoire de météorologie dynamique (LMD MGCM). Here the comparison is made between local profiles inverted from observations using the onion peeling method and profiles extracted from the model at the grid point nearest to the terminator. The largest ozone amounts are found in spring (Ls = 20-50 • ) at high northern latitudes (60 • N-90 • N). This result is consistent with the observations of the O3 column performed by SPICAM in nadir mode [START_REF] Perrier | Global distribution of total ozone on Mars from SPICAM/MEX UV measurements[END_REF]] and with the LMD MGCM simulations that predict maximum O3 amounts in the dry polar vortices in winter and spring [START_REF] Lefèvre | Heterogeneous chemistry in the atmosphere of Mars[END_REF]]. However, our results shed light on the vertical distribution of ozone during that period: we find maximum O3 densities at about 20 km (note that the high aerosol opacity prevents observations below 10 km) and the top of the O3 springtime polar layer is located at about 30 km. These altitudes are significantly higher than those simulated by the LMD MGCM, which predicts an O3 polar layer that is essentially confined below 20 km (right column). This is also the case at southern polar latitudes (60 • S-90 • S), where O3 is detected by SPICAM solar occultations in autumn and spring above 20 km. This result is consistent with the O3 profiles derived by stellar occultation at the same time [START_REF] Lebonnois | Vertical distribution of ozone on Mars as measured by SPICAM/Mars Express using stellar occultation[END_REF]]. Both techniques show that more ozone is present above 20 km in the polar vortices than calculated by the LMD MGCM, which suggests that the dehydration in high latitudes extends higher in the real atmosphere than in the model.

At low southern latitudes (0-30 • S), the SPICAM solar occultations show between 30 and 40 km a clear seasonal maximum centred on aphelion (Ls = 70 • ). This result corroborates the nighttime SPICAM measurements by stellar occultation [START_REF] Lebonnois | Vertical distribution of ozone on Mars as measured by SPICAM/Mars Express using stellar occultation[END_REF]] and is also well reproduced by the LMD MGCM (right column). The buildup of this seasonal layer is mainly caused by changes in the vertical distribution of H2O [START_REF] Clancy | Annual (perihelion-aphelion) cycles in the photochemical behavior of the global Mars atmosphere[END_REF]]. Around aphelion, the production of hydrogen radicals is much reduced above the low hygropause and so is the ozone loss. In contrast, no ozone is detected by SPICAM at low latitudes in the perihelion season. This can be explained by the dramatic increase in ozone-destroying hydrogen radicals that follow the rise of the hygropause level at that season. The y-axis is altitude between -10 and 70 km above the areoid, the x-axis is the solar longitude, and the color scale goes from 0 to 4 x 10 9 molecules/m3.

How to compare modeled and observed profiles of a photochemically reactive species at the terminator?

Due to the rapid variation of ozone due to photolysis at sunrise and sunset, a classical comparison of local density profiles is not appropriate for solar occultations that are acquired at the terminator. The spherical symmetry hypothesis made in the onion-peeling vertical inversion method is not valid for photochemically active species (e.g., ozone) around terminator (see Figure 2). Here, we are using a method commonly used in the Earth community. This method exploits the slant profiles (ozone concentration integrated along the line-of-sight of the instrument) for the comparison. For each occultation, we have modeled with the 1D version of the LMD MGCM the ozone vertical and horizontal distribution with high solar zenith angle (or local time) resolution around the terminator (Figure 2). It can be seen that around the terminator the ozone distribution is quite far from spherical symmetry. We then integrate the model results following the lines-of-sight of the occultation (calculated from the observation geometry) to construct the modeled slant profiles (Figure 3).

So far, we have studied six occultations in this manner. In four out of six studied occultations, the slant profile comparison shows quite a good agreement between the model and the observations (for an example, see Figure 3). However, clear differences are seen in some cases, where the observed ozone densities are clearly smaller than the ones extracted from the model (not shown).

We have also tested the possible improvement of the vertical profile comparison by performing the onionpeeling vertical inversion of these slant profiles. The comparison of the model with SPICAM is slightly improved in some cases, but naturally the differences seen in the occultations where SPICAM observes clearly less ozone than predicted by the model (not shown) persist. In any case, the vertical inversion method we use even for inverting the LOS-integrated model profiles includes the spherical symmetry hypothesis, and thus the comparison of the local profiles is biased.

Thus, the most pertinent comparison is the one between the slant profiles. To completely avoid the problem of the spherical symmetry hypothesis, the inhomogeneity of the atmosphere (for temperature, density/pressure and photochemically active trace gas concentrations) should be fully accounted for within the vertical inversion method. This will be achieved within the UPWARDS project (Work Package lead: BIRA, Belgium).

Conclusion

We will present a four-year climatology of the ozone vertical distribution on Mars acquired from the analysis of SPICAM UV solar occultations, and we compare it to the modeled climatology. Due to strong spatial variation of ozone across the solar terminator, a classical comparison of local density profiles is not appropriate for solar occultations. We introduce here a method of extracting and comparing the slant profiles also from model results (often used in the Earth community). When the slant profiles are compared, some occultations still present large differences when compared to the model: in two out of six studied cases SPICAM sees a lot less ozone than predicted by the model. For a fully realistic observation-model comparison for photochemical species one should compare the slant profiles or account for the heterogeneity of the atmosphere around the terminator in the vertical inversion method (dedicated Work Package in the European project UPWARDS H2020).

Figure 1 :

 1 Figure 1: Seasonal evolution of ozone in different latitude bands as observed by SPICAM and modeled by the LMD Mars GCM.The y-axis is altitude between -10 and 70 km above the areoid, the x-axis is the solar longitude, and the color scale goes from 0 to 4 x 10 9 molecules/m3.

Figure 2 :

 2 Figure 2: Inhomogeneity of ozone around the terminator and its consequences on the analysis: Modeled ozone profiles (colors) around the terminator during the occultation as a function of the solar zenith angle and examples of three LOS of SPI-CAM (symbols).

Figure 3 :

 3 Figure 3: Ozone concentrations integrated along the line-ofsight (slant profiles) from the model (dash-dot line) and the observations (solid line).

Figure 4 :

 4 Figure 4: Ozone vertical profiles from SPICAM (solid line) and the model (dash-dot line) after vertical inversion of the slant profiles from Figure 3, and comparison to the modeled local profile at the terminator (dashed line).

Acknowledgements

SPICAM is funded by the French Space Agency CNES and this work has received funding from the European Union's Horizon 2020 Programme (H2020-Compet-08-2014) under grant agreement UPWARDS-633127.