Jonatha Anselmi
email: jonatha.anselmi@inria.fr

Bruno Gaujal
email: bruno.gaujal@inria.fr

Louis Sébastien Rebuffi
email: louis-sebastien.rebuffi@univ-grenoble-alpes.fr

Optimal Speed Profile of a DVFS Processor under Soft Deadlines

Keywords: Energy minimization, Optimal speed, DVFS, Markov decision processes

We consider a Dynamic Voltage and Frequency Scaling (DVFS) processor executing jobs with obsolescence deadlines: A job becomes obsolete and is removed from the system if it is not completed before its deadline. The objective is to design a dynamic speed policy for the processor that minimizes its average energy consumption plus an obsolescence cost per deadline miss. Under Poisson arrivals and exponentially distributed deadlines and job sizes, we show that this problem can be modeled as a continuous time Markov decision process (MDP) with unbounded state space and unbounded rates. While this MDP admits a continuous time optimality equation for its average cost, the standard uniformization approach is not applicable. Inspired by the scaling method introduced by Blok and Spieksma, we first define a family of truncated MDPs and we then show that the optimal speed profiles are increasing in the number of jobs in the system and are uniformly bounded by a constant. Finally, we show that these properties are inherited from the original (infinite) system. The proposed upper bound on the optimal speed profile is tight and is used to develop an extremely simple policy that accurately approximates the optimal average cost in heavy traffic conditions.

Introduction

Minimizing the energy consumption of embedded systems with real-time execution constraints is becoming more and more important. More functionalities and better performance/cost trade-offs are expected from such systems because of the increased use of real-time applications and the fact that batteries are becoming standard power supplies. Dynamically changing the speed of the processor is a common and efficient way to reduce energy consumption and remarkable gains can be obtained when considering cacheintensive and/or CPU-bound applications as the CPU energy consumption may dominate the overall energy consumption [START_REF] Snowdon | Power Management and Dynamic Voltage Scaling: Myths and Facts[END_REF]. In fact, this is the reason why modern processors are equipped with Dynamic Voltage and Frequency Scaling (DVFS) technology [START_REF] Weiser | Scheduling for Reduced CPU Energy[END_REF].

In the deterministic case where job sizes and arrival times are known, a vast literature addressed the problem of designing both off-line and on-line algorithms to compute speed profiles that minimize the energy consumption subject to hard real-time constraints (deadlines) on job execution times; see, e.g., [START_REF] Yao | A scheduling model for reduced CPU energy[END_REF]; [START_REF] Bansal | Speed Scaling to Manage Energy and Temperature[END_REF]; [START_REF] Li | An O(n 2 ) Algorithm for Computing Optimal Continuous Voltage Schedules[END_REF] and the references therein. In a stochastic environment where only statistical information is available about job sizes and arrival times, it turns out that combining hard deadlines and energy minimization via DVFS-based techniques is much more difficult. In fact, forcing hard deadlines requires to be very conservative, i.e., to consider the worst cases. In spite of these difficulties, this problem has been investigated in [START_REF] Lorch | Improving Dynamic Voltage Scaling Algorithms with PACE[END_REF] for a single job and in [START_REF] Gaujal | Dynamic Speed Scaling Minimizing Expected Energy Consumption for Real-Time Tasks[END_REF] for multiple jobs. The former approach constructs the optimal speed profile explicitly in "closed form" while the latter relies on the numerical solution of a discrete time Markov Decision Process (MDP) [START_REF] Martin L Puterman | Markov decision processes: discrete stochastic dynamic programming[END_REF]. The latter approach has several drawbacks: i) it requires a discretization of both time and space, which introduces by itself an approximation on the optimal solution, ii) deadlines and job sizes need to be bounded, and iii) the size of the state space of the underlying MDP is exponential in the size of the maximal deadline. These issues makes this approach unusable in practice.

The approach followed in this paper circumvents the difficulties described above by replacing the hard real-time constraints, i.e., jobs have hard deadlines that must be satisfied, by soft real-time constraints, i.e., jobs may miss their deadlines, at some cost. While the hard deadline of a job must be known at the job arrival, soft deadlines allow for a different information structure: here, only the deadline distribution is known at the job arrival. In this paper, we further assume that jobs missing their deadlines become obsolete and are dropped. Obsolescence is often found in real-time systems where the information carried by jobs may not be valid any longer after their deadline as it will be replaced by fresher input coming from other jobs. Therefore, obsolete jobs become useless and can get discarded from the queue. Dropping obsolete jobs can also model impatient customers: customers wait for service for some time (deadline) and quit (are dropped) if not served before that time.

Contribution

We investigate the problem above in a Markovian setting where jobs join the system following a Poisson process and both the deadlines and sizes of jobs are exponentially distributed. Under these assumptions, our goal is to minimize the average cost, i.e., the average energy spent by the processor per second plus the penalty due to jobs missing their deadlines. We formulate this problem as an MDP in continuous time where the state is the number of jobs in the system and the action is the processor speed.

Our main result, Theorem 2.2, shows the existence of an optimal speed profile that is increasing in the number of jobs in the system and upper bounded by some constant. This constant is defined in (3) as the minimizer of a function that comes out from our analysis. Surprisingly, our bound does not depend on the deadlines and arrival rates. In other words, our bound on the optimal speed does not change upon variations of these job characteristics. In addition, it yields a simple approximation for the optimal policy and several numerical tests show that such approximation is accurate in heavy-traffic conditions. Finally, the proposed approximation is used to control the proportion of jobs that leave the system because they missed their deadline in a simple manner.

Underlying the proof of our main result, there are some technical challenges that we now discuss. The proposed MDP satisfies the regularity assumptions (stability, unichain) needed to establish an optimality equation as described in [START_REF] Guo | Continuous-time Markov decision processes[END_REF]. However, this is not enough to show structural properties of the optimal policy. In fact, the classical approach to do this is to uniformize the MDP and to investigate the properties of the corresponding discrete time value iteration operator. Unfortunately, this is not possible in our case because the transition rates are unbounded. To uniformize the MDP, a typical approach consists of truncating the state space. Indeed, this is the approach we follow. However, we notice that a naive truncation will not help because the truncation barrier has a strong impact on the structure of the optimal policy in the sense that it would not preserve any monotonicity property that it may have without truncation. This is shown in Figure 2. Instead, we use the technique proposed in [START_REF] Blok | Countable state Markov decision processes with unbounded jump rates and discounted cost: optimality equation and approximations[END_REF], which smoothly scales down the upward rates of the truncated system as a function of the size of its state space. This technique has been successfully used in Hyon and Jean-Marie (2020); [START_REF] Bhulai | On structural properties of the value function for an unbounded jump Markov process with an application to a processor sharing retrial queue[END_REF] to show structural properties of controlled queueing systems. However, these works focus on discounted costs. Here, we use the same truncation technique but we apply it to the average cost. To the best of our knowledge, this has never been done before. In our specific case, the convergence to the infinite system will be guaranteed by the monotone convergence theorem.

Our paper is organized as follows. In Section 2, our model and the corresponding MDP are described in detail. We also present our main result (Theorem 2.2) as well as some hindsight on the construction of the proof. Section 3 shows how the MDP is truncated and scaled and shows the proof of the monotonicity as well as the construction of the upper bound on the optimal speed. Section 4 focuses on the convergence when the truncation point goes to infinity. Section 5 uses Theorem 2.2 to provide an approximation on the optimal policy and estimate the deadline-miss probability. Finally, Section 6 draws the conclusions of our work and addresses further research.

System Description, Problem Statement and Main Result

The system described here is a model for the dynamics of a real-time device composed of a single computing resource (a processor) where incoming jobs need to be executed under a constraint on the amount of time that they spend in the system.

Processor. This is a DVFS processor whose speed can continuously vary in the interval [0, S max ]. We consider that speed changes are immediate and induce no energy cost. When the processor works at speed s, it processes s units of work per second while its power dissipation is w(s) watts. The classic simple model for the dynamic power dissipation of any CMOS circuit is w = KαV 2 f , see, e.g., [START_REF] Snowdon | Power Management and Dynamic Voltage Scaling: Myths and Facts[END_REF], where K is a constant, α measures the activity of the logical gates, V is the supply voltage and f is the clock frequency. The clock frequency of the gates is often linearly related to the voltage and therefore DVFS processors adjust both variables together. Within the model above, this means that w(s) is cubic in the speed s. In this paper, we just require that w(s) is continuous, increasing and strictly convex in the speed s.

Jobs. They form a stochastic point process, with Poisson arrivals with rate λ, i.i.d. deadlines exponentially distributed with rate δ and i.i.d. sizes exponentially distributed with rate µ. Without loss of generality, we assume that µ = 1.

Dynamics. At any point in time t, the processor chooses its speed s(t) and executes one of the jobs in its backlog queue. We notice that the choice of the job in execution, named active in the following description, is irrelevant here because of the memory-less property of the deadlines and of the sizes. Thus, at any point in time, at most one job can be active. As mentioned above, this induces an instantaneous energy cost of w(s(t)). Now, three events can happen in continuous time:

1. A new job may join the queue. 2. The active job is completed before its deadline. In this case, the job leaves the system. 3. One job (active or inactive) reaches its deadline. In this case, this job becomes obsolete, it is removed from the queue and an immediate cost equal to C is paid.

Cost Function. If we denote by M T the number of missed deadlines in the time interval [0, T ], the objective of this paper is to study the speed profile s(t) of the processor that minimizes the long-run average cost given by the missed deadlines plus power consumption, say J . Specifically, this is given by

J := lim sup T →∞ 1 T CM T + T 0 w(s(t)) dt .
At this point, we claim that this problem can be modeled by a continuous time Markov decision process with a discrete state space. To see this, let us consider the system at time t under the speed profile s(•) and let i(t) denote the number of jobs present in the system at time t. Each job is characterized by two values: (a k , b k ), respectively, the sojourn time of job k (how long the job has been present) and its service quantity (how many elementary operations have already been executed on this job by the processor). Also let r t be the time elapsed since the last arrival. Then the state of the system at time t is the tuple (r t , (a 1 , b 1 ), . . . , (a i(t) , b i(t) )). From time t, the cost function only depends on the expected future behavior of the system, namely, the future arrivals, future services and future drops of packets due to deadline misses. The memory-less property of the exponential distribution implies that this expected future behavior does not really depend on (r t , (a 1 .b 1 ), . . . , (a i(t) , b i(t) )), but only on i(t), the current number of jobs. As a by-product, this implies that the choice of the optimal speed at time t only depends on i(t), so that the optimal speed profile can only change when i(t) changes, i.e., when an arrival, a service or a drop occur.

Markov Decision Process. We now formulate the problem of interest as an MDP. As mentioned before, the state space is N and a state represents the number of jobs in the system. The action space is [0, S max ], i.e., the set of available speeds for the DVFS processor. Let σ = (σ i ) i∈N denote a stationary and deterministic speed policy adopted by the processor, i.e., σ i ∈ [0, S max ] is the speed used in state i. It is well known that focusing on stationary and deterministic policies can be done with no loss of optimality in our case (Guo and Hernandez-Lerma, 2009, Theorem 5.9); in other words, we may consider the broader set of history dependent randomized policies and then use that there exist stationary and deterministic policies that are optimal with respect to the cost defined below. For i, j ∈ N, the transition rates are given by:

q i,j (σ) =                λ if j = i + 1 (arrival of a new job) σ i + iδ if i ≥ 1 and j = i -1 (completion of the active job or obsolescence of one job) -λ -iδ -σ i if j = i 0 otherwise.
Under the speed policy σ = (σ i ) i∈N , the induced Markov chain, denoted by X σ , is a birth and death process that resembles an M/M/∞ queue but with an additional decreasing rate, which comes from the processing speed; see Figure 1.

By the ergodicity of the Markov chain X σ under all policies σ, the longrun cost J is equal to the long-run expected cost. Letting E σ i 0 denote the expectation given a speed policy σ and starting state i 0 , we have

J = J(i 0 , σ) := lim sup T →∞ 1 T T 0 E σ i 0 c(X σ (t), σ) dt. 0 • • • i i + 1 (i + 1)δ + σ i+1 λ • • • Figure 1: Markov chain diagram under policy σ.
In this equation, the immediate cost function c(•, •) is the expected cost incurred by the system at time t. It only depends on the current state and the current speed. Conditional on the state (X σ (t) = i), the obsolescence rate is iδ. Thus, the expected cost is:

c(i, σ) := Ciδ + w(σ i ).
With a slight abuse of notation, we will use both notations c(i, σ) or c(i, σ i ) since c only depends on the speed used in state i and not on the whole policy.

As mentioned before, for given speed policy σ, X σ is ergodic so that the MDP is unichain (all states are positive recurrent under all policies). This implies that the cost can be defined independently of the starting state. Moreover, for a given σ, there exists a unique stationary measure π σ for X σ so that we can define the cost independently of the initial state and express it as a function of the stationary measure:

∀i 0 , J(σ) := J(i 0 , σ) = E π σ c(X σ , σ) = i π σ i c(i, σ i ). (1)
Here, E π σ is the expectation with respect to the invariant measure of X σ . Stationary policies that minimize (1) are optimal speed policies for the model.

In particular, they are also optimal over all policies (history dependent and randomized) (Guo and Hernandez-Lerma, 2009, Theorem 5.9) . Also, our MDP satisfies all the conditions given in (Guo and Hernandez-Lerma, 2009, Theorem 5.9) to assert the existence of an optimal stationary deterministic policy σ * and an optimality equation of the form

J * = J(σ * ) = c(i, σ * i ) + j h * (j)q i,j (σ * i ) = min s∈[0,Smax] c(i, s) + j h * (j)q i,j (s), ∀i ∈ N, (2)
where h * is a real function defined on N, usually referred to as bias of the optimal policy.

Main result. The goal of this paper is to investigate structural properties on σ * and J * . First, let us define B as

B := arg min s∈R + (w(s) + C(λ -s)). ( 3
)
This constant is well defined in R + ∪ {+∞} because w is strictly convex. In addition, we have the following remark.

Remark 2.1. If w is super-linear, i.e., lim s→∞ w(s) s = ∞, then w(s)+C(λ-s) is also super-linear and B is finite. In practice, all models of power dissipation are super-linear in s, e.g., [START_REF] Chandrakasan | Low-power CMOS digital[END_REF]. In the simple case where w(s) = Ks 3 , we first notice that the constant K can be set to 1 without loss of generality because to compensate its effect one can adjust the missed deadline cost C accordingly. Then, in this case, we obtain B = C 3 . Our main result is the following.

Theorem 2.2. There exists a deterministic optimal policy σ * = (σ * i ) i∈N that is increasing in i and upper bounded by B.

Remark 2.3. The optimal speed policy of the processor is always bounded by a finite constant, namely min(B, S max ). By definition, B only depends on w (the power dissipation of the processor) and C the cost of each missed deadline. Thus, we remark that B is independent of the job characteristics (arrival rate, deadline and size distributions). This is both surprising and helpful in practice. Indeed, if B is finite, one can set a priori the maximal speed of the processor to S max := B. This guarantees that in all cases, no cost reduction would be possible by using a more powerful processor. Further discussion about parameter settings, in particular the link between C and the probability that jobs miss their deadline under the optimal policy, will be discussed in Section 5.

A proof of Theorem 2.2 is developed in Sections 3 and 4. Before delving into the proof, we devote the remainder of this section to explain the technical difficulties underlying our problem and the general approach that we follow.

The optimality equation ( 2) cannot be uniformized because the rates q ii (σ) are unbounded in i. Therefore, the study of structural properties of the optimal policy must be done by constructing a sequence of truncated MDPs whose optimal policies converge to σ * and for which we can prove monotonicity and boundedness. This approach has been proposed for the first time in [START_REF] Blok | Countable state Markov decision processes with unbounded jump rates and discounted cost: optimality equation and approximations[END_REF] for MDPs with discounted cost by truncating the state space and scaling the rates of all the events that take the system out of the truncated space. This has been successfully applied in, e.g., [START_REF] Bhulai | On structural properties of the value function for an unbounded jump Markov process with an application to a processor sharing retrial queue[END_REF]; Hyon and Jean-Marie (2020), to show that threshold type policies yield optimal admission control in one queue. However, all these applications consider discounted costs. To the best of our knowledge, no work has been done for the average cost. In the following, we will show that in our case, the scaling technique of [START_REF] Blok | Countable state Markov decision processes with unbounded jump rates and discounted cost: optimality equation and approximations[END_REF] also works for the average cost and is the key ingredient to show Theorem 2.2, which gives new hindsights on the optimal policy. Thus, our result is another evidence of the power of this scaling approach, though our proof is quite different from the approach used in the discounted case. In fact, the classical approach is to show that the value iteration operator preserves structural properties of the cost and of the policy (typically convexity properties of the cost and level sets of the policy), so that successive iterations of the operator will also preserve the properties and converge to the optimal cost/policy. Here, we will directly consider the fixed point optimality equation and prove monotonicity of the policy by induction on the state (see Sections 3.1, 3.2, 3.3).

Truncated Model

As mentioned before, the original MDP cannot be uniformized because the transition rates q ii (σ) grow to infinity when i goes to infinity. To construct a discrete time model, we truncate the state space following the guidelines from [START_REF] Blok | Countable state Markov decision processes with unbounded jump rates and discounted cost: optimality equation and approximations[END_REF] to construct a finite state MDP M N with linearly decaying arrival rates. The new state space is {0, • • • , N }, with the same continuous action space [0, S max ] and an average cost per second J N (σ) defined as

J N (σ) := lim sup T →∞ 1 T T 0 E σ i 0 c(X σ (t), σ) dt (4)
and with transition rates for all i, j ≤ N given by

q N i,j (σ) :=          λ N i if i < N and j = i + 1 σ i + iδ if i > 0 and j = i -1 -λ i -iδ -σ i if j = i 0 otherwise,
where the decaying arrival rate is λ N i := λ 1 -i N . Using decaying arrival rates will be a key ingredient in this paper. To illustrate this, let us also consider a naive truncated MDP, M N with fixed arrival rates in each state given by λ. The state space is {0, • • • , N }, with the same continuous action space [0, S max ] and the same cost function but with modified transition rates given by

q N i,j (σ) :=          λ if i < N and j = i + 1 σ i + iδ if i > 0 and j = i -1 -λ -iδ -σ i if j = i 0 otherwise,
The two finite MDPs M N and M N have been solved numerically using the following parameters: λ = 10, δ = 0.14, C = 200, w(v) = v 3 , N = 60. The respective optimal policies σ * and σ * are displayed in Figures 2(a)-(b). As one can see, the two models behave very differently. At the last state of M N , the speed does not need to be as high as the maximal speed, as the arrival rate drops from λ to 0, so that the speed plummets near the last state. In contrast, the optimal policy for M N is increasing from 0 to some bound slightly below C 3 . Additional numerical tests where we let N grow (reported in Figure 3) further suggest that σ * is increasing in the state i as well as in the level of truncation N . The bound B also appears to be rather tight when N is large in the example reported in Figure 3. This makes the study of M N promising and, in the remainder, we focus on this MDP with decaying arrival rates.

Since N is fixed here, we may remove it in the notation for simplicity. As the state space is finite, we can uniformize this MDP to get a discrete time MDP. Choosing

U := λ + N δ + S max (5)
as uniformization constant, we get a discrete time MDP D N with transition probabilities given by

p N i,j (σ) =          1 U λ i if i < N and j = i + 1 1 U (σ i + iδ) if i > 0 and j = i -1 1 U Ūi if j = i 0 otherwise,
where the complementary probability to stay in state i is Ūi := U -λ i -δi-σ i . Again, we focus on stationary policies σ that minimize the cost J N (σ). The long-run average cost per step for the discrete time MDP is

g(σ) := lim sup K→∞ 1 K K k=1 E σ i 0 c(X σ (k), σ).
A classical result for MDPs (see for example (Puterman, 2014, Section 11.5.3)) says that the discrete and continuous time MDPs are related by the following relations: for any policy σ, U g(σ) = J N (σ). This has the following consequences:

Both MDPs have the same optimal policies. Optimal long-run average costs coincide up to the multiplicative uniformization constant:

U g * = U g(σ * ) = J N (σ * ) = J *
N . We will first show the counterpart of Theorem 2.2 in the finite case. Let us define

B N := arg min s∈R + w(s) + C(λ -s) 1 + λ δN , (6)
which is well defined in R + ∪ {∞} and unique because w is strictly convex.

We have the following properties for the optimal speed policy.

Theorem 3.1. (i) The optimal speed policy σ * minimizing (4) is unique.

(ii) The optimal speed policy is increasing in i: ∀i ≤ N, σ * i < σ * i+1 .

(iii) The optimal speed policy is upper-bounded: ∀i ≤ N, σ * i ≤ B N . The proof of this result will be the object of Section 3.1 (monotonicity, item (ii)), 3.2 (upper-bound, item (iii) ), 3.3 (uniqueness, item (i)). Then, the relation with the original infinite MDP will be shown in Section 4.

Proof of Theorem 3.1(ii): Monotonicity of the Optimal Speed

In this subsection, we denote by σ * any optimal speed policy, X * the associated Markov chain and J * N the optimal cost induced by σ * . Thus, the optimal long-run average gain per step is

g * := lim K→∞ 1 K E σ * K k=1 1 U c X σ * (k), σ * . ( 7
)
When the state space is finite, the bias H ∈ R N +1 for the optimal speed σ * is defined up to an additive constant by

H(i) := E σ * i ∞ k=1 c X σ * (k), σ * -U g * , ∀i ≤ N. ( 8
)
To fix the value of the bias vector, we set H(0) := 0. Since the MDP is finite, unichain, the action space is compact and the costs and transition probabilities are continuous and bounded in the actions, (Puterman, 2014, Theorem 8.4.7) guarantees the existence of the optimality equations for the optimal cost and for the bias. Specifically, for any state i ∈ {0, . . . , N },

g * + H(i) = 1 U min s∈[0,Smax] c(i, s) + (δi + s)H(i-1) + Ūi H(i) + λ i H(i+1) (9)
with H(-1) = H(N + 1) = 0 by convention.

For each state i, an optimal action σ * i is the choice of a speed minimizing the right hand side term. Notice that necessarily, σ * 0 = 0 (the speed of the processor must be 0 when there is no work to do).

Using (9), for i ≥ 1, we can subtract H(i -1) from H(i) and choose s = σ * i-1 in (9) to get

U (H(i) -H(i -1)) ≤ δC + λ 1 - i N (H(i + 1) -H(i)) + Ū i - λ N (H(i) -H(i-1)) + (δ(i -1) + σ * i-1 )(H(i-1) -H(i-2)) (10) with Ū i := Ūi + σ * i -σ * i-1 .
Similarly, doing the same subtraction for s = σ * i , we obtain

U (H(i) -H(i -1)) ≥ δC + λ 1 - i N (H(i + 1) -H(i)) + Ūi - λ N (H(i) -H(i-1)) + (δ(i-1) + σ * i )(H(i-1) -H(i-2)). ( 11
)
Combining both inequalities together we get the inequality

σ * i -σ * i-1 (H(i) -2H(i -1) + H(i -2)) ≥ 0.
From this, we can deduce the following property of the model. Proposition 3.2. If H satisfies the following notion of discrete convexity:

∀i ≥ 2, (H(i) -2H(i -1) + H(i -2)) > 0, ( 12
)
then the optimal speed policy σ * is increasing.

Here are some additional properties of the finite MDP that will be needed for the proof of the main result.

Lemma 3.3. The asymptotic cost per second is upper-bounded:

U g * ≤ Cλ 1+ λ N δ .
Proof. The cost U g * is the optimal asymptotic cost per unit of time. Therefore, we have that U g * ≤ U g 0 , with U g 0 being the asymptotic cost per unit of time when the speed is 0 for each state. Let π (0) be the asymptotic distribution in that case. It is not difficult to see that π

(0) i = 1 + λ N δ -N N i λ N δ i
and thus we can then compute the associated asymptotic cost as follows:

g 0 = 1 + λ N δ -N N i=0 π (0) i δiC = δC 1 + λ N δ -N N i=0 i N i λ N δ i = δC 1 + λ N δ -N N i=1 N N -1 i -1 λ N δ i = Cλ 1 + λ N δ -N 1 + λ N δ N -1 = Cλ 1 + λ N δ
, which concludes the proof.

We now want to show by backward induction on i that H satisfies (12). For 0 ≤ i ≤ N -1, the exact property P(i) that we will show is

H(i) -H(i -1) < H(i + 1) -H(i) < C 1 + λ N δ . ( 13
)
Before investigating the initialization step, we show the following preliminary inequality on the optimal cost.

Lemma 3.4. Let σ := E π * [σ * X * ] be the average speed under the stationary measure π * . Then, Proof. Consider the optimal Markov chain X * displayed in Figure 4, where for clarity we represented the continuous Markov chain, as its behavior is similar as the uniformized discrete time Markov chain. Using Jensen's inequality, since w is strictly convex and the stationary measure π * is non trivial, we get the following strict inequality:

U g * > w(σ) + C λ -σ 1 + λ N δ . 0 • • • i i + 1 (i + 1)δ + σ * i+1 λ 1 -i N • • • N
U g * = J * N = E π * [w(σ * X * )] + CδE π * [X * ] > w(E π * [σ * X * ]) + CδE π * [X * ] = w(σ) + CδE π * [X * ]. Now, we calculate E π * [X * ].
The underlying detailed balance equations are

π * i+1 [δ(i + 1) + σ * i+1 ] = π * i λ -i λ N .
The LHS is 0 for i = -1, and the RHS is also 0 for i = N . When we sum these equations from 0 to N -1, we can therefore write

N i=0 π * i (δi + σ * i ) = N i=0 π * i λ -i λ N , which gives δE σ * [X * ] + σ = λ -E σ * [X * ] λ N , i.e., E σ * [X * ] = λ -σ δ + λ N . (14)
Thus, we finally have

U g * > w(σ) + C λ -σ 1 + λ N δ
as desired.

Remark 3.5. The computation of the expectation that gave (14) remains true for any speed policy, and in that case σ becomes the average speed of that policy. This formula of the expectation will be needed later to rewrite the cost J * N in a different way. For the initialization of the induction, we need the following lemma.

Lemma 3.6. H(N ) -H(N -1) < C 1+ λ N δ .
Proof. For any σ, the optimality equation ( 9) at i = N gives

U (g * + H(N )) ≤ w(s) + CN δ + (δN + s)H(N -1) + (U -δN -s)H(N ),
which holds if and only if

U g * + (δN + s)H(N ) ≤ w(s) + CN δ + (δN + s)H(N -1),
which holds if and only if

H(N ) -H(N -1) ≤ w(s) + CN δ -U g δN + s , (15)
with equality if the chosen σ is optimal. We therefore want to show the following inequality for some s:

h(s) := w(s) + CN δ -U g * δN + s < C 1 + λ N δ
To show this, we use the inequality (3.4) on U g * and choose s = σ. This gives

h(σ) < C N δ -λ-σ 1+ λ N δ δN + σ < C N δ + σ (δN + σ)(1 + λ N δ ) < C 1 + λ N δ .
This concludes the proof by (15). Therefore, for the initialization step we write the inequality (11) for i = N to obtain

U (H(N ) -H(N -1)) ≥ δC + U -δN -s - λ N (H(N ) -H(N -1)) + (δ(N -1) + s)(H(N -1) -H(N -2)), which implies (δ(N -1) + s) (H(N ) -H(N -1)) ≥ δC - λ N + δ (H(N ) -H(N -1)) + (δ(N -1) + s)(H(N -1) -H(N -2)),
which implies (using Lemma 3.6)

(δ(N -1) + s) (H(N ) -H(N -1)) > (δ(N -1) + s)(H(N -1) -H(N -2)),
which as desired gives H(N ) -H(N -1) > H(N -1) -H(N -2). This new inequality and Lemma 3.6 imply P(N -1). Now, let us assume that the property P(i) is true for some i ≤ N -1. We then have, writing ( 11) with s = σ * i and using the first inequality of P(i) (13),

U (H(i) -H(i -1)) ≥ δC + U -δi -s - λ N (H(i) -H(i -1)) + (δ(i 1) + s)(H(i -1) -H(i -2)) which implies δi + s + λ N (H(i) -H(i -1)) ≥ δC + (δ(i -1) + s)(H(i -1) -H(i -2)), which implies (δ(i -1) + s) (H(i) -H(i -1)) > (δ(i -1) + s)(H(i -1) -H(i -2)).
This inequality comes from the second inequality of P(i) (13), and from its first inequality we finally obtain both inequalities of P(i -1). Thus, H(i) -H(i -1) > H(i -1) -H(i -2) and the backward induction is complete. Since H is strictly convex, with Proposition 3.2 we deduce that σ * is increasing.

Proof of Theorem 3.1(iii): Upper Bound on the Optimal Speed

By definition, we recall that B N is the unique minimum of u, see ( 6). We also have

J N (σ) = i∈S π i u(σ i ),
where for all a ≥ 0,

u(a) := w(a) + C λ -a 1 + λ N δ .
Proposition 3.7. If σ * is an optimal increasing speed policy, then it is upperbounded by B N , which means σ * i ≤ B N for all i ≤ N . Proof. We show this statement by contradiction. Let σ * be an optimal increasing speed policy. As σ * is increasing, assume that σ * N > B N Let π be the associated stationary measure. As σ * N > B N ≥ 0 = σ * 0 , we can define i 0 > 0 as the smallest state such that σ * i 0 > B N . Now, we can define the following policy

σ i = σ * i if i < i 0 B N if i ≥ i 0
and let π denote its associated stationary measure. By definition, we notice that the minimum of the function u is reached at B N . Moreover, for all i < i 0 , we get from the local balance equations that π i > π i . Indeed, when i ≥ i 0 -1:

π i+1 π i = λ 1 -i N (i + 1)δ + σ * i+1 < π i+1 π i ,
with equality only when i < i 0 -1. It remains to show that the cost associated to this new speed is lower than the original one. We have

J(σ * ) -J( σ) = N i=0 π i u(σ * i ) - N i=0 π i u( σ i ) = i 0 -1 i=0 (π i -π i )u(σ * i ) + N i=i 0 π i u(σ * i ) - N i=i 0 π i u(B N ) > i 0 -1 i=0 (π i -π i )u(B N ) + N i=i 0 π i u(B N ) - N i=i 0 π i u(B N ) > 0.
This contradicts the optimality of σ * and concludes the proof.

Proof of Theorem 3.1(i): Uniqueness of the Optimal Policy

The following lemma gives a different expression for σ * i using the notion of generalized inverse of the derivative w of w. First, since w is strictly convex, w is continuous and well defined everywhere but on a countable subset D ⊆ [0, S max ]. Moreover, w is strictly increasing and diverges, so that we can correctly define the inverse

w -1 : y → inf{s ∈ R, w (s) ≥ y},
and w -1 is increasing.

Lemma 3.8. If σ * is an optimal speed policy with bias function H, we can relate the variation of the bias to the speed at a given state:

σ * i = w -1 (∆ i ), ∀i ≥ 1 with ∆ i := H(i) -H(i -1) for i ∈ [1, N ].
Proof. Let σ * be an optimal speed policy with bias function H. For ≥ 1, the optimality equation ( 9) can be written as

J * N = min s {w(s) + Ciδ -(δi + s)∆ i + λ i ∆ i+1 }, (16)
so that for any speed s:

w(σ * i ) -σ * i ∆ i ≤ w(s) -s∆ i . This yields w(σ * i )-w(s) σ * i -s ≤ ∆ i for s ≤ σ * i w(σ * i )-w(s) σ * i -s ≥ ∆ i for s ≥ σ * i . Let D = [0, S max ]\D. These inequalities on ∆ i give: w (σ * - i ) ≤ ∆ i ≤ w (σ * + i ), so that there are two possibilities. Either i) ∆ i ∈ w D and w (σ * i ) = ∆ i , or ii) ∆ i / ∈ w D , so that σ * i = w -1 (∆ i ) ∈ D.
In both cases, we have σ * i = w -1 (∆ i ). Proposition 3.9. The optimal speed policy is unique, and therefore Blackwell optimal.

Proof. Let σ, σ be two optimal speed policies, H, H their respective biases and ∆, ∆ the respective variations of the biases. We will show by induction that the variation of the speed and biases are equal.

We already have that σ(0) = σ(0) = 0. The optimality equation ( 9) for i = 0 then gives

∆ 1 = ∆ 1 = U g * λ = J * N λ .
Using Lemma 3.8, we then have

σ 1 = w -1 (∆ 1 ) = w -1 ( ∆ 1 ) = σ 1 .
Assume that for some i ≥ 1, σ i = σ i and ∆ i = ∆ i . Writing ( 16) for both optimal speeds, and using the assumption of the induction, we have ∆ i+1 = ∆ i+1 . Then by using Lemma 3.8 again,

σ i+1 = w -1 (∆ i+1 ) = w -1 ( ∆ i+1 ) = σ i+1 .
The induction is complete, so that the optimal speed policy is unique.

This completes the last part of the proof of Theorem 3.1.

Convergence of the Truncated MDPs

To show the convergence of the sequence of the truncated MDPs to the infinite one as N goes to infinity, we first show monotonicity properties in N . These guarantee the existence of the limit and allow us to invoke the monotone convergence theorem to show that this limit the optimality equation of the infinite MDP.

The Optimal Speed is Increasing in the Size of the State Space

The next proposition states that the optimal cost is increasing in the number of states.

Proposition 4.1. Let N ≥ 1. Then, J * N ≤ J * N +1 .
Proof. Let N ≥ 1 and let σ be the optimal speed policy when the state space is {0, • • • , N + 1}. Let σ be a speed policy for the MDP with state space {0, • • • , N }, defined as the truncation of σ: for i ∈≤ i, σ i = σ i . Recall that X σ is the continuous time Markov chain with speed policy σ on the reduced state space {0, • • • , N }, and therefore cannot be compared directly with X σ , which is defined on {0, • • • , N +1}. Thus, let X be the following discrete time Markov chain, with state space {0, • • • , N + 1} and transition probabilities given by, for i ≤ N + 1,

p i,j =                1 U (N +1) λ (N ) i if 1 ≤ j = i + 1 ≤ N 1 U (N +1) (σ i + iδ) if 0 ≤ j = i -1 ≤ N -1 1 U (N +1) Ū (N +1) i if 0 ≤ j = i ≤ N 1 U (N +1) (U -(N + 1)δ -σ N +1 ) if j = i = N + 1 0 otherwise,
with the notations:

U (N +1) := 2λ + (N + 1)δ + σ max , λ (N ) i := λ 1 -i N and Ū (N +1) i := U -λ (N ) i -δi -σ i .
With a slight abuse of notation, we denote by X σ and X σ the uniformized Markov chains with the same uniformization constant U (N +1) , so that we will be able to compare both Markov chains defined with the same time step. Moreover, notice that the Markov chain X is not irreducible: the use of the last state is only to extend the chain to a larger number of states while keeping the behavior of X similar to the one of X σ .

In order to effectively compare X and X σ , we will now define a coupling ( Y , Y ) such that Y and Y have the same distributions as X and X σ respectively. For each time step k ∈ N, let A(k) be a uniformly distributed random variable on [0, U +1) ]. Then,

Y (k + 1) =                                    Y (k) + 1 if A(k) ∈ [0, λ (N ) Y (k) ] and Y ≤ N Y (k) if A(k) ∈ [λ (N ) Y (k) , U (N +1) -Y (k)δ -σ Y (k) ] and Y (k) ≤ N Y (k) -1 if A(k) ∈ [U (N +1) -Y (k)δ -σ Y (k) , U (N +1) ]
and

Y (k) ≤ N Y (k) if A(k) ∈ [0, U (N +1) -Y (k)δ -σ Y (k) ] and Y (k) = N + 1 Y (k) -1 if A(k) ∈ [U (N +1) -Y (k)δ -σ Y (k) , U (N +1) ]
and

Y (k) = N + 1,
and similarly:

Y (k + 1) =      Y (k) + 1 if A(k) ∈ [0, λ (N +1) Y (k) ] Y (k) if A(k) ∈ [λ (N +1) Y (k) , U (N +1) -Y (k)δ -σ Y (k) ] Y (k) -1 if A(k) ∈ [U (N +1) -Y (k)δ -σ Y (k) , U (N +1) ].
By construction, if Y (0) = Y (0), then we show by induction that for all n, Y (k) ≤ Y (k).

We now check for all possible cases. For all i ≤ N + 1, λ

(N ) i < λ (N +1) i , therefore if Y (k) = Y (k), then Y (k + 1) ≤ Y (k + 1). If Y (k) = Y (k) -1 with Y (k + 1) = Y (k) + 1, then Y (k + 1) ≥ Y (k) by definition of U , as U -Y (k)δ -σ Y (k) ≥ λ ≥ λ (N ) Y (k)
, and therefore Y (k + 1) ≥ Y (k + 1). In the remaining cases, Y (k + 1) ≤ Y (k + 1).

Hence, for all k, Y (k) ≤ Y (k), so that X ≤ st X σ ; here, ≤ st denotes the stochastic order [START_REF] Shaked | Stochastic orders and their applications[END_REF]. As i → c(i, σ i ) is increasing, this implies:

E σ c( X, σ) ≤ E σ c(X σ , σ) = J * N +1
, and we have

J N ( σ) = E σ c(X σ , σ) = E σ c( X, σ).
Therefore,

J * N ≤ J N ( σ) ≤ J N +1 (σ) = J * N +1
as desired.

The next proposition states that the optimal speed policy is increasing in the size of the state space.

Proposition 4.2. Let N ≥ 1 and let σ (N ) be the unique optimal speed for the N -th MDP. Then, σ

(N ) i ≤ σ (N +1) i and ∆ (N ) i ≤ ∆ (N +1) i for all i ≤ N .
Proof. We use the expression of the optimal speed from Lemma 3.8 and the bias function to show by induction on the states i ≥ 1 that:

∆ (N ) i ≤ ∆ (N +1) i and σ (N ) i ≤ σ (N +1) i , (17)
where ∆

(N ) i = H(i) -H(i -1
) and H is the bias function for the MDP with state space [0, N ].

We first have that σ

(N ) 0 = σ (N +1) 0 = 0 and that ∆ (N ) 1 = J * N λ ≤ J * N +1 λ = ∆ (N +1) 1
, where the inequality comes from Proposition 4.1. Let us now assume that for some i ≥ 1, ∆

(N ) i ≤ ∆ (N +1) i and σ (N ) i ≤ σ (N +1) i . ∆ (N ) i ≤ ∆ (N +1) i
, we directly have that σ

(N ) i ≤ σ (N +1) i .
To show the first inequality, we write the optimality equation ( 16) for N and N + 1 with their respective optimal speed, and subtract one from the other to get:

J * N +1 -J * N = -(iδ + σ (N +1) i )∆ (N +1) i + (iδ + σ (N ) i )∆ (N ) i + λ i (∆ (N +1) i+1 -∆ (N ) i+1 ).
Using Proposition 4.1 and the induction assumption, we get

λ i (∆ (N +1) i+1 -∆ (N ) i+1 ) ≥ 0.
Then using Lemma 3.8, as ∆

(N ) i+1 ≤ ∆ (N +1) i+1
and w -1 is increasing,

σ (N +1) i+1 = w -1 (∆ (N +1) i+1 ) ≥ w -1 (∆ (N ) i+1 ) = σ (N ) i+1 .
This completes the induction so that the optimal speed policy is increasing.

Convergence Results and Proof of Theorem 2.2

For the truncated MDP in discrete time, let g * , σ (N ) and H (N ) be the optimal average cost, the optimal policy and its bias. They satisfy the optimality equation.

g * + H (N ) (i) = 1 U c(i, σ ) i ) + j H (N ) (j)p N i,j (σ (N ) i ) ∀i ≤ N.
This implies that for the truncated model in continuous time, N ) and H (N ) also satisfy the optimality equation

J N = U g * , σ (
J N = c(i, σ (N ) i ) + j H (N ) (j)q N i,j (σ (N ) i ) ∀i ≤ N.
Furthermore, for all N , J N ≤ Cλ by Lemma 3.3, σ (N ) ≤ B N ≤ B by Theorem 3.1 and since the function w is increasing, H (N ) (i) ≤ iw (B) for all i ≤ N by Lemma 3.8. Now, by the monotonicity of J N (Proposition 4.1) and the monotonicity of σ (N ) and H (N ) (Proposition 4.2), they all converge to finite non-negative limits when N goes to infinity, denoted respectively by J ∞ , σ (∞) and H (∞) .

As for the rates, q N ij (s) also converges monotonically to q ij (s) and is continuous in s. Finally, the immediate cost c is continuous in σ. The monotone convergence theorem implies that these limits satisfy an optimality equation,

J ∞ = c(i, σ (∞) i ) + j H (∞) (j)q i,j (σ (∞) i ) ∀i ∈ N.
This shows that these limits are respectively the optimal average cost (J ∞ ), an optimal policy (σ (∞) ) and its bias (H (∞) ) for the original MDP. This completes the proof of Theorem 2.2.

Cost and Deadline-Miss Probability Approximations

Our main result, Theorem 2.2, may suggest to consider the simple policy σ B defined by σ B i = B I {i>0} for all i ≥ 0, where B is defined in (3) and I A is the indicator function of the event A. Thus, σ B uses constant speed B whenever the system is busy. In this section, we numerically show that the dynamics induced by σ B are "close" to the ones induced by σ * in the sense that the optimal average cost J(σ * ) is very well approximated by J(σ B ).

Then, we investigate mathematical properties of the Markov chain induced by σ B and we give an upper bound on the stationary probability of missing deadlines as a function of the model parameters. In particular, by varying the cost parameter C, this bound can be used to keep such probabilities below a desired threshold.

In the following, we let π B denote the stationary probability of the Markov chain induced by policy σ B . Using the detailed balance equations, we obtain

π B i = π B 0 λ i i j=1 (δj + B) , i ≥ 1 and, using i≥0 π B i = 1, π B 0 = i≥0 λ i i j=1 (δj + B) -1 = δ B λ δ B δ e λ δ γ B δ , λ δ (18)
where γ(•, •) is the lower incomplete gamma function.

We will refer to the following proposition, which can be proven by summing the detailed balance equations of the underlying Markov chain as done in the proof of Lemma 3.4; for this reason, we omit the proof.

Proposition 5.1. Let X be distributed as π

B . Then, E[X] = 1 δ (λ -B(1 - π B 0 )).

Approximation of the Average Cost

By definition, the average cost induced by σ B , J(σ B ), is an upper bound on the optimal average cost, J(σ * ). In particular,

J(σ * ) ≤ J(σ B ) = i≥0 c(i, σ B ) π B i = Cδ i≥1 i π B i + i≥1 w(σ B i ) π B i = C(λ -B(1 -π B 0 )) + w(B)(1 -π B 0 )
. where the last equality follows by Lemma 5.1.

Deadline-Miss Probabilities

Let us consider the probability p M that a job misses its deadline under the stationary regime of policy σ B . This is defined by

p M := i≥1 π B i 1 -π B 0 δi δi + B . (20)
Our objective here is to control p M by fine-tuning the model parameter C.

In other words, we want to design C such that p M remains below a given threshold. Though the structure of π B 0 in (18) implies that the exact relation between C and p M is not trivial, this problem can be clearly solved numerically. Nonetheless, we aim at developing an upper bound on p M allowing for a simple analytic evaluation.

The following proposition provides a first upper bound on p M .

Proposition 5.2.

p M ≤ 1 1 -π B 0 1 - δπ B 0 δ + B - B λ . (21)
Proof. This proposition is proven by the following inequalities

p M ≤ i≥1 π B i 1 -π B 0 δ(i + 1) δ(i + 1) + B = δ λ(1 -π B 0 ) i≥1 (i + 1) λ i+1 i+1 j=1 (δj + B) π B 0 = 1 1 -π B 0 δ λ i≥1 (i + 1)π B i+1 = 1 1 -π B 0 1 - B λ (1 -π B 0 ) - δ δ + B π B 0 1 -π B 0 (22) = 1 1 -π B 0 - B λ - δ δ + B π B 0 1 -π B 0 .
In ( 22), we have used Proposition 5.1.

By coupling the underlying Markov chain under σ B with an auxiliary M/M/∞ queue with arrival rate λ and service rate δ + B, we notice that the state of the former is stochastically dominated by the latter. Therefore, π B 0 ≤ e -λ δ+B . Using that the mapping x → 1 1-x 1 -δx δ+B is increasing over [0, 1) and the previous inequality on π B 0 in (21), we obtain

p M ≤ p M = 1 1 -e -λ δ+B 1 - δe -λ δ+B δ + B - B λ . (23)
We notice that p M = 1 if B = 0 and that p M → 0 as B → ∞.

Now, the bound p M can also be used to adjust the cost of missed deadlines C so that the proportion of jobs that miss their deadline will stay below some acceptable level α.

Since Equation ( 23) cannot be inverted in closed form, a first order Taylor expansion gives p M ≈ δ δ+B when B goes to infinity. Using the value B = C/3 for the classical power dissipation w(s) = s 3 (see Remark 2.3), we get

C ≈ 3δ 2 1-p M p M 2 . Therefore, C = 3δ 2 1-α α 2
is the cost of missed deadlines that keeps the deadline-miss probability below α.

Accuracy Assessment

By means of numerical calculations, we now evaluate the accuracy of the bound in (19). For this purpose, we consider the two following set-ups. In the first one, we focus on one acceptable value for the probability to miss a deadline under policy σ B and we consider that such an acceptable level is p M = 0.1. We let the deadline rate δ and the arrival rate λ vary in [0.1, 1] and [0.1, 50], respectively. Each value of the couple (λ, δ) induces a unique value for B through (20) and for C through (3). Then, for each value of (λ, δ), we compute the percentage relative error E := J(σ B ) -J(σ * ) J(σ * ) 100% (24)

where J(σ B ) and J(σ * ) are computed numerically by truncating the state space to some N large enough so that increasing N does not change the average cost by more than 10 -4 ; in all cases, N ≤ 250.

Figure 5 depicts the level sets for the values of E. We observe that E decreases as the arrival rate increases and that it is small under a wide set of parameters. In particular, in heavy traffic conditions and with "large" deadlines the percentage relative error can be smaller than 3%. Thus, we conclude that our approximation is accurate within these conditions. On the other hand, in light load conditions and with "short" deadlines, E can be above 20% and in this case the optimal speed profile is more complex.

In the second set of experiments, we consider a fixed cost C instead of a fixed miss probability. Specifically, we let C = 300, which implies that the upper bound on the speed is B = C/3 = 10. For each pair (λ, δ), we compute the relative error E for policy σ B compared with the optimal policy, σ * . The corresponding results are displayed in Figure 6. Unlike for the case with fixed probability, the error decreases sharply as the arrival rate λ increases. In particular, the error is very small in heavy traffic (λ ≥ B = 10).

Conclusion and Perspectives

In a stochastic environment, it is well known that the problem of determining the load-dependent speed profile that a DVFS processor should follow to minimize the mean energy consumption under hard deadlines is difficult. In this paper, the basic idea is to follow a Lagrangian approach where hard deadlines are replaced by soft deadlines, meaning that jobs are allowed to miss their deadline, though each missed deadline comes with a penalty that can be fine-tuned to keep the proportion of jobs missing their deadline as small as desired. The resulting advantage of this approach stands in the ability of formulating the problem above as a Markov decision process for which we can establish constructive structural results (Theorem 2.2). Beyond the existence of monotone optimal policies, we have found that the optimal speed is bounded from above by some constant that does not depend on the deadline and the arrival rates, which may be quite surprising, and that such constant let us define an extremely simple policy whose average cost is close to the optimal one.

There are some open questions that we leave as future research. First, we have assumed that the available processing speeds vary continuously on a compact set. On the other hand, it may be convenient to consider the case where only a finite number of speeds is available. Second, it is interesting to investigate whether our results are insensitive to the job size and/or deadline distributions. This may be justified by the fact that the proposed queueing system is somewhat similar, when the speeds are bounded, to an M/M/∞ queue. Furthermore, if service times were to follow a phase-type distribution, we could use again an MDP formulation to model the problem but the analysis would require more work as the state space would be much larger. Finally, we wonder whether the optimal policy could be "learned" and, in this respect, how the proposed upper bound could speed up the learning process.

  M N with decaying arrival rates (b) M N with no decay of the arrival rates

Figure 2 :

 2 Figure 2: Optimal policies σ * and σ * for the two truncated MDPs, respectively, M N (a) and M N (b).

Figure 3 :

 3 Figure 3: Two optimal speed policies in M N , for N = 60 and N = 1000 and the bound B = C 3 .

Figure 4 :

 4 Figure 4: The transition rates of the Markov chain X * .

Figure 5 :

 5 Figure 5: Level sets of the percentage relative error E with a fixed miss probability p M = 0.1. The darker (lighter) the zone, the higher (lower) the error.

Figure 6 :

 6 Figure 6: Level sets of the percentage relative error E with a fixed cost per deadline, C = 300.

Supported by the ANR under the "Investissements d'avenir" program (ANR-15-IDEX-02) and the LabEx PERSYVAL (ANR-11-LABX-0025-01)


Viktoria Wagner
email: viktoria.wagner@ualberta.ca

Martin Večeřa
email: martinvec@seznam.cz

Borja Jiménez-Alfaro
email: jimenezalfaro@uniovi.esjanpergl

Jan Pergl
email: jan.pergl@ibot.cas.cz

Jonathan Lenoir
email: jonathan.lenoir@u-picardie.fr

Jens-Christian Svenning
email: svenning@bio.au.dk

Petr Pyšek
email: petr.pysek@ibot.cas.cz

Emiliano Agrillo
email: emiliano.agrillo@isprambiente.itidoiabiurrun

Idoia Biurrun
email: idoia.biurrun@ehu.eus

Juan Antonio Campos
email: juanan.campos@ehu.es

Jörg Ewald
email: joerg.ewald@hswt.de

Federico Fernández-González
email: federico.fdez@uclm.esute

Ute Jandt
email: ute.jandt@botanik.uni-halle.de

Valerijus Rašomavičius
email: valerijus.rasomavicius@gamtc.lt

Urban Šilc

Željko Škvorc
email: zskvorc@sumfak.hr

Kiril Vassilev

Thomas Wohlgemuth
email: thomas.wohlgemuth@wsl.ch

Milan Chytrý
email: chytry@sci.muni.cz

Borja Jiménez

Federico Fernández

Alien plant invasion hotspots and invasion debt in European woodlands

Keywords: Alien plants, biological invasions, Europe, European Vegetation Archive (EVA), exotic plants, forest, habitat type, invasibility, neophytes, non-native plants

Questions: European woodlands harbor at least 386 alien plant species but the factors driving local invasions remain unknown. By using a large vegetation-plot database,

Introduction

At least 386 alien plant species have been introduced to European woodlands deliberately or accidentally after the year 1500 [START_REF] Wagner | Alien plant invasions across European woodlands[END_REF]. While most alien species reach only low abundances and frequencies, some can become invasive by spreading fast. Some few species can even become transformative by changing plant species composition, suppressing native tree regeneration [START_REF] Langmaier | A systematic review of the impact of invasive alien plants on forest regeneration in European temperate forests[END_REF], and altering ecosystem functioning (e.g., [START_REF] Aerts | Invasion by the alien tree Prunus serotina alters ecosystem functions in a temperate deciduous forest[END_REF]. For example, black locust (Robinia pseudoacacia) can homogenize woodland understory composition [START_REF] Šibíková | Effect of planting alien Robinia pseudoacacia trees on homogenization of Central European forest vegetation[END_REF][START_REF] Vítková | Robinia pseudoacacia dominating vegetation types of Southern Europe: species composition, history, distribution and management[END_REF] and trigger soil acidification and nitrification [START_REF] Lazzaro | How ecosystems change following invasion by Robinia pseudoacacia: Insights from soil chemical properties and soil microbial, nematode, microarthropod and plant communities[END_REF]. Similarly, black cherry (Prunus serotina) can change nitrogen, phosphorus, and carbon cycles to its advantage, potentially altering the photosynthetic capacity of the long-lived native broadleaved species [START_REF] Aerts | Invasion by the alien tree Prunus serotina alters ecosystem functions in a temperate deciduous forest[END_REF].

To decrease possible detrimental effects, the European Union has outlined strategies for the prevention and management of invasive alien species (European Parliament, 2014). However, designing efficient strategies depends on understanding the drivers of local alien plant invasions across European regions. Such knowledge is particularly important for woodlands because they cover a third of the continent's territory and provide essential ecosystem services, including carbon sequestration, water retention, timber production, and biodiversity conservation (Forest Europe, 2015). However, the factors that drive alien plant invasions in European woodlands are not well understood. A recent Europe-wide study has shown that in general, macroclimate is an important abiotic predictor of local plant species richness in woodlands, in congruence with a north-west to south-east continental gradient [START_REF] Večeřa | Alpha diversity of vascular plants in European forests[END_REF]. However, it is not clear whether the same patterns holds true specifically for local alien plant species richness. Macroclimate was not an essential predictor of local alien plant species richness in deciduous woodlands of Catalonia [START_REF] Gassó | Regional context affects native and alien plant species richness across habitat types[END_REF] and was just one of several predictors of local alien plant species richness in woodlands of the Czech Republic [START_REF] Divíšek | High-resolution and large-extent mapping of plant species richness using vegetation-plot databases[END_REF]. This indicates that other predictors could be more important. For instance, there is strong evidence that the level of invasion generally decreases with elevation, partly due to a drop in human activities and alien plant propagule pressure, which are usually confounded [START_REF] Alexander | Assembly of nonnative floras along elevational gradients explained by directional ecological filtering[END_REF][START_REF] Pyšek | Colonization of high altitudes by alien plants over the last two centuries[END_REF]. By comparison, the role of abiotic soil conditions in structuring alien plant invasions is less clear. Alien plant invasions are often favored at sites with increased soil nutrient availability [START_REF] Davis | Fluctuating resources in plant communities: a general theory of invasibility[END_REF][START_REF] Gurevitch | Effects of experimental manipulation of light and nutrients on establishment of seedlings of native and invasive woody species in Long Island, NY forests[END_REF], yet many alien species can grow in low-nutrient ecosystems, likely because of their higher resource use efficiency compared to native species [START_REF] Funk | Resource-use efficiency and plant invasion in low-resource systems[END_REF]. Soil pH explains some variation in local plant species richness in European woodlands (e.g. [START_REF] Večeřa | Alpha diversity of vascular plants in European forests[END_REF] but its importance specifically for local alien plant richness is unknown. Human presence and disturbance at the landscape scale can accelerate alien plant invasions by enhancing propagule pressure and facilitating habitat edge effects [START_REF] Vilà | Plant invasions in the landscape[END_REF][START_REF] Conedera | Drivers of evergreen broadleaved spread into deciduous forest in the Southern Alps[END_REF]. Levels of invasion have been found to increase with landscape fragmentation [START_REF] Ohlemüller | Local vs regional factors as determinants of the invasibility of indigenous forest fragments by alien plant species[END_REF], proportion of urban and industrial land (Chytrý et al., 2008a;[START_REF] Liendo | Drivers of plant invasion at broad and fine scale in short temperate streams[END_REF], proximity to roads [START_REF] González-Moreno | Landscape context modulates alien plant invasion in Mediterranean forest edges[END_REF], and human population density (Pyšek et al., 2010b;[START_REF] Spear | Human population density explains alien species richness in protected areas[END_REF]but see Carboni et al., 2010).

Many studies have revealed marked differences in the level of invasion across habitat types (e.g., [START_REF] Chytrý | Invasions by alien plants in the Czech Republic: a quantitative assessment across habitats[END_REF]2008b;[START_REF] Campos | Assessing the level of plant invasion: A multi-scale approach based on vegetation plots[END_REF][START_REF] González-Moreno | Plant invasions are context-dependent: multiscale effects of climate, human activity and habitat[END_REF][START_REF] Wagner | Alien plant invasions across European woodlands[END_REF][START_REF] Giulio | Alien flora across European coastal dunes[END_REF]. Habitat types explained more variation (18%) in levels of invasion by neophytes (alien species that arrived in Europe after 1500 AD) than propagule pressure (3.1%) and climate (5.9%) in a study from the Czech Republic (Chytrý et al., 2008a). Habitat type was also the strongest predictor of alien plant invasions in Iberian vegetation [START_REF] González-Moreno | Plant invasions are context-dependent: multiscale effects of climate, human activity and habitat[END_REF]. This signal could reflect the integrating role of habitat type as a surrogate capturing the effects of both abiotic and biotic filters but might also reflect a confounding effect of the environmental variables that underlie habitat typologies, like elevation, landscape configuration, soil properties, or biotic interactions.

Identifying the main factors driving alien plant invasions is complex because invasions are dynamic processes and alien species pools not fully saturated due to dispersal limitation [START_REF] Brown | Diversity and invasibility of southern Appalachian plant communities[END_REF]. This could lead to invasion debt, i.e. a time lag between the introduction and spread of alien species [START_REF] Essl | Socioeconomic legacy yields an invasion debt[END_REF]2012). At a large scale, this could translate into hotspots of invasion among regions where introduction and spread have been frequent, and invasion debt among regions that are lagging in this respect [START_REF] Ronk | Observed and dark diversity of alien plant species in Europe: estimating future invasion risk[END_REF].

Hence, an analysis of alien species richness and abundance using quantiles closer to the maximum could provide better insights than the more conventional analysis of the mean response [START_REF] Brown | Diversity and invasibility of southern Appalachian plant communities[END_REF].

Here, we are building on our previous work on alien plant invasions in European woodlands [START_REF] Wagner | Alien plant invasions across European woodlands[END_REF] by analyzing a large dataset of woodland vegetation-plot data stored in the European Vegetation Archive [START_REF] Chytrý | European Vegetation Archive (EVA): An integrated database of European vegetation plots[END_REF]. Whereas our previous study focused on characterizing the alien species pool, this study aims to understand the possible drivers of alien plant invasion in European woodlands. Specifically, we asked: (1) How do relative alien species richness and sum of alien species covers change across regions, habitat types, environmental conditions and human disturbance? [START_REF] Brundu | Code of Conduct on Plantation Forestry and Invasive Alien Trees -T-PVS/Inf[END_REF] What is the importance of these predictors for explaining variation in alien plant invasions?

Methods

Vegetation data

We extracted vegetation-plot data from the European Vegetation Archive (EVA, version 13 January 2016, see Online Supplement and [START_REF] Wagner | Alien plant invasions across European woodlands[END_REF] for details on data preparation). This dataset contained plots from Western, Central and Southern Europe, as well as Lithuania (Fig. 1). For every vascular plant species in a plot, we assigned information on its status in the respective country (native, alien, or uncertain), with alien species defined as those that have been introduced after the year 1500 (i.e. neophytes sensu [START_REF] Thellung | Pflanzenwanderungen unter dem Einfluss des Menschen[END_REF]) and Raunkiaer's life-form (e.g., hemicryptophytes, phanerophytes, therophytes). The dataset encompassed data for plots surveyed in and after the year 1970 (53% collected after 2000)

and was restricted to plots with an area ranging from 100 to 1000 m 2 (or lacked plot size information, assuming that most of these plots were within this range; 77% of plots had a size of 100-500 m 2 ). We assigned plots to European Nature Information System (EUNIS) woodland habitat types via a custom cross-walk (see Wagner et al., 2017 for details). The used version of the EUNIS habitat classification corresponds to that used in the European Red List of Habitats [START_REF] Janssen | European Red List of habitats. Part 2. Terrestrial and freshwater habitats[END_REF]. A crosswalk to the recently revised EUNIS version is provided by [START_REF] Chytrý | EUNIS Habitat Classification: expert system, characteristic species combinations and distribution maps of European habitats[END_REF]. Only plots unambiguously assigned to a habitat type were included in the analysis.

Our study focused on the spontaneously established alien vascular plants rather than outcomes of deliberate planting. Hence, we excluded plots with alien tree and shrub species mostly known from planted populations (e.g., Abies alba in Belgium, Pinus strobus in Slovenia). We also removed plots with alien tree species that are naturalized but widely planted if they covered more than 20% in the tree layer or an unidentified vegetation layer. In addition, we also filtered out plots in which taxa with an uncertain status (native vs. alien) had more than 20% cover (for details, see [START_REF] Wagner | Alien plant invasions across European woodlands[END_REF].

Response variables

For every vegetation plot, we calculated two measures of the level of invasion: (i) relative alien species richness, defined as the percentage of the number of alien vascular plant species among all vascular plant species in a plot, and (ii) absolute sum of alien species covers (in %) in a plot (see Appendix S3, Fig. S3.1 for the frequency distribution of response variables). We chose relative alien species richness instead of the total number of alien species per plot because our dataset included different plot sizes. If an alien species occurred in different layers, we included its absolute sum across layers (alien species occurred in different layers only in 16% of plots). We also considered the presence/absence of any alien species in a plot as a metric but dropped it since subsequent statistical models performed poorly, likely due to the overabundance of absences. In addition, we computed response metrics for five alien species groups: (i) all species; (ii) all species except Impatiens parviflora; (iii) phanerophytes;

(iv) hemicryptophytes; and (v) therophytes (including I. parviflora). We removed Impatiens parviflora from one subset since it has a much higher frequency in European woodlands than any other alien species [START_REF] Wagner | Alien plant invasions across European woodlands[END_REF] and could dominate the results. We focused on phanerophytes, hemicryptophytes, and therophytes because they represent the three most frequent life-forms among alien vascular plant species in European woodlands [START_REF] Wagner | Alien plant invasions across European woodlands[END_REF].

Predictor variables

We considered 75 predictor variables for our analysis: habitat type, region, elevation, six climate variables and two Principal Component (PC) axes based on a PC analysis (PCA) of climate variables, four soil variables and two PC axes based on a PCA of these soil variables, and 50 human disturbance variables (including variables extracted for different spatial buffers) and eight PC axes based on PCAs of these human disturbance variables (two axes for each buffer category of 200m, 500m, 1000m and 2000m) (Table 1). We hypothesized that patterns in our data could vary by region, possibly due to regional differences in alien plant introduction and planting history. To test this hypothesis, we assigned plots to a grid of 250 km × 250 km, with grid cell ID as a predictor factor with 41 levels. We chose this specific resolution because it was a compromise between capturing regional effects and limiting the total number of levels in this variable. Information on elevation was derived from the original records or, if missing, extracted with the raster package [START_REF] Hijmans | raster: Geographic Data Analysis and Modeling[END_REF] in R (R Core Team, 2020) following [START_REF] Jarvis | Hole-filled SRTM for the globe. Version 4[END_REF]. We extracted climate variables that could potentially affect alien plant species distribution and abundance (BIO1: Mean annual temperature; BIO4: Temperature seasonality; BIO6: Minimum temperature of the coldest month; BIO10: Mean temperature of the warmest quarter; BIO12: Annual precipitation sum; BIO15: Precipitation seasonality) from the CHELSA dataset, which has a resolution of 30 arc seconds (~1 km × 1 km at the equator; [START_REF] Karger | Climatologies at high resolution for the earth's land surface areas[END_REF]. Topsoil properties (pH and fraction of sand, organic carbon, and gravel) were retrieved from the Harmonized World Soil Database, which offers data at a 0.05 resolution (~5.5 km × 3.9 km at the equator; [START_REF] Wieder | Regridded Harmonized World Soil Database v1.2. Data set[END_REF].

We extracted several measures to capture human disturbance. Imperviousness (the percentage of sealed soil, including roads, buildings, parking lots and industry areas) was prepared from a raster dataset of the Copernicus satellite program, which is based on an automatic derivation of calibrated Normalized Difference Vegetation Index (NDVI) in the year 2006 (European Environment Agency, 2013; 100 m × 100 m resolution). Estimates of human population density were retrieved from a dataset of the GEOSTAT initiative of the European Union (Eurostat, 2012), based on a 1 km × 1 km grid for the year 2011. We calculated the level of woodland fragmentation and woodland edge complexity based on Copernicus satellite products (European Environment Agency, 2016; original resolution: 20 m × 20 m). To quantify fragmentation, we used the area-weighted average patch size index (Sa; [START_REF] Turner | Landscape Ecology in Theory and Practice: Pattern and Process[END_REF], computed as Sa = ∑(Sk 2 ) / ∑ (Sk), where Sk is the size of the k-th woodland patch within a spatial buffer. The larger the index, the larger the patch sizes within a buffer and the less fragmentation. The complexity of woodland edges was quantified using Patton's diversity index [START_REF] Patton | A diversity index for quantifying habitat "edge[END_REF], calculated as the perimeter of an individual woodland patch divided by the perimeter of a circle covering the same area as that of the focal patch. Values > 1 indicate more complex patch edges and less interior woodland environment compared to the circular patch area. We estimated the proximity to the next propagule dispersal corridor by calculating the distance between plots and the closest road or railroad, respectively, in OpenStreetMap (OpenStreetMap Contributors, 2019). We considered all roads (motorways, trunks, primary, secondary, tertiary roads, residential streets) and all railroads (regular railroads, light rails, monorails, narrow gauges, racks, tram lines). Road density was extracted from global road density raster data (9 km × 9 km resolution; [START_REF] Meijer | Global patterns of current and future road infrastructure[END_REF]. We estimated the cover of human-altered landscapes based on CORINE land-cover types [START_REF] Bossard | CORINE Land Cover technical guide -Addendum 2000[END_REF], which we aggregated into four groups: (i) Intensive agricultural land cover, (ii) Artificial land-cover types I: Artificial land cover potentially related to the spread of alien species, (iii) Artificial land-cover types II: Theoretically the highest affinity to spread of alien species, and (iv) All human-dominated land-cover types combined (see Appendix S2, Table S2.1 for more details). We hypothesized that human disturbance has an effect at the landscape level by influencing the area surrounding a plot. Hence, we used spatial buffers of 200 m, 500 m, 1000 m, and 2000 m when extracting mean values for human population density, imperviousness and road density (using QGIS; QGIS Development Team, 2019), and woodland fragmentation, woodland edge complexity, and the cover of human-dominated land-cover types (using ArgGIS; ESRI, 2011). For comparison purposes, we also extracted maximum values for imperviousness, human population density, and road density within the buffers in QGIS.

In addition, we constructed simplified and uncorrelated variables for three predictor groups (climate, soil, human disturbance) by calculating two ordination axes using PCA in the vegan package [START_REF] Oksanen | vegan: Community Ecology Package. Version 2.5-3[END_REF] in R (R Core Team, 2020; see details in Appendix S1). These single predictor variables and the PC axes were subsequently used as predictor variables in our bivariate analyses.

Data filtering

After preparing our response and predictor variables, we took several data filtering steps.

First, we included only plots with geographic coordinates and a location uncertainty of ≤ 1000 m (44% of plots in the final dataset had an uncertainty of ≤ 500 m). However, plots that lacked information on location uncertainty but were surveyed after the year 2000 were included (33%), assuming that by that time, global positioning system (GPS) devices or accurate georeferencing were used to record plot locations. We reduced spatial clustering and sampling bias at the landscape and regional scales by randomly resampling up to 10 plots within a combination of 10 km × 10 km grid cell and habitat type. Next, to balance out regional unevenness, we randomly selected up to 200 plots for a combination of 250 km × 250 km grid cell and habitat type. Regional grid cells with < 100 plots were excluded from the analysis. Fagus woodlands on non-acid soils were much more common than other habitat types. To balance out this overrepresentation, we randomly resampled up to 200 plots of this habitat type per country. We finally excluded rare habitat types with < 400 plots in the whole dataset after resampling and deleted outliers (see Appendix S1 in Supplementary Information). The final dataset included 16,211 plots (Fig. 1) and encompassed 15 EUNIS habitat types (see Appendix S2, Tables S2.2, S2.3 and S2.4 for an overview of databases, regions and habitats represented in the final dataset).

Bivariate relationships

We fitted standard generalized additive models (GAMs; [START_REF] Wood | Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models[END_REF] and quantile GAMs (qGAMs; [START_REF] Fasiolo | Fast calibrated additive quantile regression[END_REF] to bivariate relationships between each of the two studied response variables on the one hand, and predictor variables on the other hand (both model types were run with a Gaussian error family). GAMs and qGAMs allowed us to explore trends using the original, untransformed variables. We used GAMs to inspect trends in average response rate and qGAMs to analyze values closer to the maximum response rates (𝜏 = 0.90, and 0.99). We included the qGAMs approach because our data were zero-inflated (alien plants were present only in 17% of plots), and we hypothesized that signals could be fitted better at higher quantiles rather than as averages. We ran GAMs and qGAMs using the R packages mgcv [START_REF] Wood | Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models[END_REF] and qgam [START_REF] Fasiolo | Fast calibrated additive quantile regression[END_REF], respectively. Smoothing parameters were constructed using cubic regression splines with k = 10 dimensions and estimated with restricted maximum likelihood estimation.

Variable importance

We used Random Forest (RF) analyses in the randomForest package [START_REF] Liaw | Classification and regression by randomForest[END_REF] in R to assess the importance of predictor variables for structuring relative alien species richness and the sum of alien species covers. RF analysis is a machine learning algorithm that averages predictions over several recursively partitioned decision trees based on bootstrapped data ("bagging"), creating a higher accuracy than decision trees. It uses a random subset of predictors for each tree to assess predictor importance in the presence of multicollinearity [START_REF] Breiman | Random forests[END_REF]. Nonetheless, collinearity was present when we included all predictor variables, as suggested by generalized variance inflation factors (GVIF, Appendix S2, Table S2.6), which could bias RF results [START_REF] Dormann | Collinearity: a review of methods to deal with it and a simulation study evaluating their performance[END_REF]. Hence, we considered all predictor variables but only included those with a GVIF ≤ 2 in the RFs.

Each RF was based on 10,000 grown trees, with four randomly chosen predictors considered at each split, minimum node sizes of 40 plots, and randomization of cases without replacement. Variable importance was assessed based on out of bag (OOB) data, as the average increase in squared OOB residuals when a variable in the OOB data is permuted (%IncMSE). The latter metric is unbiased compared to the increase in node purity, another metric reported by the package (Grömping, 2009 and references therein). To compare the importance of predictor variables for our two response variables, we normalized %IncMSE values across predictors to sum to 100% [START_REF] Grömping | Variable Importance Assessment in Regression: Linear Regression versus Random Forest[END_REF].

Results

Bivariate relationships

Both the relative alien species richness and the sum of alien species covers differed by habitat type (Fig. 2). Differences were apparent in GAM estimates but amplified when using qGAMs at the maximum response rates (𝜏 = 0.90 and 0.99). According to GAM estimates, relative alien species richness was highest in temperate and boreal mountain Betula and Populus tremula woodland on mineral soils (G1.9a), acidophilous Quercus woodland (G1.8), Fagus woodland on acid soils (G1.6b), and temperate and boreal softwood riparian woodland (G1.1) (Fig. 2a). In general, habitats with low relative alien species richness had also low sum of alien species covers. However, the ranking based on the sum cover was more nuanced for habitats with higher relative alien species richness, with Temperate and boreal softwood riparian woodland (G1.1 ) having the highest sum of alien species covers (Fig. 2b). In general, woodland habitat types at temperate latitudes, such as those listed above, had higher levels of invasion than Mediterranean or montane habitat types like Mediterranean evergreen Quercus woodland (G2.1) and Temperate and submediterranean montane Pinus sylvestris-Pinus nigra woodland (G3.4b).

The relative alien species richness and the sum of alien species covers differed across regional grid cells (Fig. 3). In general, it was higher in northern than in southern regional grid cells and peaked in north-western and central Europe. These regional patterns were consistent for relative alien species richness and the sum of alien species covers and for responses at the mean and 𝜏 = 0.99 quantiles (Fig. 3).

Relative alien species richness and the sum of alien species covers declined with increasing elevation, distance to the nearest road, and distance to the nearest railroad and increased with maximum imperviousness (percentage of sealed soil cover) in 2000m buffers (Fig. 4).

However, these relationships were only pronounced when values closer to the maximum were analyzed (𝜏 = 0.99 and 0.90). None or only weak relationships were found in the analysis of average trends (𝜏 = 0.5; Fig. 4). With respect to elevation, relative alien species richness declined at elevations 0-750 m, whereas the sum of alien species covers showed a slight increase and then a decrease in this elevation range (Fig. 4). Neither the relative alien species richness nor the sum of alien species covers showed a clear trend with climate and soil variables but displayed only local maxima along these gradients (Appendix S3, Figs. S3.2 and

S3.3).

A closer inspection of bivariate relationships between relative alien species richness and the sum of alien species covers on the one hand and elevation, distance to the nearest road, distance to the nearest railroad, and the mean annual temperature on the other, revealed that peaks were driven by regional signals. For all of these predictor variables, most plots that matched the maxima in relative alien species richness at 𝜏 = 0.99 and 0.90 quantiles were from the Netherlands, Belgium, Poland and the Czech Republic (Table 2). Plots that were associated with maxima in the sum of alien species covers at these quantiles were from the Czech Republic, the Netherlands, and France (Table 3).

Variable importance

The variables that we considered in our RFs models explained together 39.6% of the total variation in relative alien species richness. By comparison, the total amount of variation explained for the sum of alien species covers was only 20.9% (Fig. 5a,b). Habitat type and regional grid cells were the most important predictors (> 10%) for both response variables.

Variable importance and percentage of explained variation changed little when we removed

Impatiens parviflora (n = 15,214 plots;Appendix S3,Fig. S3.4a,b) or inspected patterns only among alien phanerophytes (Fig. S3.4c,d). By contrast, when we analyzed only alien hemicryptophytes (Fig. S3.4e,f) or alien therophytes (Fig. S3.4 g,h), the percentage of explained variation decreased.

Discussion

Several key findings suggest that European woodlands display regional hotspots of local alien plant invasion and that woodlands in some regions might not be fully saturated with alien species but rather experience invasion debt [START_REF] Brown | Diversity and invasibility of southern Appalachian plant communities[END_REF][START_REF] Essl | Ancient and recent alien species in temperate forests: steady state and time lags[END_REF]. First, although macroclimate is an important predictor of local plant species richness in European woodlands [START_REF] Večeřa | Alpha diversity of vascular plants in European forests[END_REF], it did not emerge as a significant driver of local relative alien species richness across European woodlands. Second, local levels of alien plant invasions differed strongly by region, with regional hotspots in north-western and central Europe. Third, for continuous predictors that showed relationships with relative alien species richness and the sum of alien species covers (i.e. elevation, distance to the nearest road or nearest railroad, imperviousness), effects were only pronounced at higher quantiles and non-detectable at the average response of our response variables. When we looked closer at peaks in these relationships, the signal at the maximum quantiles was driven by vegetation plots located within regional invasion hotspots.

Interestingly, the regional hotspots that we detected broadly match hotspots of regional alien species richness described by [START_REF] Ronk | Observed and dark diversity of alien plant species in Europe: estimating future invasion risk[END_REF]. In contrast to our work, this study assessed richness at a coarser spatial resolution, based on species lists for 50 km × 50 km grid cells in Europe. This may imply that the magnitude of alien plant invasions at the local scale in European woodlands is influenced by regional differences in alien species pool size. In theory, the high estimates for the Czech Republic and France could be due to more intense sampling in these regions (% of plots in the final dataset: 15.1%, 12.1%, respectively). Plots in the Netherlands and Belgium made up only 5.9% and 2.2% of the final dataset, respectively (Appendix S2, Table S2.3), albeit these countries are also smaller. Here, the high levels of invasion could also reflect the region's legacy as an invasion hub, especially the long history of international trade at its seaports and high density of traffic corridors (Eurostat, 2018). Another explanation is the lower richness in the native woodland flora in northwestern Europe [START_REF] Večeřa | Alpha diversity of vascular plants in European forests[END_REF]. Since we used a relative measure of alien species richness (to account for different plot sizes), calculated as the number of alien species divided by all species in a plot, naturally species-poor forests could appear to be relatively more invaded. The lower magnitude of alien plant invasions in southern Europe is in line with patterns observed at the higher spatial scale by [START_REF] Ronk | Observed and dark diversity of alien plant species in Europe: estimating future invasion risk[END_REF]. The prevalence of summer drought in this region could make its woodlands more resistant to invasion. Alternatively, given the high rates of alien plant invasions in other regions with summer drought, like California (van Kleunen et al., 2015), our results could also imply that invasion debt is more prominent in southern European woodlands [START_REF] Rouget | Invasion debtquantifying future biological invasions[END_REF]. These hypotheses deserve more attention by researchers in the future. Furthermore, our study did not cover all of Europe and might have missed other regional patterns of local alien plant invasions in European woodlands.

Our results also join an increasing number of studies that found habitat type to be an important predictor of the magnitude of invasions (Chytrý et al., 2008a;Pyšek et al., 2010a;[START_REF] Pyšek | Habitat invasion research: where vegetation science and invasion ecology meet[END_REF][START_REF] Dyderski | Context-dependence of urban forest vegetation invasion level and alien species' ecological success[END_REF]. Our study encompassed a smaller set of habitat types than [START_REF] Wagner | Alien plant invasions across European woodlands[END_REF] due to stronger data filtering. However, the trends in mean relative alien species richness among habitat types were broadly similar. For instance, the three highest-ranked habitats according to the mean (G1.1, G1.8, G1.9a), as shown by GAMs, were among the five highest-ranked habitat types in [START_REF] Wagner | Alien plant invasions across European woodlands[END_REF].

The high levels of invasion in temperate and boreal softwood riparian woodlands are well known [START_REF] Schnitzler | Examining native and exotic species diversity in European riparian forests[END_REF]. They can be linked to several habitat characteristics, including a high propagule pressure at low elevations, their role as dispersal corridors and intersections by roads and railroads, and their exposure to nutrient pulses, frequent disturbance, and light gaps due to flooding. The high rankings for Acidophilous Quercus woodland and Temperate and boreal mountain Betula and Populus tremula woodland on mineral soils are more surprising. In theory, they could be explained by native species richness per plot, which tends to be low in these habitat types (Wagner et al., 2017, their Appendix S5). However, in general, we did not find a correlation between mean alien species richness and mean native species richness per habitat type (Pearson's r = -0.41, d.f. = 13, p = 0.126; calculated for the most frequent plot size of 100 m 2 ). In addition, the high levels of invasion in these habitat types may be due to their relatively open canopy cover, occurrence in regional invasion hotspots, higher local disturbance through management, and the ability of frequent alien species to colonize acidic soils and forest gaps, as displayed for example by Prunus serotina [START_REF] Godefroid | Ecological factors controlling the abundance of non-native invasive black cherry (Prunus serotina) in deciduous forest understory in Belgium[END_REF] or Impatiens parviflora [START_REF] Chmura | Autoecology of Impatiens parviflora DC. in natural forest communities[END_REF].

The magnitude of alien plant invasions in European woodlands decreased with increasing elevation. Similar signals at the local scale have been explained by declining propagule pressure or directional filtering. The latter process refers to an overrepresentation of taxa with a niche optimum at lower elevations, within the low elevation species pool (e.g., [START_REF] Alexander | Assembly of nonnative floras along elevational gradients explained by directional ecological filtering[END_REF][START_REF] Pyšek | Colonization of high altitudes by alien plants over the last two centuries[END_REF]. However, in our study, the decline of the levels of alien plant invasions with elevation could also be partly due to the confounding effect of region, with maxima reported from regional hotspots at low elevations, like the Netherlands.

Distances to the nearest road and railroad and cover of impervious surfaces were the only clear indicators of human-induced disturbance associated with the magnitude of alien plant invasions in European woodlands. In general, the decline of alien species richness and cover with increasing distances from these corridors supports the role of road and railroad traffic as a pathway for alien plant invasions [START_REF] Von Der Lippe | Do cities export biodiversity? Traffic as dispersal vector across urban-rural gradients[END_REF]. Imperviousness showed a positive relationship with our response variable, but only when we extracted maximum imperviousness values in a 2000-m buffer around the plot. This could reflect the importance of urban areas as propagule sources and as dispersal corridors [START_REF] Klotz | Urbanisation and alien invasion[END_REF] in the wider perimeter of woodlands. In our study, the effect of distances to the nearest road and railroad, and imperviousness, may also reflect a confounding effect with region.

North-western Europe, where the magnitude of invasion was particularly high, has one of the highest road and imperviousness densities in Europe (Eurostat, 2018;European Environmental Agency, 2021). Surprisingly, our study did not detect a link between alien plant invasions and other human disturbance predictors, including population density, woodland fragmentation and edge complexity, and the cover of human-dominated land. We hypothesize that other factors, like introduction and planting history, could over-ride the importance of human disturbance at the landscape scale. For instance, given that trees and shrubs are over-represented in the European alien woodland flora compared to natives [START_REF] Wagner | Alien plant invasions across European woodlands[END_REF], our finding could echo the large-scale planting of this group in the mid-18 th to early 20 th century [START_REF] Bucharová | Introduction history and species characteristics partly explain naturalization success of North American woody species in Europe[END_REF][START_REF] Nyssen | The history of introduced tree species in Europe in a nutshell[END_REF]. Well-known examples are the planting of Prunus serotina [START_REF] Schrader | Risk analysis for alien plants in European forests, illustrated by the example of Prunus serotina[END_REF][START_REF] Nyssen | The history of introduced tree species in Europe in a nutshell[END_REF] and Robinia pseudoacacia [START_REF] Vítková | Black locust (Robinia pseudoacacia) beloved and despised: A story of an invasive tree in Central Europe[END_REF] in north-western and central Europe. This hypothesis is supported by [START_REF] Bucharová | Introduction history and species characteristics partly explain naturalization success of North American woody species in Europe[END_REF] and [START_REF] Pyšek | Planting intensity, residence time, and species traits determine invasion success of alien woody species[END_REF], who found planting frequency to be the most important predictor of the success of alien trees in Europe. Planting could have removed landscape filters and facilitated the spread of alien phanerophytes into European woodlands. Furthermore, our proxies of human disturbance at the continental scale were available only at the turn of the 21 st century, roughly when 50% of our vegetation plots have been sampled. However, invasions may be the result of much earlier human disturbance activities. Linking our data to historic disturbance and planting data could have explained more variation, but such data are more difficult to gather at the European extent.

Soil properties explained only a small amount of variation in our data. It is possible that their spatial resolution (~ 5.5 × 3.9 km) was too coarse to detect any signal in the data. Linking our data to soil information at a finer spatial resolution might have improved models, but such data were not available and are generally underrepresented or recorded inconsistently in vegetation databases. Future studies that couple vegetation surveys with detailed analyses of soil abiotic and biotic properties could shed more light in this respect.

RFs predicted 39.6% and 20.9% of the variation in relative alien species richness and the sum of alien species covers, respectively. The lower variance explained in the latter response variable is likely due a higher imprecision in cover estimates compared to species richness estimates, and to a stronger influence of effects not captured in our models, such as local disturbance [START_REF] Chabrerie | Disentangling relationships between habitat conditions, disturbance history, plant diversity, and American black cherry (Prunus serotina Ehrh.) invasion in a European temperate forest[END_REF]. In addition, we considered only single effects in driving alien plant invasions. However, invasion patterns are driven by interactions of multiple effects [START_REF] Thuiller | Interactions between environment, species traits, and human uses describe patterns of plant invasions[END_REF][START_REF] Pyšek | Naturalization of central European plants in North America: species traits, habitats, propagule pressure, residence time[END_REF]. An analysis of interactions between habitat type and human disturbance, in particular, could reveal whether some habitat types are more prone to invasions under increased human disturbance.

Conclusions

European woodlands show clear signs of invasion debt for alien plants, with regional hotspots of invasions, links to predictors representing propagule pathways, and the strongest environmental predictions for maximum levels of invasion. Looking into the future, our results imply that alien plant invasions will continue to expand in European woodlands, likely driven by the complex effects of regional introduction history, factors associated with differences between habitat types, and the prevalence of dispersal corridors and urban areas.

A challenge for future research is to clarify the role of local and historic factors as drivers of alien plant invasions. Avoidance strategies, such as regulations of alien plant introductions and plantings [START_REF] Brundu | Code of Conduct on Plantation Forestry and Invasive Alien Trees -T-PVS/Inf[END_REF][START_REF] Pötzelsberger | Mapping the patchy legislative landscape of non-native tree species in Europe[END_REF], must be key components for the protection of European woodlands from future invasions. Given that alien plant invasions are dynamic processes and will likely proceed further, we advocate for continuous surveys in European woodlands at the local and regional scales.

Fig. 1 .

 1 Fig. 1. Location and relative alien species richness of vegetation plots included in the final

Fig. 2 .

 2 Fig. 2. Differences in (a) the relative alien species richness and (b) the sum of alien species

Fig. 3 .

 3 Fig. 3. Differences in relative alien species richness (a, b) and the sum of alien species covers

Fig. 4 .

 4 Fig. 4. Bivariate relationships between (a-d) relative alien species richness and (e-h) sum of

Fig. 5 .

 5 Fig. 5. Importance of predictors for explaining (a) relative alien species richness and (b) sum

Acknowledgements

We would like to thank all persons involved in collecting, curating and checking the data that were mined for this project, especially Christian Berg, Henry Brisse, Andraž Čarni, Laura Casella, János Csiky, Mirjana Ćuk, Renata Ćušterevska, Els De Bie, Michele De Sanctis, Panayotis Dimopoulos, Anton Drescher, Tetiana Dziuba, Úna FitzPatrick, Xavier Font, Valentin Golub, Stephan Hennekens, Dana Holubová, Adrian Indreica, Florian Jansen, Zygmunt Kącki, Ilona Knollová, Flavia Landucci, Vassiliy Martynenko, Viktor Onyshchenko, Tomáš Peterka, Vadim Prokhorov, John Rodwell, Solvita Rūsiņa, Joop Schaminée, Ioannis Tsiripidis, Milan Valachovič, Roberto Venanzoni, Risto Virtanen, Lynda Weekes and Wolfgang Willner. We thank Jan Divíšek for helpful comments on the methods and preparation of predictor variables.

Funding information VW initiated this study as a part of the InvasEVe project financed by the SoMoPro II programme and considers it as part of her Discovery Grant from the Canadian Natural Science and Engineering Research Council. VW's research leading to these results has acquired a financial grant from the People Programme (Marie Curie Action) of the Seventh Framework Programme of the EU according to REA Grant Agreement No. 291782, and it was further co-financed by the South-Moravian Region. MV and MC were supported by the Czech Science Foundation (19-28491X). PP and JP were supported by the Czech Science Foundation (19-28807X) and The Czech Academy of Sciences (RVO 67985939). JCS considers this work a contribution to his VILLUM Investigator project "Biodiversity Dynamics in a Changing World" funded by VILLUM FONDEN (grant 16549). IB and JAC were funded by the Basque Government (IT936-16). FFG was funded by the Plan Propio of the UCLM (2020-GRIN-29214).

Data availability statement

Vegetation-plot data used in this study are stored in European Vegetation Archive with reference to project no. 1 (2014-05-31). Accompanying data and R scripts can be accessed at 10.6084/m9.figshare.13697329.

Author contributions

VW and MC conceived the idea and developed the study design. VW and MV prepared the data. VW carried out the statistical analysis, with feedback from co-authors, and led the writing. All co-authors supplied data, species status and life-form assignments, and helped develop the manuscript.

SUPPORTING INFORMATION

Additional supporting information may be found online in the Supporting Information section.

Appendix S1. Methods supplement

Appendix S2. Supplementary tables

Appendix. S3. Supplementary figures