Joke Blom
email: joke.blom@cwi.nl

Koen De Turck
email: kdeturck@telin.ugent.be

Michel Mandjes •

M Mandjes
email: m.r.h.mandjes@uva.nl

Analysis of Markov-modulated infinite-server queues in the central-limit regime

Keywords:

come

Introduction

The infinite-server queue has been intensively studied, perhaps owing to its wide applicability and attractive computational features. In these systems jobs arrive according to a given arrival process, go into service immediately, are served in parallel, and leave when their service is completed.

An important feature of this model is that there is no waiting: jobs do not interfere with each other. The infinite-server queue was originally developed to analyze the probabilistic properties of the number of calls in progress in a trunk in a communication network, as an approximation of the corresponding system with many servers. More recently, however, various other application domains have been identified, such as road traffic [START_REF] Van Woensel | Modeling traffic flows with queueing models: a review[END_REF] and biology [START_REF] Schwabe | Transcription stochasticity of complex gene regulation models[END_REF].

In the most standard variant of the infinite-server model, known as the M/M/∞ model, jobs arrive according to a Poisson process with a fixed rate λ, where the service times are i.i.d. samples from an exponential distribution with mean µ -1 (independent of the job arrival process). A classical result states that the stationary number of jobs in the system has a Poisson distribution with mean λ/µ. Also the transient behavior of this queueing system is well understood.

In many practical situations, however, the assumptions underlying the standard infinite-server model are not valid. The arrivals often tend to be 'clustered' (so that the assumption of a fixed arrival rate does not apply), while also the service distribution may vary over time. This explains the interest in Markov-modulated infinite-server queues, so as to incorporate 'burstiness' into the queue's input process. In such queues, the input process is modulated by a finite-state (of dimension d ∈ N) irreducible continuous-time Markov process (J(t)) t∈R , often referred to as the background process or modulating process, with transition rate matrix Q ≡ (q ij ) d i,j=1 . If J(t) is in state i, the arrival process is (locally) a Poisson process with rate λ i and the service times are exponential with mean µ -1 i (while the obvious independence assumptions are assumed to be fulfilled).

The Markov-modulated infinite-server queue has attracted some attention over the past decades (but the number of papers on this type of system is relatively modest, compared to the vast literature on Markov-modulated single-server queues). The main focus in the literature so far has been on characterizing the steady-state number of jobs in the system; see e.g. [START_REF] Auria | M/M/∞ queues in semi-Markovian random environment[END_REF][START_REF] Falin | The M/M/∞ queue in a random environment[END_REF][START_REF] Fralix | An infinite-server queue influenced by a semi-Markovian environment[END_REF][START_REF] Keilson | The matrix M/M/∞ system: retrial models and Markov modulated sources[END_REF][START_REF] O'cinneide | The M/M/∞ queue in a random environment[END_REF] and references therein. Interestingly, there are hardly any explicit results on the probability distribution of the (transient or stationary) number of jobs present: the results are in terms of recursive schemes to determine all moments, and implicit characterizations of the probability generating function.

An idea to obtain more explicit results for the distribution of the number of jobs in the system, is by applying specific time-scalings. In [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF][START_REF] Hellings | Semi-Markov-modulated infinite-server queues: approximations by time-scaling[END_REF] a time-scaling is studied in which the transitions of the background process occur at a faster rate than the Poisson arrivals. As a consequence, the limiting input process becomes essentially Poisson (with an arrival rate being an average of the λ i s); a similar property applies for the service times. Under this scaling, one gets in the limit the Poisson distribution for the stationary number of jobs present. Recently, related transient results have been obtained as well, under specific scalings of the arrival rates and transition times of the background process [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF][START_REF] Blom | Time-scaling limits for Markov-modulated infinite-server queues[END_REF].

Contribution. Our work considers a time-scaling featuring in [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF][START_REF] Blom | Time-scaling limits for Markov-modulated infinite-server queues[END_REF] as well. In this scaling, the arrival rates λ i are inflated by a factor N , while the background process (J(t)) t∈R is sped up by a factor N α , for some α ∈ (0, ∞). The primary focus is on the regime in which N grows large.

The object of study is the number of jobs in the scaled system at time t, in the sequel denoted by M (N ) (t). More specifically, we aim at deriving a central limit theorem (clt) for M (N ) (t), as well as for its stationary counterpart M (N ) . Interestingly, we find different scaling regimes, based on the value of α. The rationale behind these different regimes lies in the fact that for α > 1 the variances of M (N ) (t) and M (N ) grow essentially linearly in N , while for α < 1 they grow as N 2-α .

It is important to notice that there are actually two variants of this Markov-modulated infiniteserver queue. In the first (to be referred to as 'Model i') the service times of jobs present at time t are subject to a hazard rate that is determined by the state J(t) of the background process at time t. In the second variant (referred to as 'Model ii') the service times are determined by the state of the modulating process at the job's arrival epoch (and hence can be sampled upon arrival).

The main contribution of our work is that we develop a unified approach to prove the clt s for both Model i and Model ii for the scalings given above, for arbitrary α ∈ (0, ∞), and for both the transient and stationary regimes. The technique used can be summarized as follows. We first derive differential equations for the probability generating functions (pgf s) of the transient number of jobs in the system M (N ) (t) as well as its stationary counterpart M (N ) (for both models). The next step is to establish laws of large numbers: we identify (t) ( , respectively) to which

N -1 M (N ) (t) (N -1 M (N ) , respectively) converges as N → ∞. This result indicates how M (N ) (t)
and M (N ) should be centered in a clt. The thus obtained centered random variables are then normalized (that is, divided by N γ , for an appropriately chosen γ), so as to obtain the clt. As suggested by the asymptotic behavior of the variance of M (N ) (t) and M (N ) , as we pointed out above, the appropriate choice of the parameter γ in the normalization is γ = 1 2 for α > 1, and γ = 1 -α 2 for α < 1. The proofs rely on (non-trivial) manipulations of the differential equations that underly the pgf s. For α < 1 the deviation matrix [START_REF] Coolen-Schrijner | The deviation matrix of a continuous-time Markov chain[END_REF] appears in the clt in the expression for the variance.

Relation to previous work. In our preliminary conference paper [START_REF] Blom | A central limit theorem for Markov-modulated infinite-server queues[END_REF] we just covered Model i, with an approach similar to the one featuring in the present paper. In [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF] the transient regime of Model ii is analyzed, but just for α > 1, relying on a different and more elaborate methodology. New results of this paper are: (i) Model ii for α ≤ 1, (ii) the clt for the stationary number of jobs M (N ) in Model ii, (iii) results on the correlation across time. The main contribution, however, concerns the unified approach: where earlier work has been using ad hoc solutions for the scenario at hand, we now have a general 'recipe' to derive clt s of this kind. Current work in progress aims at functional versions of the clt s for the process (M (N ) (t)) t∈R ; [START_REF] Anderson | A functional central limit theorem for a Markov-modulated infinite-server queue[END_REF] covers the special case of uniform service rates, which constitutes the intersection of Model i and Model ii.

Organization. The organization of this paper is as follows. In Section 2, we explain the model in detail and introduce the notations used throughout the paper. Section 3 provides a systematic explanation of our technique for proving this kind of clt s; in addition, we demonstrate the approach for a special case, viz. the transient analysis for the model with uniform service rates (in which Model i and Model ii coincide). In Section 4 we recall the results for Model i as derived in the precursor paper [START_REF] Blom | A central limit theorem for Markov-modulated infinite-server queues[END_REF]. Then in Section 5, we state and prove for Model ii the clt s, both for the stationary and transient distribution. The single-dimensional convergence can be extended to convergence of the finite-dimensional distributions (viz. at different points in time); see Section 6. In Section 7, we provide some numerical examples so as to get insight into the speed of convergence to the various limiting regimes. The final section of the paper, Section 8, contains a discussion and concluding remarks.

Model description and preliminaries

In this section, we first provide a detailed model description. We then give a number of explicit calculations for the mean and variance of M (N ) (t), that indicate how this random variable should be centered and normalized so as to obtain a clt. We conclude by presenting a number of preliminary results (e.g., a number of standard results on deviation matrices).

2.1.

Model description, scaling. The main objective of this paper is to study an infinite-server queue with Markov-modulated Poisson arrivals and exponential service times. In full detail, the model is described as follows.

Model. Consider an irreducible continuous-time Markov process (J(t)) t∈R on a finite state space {1, . . . , d}, with d ∈ N. Let its transition rate matrix be given by Q ≡ (q ij ) d i,j=1 ; here the rates q ij are nonnegative if i = j, whereas q ii =j =i q ij (so that the row sums are 0). Let π i be the stationary probability that the background process is in state i, for i = 1, . . . , d; due to the irreducibility assumption there is a unique stationary distribution. The time spent in state i (often referred to as the transition time) has an exponential distribution with mean 1/q i , where q i := -q ii . Let M (t) denote the number of jobs in the system at time t, and M its steady-state counterpart.

The dynamics of the process (M (t)) t∈R can be described as follows. While the process (J(t)) t∈R , usually referred to as the background process or modulating process, is in state i ∈ {1, . . . , d}, jobs arrive at the queue according to a Poisson process with rate λ i ≥ 0. The service times are assumed to be exponentially distributed with rate µ i , however, more importantly this statement can be interpreted in two ways: Model i: In the first variant of our model, the service times of all jobs present at a certain time instant t are subject to a hazard rate determined by the state J(t) of background chain at time t, regardless of when they arrived. Informally, if the system is in state i, then the probability of an arbitrary job leaving the system in the next ∆t time units is µ i ∆t.

Model ii: In the second variant the service rate is determined by the background state as seen by the job upon its arrival. If the background process was in state i, the service time is sampled from an exponential distribution with mean µ -1 i .

The difference between the two models is nicely illustrated by the following alternative representation [START_REF] Auria | M/M/∞ queues in semi-Markovian random environment[END_REF]. In Model i M (t) has a Poisson distribution with random parameter ψ(J), while in Model

ii it is Poisson with random parameter ϕ(J), where J ≡ (J(s)) s∈[0,t] , and

(1)

ψ(f ) := t 0 λ f (s) e -t s µ f (r) dr ds, ϕ(f ) := t 0 λ f (s) e -µ f (s) (t-s) ds,
with f : [0, t] → {1, . . . , d}.

Scaling. In this paper, we consider a scaling in which both (i) the arrival process, and (ii) the background process are sped up, at a possibly distinct rate. More specifically, the arrival rates are scaled linearly, that is, as λ i → N λ i , whereas the background chain is scaled as q ij → N α q ij , for some positive α. We call the resulting process (M (N ) (t)) t∈R , to stress the dependence on the scaling parameter N ; the corresponding background process is denoted by (J (N ) (t)) t∈R .

The main objective of this paper is the derivation of clt s for the number of jobs in the system, as N grows large. As mentioned in the introduction, the parameter α plays an important role here:

it turns out to matter whether α is assumed smaller than, equal to, or larger than 1. Letting the system start off empty at time 0, we consider the number of jobs present at time t, denoted by M (N ) (t); we write M (N ) for its stationary counterpart.

Our main result is a 'non-standard clt': for a deterministic function (t),

(2)

M (N ) (t) -N (t) N γ
converges in distribution to a zero-mean Normal distribution with a certain variance, say, σ 2 (t).

It is important to note that in the case α > 1 we have that the parameter γ equals the usual 1 2 , while for α ≤ 1 it has the uncommon value 1 -α 2 . A similar dichotomy holds for the stationary counterpart M (N ) . In the next subsection, we present explicit calculations for the mean and variance of M (N ) (t) and M (N ) that explain the reason behind this dichotomy.

2.2. Explicit calculations for the mean and variance. We now present a number of explicit calculations for the mean and variance of the number of jobs present; for ease we consider the case that µ i = µ for all i ∈ {1, . . . , d}, so that Models i and ii coincide. We assume J(0) is distributed according to the stationary distribution of the Markov chain J(t). Directly from, e.g., [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF], for any

N ∈ N, EM (N ) (t) N = (t) := 1 -e -µt µ d i=1 π i λ i , EM (N ) N = := 1 µ d i=1 π i λ i .
We now concentrate on the corresponding variance; we first consider the non-scaled system, to later explore the effect of the time-scaling. In the sequel we use the notation p ij (t) := P(J(t) = j | J(0) = i). The 'law of total variance', with J ≡ (J(s)) t s=0 , entails that

(3) Var M (t) = E Var(M (t) | J) + Var E(M (t) | J).
We first recall from (1) that M (t) obeys a Poisson distribution with the random parameter ϕ(J).

As a result, the second term on the right of (3) can be written as

Varϕ(J) = Var t 0 λ J(s) e -µ (t-s) ds = t 0 t 0 Cov λ J(u) e -µ (t-u) λ J(v) e -µ (t-v) dudv,
which can be decomposed into I 1 + I 2 , where

I 1 := d i=1 d j=1 λ i λ j K ij , with K ij := t 0 v 0 e -µ(t-u) e -µ(t-v) π i (p ij (v -u) -π j ) dudv, I 2 := d i=1 d j=1 λ i λ j L ij , with L ij := t 0 t v e -µ(t-u) e -µ(t-v) π j (p ji (u -v) -π i ) dudv.
Let us first evaluate K ij . To this end, substitute w := v -u (i.e., replace u by v -w), and then interchange the order of integration, so as to obtain

K ij = e -µt π i t 0 t w e 2µv dv e -µ(t+w) (p ij (w) -π j ) dw.
Performing the inner integral (i.e., the one over v) leads to

K ij = 1 2µ e -µt π i t 0 e µ(t-w) -e -µ(t-w) (p ij (w) -π j ) dw.
The integral L ij can be evaluated similarly:

L ij = e -µt π j t 0 t-w 0 e 2µv dv e -µ(t-w) (p ji (w) -π i ) dw = 1 2µ e -µt π j t 0 e µ(t-w) -e -µ(t-w) (p ji (w) -π i ) dw = K ji .
The first term in the right hand side of ( 3) is easily evaluated, again relying on the fact that M (t) has a Poisson distribution, conditional on J:

E Var(M (t) | J) = d i=1 π i λ i t 0 e -µs ds = 1 -e -µt µ d i=1 π i λ i = (t).
Now we study the effect of the time-scaling: we replace λ i by N λ i (for i = 1, . . . , d) and p ij (w) by

p ij (N α w) (for i, j = 1, . . . , d). Introduce the deviation matrix D, by [D] ij := ∞ 0 (p ij (t) -π j ) dt;
see e.g. [START_REF] Coolen-Schrijner | The deviation matrix of a continuous-time Markov chain[END_REF]. Combining the above results, it is a matter of some elementary algebra to verify that, in obvious notation,

Var M (N ) (t) ∼ N (t) + N 2-α 1 -e -2µt µ d i=1 d j=1 π i λ i λ j [D] ij .
From this relation, the above mentioned dichotomy becomes clear. It is observed that for α > 1 the variance of M (N ) (t) grows linearly in N , and is essentially equal to the corresponding mean, viz. N (t). The intuition here is that in this regime the background process jumps faster than the arrival process, so that the arrival stream is nearly Poisson with parameter d i=1 π i λ i . The resulting system behaves therefore, as N → ∞, essentially as an M/M/∞. If α < 1 the background process is slower than the arrival process. The variance of M (N ) (t) now grows like N 2-α , proportionally to a constant that is a linear combination of the entries of the deviation matrix D.

The above computations were done for the transient number of jobs M (N ) (t), but obviously an analogous reasoning applies to its stationary counterpart M (N ) .

2.3.

Preliminaries on deviation matrices, additional notation. In this subsection, we recall a number of key properties of deviation matrices; for more detailed treatments we refer to e.g. the standard texts [START_REF] Keilson | Markov Chain Models: Rarity and Exponentiality[END_REF][START_REF] Kemeny | Finite Markov chains[END_REF][START_REF] Syski | Ergodic potential[END_REF], as well as the compact survey [START_REF] Coolen-Schrijner | The deviation matrix of a continuous-time Markov chain[END_REF]. We also introduce additional notation, which is intensively used later on.

We define the diagonal matrices Λ and M, where [Λ] ii = λ i and [M] ii = µ i . We denote the invariant distribution corresponding to the transition matrix Q by the vector π; we follow the convention that vectors are column vectors unless stated otherwise, and that they are written in bold fonts.

As π denotes the invariant distribution, we have π T Q = 0 T and π T 1 = 1, where 0 and 1 denote vectors of zeros and ones, respectively. In the sequel we frequently use the 'time-average arrival rate' λ ∞ := d i=1 π i λ i = π T Λ1, and the 'time average departure rate' µ ∞ := d i=1 π i µ i = π T M1. We recall some concepts pertaining to the theory of deviation matrices of Markov processes; see e.g. [START_REF] Coolen-Schrijner | The deviation matrix of a continuous-time Markov chain[END_REF]. In particular, we let Π := 1π T denote the ergodic matrix. We also define the fundamental matrix F := (Π -Q) -1 . It turns out that the deviation matrix D, introduced above, satisfies D = F -Π. We will frequently use the identities QF = F Q = Π -I, as well as the facts that ΠD = DΠ = 0 (here 0 is to be read as an all-zeros d × d matrix) and F 1 = 1.

We use the following three vector-valued generating functions throughout the paper: p denotes the unscaled probability generating function (pgf); p ≡ p(N) denotes the corresponding moment generating function (mgf) under the law-of-large-numbers scaling; and p ≡ p(N) denotes the mgf centered and normalized appropriately for the central limit theorem at hand. For the transient cases, these generating functions involve an extra argument t to incorporate time. Importantly, all three generating functions are vectors of dimension d as we consider distributions jointly with the state of the background process; to make the notation easier, we assume that these vectors are row vectors. Lastly, φ ≡ φ (N ) denotes the scalar mgf under the centering and normalization (obtained by summing the elements of p).

Outline of clt proofs

In this section we point out how we set up our clt proofs. In the next two sections this 'recipe' is then applied to analyze Model i and Model ii, covering both the transient and stationary number of jobs in the system. We use a fairly classical approach to proving the clt s for centered and normalized sequences of random variables of the type (2). More specifically, our objective is to show that under the appropriate normalization (i.e., an appropriate choice of γ), the moment generating function of (2) converges to that of the Normal distribution; the same is done for the stationary counterpart of (2).

Our technique consists of the following steps. (e) Discard asymptotically vanishing terms, so as to obtain a unique limit solution, viz., φ(ϑ) = exp(ϑ 2 σ 2 (t)) in the transient case and φ(ϑ) = exp(ϑ 2 σ 2 ) in the stationary case. We explicitly identify σ 2 (t) and σ 2 . This limit solution resulting from the last step corresponds to a zero-mean Normal distribution.

Due to Lévy's continuity theorem, this pointwise convergence of characteristic functions implies convergence in distribution to the zero-mean Normal random variable, so that we have derived the clt.

Issues related to the uniqueness of the solution of the differential equation are dealt with in Appendix A. Below we demonstrate this proof technique for the special case that the service rates in each of the states are identical, i.e., M = µI for some µ > 0, so that Models i and ii coincide. Importantly, Prop. 1 in Section 3.1 holds for general M.

3.1. Differential equations for the pgf p. First we derive a system of differential equations for the pgf of the number of jobs in the system, jointly with the background state. We consider the bivariate process (M (t), J(t)) t∈R , which is an ergodic Markov process on the state space {1, . . . , d}× N. With the states of this process enumerated in the obvious way, it has the (infinite-dimensional) transition rate matrix

           Q -Λ Λ M Q -M -Λ Λ 2M Q -2M -Λ Λ 3M Q -3M -Λ Λ . . . . . . . . .            .
We set out to find the transient distribution (p k (t)) ∞ k=0 , where p k (t) is a d-dimensional row-vector whose entries are defined by [p k (t)] j := P(M (t) = k, J(t) = j). The (row-vector-)pgf p(t, z) is then defined through

p(t, z) := ∞ k=0 p k (t)z k , such that [p(t, z)] j = E z M (t) 1 {J(t)=j} .
Proposition 1. The pgf p(t, z) satisfies the following differential equation:

∂p(t, z) ∂t = p(t, z) Q + (z -1) p(t, z) Λ - ∂p(t, z) ∂z M .
Proof. The result follows from classical arguments. By virtue of the Chapman-Kolgomorov equation, we have that

(4) dp k (t) dt = p k-1 (t)Λ + p k (t)(Q -Λ -kM) + (k + 1)p k+1 (t)M,
for all k ∈ N, where we put p -1 (t) := 0 for all t ≥ 0.

From the standard relations

∞ k=0 (k + 1)p k+1 (t)z k = ∂p(t, z) ∂z , and
∞ k=0 kp k (t)z k = z ∂p(t, z) ∂z ,
we obtain by multiplying both sides of (4) by z k and summing over k ∈ N,

∂p(t, z) ∂t = zp(t, z)Λ + p(t, z)(Q -Λ) -z ∂p(t, z) ∂z M + ∂p(t, z) ∂z M.
The claim follows directly.

We assume that at time 0 the system starts off empty. Under the scaling Λ → N Λ and Q → N α Q, Prop. 1 implies that we have the following system of partial differential equations governing (M (N ) (t), J (N ) (t)):

(5)

∂p (N ) (t, z) ∂t = N α p (N ) (t, z) Q + (z -1) N p (N ) (t, z)Λ - ∂p (N ) (t, z) ∂z M
describing the pgf p (N ) of the number of jobs in the scaled system.

Mean behavior.

In the remainder of this section, we assume that all service rates are identical: M = µI. To obtain the limiting behavior of N -1 M (N ) (t) when N grows large, it turns out to be convenient to take the following steps.

(i) Rewrite the differential equation ( 5) as a recurrence relation for p (N ) involving the fundamental matrix F ; recall from Section 2.3 the relation QF = Π -I.

(ii) Translate this into a recurrence relation in terms of the mgf p(N) of N -1 M (N ) (t), using a

Taylor expansion for z = exp(ϑ/N ).

(iii) Sum over the possible background states by postmultiplying with 1, so as to obtain a scalar mgf; in this step we make use of the identity F 1 = 1.

(iv) Obtain the limiting differential equation by taking the limit for N → ∞. This equation has a closed solution. This is the mgf of the limiting constant (t).

In this way we have proven the convergence in distribution of N -1 M (N ) (t) to (t); as this limit is a constant, convergence in probability follows immediately.

Let us go through the procedure in full detail now. Postmultiplication of Eqn. ( 5) with F and N -α , using QF = Π -I, results in the recurrence relation

p (N ) (t, z) = p (N ) (t, z) Π + N -α (z -1) N p (N ) (t, z)Λ - ∂p (N ) (t, z) ∂z M F (6) -N -α ∂p (N ) (t, z) ∂t F.
We are now set to state and prove the mean behavior of M (N ) (t). Define (t) := (1 -e -µt ), with

= λ ∞ /µ. Lemma 1. N -1 M (N ) (t) converges in probability to (t), as N → ∞.
Proof. We introduce the transient scaled moment generating function p(N) (t, ϑ):

p(N) (t, ϑ) := p (N ) (t, z), with z ≡ z (N ) (ϑ) = exp(ϑ/N ). Evidently, ∂ p(N) (t, ϑ) ∂t = ∂p (N ) (t, z) ∂t , ∂ p(N) (t, ϑ) ∂ϑ = ∂p (N ) (t, z) ∂z dz dϑ = z N ∂p (N ) (t, z) ∂z .
Substituting these expressions in Eqn. [START_REF] Auria | M/M/∞ queues in semi-Markovian random environment[END_REF] and noting that z

±1 = 1 ± ϑN -1 + O(N -2 ), we obtain p(N) (t, ϑ) = p(N) (t, ϑ)Π + N -α ϑ p(N) (t, ϑ) Λ -ϑ ∂ p(N) (t, ϑ) ∂ϑ µI - ∂ p(N) (t, ϑ) ∂t F + o(N -α ).
The above implies that p(N) (t, ϑ) = p(N) (t, ϑ)Π + O(N -α ), and the same holds for the partial derivatives of p(N) (t, ϑ), so all p(N) (t, ϑ) between the brackets can be replaced by p(N) (t, ϑ) Π.

Postmultiplying by 1N α and using the identities Π1 = 1 and It is straightforward to check that p(t, ϑ)1 = exp(ϑ (t)) satisfies the equation as well as the boundary conditions p(t, 0)1 = 1 and p(0, ϑ)1 = 1. Now the stated follows directly.

F 1 = 1, yields 0 = ϑλ ∞ p(N) (t, ϑ) 1 -µϑ ∂ p(N) (t, ϑ) ∂ϑ 1 - ∂ p(N) (t, ϑ) ∂t 1 + o( 1

3.3.

Recurrence relations for the centered and normalized mgf p(N) . Now that we have derived the weak law of large numbers, we introduce in the next step the centered and normalized mgf p(N) (t, ϑ), that is, centered around N (t) and normalized by N -γ , with the scalar γ yet to be determined. We perform a change of variables in the recurrence relation for p (N ) , Eqn. [START_REF] Auria | M/M/∞ queues in semi-Markovian random environment[END_REF], so as to obtain the recurrence relation for the centered and normalized mgf p(N) .

The pgf p (N ) can be expressed in the normalized and centered mgf p(N) using

p(N) (t, ϑ) = exp(-N (t)ϑ/N γ ) p (N ) (t, exp(ϑ/N γ )) ,
which can be written as

p (N ) (t, z) = exp( (t)ϑN 1-γ ) p(N) (t, ϑ), with z ≡ z (N ) (ϑ) = exp(ϑN -γ ). It is readily verified that ∂p (N ) (t, z) ∂z dz dϑ = exp( (t)ϑN 1-γ ) (t)N 1-γ p(N) (t, ϑ) + ∂ p(N) (ϑ) ∂ϑ ; dz dϑ = N -γ exp(ϑN γ ) = N -γ z,
so the derivatives of p (N ) can be expressed in terms of the corresponding derivatives of p(N) :

∂p (N ) (t, z) ∂t = exp( (t)ϑN 1-γ ) (t)ϑN 1-γ p(N) (t, ϑ) + ∂ p(N) (t, ϑ) ∂t , ∂p (N ) (t, z) ∂z = 1 z exp( (t)ϑN 1-γ ) N (t) p(N) (t, ϑ) + N γ ∂ p(N) (t, ϑ) ∂ϑ .
Now perform the change of variables and substitute the expressions for p (N ) (t, z) and its partial derivatives into Eqn. [START_REF] Auria | M/M/∞ queues in semi-Markovian random environment[END_REF]. Dividing by exp( (t)ϑN 1-γ ) yields the following recurrence relation for p(N) :

p(N) (t, ϑ) = p(N) (t, ϑ)Π + N 1-α z (N ) (ϑ) -1 p(N) (t, ϑ)ΛF -N 1-α 1 - 1 z (N ) (ϑ) (t) p(N) (t, ϑ)MF -N γ-α 1 - 1 z (N ) (ϑ) ∂ p(N) (t, ϑ) ∂ϑ MF -N 1-α-γ (t)ϑp (N ) (t, ϑ)F -N -α ∂ p(N) (t, ϑ) ∂t F. (7)
3.4. Differential equation for the scalar, centered and normalized, mgf φ (N ) . The next step is to expand z in a Taylor series. Assuming certain restrictions on γ (that we later justify) we delete all terms of order smaller than N -α . The resulting recurrence relation is iterated and manipulated until all terms in the right-hand side contain p(N) Π. Next we postmultiply this system of partial differential equations by 1, so as to obtain a scalar partial differential equation in terms of

φ (N ) (t, ϑ) := p(N) (t, ϑ)1.
In this step we make use of the definition of Π := 1π T and the identities Π1 = 1 and F 1 = 1.

The Taylor expansions of z and z -1 are

z ±1 = 1 ± ϑN -γ + 1 2 ϑ 2 N -2γ + O(N -3γ ),
Applying these to Eqn. [START_REF] Coolen-Schrijner | The deviation matrix of a continuous-time Markov chain[END_REF] results in

p(N) (t, ϑ) = p(N) (t, ϑ)Π + ϑN 1-α-γ p(N) (t, ϑ)(Λ -(t)M -(t)I)F + ϑ 2 2 N 1-α-2γ p(N) (t, ϑ)(Λ + (t)M)F -ϑN -α ∂ p(N) (t, ϑ) ∂ϑ MF -N -α ∂ p(N) (t, ϑ) ∂t F + O(N 1-α-3γ ) + O(N -α-γ ). ( 8
)
Under the assumption that γ > 1/3 (to be justified later) the order terms can be replaced by o(N -α ).

Next we iterate Eqn. ( 8) until all terms in the right-hand side either contain p(N) (ϑ)Π or are of O(N -α ). For the latter we assume a second restriction, viz., γ ≥ 1 -α/2 (also justified later). We

thus obtain p(N) (t, ϑ) = p(N) (t, ϑ)Π +ϑN 1-α-γ p(N) (t, ϑ)Π + ϑN 1-α-γ p(N) (t, ϑ)(Λ -(t)M -(t)I)F + O(N 1-α-2γ ) + O(N -α ) (Λ -(t)M -(t)I)F + ϑ 2 2 N 1-α-2γ p(N) (t, ϑ)Π + O(N 1-α-γ ) + O(N -α ) (Λ + (t)M)F -ϑN -α ∂ p(N) (t, ϑ) ∂ϑ MF -N -α ∂ p(N) (t, ϑ) ∂t F + o(N -α );
here we remark that in the O(N -α )-terms p(N) can be replaced by p(N) Π as an immediate consequence of the fact that Eqn. (8) implies p(N) = p(N) Π + o(1), while the same applies to its derivatives. The above equation can be rewritten as

p(N) (t, ϑ) = p(N) (t, ϑ)Π + ϑN 1-α-γ p(N) (t, ϑ) Π(Λ -(t)M -(t)I)F +ϑ 2 N 2-2α-2γ p(N) (t, ϑ) Π(Λ -(t)M -(t)I)F (Λ -(t)M -(t)I)F + ϑ 2 2 N 1-α-2γ p(N) (t, ϑ) Π(Λ + (t)M)F -ϑN -α ∂ p(N) (t, ϑ) ∂ϑ ΠMF -N -α ∂ p(N) (t, ϑ) ∂t ΠF + o(N -α ). (9)
Now postmultiply Eqn. ( 9) by 1 N α ; using the identities Π1 = 1 and F 1 = 1, and the definition Π := 1π T . We obtain

0 = ϑN 1-γ φ (N ) (t, ϑ) π T Λ -(t)M -(t)I 1 +ϑ 2 N 2-α-2γ φ (N ) (t, ϑ) π T (Λ -(t)M -(t)I)F (Λ -(t)M -(t)I)1 + ϑ 2 2 N 1-2γ φ (N ) (t, ϑ) π T (Λ + (t)M)1 -ϑµ ∂φ (N ) (t, ϑ) ∂ϑ - ∂φ (N ) (t, ϑ) ∂t + o(1).
Directly from the definition of (t), it is seen that the first term on the right-hand side vanishes.

In addition, it takes some elementary algebra to check that

π T (Λ -( (t)µ + (t))I)F (Λ -( (t)µ + (t)I)1 = π T ΛDΛ1 =: U,
where we used F = D + Π = D + 1π T , and

1 2 π T (Λ + (t)µI)1 = λ ∞ 1 - e -µt 2 .
This results in the partial differential equation

∂φ (N ) (t, ϑ) ∂t + ϑµ ∂φ (N ) (t, ϑ) ∂ϑ = ϑ 2 φ (N ) (t, ϑ) N 2-α-2γ U + 1 2 N 1-2γ λ ∞ (1 - 1 2 e -µt ) + o(1). (10)
3.5. Limit solution. The last step in our proof is to obtain the limiting differential equation for φ(t, ϑ), being the limit of φ (N ) (ϑ, t). Its unique solution corresponds to a normal distribution

N (0, σ 2 (t)).
First, note that if we choose γ larger than both 1 -α/2 and 1/2, we do not obtain a clt, but rather that the random variable under study converges in distribution to the constant 0. Hence, we take γ = max{1 -α/2, 1/2}, in which case the largest term dominates in [START_REF] Hellings | Semi-Markov-modulated infinite-server queues: approximations by time-scaling[END_REF], with both terms contributing if α = 1. Note that this choice is consistent with the restrictions on γ we used during our proof. We obtain by sending N → ∞, [START_REF] Keilson | Markov Chain Models: Rarity and Exponentiality[END_REF] ∂φ(t, ϑ) ∂t + ϑµ ∂φ(t, ϑ) ∂ϑ = ϑ 2 φ(t, ϑ) g(t),

with g(t) := U 1 {α≤1} + (λ ∞ (1 -e -µt /2)) 1 {α≥1} .
We propose the ansatz

φ(t, ϑ) = exp 1 2 ϑ 2 e -2µt f (t) ,
for some unknown function f (t); recognize the mgf associated with the Normal distribution. This leads to the following ordinary differential equation for f (t):

f (t) = 2e 2µt g(t),
which is obviously solved by integrating the right-hand side. From this we immediately find the expression for the variance σ 2 (t) of the Normal distribution.

With this last step we have proven our claim. It is instructive to compare the findings with the expressions obtained in Section 2.2.

Theorem 1. Consider Model i or ii with µ i = µ for all i ∈ {1, . . . , d}. The random variable

M (N ) (t) -N (t) N γ
converges to a Normal distribution with zero mean and variance σ 2 (t) as N → ∞; here the parameter γ equals max{1 -α/2, 1/2}, and σ 2 (t) := σ 2 m (t)1 {α≤1} + (t)1 {α≥1} , with σ 2 m (t) := µ -1 (1 -e -2µt )U . Corollary 1. Consider Model i or ii with µ i = µ for all i ∈ {1, . . . , d}. The random variable

M (N ) -N N γ
converges to a Normal distribution with zero mean and variance σ 2 as N → ∞; here the parameter γ equals max{1 -α/2, 1/2}, and σ 2 := σ 2 m 1 {α≤1} + 1 {α≥1} , with σ 2 m := µ -1 U .

Model i: Stationary and transient distribution

In this section we briefly recall the results of the steps for Model i, both for the stationary and time-dependent behavior. The proofs are analogous to those in the previous section. Comparing the results of Lemma 2 and Thm. 2 with those of Lemma 1 and Thm. 1, respectively, the effect of heterogeneous service rates becomes visible. Proposition 2. Consider Model i. In the stationary case the pgf p(z) satisfies the following differential equation:

p(z)Q = (z -1) dp(z) dz M -p(z)Λ .
In the transient case the pgf p(t, z) satisfies the following differential equation:

∂p(t, z) ∂t = p(t, z) Q + (z -1) p(t, z) Λ - ∂p(t, z) ∂z M . Define (i) := λ ∞ /µ ∞ , and (i) (t) = (i) (1 -e -µ∞t ). Lemma 2. Consider Model i. As N → ∞, (1) N -1 M (N ) (t) converges in probability to (i) (t).
(2) N -1 M (N ) converges in probability to (i) .

Theorem 2. Consider Model i. The random variable

M (N ) (t) -N (i) (t) N γ
converges to a Normal distribution with zero mean and variance σ 2 (t) as N → ∞; here

σ 2 (t) := σ 2 m (t)1 {α≤1} + (i) (t)1 {α≥1} , with σ 2 m (t) := 2e -2µ∞t t 0 e 2µ∞s π T (Λ -(i) (s)M)D(Λ -(i) (s)M)1 ds.
The random variable

M (N ) -N (i) N γ
converges to a Normal distribution with zero mean and variance σ 2 as N → ∞; here

σ 2 := σ 2 m 1 {α≤1} + (i) 1 {α≥1} , with σ 2 m := µ -1 ∞ π T (Λ -(i) M)D(Λ -(i) M)1.
In both cases the parameter γ equals max{1 -α/2, 1/2}.

The formula for σ 2 m (t) can be evaluated more explicitly. Define G m,n (t) := e -mµ∞t -e -nµ∞t , for m, n ∈ N. Direct computations yield that σ 2 m (t) equals

U 1 µ ∞ G 0,2 (t) + Û (i) µ ∞ (2G 1,2 (t) -G 0,2 (t)) + Ǔ ( (i) ) 2 µ ∞ G 0,2 (t) -4G 1,2 (t) + 2µ ∞ te -2µ∞t ,
with Û := π T MDΛ1 + π T ΛDM1 and Ǔ := π T MDM1. It is readily verified that σ 2 m (t) → σ 2 m as t → ∞, as expected.

Results for Model ii

In this section we study Model ii: the service times are now determined by the background state as seen by the job upon arrival. The approach is as before: we first derive a system of differential equations (Section 5.1), then establish the mean behavior by means of laws of large numbers (Section 5.2), and finally derive the clt s (Section 5.3). 5.1. Differential equations for the pgf p. For the transient distribution, a system of differential equations was previously derived in [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF]. It is based on the observation that M (t) has a Poisson distribution with (random) parameter ϕ(J), see [START_REF] Anderson | A functional central limit theorem for a Markov-modulated infinite-server queue[END_REF]. The intuition behind this formula is that a job arriving at time s survives in the system until time t with probability e -µ i (t-s) (assuming that the background process is in state i), which is distributionally equivalent with 'thinning' the Poisson parameter with exactly this fraction. This description yields, after some manipulations, the following differential equation for the pgf, the row vector p(t, z):

(12) ∂p(t, z) ∂t = p(t, z) Q + (z -1)p(t, z)∆(t),
where Q = (q ij ) d i,j=1 is the transition rate matrix of the time-reversed version of J(•) (i.e., qij := q ji π j /π i ), and ∆(t) denotes a diagonal matrix with entries [∆(t)] ii := λ i exp(-µ i t).

Remark 1. It is noted that the definition of p is slightly different from the one used in [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF]. In the present paper we consider the generating function of the number of jobs present at time t jointly with the state of the background process at time t, whereas [2, Prop. 2] considers the generating function of the number of jobs present at time t conditioned on the background state at time 0. As a consequence, we obtain a slightly different equation, but it is easy to translate them into each other. ♦ Our objective is to set up our proof such that it facilitates proving both the transient and stationary

clt. Naïvely, one could try to obtain a differential equation for the stationary behavior by sending t → ∞ in ( 12), but it is readily checked that this yields a trivial relation only: 0 = 0. A second naïve approach would be to establish the clt for M (N ) (t), and to send then t to ∞; it is clear, however, that this procedure relies on interchanging two limits (N → ∞ and t → ∞), of which a formal justification is lacking.

We therefore resort to an alternative approach. It relies on a description based on a more general state space: we do not only keep track of the number of jobs present, but we rather record the numbers of jobs present of each type, where 'type' refers to the state of the background process upon arrival. To this end, we introduce the d-dimensional stochastic process

M (t) = (M 1 (t), . . . , M d (t)) t∈R ,
where the k-th entry denotes the number of particles of type k in the system at time t. The transient and stationary total numbers of jobs present are denoted by

M (t) := d k=1 M k (t), M := d k=1 M k ,
respectively. As usual, we add a superscript (N ) when working with the model in which imposed our scaling on the arrival rates and the transition rates of the background process.

As before, we first derive a differential equation for the unscaled model. The generating function p(t, z) is defined as follows:

[p(t, z)] j = E d k=1 z M k (t) k 1 {J(t)=j} .
In addition, E k is a matrix for which [E k ] kk = 1, and whose other entries are zero. For a row vector q, the multiplication q E k thus results in a (row) vector which leaves the k-th entry of q unchanged while the other entries become zero. The following result covers the transient case.

Proposition 3. Consider Model ii. The pgf p(t, z) satisfies the following differential equation:

∂p(t, z) ∂t = p(t, z)Q + d k=1 (z k -1) λ k p(t, z) E k -µ k ∂p(t, z) ∂z k .
With the pgf p(z 1 , . . . , z d ) defined in the obvious way, the differential equation for the stationary case is the following.

Proposition 4. Consider Model ii. The pgf p(z) satisfies the following differential equation:

0 = p(z)Q + d k=1 (z k -1) λ k p(z) E k -µ k ∂p(z) ∂z k .
The proofs of these propositions are straightforward, and follow the same lines as before: we consider the generator of the Markov process, and transform the Kolmogorov equation (for the transient case) and the invariance equation (for the stationary case).

The partial differential equation for the transient scaled model follows directly from Prop. 3, by replacing λ k by N λ k , and

Q by N α Q. It results in (13) ∂p (N ) (t, z) ∂t = N α p (N ) (t, z) Q + d k=1 (z k -1) N λ k p (N ) (t, z) E k -µ k ∂p (N ) (t, z) ∂z k .
The stationary case can be dealt with analogously, relying on Prop. 4.

Our objective is to derive the clt for both the transient and stationary case. We do so by presenting the full analysis for the transient case; in the stationary case we can leave out one term. Importantly, this approach does not have the problem of illegitimately interchanging two limits.

Mean behavior.

As before, we first derive the law of large numbers. Again we rewrite the differential equations (13) as a recurrence relation for p (N ) that involves the fundamental matrix

F : p (N ) (t, z) = p (N ) (t, z)Π + N -α d k=1 (z k -1) N λ k p (N ) (t, z) E k -µ k ∂p (N ) (t, z) ∂z k F -N -α ∂p (N ) (t, z) ∂t F (14)
for the transient case, and likewise for the stationary case.

The following lemma establishes weak laws of large numbers for M (N ) (t) and M (N ) (t), as well as their steady-state counterparts M (N ) and M (N ) . We first define

(ii) k (t) := π k λ k µ k (1 -e -µ k t ), (ii)
k := π k λ k µ k . Also, (ii) (t) := k (ii) k (t) and (ii) := k (ii) k . Lemma 3. Consider Model ii. As N → ∞, (1)
N -1 M (N ) (t) converges in probability to (ii) (t).

(2) N -1 M (N ) converges in probability to (ii) .

(3) N -1 M (N ) (t) converges in probability to (ii) (t), and N -1 M (N ) to (ii) .

Proof. Similarly to the proof of Lemma 1, we first introduce the scaled moment generating function

p(N) (t, ϑ) := p (N ) (t, z), with z k ≡ z (N ) k (ϑ k ) = exp(ϑ k /N ), for k = 1, . . . , d. We see immediately that ∂ p(N) (t, ϑ) ∂t = ∂p (N ) (t, z) ∂t , ∂ p(N) (t, ϑ) ∂ϑ k = ∂p (N ) (t, z) ∂z k dz k dϑ k = z k N ∂p(t, z) ∂z k .
Now we substitute these expressions in Eqn. [START_REF] O'cinneide | The M/M/∞ queue in a random environment[END_REF], and note that

z ±1 k = 1 ± ϑ k N -1 + O(N -2 ). As a consequence, p(N) (t, ϑ) = p(N) (t, ϑ)Π + N -α d k=1 ϑ k λ k p(N) (t, ϑ) E k -µ k ∂ p(N) (t, ϑ) ∂ϑ k F -N -α ∂ p(N) (t, ϑ) ∂t F + o(N -α ).
It directly follows that p(N) (t, ϑ) = p(N) (t, ϑ)Π + O(N -α ), and hence also

∂ p(N) (t, ϑ) ∂t = ∂ p(N) (t, ϑ) ∂t Π + O(N -α ), ∂ p(N) (t, ϑ) ∂ϑ k = ∂ p(N) (t, ϑ) ∂ϑ k Π + O(N -α ).
The next step is to postmultiply the previous display by 1 N α , and after some elementary steps we obtain the following scalar partial differential equation in p(N) (t, ϑ)1:

∂(p (N ) (t, ϑ)1) ∂t = d k=1 ϑ k π k λ k (p (N ) (t, ϑ)1) -µ k ∂(p (N ) (t, ϑ)1) ∂ϑ k + o(1). Now let N → ∞; define p(t, ϑ)1 := lim N →∞ p(N) (t, ϑ)1.
We propose the following form for the limiting function p(t, ϑ)1:

p(t, ϑ)1 = exp d k=1 ϑ k ¯ k (t) ,
for specific functions ¯ k (•) (to be determined later). Plugging this form into the differential equation, it means that the following equation must be fulfilled by the ¯ k (•):

d k=1 ϑ k ¯ k (t) -π k λ k + µ k ¯ k (t) = 0.
As this must hold for any ϑ k , this equation leads to a separate differential equation for every ¯ k (t), which moreover agrees with the one in the first part of the claim (¯ k (t) =

(ii) k (t), that is). We conclude that we have established the claim for the transient case:

N -1 M (N ) (t) converges in probability to (ii) (t) as N → ∞.
For the stationary case, we can follow precisely the same procedure, but without the partial derivative with respect to time, so that we now end up with a differential equation in p(ϑ)1 as follows:

0 = d k=1 ϑ k π k λ k (p(ϑ)1) -µ k ∂(p(ϑ)1) ∂ϑ k , for which p(ϑ)1 = exp( d k=1 ϑ k (ii)
k ) forms a solution. This completes the proof of the second claim. The third claim follows trivially. 5.3. Central limit theorems. Next, we state and prove the clt result for Model ii. To this end, we first define the (symmetric) matrices V (t) and V := lim t→∞ V (t) with entries

[V (t)] jk := λ j λ k [ D] jk µ j + µ k (1 -e -(µ j +µ k )t ), [V ] jk = λ j λ k [ D] jk µ j + µ k ;
here D denotes the (symmetric) matrix defined by

[ D] jk = (π j [D] jk + π k [D] kj ). Also, C := lim t→∞ C(t), where [C(t)] jk := [V (t)] jk 1 {α≤1} + (ii) j (t)1 {α≥1} 1 {j=k} .
The following theorem is the main result of this section.

Theorem 3. Consider Model ii. The random vector

M (N ) (t) -N (ii) (t) N γ
p(N) (t, ϑ) = exp -N 1-γ d k=1 ϑ k (ii) k (t) p (N ) (t, z) .
We wish to perform a change of variables in Eqn. [START_REF] O'cinneide | The M/M/∞ queue in a random environment[END_REF] to obtain a recurrence relation in p(N) (t, ϑ).

To this end, note that

∂p (N ) (t, z) ∂z k dz k dϑ k = exp N 1-γ d k=1 ϑ k (ii) k (t) (ii) k (t)N 1-γ p(N) (t, ϑ) + ∂ p(N) (t, ϑ) ∂ϑ k , where dz k dϑ k = N -γ exp(ϑ k N -γ ) = N -γ z k .
Also,

∂p (N ) (t, z) ∂t = exp N 1-γ d k=1 ϑ k (ii) k (t) k ϑ k d (ii) k (t) dt N 1-γ p(N) (t, ϑ) + ∂ p(N) (t, ϑ) ∂t .
Now perform the change of variables, and substitute the expressions for the partial derivatives of p (N ) (t, z) into Eqn. [START_REF] O'cinneide | The M/M/∞ queue in a random environment[END_REF]. Dividing the equation by exp(N 1-γ d k=1 ϑ k (ii) k (t)) gives the following recurrence relation for p(N) (t, z):

p(N) (t, ϑ) = p(N) (t, ϑ)Π + N 1-α d k=1 (z k -1)λ k p(N) (t, ϑ) E k F -N -α d k=1 1 - 1 z k N γ µ k N 1-γ (ii) k (t)p (N ) (t, ϑ) + ∂ p(N) (t, ϑ) ∂ϑ k F -N 1-α-γ d k=1 ϑ k d (ii) k (t) dt p(N) (t, ϑ)F -N -α ∂ p(N) (t, ϑ) ∂t F.
The next step is to introduce the second order Taylor expansions for z k and z -1 k :

z ±1 k = 1 ± ϑ k N -γ + 1 2 ϑ 2 k N -2γ + O(N -3γ ).
Ignoring all terms that are provably smaller than N -α under the assumption that γ > 1/3 (justified later), and combining terms of the same order, we obtain

p(N) (t, ϑ) = p(N) (t, ϑ) Π + N 1-α-γ d k=1 ϑ k p(N) (t, ϑ) λ k E k -µ k (ii) k (t)I - d (ii) k (t) dt I F + N 1-α-2γ d k=1 ϑ 2 k 2 p(N) (t, ϑ) λ k E k + µ k (ii) k (t)I F -N -α d k=1 ϑ k µ k ∂ p(N) (t, ϑ) ∂ϑ k F -N -α ∂ p(N) (t, ϑ) ∂t F,
up to an error term that is o(N -α ). As we did in the proof of the clt in Section 3 with Eqn. (8), we iterate and manipulate this relation, under the assumption that γ ≥ 1 -α/2 (justified later), until all terms in the right-hand side contain p(N) Π. Then we postmultiply with 1 N α , and develop a differential equation in terms of φ (N ) (t, ϑ) := p(N) (t, ϑ) 1. After some (by now quite familiar) manipulations, we obtain the following partial differential equation in φ (N ) (t, ϑ):

∂φ (N ) (t, ϑ) ∂t + d k=1 ϑ k µ k ∂φ (N ) (t, ϑ) ∂ϑ k = 1 2 φ (N ) (t, ϑ)   N 2-α-2γ d j=1 d k=1 ϑ j ϑ k λ j λ k [ D] jk + N 1-2γ d k=1 ϑ 2 k π k (λ k + µ k (ii) k (t)) + o(1),
where we have used that

π T   d j=1 ϑ j (λ j E j -µ j (ii) j (t)I - d (ii) k (t) dt I)   F d k=1 ϑ k (λ k E k -µ k (ii) k (t)I - d (ii) k (t) dt I) 1 = d j=1 d k=1 ϑ j ϑ k λ j λ k π T E j DE k 1 = 1 2 d j=1 d k=1 ϑ j ϑ k λ j λ k [ D] jk .
The last part of the proof concerns the limiting behavior as N → ∞. Pick, as before, γ = max{1α/2, 1/2}, to obtain the following partial differential equation:

∂φ(t, ϑ) ∂t + d k=1 ϑ k µ k ∂φ(t, ϑ) ∂ϑ k = 1 2 φ(t, ϑ)   d j=1 d k=1 ϑ j ϑ k λ j λ k [ D] jk 1 {α≤1} + d k=1 ϑ 2 k (π k λ k + µ k (ii) k (t))1 {α≥1}   .
It is straightforward to verify that the following expression constitutes a solution for this differential equation:

φ(t, ϑ) = exp   1 2 d j=1 d k=1 ϑ j ϑ k [V (t)] jk 1 {α≤1} + 1 2 d k=1 ϑ 2 k (ii) k (t)1 {α≥1}   .
If we redo the derivation for the stationary case (i.e., we now discard the terms originating from the derivative with respect to t in the original partial differential equation), we end up with

φ(ϑ) = exp   1 2 d j=1 d k=1 ϑ j ϑ k [V ] jk 1 {α≤1} + 1 2 d k=1 ϑ 2 k (ii) k 1 {α≥1}   .
This completes the proof.

Corollary 2. Consider Model ii. An immediate consequence of Thm. 3 is that, with γ as defined before, the random variables

M (N ) -N (ii) N γ and M (N ) (t) -N (ii) (t) N γ
converge to Normal distributions with zero mean and variances

d j=1 d k=1 [V ] jk 1 {α≤1} + (ii) 1 {α≥1} and d j=1 d k=1 [V (t)] jk 1 {α≤1} + (ii) (t)1 {α≥1} ,
respectively, as N → ∞.

Correlation across time

Above we analyzed the joint distribution of the two queues at a given point in time. A related question, to be covered in this section, concerns the joint distribution at distinct time epochs. For ease we assume that the service rates are identical (and equal to µ), so that Model i and Model ii coincide.

6.1. Differential equation. We follow the line of reasoning of [2, Prop. 2]; we consider again the non-scaled model, but, as before, these results can be trivially translated in terms of the N -scaled model. Fix time epochs 0 ≡ s 1 ≤ s 2 ≤ • • • ≤ s K for some K ∈ N. The goal of this subsection is to characterize the joint transform, for j = 1, . . . , d,

Ψ j (t, z) := E K k=1 z M (t+s k ) k J(0) = j .
Assume a job arrives between 0 and ∆t, for an infinitesimally small ∆t. Then it is still in the system at time t+s k , but not anymore at t+s k+1 with probability f k (t)-f k+1 (t), where f k (t) := e -µ(t+s k ) .

As a consequence, we obtain the following relation:

Ψ j (t, z) = λ j ∆t b(t, z) Ψ j (t -∆t, z) + i =j q ji ∆t Ψ i (t -∆t, z) +   1 -λ j ∆t - i =j q ji ∆t   Ψ j (t -∆t, z) + o(∆t), where b(t, z) := (1 -f 1 (t)) + z 1 (f 1 (t) -f 2 (t)) + • • • + (z 1 • • • z K-1 )(f K-1 (t) -f K (t)) + (z 1 • • • z K )f K (t).
With elementary manipulations, we obtain

Ψ j (t, z) -Ψ j (t -∆t, z) ∆t = d i=1 q ji Ψ i (t -∆t, z) + a j (t, z)Ψ j (t -∆t, z) + o(1),
where a j (t, z) := λ j (b(t, z) -1) . Now letting ∆t ↓ 0, and defining A(t, z) := diag{a(t, z)}, we obtain the differential equation, in vector notation,

∂ ∂t Ψ(t, z) = (Q + A(t, z))Ψ(t, z).
6.2. Covariance. We now explicitly compute Cov(M (s), M (t)), assuming, without loss of generality, that s ≤ t; the computations are similar to the ones in Section 2.2 (and therefore some steps are left out). The 'law of total covariance', with J ≡ (J(r)) t r=0 , entails that ( 15)

Cov(M (s), M (t)) = E Cov(M (s), M (t) | J) + Cov(E(M (s) | J), E(M (t) | J)).
Due to the fact that M (s) obeys a Poisson distribution with the random parameter ϕ(J), the second term in the right hand side of ( 15) can be written as I 1 + I 2 , where

I 1 := d i=1 d j=1 λ i λ j K ij , where K ij := s 0 v 0 e -µ(s-u) e -µ(t-v) π i (p ij (v -u) -π j ) dudv, I 2 := d i=1 d j=1 λ i λ j L ij , where L ij := s 0 t v e -µ(s-u) e -µ(t-v) π j (p ji (u -v) -π i ) dudv.
It takes some standard algebra to obtain

K ij = e -µt π i s 0 s w e 2µv dv e -µ(s+w) (p ij (w) -π j ) dw = 1 2µ e -µt π i s 0 e µ(s-w) -e -µ(s-w) (p ij (w) -π j ) dw. Similarly, L ij = L (1) ij + L (2)
ij , where

L (1) ij := 1 2µ e -µt π j e µs -e -µs t-s 0 e µw (p ji (w) -π i ) dw, L (2) ij :
= 1 2µ e -µs π j t t-s e µ(t-w) -e -µ(t-w) (p ji (w) -π i ) dw.
Now concentrate on the first term in the right hand side of [START_REF] Schwabe | Transcription stochasticity of complex gene regulation models[END_REF]. To this end, consider the following decomposition:

M (s) := M (1) (s, t) + M (2) (s, t), M (t) := M (2) (s, t) + M (3) (s, t),
where M (1) (s, t) are the jobs that arrived in [0, s) that are still present at time s but have left at time t, M (2) (s, t) the jobs that have arrived in [0, s) that are still present at time t, and M (3) (s, t)

the jobs that have arrived in [s, t) that are still present at time t. Observe that, conditional on J, these three random quantities are independent. As a result,

E Cov(M (s), M (t) | J) = E Var(M (2) (s, t) | J).
Mimicking the arguments used in [START_REF] Auria | M/M/∞ queues in semi-Markovian random environment[END_REF], it is immediate that M (2) (s, t) has a Poisson distribution with random parameter ξ(J), where

ξ(f ) := s 0 λ f (r) e -µ f (r) (t-r) dr.
We conclude that

E Cov(M (s), M (t) | J) = Eξ(J) = d i=1 π i λ i s 0
e -µ(t-r) dr = (s) e -µ(t-s) .

On the other hand, Theorem 2 suggests the following asymptotic expression for the variance of

M (N ) : V 1 (N ) := N (i) 1 {α≥1} + N 2-α σ 2 m 1 {α≤1} .
This expression discards one of the two contributions to the variance, and may therefore be less accurate when both terms are of comparable size. To remedy this effect, we propose the following simple alternative

V 2 (N ) := N (i) + N 2-α σ 2 m ,
which is asymptotically equivalent with V 1 (N ) as N grows large.

In Fig. 1, we illustrate these approximations in the three different regimes σ 2 m , ≈ σ 2 m , and σ 2 m , for a two-state Markov process with generator Q and varying values of α. The parameter values for the three cases are

Q =   -1 1 3 -3   ,   -2 2 1 -1   ,   -1 1 3 -3   ,
and λ = [START_REF] Anderson | A functional central limit theorem for a Markov-modulated infinite-server queue[END_REF][START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF], [START_REF] Anderson | A functional central limit theorem for a Markov-modulated infinite-server queue[END_REF][START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF], [START_REF] Anderson | A functional central limit theorem for a Markov-modulated infinite-server queue[END_REF]50], µ = [2, 1], [100, 1], [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF][START_REF] Anderson | A functional central limit theorem for a Markov-modulated infinite-server queue[END_REF]. We observe that in all cases both approximations tend to the exact values as N gets larger, but the errors are dependent on the specific choices of the parameters of the Markov process. As to be expected, V 2 (N ) is the more accurate one. The contourplots in the middle row give the relative error in the approximation V 1 (N ). They nicely show the effect of the absence of one of the terms in the approximation: for σ 2 m the relative error is almost one if α = 1 -ε, whereas for σ 2 m this is the case for α = 1 + ε. If the two terms are in balance ( ≈ σ 2 m ), we see an increase of the relative error around α ≈ 1, which is absent in approximation V 2 (N ), plotted in the bottom row.

Discussion and conclusion

In this paper we derived central limit theorems (clt s) for infinite-server queues with Markovmodulated input. In our approach the modulating Markov chain is sped up by a factor N α (for some positive α), while the arrival process is sped up by N . Interestingly, there is a phase transition in the sense that the normalization to be used in the clt depends on the value of α: rather than the standard normalization by √ N , it turned out that the centered process should be divided by N γ , with γ equal to max{1 -α/2, 1/2}. We have proved this by first establishing systems of differential equations for the (transient and stationary) distribution of the number of jobs in the system, and then studying their behavior under the scaling described above.

We have also derived a clt for the multivariate distribution of the number of jobs present at different time instants, complementing the analysis for just α > 1 in [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF]. We anticipate weak with the boundary condition φ(0, ϑ 1 , . . . , ϑ d ) = 1 (i.e., an empty system at t = 0) for which the uniqueness question can be tackled in a similar but slightly different fashion (as t is now an explicit variable of the problem). This form occurs in the proofs of Thms. 1 and 3 (as well as in the proofs Lemmas 1 and 3 with the slight difference that there is a negative sign in the ∂/∂t-term, which hardly changes our argument). Indeed, we consider the parametric curve:

(t, ϑ 1 (t), • • • , ϑ d (t), φ(t)) , with the same ode s imposed on ϑ k (t) (and hence having the same solution as well), while dφ(t) dt = g(t, ϑ 1 (t), . . . , ϑ d (t)) φ(t)

has again a solution with one undetermined constant. In order to find the solution at (t, ϑ 1 , . . . , ϑ d ),

we put ϑ k (t) = ϑ k , from which we find ϑ k (0) = ϑ k exp(-µ k t). These relations together with φ(0) = 1 ensure that each ode has a unique solution, and hence the original pde has a unique solution as well.

  (a) Derive a differential equation for the pgf p of the random quantities M (t) and M . (b) Establish the 'mean behavior' (t) ( , respectively) of M (N ) (t) (M (N ) , respectively). This law of large numbers follows by manipulating the mgf p ≡ p(N) , obtaining a scalar limit solution exp(ϑ (t)) in the transient case, and exp(ϑ ) in the stationary case. (c) Reformulate the differential equation for the uncentered and unnormalized pgf p into a recurrence relation for the centered and normalized mgf p ≡ p(N) .

  (d) Manipulate and iterate this equation, approximate by suitable Taylor expansions, to obtain a differential equation for the scalar mgf φ under the chosen centering and normalization.

  ); recall the definitions Π := 1π T and λ ∞ := π T Λ1. Define p(t, ϑ)1 as the limit of p(N) (t, ϑ)1 as N → ∞. Now multiply the differential equation with N α and let N → ∞. We thus obtain a scalar partial differential equation in p(t, ϑ)1 ∂(p(t, ϑ)1) ∂t = ϑλ ∞ (p(t, ϑ)1) -µϑ ∂(p(t, ϑ)1) ∂ϑ .

  converges to a d-dimensional Normal distribution with zero mean and covariance matrix C(t) as N → ∞. In both cases the parameter γ equals max{1 -α/2, 1/2}. The random vectorM (N ) -N (ii)N γ converges to a d-dimensional Normal distribution with zero mean and covariance matrix C asN → ∞.Proof. Mimicking the proof of the clt in Section 3, we start again with setting up a recurrence relation for the centered and normalized mgf p(N) . Define z by z k ≡ z(N ) k (ϑ k ) := exp(ϑ k N -γ), for k = 1, . . . , d, with the value of γ to be determined later on. We first concentrate on the transient case and introduce the centered and normalized mgf p(t, ϑ):

Figure 1 .

 1 Figure 1. Illustration of the behaviour of the approximation in the three different regimes: (left) σ 2 m , (middle) ≈ σ 2 m , (right) σ 2 m . Top row: Plots of the variance of M (N ) along with two approximations; black: α = 0.7; blue: α = 1.0; red: α = 1.3. Full lines represent the exact values, dashed lines represent the first approximation V 1 (N ) and dash-dotted lines represent V 2 (N ). Middle row: Contourplots of the relative error in the approximation V 1 (N ) for varying α and N . Bottom row: Same for approximation V 2 (N ).

Work done while K. de Turck was visiting Korteweg-de Vries with greatly appreciated financial support from Fonds Wetenschappelijk Onderzoek / Research Foundation -Flanders. He is also a Postdoctoral Fellow of the same foundation.

When scaling λ → N λ and Q → N α Q, for α > 0, it is readily verified that for N large, Cov(M (N ) (s), M (N ) (t)) ∼ N (s) e -µ(t-s) + N 2-α e -µt µ e µs -e -µs U, recalling that U := π T ΛDΛ1. When taking s = t, we obtain formulae for the variance that are in line with our findings of Section 2.2.

6.3. Limit results. We again consider the situation in which the modulating Markov chain J(•)

is sped up by a factor N α (for some positive α), while the arrival rates λ i are sped up by N . In this subsection we consider the (multivariate) distribution of the number of jobs in the system at different points in time. While in [START_REF] Blom | Markov-modulated infinite server queues with general service times[END_REF] we just covered the case of α > 1, we now establish a clt for general α.

As the techniques used are precisely the same as before, we just state the result. We first introduce some notation. Define [

Theorem 4. The random vector

converges to a K-dimensional Normal distribution with zero mean and covariance matrix Č(t) as

We observe that the limiting centered and scaled process, as t → ∞, has the correlation structure of an Ornstein-Uhlenbeck process S(t) (at the level of finite-dimensional distributions), that is, the solution to the stochastic differential equation

with W (•) standard Brownian motion.

Numerical illustration

In this section, we briefly illustrate the accuracy of the approximations that are suggested by the limit theorems of this paper. In particular, we consider the variance of the queue content M (N ) under stationarity for Model i. In this case, there is an exact expression for the variance [START_REF] O'cinneide | The M/M/∞ queue in a random environment[END_REF]:

convergence to an Ornstein-Uhlenbeck process with appropriate parameters, but establishing such a claim will require different techniques.

Appendix A. Uniqueness of solutions of the PDEs

In the various proofs of this article, we have 'solved' the differential equations by guessing a solution and establishing that it satisfies both the differential equation itself and the boundary conditions.

We now show that the solutions are indeed unique by relying on the method of characteristics [START_REF] Hilbert | Methoden der mathematischen Physik[END_REF].

The method consists of rewriting the partial differential equation (pde) as a system of ordinary differential equations along so-called characteristic curves, for which the theory of existence and uniqueness is well-developed.

As all occurring pde s are of a similar form and moreover quasi-linear, we can suffice by establishing uniqueness for the two types of pde s, the first of which is as follows:

for some function g(•) with boundary condition φ(0, . . . , 0) = 1. This pertains to differential equations in the proofs of Lemma 3 and Thm. The ode s in ϑ k (t) have the following solution:

while the ode for φ is also quasi-linear with a continuous function g(•), such that a general solution can be found with one undetermined constant. In order to construct the solution at an arbitrary point (ϑ 1 , . . . , ϑ d ), one puts ϑ k (0) = ϑ k and then combines this with the boundary condition 1 = φ(0, • • • , 0), which indeed gives us the condition to make the solution of the ode in φ(t)

unique.

Next, we consider the pde:

µ k ϑ k ∂φ ∂ϑ k = g(t, ϑ) φ(t, ϑ 1 , . . . , ϑ d ),


O Deforges
email: odeforge@insa-rennes.fr

M Babel
email: mbabel@insa-rennes.fr

J Ronsin
email: jronsin@insa-rennes.fr

Codage d'Objets de type VOP Par Représentation en Régions Auto-Extractible

Keywords: Niveau bloc, U/V codés régions : 0.17 bpp Codage interne VOP à 0.45 bpp Y: texture, U/V : niveau bloc

Le nouveau standard multimédia Mpeg-4 traite aussi bien des aspects compression que de fonctionalités de haut niveau telles que la manipulation d'objets. La méthode LAR offre un cadre unique permettant une compression efficace tout en offrant une segmantation à coût « quasi-nul » à travers le principe de représentation en régions auto-extractibles. Aussi montronsnous ici que notre méthode peut constituer une alternative au mode classique de codage intra de Mpeg-4. En particulier, l'approche réussi à unifier le codage de la forme et de la texture d'un VOP.

Abstract -The new multimedia standard Mpeg-4 addresses compression aspects and high level functionalities such as object manipulation.

The LAR method provides an unique framework achieving an efficient compression while allowing a nearly free reliable segmentation map of the image through the principle of self-extracting region representation. Thus, we demonstrate that our method can be a profitable alternative to the traditional intra mode coding of natural video in Mpeg-4. In particular, the approach succeeds to unify the coding of the shape and the texture of a VOP.

Introduction

Mpeg-4 est le nouveau standard multimédia destiné aussi bien au codage de l'audio et de la vidéo, qu'à la synchronisation de ces flux et de leurs combinaisons [START_REF]Signal Processing : Image communication[END_REF]. L'introduction de fonctionnalités au niveau objet telles que la composition de scènes, l'interaction avec l'utilisateur, a nécessité la définition de descripteurs de forme et de codage de contenu associé. Un VOP (Virtual Object Plan) dans Mpeg-4 peut représenter toute une image ou seulement un objet dans celle-ci. La description du VOP concerne le codage de sa forme. Pour compresser son contenu (texture), les techniques classiques de transformée par bloc ont dû être adaptées aux formes quelconques.

La méthode LAR (Locally Adpative Resolution) a été introduite initialement pour le codage d'image en multiniveaux de gris. Des extensions ont ensuite été apportées pour les images couleur et le codage par région d'intérêt. Des travaux récents ont permis d'adapter avantageusement l'approche à du codage sans perte [START_REF] Babel | Décomposition Pyramidale à Redondance Minimale pour Compression d'Image sans Perte[END_REF], ainsi que d'améliorer les performances du codeur basse résolution. Cet article s'attache à montrer plus généralement le schéma de codage complet homogène, depuis le codage très bas débit, jusqu'à la représentation en régions pouvant être utilisée à la fois pour la compression des images de chrominance, et la description de région d'intérêt. En particulier, nous montrons ici que notre approche peut globalement constituer une alternative au mode standard de Mpeg-4, pour le codage d'image intra, mais aussi intégrer de manière plus naturelle la notion de VOP.

Récemment, les résultats obtenus en compression par H.264 ont montré leur supériorité sur un codage basé bloc classique, et même sur des approches ondelettes [START_REF]Performance comparison: H.26L intra coding vs. JPEG2000[END_REF]. L'atout essentiel en codage intra a été apporté par la possibilité de disposer de tailles de blocs différentes en fonction du contenu de l'image. Notre approche suit cette philosophie mais la décline beaucoup plus en avant. Des tests comparatifs avec H.264 pourront être trouvés dans [START_REF] Babel | Adaptive Multiresolution Scheme for Efficient Image Compression[END_REF].

Ce papier est organisé comme suit : Le principe du codeur LAR de base est rappelé section 2. La partie suivante présente l'approche dite de représentation en région autoextractible, et son utilisation pour la compression. L'intégration de VOP est introduite section 4.

Codeur LAR : principes généraux

Le schéma complet du codeur apparaît figure 1. Le codeur LAR pour la compression d'image en multi-niveaux de gris est composé de deux étages : un codeur dit spatial pour les forts de taux de compression, et un codeur spectral pour coder l'image d'erreurs [START_REF] Babel | Adaptive Multiresolution Scheme for Efficient Image Compression[END_REF].

Codeur Spatial

Le principe repose sur la notion de résolution locale (taille de pixel) variable. L'image est découpée en macro-blocs 16x16, eux-mêmes divisés sur une structure de type quad-tree (taille minimale 2x2) selon l'activité locale. En choisissant comme critère une mesure de gradient dans le bloc, une image basse résolution est obtenue en remplissant chaque bloc par sa valeur moyenne. Cette image respecte les contours et la représentation à taille de bloc variable fournit une carte de segmentation grossière de l'image. Une "quantification psychovisuelle" simple peut ainsi être mise en oeuvre en appliquant une quantification forte pour les petits blocs situés sur les contours, et une plus fine pour les grands blocs localisés dans les zones uniformes. Initialement codés par un schéma prédictif de type MICD à simple balayage raster, nous avons développé une nouvelle forme de décomposition pyramidale de type quad-tree à redondance minimale et basé sur des contextes enrichis de prédiction à 360° [START_REF] Babel | Décomposition Pyramidale à Redondance Minimale pour Compression d'Image sans Perte[END_REF]. La méthode offre ainsi un schéma de codage "scalable" basé sur un souséchantillonnage non uniforme spatial, où l'image est transmise progressivement par raffinement local de la résolution. L'image finale des blocs est alors souséchantillonnée d'un facteur 2 (bloc 2x2 pleine résolution -> bloc 1x1). Un post-traitement simple est appliqué pour lisser les zones homogènes (blocs supérieurs à 1), puis un filtre d'interpolation est utilisé pour étendre l'image à la pleine résolution.

Codeur Spectral

Le codage de l'image d'erreur correspond réellement ici à la transmission de la texture locale. Nous appliquons une approche de type DCT à taille de bloc variable, où à la fois la taille et la composante DC sont déjà fournies par le codeur spatial. La nature des blocs permet ici un codage progressif dépendant du contenu : restauration de la texture des zones uniformes à travers les grands blocs et/ou amélioration des contours à travers les blocs 2x2.

Représentation en Régions Autoextractibles

La particularité majeure du codeur spatiale est de ne pas engendrer de fortes distorsions et de respecter les contours. Aussi est-il possible de réaliser une segmentation directement sur cette image. La représentation en régions auto-extractible consiste donc à transmettre l'image basse résolution, puis à effectuer au codeur et au décodeur la phase de segmentation. Si le découpage en blocs peut être vu comme une phase de pré-segmentation (« split »), la phase suivante est donc une phase de fusion effectuée sur les blocs. Notre approche offre ainsi une solution unique de représentation en régions, avec :

un coût nul (pas d'information à transmettre pour la description des contours), -la forme des régions et leur contenu partagent une unique grille de représentation : les blocs. Le but n'est pas ici d'obtenir le meilleur algorithme de segmentation, mais bien de permettre des fonctionalités au niveau région pour les images bas débit.

Notre méthode de segmentation est basée sur les hypergraphes. Elle est présentée brièvement dans le paragraphe suivant. (régions) i, j ∈ V il existe un arc (i,j) si les régions sont adjacentes. Un unique coût est généralement attribué à chaque arc. Le processus de fusion optimal va alors consister à trouver les deux régions adjacentes les plus similaires, puis les fusionner. Si cette algorithme donne de bons résultats, en revanche il est très coûteux en termes de temps de calcul. Afin de faire converger l'algoritme plus rapidemment, la fusion entre deux régions est ici simplement effectuée si leur fonction de coût est inférieure à un seuil. Le processus est réitéré jusqu'à ce qu'aucune nouvelle région ne soit créée. Nous avons introduit une fonction de coût COST (R i , R j ) entre deux régions qui intègre un calcul de distance à la fois sur les luminances moyennes et sur le gradient. Cette fonction est également pondérée par un paramètre lié à la taille de la partition (DIM(R i )), de sorte que les petites régions soient plus rapidement absorbées

Segmentation

Υ K k k R S 1 = = , R i ∩ R j = ∅, ∀i, j ∈ {1, …, K}, pour i ≠ j,
COST" (R i , R j ) = COST (R i , R j ) × log 10 ( DIM(R i ) ) ⇒ COST" (R i , R j ) < COST" (R j , R i ) si DIM(R i ) < DIM(R j )
Le RAG est alors un graphe orienté. Le nombre total de fusion est borné à (K 0 -K F ), si K F est la partition finale. En faisant croître le seuil de fusion afin d'obtenir successivement des partitions finales décroissantes, la représentation en régions devient hiérarchique (voir fig. 3).

Codage de Région d'Intérêt (ROI)

La représentation en régions permet une sélection semiautomatique d'une ROI qui sera déterminée comme un ensemble de régions élémentaires, autrement dit décrite par uniquement les étiquettes des régions concernées. Les régions étant issues des blocs élémentaires, l'amélioration de la qualité de la ROI est immédiate, puisqu'il suffit d'activer le codeur spectral uniquement pour les blocs concernés, c'est à dire appartenant à la ROI.

Codage basé région des composantes U/V

Dans les méthodes classiques, le passage au codage des images couleur consiste essentiellement à dupliquer le schéma général pour les trois composantes. Dans notre approche, nous tirons à la fois profit de l'information couleur pour la phase de segmentation, et utilisons la représentation en régions pour coder à très bas débit l'information de chrominance.

Les images sont traitées dans un espace YUV, et les tailles de bloc sont estimées à partir des trois images. Les résultats ont montré que le codeur spatial seul appliqué aux composantes UV permet d'obtenir une très bonne qualité d'image d'un point de vue couleur. Il est toutefois possible de compresser plus fortement l'information en fournissant uniquement les valeurs de chrominance moyenne de chaque région issue de la segmentation. Afin de garantir des régions cohérentes d'un point de vue chrominance, nous avons introduit le principe de "contrôle de la chrominance" dans le processus de segmentation : la segmentation reste basée sur la seule luminance pour la recherche des plus proches voisins, mais la fusion finale est validée au codeur par le calcul d'une distance chromatique entre les deux candidats. Un bit est donc transmis pour chaque tentative de fusion, mais les composantes YUV étant fortement corrélées, seuls quelques rejets seront effectifs, engendrant une entropie de l'information très faible. Le surcoût apporté par ce contrôle est typiquement compris entre 0.01 et 0.04 bpp, mais permet d'obtenir au décodeur des régions globalement cohérentes d'un point de vue couleur.

Le réhaussement des composantes UV consiste éventuellement à utiliser le codeur spatial pour la compression de l'erreur de chrominance par bloc, de manière globale sur toute l'image ou localement dans une ROI. Le codage basé région des images UV et la sélection d'une ROI peuvent s'appuyer sur des niveaux hiérarchiques différents de la représentation en régions. Des résultats de segmentation sont donnés figure 3. Des tests comparatifs uniquement en termes de compression entre Jpeg2000, Mpeg-4 mode intra (issue de MoMuSys) et le LAR , apparaissent figure 4.

Codage de VOP

Pour un VOP de forme quelconque, il est nécessaire de transmettre en général ses contours, puis de coder sa texture interne. Pour les blocs contours, il faut recourir à des transformées spécifiques de type SA-DCT.

Considérant une image intra (I) codée en LAR, nous proposons plusieurs alternatives pour la description et le codage de la texture du VOP.

• Transmission préalable de l'image basse résolution, VOP résultant de la segmentation Le VOP est décrit uniquement pour des étiquettes de régions contenues dans celui-ci (= ROI). La représentation hiérarchique en régions permet localement un choix fin de résolution. Le choix du VOP peut être effectué aussi bien au codeur qu'au décodeur. Ce mode est le moins coûteux d'un point de vue description du VOP.

Les deux modes suivants considèrent un choix de VOP indépendant de notre segmentation. Afin que la représentation en blocs soit adaptée au VOP, le masque binaire du VOP est également utilisé dans l'estimation des tailles.

• Transmission préalable de la forme du VOP Schéma classique engendrant un coût de description. L'utilisation du codeur LAR permet néanmoins une dissociation spatiale des blocs entièrement extérieurs ou intérieurs au VOP.

• Transmission préalable de l'image basse résolution, VOP issu d'une phase de segmentation quelconque Notre processus de segmentations en région peut alors être effectué en incorporant un mécanisme de contrôle par la forme du VOP identique à celui de la chrominance. Cela garantit ainsi des régions entièrement externes ou internes au VOP. La figure 5 illustre un tel schéma de codage, avec Y compressée par le codeur spatial (0.127 bpp), U/V codés régions (119 régions, contrôle chromatique : 0.034 bpp + valeur par région : 0.0067 bpp). Le coût additionnel du contrôle du VOP est de 302 bits (0.003 bpp), alors que le coup de codage de cette même forme est de 452 bits (source MoMuSys). Le VOP global est ici reconstruit par 12 régions élémentaires, et réhaussé par correction de la luminance et chrominance.

Conclusion et perspectives

Nous montrons ici un schéma de codage général alliant à la fois progressivité, efficacité dans la compression, et fonctionnalités au niveau région disponibles de manière automatique et à très faible coût. L'approche à plusieurs couches de codage conserve la même structure d'échantillonnage, permettant ainsi une totale concordance entre forme des objets et codage de leur texture. Nous montrons également que le principe de VOP introduit dans Mpeg-4 peut s'intégrer de manière plus naturelle dans notre schéma.

FIG. 2 :

 2 FIG. 2: codeur spatial , a) image source, b) grille à taille variable (0.003 bpp), c) Image des blocs (0.13 bpp), d) après postraitement.

Soit

  S={ (x,y) : 1 ≤ x ≤ N c , 1 ≤ y ≤ N r } les coordonnées spatiales d'un pixel dans une image de N r lignes et N c colonnes. La segmentation d'une image en K régions consiste à trouver la partition ∆ K de S telle que :

FIG. 3 :

 3 FIG. 3 : Exemples de représentation en régionobtenues par contrôle de chrominance

CODEUR SPATIAL CODEUR REGIONS UV Erreurs de Blocs CODEUR CHROMINANCE Blocs Y

		CODEUR LUMINANCE
			Image source (Y)		Images U/V
			Sous-échantil.
			Non-uniforme
		CODEUR SPATIAL	Moyennage blocs Quantification blocs Grille	Region Représentation en régions	Moyennage régions Segmentation en régions
				Post-		Quantification	Moyennage
			traitement		régions	blocs
			+ Image basse résolution -	Sélect. ROI	+ UV basées région -
	CODEUR SPECTRAL	FDCT blocs 2x2 Coef AC quant./coding	FDCT blocs 16x16 Coef AC quant./coding Image texture	Blocs validés	Quantification blocs traitement Post-
		Texture blocs 2x2	Texture blocs 16x16		Image basse résolution UV
				FIG. 1 : Schéma de principe général