Eitan Altman

Alberto Avritzer

Rachid El-Azouzi

Leandro Pfleger De Aguiar

Daniel Sadoc Menasche

Leandro Pfleger

Rejuvenation and the Spread of Epidemics in General Topologies

Keywords: Performance analysis, Reliability, Security

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Epidemic models have received significant attention in the past few decades to study the propagation of viruses, worms and ideas in computer and social networks [START_REF] Newman | Networks: an introduction[END_REF], [START_REF] Huang | Software rejuvenation: Analysis, module and applications[END_REF]. In the case of viruses, the goal is to understand how the topology of the network and the properties of the nodes that comprise the network, together, impact the spread of the epidemics [START_REF] Newman | Networks: an introduction[END_REF], [START_REF] Avritzer | Using performance signatures and software rejuvenation for worm mitigation in tactical manets[END_REF]. If different treatments are available, at different costs, it is important to understand their impact on each individual node and on the network as a whole.

In this paper, we propose rejuvenation as a way to cope with epidemics. Rejuvenation consists of periodically sampling the nodes for early detection of viruses and preemptively running applications to prevent failures in the future. Full repair, in contrast, occurs after failure and requires significantly more effort.

We distinguish between a state in which the virus is incubating and in which symptoms might not be visible and yet they may be contagious and infecting other nodes, and a state of failure where symptoms are clear. One of the visible symptoms might be the decrease of capacity of nodes, also referred to as node degradation or failure. Aging of infected nodes produces degradation and allows the identification of failed nodes. Alternatively, nodes that carry the virus might be identified through sampling.

Sampling a node one can detect viruses and intrusions through their signatures [START_REF] Avritzer | Monitoring for security intrusion using performance signatures[END_REF]. In particular, viruses and security intrusions may be detected through identification of subtle changes in performance metrics (CPU, Memory, IO), also referred to as performance signatures, and the performance impact might depend on the sophistication of the attack. In some cases, it might be possible to detect a virus by analyzing performance signatures, even though the service capacity of the node as a whole does not show any level of degradation. In any case, sampling costs must be incurred if one wants to examine nodes to search for viruses and intrusions at an early stage. Using the proposed model, we numerically show the impact of sampling and rejuvenation on the steady-state population dynamics.

The two key questions posed in this paper are the following, 1) what is the impact of rejuvenation on the spread of epidemics and how does the effect of rejuvenation compare to other strategies such as relying only on full restoration after failures? 2) what is the effect of incubation time and how to quantify the effort for detecting viruses at an early stage?

Note that there is a vast literature on these topics [START_REF] Huang | Software rejuvenation: Analysis, module and applications[END_REF], but not on general topologies (graphs). Our objective is to understand the impact of topology on the spread of epidemics.

To answer these questions, we propose an analytical model whose main components are a graph, wherein nodes are connected to their neighbors through edges, and each node can be in one of five states: a) vaccinated, b) healthy, c) prone to failure (infected), d) under rejuvenation and e) failed. Infected nodes have the virus and can infect others, without showing immediate symptoms. Infected nodes that are not sampled remain infecting other nodes until one discovers that they are infected and a treatment strategy is applied. The discovery that a node is infected might occur either because the node is sampled or because it fails, in which case its performance degradation signalizes that the node requires treatment.

Using the proposed model, we numerically show the impact of sampling, rejuvenation, vaccination and full repair on the steady-state population dynamics. Our key contributions are the following:

Comparison of rejuvenation and full repair: we show how rejuvenation and full repair can be used to cope with the spread of epidemics, indicating their advantages and disadvantages. We account for sampling costs and benefits, in order to identify whether nodes are infected at an early stage. Although sampling incurs costs, it might prevent even larger healing costs, as well as avoid the spread of the virus while infected nodes incubate the virus without symptoms.

Analytical model: we present an analytical model which can be used to study the impact of rejuvenation on the spread of epidemics. Our model is tractable and allows to compute the transient and steady-state number of nodes in each of the possible states as a function of different system parameters. Using the proposed model, we establish conditions for the uniqueness of the optimal sampling rate, and we numerically indicate how different system parameters affect the optimal sampling rate.

The rest of this paper is organized as follows. In Section II we present a simple instance of our model, which already gives insights about the methodology used in the general model introduced in Section III. In Section IV we present numerical results, obtained with the proposed model, which relate the optimal sampling rate to different system parameters. Section V contains related work and Section VI concludes.

II. AGING AND INFECTION

To appreciate the key features of the proposed model, in this section we illustrate the proposed model in a simple setup accounting for aging, infection and full repair. Sampling costs and rejuvenation are not taken into account and the only treatment available is to fully repair failed nodes. In the next section, the simplifying assumptions are removed. The model presented in this section extends the model introduced in [2, Figure 1] to a network setting and general topologies.

Figure 1 illustrates the key components of our model. In state P the node is prone to failure but its performance remains the same as that of a healthy node, hence an external observer cannot distinguish between a healthy node and an infected node (except after careful examination also referred to as sampling). In state F , in contrast, the node has failed.

Each infected or failed node infects each of its neighbors at rate β. Let A ij be equal to 1 if node j is a neighbor of i and 0 otherwise. A healthy node i becomes infected by its neighbors at rate r i,2 (t), and by exogenous threats at rate κ i . An infected node further develops the infection and fails at rate λ. A failed node submitted to treatment heals at rate r 1 (Table I contains the notation). We assume that all times between events are exponentially distributed.

Let X i,0 (t), X i,P (t) and X i,F (t) be three indicator random variables equal to 1 if node i is healthy, prone to failure and failed, respectively, and 0 otherwise. Accordingly, let π i,0 (t), π i,P (t) and π i,F (t) be the probabilities that node i is healthy, prone to failure (infected) and failed at time t, respectively. Note that

E[X i,0 (t)] = π i,0 (t), E[X i,P (t)] = π i,P (t) and E[X i,F (t)] = π i,F (t).
The time change of E[X i,P (t)] is given by

dE[X i,P (t)] dt = E[-λX i,P (t)] + + E   X i,0 (t)   N j=1 A ji β(1 -X j,0 (t)) + κ i     (1)
The direct and exact solution of (1) requires the characterization of the expected value of X i,0 conditioned on all possible combinations of states of neighbors of node i, for i = 1, 2, . . . , N . The number of variables to be determined grows exponentially with respect to N . For this reason, in the remainder of this paper we rely on a mean-field (MF) approximation referred to as N-interwinded MF approximation (NIMFA) [START_REF] Van Mieghem | Virus spread in networks[END_REF] to compute the fraction of infected nodes. The approximation consists of replacing the expectation of the product of random variables by the product of their expectations. Let the endogenous infection rate of node i by its neighbors, at time t, be

r i,2 (t) = N j=1 A ji β(1 -π j,0 (t)) (2)
Applying the NIMFA approximation to (1) and making use of ( 2) we obtain

dπ i,P (t) dt = -λπ i,P (t) + π i,0 (t) (r i,2 (t) + κ i ) (3) dπ i,F (t) dt = λπ i,P (t) -r 1 π i,F (t) (4)
dπ i,0 (t) dt = - dπ i,P (t) dt + dπ i,F (t) dt (5)
π i,0 (t) = 1 -π i,P (t) -π i,F (t) (6)
Equation ( 3) characterizes the probability that node i is infected at time t, which decreases at failure rate λ in case node i is infected and increases at infection rate r i,2 (t) + κ i in case node i is healthy. Equation ( 4) characterizes the probability that node i is failed, which increases at rate λ in case node i is infected and decreases at repair rate r 1 in case node i is failed. Equations ( 5) and ( 6) capture the fact that each node must be either healthy, infected or failed. Let π i,0 , π i,P and π i,F be the steady state probabilities that node i is susceptible, infected and failed, respectively. Next, consider the special case of a connected graph where all nodes have the same degree k. In addition, assume κ i = κ for i = 1, . . . , N . Then, by symmetry we have π i,F = π F , π i,0 = π 0 and π i,P = π P for i = 1, . . . , N , where π 0 , π P and π F are the steady state probabilities that a node is susceptible, infected and failed, respectively, in the symmetric topology. It follows from ( 3)-( 6) that the nontrivial equilibrium for the symmetric topology is given by the value of π F which solves the following equation,

1 + r 1 λ 2 βkπ 2 F + r 1 + (κ -βk) r 1 λ + 1 π F = κ (7)
In the absence of exogenous infection rates (κ = 0) and in case infected nodes immediately degrade (λ = ∞), the solution of ( 7) is given by [6, equation ( 2)].

III. SAMPLING AND REJUVENATION

In this section we generalize the model presented in the previous section to account for sampling costs as well as rejuvenation. Our goal is to show that a simple adaptation of the previously presented model allows us to account for these additional features. To this aim, we extend [2, Figure 2] to general topologies.

We start presenting the basic balance equations which capture the dynamics presented in Figure 2. According to the NIMFA approximation [START_REF] Van Mieghem | Virus spread in networks[END_REF] we have

dπ i,P (t) dt = -(λ + r 4 )π i,P (t) + (8)
+ π i,0 (t)   N j=1 A ji β(1 -π j,0 (t)) + κ i   dπ i,R (t) dt = r 4 π i,P (t) -π i,R (t)(r 3 + γ) ( 9
)
dπ i,F (t) dt = λπ i,P (t) + γπ i,R (t) -r 1 π i,F (t) (10)
dπ i,0 (t) dt = - dπ i,P (t) dt + dπ i,R (t) dt + dπ i,F (t) dt ( 11
)
π i,0 (t) = 1 -π i,P (t) -π i,R (t) -π i,F (t) -π i,V (t) (12)
Comparing Figures 1 and2, the two new states represented in Figure 2 are R and V which correspond to nodes under rejuvenation that were sampled, wherein infection was detected at an early stage, and vaccinated nodes, respectively. We assume that sampling for early infections occurs according to a Poisson process with rate r 4 , and failures after

0 P R F r 3 r 1 λ γ r 4 r 2 V Figure 2.
A node of the complete model with rejuvenation. Edges representing vaccination, rejuvenation, full repair and aging are represented by dotted, broken-dotted, broken and solid lines, respectively infection occur with rate λ. Therefore, we have transitions from P to R and F with rates r 4 and λ, respectively. We have a transition from R to F with rate γ, characterizing failures that might occur during rejuvenation. The effect of γ is subject for future work, and in this paper we assume γ = 0 as in [START_REF] Huang | Software rejuvenation: Analysis, module and applications[END_REF]. Under normal conditions, a node that is prone to failure can be rejuvenated. Once started, the rejuvenation takes mean time 1/r 3 . Finally, full repairs occur from state F , and yield a transition to state 0.

Vaccinated nodes remain in state V , which is an absorbing state, as we assume that the vaccine completely imunizes the nodes and they are removed from the graph. In the remainder of this paper, except otherwise noted we focus on nodes that are not vaccinated.

Note that if we do not allow for sampling and rejuvenation (r 4 = r 3 = 0) the model introduced in this section degenerates to the one presented in Section II. The proposed model is general enough to allow for different recovery strategies at different nodes (i.e., some nodes might be amenable to rejuvenation and early sampling whereas others might count only with rejuvenation). The methodology to solve the model above is the same as the one used in the previous section.

A. Regular Topology

Next, we let κ i = 0 and consider the special case of a connected graph where all nodes have the same degree k. In addition, we restrict ourselves to nodes that are not vaccinated (recall that vaccinated nodes are removed from the infection network). Due to symmetry, we drop the subscript i from the state variables to obtain,

π 0 = λ + r 4 βk( λ r1 + 1 + r4 r3 ) (13)
π P = 1 -π 0 λ r1 + 1 + r4 r3 ( 14
) π R = r 4 r 3 π P , π F = λ r 1 π P ( 15
)
where we assume that βk > (λ+r 4 )/((λ/r 1 )+1+(r 4 /r 3 )).

Otherwise, only the trivial solution (π 0 , π P , π R , π F ) = (1, 0, 0, 0) exists.

B. Costs: Infection Probabilities and Sampling Rates

Next, we consider the infection and sampling costs incurred by the nodes. The setup is the one considered in Section III-A. The infection cost is the cost due to a machine having a virus. Let c v be the average cost per time unit incurred when a machine is infected. The cost c v is proportional to the severity of the infection. The average infection cost per time unit incurred by a node is proportional to c v and to the probability that a node is infected, 1 -π 0 , and is given by

C v (r 4 ) = c v (1 -π 0 ).
The sampling cost is the cost to examine the machines (e.g., to run a comprehensive checkup). Let c s be the sampling cost (per sample) to identify an infection at an early stage. Samples are taken at states 0 and P . The fraction of time spent at these states is (π 0 +π P ). Therefore, the expected number of samples taken per time unit is (π 0 + π P )r 4 and the expected sampling cost per time unit is C s (r 4 ) = (π 0 + π P )r 4 c s .

Let C(r 4 ) be the mean total cost per time unit. It is given by the sum of the mean sampling and infection costs,

C(r 4 ) = C s (r 4 ) + C v (r 4 ) = c s r 4 (π 0 + π P ) + c v (1 -π 0 ) (16
) Theorem 1: If sampling costs are negligible (c s = 0), the optimal sampling strategy consists of sampling the system at the highest possible rate (i.e., setting

r 4 = ∞) if 1 λ + 1 r 1 > 1 r 3 ( 17
)
and not sampling (i.e., setting r 4 = 0) otherwise.

Proof: Increasing the sampling rate is beneficial if the derivative of ( 13) with respect to r 4 is positive, i.e., dπ 0 /dr 4 ≥ 0. It is easily verified that dπ 0 /dr 4 ≥ 0 if ( 17) is satisfied.

Let r ⋆ 4 be the optimal sampling rate. In the presence of sampling costs (c s > 0), the following theorem establishes that either r ⋆ 4 = 0 or r ⋆ 4 is the unique local minimum of C(r 4 ) in (0, ∞).

Theorem 2: If c s > 0, the mean total cost C(r 4 ) admits at most one local minimum in (0, ∞). Proof: Replacing ( 13) and ( 14) into ( 16) we obtain an expression for C(r 4 ) that depends only on system parameters. Taking the second derivative of the resulting expression with respect to r 4 , we note that = c s r 3 /(kβ) is strictly positive, meaning that increasing r 4 beyond a threshold increases the mean total cost unboundedly. Therefore, i) C(r 4 ) cannot be strictly concave, ii) if C(r 4 ) is strictly convex it admits only one local minimum in [0, ∞), which is also the global minimum, and iii) if C(r 4 ) has a single inflection point it has at most two local minima, one at r 4 = 0 and another in (0, ∞). The result follows from i), ii) and iii).

IV. NUMERICAL RESULTS

In this section we numerically illustrate some of the conclusions derived in the previous sections. Our goals are a) to quantify the potential gains of rejuvenation, b) to show how different system parameters impact the mean total cost and the optimal sampling rate. Our base setup is the one considered in [2, Figure 3]: λ = 1/(12 × 30 × 24), r 1 = 2 and r 3 = 3. We consider a full mesh topology (complete graph) and let k = 10, β = 0.08, c s = 9 and c v = 10. Given the initial setup, we vary the system parameters according to our experimental goals.

Figure 3 shows the tradeoff involved in the choice of the sampling rate. Initially, when r 4 = 0 (no rejuvenation), the mean total cost equals infection costs (C(0) ≈ 10.00) and it is beneficial to increase the sampling rate. When r 4 = 0.60, the mean total cost reaches its minimum of C(0.60) = 8.8. Further increasing r 4 increases the mean total cost. When r 4 > 0.6 sampling costs dominate the total cost.

Figure 4 shows how the sampling cost per sample, c s , affects the optimal sampling rate. As c s increases, the optimal sampling rate decreases. The optimal sampling rate r 4 is roughly 0.85, 0.60, 0.35 for c s = 7, 9, 11 respectively.

Figure 5 shows how the infection rate per contact, β, affects the optimal sampling rate. As β increases, the optimal sampling rate decreases. The optimal sampling rate r 4 is roughly 0.7, 0.6, 0.5 for β = 0.07, 0.08, 0.09 respectively. This is because as nodes are more prone to infection, sampling must occur at a higher rate to cope with greater infection costs.

Figure 6 shows how the mean time to complete rejuvenation, 1/r 3 , affects the optimal sampling rate. As 1/r 3 decreases, the mean total cost decreases as the system repair time is reduced. In addition, the optimal sampling rate increases as r 3 increases. As the time to execute rejuvenation decreases, it becomes more advantageous to detect infections at very early stages and immediately treat the nodes as opposed to waiting for them to fail and then execute a full repair. The optimal sampling rate r 4 is roughly 0.35, 0.5, 0.6 for r 3 = 1, 2, 3 respectively.

Infection costs might differ between nodes, as more Figure 6. Varying the mean time to complete rejuvenation, 1/r 3 .

Figure 7. Asymmetric setup with c

(1)

v = 7, c (2)
v = 14.

critical nodes might incur higher costs than others. Consider an asymmetric setup where nodes are divided in two classes.

The population is equally divided between the two classes. Nodes of classes 1 and 2 incur costs c

(1) v = 7 and c

(2) v = 14, respectively, due to infection, when infected by viruses. The model presented in the previous section naturally extends to multiple classes. We use superscripts (1) and (2) to denote variables associated to classes 1 and 2. Figure 7 shows the mean total cost C(r

(1) 4 , r (2)
4 ) as a function of the sampling rates of nodes in the two classes. The minimum mean total cost is 8.9 and is attained when (r

(1) 4 , r (2)
4 ) = (0.6, 0.7). Nodes of class 2 need to be sampled at higher rates than those of class 1 because c

(2) v > c (1) v . In contrast, if c (1) v = c (

V. RELATED WORK

There is a vast literature on epidemic models to analyze the spread of diseases [START_REF] Newman | Networks: an introduction[END_REF], viruses in networks [START_REF] Wang | Epidemic spreading in real networks: An eigenvalue viewpoint[END_REF] and ideas in social networks [START_REF] Netrapalli | Learning the graph of epidemic cascades[END_REF]. Rejuvenation as a possible strategy to mitigate the effects of worm epidemics was introduced in [START_REF] Avritzer | Using performance signatures and software rejuvenation for worm mitigation in tactical manets[END_REF]. In this paper, we propose a novel model to account for the effects of rejuvenation in general topologies.

From an industrial viewpoint, the use of rejuvenation in telecommunication systems was discussed in [START_REF] Alonso | Software rejuvenation: Do IT and telco industries use it?[END_REF], [START_REF] Avritzer | Methods and opportunities for rejuvenation in aging distributed software systems[END_REF]. The typical approach consists of monitoring OS-level indicators to detect performance degradation and anomalies [START_REF] Bovenzi | An oslevel framework for anomaly detection in complex software systems[END_REF], [START_REF] Garg | Analysis of software rejuvenation using Markov regenerative stochastic Petri net[END_REF]. Avritzer and Weyuker [START_REF] Avritzer | Monitoring smoothly degrading systems for increased dependability[END_REF] propose heuristics to decide if a node requires treatment as a function of its performance degradation. In this paper, in contrast, we assume that the performance degradation of infected but non-failed nodes is negligible, and the performance degradation of failed nodes is immediately observable. Future work consists of considering multiple performance degradation levels and the effort to detect them.

We show that approximations to compute steady-state metrics on networks, such as the NIMFA model [START_REF] Van Mieghem | Virus spread in networks[END_REF], can be extended to account for rejuvenation and sampling. Our model extends the standard susceptible-infected (SI) model in a number of ways. In the standard SI model, nodes are either susceptible or infected, and transition among these two states until they all reach the susceptible state, which is an absorbing state. Previous works considered the time till absorption [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF] or the quasi-stationary regime [START_REF] Méléard | Quasi-stationary distributions and population processes[END_REF], [START_REF] Wang | Epidemic spreading in real networks: An eigenvalue viewpoint[END_REF], [START_REF] Van Mieghem | Virus spread in networks[END_REF], [START_REF] Hayel | Complete game-theoretic characterization of sis epidemics protection strategies[END_REF], [START_REF] Auletta | Metastability of logit dynamics for coordination games[END_REF]. As we show in Section II, some of the results presented in [START_REF] Van Mieghem | Virus spread in networks[END_REF], [START_REF] Hayel | Complete game-theoretic characterization of sis epidemics protection strategies[END_REF] are special cases of ours, obtained by setting the exogenous rate to zero.

Newman [START_REF] Newman | Networks: an introduction[END_REF]Section 18.3] studies dynamics on graphs, with more than one state variable per vertex. Our work extends [START_REF] Newman | Networks: an introduction[END_REF] by making use of the N-Interwined Mean-Field Approximation (NIMFA) [START_REF] Van Mieghem | Virus spread in networks[END_REF]. In addition, we also discuss sampling and rejuvenation strategies which are out of the scope of [START_REF] Newman | Networks: an introduction[END_REF].

Asavathiratham proposes the influence model [START_REF] Asavathiratham | The influence model: A tractable representation for the dynamics of networked Markov chains[END_REF] as a tractable way to handle pairwise interactions of nodes, each node having its associated Markov chain to represent its internal state. Garetto et al. [START_REF] Garetto | Modeling malware spreading dynamics[END_REF] study the spread of viruses in a network under this model. According to the influence model, a susceptible node that has one single neighbor, which turns out to be infected, will become infected at the same rate as a node with multiple neighbors all of which are infected. This is because interactions occur in a pairwise manner. In the model proposed in this paper, in contrast, the larger the number of infected neighbors incident to a susceptible node, the larger the infection rate of the susceptible node [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF], [START_REF] Van Mieghem | Virus spread in networks[END_REF].

VI. CONCLUSION

In light of security threats posed by viruses and worms to cyber infrastructures, the study of the spread of epidemics has gained significant attention in the past few years. In this paper, we discussed the role of sampling and rejuvenation as solutions to mitigate epidemics. We showed when and how rejuvenation is beneficial, and we established conditions under which the sampling rate admits one global minimum. Then, we numerically illustrated how different system parameters impact the optimal sampling rate. We believe that the models proposed in this paper are the first step towards a better understanding of how to sample and rejuvenate nodes accounting for security threats. This work opens up several avenues for future exploration, including control and game theoretical aspects involved in the choice of when and how to sample, vaccinate, rejuvenate and heal nodes.

Figure 1 .

 1 Figure 1. Simple model with 4 nodes without rejuvenation.

Figure 3 .

 3 Figure 3. Finding the optimal sampling rate

4 = 0

 40 d 2 C(r4) dr 2 admits at most one real root. It follows that C(r 4 ) is either i) strictly concave, ii) strictly convex or iii) has one single inflection point. But lim r4→∞ dC(r4) dr4

Figure 4 .

 4 Figure 4. Varying the sampling cost per sample, cs.

Figure 5 .

 5 Figure 5. Varying the infection rate per contact, β.

  2) v = 10 (not shown in the figure) the minimum cost is 8.8 and is attained when (r (0.6, 0.6) (in accordance to Figure 4, for c s = 9).

Table I TABLE

 I OF NOTATION. NOTE THAT WE DROP THE DEPENDENCY ON t WHEN CONSIDERING STEADY-STATE, AND THE DEPENDENCY ON i

	variable	description
	π i,0 (t)	probability that node i is healthy at time t
	π i,P (t)	probability that node i is prone to failure (infected) at time t
	π i,R (t)	probability that node i is under rejuvenation at time t
	π i,F (t)	probability that node i is failed at time t
	κ i	exogenous infection rate at node i
	r 1	full repair rate
	r i,2	endogenous infection rate at node i
	r 3	repair rate (rate of recovery once rejuvenation started)
	r 4	sampling rate
	λ	failure rate
	β	rate of infection of infected node towards each neighbor
	A ij	equals 1 if node j is a neighbor of i and 0 otherwise
	k	number of neighbors of each node in a regular graph
	N	number of nodes in the network
	cv	cost of infection per time unit
	cs	cost of sampling per sample
	C(r 4 )	mean total cost

WHEN CONSIDERING SYMMETRIC TOPOLOGIES.

Acknowledgment: Research conducted within the context of the THANES Associate Team, jointly supported by Inria (France) and FAPERJ (Brazil). This work has been partially supported by the European Commission within the framework of the CONGAS project FP7-ICT-2011-8-317672 and by Siemens Corporation, Corporate Technology.

PROOF OF THEOREM 1

The probability that a node is healthy is π 0 , given by [START_REF] Avritzer | Monitoring smoothly degrading systems for increased dependability[END_REF], π 0 = r 1 r 3 (λ + r 4 ) βk(λr 3 + r 1 r 3 + r 1 r 4 )

Therefore, dπ 0 dr 4 = r 1 r 3 (λr 3 -λr 1 + r 1 r 3 ) βk(λr 3 + r 1 r 3 + r 1 r 4 ) 2

As the denominator of ( 19) is always positive, dπ0 dr4 > 0 if

and in this case it is worth increasing r 4 unboundedly (r 4 = ∞). Otherwise, set r 4 = 0.

APPENDIX PROOF OF THEOREM 2

Follows below the Matlab code used to symbolically obtain the results discussed in the proof of theorem 2. % unfortunately, the derivative of C with % respect to r_4 is too complicate % so we cannot find its real roots directly % (there are 3 roots, that admit % closed form expressions, but they % are not simple, and it is not easy % to identify which ones are imaginary) simpleroots=simple(solve(... derivtotalcost==0,sr4)) pretty(simpleroots) % therefore, we compute the following limit % instead mylimit=limit(derivtotalcost,sr4,inf) % mylimit = c_s r_3 / (k beta) Note that the last limit is reported and used in the proof of Theorem 2.


V Baldoni
email: baldoni@mat.uniroma2.it

Nicole Berline
email: nicole.berline@math.polytechnique.fr

M K Öppe

Michèle Vergne
email: vergne@math.jussieu.fr

Matthias Köppe

INTERMEDIATE SUMS ON POLYHEDRA: COMPUTATION AND REAL EHRHART THEORY

Keywords: Primary); 52C07, 68R05, 68U05, 52B20 (Secondary)

We study intermediate sums, interpolating between integrals and discrete sums, which were introduced by A. Barvinok [Computing the Ehrhart quasi-polynomial of a rational simplex, Math. Comp. 75 (2006Comp. 75 ( ), 1449Comp. 75 ( -1466]]. For a given semi-rational polytope p and a rational subspace L, we integrate a given polynomial function h over all lattice slices of the polytope p parallel to the subspace L and sum up the integrals. We first develop an algorithmic theory of parametric intermediate generating functions. Then we study the Ehrhart theory of these intermediate sums, that is, the dependence of the result as a function of a dilation of the polytope. We provide an algorithm to compute the resulting Ehrhart quasi-polynomials in the form of explicit step polynomials. These formulas are naturally valid for real (not just integer) dilations and thus provide a direct approach to real Ehrhart theory.

Introduction

Let p be a rational polytope in V = R d and h(x) a polynomial function on V . A classical problem is to compute the sum of values of h(x) over the set of integral points of p,

S(p, h) = x∈p∩Z d h(x).
The sum S(p, h) has generalizations, the intermediate sums S L (p, h), where L ⊆ V is a rational vector subspace, introduced by Barvinok [START_REF]Computing the Ehrhart quasi-polynomial of a rational simplex[END_REF]. They interpolate between the discrete sum S(p, h) and the integral p h(x) dx. For a polytope p ⊂ V and a polynomial h(x)

S L (p, h) = x p∩(x+L)
h(y) dy, where the summation index x runs over the projected lattice in V /L. In other words, the polytope p is sliced along affine subspaces parallel to L through lattice points and the integrals of h over the slices are added up. For L = V , there is only one term and S V (p, h) is just the integral of h(y) over p, while, for L = {0}, we recover the discrete sum S(p, h). As in the discrete case, a powerful method is to consider the intermediate generating function

S L (p)(ξ) = x p∩(x+L)
e ξ,y dy,

for ξ ∈ V * . If p is a polyhedron, not necessarily compact, the generating function still makes sense as a meromorphic function. By writing a polyhedron as the sum of its cones at vertices (Brion's theorem), we need only study the case where p is an affine cone.

It is then natural to turn to the Ehrhart theory of the intermediate sums S L , that is, the study of the intermediate sum S L (tp, h) of a dilated polytope tp as a function of the dilation parameter t. It turns out that, just like in the classical case, the Ehrhart function t → S L (tp, h) is a quasi-polynomial, that is, a function of the form

S L (tp, h) = d+M m=0 E L m (p, h; t) t m , (2)
where the coefficients E L m (p, h; t) depend only on t mod q, where q is an integer such that qp has lattice vertices.

The main results of this article are:

1. a polynomial time algorithm for the computation of the intermediate generating function of a simplicial affine cone, when the slicing space has fixed codimension, Theorem 24, 2. a polynomial time algorithm for the computation of the weighted intermediate sum S L (p, h) of a simple polytope p (given by its vertices), and the corresponding Ehrhart quasipolynomial t → S L (tp, h), when the slicing space has fixed codimension and the weight depends only on a fixed number of variables, or has fixed degree, Theorem 28.

The main new feature of this article is the use of a decomposition of any simplicial cone into simplicial cones which have a face parallel to the given space L, the decomposition being given by a closed formula, Theorem 18. This decomposition is borrowed from [START_REF] Brion | Residue formulae, vector partition functions and lattice points in rational polytopes[END_REF] where it plays a key role in the study of partition functions. We prove that this decomposition is done by a polynomial algorithm, if the codimension of L is fixed. Thus, the computation of the intermediate generating function of a simplicial affine cone is reduced to the case where L is parallel to a face of the cone, and thus the generating function factors into a discrete sum and an integral. This case has been already studied in [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF]. The output of the algorithm is a "short formula" where both ξ and the vertex s of the cone appear as symbolic variables (see Examples 6,7 and 26). Once a short formula for the generating function is available, we follow the method of [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF] for the computation of the weighted Ehrhart quasipolynomial of a simple polytope. The method uses the fact that the intermediate generating function decomposes into meromorphic terms of homogeneous ξ-degree. The output is again a short formula where the dilating parameter t appears as a symbolic variable inside step polynomials, see Example 30.

The results remains valid for non-negative real values of t. Thus our method provides a direct proof of a real Ehrhart theorem (Theorem 31). This extension of Ehrhart theory to real dilation factors has been studied by E. Linke [START_REF] Linke | Rational Ehrhart quasi-polynomials[END_REF] for the classical case. It is even more natural for the intermediate valuations, as one of the cases is L = V . Then S L (tp, h) is the integral of h over the tp, the result being clearly a polynomial formula valid for any non-negative real value of t.

The present manuscript together with the article [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF] also provide the foundation for a future work, in which the weighted Barvinok patched sums L λ(L)S L (p, h) for a simple polytope will be studied.

Notations and basic facts

2.1. Rational and semi-rational polyhedra. We consider a rational vector space V of dimension d, that is to say a finite dimensional real vector space with a lattice denoted by Λ. We will need to consider subspaces and quotient spaces of V , this is why we cannot just let

V = R d and Λ = Z d . An element v ∈ V is called rational if nv ∈ Λ for some non-zero integer n. The set of rational points in V is denoted by V Q . A subspace L of V is called rational if L ∩ Λ is a lattice in L. If
L is a rational subspace, the image of Λ in V /L is a lattice in V /L, so that V /L is a rational vector space. The image of Λ in V /L is called the projected lattice. It is denoted by Λ V /L . A rational space V , with lattice Λ, has a canonical Lebesgue measure dx = dm Λ (x), for which V /Λ has measure 1.

A convex polyhedron p in V (we will simply say polyhedron) is, by definition, the intersection of a finite number of closed half spaces bounded by affine hyperplanes. If the hyperplanes are rational, we say that the polyhedron is rational. If the hyperplanes have rational directions, we say that the polyhedron is semi-rational. For instance, if p ⊂ V is a rational polyhedron, t is a real number and s is any point in V , then the dilated polyhedron tp and the translated polyhedron s+p are semirational. Unless stated otherwise, a polyhedron will be assumed to be semi-rational.

In this article, a cone is a polyhedral cone (with vertex 0) and an affine cone is a translated set s + c of a cone c. A cone c is called simplicial if it is generated by independent elements of V . A simplicial cone c is called unimodular if it is generated by independent integral vectors v 1 , . . . , v k such that {v 1 , . . . , v k } can be completed to an integral basis of V . An affine cone a is called simplicial (resp. simplicial unimodular ) if the associated cone is.

A polytope p is a compact polyhedron. The set of vertices of p is denoted by V(p). For each vertex s, the cone of feasible directions at s is denoted by c s .

We denote by L(V ) the space of functions on V which is generated by indicator functions [q] of polyhedra q which contain lines.

Valuations and generating functions.

Definition 1. Let M be a vector space. An M -valued valuation Φ is a map from the set of polyhedra p ⊂ V to the vector space M such that whenever the indicator functions [p i ] of a family of polyhedra p i satisfy a linear relation i r i [p i ] = 0, then the elements Φ(p i ) satisfy the same

relation i r i Φ(p i ) = 0.
Definition 2. We denote by H(V * ) the ring of holomorphic functions defined around 0 ∈ V * . We denote by M(V * ) the ring of meromorphic functions defined around 0 ∈ V * and by M (V * ) ⊂ M(V * ) the subring consisting of meromorphic functions φ(ξ) which can be written as a quotient of a holomorphic function by a product of linear forms.

The intermediate generating functions of polyhedra which are studied here are important examples of functions in M (V * ). Proposition 3. Let L ⊆ V be a rational subspace. There exists a unique valuation S L which to every semi-rational polyhedron p ⊂ V associates a meromorphic function S L (p) ∈ M(V * ) so that the following properties hold:

(a) If p contains a line, then S L (p) = 0. (b) S L (p)(ξ) = x∈Λ V /L p∩(x+L) e ξ,y dm L (y), (3)
for every ξ ∈ V * such that the above sum converges. (c) For every point s ∈ Λ + L, we have S L (s + p)(ξ) = e ξ,s S L (p)(ξ).

When necessary we will indicate the lattice and use the notation S L (p, Λ).

Remark 4. The fact that S L (p) is actually an element of M (V * ) will be a consequence of the explicit computations of the next section.

For L = {0}, we recover the valuation S given by

S(p)(ξ) = x∈p∩Λ e ξ,s , provided this sum is convergent. For L = V , S V (p)(ξ) is the integral I(p)(ξ) = p e ξ,x
dx, if p is full dimensional, and S V (p) = 0 otherwise. If p is compact, the meromorphic function S L (p)(ξ) is actually regular at ξ = 0, and its value for ξ = 0 is the valuation E L ⊥ (p) considered by Barvinok [6]. The proof is entirely analogous to the cases L = {0} and L = V , see for instance Theorem 3.1 in [START_REF] Barvinok | An algorithmic theory of lattice points in polyhedra[END_REF] and we omit it.

Remark 5 (Rational vs. semi-rational polyhedra). Although [START_REF] Barvinok | An algorithmic theory of lattice points in polyhedra[END_REF] and other classical texts deal only with rational polyhedra, their proofs extend easily to semi-rational ones. For instance, by summing geometric series, one obtains immediately the generating function of a unimodular cone with a real vertex.

Here are some simple examples. For t ∈ R we denote by {t} the fractional part of t, i.e., 0 ≤ {t} < 1 and t -{t} ∈ Z, and by t the ceiling of t, i.e., t = t + {-t}.

Example 6. For any t ∈ R, we have

S(t + R ≥0 )(ξ) = e t ξ
1 -e ξ = e tξ e {-t}ξ 1 -e ξ .

More generally, let c ⊂ R d be a unimodular cone with primitive edge generators v 1 , . . . , v d . Let s be any point in R d , with s = d i=1 s i v i . Then

S(s + c)(ξ) = d i=1 e s i ξ,v i 1 -e ξ,v i = e s,ξ d i=1 e {-s i } ξ,v i 1 -e ξ,v i . Example 7. Let c ⊆ R 2 be the first quadrant and L = Re 1 . We write s = (a, b) ∈ R 2 .
We have

S L (s + c)(ξ 1 , ξ 2 ) = ∞ n= b +∞ a e nξ 2 +xξ 1 dx = ∞ n= b e nξ 2 • -e aξ 1 ξ 1 = e b ξ 2 1 -e ξ 2 • -e aξ 1 ξ 1 = e s,ξ e {-b}ξ 2 1 -e ξ 2 • -1 ξ 1 .
Example 8. Let p ⊂ R 2 be the triangle with vertices (0, 0), (1, 0), (0, 1).

Let L = R(1, 0). Straightforward calculations give S L (p) = 1-e ξ 1 ξ 1 , I(p)(ξ) = 1 ξ 1 -ξ 2 1-e ξ 1 ξ 1 -1-e ξ 2 ξ 2
, and S(p)(ξ) = 1 + e ξ 1 + e ξ 2 . These three functions are indeed analytic.

A consequence of the valuation property is the following fundamental theorem. It follows from the Brion-Lawrence-Varchenko decomposition of a polytope into the supporting cones at its vertices [START_REF] Brion | Points entiers dans les polyèdres convexes[END_REF].

Theorem 9 (Brion's theorem). Let p be a polytope with set of vertices V(p). For each vertex s, let c s be the cone of feasible directions at s. Let L ⊆ V be a rational subspace. Then Where the slicing space is parallel to a face. Let c ⊂ V be a simplicial cone with integral generators (v j , j = 1, . . . , d), and let s ∈ V . Let us recall the expression for the intermediate generating function S L (s+c)(ξ), when the slicing subspace L is generated by some of the v j 's. For I ⊆ {1, . . . , d}, we denote by L I the linear span of the vectors (v i , i ∈ I). We denote the complement of I in {1, . . . , d} by I c . We have V = L I ⊕ L I c . For x ∈ V we denote the components by

S L (p) = s∈V(p) S L (s + c s ).
x = x I + x I c .
Thus we identify the quotient V /L I with L I c and we denote the projected lattice by Λ

I c ⊂ L I c . Note that L I c ∩ Λ ⊆ Λ I c , but the inclusion is strict in general. Example 10. Let v 1 = (1, 0), v 2 = (1, 2), I = {2}. The projected lattice Λ I c on L I c = Rv 1 is Z v 1 2 . See Figure 1.
We denote by c I the cone generated by the vectors v j , for j ∈ I and by b I the parallelepiped b I = i∈I [0, 1[ v i . Similarly we denote by c I c the cone generated by the vectors v j , for j ∈ I c and b

I c = i∈I c [0, 1[ v i .
The projection of the cone c on V /L I = L I c identifies with c I c . Note that a generator v i , i ∈ I c , may be non-primitive for the projected lattice Λ I c , even if it is primitive for Λ, as we see in the previous example. We write s = s I + s I c .

First let us recall from [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF] how the intermediate generating function S L I (s + c, Λ) decomposes as a product. This already appeared in Example 7.

Proposition 11. [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF]. The intermediate sum for the full cone s + c breaks up into the product

S L I (s + c, Λ)(ξ) = S(s I c + c I c , Λ I c )(ξ) I(s I + c I , L I ∩ Λ)(ξ). ( 4
)
We remark that in the above equation, the function S(s

I c +c I c , Λ I c )(ξ) is a meromorphic function on the space (L I c ) * . The integral I(s I + c I , L I ∩ Λ)(ξ)
is a meromorphic function on the space (L I ) * . We consider both as functions on V * through the decomposition V = L I ⊕L I c .

Cone decomposition with respect to a linear subspace.

Let L be a linear subspace of V . Let c be a full dimensional simplicial cone in V with generators w 1 , . . . , w d . The main theorem of this section, Theorem 18, gives a signed decomposition of c into a family of full dimensional simplicial cones c σ , where each c σ has a face parallel to L. More precisely, this decomposition is a relation between indicator functions of cones, modulo the space L(V ) of functions on V which is generated by indicator functions [q] of polyhedra q which contain lines. This decomposition is borrowed from [START_REF] Brion | Residue formulae, vector partition functions and lattice points in rational polytopes[END_REF]. The proof relies on the following simple lemma which has its own interest. The extreme cases, dim L = 1 or codim L = 1, of Theorem 18 follow easily from this lemma (Proposition 14). The general case is proven by induction on either dim L or codim L.

Given vectors u 1 , . . . , u p ∈ R d , we denote the cone p i=1 R ≥0 u i by c(u 1 , . . . , u p ) and its relative interior p i=1 R >0 u i by c(u 1 , . . . , u p ). Lemma 12. Let u 1 , . . . , u d be a basis of V and let

u d+1 = -d k=1 u k . For i = 1, . . . , d + 1, define semi-open cones k i = c(u 1 , . . . , u i-1 ) + c(u i+1 , . . . , u d+1 ).
Then V is the disjoint union of the k i for i = 1, . . . , d + 1.

Proof. Let x = t 1 u 1 + . . . + t d u d . If t k ≥ 0 for all k, then x ∈ k d+1 . Otherwise, let i be the largest index with 1 ≤ i ≤ d, for which the minimum min 1≤k≤d t k is reached. Writing u i = -1≤k≤d+1,k =i u k , we obtain x = 1≤k≤d,k =i (t k -t i )u k -t i u d+1 . As t k -t i ≥ 0 if 1 ≤ j ≤ d and t j -t i > 0 if i < j ≤ d, and t i < 0, we see that x ∈ k i . Let us show that k i ∩ k j is empty if i < j. Let x ∈ k i , x = 1≤k≤d+1,k =i t k u k . Writing u j = -1≤k≤d+1,k =j u k , we obtain x = 1≤k≤d, k / ∈{i,j} (t k -t j )u k -t j u i .
As t j > 0, it follows that x / ∈ k i .

From the lemma, we deduce other interesting cone decompositions. The first one is the well-known stellar decomposition.

Proposition 13 (Stellar decomposition). Let w 1 , . . . , w d be a basis of V . Let r ≤ s and

v = w 1 + • • • + w r -(w s+1 + • • • + w d ).
Let c be the cone generated by w i 's for 1 ≤ i ≤ d and let k be the cone generated by c and the vector v.

k i = c(w 1 , . . . , w i-1 ) + c(w i+1 , . . . , w d , v) for 1 ≤ i ≤ r, if r ≥ 1, k i = c(w 1 , . . . , w i-1 ) + c(w i+1 , . . . , w d , v) for s + 1 ≤ i ≤ d, if s < d. Then: (a.1) If r ≥ 1, k is the disjoint union of the k i for 1 ≤ i ≤ r. (a.2) k \ c is the disjoint union of the k i for s + 1 ≤ i ≤ d. (If s = d, then k \ c = ∅.) (b) For i ∈ {1, . . . , r} ∪ {s + 1, . . . , d}, let c i = c(w 1 , . . . , w i-1 , w i+1 , . . . , w d , v). Then [c] - r i=1 [c i ] + d i=s+1 [c i ] ≡ 0 modulo indicator functions of lower dimensional cones.
Proof. Let us prove (a.2) first. By quotienting out by the subspace generated by w 1 , . . . , w s , we can assume that s = 0, thus v = -(w 1 + . . . + w d ). In this case, we just apply Lemma 12 with u k = w k and u d+1 = v. As the cone k d+1 of Lemma 12 is just c, we obtain (a.2).

Let us now assume that l ≥ 1. We write

w 1 = v + w d + • • • + w s+1 -(w r + • • • + w 2 ).
We apply (a. Next, we obtain decompositions of any simplicial cone into fulldimensional cones which have a given edge, or a given facet. Recall that L(V ) is the linear span of indicator functions of cones which contain straight lines. Proposition 14. Let w 1 , . . . , w d be a basis of V and let c be the cone generated by w i 's for

1 ≤ i ≤ d. Let 0 ≤ r ≤ s ≤ d + 1. a) Let v = w 1 + • • • + w r -(w s+1 + • • • + w d ).
(5) Then we have the following relation between indicator functions of cones.

[c] ≡ r i=1 (-1) i-1 [c(-w 1 , . . . , -w i-1 , w i+1 , . . . , w d , v)] + d i=s+1 (-1) d-i [c(w 1 , . . . , w i-1 , -w i+1 , . . . , -w d , -v)] mod L(V ). (6)
b) Let L ⊂ V be a hyperplane. Assume that w i ∈ L if and only if r + 1 ≤ i ≤ s, that the w i 's lie on one side of L for 1 ≤ i ≤ r and on the other side for s + 1 ≤ i ≤ d. For i ∈ {1, . . . , r} ∪ {s + 1, . . . , d}, denote by ρ i the projection V → L parallel to w i . Then we have the following relation between indicator functions of cones.

[c] ≡ r i=1 [R + w i + ρ i (c)] - d i=s+1 [R + (-w i ) + ρ i (c)] mod L(V ). (7) Example 15. c = c(w 1 , w 2 ), v = w 1 + w 2 . Then [c] = -[c(-w 2 , v)] + [c(w 1 , v)] + [c(w 2 , -w 2 , w 1 )] . Indeed [c(w 2 , -w 2 , w 1 )]-[c] is the indicator function of the quadrant c(w 1 , -w 2 ) minus that of the line Rw 1 . This is also the case for [c(-w 2 , v)] + [c(w 1 , v)].
See also Figure 2.

Example 16. We consider the 3-dimensional cone c = c((-1, 0, 0), (-1, 2, 0), (-1, 0, 3)). i) First let L be the subspace with basis (0, 1, 1). Formula (6) gives

[c] ≡ [c((-1, 0, 0), (-1, 2, 0), (0, 1, 1))] -[c((-1, 0, 0), (1, 0, -3), (0, 1, 1))] + [c((-1, 2, 0), (-1, 0, 3), (0, -1, -1))] .

ii) Now, let L be the subspace with basis (1, 0, 0), (0, 1, 1). Formula (7) gives Proof of Proposition 14. Let us prove a). By quotienting out by the subspace generated by w r+1 , . . . , w s , we can assume that r = s. We write [START_REF] Barvinok | Polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed[END_REF] as

w r + • • • + w 1 -v -w d -• • • -w r+1 = 0.
By applying Lemma 12 with u 1 = -w r , . . . , u r+1 = -v, . . . , u d+1 = w r+1 , we see that the whole space V is the disjoint union of the semiopen cones

k 1 = c(w r-1 , . . . , -v, -w d , . . . , -w r+1 ), . . . k r = c(w r , . . . , w 2 ) + c(-v, -w d , . . . , -w r+1 ), k r+1 = c(w r , . . . , w 1 ) + c(-w d , . . . , -w r+1 ), k r+2 = c(w r , . . . , w 1 , -v) + c(-w d-1 , . . . , -w r+1 ), . . . k d+1 = c(w r , . . . , w 1 , -v, -w d , . . . , -w r+2 ).
For any basis u 1 , . . . , u d , we have

[(c(u 1 , . . . , u k ) + c(u k+1 , . . . , u d ))] ≡ (-1) d-k [c(u 1 , . . . , u k , -u k+1 , . . . , -u d )] mod L(V ).
Therefore we obtain

[k 1 ] ≡ (-1) d [c(-w r-1 , . . . , -w 1 , v, w d , . . . , w r+1 )] . . . [k r ] ≡ (-1) d-r [c(w r , . . . , w 2 , v, w d , . . . , w r+1 )] [k r+1 ] ≡ (-1) d-r-1 [c(w r , . . . , w 1 , w d , . . . , w r+1 )] . . . [k d+1 ] = [c(w r , . . . , w 1 , -v, -w d , . . . , -w r+2 )] .

Moreover these indicator functions add up to [V ] ≡ 0 mod L(V ). Thus we have proven a).

Let us show that b) is just the dual of a stellar decomposition. Let α ∈ V * be such that ker α = L, and such that α(w i ) > 0 for 1 ≤ i ≤ r and α(w i ) < 0 for s + 1 ≤ i ≤ d. Denote the dual cone by c

• ⊂ V * , c • = { ξ ∈ V * : ξ, x ≥ 0 for all x ∈ c }.
Let (γ 1 , . . . γ d ) be the dual basis of (v 1 , . . . , v d ). Thus c • = c(γ 1 , . . . , γ d ).

From the stellar decomposition of the cone c • with respect to the extra edge α, we see that the expression

[c • ] - i∈{1,...,r}∪{s+1,...,d} [c(γ 1 , . . . , γi , . . . , γ d , α)]
is a linear combination of indicator functions of lower dimensional cones. If a linear identity holds for indicator functions of cones, the same identity holds for the indicator functions of their duals, (see for instance [START_REF] Barvinok | An algorithmic theory of lattice points in polyhedra[END_REF], Corollary 2.8). Moreover, the dual of a lower dimensional cone contains a line, therefore by taking indicator functions of dual cones we obtain an equality mod L(V ). (This "duality trick" goes back to [START_REF] Brion | Points entiers dans les polyèdres convexes[END_REF]). There remains to compute the dual cone (c(γ 1 , . . . , γi , . . . , γ d , α)) • . We see easily that if α, w i > 0, then (c(γ 1 , . . . , γi , . . . , γ d , α))

• = R + w i + ρ i (c), while if α, w i < 0, then (c(γ 1 , . . . , γi , . . . , γ d , α)) • = R + (-w i ) + ρ i (c).

Thus we have proven b).

As we will see now, Proposition 14 consists of particular cases of the next result, Brion-Vergne decomposition. Let L be a linear subspace of V . Let c be a simplicial cone in V with generators w 1 , . . . , w d . We consider the subsets σ ⊂ {1, . . . , d} such that the vectors w j , for j ∈ σ, form a supplementary basis to L in V . In other words we have a direct sum decomposition

V = j∈σ Rw j ⊕ L.
The corresponding projection on L is denoted by ρ σ : V → L. The set of these σ is denoted by B(c, L).

Remark 17. At most d k 0

= O(d k 0 ) such sets σ exist. Thus, if codim L is a fixed constant k 0 , this is a polynomial quantity, and the sets σ ∈ B(c, L) can be enumerated by a straightforward algorithm in polynomial time. (In practice, reverse search [START_REF] Avis | Reverse search for enumeration[END_REF] would be used for this enumeration.)

A vector a in V /L is said to be generic with respect to the pair (c, L) if a does not belong to any hyperplane generated by vectors among the projections of w 1 , . . . , w d . In other words, for any σ ∈ B(c, L), we have a = j∈σ a σ,j (w j mod L), where a σ,j = 0 for every j ∈ σ.

a c L w 2 mod L w 1 mod L V /L v w 2 w 1 ≡ a L w 1 mod L σ = {1} c σ w 2 ρ {1} (w 2 ) {1},1 w 1 - a L c σ σ = {2} w 2 mod L {2},2 w 2 w 2 ρ {2} (w 1 )
Figure 2. Brion-Vergne decomposition of a cone c into cones with a face parallel to the subspace L, modulo cones with lines. The vectors in the quotient V /L determine the signs σ,j .

Theorem 18 (Brion-Vergne decomposition, [START_REF] Brion | Residue formulae, vector partition functions and lattice points in rational polytopes[END_REF], Theorem 1.2). Let L be a linear subspace of V . Let c be a full dimensional simplicial cone in V with generators w 1 , . . . , w d . Fix a vector a ∈ V /L that is generic with respect to the pair (c, L) and belongs to the projection of c on V /L. For σ ∈ B(c, L), let a = j∈σ a σ,j (w j mod L).

Let σ,j be the sign of a σ,j and (σ) = j∈σ σ,j . Denote by c σ ⊂ V the cone with edge generators σ,j w j for j ∈ σ, and ρ σ (w k ) for k / ∈ σ.

Then we have the following relation between indicator functions of cones.

[c] ≡

σ∈B(c,L) (σ) [c σ ] mod L(V ). ( 8
)
The theorem is illustrated in Figure 2.

Remark 19. As explained in [START_REF] Brion | Residue formulae, vector partition functions and lattice points in rational polytopes[END_REF], the case where dim L = 1 reduces to case a) of Proposition 14. Let us show that the case where codim L = 1 reduces to case b). A vector a in V /L is generic with respect to (c, L) if and only if a = 0. Up to renumbering, we can assume that w 1 , . . . , w r lie on the same side of L as a, while w s+1 , . . . , w d lie on the other side.

The set B(c, L) consists of the singletons σ i = {i} such that w i / ∈ L, ie i ∈ {1, . . . , r} ∪ {s + 1, . . . , d}. The sign σ i ,i is +1 if i ∈ {1, . . . , r} and -1 if i ∈ {s + 1, . . . , d}, therefore Equation (8) becomes precisely Equation [START_REF] Barvinok | An algorithmic theory of lattice points in polyhedra[END_REF].

Theorem 20 (Efficient construction of a Brion-Vergne decomposition). Fix a number k 0 . There exists a polynomial-time algorithm for the following problem. Given a rational simplicial cone c and a rational subspace L of V of codimension at most k 0 , compute a decomposition [START_REF] Brion | Points entiers dans les polyèdres convexes[END_REF].

Proof. By Remark 17, we can enumerate the bases σ ∈ B(c, L) in polynomial time.

A rational generic vector a can be constructed in polynomial time as follows; this technique was already used in [START_REF] Barvinok | Polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed[END_REF]. Fix one of the bases σ ∈ B(c, L) as a basis of V /L. Then the projections of w 1 , . . . , w d have rational coordinate vectors w1 , . . . , wd ∈ Q k 0 with respect to this basis.

Consider the moment curve ā(t) = (1, t, t 2 , . . . , t k 0 -1 ) ∈ R k 0 + for t ≥ 0. For any σ ∈ B(c, L) and j ∈ σ, we have by Cramer's rule that the coefficient a σ,j (t) for a = a(t) is nonzero if and only if f σ,j (t) = det(a(t); wi : i ∈ σ, i = j) = 0. The functions f σ,j (t) are polynomial functions in t of degree at most k 0 (k 0 -1), which are not identically zero. Hence each has at most k 0 (k 0 -1) zeros. Thus, by computing the coefficients a σ,j (t) for M = k 0 (k 0 -1)|B(c, L)| + 1 different values for t, we will find one t = t that is not a zero of any f σ,j . Clearly this search can be implemented in polynomial time.

Then the vector a = i∈σ āi (t)(w i mod L) ∈ V /L is generic with respect to the pair (c, L), and since āi (t) ≥ 0, it also lies in the projection of c on V /L. (We could also use explicit deterministic constructions of generic vectors of polynomial encoding size, without search, similar to [START_REF] Köppe | A primal Barvinok algorithm based on irrational decompositions[END_REF][START_REF] Köppe | Computing parametric rational generating functions with a primal Barvinok algorithm[END_REF].)

Remark 21 (Relation to variants of Barvinok's decomposition). Barvinok's method to efficiently compute the discrete generating function S(s + c)(ξ) of a cone in fixed dimension decomposes a simplicial cone into cones with smaller indices (and ultimately into unimodular cones). The original method by Barvinok [5] used a primal decomposition and inclusion-exclusion to take care of lower-dimensional cones. The dual of the stellar decomposition is used in the variant of Barvinok's method that was presented in [START_REF] Barvinok | An algorithmic theory of lattice points in polyhedra[END_REF] and later implemented in LattE [START_REF] De Loera | LattE, version 1.2[END_REF]; this avoids having to take care of lower dimensional cones. Instead of dualizing, one can directly uses Formula (6), the case where dim L = 1. We use this primal variant of Barvinok's decomposition in our implementation. Other primal variants were studied in [START_REF] Köppe | A primal Barvinok algorithm based on irrational decompositions[END_REF] and [START_REF] Köppe | Computing parametric rational generating functions with a primal Barvinok algorithm[END_REF] and implemented in LattE macchiato.

Complexity of the intermediate generating function.

We will use the following notations.

Definition 22. Let T (t, x) = e tx x 1 -e x = - ∞ n=0 B n (t) x n n! , (9)
where B n (t) are the Bernoulli polynomials.

Definition 23. For a and q a positive integer, we denote n q := q 1 q n ∈ qZ, {n} q := (n mod q) ∈ [0, q), which give the unique decomposition

n = n q + {n} q .
Thus a 1 = a and {a} 1 = {a} will be the ordinary "floor" and "fractional part" of a.

Theorem 24 (Short formula for S L (s + c)(ξ)). Fix a non-negative integer k 0 . There exists a polynomial time algorithm for the following problem. Given the following input: (O 1 ) a finite set Γ, (O 2 ) for every γ in Γ, integers α (γ) and lattice vectors w

(γ) i ∈ Z d for i = 1, . . . , d,
which have the following properties.

(1) For every γ, the family (w (γ) i ), for i = 1, . . . , d, is a basis of R d and if we denote the dual basis by (η

(γ) i ), for i = 1, . . . , d, then η (γ) i ∈ Z d for i = 1, . . . , k 0 .
(2) For every s ∈ R d , we have the following equality of meromorphic functions of ξ:

S L (s + c, Z d )(ξ) = e ξ,s γ∈Γ α (γ) k 0 i=1 T -η (γ) i , s , ξ, w (γ) i
•

1 d i=1 ξ, w (γ) i
.

Remark 25. Consider the term corresponding to γ ∈ Γ in [START_REF] De Loera | LattE, version 1.2[END_REF]. Define

s (γ) = k 0 i=1 -η (γ) i , s w (γ)
i .

The point s + s (γ) is in k 0 i=1 Zw (γ) i ⊕ d i=k 0 +1 Rw (γ)
i , and formula (10) reads also

S L (s + c, Z d )(ξ) = γ∈Γ α (γ) e ξ,s+s (γ) 1 k 0 i=1 1 -e ξ,w (γ) i 1 d i=k 0 +1 ξ, w (γ) i
.

Formula ( 11) is easier to grasp, while [START_REF] De Loera | LattE, version 1.2[END_REF] is used in the program, where we avoid all the singular hyperplanes ξ, w

= 0, for γ ∈ Γ and i = 1, . . . , d, by deforming the linear form ξ into ξ + .

Proof. Using Theorem 18, we construct a signed decomposition of c into cones c σ which have a face parallel to L, so that for any s ∈ R d we have

[s + c] ≡ σ σ [s + c σ ] (modulo cones containing lines), (12)
where σ ∈ {±1}. Here σ runs through the set B(c, L) of Theorem 18, which can be enumerated in polynomial time by Remark 17. This makes the construction of the signed decomposition a polynomial-time algorithm. As the valuation S L vanishes on cones which contain lines, we get

S L (s + c)(ξ) = σ∈B(c,L) σ S L (s + c σ )(ξ).
Each cone c σ is given by its edge generators w σ i , i = 1, . . . , d, which are such that (w σ i ), for i = k 0 + 1, . . . , d), is a basis of L. Then we apply Theorem 28 [Short formula for S L I (s + c, Z d )(ξ)] of [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF] with input c σ and I = {k 0 + 1, . . . , d}. See Examples 6 and 7.

Example 26 (Continuation of Example 16). Again c is the cone c((0, -2, 0), (1, -2, 0), (0, -2, 3)).

Here is the output of our Maple program computing the generating functions S L (s + c)(ξ) (Formula (11)). We set s = (a, b, c) and ξ = (x, y, z). i) Let L be the subspace with basis (0, 1, 1). In this case, the three cones which occur in the decomposition of Example 16 have non-unimodular projections in V /L, so that the projections are further decomposed into unimodular cones. This is why we have six terms in the expression of S L (s + c):

S L (s + c)(x,
- e -{3 a-b+c}z+{-a}(x-3z) (1 -e -z )(1 -e x-3 z )(-y -z) + e -{-b+c}z+{-a}x (1 -e -z )(1 -e x )(-y -z) + e {-b+c}y-{a}x (1 -e y )(1 -e -x )(-y -z) + e {a}(-x+2y)-{2 a+b-c}y (1 -e -x+2 y )(1 -e -y )(-y -z) .
ii) Let L be the subspace with basis (1, 0, 0), (0, 1, 1). Then:

S L (s + c)(x, y, z) = e ax+by+cz -6 e {-b+c}(-1 2 x+y)
1 -e - Example 27. We consider the four triangles T i of Figure 3 which subdivide the square q = [0, 4] × [0, 4]. Let L = R(0, 1). Let us apply the inclusion-exclusion formula together with the valuation property of S L . As S L vanishes on the edges and vertices which occur in the intersections of the triangles, it follows that the sum 4 i=1 S L (T i ) must be equal to S L (q). The example has been chosen so that the intermediate valuation S L (q) can also be computed directly by elementary means. This allows us to verify our methods on the example. A direct computation gives S L (q) = e 4y -1 y (1 + e x + e 2x + e 3x + e 4x ). Table 1 shows the output of our Maple program that implements the algorithm of the theorem for S L (T i ), i = 1, . . . , 4. The equality 4 i=1 S L (T i ) = S L (q) is readily checked on these expressions. (1-e x+2/15 y )y q ( 0 0 ) , ( 0 4 ) , ( 4 4 ) , (

e 4y -1 y (1 + e x + e 2x + e 3x + e 4x )

Ehrhart quasipolynomials for intermediate valuations

Let us fix a rational polytope p and a weight h(x). It is well known ("Ehrhart's Theorem") that when we dilate p by an integer t ∈ N, the function t → S(tp, h) is a quasi-polynomial function of t of degree d + M where M is the degree of h,

S(tp, h) = d+M m=0 E m (p, h; t) t m . ( 13
)
The coefficients E m (p, h; t) depend only on {t} q where q is an integer such that qp has lattice vertices. When p is a lattice polytope, then the coefficients E m (p, h; t) do not depend on t, and t → S(tp, h) becomes a polynomial. E. Linke [START_REF] Linke | Rational Ehrhart quasi-polynomials[END_REF] has proved that the expression ( 13) is still valid for real dilations t ∈ R. This is the generalized setting that we are going to use in this section. Note that if we allow real dilations, then even for lattice polytopes p we obtain a quasi-polynomial rather than a polynomial; the fractional part function {•} q with q = 1 appears in the expressions for the coefficients E m (p, h; t). For instance, the number of integers in the interval [0

, t] = t[0, 1] is t + 1 -{t} 1 .
We are going to extend Ehrhart's Theorem to the intermediate valuation S L (p). Indeed, one still has

S L (tp, h) = d+M m=0 E L m (p, h; t) t m (14)
where the coefficients E L m (p, h; t) depend only on {t} q . Note that when L = V , then S L (tp, h) is polynomial, and coincides with the integral of h on tp.

More precisely, we will show that the Ehrhart coefficients are steppolynomial functions of t, in the sense of [START_REF] Verdoolaege | Counting with rational generating functions[END_REF]. In the next theorem, we prove these results and furthermore, we show that the coefficients E L m (p, h; t) are given by short formulas, provided that in the input (L, p, h), the subspace L has fixed codimension, p is a simple polytope given by its vertices, and the weight is a power of a linear form h(x) = , x M .

Theorem 28. For every fixed number k 0 ∈ N, there exists a polynomialtime algorithm for the following problem. Input:

(I 1 ) a number d ∈ N in unary encoding, with d ≥ k, (I 2 ) a rational subspace L ⊆ R d of codimension k 0 , represented by d -k 0 linearly independent vectors b 1 , . . . , b d-k 0 ∈ Q d in binary
encoding, (I 3 ) a finite index set V, (I 4 ) a simple rational polytope p, given by the set V of its vertices, rational vectors s i ∈ Q d in binary encoding, (I 5 ) a rational vector ∈ Q d in binary encoding, (I 6 ) a number M ∈ N in unary encoding.

Output:

(O 1 ) a finite index set Γ, (O 2 ) polynomials f γ,m ∈ Q[r 1 , . . . , r k 0 ] and numbers ζ γ,m i ∈ Z, q γ,m i ∈ N for γ ∈ Γ, i = 1, . . . , k 0 and m = 1, . . . , d + M , such that for t ∈ R we have S L (tp, , x M ) = d+M m=0 E m (t) t m
, where E m (t) is given by the step-polynomial

E m (t) = γ∈Γ f γ,m {ζ γ,m 1 t} q γ,m 1 , . . . , {ζ γ,m k 0 t} q γ,m k 0 .
Remark 29. Similar results hold for a more general weight h(x), as in [START_REF]How to integrate a polynomial over a simplex[END_REF] for integrals over a simplex. One can assume that the weight is given as a polynomial in a fixed number D of linear forms, h(x) = f ( 1 , x , . . . , D , x ), or has a fixed degree D.

Proof of Theorem 28. As first observed in [START_REF] Baldoni | Local Euler-Maclaurin expansion of Barvinok valuations and Ehrhart coefficients of rational polytopes[END_REF], the path from generating functions to Ehrhart quasi-polynomials relies on the following important property of functions φ(ξ) ∈ M (V * ): such a function has a unique expansion into homogeneous rational functions

φ(ξ) = m≥m 0 φ [m] (ξ).
If P is a homogeneous polynomial on V * of degree p, and D a product of r linear forms, then P D is an element in M (V * ) homogeneous of degree m = p -r. For instance, ξ 1 ξ 2 is homogeneous of degree 0. On this example, we observe that a function in M (V * ) which has no negative degree terms need not be analytic.

For a fixed t ∈ R, the polynomial function ξ → S L (tp, ξ,x M M ! ) is the term of degree M of of the holomorphic function S L (tp)(ξ). Let c s be the supporting cone of p at the vertex s. By Brion's theorem applied to the semi-rational polytope tp, we write S L (tp)(ξ) as the sum of the intermediate generating functions of the cones at the vertices ts, s ∈ V, of the dilated polytope tp. The crucial point is that the cone c s does not change when the polytope is dilated.

S L (tp)(ξ) = s∈V S L (ts + c s )(ξ). (15)
Let q s be the smallest integer such that q s s is a lattice point. Then (t -{t} qs )s is a lattice point, therefore we have S L (ts + c s )(ξ) = e ξ,(t-{t}q s )s S L ({t} qs s + c s )(ξ) = e t ξ,s e -ξ,{t}q s s S L ({t} qs s + c s )(ξ).

Expanding e t ξ,s in powers of t and looking for the coefficient of t i , we obtain

S L tp, ξ,x M M ! = i≥0 t i s∈V ξ,s i i!
e -ξ,{t}q s s S L ({t} qs s + c s )(ξ) [M -i] .

Thus we must analyze the term of ξ-degree M -i of e -ξ,{t}q s s S L ({t} qs s + c s )(ξ).

We use Formula [START_REF] De Loera | LattE, version 1.2[END_REF] and consider the summands indexed by γ ∈ Γ, one at a time. From then on, the proof is entirely similar to the proof of Theorem 37 in [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF], and we omit it.

Example 30 (Example 27, continued). Table 2 shows the output of our Maple program for the Ehrhart quasi-polynomials of the four triangles of Example 27 with respect to the weight h = 1. Note that the dilating parameter t is real. We obtain Setting aside considerations of efficient computation, we can prove the following theorem.

Theorem 31 (Real Ehrhart Theorem). Let p ⊆ R d be a rational polytope, h be a polynomial function of degree M , and L ⊆ R d be a rational subspace of codimension k 0 . Then t → S L (tp, h) is a quasi-polynomial function of t ∈ R >0 . More precisely,

S L (tp, h) = d+M m=0 E L m (p, h; t) t m , (16)
where E L m (p, h; t) is a step-polynomial of the form E m (t) = γ∈Γ f γ,m {ζ γ,m 1 t} q γ,m 1 , . . . , {ζ γ,m k 0 t} q γ,m k 0

for t ∈ R >0 ,
where Γ is a finite index set, with polynomials f γ,m ∈ Q[r 1 , . . . , r k 0 ] and numbers ζ γ,m i ∈ Z, q γ,m i ∈ N for γ ∈ Γ, i = 1, . . . , k 0 and m = 1, . . . , d + M .

Proof. The polytope p is no longer assumed to be simple, so we compute a triangulation of each of the vertex cones and use inclusion-exclusion to avoid overcounting. Using a decomposition of h into powers of linear forms, then the method of the previous theorem gives the result.

v 2 v 1 Figure 1 .

 11 Figure 1. The projected lattice Λ {1}

  (a) Define the semi-open cones k i by

  2) with w1 = v, w2 = w d etc., and v = w 1 . The corresponding semi-open cones are k1 = c(v, w d , . . . , w r+1 ) + c(w r-1 , . . . , w 2 , w 1 ) = k r , . . . , kr = c(v, w d , . . . , w 2 ) = k 1 , thus we obtain (a.1).For i ∈ {1, . . . , r} ∪ {s + 1, . . . , d}, the (closed) cone c i differs from the semi-open cone k i by the union of some faces of lower dimension, therefore (b) follows from (a) and inclusion-exclusion relations.

  [c] ≡ [c((-1, 2, 0), (-1, 0, 0),(-5, 6, 6))] -[c((1, 0, -3), (-1, 0, 0), (-5, 6, 6))] . This example will be continued in Example 26.

(I 1 )

 1 a number d in unary encoding, (I 2 ) a simplicial cone c = c(v 1 , . . . , v d ) ⊂ R d , represented by the primitive vectors v 1 , . . . , v d ∈ Z d in binary encoding, (I 3 ) a rational subspace L ⊆ R d of codimension k 0 , represented by d -k 0 linearly independent vectors b 1 , . . . , b d-k 0 ∈ Z d in binary encoding, compute the following output in binary encoding:

3 5( 1 - 2 5 y+ 3 5

 3123 z )(1 -e -x+3 z )(y + z) + e {a}(-x+2 y)+{2 a+b-c}(-e -x+2 y )(1 -e -z )(y + z)

1 2 1 3

 11 x+y x (-5 x + 6 y + 6 z) x-z x (-5 x + 6 y + 6 z) .

Figure 3 .Table 1 .

 31 Figure3. The square q subdivided into 4 triangles T i . Lattice subspaces parallel to the subspace L = R(0, 1) intersect the walls between the triangles at most in a point, and thus S L vanishes on the intersections of the triangles.Table1. Intermediate valuations of the square and the 4 triangles

4 i=1S

 4 L (tT i ) = 16 t 2 + 4 t -4 t{4 t} 1 -1 30 {-5 t} 2 + 1 30 {5 t} 2 + 1 60 ({-5 t} 2 ) 2 -1 60 ({5 t} 2 ) 2 = 16 t 2 + 4 t -4 t{4 t} 1 ,using a relation between fractional parts, {x} 2 -{-x} 2 = {x} -{-x}, for the simplification. Indeed, a direct calculation gives the same answer.

Table 2 .

 2 Intermediate Ehrhart quasi-polynomials of the square and the 4 triangles p Ehrhart quasi-polynomial S L (tp, 1) -7 8 ({4 t} 1 ) 2 -7 10 {-5 t} 2 + 7 8 {4 t} 1 + 7 20 ({-5 t} 2 ) 2 + 11 6 t -11 24 ({4 t} 1 ) 2 + 11 24 {4 t} 1T 3 3 t 2 + 2 -4 {4 t} 1 t -1 3 ({-5 t} 2 ) 2 + 2 3 {-5 t} 2 + 4 3 ({4 t} 1 ) 2 -4 3 {4 t} 1

	T 1 4 t 2 T 2 21 22 3 t 2 T 4 5 12 t 2 + 1 6 t -1 60 ({5 t} 2 ) 2 + 1 30 {5 t} 2

q 16 t 2 + 4 -4 {4 t} 1 t

Acknowledgments References Acknowledgments This article is part of a research which was made possible by several meetings of the authors, at the Centro di Ricerca Matematica Ennio De Giorgi of the Scuola Normale Superiore, Pisa in 2009, in a SQuaRE program at the American Institute of Mathematics, Palo Alto, in July 2009 and September 2010, and in the Research in Pairs program at Mathematisches Forschungsinstitut Oberwolfach in March/April 2010. The support of all three institutions is gratefully acknowledged. V. Baldoni was partially supported by the Cofin 40%, MIUR. M. Köppe was partially supported by grant DMS-0914873 of the National Science Foundation.


Serge Guelton

Mehdi Amini

Ronan Keryell

Béatrice Creusillet

Pyps

PyPS a programmable pass manager

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Surasak Phoemsapthawee

Marc Le Boulluec

François Jean-Marc Laurens

A potential-flow

come

× × × ∇ ρ φ φ θ θ ψ - Σ Σ -Σ Σ - - - × × × × × × × × × ρ -∇ ρ × -∇ × σ μ υ - υ × ρ ρ π α Λ Λ π Λ α Λ π ρ ω - ω
produce very similar results as expected. This will be the case for most trajectories simulations but in the next section, we will present a case where both approaches yield very different results.

Helical trajectory

In this section, an examination of the helical trajectory is explored and tested by commanding the glider to dive with positive trim and heel command angles with the intention to change the glider direction. In the example presented here, trim and heel command angles are both equal to 10 degrees. Three different geometries of the stabilizer are considered. The dimensions of the three stabilizers are illustrated in Fig. 10.

Fig.10 Dimensions of stabilizers

The glider is expected to turn because of the lateral component of the main wing lift force when the glider is heeling. In the case of a positive heel angle as in these simulations, the glider is then expected to turn to the starboard direction. The obtained trajectories for both approaches are presented in Fig. 11. Unlike the sawtooth simulations, the two approaches produce very different results.

According to the numerical results of the simple approach, the glider moves to starboard as expected. The stabilizer size is only of little effect; the three trajectories are almost the same.

The glider turns with the steering force generated by the main wing and the stabilizer generates the hydrodynamic moment to adjust the glider orientation to the incident inflow. However, this behavior is not always observed in reality. The potential flow results retrieve the experimentally observed counter-steering behavior. The turning equilibrium conditions of the three stabilizers are compared in Table 1. These numerical results confirm that the stabilizer size plays an important role in this behavior. The smallest stabilizer (stabilizer I) causes the glider to turn to the counter-steering direction (the port direction in this case) while the other two do not. When the stabilizer is large enough (stabilizer III), the glider behaves as expected and predicted by the simpler approach. In the case of the intermediate size (stabilizer II), the stabilizer is not large enough to steer the glider properly. When the stabilizer is too small, it cannot produce enough hydrodynamic moment to counteract the counter-steering hydrodynamic moment. The counter-steering hydrodynamic moment is generated by the main wing lift lateral component pointing to starboard (see Fig. 12) and the lateral force due to the body drifting pointing to port. If not sufficiently countered, the moment will make the glider Sterne turn to port because the body drifting drag force center is situated well in front of the wing. In this case the source of the counter-steering moment was hence generated by the lateral friction force on the drifting body. The simple hydrodynamic model does not have this lateral force because as described in section 2.2, in this simple parametric model, the friction is only applied along the main axis. Once it starts to turn to port, the stabilizer sees a flow from the port direction and its lift points to starboard generating a moment at port helping the glider to find a stable route in the port direction.

With higher accelerations, the added mass tensor could take over this effect. The way the body added mass tensor is taken into account is the same for both models. In such cases, both models would see the same counter steering moment but it does not necessarily imply that they would predict the same trajectory because of the fluid inertia effect on the lifting bodies. In the simple model, the hydrodynamic moment generated by the stabilizer is fully perceived as soon as the glider heels without any delay. In the potential flow simulations, there is a delay between the geometric position and the hydrodynamic response. This delay exists because of the fluid inertia applied by the stabilizer. Because of this delay, if the stabilizer counteracting moment is not strong enough, the glider could find the time to position itself in the other equilibrium state causing the counter-steering behavior.

In principle, it is possible to add such effects and other effects in a parametric model. For instance, we could add the interaction between the glider body and its appendages using a correction factor as suggested by [START_REF] Caldeira-Savaira | The active control of swath motions[END_REF]. Adding some of these parameters would not significantly increase CPU time, but it would rapidly make the parametric model look like a patchwork difficult to control and maintain. Before any additional parameter is included within the parametric model, its importance and effect have to be clearly identified. Furthermore, not all effects are easily replaced by a simple model involving a single parameter. The lateral friction drag for instance depends on the Reynolds number but also on the body geometry. Although the potential flow simulator does not take all the physics into account, it can be used as a numerical towing tank facility to improve the simulator, i.e. the parametric model. Its CPU time is not negligible compared to the simple hydrodynamic code but it consists of a stage between the parametric model and the experimental trials.

Computation resource consumption

We now consider the underwater glider equipped with the stabilizer III to simulate a non trivial trajectory. The objective is to simulate a helical diving trajectory followed by a surfacing contra-rotating helical trajectory. As in the previous simulations, the glider is launched with a 0.2 m/s horizontal velocity in the b x -positive direction. The heel command angle is fixed at 20 degrees. The trim command and the ballast are varied alternatively every 1100 seconds. For the first 1100 seconds, the trim command angle is 10 degrees and the ballast takes the water in for diving. For the next 1100 seconds, the trim command angle is -10 degrees and the ballast drains the water off for surfacing. The total duration of the simulation is 2332 seconds and takes 20000 time steps. This simulation involves about 24 hours of computation time on a standard workstation (CPU 4 cores with 2.66 GHz).

Conclusion and perspective

An Euler-Newton equations solver is coupled with a potential flow code to simulate 6-DOF trajectories of underwater gliders. This simulator can be used to study the hydrodynamic behavior of gliders in order to improve the automatic flight control and to optimize the geometry of the glider. A numerical study of the hydrodynamic behavior of an underwater glider has been conducted. All results are compared with a simple parametric simulator. A series of simulations considering sawtooth trajectories was first conducted. In this case both simulators presented very similar trajectories. As expected, the glider velocity varies as a function of the trim angle. To confirm a behavior observed experimentally, a second series of simulations examining the glider steering was conducted. The potential flow simulator shows that the stabilizer geometry plays an important role in steering control. Inappropriate stabilizer geometry can cause counter-steering behavior that the parametric simulator cannot anticipate. A non trivial case was finally presented to demonstrate the capabilities of the potential flow based glider simulator. Enhancing the parametric model with additional features is easy, but simple models are the best way to determine which physical effect is really important. The potential flow code is a good practical tool to simulate complex trajectory patterns and to verify whether the parametric code predicts the same trajectory. However, further validations utilizing real data and experimental model testing are necessary to increase confidence in the potential simulator before it can be adopted as a template for the parametric model enhancements.

Fig. 13 Example of trajectory simulation for a real-time duration of 2332 seconds, of which the CPU time is about 24 hours on a standard PC (the glider appears 10 times bigger than it is) Good quality experimental results will be extremely difficult to obtain. The behavior of the Sterne glider has been obtained from sea trials. The onboard instruments record the velocities, depth, heel and pitch angles. The trajectory is deduced from these data. There is no information concerning the environmental conditions at sea. To conduct a proper experimental investigation at scale 1, we would need a calm pool of water with very large dimensions. Although the glider is only 2 meters long, its gyration radius is about 30 meters. In Fig. 13, the glider appears much larger than in reality as stated in the caption. We plan to perform some experimental trials on a reduced size model at the Ifremer tank facility in Brest.

Since the scale cannot be too small, we will not be able to conduct anything resembling the case of Fig. 13 in the tank. However, a full gyration is not necessary to study the counter-steering behavior. The scale of the experimental model is the key issue since it will be a compromise between several factors including the Reynolds number, manufacture and the observable trajectories in the Ifremer towing tank facility.

Fig. 11

 11 Fig.11 Top view of helical trajectory during 1000 seconds for 10° trim and 10° heel command angles; top with potential flow approach, bottom with simple approach

Table 1 Turning equilibrium conditions from the potential flow approach results; the gliding diving with 10° trim and 10° heel command angle

 1

	Stabilizer	I	II	III
	Turning radius (m)	22	171	39
	Turning rate (deg/s)	0.88	-0.11	-0.48
	Advance speed (m/s)	0.334	0.339	0.329
	Diving speed (m/s)	0.062	0.058	0.056
	AOA of stabilizer (deg)	8.85	3.95	2.35
	AOA of main wing (deg)	1.68	1.55	1.63
	Heel angle (deg)	12.56	9.68	8.51
	Trim angle (deg)	7.08	7.51	7.67

Fig.12 Force diagram in steering mechanics

Acknowledgement

This work is partially supported by the Europole Mer (http://www.europolemer.eu).

Laurens Jean-Marc is a senior lecturer in naval hydrodynamics at Ensta Bretagne, Brest, France. His main research fields include resistance and propulsion, ship and underwater vehicles maneuverability.

Deniset François is an associate professor in naval hydrodynamics at french naval academy, Brest, France. His main research fields include performance of hydrofoils and marine propulsion systems, fluid-structure interactions and cavitation.