Pierre Martinon

Joseph Gergaud
email: gergaud@enseeiht.fr

Using switching detection and variational equations for the shooting method

Keywords: Shooting method, Jacobian evaluation, switching detection, variational equations

We study in this paper the resolution by single shooting of an optimal control problem with a bang-bang control involving a large number of commutations. We focus on the handling of these commutations regarding the precise computation of the shooting function and its Jacobian. We first observe the impact of a switching detection algorithm on the shooting method results. Then, we study the computation of the Jacobian of the shooting function, by comparing classical finite differences to a formulation using the variational equations. We consider as an application a low thrust orbital transfer with payload maximization. This kind of problem presents a discontinuous optimal control, and involves up to 1800 commutations for the lowest thrust.

Context

The context of this work is the numerical resolution, by a single shooting method, of an optimal control problem with a discontinuous (bang-bang) optimal control. This large number of commutations causes two main problems. First, how to determine the correct control structure, i.e. the number and location of commutations. This matter, which is not the object of the present paper, was solved in our case by the use of continuation methods, as recalled briefly 1 in Section 1.2.

Then, these discontinuities cause numerical difficulties with the shooting function, and require careful handling if one wants to solve the problem precisely. We focus here on two points, namely the computation of the shooting function itself, and of its Jacobian required for the shooting method. The first point involves the detection of commutations during the integration of the initial value problem needed to compute the shooting function. The second point deals with the use of variational equations to compute the Jacobian, compared to traditional finite differences.

These points are illustrated with numerical experiments on a low thrust orbital transfer with minimization of the fuel consumption, which can present up to 1800 commutations.

Shooting method

We begin with a brief presentation of the single shooting method, which is part of the indirect methods, and is based on Pontryagin's Maximum Principle (we refer readers interested in these methods to [START_REF] Pontriaguine | Théorie Mathématique des Processus Optimaux[END_REF][START_REF] Cesari | Optimization theory and application[END_REF][START_REF] Ascher | Numerical Solution of Boundary Value Problems for Ordinary Differential Equations[END_REF][START_REF] Stoer | Introduction to Numerical Analysis[END_REF][START_REF] Bonnans | The shooting algorithm for optimal control problems: a review of some theoretical and numerical aspects[END_REF] for instance). We recall that direct methods, on the other hand, typically involve the partial or total discretization of the problem, and then use various approaches (SQP and interior point techniques for instance) to solve the resulting optimization problem. Direct methods are thus supposed to be robust, but the counterpart is a relatively low precision, and a huge problem size depending on the discretization stepsize used. This makes these methods ill-suited to some particular cases, such as the problems studied here, which present a bang-bang control structure with a huge number of commutations.

Back to the indirect methods, single shooting consists in finding a zero of the shooting function associated with the original problem. There is no discretization, even if the method still involves an integration of the system in some way. It is a fast and high precision method, but often requires a good initial guess: as they typically consist in applying a Newton solver to the shooting function, the convergence radius may be quite small, depending on the problem. This is particularly true for problems involving commutations, or worse, state constraints or singular arcs (both are not in the scope of this paper), as it is often necessary to have a good knowledge of the structure of the optimal control. In our case, we use a continuation approach to obtain a suitable initial point, as described later.

We consider a general optimal control problem in the Bolza form (P )

           M in g(t 0 , x(t 0 ), t f , x(t f )) + t f t0 l(t, x, u) dt Objective ẋ = f (t, x, u) Dynamics u ∈ U Admissible Controls ψ 0 (t 0 , x(t 0 )) = 0 Initial Conditions ψ 1 (t f , x(t f )) = 0
T erminal Conditions Notation/Remark: for clarity, the time t will often be omitted in the formulas, except in the ambiguous cases.

We use here and in all the following the notations: x ∈ R n for the state, u ∈ R m for the control, U the compact convex set of admissible controls, and f for the state dynamics. We assume that U is a compact convex subset of R m and that the initial and final times t 0 and t f are fixed. We introduce the costate p, of same dimension as the state x, and define the Hamiltonian by H(t, x, p, u) = l(t, x, u) + (p|f (t, x, u)). Remark: we thus assume in all the following that we are in the normal case, meaning that the costate p 0 associated to the integral objective l is non zero, and can be made equal to 1, by dividing the whole costate p by p 0 ... The Pontryagin's Maximum Principle then states that, under the assumptions: • ∃ (x, u) feasible for (P ), with x absolutely continuous and u measurable.

• f and l are continuous with respect to u and C 1 with respect to t and x.

• g, ψ 0 , ψ 1 are C 1 with respect to x.

Let (x, u) be an optimal pair for (P ), then (i) ∃ p = 0 absolutely continuous such that we have the Hamiltonian system

   ẋ = ∂H ∂p (t, x, p, u) ṗ = -∂H ∂x (t, x, p, u) (ii) u is solution of M in w∈U H(t, x, p, w) ae in [t 0 , t f ]. (iii) "Transversality conditions": ∃ (µ 0 , µ 1 ) such that (T C)            ψ 0 (x(t 0 )) = 0 p(t 0 ) = -∂Φ ∂x0 (t 0 , x(t 0 ), t f , x(t f ), µ 0 , µ 1 ) ψ 1 (x(t f )) = 0 p(t f ) = ∂Φ ∂x f (t 0 , x(t 0 ), t f , x(t f ), µ 0 , µ 1 ) with Φ : (t 0 , x 0 , t f , x f , µ 0 , µ 1 ) → g(t 0 , x 0 , t f , x f ) + (ψ 0 (t 0 , x 0 )|µ 0 ) + (ψ 1 (t f , x f )|µ 1 )
Remark: this is why indirect methods are sometimes referred to as necessary condition methods.

Now we denote y = (x, p) and ϕ the state-costate dynamics derived from the Hamiltonian system. We assume here that the expression of the optimal control given by the necessary conditions (Hamiltonian minimization) is actually a function, noted γ, ie u(t) = ArgM in w∈U H(t, x, p, w) = γ(t, y). Solving (P ) is equivalent to solving the following Boundary Value Problem ‡ (BV P )

   ẏ = ϕ(t, y, γ(t, y)) ae in [t 0 , t f ] c 0 (t 0 , y(t 0 )) = 0 Boundary Conditions at t 0 c 1 (t f , y(t f )) = 0
Boundary Conditions at t f Note: these Boundary Conditions c 0 and c 1 correspond to the Transversality Conditions mentioned above, that contain the Initial and Terminal conditions of (P ) in addition to the constraints on the costate p.

It is possible to integrate y = (x, p) if we set the value of y(t 0 ), and we then obtain the following Initial Value Problem

(IV P ) ẏ = ϕ(t, y, γ(t, y)) ae in [t 0 , t f ] y(t 0 ) = ζ Initial V alue
We introduce now an application called the shooting function, which basically maps the initial value ζ to the value of the Boundary conditions at t f for the corresponding solution of (IV P ). In practice, the initial conditions ψ 0 of the problem (P ) already give a part of y(t 0 ), so the unknown of the shooting function is reduced to the "missing" part, that we note z. A frequent situation is when the initial conditions determine the initial state x(t 0 ), therefore z is actually the initial costate p(t 0 ). Then the value of the shooting function is given by the boundary conditions at t f for the solution y(•, z) of (IV P ) corresponding to the initial value y(t 0 , z) = (x 0 , z)

(Shooting f unction) S : z → c 1 (y(t f , z)).
Finding a zero of the shooting function S is then equivalent to the resolution of (BV P ), and therefore also gives a solution of (P ). The "shooting method" thus consists in solving the equation S(z) = 0, as summarized below:

Optimal Control Problem (P ) -→ Boundary Value Problem (BV P ) -→ Initial Value Problem (IV P )
and Shooting function S

We consider now the discontinuous case, when the optimal control presents switchings (or commutations). More precisely, the Hamiltonian minimization gives

u = γ 1 (t, y) , if ψ(t, y) < 0 u = γ 2 (t, y) , if ψ(t, y) > 0 u = Γ(t, y)
, if ψ(t, y) = 0 with Γ(t, y) a subset of U , and ψ the switching function, whose zeros correspond to the commutations of the optimal control. We assume that the set of switching times t such that ψ(t, y) = 0 is finite. We then obtain y(•, z) as solution of the initial value problem

(IV P ) disc    ẏ = ϕ 1 (t, y) , if ψ(t, y) < 0 ẏ = ϕ 2 (t, y) , if ψ(t, y) > 0, y(t 0 , z) = (x 0 , z)
We would like here to highlight some difficulties related to the single shooting method in presence of a large number of commutations, namely the evaluation of the shooting function and its Jacobian.

The orbital transfer problem

We study here an orbital transfer problem, originally submitted by CNES § . We consider a satellite with a mass of 1500 kg and low thrust electro-ionic propulsion (with thrusts ranging from 10 Newtons to 0.1 Newton). We want to transfer it from a strongly elliptic, slightly inclined orbit, to a circular, equatorial, geostationary orbit. The objective is to maximize the payload, i.e., to minimize the fuel consumption during the transfer (NB: this is not a minimumtime problem). This kind of problems was studied for instance in [START_REF] Bartholomew-Biggs | The solution of some difficult problems in low-thrust interplanetary trajectory optimization[END_REF][START_REF] Oberle | Existence and Multiple Solutions of the Minimum-Fuel Orbit Transfer Problem[END_REF], and this specific family of problems in [START_REF] Gergaud | Homotopy method for minimum consumption orbit transfer problem[END_REF][START_REF] Gergaud | Low thrust minimum-fuel orbital transfer: an homotopic approach[END_REF][START_REF] Martinon | Numerical resolution of optimal control problems by a Piecewise Linear continuation method[END_REF][START_REF] Haberkorn | Transfert orbital à poussée faible avec minimisation de la consommation: résolution par homotopie différentielle[END_REF].

Model

We consider that the forces applied to the satellite are the Earth attraction (central force) and the engine thrust: r = -µ r |r| 3 + T m with r the position vector in R 3 , T the thrust also in R 3 (ie the control), m the mass of the satellite, µ the gravitational constant of the Earth (µ = Gm T , with G the universal gravitational constant and m T the mass of the Earth), and | • | the Euclidean norm.

State

For a better numerical stability, we use orbital parameters instead of Cartesian coordinates to define the state variables, to avoid the strong oscillations that would result from the high number of revolutions. These parameters define the shape of the orbit and the position of the satellite on it, as illustrated on Figure 1 (where w is the true anomaly, Ω n the ascending node longitude, ω the argument of perigee and i the inclination with respect to the equatorial plane). More precisely, the new state variables are:

• Orbit parameter P • Eccentricity vector (e x , e y ), in the orbit plane, oriented towards perigee • Rotation vector (h x , h y ), in the equatorial plane, collinear to the intersection of orbit and equatorial planes

• True longitude L = Ω n + ω + w • Mass m
The eccentricity and rotation vectors are given by (e is the orbit eccentricity): e x = e cos(Ω + ω) , e y = e sin(Ω + ω) h x = tan(i/2) cosΩ , h y = tan(i/2) sinΩ

Control

The normalized three-dimensional control u such that T = T max u is expressed in the moving reference frame attached to the satellite, as radial thrust q, transverse thrust s and normal § Centre National d'Etudes Spatiales Copyright c 2000 John Wiley & Sons, Ltd.

Optim. Control Appl. Meth. 2000; 0:0-0 Prepared using ocaauth.cls

Optimal Control Problem

If we note x = (P, e x , e y , h x , h y , L), with T M ax the maximal thrust we have the following formulation of the problem:

(P )                    M in t f t0 |u(t)| dt ẋ = f 0 (x) + TMax m B(x)u ṁ = -βT M ax |u| |u| ≤ 1 x(t 0 ) = (11625, 0.75, 0, 0.0612, 0, π) m(t 0 ) = 1500 x(t f ) = (42165, 0, 0, 0, 0, f ree) m(t f ) is f ree t 0 = 0 t f is f ixed with f 0 (x) = µ P 0 0 0 0 0 W 2 P T
Copyright c 2000 John Wiley & Sons, Ltd.

Optim. Control Appl. Meth. 2000; 0:0-0 Prepared using ocaauth.cls and

B(x) = P µ         0 2P/W 0 sinL cosL + (e x + cosL)/W -Ze y /W -cosL sinL + (e y + sinL)/W Ze x /W 0 0 (CcosL)/(2W ) 0 0 (CsinL)/(2W ) 0 0 Z/W        
where

W = 1 + e x cosL + e y sinL , Z = h x sinL -h y cosL , C = 1 + h 2 x + h 2 y
The physical constants µ and β are respectively 398600.47km 3 s -2 and 0.05112km -1 s.

Remark: the presence of the norm of the control in the mass dynamic would lead to a non differentiable problem with a direct method.

Switching function and Control

We note p the costate associated to x, and p m the costate associated to m, and define the switching function ψ:

ψ(x, m, p, p m ) = 1 -βT max p m - T M ax m |B(x) t p|.
We now make two assumptions (that we numerically check, see [START_REF] Haberkorn | Transfert orbital à poussée faible avec minimisation de la consommation: résolution par homotopie différentielle[END_REF]):

• B(x) t p is non-zero on [t 0 , t f ]
• There is no singular arc, ie ψ does not vanish on any finite interval.

The application of Pontryagin's Maximum Principle then leads to the following expression of the optimal control

u = -B(x) t p |B(x) t p| if ψ(x, m, p, p m ) < 0 u = 0 if ψ(x, m, p, p m ) > 0
We can see that this control has a bang-bang structure, as its norm switches between 0 and 1 at zeros of the switching function ψ. The two cases define the two dynamics ϕ 1 and ϕ 2 .

Shooting function

In this context, we can prove that the shooting function S is smooth on a suitable open subset of R n . More precisely, if we consider a security zone, corresponding to the fact that the satellite does not come too close to Earth, and does not use all of its fuel (note: the following definitions use the Cartesian coordinates, the state being decomposed as (position,velocity,mass)):

Let A = {x = (r, v, m) ∈ R 7 |r > ρ 0 > 0 , m > M 0 > 0}
where ρ 0 is a safety minimal altitude and M 0 the minimal final mass (for instance the mass of the satellite without fuel).

Let us now define

Ω = {z ∈ R n | x(t, z) ∈ A, ∀t p v (t, z) = 0, ∀t ψ(y(t, z)) 2 + (p r (t, z)|p v (t, z)) 2 = 0, ∀t ψ(y(0, z)) = 0 and ψ(y(t f , z)) = 0 } Then we have i) for z ∈ Ω, t → ψ(y(t, z)
) is continuously differentiable, and finitely many switching times t i exist such that ψ(y(t i , z)) = 0. ii) for every commutation, the flow is transverse to the commutation surface, and goes through the surface iii

) Ω is an open subset of R n ans S is C ∞ on Ω.
Proof. see [START_REF] Gergaud | Homotopy method for minimum consumption orbit transfer problem[END_REF].

Continuation approach

It is well known that shooting methods are quite sensible to the initial guess, and it is indeed the case for this family of problems. We therefore use continuation approaches to obtain a suitable starting point (we refer readers interested in continuation methods to [START_REF] Allgower | Simplicial and continuation methods for approximating fixed points and solutions to systems of equations[END_REF][START_REF] Allgower | Numerical Continuation Methods[END_REF][START_REF] Allgower | Piecewise linear methods for nonlinear equations and optimization[END_REF]). The principle is to find a suitable homotopy that connects the problem to an easier one, and then to follow the zero path of this homotopy, from the easy problem to the original one. In our case, we choose for the "easy" problem the same transfer with minimization of an "energy" criterion, namely

t f t0 |u(t)| 2 dt.
For our homotopy, we take the shooting function corresponding to the orbital transfer with the following objective, parametrized by λ ∈ [0, 1]:

J λ = t f t0 λ|u(t)| + (1 -λ)|u(t)| 2 dt.
The resulting perturbed problem (P λ ) has a strictly convex Hamiltonian (with respect to u), with a continuous optimal control, and is much easier to solve than (P ) = (P 1 ). The path following requires some precautions, but this approach provided us with a good initial point for the single shooting method, with no a priori assumption on the structure of the optimal control, see [START_REF] Gergaud | Homotopy method for minimum consumption orbit transfer problem[END_REF][START_REF] Gergaud | Low thrust minimum-fuel orbital transfer: an homotopic approach[END_REF][START_REF] Haberkorn | Transfert orbital à poussée faible avec minimisation de la consommation: résolution par homotopie différentielle[END_REF][START_REF] Martinon | Numerical resolution of optimal control problems by a Piecewise Linear continuation method[END_REF].

Remark: For a better numerical behavior, we also use longitude instead of time as the integration variable. The resulting problem is non-autonomous, with a fixed final longitude (ie "rendez-vous") and free final time. Along with the choice of the orbital parameters, this leads to a very good numerical stability of the system. This allows us to solve the problem with a single shooting method, without having to resort to multiple shooting.

Methods

Now let us go back to some specific difficulties for the shooting method in this context. First, the integration of (IV P ) disc itself can be tricky, and typically requires a good variable step integrator. Then, the use of variable step can hinder the computation of the Jacobian by classical finite differences (cf [START_REF] Bock | Numerical treatment of inverse problems in chemical reaction kinetics[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations I: Non-stiff Problems[END_REF]). We propose to study the behavior of some practical ways to deal with these difficulties.

Switching detection

First, instead of letting the variable step integrator deal with the discontinuities of the right hand side on its own, we try to detect the switchings during the integration of (IV P ) disc . The principle of the detection is described below, and assumes that the integration method used is able to provide a "dense output", ie a cheap approximation of y on each integration interval. This is typically done by some polynomial interpolation, see [START_REF] Hairer | Solving Ordinary Differential Equations I: Non-stiff Problems[END_REF] for a thorough description of dense outputs available for many integration methods.

At integration time t, compute the switching function ψ(t, y).

If the sign of ψ has changed since the previous step, a commutation has occurred (actually, an odd number of commutations, but most often one in practice).

Then we locate the commutation by dichotomy, by using the dense output of the integrator to provide y and ψ. Then we perform the control switch, store the sign of ψ, and restart the integration back from the commutation time.

Else the integration goes on normally. This algorithm is quite simple, but effective. The main issue is the risk of skipping one or several commutations during an integration step. An easy improvement of the detection method is to perform additional sign checks at a fixed number of intermediate points, and not only at the end of the step. If these additional checks bring some new commutations, then the number of intermediate checkpoints should be increased until no new commutations appear. Besides, with a variable step integrator, we can also expect the stepsize control mechanism to reject a step that would have missed commutations. In practice, for the orbital transfer studied in this paper, setting 10 intermediate checkpoints appears to be sufficient to catch all commutations (up to 1800 for the longest transfer with the lowest thrust).

Jacobian evaluation -variational equations

The first way to compute the Jacobian of the shooting function is simply to use finite differences. We use for the shooting method the nonlinear solver hybrd (see code references on page 20), which performs finite differences with a step of

h j = √ |x j |,
where is the error on the shooting function evaluation (for the numerical experiments, we set to the relative error of the integrator). This method will be noted "external differentiation" (END) in the following. The drawback of this method is that when used with a variable step integrator, the integration steps can vary at the points where S is evaluated for the finite differences, which can impair the approximation of the Jacobian (see [START_REF] Bock | Numerical treatment of inverse problems in chemical reaction kinetics[END_REF], [START_REF] Hairer | Solving Ordinary Differential Equations I: Non-stiff Problems[END_REF] p. 200).

Then, another possibility is to use the variational equations to compute the Jacobian of the shooting function. We first recall the smooth case, when we consider the derivative with respect to the initial condition of the system

(IV P )    ẏ(t) = ϕ(t, y(t)) y i (t 0 ) = y 0 i for i = 1, . . . , n y i+n = z i for i = 1, . . . , n.
If we note y(•, z) the solution of (IV P ), we know that ∂y ∂zi (t f , z) is solution of the variational system

(V AR) j Ẏj (t) = ∂ϕ ∂y (t, y(t))Y j (t) Y j (t 0 ) = 0 • • • 1 • • • 0 T
It is then possible to compute the Jacobian of the shooting function by integrating (V AR) along with (IV P ). We will use two methods to compute the right hand side of (V AR). First, we can just obtain ∂ϕ ∂y by automatic differentiation (see code references on page 20). Another possibility is to approximate the right hand side by finite differences (cf [START_REF] Hairer | Solving Ordinary Differential Equations I: Non-stiff Problems[END_REF] p.201):

Ẏ (:, j) ≈ 1 h (ϕ(y + hY (:, j)) -ϕ(y))
We will use the notations VARad and VARfd to distinguish these two methods.

Remark on the error orders (see [START_REF] Hairer | Solving Ordinary Differential Equations I: Non-stiff Problems[END_REF]): For the external differentiation (END), the error of the finite differences with a step of h on the shooting function is of size O(T ol/h) + O(h), where T ol is the tolerance of the integrator used for (IV P ). Therefore, it seems natural to chose h ≈ √ T ol, to obtain an error of size O( √ T ol).

For the variational formulations, the Jacobian is obtained from the integration of the system (V AR). If we use automatic differentiation for the right hand side (VARad), then the error is of size T ol. If finite differences are used (VARfd), the error is of size O(h) + O(macheps/h), where macheps is the machine precision. Choosing h ≈ √ macheps gives an error of size O( √ macheps) on the right hand side of (V AR), for a final error on the Jacobian of size O(T ol + √ macheps). Thus for T ol ≥ √ macheps, the method VARfd is as precise as VARad.

A simple example: We illustrate these different methods on a simple example:

(P ε )            M in 2 0 (|u(t)| -ε(ln |u(t)| + ln(1 -|u(t)|)))dt ẋ(t) = -x(t) + u(t) |u(t)| < 1 x(0) = x 0 = 0 x(2) = x f = 0.5
We consider ε = 0.01, and use Matlab's ode45 to integrate the corresponding (IV P ), with a tolerance of T ol = 10 -4 . The solution is z * = -0.2715, and the corresponding trajectory and We trace on Figure 4 the derivative of the shooting function ∂S ∂z , as computed by the external differentiation (with a finite difference step of T ol and √ T ol), and the two variational formulations VARad and VARfd. According to the formulas recalled above, we expect the error on the Jacobian to be of size O(1) and O(10 -2 ) for the two END, and O(10 -4 ) for both VARad and VARfd. The graphs confirm that the first END (h = T ol) is the worst of all, followed by the second END (h = √ T ol), the two variational methods VARad and VARfd being the most accurate, and of similar precision. Note: The two peaks observed on S at z ≈ +/ -0.14 correspond to a strong variation (with respect to z) of the control values near t f = 2, which result in a strong variation of S(z).

Variational system -Discontinuous case

Let us now describe the use of the variational equations in presence of commutations. For a fixed z, let us denote by t 0 < t 1 < . . . < t q < t f the q switching times, and decompose the shooting function S(z) = S q • . . . • S 0 (z), with the applications

S 0 : Ω ⊂ R n -→ Ω 1 ⊂ R × R 2n z -→ (t 1 (z), y(t 1 (z), z)), S q : Ω q ⊂ R × R 2n -→ Ω q+1 ⊂ R 2n (t i , y i ) -→ y(t i+1 (t i , y i ), t i , y i ). and S i : Ω i ⊂ R × R 2n -→ Ω i+1 ⊂ R × R 2n (t i , y i ) -→ (t i+1 (t i , y i ), y(t i+1 (t i , y i ), t i , y i )), i = 1, . . . , q -1
The implicit function theorem gives the derivatives of the switching time t i+1 (t i , y i ) Here Y ti (.) = ∂y ∂ti (., t i , y i ) and Y yi (.) = ∂y ∂yi (., t i , y i ) are solutions of the variational systems (we note in the following ϕ i the right hand side on ]t i , t i+1 [, which can be either ϕ 1 or ϕ 2 .)

t i+1 (t i , y i ) = ∂ψ ∂t (t,
(V AR) ti Ẏ (t) = ∂ϕi ∂y (t, y(t))Y (t) Y (t i ) = -ϕ i (t i , y i ) and (V AR) yi Ẏ (t) = ∂ϕi ∂y (t, y(t))Y (t) Y (t i ) = I
Note: as the variational system is linear with respect to Y, we have Y ti = -Y yi ϕ i (t i , y i ), so only Y ti needs to be integrated in practice.

We then have for

S i S i (t i , y i ) = ∂ti+1 ∂ti (t i , y i ) ∂ti+1 ∂yi (t i , y i ) ∂ti+1 ∂ti (t i , y i )ϕ i (t i+1 , y i+1 ) + Y ti (t i+1 ) ϕ i (t i+1 , y i+1 ) ∂ti+1 ∂yi (t i , y i ) + Y yi (t i+1
) And for S 0 and S q , we have S 0 (z) = ∂t1 ∂z (z) ϕ 0 (t 1 (z), y 1 (t 1 (z), t 0 , z)) ∂t1 ∂z (z) + ∂y ∂z (t 1 (z), t 0 , z) and S q (t q , y q ) = ∂y ∂tq (t f , t q , y q ) ∂y ∂yq (t f , t q , y q ) . This gives the derivatives of the shooting function S (z) = S q • . . . • S 0 (z) (note: a similar approach can be found for instance in [START_REF] Kaya | Computations and Time-Optimal Controls[END_REF]). As in the smooth case, the partial derivatives in the formulas can be computed either by automatic differentiation, or by finite differences, and we will keep the notations VARad and VARfd in the following.

An improved formulation: However, one might notice that the above formula actually computes Y yi (.) = ∂y ∂yi (., t i , y i ) on each sub-arc, which implies integrating a variational system of dimension 2n × 2n. The point is, we only need the derivatives with respect to the unknown z ∈ R n for the Jacobian of the shooting function. We will see that it is not necessary to compute the full derivatives with respect to y i , but only half of it, corresponding to the components of z.

Copyright c 2000 John Wiley & Sons, Ltd.

Optim. Control Appl. Meth. 2000; 0:0-0 Prepared using ocaauth.cls

We will now assume for the reasoning that the unknown z corresponds to the n last components of y = (x, p). This means that the unknown of the shooting function is the initial costate, which is often the case from the application of the PMP. In other cases, it would be sufficient to rearrange the components of y to put the n components corresponding to z at the end.

We consider the integration of y on the interval [0, T ] from the initial value y 0 , with a commutation occurring at τ ∈]0, T [. We assume here that the system is autonomous (it is always possible to transform a non-autonomous system into an autonomous one by adding a state component corresponding to the independent variable, of derivative equal to 1), and therefore replace the notation y(t, t 0 , y 0 ) by y(t -t 0 , y 0 ).

Assuming the dynamic is ϕ 1 before τ and ϕ 2 afterward, the integrated value of y at T , denoted I(y 0 ), is given by: I(y 0 ) = y 2 (T -τ, y 1 (τ, y 0 )) , with τ = τ (y 0 ) determined by ψ(y 1 (τ, y 0 )) = 0.

With the notations y

τ = y 1 (τ, y 0 ), y T = y 2 (T -τ, y τ ), Y 1 (t, y 0 ) = ∂y1(t,y0) ∂y0 and Y 2 (t, y 0 ) = ∂y2(t,y0) ∂y0
, we can write the derivatives as

I (y 0 ) = Y 2 (T -τ, y τ )(Y 1 (τ, y 0 ) + ϕ 1 (y τ )τ (y 0 )) -ϕ 2 (y T )τ (y 0 )
where τ (y 0 ) can be obtained like before by the implicit function theorem:

τ (y 0 ) = - ∂ψ ∂y (y τ )Y 1 (τ, y 0 ) ∂ψ ∂y (y τ )ϕ 1 (y τ ) ( * )
If we integrate the variational system until the commutation, taking T = τ + in the formula above gives the derivatives at the commutation (as Y 2 (τ + -τ, y τ ) = I)

I (y 0 ) = Y 1 (τ, y 0 ) + (ϕ 1 (y τ ) -ϕ 2 (y τ + ))τ (y 0 )
Instead of reinitializing the integrated Y to I, we perform the update

Y ← Y + (ϕ 1 (y τ ) -ϕ 2 (y τ + ))τ (y 0 ).
Thus Y is actually equal to ∂y ∂y0 (τ + , y 0 ), and we have crossed the commutation while keeping the correct derivatives. We can then continue the integration of the variational system until the next commutation, and so on until we reach the prescribed final time. At this point the integrated Y holds the value of ∂y ∂y0 (t f , y 0 ), which immediately gives the Jacobian of the shooting function S.

A first advantage of this second formulation is to get rid of the composition matrix products that appear in S (z) = S q • . . . • S 0 (z). Moreover, it is now easy to see that the update formula above is still valid if we replace Y = ∂y ∂y0 by Y = ∂y ∂p0 and τ (y 0 ) by ∂τ ∂p0 . And it is clear that taking Y = ∂y ∂p0 in ( * ) precisely gives the value of ∂τ ∂p0 instead of τ (y 0 ). Recalling the assumption that the unknown of the shooting function is the initial costate, we therefore only need to integrate the variational system for Y = ∂y ∂p0 , which is of dimension 2n × n, half than for the first formulation. This improved formulation gives similar results to the previous one, but is about twice faster, which could be expected. The numerical results presented in the following correspond to this faster algorithm.

Remark on the complexity: With this second formulation, using finite differences to approximate the right hand side of the variational system has the same complexity as using finite differences to compute the Jacobian of the shooting function.

Finite differences to compute S (z) require n+1 evaluations of S, which corresponds roughly to (n + 1) × N steps evaluations of the dynamic ϕ, if we assume that the n + 1 calls to S nearly use the same number of integration steps.

Finite differences to compute the right hand side of the variational system require n + 1 evaluations of ϕ, thus leading to the same total complexity if the variational system nearly requires the same number of integration steps then the evaluation of S (ie the integration of (IV P ) disc ). This is obviously true with a fixed step integrator, but must be taken into account with a variable step integrator, as integrating the larger variational system will usually require more steps than the simple initial value problem for S. This is problematic, as we compute the derivatives of a function that numerically slightly differs from S, the y component being integrated with more steps than when we just compute the value of S. A way to correct this behavior is to force the stepsize control algorithm of the integrator to ignore all the components corresponding to Y , and take only into account the 2n components of y. This is easily done in practice if the integrator accepts vector (instead of scalar) tolerances.

Numerical results

We examine now the numerical results of the various methods described above, applied to the orbital transfer problem. The numerical experiments in the following use the dopri5 integrator from E. Hairer and G. Wanner (see [START_REF] Hairer | Solving Ordinary Differential Equations I: Non-stiff Problems[END_REF]), a 5th order Dormand-Prince embedded Runge-Kutta method with a 4th order dense output.

To begin with, we show on Figure 5 the solution for the transfer with a maximal thrust of 10 Newtons (Notes: The "+" signs on the state/costate graphs indicate the final conditions). From left to right, we have the following graphs, with the time as abscissa: -the state x ∈ R 8 : orbital parameters P, e x , e y , h x , h y , L, mass m and final time t f ¶ . -in the middle column, the corresponding costate -in the right column, the three components (q, s, w) and norm |u| of the control.

It can be seen that the evolution of the state variables is quite smooth, due to the use of orbital parameters instead of Cartesian coordinates, and also to the low thrust of the propulsion. As brutal corrections (like in the impulsional case) are not possible, the cheapest way to move the satellite to its final orbit seems to slowly and regularly correct the orbit along each revolution (more precisely at the apogee, see the trajectories on Figures 6 and7 below). ¶ The original problem with free final time was reformulated as a fixed final time one via the usual transformation t = t f .s, with s ∈ [0, 1]. This adds the components t f and pt f to the state and costate vectors respectively. In particular, the inclination regularly decreases until the end of the transfer, so it does not seem optimal to correct the inclination first and then proceed with a coplanar transfer, for instance. As expected, the control is bang-bang, with a succession of full thrust and drifting arcs, the commutations being clearly visible on the graph showing the control norm.

We see now on Figure 6 the trajectory projected on the equatorial plane, with the thrust and drifting arcs shown in red (darker) and green respectively. We observe that the thrust arcs are located at the apogees, and also at the last perigees, which may be linked to the fact that at the end of the transfer the orbit is nearly circular, so the distinction between apogee and perigee becomes thin. According to the criterion (consumption minimization), we can assume that thrust arcs occur where corrections to the orbit are the most effective.

For this thrust of 10N, the optimal control presents only 18 switchings, whereas the transfer for 0.1N (on Figure 7) has about 1800 switchings, but with a qualitatively similar trajectory and control structure. It actually happens that this family of problems is extremely regular with respect to the maximal thrust, with several empirical laws that are extremely well verified, most importantly T hrust(N ) × t f ≈ C te , see [START_REF] Haberkorn | Transfert orbital à poussée faible avec minimisation de la consommation: résolution par homotopie différentielle[END_REF] (a similar law was first observed for the minimum time problems, see [START_REF] Cong | Contrôle optimal et transfert en temps minimal[END_REF]). A similar law holds for the total number of switchings, and also for the total number of integration steps (for a given integration method). In the following, we will show the numerical results of the switching detection and variational equations methods on this family of problems. Integration is performed by the dopri5 code, with absolute and relative tolerances set to 10 -8 /10 -6 . The same tests, conducted with lower (10 -6 /10 -4 ) or higher (10 -12 /10 -10 ) precision, give qualitatively similar results. The only notable difference is that the basic external differentiation, relatively to the other methods, performs less badly with a high precision, and much worse with a low precision.

Switching detection

We begin the tests with the switching detection method, by comparing the results of the shooting algorithms with normal integration, and with the detection algorithm described on page 8. In both cases, the Jacobian of the shooting function is computed by the hybrd solver, by finite differences (external differentiation). Tables I and II show the results of the shootings with and without switching detection. The legend of the columns is: -T M ax : maximal thrust in Newtons (ranges from 10N to 0.1N) -Nf (Njac): number of shooting function S (resp. Jacobian) evaluations by the hybrd solver -Time: computational time in seconds (on a 2.4GHz Celeron) -|S|: norm of the shooting function at the solution -Objective: objective value (ie fuel consumption) at the solution -Steps, Rejected, ratio: total integration steps, rejected steps and total/rejected ratio First, we observe a better convergence with the detection, as indicated by the norm of the shooting function at the solution, and notice a strong decrease of the execution time, the detection method being two to four times faster than the standard integration. The objective value at the solution is almost identical, and the solution themselves are quite close (the scaled solution vectors are identical up to four digits at 10N, and two digits at 0.1N). We note that the number of function and Jacobian evaluations does not seem to be really affected by the detection.

The improvement is actually due to the great impact of the detection on the integration steps, with a reduction of the total steps by a factor 5 roughly, and of the rejected steps by a factor 10, the rejected step ratio falling from 50% (!) to 25%. These values are shown on Figure 8, with a logarithmic scale for the number of total steps, and in abscissa the thrust with an inverted logarithmic scale, the shortest transfer (10N) being on the left and the longest (0.1N) on the right. We can check that the empirical relation (for this problem family) T hrust(N ) × Steps ≈ C te is numerically well verified. This relation is a consequence of the T hrust(N ) × t f ≈ C te law mentioned previously, and holds for different integration methods (eg the well known rkf45), of course with different constants. Figure 9 illustrates the integrated trajectory (ie the accepted steps) for the 0.1N transfer, and the reduction of the number of steps is clearly visible. We notice that the steps are more clustered on the "left" and "right" parts of the orbit, which corresponds to the fact that thrust arcs are located around the apogees, and perigees at the end of the transfer. It seems natural that less steps are needed to integrate the drifting arcs in the "middle" of the graphs. This is confirmed by looking at the integration stepsize for the 10N transfer, as shown on Figure 10 (plain line for standard integration, bold line with switching detection, with a logarithmic scale for the stepsize). In the standard case, we observe three stepsize ranges: -relatively large steps ([0.5, 1] roughly), corresponding to the drifting (no thrust) arcs.

-smaller steps (about 10 -1 ), corresponding to the thrust arcs at the apogees and last perigees.

-extremely small steps (near 10 -6 ), just before the commutations (a lot of rejected steps also occur here, as the integrator reduces the stepsize until the error due to the commutation crossing is acceptable).

With the switching detection method, we see that the stepsize on the thrust and drifting arcs is nearly unaffected, whereas the very small steps at the commutations have almost disappeared. The switching detection therefore effectively does its job, and clearly improves the integration of (IV P ) disc for the computation of the shooting function. It should be kept in mind that this method's performance probably depends on the dense output of the integrator, which should be cheap and reasonably accurate (for the description of the dense output of dopri5 see [START_REF] Hairer | Solving Ordinary Differential Equations I: Non-stiff Problems[END_REF], page 192).

Variational equations

We compare now the two variational formulations VARad and VARfd to the switching detection method ENDc. We recall that both variational formulations actually also include the switching detection algorithm. The initial point and settings are the same as before. For the two variational formulations VARad and VARfd, we use a variant of hybrd (namely hybrj), whose only difference with hybrd is that it takes the user-supplied Jacobian instead of approximating it by finite differences. We indicate here on Table III function at the solution and the execution times.

The three formulations converge in all cases, and we note that the convergence is qualitatively similar for ENDc and the VAR methods. At the solutions, the objective value and the rejected step ratios are very close for the three formulations. These tests do not show any significant difference in the number of Jacobian evaluations between the three methods. Concerning the execution times, we see that the automatic differentiation (VARad) is slower than the other two. The times for ENDc and VARfd are close, in accordance with the remark on complexity on page 13, and the variational formulation seems a bit faster, which might indicate a better Jacobian approximation.

We now have a look at the Jacobian computed by the various formulations at the initial point of the shooting. We show on Figure 11 the Froebinius norm of the difference between the Jacobian of each formulation. More specifically, we trace the following distances: END to ENDc, END to VARad, ENDc to VARad, and VARfd to VARad.

We see that the two variational formulations are the closest ones, which is not surprising. Their distance ranges from 10 -4 to 10 -2 , and grows when the thrust decreases, probably because the transfer time and number of commutations increase. The next shorter distance is between ENDc and VARad, between 10 -1 and 1 roughly. The distance between ENDc and VARfd is obviously of the same order, which is why it was omitted on the graph. On the other hand, the greatest distance is between END and the other three formulations (same remark as before, the END-VARfd has been omitted, and END-ENDc is very close to END-VARad...).

Remark: taking the formulas for the continuous case recalled on page 9, they would give an error size of O( √ T ol) for the two END methods and O(T ol) for the two VAR methods, with T ol = 10 -6 . Although these results do not hold for the discontinuous case, it is interesting to note that the difference between the Jacobians of the END and VAR methods seems actually close to √ T ol = 10 -3 .

Concerning the Jacobian numerical computation, it seems reasonable to assume that the variational methods provide the most accurate value, followed by the external differentiation with switching detection (ENDc), the basic external differentiation (END) being the furthest one. But while the improvement given by the switching detection is quite noticeable on the shooting results, the use of the variational equations does not seem to further improve the convergence. The approximation of the Jacobian by finite differences, when coupled with the switching detection, is here sufficient to ensure a good convergence of the quasi-Newton solver.

Conclusion

As a conclusion, we observe that on this family of problems with a large number of commutations, a switching detection method relying on the integrator's dense output can significantly improve both the convergence and the execution time of the shooting method. Besides, we see that using the variational equations to obtain the Jacobian of the shooting function achieves the same qualitative results in terms of convergence, with a slight improvement in terms of execution times for the finite differences formulation (VARfd).

For the automatic differentiation, we used tapenade from INRIA (http://wwwsop.inria.fr/tropics/tapenade.html, see [START_REF] Hascoët | TAPENADE: a tool for Automatic Differentiation of programs[END_REF][START_REF] Hascoët | TAPENADE 2.1 user's guide[END_REF] for more details).

Figure 1 .

 1 Figure 1. Orbital parameters.

Figure 2 .

 2 Figure 2. Three-dimensional control.

Figure 3 .

 3 Figure 3. Trajectory and control at the solution z * = -0.2715, and shooting function S on [-1, 1] .

Figure 4 .

 4 Figure 4. Illustration of the various Jacobian computation methods.

Copyright c 2000 6 Figure 5 .

 200065 Figure 5. Orbital transfer: solution for a 10N thrust (abscissa is time for all graphs).

Figure 6 .

 6 Figure 6. Orbital transfer: trajectory for a 10N thrust.

Figure 7 .

 7 Figure 7. Orbital transfer: trajectory for a 0.1N thrust.

Figure 8 .

 8 Figure 8. Total integration steps, and rejected steps ratio.

Figure 9 .

 9 Figure 9. Trajectories for the 0.1N transfer.

Figure 10 .

 10 Figure 10. Integration stepsize for the 10N transfer.

Copyright c 2000

 2000 John Wiley & Sons, Ltd. Optim. Control Appl. Meth. 2000; 0:0-0 Prepared using ocaauth.cls

Figure 11 .

 11 Figure 11. Distance between the Jacobian -10N to 0,1N transfers.

  i+1 , y(t i+1 , t i , y i )) Y ti (t i+1 ) ∂ψ ∂y (t i+1 , y(t i+1 , t i , y i )) Y yi (t i+1 )

	y(t, t i , y i )) +	∂y ∂ψ	(t, y(t, t i , y i )). ẏ(t, y(t, t i , y i ))	-1
	∂ψ ∂y (t

Table I .

 I Shooting results without switching detection

	TMax Nf (Njac) Time	|S| Objective Steps Rejected ratio
	10N	61 (5)	2 1.87 10 -6	121.21	846	411	0.49
	5N	72 (5)	7 4.20 10 -6	121.58	1683	813	0.48
	1N	59 (4)	30 1.06 10 -5	121.77	8383	4136	0.49
	0.5N	102 (7)	92 5.53 10 -4	121.69 16506	8083	0.49
	0.2N	48 (3)	116 2.49 10 -4	121.71 41360	20274	0.49
	0.1N	75 (5)	327 2.42 10 -4	121.70 82684	40435	0.49
		Table II. Shooting results with switching detection
	TMax Nf (Njac) Time	|S| Objective	Steps Rejected ratio
	10N	24 (1)	1 1.71 10 -14	121.21	167	33	0.20
	5N	121 (8)	5	1.55 10 -7	121.58	310	64	0.21
	1N	47 (3)	10	1.59 10 -8	121.78	1523	396	0.26
	0.5N	98 (7)	41	9.35 10 -7	121.69	2982	765	0.26
	0.2N	79 (5)	78	5.58 10 -5	121.71	7615	1909	0.25
	0.1N	33 (2)	65	3.61 10 -7	121.70 15455	3995	0.26

Table III .

 III the norm of the shooting Shooting results -Switching detection alone and variational equations

	10 -8/6		EN Dc		V ARad		V ARf d
	TMax	Time	|S| Time	|S| Time	|S|
	10N	1 1.71 10 -14	1 3.26 10 -14	1 2.73 10 -14
	5N	5	1.55 10 -7	5 6.13 10 -14	2 3.57 10 -13
	1N	10	1.59 10 -8	12	2.84 10 -8	19	2.50 10 -8
	0.5N	41	9.35 10 -7	93	3.63 10 -6	22	6.89 10 -7
	0.2N	79	5.58 10 -5	85	1.21 10 -4	56	6.29 10 -6
	0.1N	65	3.61 10 -5	114	4.16 10 -5	57	4.00 10 -5
	Copyright c 2000 John Wiley & Sons, Ltd.		Optim. Control Appl. Meth. 2000; 0:0-0
	Prepared using ocaauth.cls

Contract/grant sponsor: This work was partially funded by the CNES † , contract R&T 02/CNES/0257/

Codes references:

The code (shooting method and homotopy) used for all the presented numerical experiments is available at http://enseeiht.fr/~martinon/simplicial. See [START_REF] Martinon | Numerical resolution of optimal control problems by a Piecewise Linear continuation method[END_REF] for a complete description. The package includes the third party codes described below.

The integrator used for the experiments of this paper is dopri5 from E.Hairer and G.Wanner (see [START_REF] Hairer | Solving Ordinary Differential Equations I: Non-stiff Problems[END_REF]), available at http://www.unige.ch/~hairer/software.html.

The nonlinear solver for the shooting method is hybrd/hybrj from Burton s. Garbow, Kenneth e. Hillstrom, Jorge j. More (see [START_REF] Powell | A Hybrid Method for Nonlinear Equations[END_REF][START_REF] Garbow | Implementation Guide for MINPACK-1[END_REF][START_REF] More | User Guide for MINPACK-1[END_REF]), available at http://netlib.bell-labs.com/netlib/minpack/hybrd.f.gz.


Long Li

Werner Bauer

Etienne Mémin

Stochastic modeling of mesoscale eddies in barotropic wind-driven circulation

. The numerical assement shows that the proposed random model can well capture on a coarse mesh the correct four-gyre time-averaged circulation structure, as predicted by a direct numerical simulation at a much finer resolution.

dX t dt = u + ηt , ηt (x) = ∞ k λ k Φ k (x)ξ k (1)
where the mesoscale eddy is represented by a timeuncorrelated noise with ξ k ∼ N (0, 1) i.i.d., yet correlated in space with its covariance operator defined as:

Q t (x, y) = k λ k Φ k (x)Φ T k (y) (2)
2. Stochastic Reynolds Transport Theorem (SRRT). Following the framework of [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF], the rate of change of a scalar process θ within a volume, transported by the random flow ( 1) is given by

d dt V(t) θdV = V(t) Dθ Dt + θ∇•(u -∇•a) dV (3) Dθ Dt = ∂θ ∂t + (u + ηt )•∇θ - 3 i,j=1 ∂ 2 ij (a ij θ) (4)
where a(x, t) = dt 2 Q t (x, x) is defined as a one-point one-time variance tensor and the eddies are assumed to be incompressible, i.e. ∇• ηt = 0. 3. Derivation of QG Model.

Stochastic Barotropic Vorticity Equation (SBVE).

The governing equations for SBVE can be written in dimensionless form as:

∂q ∂t + J(ψ, q) = F + D + M (5) M = e 3 •∇× ( ηt •∇)u- 2 i,j=1 ∂ 2 ij (a ij u) +e 2 • ηt (6)
where q = Ro∇ 2 ψ + y is the potential vorticity in the beta plane with Ro the Rossby number, ψ is the mean velocity streamfunction, J is the Jacobian operator, F = sin(πy) is the symmetric double-gyre forcing, and D is an enstrophy-dissipation operator which can be given by ( δ 2 L ) 3 ∇ 4 ψ or ( δ 4 L ) 5 ∇ 6 ψ, with δ 2 , δ 4 the width of the Munk boundary layer. More detais on the physical mechanism is referred to [START_REF] Holm | Modeling mesoscale turbulence in the barotropic double-gyre circulation[END_REF]. We denote that, without any forcing and dissipation, SBVE conserves the mean kinetic energy along each realization:

d dt Ω 1 2 ∇ψ 2 dx = 0. (7)
On-line learning of eddy's EOFs from mean flow Method: 1) Sliding window to draw N o pseudoobservations at each grid point, removing average to substract fluctuations within each window size (L):

u L = u L -u L ;
2) Apply singular value decomposition (SVD) of fluctuations to get a set of empirical orthogonal functions (EOFs): η(L)

t = N o -1 k=1 √ λ k Φ (L) k ξ k ; 3) Adapt variance to grid scale (l) using a turbulence power-law [2]: η(l) t = ( l L ) 1/3 η(L) t .
Benefits: parameter-free, no need of data, timedependent EOFs, scaling auto-adapted.

Time-statistical nature of double-gyre circulation

Conclusion

The proposed random model has an explicit eddy representation which respects a set of physical consistency laws, such as SRRT (3) and energy conservation (7). From the numerical results of the time-statistics, we understand that the random model helps us to correct the wrong numerical dissipation of energy due to the Laplacian/Hyper-viscosity. It results that the physical four-gyre pattern in time average are recovered on a coarse-grid computation. In another work, the performance of our random model has been evaluated and analyzed in terms of uncertainty quantification and ensemble forecasting. Therein, results confirm that the spread is more accurate compared to a deterministic model with a pertubation of the initial condition. This ability is in particular essential for data assimilation applications. For future work, we attempt to move on to a multilayer model including baroclinicity. The objective will be to understand more about the isopycnal mixing through a specific parametrization.

Model 1 .

 1 Decomposition of Eddy-Mean Flow.

1 :

 1 Test of convergences for the stationary states -From left to right: time serie of kinetic energy; RMSEs of statistics for streamfunction by perturbation of time intervals; RMSEs of statistics for potential vorticity.

Figure 2 :

 2 Figure 2: Comparison of contour plots for different statistics: left top three -time-averaged streamfunctions; right top three -time-averaged potential vorticity; left bottom three -eddy kinetic energy (EKE); right top three -eddy potential enstrophy (EPE).


Marcus Denker

Mathieu Suen

Stéphane Ducasse

The Meta in Meta-object Architectures

est destinée au dépôt

Introduction

Reflection, a feature common to many modern programming languages, is the ability to query and change a system at runtime [START_REF] Kiczales | The Art of the Metaobject Protocol[END_REF]. Reflection is a very desirable feature, particularly suited for the long-living and highly dynamic systems of today.

The original model of reflection as defined by Smith [START_REF] Cantwell | Reflection and semantics in a procedural language[END_REF] is based on meta-level interpretation. The program is interpreted by an interpreter, such interpreter is interpreted by a metainterpreter leading to a tower of interpreters each defining the semantics of the program (the interpreter) it interprets.

However, a tower of interpreters is too slow in practice. To enable reflection in mainstream languages such as CLOS [START_REF] Kiczales | The Art of the Metaobject Protocol[END_REF], Smalltalk [START_REF] Ducasse | Evaluating message passing control techniques in Smalltalk[END_REF][START_REF] Rivard | Pour un lien d'instanciation dynamique dans les langages à classes[END_REF][START_REF] Röthlisberger | Unanticipated partial behavioral reflection: Adapting applications at runtime[END_REF] or Java [START_REF] Tanter | Partial behavioral reflection: Spatial and temporal selection of reification[END_REF], the tower of interpreters is replaced with a reflective architecture [START_REF] Maes | Concepts and experiments in computational reflection[END_REF] where meta-objects control the different aspects of reflection offered by the language. Meta-objects define the new or modified behavior and describe where this new behavior is active. For example, in systems that use metaclasses like CLOS [START_REF] Kiczales | The Art of the Metaobject Protocol[END_REF], Neoclasstalk [START_REF] Bouraqadi | Safe metaclass programming[END_REF], or MetaClassTalk [START_REF] Bouraqadi | Un MOP Smalltalk pour l'étude de la composition et de la compatibilité des métaclasses[END_REF], the metaclass of a class defines both the new behavior and which classes are affected by the new behavior. Recently, more elaborated schemes have been proposed (e.g., partial behavioral reflection [START_REF] Röthlisberger | Unanticipated partial behavioral reflection: Adapting applications at runtime[END_REF][START_REF] Tanter | Partial behavioral reflection: Spatial and temporal selection of reification[END_REF]) that provide a more flexible and fine-grained way to specify both the location been reflected and the meta-object invoked.

In general, meta-objects provide the implementation of new behavior which is called at certain defined places from the base system. It is important to note that meta-objects are not special objects, the execution of code as part of a metaobject is not different to any execution occuring in the base level application. As both base and metacomputation are handled the same, we are free to call any part of the base-system in the meta-level code.

As a matter of fact, this means that meta-level code can actually trigger again the execution of meta-level functionality. There is nothing to prevent the metalevel code to request the same code to be executed again, leading to an endless loop resulting in a system crash. This is especially a problem when reflecting on core language features (e.g., the implementation of Arrays or Numbers) since the chances are high that such features are used to implement reflective features themselves. These cases of spurious endless recursion of meta-object calls have been noted in the past in the context of CLOS [START_REF] Chiba | Avoiding confusion in metacircularity: The meta-helix[END_REF].

The ability to reflect on system classes is especially important when using reflection for dynamic analysis. A tracing tool that is realized with reflection should be able to trace a complete run of the system, not only the application code. In addition to the problem of recursion, such a tracer has the problem of recording the execution of trace code itself.

If we go back to the infinite tower (the origin of meta-level architectures) we can see that here these problems do not exist by construction: going meta means jumping up to another interpreter. A reflective function is always specific to one interpreter. As a function that is reflective at a meta-level I n is not necessarily reflective in I n+1 , the problem of infinite recursion does not happen.

An important question then is the difference between the meta-object and interpreter/infinite tower approach. The metaness in the case of the tower is defined by the interpreter used. The interpreter forms a context that defines if we are executing at the base level or at the meta-level. Calling reflective functionality (so called reification) is always specific to one interpreter. The meta-object approach now in contrast is lacking any mechanism to specify this contextual information: when executing a meta-level program, in a meta-object based reflective system, we lack the information that this program is executing at the meta-level. In addition, all reifications are globally active: we can not define to only trigger meta-object activation when executing base level code. The research question is then how can we incorporate the infinite tower property of explicitly representing the execution context into meta-object based architectures.

Our solution to this problem is to extend meta-object based reflective systems with a first class notion of meta-level execution and the possibility to query at any point in the execution whether we are executing at the meta-level or not.

To model meta-level execution we propose the notion of a first class context and context-aware reifications.

This paper is organized as follows. First we present a simple example to illustrate our problem and elaborate on the system used for evaluating the solution. Section 3 discusses the problem in detail, the next section then provides an overview of context and contextual reflection which can solve the presented problem (Section 4). We present an implementation in Section 5 followed by an evaluation (Section 6). After an overview of related work (Section 7) we conclude with a discussion of future work in Section 8

Context and Example

For the rest of the paper, we discuss the problem of meta-object call recursion in the context of Reflectivity, a reflective library implemented in Squeak [START_REF] Ingalls | Back to the future: The story of Squeak, a practical Smalltalk written in itself[END_REF]. We use Reflectivity to validate our solution. The ideas we discuss in this paper are the outcome of a larger project whose goal is to provide better control of reflection. Thus we started with a model of behavioral reflection that already allows for fine-grained spatial and temporal control of reification [START_REF] Röthlisberger | Unanticipated partial behavioral reflection: Adapting applications at runtime[END_REF][START_REF] Tanter | Partial behavioral reflection: Spatial and temporal selection of reification[END_REF]. Nevertheless it should be noted that the problem presented in this paper represents a universal problem: we are not just fixing a bug in our Reflectivity framework.

We first provide a short overview of partial behavioral reflection and discuss how to implement a simple example, which we use in the rest of the paper.

Partial Behavioral Reflection

With Reflex, Tanter introduced partial behavioral reflection [START_REF] Tanter | Partial behavioral reflection: Spatial and temporal selection of reification[END_REF] in the context of Java, the model was later applied to dynamic languages [START_REF] Röthlisberger | Unanticipated partial behavioral reflection: Adapting applications at runtime[END_REF]. Here metaobjects are associated not per object (as in 3-KRS [START_REF] Maes | Concepts and experiments in computational reflection[END_REF]) or per metaclass (as in CLOS [START_REF] Kiczales | The Art of the Metaobject Protocol[END_REF]), but per instruction. The core idea of the Reflex model is to bind meta-objects to operations using a link (see Figure 1). One can think about the link as the jump to the meta-level reified as an object. A link thus conceptually invokes messages on a meta-object at occurrences of selected operations.

Link attributes enable further control of the exact message sent to the metaobject. For example, we can control if the meta-object is supposed to be invoked before, after or instead of the original operation, an activation condition link attribute controls if the link is really invoked.

For our experiment, we use an implementation of partial behavioral reflection for Smalltalk that uses an abstract syntax tree for selecting which instructions to reflected on. Before execution, the AST is compiled on demand to a low-level representation that is executable, in our case to byte-codes executable by the Squeak virtual machine. More in-depth information about this system and its implementation can be found in the paper on sub-method reflection [START_REF] Denker | Submethod reflection[END_REF].

Annotating a node of the AST with a link thus results in code that when executed will call the specified meta-object.

A Simple Example

The problem of meta-object call recursion is a known problem [START_REF] Chiba | Avoiding confusion in metacircularity: The meta-helix[END_REF]. To show that it is relevant in practice, we have decided to keep our example as simple as possible. In Section 6.2 we discuss as more complex scenario how our solution is useful for dynamic analysis in general.

Imagine that we want to trace system activity: we want the system to beep when it executes a certain method. This audio based debugging is an interesting technique to determine if a certain piece of code is executed or not. In Squeak, there is a class Beeper that provides all the beeping functionality. When calling beep, the beep sound is played via the audio subsystem. The following example shows how to create a link that will invoke the message beep on the Beeper class.

beepLink := Link new metaObject: Beeper. beepLink selector: #beep Now we can install this beepLink on a method that is part of the system libraries. We take on purpose the method add: in OrderedCollection, a central class part of the collection libraries that is heavily used by the system. To set the link, we invoke the method link: on the AST-Node that stands for the whole method:

(OrderedCollection>>#add:) methodNode link: beepLink.

As result, a sound should be emited each time the method OrderedCollection»#add: is called. But as soon as we install this link, the system freezes. This clearly was not the intended outcome.

Infinite Meta-object Call Recursion

Let's analyze the cause for the problem presented above. After a discussion of some ad-hoc solutions, we show that the problem is caused by a missing model for the concept of the meta-level execution. The Problem. To ensure that the problem is not caused by our framework, we modify the example to call a different method at the meta-object, the method beepPrimitive which directly executes functionality defined in the virtual machine.

beepLink := Link new metaObject: Beeper. beepLink selector: #beepPrimitive.

When installing this link, we can see that it works as expected: we can hear a beep for all calls to the add: method, for example when typing characters in the Squeak IDE.

The problem thus lies in the code executed by the meta-object. The Squeak sound subsystem uses the method add: of OrderedCollection at some place to emit the beep sound. Thus, we call the same method add: from the meta-object that triggered the call to the meta in the first place. Therefore we end up calling the meta again and again as shown in Figure 2. This is clearly not a suitable semantics for behavioral reflection: it should be possible to use behavioral reflection even on system libraries that are used to implement the meta object functionality themselves. We present now two ad-hoc solutions.

Code Duplication. As the problem is caused by calling base level code from the meta-object, one solution would be to never call base level code from the meta-object, but instead provide a renamed copy for all classes and use these from the meta-level. Duplicating the complete system has, of course, all the standard problems of code duplication: space is wasted. In addition, the copy can easily become out of sync with the original. The problems could be minimized by just copying those methods that are really needed.

In practice, it is not easy to identify these methods, especially in dynamic languages. In addition this would cause changes in the reflective layer to become fragile because any change would require the programmer to update the copied version of the base level. This is clearly not a good solution. Adding Special Tests. Another solution could be to add special code to check if a recursion happens. The problem with this solution is that it is ad-hoc, the codebase becomes cluttered with checking instructions. It thus just patches the symptoms, the recursive call, and does not address the real problem.

The Real Problem: Modeling Meta-level Execution. The ad-hoc solutions are not satisfactory. The real problem exemplified by a recursive meta-call is that when going from infinite towers to meta-objects, the awareness that an execution occurs at the meta-level has been lost. Normally activation of a meta-object means a jump to the meta-level. The problem now is that this meta-level does not really exist in meta-object-based architectures: there is no way to query the system to know whether we are at the meta-level or not.

It should be noted again that the problem we have seen is not specific to a particular behavioral reflection framework. We observe the same problem when applying MethodWrappers [START_REF] Brant | Wrappers to the rescue[END_REF] to system classes. Method wrappers wrap a method with before/after behavior. MethodWrappers are reflectively implemented in Smalltalk and thus use lots of system library code during the execution of the wrapped methods. The same problem was identified for CLOS [START_REF] Chiba | Avoiding confusion in metacircularity: The meta-helix[END_REF] and is thus present in other meta-class based systems like for example MetaClassTalk [START_REF] Bouraqadi | Un MOP Smalltalk pour l'étude de la composition et de la compatibilité des métaclasses[END_REF].

Solution: A Metacontext

We have seen that the real cause for the problem of endless recursion lies in the absence of a model for meta-level execution: the fact that the system is executing meta-level code does not have a representation.

Modeling Context

At any point in the execution of some piece of code we should be able to query whether we are executing at the meta or at the base level. Such a property can be nicely modeled with the concept of context: the meta-level is a context that is active or inactive.

Such a context thus provides control-flow information: we want to know if we are in the control-flow of a call to a meta-object. But in addition, the context actually provides a reification: we have an object representing a specific controlflow, which in our case represents meta-level execution.

With a way to model meta-level execution, it is possible to solve the problem of recursive meta-object calls. A call to the meta-object can be scoped towards the base level: a meta-call should only occur when we are executing at the base level. If we are already at the meta-level, the calls should be ignored. This way, meta-object calls are only triggered by the base level computation, not the metalevel computation itself, thereby eliminating the recursion.

We will first describe a simplified model that only provides two levels of execution (base and meta) and does not allow any metameta-objects to be defined that are activated only from a meta-object execution. We describe later how to extend our model to support calls to these metameta-objects. The Metacontext. To model the meta-level, we introduce a special context, the MetaContext. This context is inactive for a normal execution of a program. MetaContext will be activated when we enter in the meta-level computation and deactivate when we leave it (Figure 3). The meta-context thus models meta-level execution.

Base Level

Meta Level

MetaContext activation

MetaContext deactivation

Fig. 3. The MetaContext activation

A simple model with just one meta-context is enough to distinguish the meta from the base level. We will see later that it makes sense to extend the metacontext to a possibly infinite tower of meta-contexts in Section 4.3. Controlling meta-object activation. Just having a way to model meta-level execution via the meta-context is not enough to solve the problem of recursion, it is just the prerequisite to be able to detect it. We need to make sure that a call to the meta-object does not occur again if we are already executing at the meta-level. Thus, the call to the meta-level needs to be guarded so it is not executed if the execution is already occurring at the meta-level. In the context of behavioral partial reflection (i.e., in the link-meta-object model that we used to show the problem), this means that the links are parameterized by the contexts in which they are active or not-active.

The Problem Revisited

With both the meta-context and the contextual controlled meta-object calls, we now can return to our example and see how our technique solves the problem of recursion. In our example, we defined the Beeper as a meta-object to be called when executing the add: method of OrderedCollection. The following steps occur (see Figure 4):

1. The add: method is executed from a base level program.

A call to the meta-object Beeper is requested:

-We first check if we are at the meta-level. As we are not, we continue with the call. -We enable the MetaContext.

-We call the meta-object.

Base Level Meta Level

Stop meta-level call Fig. 4. Stopping infinite meta-call recursion 3. Meta-object executes the beep method. 4. Meta-object calls the method add: method again 5. A Call to the meta-object is requested -We first check if we are at the meta-level. As we are executing meta-level code, the call is aborted. 6. Meta-object execution continues until it is finished. 7. On return to the base level, we deactivate the MetaContext.

Thus the recursive meta-call is aborted and the danger of recursion is eliminated. The model we described up to now with just one MetaContext is thus enough to solve the problem, but it not complete: it does not, for example, allow any calls to metameta-objects while already executing at the meta-level, which would make it impossible to observe or reason about metabehavior. In the next section, we therefore extend the model.

The Contextual Tower

As with the tower of interpreters, we can generalize the meta-context to form an infinite tower of meta-contexts. With the infinite tower of reflective interpreters, a reification is always bound to a specific interpreter. Normally, a jump from the base to the meta level means executing a reflective function that is defined as part of the interpreter I 1 . But it is possible to define a reflective function one level up: this then is only triggered by the interpreter I 2 that interprets I 1 , thus allowing to reflect on the interpreter I 1 itself. Figure 5 shows the reflective tower as visualized in the work of Smith [START_REF] Cantwell | Reflection and semantics in a procedural language[END_REF].

Transposed to our contextual model, it follows that having just one context (the meta-context) is not enough. We need more contexts for meta-meta, meta-3 and so on. If we have this contextual tower, we can for example define a metameta object that is only called when we are executing at meta-1. Meta-object calls need thus not only be defined to be active for the base level, but they can in the previous paragraph leads one to think of an infinite er of levels of reflective processors, each implementing the elow. 7 On such a view it is not coherent either to ask at level the tower is running, or to ask how many retlective s are running: in some sense they are all running at once. tly the same situation obtains when you use an editor ment, ed in APL. It is not as if the editor and the APL preter are both running together, either side-by-side or endently; rather, the one, being interior to the other, lies the anima or agency of /.he outer one. To put this er way, when you implement one process in another ss, you might want to say that you have two different sses, but you don't have concurrency; it is more a whole kind of relation. It is just this sense in which the r levels in our rcllective hierarchy are always running: of them is in some sense within the processor at the level , so that it can thereby engender it. We will not take a ipled view on which account --a single locus of agency ing between levels, or an infinite hierarchy of ltaneous processors --is correct, since they turn out to be viourally equivalent. (The simultaneous infinite tower of s is often the better way to understand processes, whereas i|!,ing-level viewpoint is sometimes the better way to rstand programs.) 3-Lisp, as we said, is an infinite reflective tower based on p. The cede at each level is like; the continuation-passing 2-MCP of Figure 14, but extended to provide a mechanism eby the user's program can gain access to fully articulated iptions of that program's operations and structures (thus ded, and located in a reflective tower, we call this code the p reflective processor). One gains this access by using what called reflective prncedures ~ procedures that, when ed, arc run not at the level at which the invocation red, but one level higher, at the level of the reflective ssor running the program, given as arguments those tures being passed around in the reflective processor.

(lambda simple [proc! ergs! cent]

(normalise (body procl) (bind (pattern proc!) args! (one procI)) cent))) Reflective procedures are essentially analogues of subroutines b be run "in tile implementation", except that they are in the same dialect as that being implemented, and can use all the power o(' the implemented language in carrying out their function (e.g., reflective procedures can themselves use reflective procedures, without limit). There is not a tower of different languages --there is a single dialect (3-Lisp) all the way up. optionally defined to be active for any of the meta-levels. This allows us to define meta*objects that reason about the system executing at any level, similar to the endless tower of interpreters.

L ve,,co . l''l J

As with all infinite structures, the most important question is how to realize it in practice. For the case of the infinite meta-context tower, there is an easy solution: contexts are objects, they can have state. We can parameterize the meta-context object with a number describing the meta-level that it encodes. Shifting to the meta-level means shifting from a context n to a context n + 1. Fig. 6. The Contextual Tower up to the second level Figure 6 shows an example. We have three contexts: the base level, the metalevel and meta-2. We see two links that are active at the base level. When called, they activate the meta-context with level 1. Then in the code of the meta-object (either the code of the meta-object itself or any library code executeted), we have a third link. This link is defined to be active only on meta-level 1, thus it will on execution enable meta-2. We show an example of such a level 1 link in Section 6.2.

An interesting and very nice property of the parameterized context is that the meta-contexts are only created on demand, they do not exist as an endless tower. This means, we have a form of partial reflection for providing a potential endless tower. If a meta-level is not needed, it does not cost anything.

The MetaContext Revised

The simple version of our idea with just one MetaContext is not enough to allow us to encode the Contextual Tower. We need a slightly modified model of MetaContext where the MetaContext is parameterized by the level. Such a parameterized MetaContext is not simply active or inactive, it is active for the level that it currently is set to. Thus when querying such a MetaContext, we give as a parameter a number denoting a meta-level. We will in the next section see how to realize such a context in practice.

Implementation

Now that we described the solution in general, we present an implementation of that model for Reflectivity, our reflection framework. We first show how the context is implemented and then discuss contextual links.

Implementation of MetaContext

The MetaContext is a class that has one instance per thread (a thread-specific singleton). The instances are created on demand and are stored per thread. As threads are objects in Squeak, we extended them to be able to store additional state directly in an associated dictionary. This mechanism is then used to store the MetaContext instance. Querying context. The MetaContext needs to model the meta-level. For that, it has one instance variable named level. We can increase the level by calling shiftLevelUp or decrease by shiftLevelDown. To test if the MetaContext is active for a certain level n we can call isActive: with a parameter denoting the level:

MetaLevel current isActive: 0

Sending current to the class MetaContext will retrieve the MetaContext singleton from the current process. If there is none yet, it will lazily create a Meta-Context with the level set to 0. Thus for a normal base level execution of code the expression above will return true.

Executing code in the MetaContext. To change the meta-level that code is executed at, we provide a way to run a block (an anonymous higher order function) one meta-level higher than the code outside the block. For example, to execute at meta-1, evaluate:

. ] valueWithMetaContext] valueWithMetaContext

The method valueWithMetaContext is implemented on the BlockClosure object, it will first shift the level of the current MetaContext up, then the block is executed. At the end the level of the context is shifted down to the previous value. We make sure that the downshift happens even in case of abnormal termination by evaluating the block using the exception handling mechanisms of Smalltalk: To make sure that the execution of the context handling code itself does not result in endless loops, we do not call any code from the sytem libraries in this method. Instead, we carefully copy all methods executed by the context setup code. The copied methods reside in the classes of the system library, but they are prefixed with ctxt and edited to call only prefixed methods. One example is the call to ctxtEnsure: seen in the method above. Concurrent meta-objects. As the MetaContext is represented by a threadspecific singleton, forking a new thread from the meta-level would mean that this thread has its own MetaContext object associated which is initialized to be at level 0. We have solved this problem by changing the implementation of threads to actually copy the meta-level information to the newly created thread. A thread created at the meta-level thus continues to run at the meta-level:

[ [ self assert: (MetaContext current isActive: 1) ] fork ] valueWithMetaContext
As the context information is not shared with the parent thread, the metalevel of the new thread is independed, it can continue to run at the meta-level even when the parent thread already returned to the base level.

Realizing Contextual Links

Now that we have a suitable way to represent the meta context, we need to make the links and the code that is generated to call them context-aware. For that, we need to solve three problems:

1. the link needs to be defined to be specific to a certain meta level.

2. link activation should occur only when code is executing on the right metalevel.

link activation should increase the meta-level.

The first problem is solved by a simple extension of the Link class, whereas the other two are concerned with the code that our system generates for link activation. We will now show the required changes in detail. Meta-level Specific Links. To allow the programmer to specify that a link is specific to a certain meta-level, we extend the link with a parameter called level. If level is not set, the link is globally active over all links (the standard behavior). The level can be set to any integer to define a link to only be active on that specific meta level. For our example, a link that is active only when executing base level code looks like this: beepLink := Link new metaObject: Beeper. beepLink selector: #beep; beepLink level: 0.

Context-Aware Link Activation. To jump up one meta-level on link activation, we make sure that the code generated for a link is wrapped in a call to valueWithMetaContext. The resulting code will look like this:

[ ...

code of the link ... ] valueWithMetaContext

In addition to that, we need to make sure that the link is only called when we are on the correct meta level. This is done by checking if the current MetaContext is executing on the same level as the level that the link is defined to be active. Only if this is true, we activate the link and call the meta-object. The code we need to generate looks like this: We will not go into the detail of how exactly the code is generated. The code can be found in the class GPTransformer of the Reflectivity distribution 2 .

Evaluation and Benchmarks

We first show that our solution solves the recursion problem, then we discuss how meta-context is useful for dynamic analysis. We describe dynamic meta-level analysis and realize an example. We present benchmarks to show the practicability of our approach.

The Problem is Solved

We show that we can really solve the practical problem. For that, we define the link that activates the Beeper to be specific to level 0: beepLink := Link new metaObject: Beeper. beepLink selector: #beep; beepLink level: 0.

Now we can install the link:

(OrderedCollection>>#add:) methodNode link: beepLink.

As soon as the link is installed, the next call to the add: method will trigger code-generation for that method. The code-generator will take the link into account and generate code as described earlier. Thus the recursive call to the add: method will not occur. We will thus hear a beep for every call to the add: method from the base level only.

Benefits for Dynamic Analysis

Our initial reason for modeling meta-level execution was to solve the problem of meta-object call recursion, and thus making reflection easier to use. But the solution we presented, modeling meta-level execution with the help of a metacontext, is useful far beyond just solving the problem of recursion. In this section, we will discuss what it means for dynamic analysis. We show how it allows the programmer to analyze the code executed by meta-objects by enabling dynamic meta-level analysis. Dynamic Analysis. One of the applications for behavioral reflection is dynamic analysis. For example, with reflection, it is fairly easy to introduce tracers or profilers into a running program [START_REF] Denker | Higher abstractions for dynamic analysis[END_REF] that normally require changes at the level of the virtual machine. One problem, though, with using reflection to introduce analysis code is that it is not clear which of the recorded events are resulting from the base level program execution and which from the code of the tracer itself. As long as we only trace application code, we can easily restrict the reflective tracer to the code of the application. But as soon as we want a complete system trace, we start to get into problems: recursion can occur easily (the problem we solved earlier), but even after working around recursion, we face another problem: how do we distinguish events that originated from the application from those that only occur due to the code executed by the tracer itself?

With our meta-level execution context, the problem described does not occur at all. The tracer (or any other tool performing dynamic analysis) is actually a meta-level program. A simple tracer would be the meta-object itself. More complex analysis tools would be called from a meta-object reifying the runtime event we are interested in. Thus, the code that performs the dynamic analysis is executing at the meta-level, while the links that trigger it are only active when executing a base level computation. This way we make sure that the infrastructure performing our analysis never affects the trace itself.

Dynamic Meta-level Analysis. An interesting challenge for dynamic analysis is that of analyzing the meta-level execution itself. Meta-level code should be lean and fast to not slow down the base computation more then really necessary. We thus are very interested in both tracing and profiling meta-level execution.

Our explicitly modeled meta-level execution makes this easy: we can define a link to be only active when executing at the meta-level. Therefore, we can install for example a trace-tool that exactly provides a trace of only those methods executed by a meta-object, even though the same methods are are called by the base level at the same time. We can thus as easily restrict the tracer towards meta-level execution as we restrict it to trace base level programs.

For our example, this means that we can use dynamic meta-level analysis to find the place in the sound-subsystem where the recursion problem happens when not using contextual links. In Section 3 we discuss that recursion happens, but we do not know where exactly the recursive call happens.

We define a link that is only active when we are executing code at level 1:

loggerLink := GPLink new metaObject: logger; selector: #log:; arguments: #(context); level: 1.

This link sends the message log: to the logger. The logger is an instance of class MyLogger:

logger := MyLogger new.
The method log: records the stack-frame that called the method where the link is installed on:

MyLogger>>log: aContext contexts add: aContext home sender copy

We install both a link calling the Beeper that is specific to level 0 and our link that is specific to level 1 on the method add:. beepLink := Link new metaObject: Beeper. beepLink selector: #beep; beepLink level: 0.

(OrderedCollection>>#add:) methodNode link: beepLink. (OrderedCollection>>#add:) reflectiveMethod methodNode link: loggerLink.

We can now inspect the logger object and see that it is recording the execution of SoundPlayer class»startPlayingImmediately: for every beep. Looking at this method, we find the code ActiveSounds add: aSound., which is the one spot in the sound system that calls the method add: of OrderedCollection. Thus we found with the help of dynamic meta-level analysis the exact call that causes the recursion problem as shown in Section 3.

Benchmarks

To assess if the system as presented is practically usable, we have carried out some benchmarks. Without the additional code for context activation, there is no overhead at all for calling a meta-object besides the call itself. The link specifies the meta-object and which method to call. The system generates from that information code that just calls the meta as specified. The problem now for analyzing the new context-enabled system is that the context code will, compared to the standard link activation, take a lot of time. In practice, though, meta-objects are usually there to do something: there is always code executed at the meta-level. So to make a practical comparison of the additional percentage of slowdown introduced by the context code, we need to compare the contextsetup code not only to the link activation of an empty meta-object, but to a meta-object that actually executes code.

We will benchmark the slowdown of the context handling for the execution of different meta-objects. The variation between the meta-objects is the number of sends done at the meta-level. For the benchmark, we create a class Base with just one empty method bench. To play the role of a meta-object, we create a class Meta with one method in which we call an empty method in a loop. This method thus simulates a meta-object doing some work. We can easily change the amount of work done by changing the loop. To know how this simple benchmark compares to real code executed, we added a meta-object calling the Beeper and one converting a float to a string.

We install a link on the method in Base to call the method in Meta. We call the base method now in a loop and measure the time: Table 1 shows the result when comparing both the original Geppetto and the context-enabled Geppetto for different meta-objects3 :

As expected, the slowdown in the case of an empty meta-object is substantial. But as soon as the meta-object itself executes code, the overhead starts to be acceptable. For calling the Beeper (from our running example), we have found an overhead of 63%. We are down to 17% on when executing 500 empty methods, and at 6.4% on 1000 methods.

It should be noted that this does not mean that the overhead observed for meta-object calls translate directly into a slowdown of a program using reflection. In a real program, the overall slowdown depends on how often meta-objects are called and how much time the program spends at the base level and meta-level compared to switching meta-levels. [START_REF] Maes | Concepts and experiments in computational reflection[END_REF], the CLOS MOP [START_REF] Kiczales | The Art of the Metaobject Protocol[END_REF] and other metaclass based systems like Neoclasstalk [START_REF] Bouraqadi | Safe metaclass programming[END_REF], or MetaClassTalk [START_REF] Bouraqadi | Un MOP Smalltalk pour l'étude de la composition et de la compatibilité des métaclasses[END_REF]. A more recent example is partial behavioral reflection for Java [START_REF] Tanter | Partial behavioral reflection: Spatial and temporal selection of reification[END_REF] and Smalltalk [START_REF] Röthlisberger | Unanticipated partial behavioral reflection: Adapting applications at runtime[END_REF]. None of these systems provide a representation for meta-level execution.

One could even argue that all these systems are not really reflective (as they can not be used to reflect on the system itself) and it is debatable if the metaobjects in these systems are really meta. It has already noted by Ferber [START_REF] Ferber | Computational reflection in class-based object-oriented languages[END_REF] that metaclasses are not meta in the computational sense, even though they are meta in the structural sense. For any computationally reflective system, we need to be able to decide if a computation is executing at the meta-level or not. Without this, there is no reason to talk about meta-objects at all. Aspect Oriented Programming. There is an ongoing debate about the relationship of aspect oriented rogramming (AOP) [START_REF] Kiczales | Aspect-Oriented Programming[END_REF] and reflection. Proponents of reflection claim that AOP is just a use-case for behavioral reflection, a pattern that can be easily implemented. The AOP community, though, would claim that reflection in turn is just the case of an aspect system where the domain is the language model itself, with the core-concepts (e.g. message sending) factored nicely into aspects and thus easily modifiable.

The ideas from this paper can contribute to this discussion. In AOP, a pointcut is globally visible: it matches even in the advice code itself by default. Conversely, we can say that the problem we noted in this paper is made an explicit and well known property in AOP, whereas it is a bug to be fixed in meta-object based reflection.

We claim that exactly here lies the difference between reflection and AOP: reflection needs the distinction between the base and the meta. Aspects, however, are a pure base level abstraction. Invoking an advice does not constitute a metacomputation, a level shift does not occur.

Contrary the notion that aspects are pure base level constructs, Stratified Aspects [START_REF] Bodden | Avoiding infinite recursion with stratified aspects[END_REF] define an extension that identify spurious recursion to be a problem and present the concept of metaadvice and metaaspects as a solution. This idea has some similarities to contextual reifications discussed in this paper, with the exception that meta-level execution is not modeled explicitly. Context-Oriented Programming. ContextL [START_REF] Costanza | Language constructs for context-oriented programming: An overview of ContextL[END_REF][START_REF] Hirschfeld | Context-oriented programming[END_REF] is a language to support Context-Oriented Programming (COP). The language provides a notion of layers, which package context-dependent behavioral variations. In practice, the variations consist of method definitions, mixins and before and after specifications. COP has no first-class notion of context, it is implicitly defined by the layers that are activated. The topic of reflection has been discussed for COP [START_REF] Costanza | Reflective layer activation in ContextL[END_REF]. But the paper looks at reflective layer activation, not a reflective model of context nor the use of context to structure or control reflection itself. Context-aware Aspects. The concept of context has seen some use in AOP [START_REF] Tanter | Context-aware aspects[END_REF]. As context specific behavior tends to crosscut base programs, it can advantageously be implemented as aspects. This leads to the notion of context-aware aspects, i.e., aspects whose behavior depends on context. The work has been continued in the direction of supporting reasoning on contexts and context history on the level of the pointcuts [START_REF] Herzeel | A temporal logic language for context awareness in pointcuts[END_REF][START_REF] Herzeel | Modularizing crosscuts in an e-commerce application in lisp using halo[END_REF].

Deployment strategies [START_REF] Tanter | Expressive scoping of dynamically-deployed aspects[END_REF] provide full control over dynamic aspect deployment both related to the call stack and to created objects and functions. As-pectBoxes [START_REF] Bergel | Aspectboxes -controlling the visibility of aspects[END_REF] are another example where aspects are controlled via a form of a context, in this case defined by the classbox model [START_REF] Bergel | Classboxes: Controlling visibility of class extensions[END_REF]. A good overview and discussion on the many mechanisms for dynamically scoping crosscutting in general is described in the work of Tanter [START_REF] Tanter | On dynamically-scoped crosscutting mechanisms[END_REF]. The MetaHelix. The problem of unwanted meta-level call recursion has been mentioned by Chiba and Kiczales [START_REF] Chiba | Avoiding confusion in metacircularity: The meta-helix[END_REF]. The problem discussed is first the structural problem that e.g., fields added by reflection to implement changed behavior show through to any user of introspection. The other problem mentioned is recursion, as any use of changed behavior can trigger a reification again. As a solution, the authors present the MetaHelix. All meta-objects have a field implemented-by that points to a version of the code that is not reflectively changed.

This approach is both more general and restrictive than our context based solution. It is more general, as it tries to solve the problem of the visibility of structural change. And it is more restrictive, as it does not model meta-level execution. The programmer has to call the right code explicitly, thus it can be seen as a controlled way to support the code copying solution presented in Section 3. The problem of structural changes is a very interesting one, as future work we plan to apply the ideas of the meta-context to structural reflection. Subjective Programming. Us [START_REF] Smith | A simple and unifying approach to subjective objects[END_REF] is a system based on Self that supports subject-oriented programming [START_REF] Harrison | Subject-oriented programming (a critique of pure objects)[END_REF]. Message lookup depends not only on the receiver of a message, but also on a second object, called the perspective. The perspective allows for layer activation similar to ContexL. The paper discusses the usefulness of subjectivity for controlling the access to reflection APIs, it does not go as far as using subjectivity for controlling behavioral reflection.

Conclusion and Future Work

In this paper we have analyzed the problem of the missing representation of metalevel execution in meta-object architectures. We have shown that the problem of infinite meta-object call recursion can be solved by introducing a representation for meta-level execution. We proposed to model the execution at the meta-level as a first class context and presented an implementation. Benchmarks show that the implementation can be realized in a practical manner.

For now, we have used the concept of context just to make the meta computation distinguishable from the base computation. We plan to extend the notion of contextual control of reification to other kinds of contexts then the meta-context.

We have experimented in the past with the idea of a first class model of change for programming languages [START_REF] Denker | Encapsulating and exploiting change with Changeboxes[END_REF]. We will explore the idea of context for structural reflection to model change. Virtual machine support for meta-contexts is interesting for two reasons. First, we hope to be able to improve performance by realizing all context setup code in the virtual machine. Second, as we explained in Section 5, the setup code executed when dealing with contexts has to be managed specially: we provide copies of all that code. We plan to move all this special code into the virtual machine.

An interesting question is how a context-aware reflective language kernel would look like and what the consequences for the language runtime and especially the reflective model would be. We plan to explore such a new reflective language kernel in the future.

Fig. 1 .

 1 Fig. 1. Partial Behavioral Reflection

Fig. 2 .

 2 Fig. 2. Infinite recursive meta-object call

Figure 14 :

 14 Figure 14: A Continaation-Passing 2-LISP MCP

Figure 15 :Fig. 5 .

 155 Figure 15: The 3-LISP Reflective Tower

[

  ... code executing on meta-1 ... ] valueWithMetaContext If code is already executing at meta-1, calling valueWithMetaContext again will execute at meta-2: [[ ... code executing on meta-2 ..

valueWithMetaContext

  MetaContext current shiftLevelUp. ˆself ctxtEnsure: [MetaContext current shiftLevelDown]

(

  MetaLevel current isActive: link level) ifTrue: [ [ ... code of the link ... ] valueWithMetaContext ].

[

  100000 timesRepeat: [Base new bench ]] timeToRun

Table 1 .

 1 Slowdown of meta-object calls with context

	meta-object	context (msecs) standard (msecs) slowdown
	0 message sends	614	34 1'705.88%
	10 message sends	723	165 338.18%
	50 message sends	1040	470 121.28%
	Beeper	1543	942	63.80%
	100 message sends	1406	856	64.25%
	200 message sends	2236	1621	37.94%
	1234.345 printString	2534	1920	31.98%
	500 message sends	4580	3907	17.23%
	1000 message sends	8543	8029	6.40%
	7 Related Work

Meta-object

Architectures. There are many examples for meta-object based reflective systems. Examples are 3-KRS

⋆ R.F. Paige and B. Meyer (Eds.): TOOLS EUROPE 2008, LNBIP 11, pp. 218-237, 2008.

The benchmark was run on an Apple MacBook Pro, 2.4Ghz Intel Core 2 Duo with 2GB RAM on Squeak Version 3.9

Acknowledgments. We thank Eric Tanter, Robert Hirschfeld, Orla Greevy, Oscar Nierstrasz, Martin von Löwis, Toon Verwaest, Adrian Kuhn, Adrian Lienhard and Lukas Renggli for discussing various aspects of the concepts presented in this paper. We acknowledge the financial support of the Swiss National Science Foundation for the project "Analyzing, capturing and taming software change" (SNF Project No. 200020-113342, Oct. 2006 -Sept. 2008) and COOK (JC05 42872) funded by the french Agence Nationale de la Recherche.