A Versatile Model for VoD Buzz Workload: Identification, Numerical Validation and Applications in Dynamic Resource Management

Keywords: Cloud Networking, Probabilistic Resource Management, Epidemic Model, Workload Generator, Statistical Estimation, Large Deviation Principle, Service Level Agreements, Video on Demand, Buzz

Dynamic resource management has become an active area of research in the Cloud Computing paradigm. Cost of resources varies significantly depending on configuration for using them. Hence efficient management of resources is of prime interest to both Cloud Providers and Cloud Users. In this report we suggest a probabilistic resource provisioning approach that can be exploited as the input of a dynamic resource management scheme. Using a Video on Demand use case to justify our claims, we propose an analytical model inspired from standard models developed for epidemiology spreading, to represent sudden and intense workload variations. As an essential step we also derive a heuristic identification procedure to calibrate all the model parameters and evaluate the performance of our estimator on synthetic time series. We show how good can our model fit to real workload traces with respect to the stationary case in terms of steady-state probability and autocorrelation structure. We find that the resulting model verifies a Large Deviation Principle that statistically characterizes extreme rare events, such as the ones produced by "buzz effects" that may cause workload overflow in the VoD context. This analysis provides valuable insight on expectable abnormal behaviors of systems. We exploit the information obtained using the Large Deviation Principle for the proposed Video on Demand use-case for defining policies (Service Level Agreements). We believe these policies for elastic resource provisioning and usage may be of some interest to all stakeholders in the emerging context of cloud networking.

Résumé : La gestion dynamique de ressources est un élément clé du paradigme de cloud computing et plus récemment de celui de cloud networking. Dans ce contexte d'infrastructures virtualisées, la réduction des coûts associés à l'utilisation et à la ré-allocation des ressources contraint les opérateurs et les utilisateurs de clouds à une gestion rationnelle de celles-ci. Dans ce travail nous proposons une description probabiliste des besoins liée à la volatilité de la charge d'un service de distribution de vidéos à la demande. Cette description peut alors servir de consigne (input) à la provision et à l'allocation dynamique des ressources nécessaires. Notre approche repose sur la construction d'un modèle stochastique inspiré des modèles de Markov standards de propagation épidémiologique, capable de reproduire des variations soudaines et intenses d'activité (buzz). Nous proposons alors une procédure heuristique d'identification du modèle à partir de séries temporelles du nombre d'utilisateurs connectés au serveur. Les performances d'estimation de chacun des paramètres du modèle sont évaluées numériquement, et nous vérifions l'adéquation du modèle aux données en comparant les distributions des états stationnaires ainsi que les fonctions d'autocorrélation des processus. Les propriétés markoviennes de notre modèle garantissent qu'il vérifie un principe de grandes déviations permettant de caractériser statistiquement l'ampleur et la durée d'évènements extrêmes et rares tels que ceux produits par les buzzs. C'est cette propriété que nous exploitons pour dimensionner le volume de ressources (e.g. bande-passante, nombre de serveurs, taille de buffers) à prévoir pour réaliser un bon compromis entre coût de re-déploiement des infrastructures et qualité de service. Cette approche probabiliste de la gestion des ressources ouvre des perspectives sur les politiques de Service Level Agreement adaptées aux clouds et servant au mieux les intérêts des opérateurs de réseaux, de services et de leurs clients.

Mots-clés : Réseaux, Cloud, Gestion probabiliste des Ressources, Modèles Epidémiques , Générateur de Charge, Estimation Statistique, Principe de Grandes Déviations, Service Level Agreement, Vidéo à la Demande, Buzz

Introduction

In recent trend of data-intensive applications with pay-as-you-go execution in a cloud environment, there are new challenges in system management and design to optimize the resource utilization. Types of the application, deployed in a cloud, can be very diverse. Some applications exhibit highly varying demand in resources. In this paper we consider a Video on Demand (VoD) system as a relevant example of a data-intensive application where bandwidth usage varies rapidly over time. A VoD service delivers video contents to consumers on request. According to Internet usage trends, users are increasingly getting more involved in the VoD and this enthusiasm is likely to grow. According to 2010 statistics a popular VoD provider like Netflix accounts for around 30 percent of the peak downstream traffic in the North America and is the "largest source of Internet traffic overall" [START_REF] Sandvine | Sandvine's spring 2011 global internet phenomena report reveals new internet trends[END_REF]. Since VoD has stringent streaming rate requirements, each VoD provider needs to reserve a sufficient amount of server outgoing bandwidth to sustain continuous media delivery (we are not considering IP multicast here). However, resource reservation is very challenging in a situation, when a video becomes popular very quickly leading to a flood of user requests on the VoD servers. This situation, also known as a "buzz", demands an adaptive resource allocation strategy to cope with the sudden (and significant) variation of workload. Following is one example of "buzz" (see Figure 1) where interest over a video "Star Wars Kid" [START_REF] Andy | Star kids the data dump[END_REF] grew very quickly within a short timespan. According to [START_REF] Bbc | Star wars kid is top viral video[END_REF] it was viewed more than 900 millions times within a short interval of time making it one of the top viral videos. Such bandwidth volatility creates significant challenges to meet, namely, both the desired QoS and efficient resource allocation. A sensible approach to this problem is to help the providers in better understanding and capturing the underlying characteristics of their applications. For example, if the information diffusion process follows a gossip-(or epidemic-) behavior, the rate at which a viewer keeps gossiping about a video and for how long (in average).

Inria

Number of downloads / day Time (days)

Figure 1 -Video server workload : time series displaying a characteristic pattern of flash crowd (buzz effect). Trace obtained from [START_REF] Andy | Star kids the data dump[END_REF].

In this report we follow a constructive approach to propose a stochastic epidemic workload generator for a VoD system based on a Markov model. We show that it succeeds to reproduce the traffic volatility, as exhibited in a real trace, including the buzz occurrence. But the principal interest of our model is that it verifies a Large Deviation Principle (LDP) that gives a probabilistic description of the mean workload of the system over different time scales. It thus adequately allows for statistically characterizing extreme rare events such as the ones produced by buzz transients. Our ultimate objective is to exploit this large deviation information as an input of an probabilistic resource management scheme. However, in order the proposed model to conform with this objective, it needs to be "identifiable" and easily calibrated on real data. The corresponding estimation procedure may not be trivial, since the VoD model is a non-parsimonious model and accounts for complex dynamics. In this report we propose a complete framework for the operators to identify the VoD model parameters based on a server workload trace After parameter estimation we devise two possible and generic ways to exploit the large deviation information in the context of probabilistic resource provisioning. They can serve as the input of resource management functionalities of the Cloud environment. It is evident that we can not define elasticity without the notion of a time scale ; the Large Deviation Principle (LDP) is capable of automatically integrating the time resolution in automatic description of the system. It is to be noted that Markovian processes do satisfy the LDP, but so do some other models as well. Hence, our proposed probabilistic approach is very generic and can adapt to address any provisioning issues, provided the resource volatility can be resiliently represented by a stochastic process for which the LDP holds true. In a nutshell our contributions in this report include :

-A Markov based versatile model to generate VoD workload, -A heuristic identification procedure for the proposed workload model, -A numerical evaluation of the estimator for each parameter of the model, -A real case study to assess the adequacy of our model to fit video workload traces, -An analysis of the Large Deviation property of the proposed Markovian model, -A discussion on the generic ways to exploit the large deviation information in the context of probabilistic resource provisioning.

Moreover, since we followed a constructive approach, each parameter of the model accounts for a specific component of the system, and so, its estimated value also permits to quantify the importance of the corresponding dynamic effect.

Rest of the paper is organized as follows. In Section 2 we discuss the related works. We describe our model and further analyze it in Section 3. Section 4 outlines the parameter estimation procedure and validates the procedure against synthetic workload traces. In Section 5 we validate both our model and the estimation procedure against the real workload traces. Section 6 presents Large Deviation Principle and numerical interpretations of the Large Deviation Spectrum. Section 7 deals with the probabilistic provisioning scheme, derived from the Large Deviation Spectrum for our use case. Finally we conclude and discuss future works in Section 8.

Inria 2 Related Work

Information dissemination in a VoD system has been an active area of research. In [START_REF] Eugster | Epidemic information dissemination in distributed systems[END_REF], it has been already demonstrated that the epidemic algorithms can be used as an effective solution for information dissemination in a VoD like P2P systems. However, in this model an individual process must have a precise idea about the total number of processes in the system. Scalability is also another challenge that the authors addressed in this work. The authors of [START_REF] Bonald | Epidemic live streaming : Optimal performance trade-offs[END_REF] studied random epidemic strategies like the random peer, latest useful chunk algorithm to achieve optimal information dissemination. But main objective of this work is to demonstrate ways to achieve performance trade-offs using unstructured, epidemic live streaming systems. However, it does not bring any information about the underlying dynamics of the streaming system. Authors of [START_REF] Bakhshi | A modeling framework for gossipbased information spread[END_REF] similarly discussed an analytical framework for gossip protocols based on the pairwise information exchange between interacting nodes. However, this model only provides an analysis of networks with lossy channels. Another relevant work to our study is derived in [START_REF] Caron | Pattern matching based forecast of non-periodic repetitive behavior for cloud clients[END_REF] where the authors proposed an approach to predict workload for cloud clients. They considered an autoscaling approach for resource provisioning and validated the result with real-world cloud client application traces. However, this work depends on similar past occurrences of the current shortterm workload history and is not appropriate to deal with sudden and short large variations of workload, as the ones produced by buzz effects. Authors of [START_REF] Garcia | Statistical characterization of a real video on demand service : User behaviour and streaming-media workload analysis[END_REF] show a statistical study of streaming traffic. They analyzed VoD session characteristics, amount and types of media delivered, popularity and daily access profile in order to develop a workload generator. However, the model does not involve the dynamics of the process itself, ergo it is not naturally adapted to infer dynamic resource allocation strategies. Authors of [START_REF] Li | Performance model of interactive videoon-demand systems[END_REF], [START_REF] Melendi | Multivariate distributions for workload generation in video on demand systems[END_REF] and [START_REF] Kanrar | Analysis and implementation of the large scale video-on-demand system[END_REF] also develop user activity models to describe the usage of system resources. Limitation of these models are that they only give average results. However, dealing with mean workloads might not be sufficient to clearly describe applications because of their potential volatility. In [START_REF] Pacheco-Sanchez | Markovian workload characterization for qos prediction in the cloud[END_REF] authors proposed a maximum likelihood method for fitting a Markov arrival process (MAP) to the web traffic measurements, collected in commonly available HTTP web server traces. This method achieves reasonable accuracy in predictive models for web workloads but lacks intuitive nature to describe users behavior like a gossip based method. In [START_REF] Niu | Understanding demand volatility in large vod systems[END_REF] the authors statistically model traffic volatility in large scale VoD systems using GARCH (generalized autoregressive conditional heteroscedasticity) process. Amazon Cloud-Watch follows this approach and provides a free resource monitoring service to Amazon Web Service customers for a given frequency. Based on such estimates of future demand, each VoD provider can individually reserve a sufficient amount of bandwidth to satisfy in average its random future demand within a reasonable confidence. However, according to the authors, this technique only models and forecasts the mean demand, or the expected demand whereas the real demand might vary around this predicted mean. They suggested to provision an additional "risk premium" to the service providers for tolerating the demand fluctuation. In another workload model the authors of [14] [15] proposed a Markov Modulated Poisson Process (MMPP) based approach for buzz modeling and then parameter estimation using the index of dispersion. However, the MMPP model includes only short-term memory in the system and the obtained statistics is not physically interpretable to draw inference about the system dynamics. The model we derive in section 3 of this report has the following advantages :

-It follows a constructive approach, based on a Markov model, -It is identifiable and succeeds to capture workload volatility, -It satisfies the large deviation properties, that can be exploited to frame dynamic resource allocation strategies.

A VoD system and its modeling

A VoD service delivers video contents to consumers on request. According to Internet usage trends, users are increasingly getting more involved in the VoD and this enthusiasm is likely to grow. A popular VoD provider like Netflix accounts for around 30 percent of the peak downstream traffic in the North America and is the "largest source of Internet traffic overall" [START_REF] Sandvine | Sandvine's spring 2011 global internet phenomena report reveals new internet trends[END_REF]. In a VoD system, consumers are video clients who are connected to a Network Provider. The source video content is managed and distributed by a Service Provider from a central data centre. With the evolution of Cloud Computing and Networking, the service in a VoD system can be made more scalable by dynamically distributing the caching/transcoding servers across the network providers. Video service providers interact with the network service providers and describe the virtual infrastructures required to implement the service (like the number of servers required, their placements and clustering of resources). The resource provider reserves resource for certain time period and may change it dynamically depending on resource requirement. Such a dynamic approach brings benefits of cost saving in the system through dynamic resource provisioning which is important for service providers as VoD workload is highly variable by nature. However, since the virtual resources used by Cloud Networking have a set-up time which is not negligible, analysis and provisioning of such a system can be very critical from the operators perspective (capex versus opex trade-off). Figure 2 shows a VoD schematic where the back-end server is connected to the data centre and the transcoding (caching) servers are placed across the network providers. Since VoD has stringent streaming rate requirements, each VoD provider needs to When multiple VoD providers (such as Netflix) are on board to use cloud services from cloud providers, there will be a market between VoD providers and cloud providers, and commodities to be traded in such a market consist of bandwidth reservations, so that VoD streaming performance can be guaranteed.

As a buyer in such a market, each VoD provider can periodically make reservations for bandwidth capacity to satisfy its random future demand. A simple way to achieve this is to estimate expectation and variance of its future demand using historical demand information, which can easily be obtained from cloud monitoring services. As an example, Amazon Cloud-Watch provides a free resource monitoring service to Amazon Web Service customers for a Inria given frequency. Based on such estimates of future demand, each VoD provider can individually reserve a sufficient amount of bandwidth to satisfy in average its random future demand within a reasonable confidence. However, this information is not helpful in case of a "buzz" or a "flash crowd" when a video becomes popular very quickly leading to a flood of user requests on the VoD servers. In situations like the one described in Figure 1, variance estimation or more generally steady state distribution can not explain burstiness of such event as time resolution is excluded from the description. The LDP, by virtue of its multi-resolution extension of the classical steadystate distribution, can describe the dynamics of rare events like this, which we believe can be of some interest for the VoD service providers.

Markov Model to describe the VoD user behavior

Epidemic models commonly subdivide a population into several compartments : susceptible (noted S) to designate the persons who can get infected, and contagious (noted C) for the persons who have contracted the disease. This contagious class can further be categorized into two parts : the infected subclass (I) corresponding to the persons who are currently suffering from the disease and can spread it, and the recovered class (R) for those who got cured and do not spread the disease anymore [START_REF] Barrat | Dynamical Processes on Complex Networks[END_REF]. There can be more categories that fall outside the scope of our current work. In these models (N S (t)) t≥0 , (N I (t)) t≥0 and (N R (t)) t≥0 are stochastic processes representing the time evolution of susceptible, infected and recovered populations respectively. Similarly, information dissemination in a social network can be viewed as an epidemic spreading (through gossip), where the "buzz" is a special event where interest for some particular information increases drastically within a very short period of time. Following the lines of related works, we claim that the above mentioned epidemic models can appropriately be adapted to represent the way information spreads among the users in a VoD system. In the case of a VoD system, infected I refers to the people who are currently watching the video and can spread the information about it. In our setting, I directly represents the current workload which is the current aggregated video requests from the users. Here, we consider the workload as the total number of current viewers, but it can also refer to total bandwidth requested at the moment. The class R refers to the past viewers. In contrast to the classical epidemic case, we introduce a memory effect in our model, assuming that the R compartment can still propagate the gossip during a certain random latency period. Then, we define the probability within a small time interval dt, for a susceptible individual to turn into an active viewer, as follows :

P S→C = (l + (N I (t) + N R (t)) β)dt + o(dt) (1)
where β > 0 is the rate of information dissemination per unit time and l > 0 fixes the rate of spontaneous viewers. The instantaneous rate of newly active viewers in the system at time t is thus :

λ(t) = l + (N I (t) + N R (t))β. (2)
Equation ( 2) corresponds to the arrival rate λ(t) of a non-homogeneous (state dependant) Poisson process. This rate varies linearly with N I (t) and N R (t).

To complete our model we assume that the watch time of a video is exponentially distributed with rate γ. As already mentioned, it also deems reasonable to consider that a past viewer will not keep propagating the gossip about a video indefinitely, but remains active only for a latency random period that we also assume exponentially distributed with rate µ (in general µ ≪ γ). Another important consideration of the model is the maximum allowable viewers (I max ) at any instant of time. This assumption conforms to the fact that the resources in the system are physically limited. For the sake of numerical tractability and without loss of generality, we also assume the number of past (but spreading rumour) viewers at a given instant to be bounded by a maximum value (R max ). With these assumptions, and posing (N I (t) = i, N R (t) = r) the current state of the Markov processes, the probability that the process reaches a different state (i ′ < I max , r ′ < R max ) at time t + dt (dt being small) reads :

P(i ′ , r ′ |i, r) (3)
= (l + (i + r)β)dt + o(dt) for (i ′ = i + 1, r ′ = r), 1 = (γi)dt + o(dt) for (r ′ = r + 1, i ′ = i -1), = (µr)dt + o(dt) for (r ′ = r -1, i ′ = i), = o(dt) otherwise.
This process defining the evolution of the current viewer and past viewer populations is a finite and irreducible Markov chain. It is to be noted that l > 0 precludes the process to reach an absorbing state. This chain is ergodic and admits a stationary regime. Above mentioned descriptions define the mechanism of information dissemination in the community in normal situations. A buzz event differs from this situation by a sudden increase of the dissemination rate β. In order to adapt the model to buzz we resort to Hidden Markov Model (HMM) to be able to reproduce the change in β. Without loss of generality we consider only two states. One with dissemination rate β = β 1 corresponds to the buzz-free case described above, and another hidden state corresponding to the buzz situation, where the value of β increases significantly and takes on a value β 2 ≫ β 1 . Transitions between these two hidden and memoryless Markov states occur with rates a 1 and a 2 respectively (see Figure 3). These rates characterize the buzz in terms of frequency, magnitude and duration. A closed-form expression for the steady state distribution of the workload (i) of this model is not trivial to derive. However, we could easily express the analytic mean workload of the system solving the flow balance equation, i.e. equaling the incoming and outgoing flow rates in steady regime. For ease, we start with β = β 1 = β 2 and generalize the result to β 1 = β 2 thereafter. We get :

E(i) = µl µγ -µβ -γβ , (4)
which, to be a positive and finite quantity, yields the stability criterion in buzz-free regime :

β -1 > µ -1 + γ -1 . ( 5
)
1. In a closed system, where the total number of viewers (susceptible, current and past) is constant, say N , the transition probability for (i ′ = i+1,r ′ = r) needs to be modified, since it would then depend on the number of susceptible viewers as well" i.e (Nir). The transition probability in this case would be (l + (i + r) β N )(Nir)dt + o(dt). Therefore, Eq. 4 and 5 need to be modified accordingly.

Inria

We now extend these results to the case where the model may exhibit a buzz activity. As β alternates between the hidden states β = β 1 and β = β 2 , with respective state probabilities a 2 /(a 1 + a 2 ) and a 1 /(a 1 + a 2 ), one can simply replace β in Eq. ( 4) and ( 5) with the equivalent average value :

β = β 1 a 2 a 1 + a 2 + β 2 a 1 a 1 + a 2 . (6)
In order to illustrate the flexibility of our workload model and to validate Eq. ( 4), we generate three synthetic traces corresponding to the different sets of parameters verifying the stability condition of relation ( 5) and reported in Table 1. Particular realizations of these processes generated over 2 21 points are displayed in Figure 4. While the synthetic traces corresponding to cases (b) and (c) reproduce distinct and easily identifiable buzz regimes, the parameter set of case (a) leads to a workload variation distinct from the typical shape of Figure 1. Nonetheless, for all 3 configurations, the empirical means estimated from the 2 21 samples of the traces are in good agreement with the expected values of Eq. ( 4). Finally, let us notice that even though we consider exponentially distributed random variables in our model, any other distributions could be used, which, according to the same balance principle, would lead to a mean workload and to a stability condition of the same kind as (4) and [START_REF] Bonald | Epidemic live streaming : Optimal performance trade-offs[END_REF]. However, the estimation procedure we derive in the next section strongly relies on the exponential assumption and it would need to be thoroughly reworked to adapt to different hypotheses.

Estimation procedure

In this section, we address the identifiability of our model and design a calibration algorithm to fit workload data. We start constructing empirical estimators for each parameter of the model and we numerically evaluate their performance on synthetic traces.

Parameters estimation

Considering a standard epidemic process X with propagation rate θ, the maximum likelihood estimate θ MLE is derived in [START_REF] Barrat | Dynamical Processes on Complex Networks[END_REF], [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF] and reads : Table 1 for the parameter values corresponding to each of these three cases. The X-axis corresponds to time (in hours unit) while the Y -axis indicates the number of active viewers.

θ MLE = n • T 0 X(t) dt -1 , (7) case (a
where n is the number of contaminations (i.e. number of increments of X) occurring within the time interval T .

Very often, maximum likelihood approach yields optimal results (in terms of estimate variance and or bias) but it is not always possible to get a closed-form expression for the estimated parameters. This can either be due to the likelihood function that is impossible to derive analytically, or to missing data that preclude straightforward application of the maximum likelihood principle. Nonetheless, solutions, such as the Expectation-Maximization (EM) or the Monte Carlo Markov Chain (MCMC) algorithms exist, which in some cases can approximate maximum likelihood estimators.

Returning to our model depicted in Figure 3, each parameter needs to be empirically estimated, assuming that the instantaneous workload time series is the only available observation.

Watching parameter γ. As γ is the departure rate of users that leave the infected state after they finished watching a video, it can directly be inferred from the number n of decrements of the observable process I(t). Therefore, the MLE of Eq. ( 7) straightforwardly applies and leads to :

γ MLE = n • T 0 I(t) dt -1 . ( 8
)
Memory parameter µ. This rate at which past viewers leave the recovery compartment and stop propagating the virus (gossip), relates to the decrement density of the non-observed process R(t).

It is thus impossible to simply apply the MLE of Eq. ( 7) unless we first construct a substitute R(t) to the missing data from the observable data set I(t). Let us recall that in our model, all current viewers turn and remain contagious for a mean period of time γ -1 + µ -1 . Then, in first approximation, we can consider that R(t) derives from the finite memory cumulative process :

R(t) = t t-(γ -1 +µ -1 ) I(u) du, (9)
which itself, depends on the parameter to be estimated µ. We propose an estimation procedure based on the inherent exponential property of the model. From the Poisson assumption, the inter-arrival time w between the consecutive arrivals of two new viewers is an exponentially distributed random variable such that E (w|

I(t) + R(t) = x) = (β x + l) -1 .
Its means that, for x fixed, the normalized random variable w = w/E(w|x) is exponentially distributed with unitary parameter and becomes independent of x. Ideally then, for each value of R(t) + I(t) = x, all the Inria sub-series w x = {w n : R(t n )+I(t n ) = x}, after normalization by their own empirical mean, yield independent and identically distributed realizations of a unitary exponential random variable.

In practice though, as R(t) is not observable, only if R(t) is accurately estimated, should this unitary exponential i.i.d. assumption hold true. From there, we propose the following algorithm : for different values of µ spanning a reasonable interval, we use R µ (t) estimated from Eq. ( 9) to build the normalized series w µ . A statistical test applied to each w µ allows for assessing the exponential i.i.d. hypothesis and then to select the value of µ that yield the best score.

More concretely, we apply to w µ = ( w n ) n=1,...,N the statistical exponentially test derived in [START_REF] Jammalamadaka | Testing exponentiality by comparing the empirical distribution function of the normalized spacings with that of the original data[END_REF] : Form the normalized spacings

v µ = v (n) = (N -n + 1)( w (n) -w (n- 1
) ) n=1,...,N where w (n) n=1,...,N stands for w µ rearranged in ascending order. Let F and G denote the cumulative distribution functions of w µ and v µ respectively, and compute the classical Kolmogorov-Smirnov distance :

T µ = 1 N sup 1≤k≤N |F (k) -G(k)|. ( 10
)
As F and G are identical for an exponentially i.i.d. random series, we then expect T µ to reach its minimum for the value of µ that gives the best estimate R µ (t) of R(t) :

µ = argmin µ T µ R = R µ . (11)
Plots of Figure 5 show the evolution of the Kolmogorov-Smirnov distance corresponding to the traces displayed in Figure 4. In the 3 cases, T µ clearly attains its minimum bound for µ close to to the actual value. The corresponding estimated processes R(t) derived from Eq. ( 11) match fairly well the real evolution of the (R) class in our model (see Figure 6). Propagation parameters β and l. According to our model, the arrival rate λ(t) of new viewers is given by Eq. ( 2). It linearly depends on the current number of active and past viewers. So, from the observation I(t) and the reconstructed process R(t) of Eq. ( 11), we could formally apply the maximum likelihood Eq. ( 7) to estimate β. In practice however, we have to bear in mind that : (i) the arrival process of rate λ(t) comprises a spontaneous viewers ingress that is governed by parameter l and which is independent of the current state of the system ; (ii) depending on the current hidden state of the model (buzz-free versus buzz state), it is alternately β = β 1 and β = β 2 that fix the propagation rate in Eq. ( 2). We designed an estimation procedure based on a weighted linear regression, that simultaneously addresses these two issues. We decompose our process R(t) (blue curve) with the estimated process R(t) (red curve) derived from expression [START_REF] Li | Performance model of interactive videoon-demand systems[END_REF]. rationale in two steps : First, let us consider the buzz-free state only and β = β 1 . As discussed in the estimation of µ the inter-arrival time w between the consecutive arrivals of two new viewers is an exponentially distributed random variable such that E (w|

R(t), R(t)
I(t) + R(t) = x) = (β x + l) -1 .
Concretely then, for different values of the sum I(t)+ R(t), we calculate the conditional empirical mean :

Ω(x) = 1 |I(x)| tn∈I(x) w n : I(x) = {t n : I(t n ) + R(t n ) = x}. ( 12
)
The linear regression of (Ω(x)) -1 against x yields at one go, both parameters estimation β (slope) and l (intercept).

Let us now return to the general form of our model with alternation of buzz and buzz-free periods. In the buzz-free case, β = β 1 corresponds to a normal workload activity, meaning that the sum I(t) + R(t) takes on rather moderate values. Conversely, when the system undergoes a buzz, β = β 2 and the population I(t) + R(t) suddenly increases to reach significantly larger values. Yet, in both cases, the quantity Ω -1 defined in Eq. ( 12) remains linear with x but with two different regimes (slopes) depending on the amplitude of I(t)+ R(t) = x. As a result, it is possible to reduce the bias that β 2 causes on the estimation of β 1 , using a weighted linear regression of Ω -1 vs x where the weights p(x) are proportional to the cardinal of the indicator sets I(x). Indeed, |I(x)| should be smaller for larger values of x because buzz episodes are expected to be less frequent than nominal activity periods. Figure 7 confirms the claim : the plots (x, Ω -1 ) show a manifest linear trend with higher variability at x's large, meaning a fewer terms entered the sum of Eq. ( 12). Formally, we can apply the exact same procedure to estimate β 2 , but considering opposite weights to favor the large values of x's. However, due to the large fluctuations of (Ω(x)) -1 in the corresponding region, the slope β 2 is subject to a very poor estimation variance. Instead, we propose to apply the ML estimator described in Eq. ( 7) on the restriction of I(t) to the buzz periods only. Strictly speaking, we should consider R(t) as well, but since a buzz event normally occurs on very small interval of time, we assume that R(t) (resp. R(t)) remains constant in the meanwhile (flash crowd viewers will enter in R compartment only after the visualization time). In practice, to automatically identify the buzz periods, we threshold I(t) and consider only the persistent increasing parts that remain above the threshold.

Transition rates a 1 and a 2 . As we already said, at time t, the inter-arrival time w separating to new incomers is a random variable drawn from an exponential law of parameter λ = β(i + r) + l, where I(t) + R(t) = i + r and β is either equal to β 1 or to β 2 . We denote f 1 (w) and f 2 (w) the corresponding densities built upon the reconstructed process R(t) and the estimated parameters ( β 1 , l) and ( β 2 , l) respectively. For a given inter-arrival time w = w n observed at time t n , we form the likelihood ratio f 2 (w n )/f 1 (w n ) to determine whether the system is in buzz or in buzzfree state. Moreover, in order to avoid non-significant state transitions we resort to a restoration method inspired by the Viterbi algorithm [START_REF] Kleinberg | Bursty and hierarchical structure in streams[END_REF]. Once we have identified the hidden states of the process, we estimate the transitions rates a 1 and a 2 from the average times spent in each state.

Numerical Validation

To evaluate the statistical performance of our estimation procedure, we resort to numerical experiments to empirically get the first and the second order moments of each parameter estimator. Owing to the versatility of our model, we must ensure that the proposed calibration algorithm performs well for a variety of workload dynamics. To this end, we systematically reproduce the experiments considering the 3 sets of parameters reported in Table 1. For each configuration, we generate 10 independent realizations of processes similar to the ones depicted in Figure 4, and use these to derive descriptive statistics.

The box-and-whisker plots of Figure 8 indicate for each estimated parameter (centered and normalized by the corresponding actual value) the sample median (red line), the inter-quartile range (blue box height) along with the extreme samples (whiskers) obtained from time series of length 2 21 points each. As expected (owing to the maximum likelihood procedure), estimation of γ shows to be the most accurate, both in terms of bias and variance. But more surprisingly though, although the estimation β 1 derives from a heuristic procedure that itself depends on the raw approximation R(t) of Eq. ( 9), the resulting performance is remarkably good : bias is always negligible (less than 5% in the worst case (c)) and the variance always confines to 10% of the actual value of β 1 . Notice also that the estimation of β 1 goes from a slight underestimation in case (a) to a slight overestimation in case (c), as the buzz effect, i.e. the value of β 2 , grows from traces (a) to (b). Compared to β 1 , the estimation of β 2 behaves more poorly and proves to be the most difficult parameter to estimate. But we have to keep in mind that this latter is only based on buzz periods which represent only a small fraction of the entire time series. Regarding the parameter µ, its estimation remains within a 20% inter-quartile range but cases (a) and (c) show a systematic bias (median hits the lower quartile bound). Let us then recall that the procedure, described by Eq. [START_REF] Kanrar | Analysis and implementation of the large scale video-on-demand system[END_REF] to determine µ selects within some discretized interval, the value of µ that yields the best T µ score. It is then very likely that the true value does not coincide with any sampled point of the interval and therefore, the procedure picks the closest one that systematically lies beneath or above. Finally, estimation of the transition parameters a 1 and a 2 between the two hidden states relies on all other parameters estimation, cumulating so all relative inaccuracies. Nonetheless and despite a systematic underestimating trend, precision remains within a very acceptable confidence interval. Convergence rate of the empirical estimators is another important feature that binds case (a) case (b) case (c) the estimate precision to the amount of available data. Using the same data set, the bar plots of Figure 9 depicts the evolution of the mean square error MSE( θ) = E{( θ -θ) 2 } -where generic θ stands for any parameter of the model -with the length N of the observable time series I(t).

-0.2 -0.1 0 0.1 0.2 0.3 β 2 µ l γ a 2 a 1 β 1 -0.3 -0.2 -0.1 0 0.1 0.2 γ l a 1 a 2 µ β 1 β 2 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 β 1 β 2 γ µ l
As our purpose is to stress the rate of convergence of these quantities towards zero, to ease the comparison, we normalize the MSE of each parameter by its particular value at maximum data length (i.e 2 21 points here). Then, the estimator rate of convergence α θ corresponds to the decaying slope of the MSE with respect to N in a log-log plot, i.e. MSE( θ) ∼ O(N -α θ ). For the different parameters of our model we obtain convergence rates that lie between α β1 = 0.9 and α a2 = 0.2, leading each time to sub-optimal convergence (α θ < 1). It is worth noticing that, despite its relatively ad hoc construction, the estimator of β 1 has an almost optimal convergence rate, which proves the rationality of our approach. 1.

Validation of the estimation procedure against a real workload trace

We now apply the calibration procedure detailed in the previous section, to fit our model on real data and to assess the data-model adequacy in the context of highly volatile workloads. As a paradigm for variable demand applications, we use a VoD trace, released by the Greek Research and Technology Network (GRNET) VoD servers [START_REF] Grnet | Video traces obtained from grnet[END_REF]. Since the trace shows modest activity with a workload level that is not sufficient to stress-test our system, we scale up the data, shrinking all inter-arrival times by a factor of 10. The resulting workload time series displayed in Figure 10, clearly shows two distinct periods of steady activity before and after the time index t = 200. We consider the sub-series on both sides of this cutoff time, as two individual workload processes, referred to as trace I and trace II respectively, and we calibrate our model on each of them separately.

Results of the parameters estimation are reported in Table 2, and we verified that in both cases the stability condition of Eq. ( 5) was satisfied. In the same vein, we also compared the empirical mean workload of each trace with its corresponding theoretical value given by the formula (4). We obtain for trace I a relative difference of 12% (E(i) = 5.59 compared to i = 4.99), and of 12.5% for trace II (E(i) = 0.621 compared to i = 0.71). Naturally, the correspondence here is not as striking as it is with the synthetic traces of Section 3. But we must bear in mind that first, ab initio nothing guarantees that the underlying system matches our model dynamics and, second, traces I and II can possibly encompass short scale non-stationary periods (e.g. day versus night activity) which are not accounted for in our model. Notwithstanding this, the match we observe is quite satisfactory and we now focus on higher order statistics to further characterize the data model adequacy. As we do not have a closed-form for the steady state distribution of our Markov process model, nor we have an analytic expression for its autocorrelation function, we use the two sets of estimated parameters of Table 2 to synthesize two time series that we compare to the real workload traces I and II. We refer to those synthetic traces as to the fitted traces I and II. The plots in Figure 11 of current viewers along time. But also, and perhaps more importantly, it demonstrates that the dynamical mechanism underlying our constructive model is able to perfectly reproduce the temporal structure of the real traces, by imposing the correct statistical dependencies between distant samples I(t) and I(t + τ ) of the process.

In addition to serve as a mean to evaluate the goodness-of-fit of our model, the estimated parameters bring on their own, a valuable insight about the system itself. For instance let us compare the propagation rates β 1 and l estimated from traces I and II, successively. In the first case, β 1 < l, meaning that arrival of new viewers is dominated by spontaneous incomers and is not so much due to information propagation through gossip. Conversely, β 1 > l for the second workload regime, indicating that the spontaneous attraction of the server has severely dropped whereas the peer-to-peer diffusion component significantly increased but not sufficiently to sustain the mean workload activity. At the same time, the index µ tripled, meaning that the mean memory period for propagation shrank by a factor of 3. This parameter could then be used as an indicator of the content interest delivered by the server, and of its lifetime in users mind.

Large Deviation Principle and its interpretation

Consider a continuous-time Markov process (X t ) t≥0 , taking values in a finite state space S, of rate matrix A = (A ij ) i∈S,j∈S . In our case X is a vectorial process X(t) = (N I (t), N R (t)) , ∀t ≥ 0, and

S = {0, • • • , I max } × {0, • • • , R max }.
If the rate matrix A is irreducible, then the process X admits a unique steady-state distribution π satisfying πA = 0. Moreover, by Birkhoff ergodic theorem, it is known that for any mapping Φ : S → R, the sample mean of Φ(X) at scale τ , i.e. 1/τ • τ 0 Φ(X s )ds converges almost-surely towards the mean of Φ(X) under the steady-state distribution, as τ tends to infinity. The function Φ is often called the observable. In our case, as we are interested in the variations of the current number of users N I (t), Φ will simply be the function that selects the first component : Φ(N I (t), N R (t)) = N I (t). The large deviations principle (LDP), which holds for irreducible Markov processes on a finite state space [START_REF] Varadhan | Large deviations[END_REF], gives a efficient way to estimate the probability for the sample mean calculated over a large period of time τ to be around a value α ∈ R that deviates from the almost-sure mean :

lim ǫ→0 lim τ →∞ 1 τ log P τ 0 Φ(X s )ds ∈ [α -ǫ, α + ǫ] = f (α). (13)
The mapping α → f (α) is called the large deviations spectrum (or the rate function). For a given function Φ, it is possible to compute the theoretical large deviations spectrum from the rate matrix A as follows. One first computes, for each values of q ∈ R, the quantity Λ(q) defined as the principal eigenvalue (i.e., the largest) of the matrix with elements A ij + qδ ij Φ(j) (δ ij = 1 if i = j and 0 otherwise). Then the large deviations spectrum can be computed as the Legendre transform of Λ :

f (α) = sup q∈R {qα -Λ(q)} , ∀α ∈ R. (14)
As described in Equation [START_REF] Niu | Understanding demand volatility in large vod systems[END_REF], α τ = i τ corresponds in our study case, to the mean number of users i observable over a period of time of length τ and f (α) relates to the probability of its occurrence as follows :

P{ i τ ≈ α} ∼ e τ •f (α) . ( 15
)
Interestingly also, if the process is strictly stationary (i.e. the initial distribution is invariant) the same large deviation spectrum f (•) can be estimated from a single trace, provided that it is "long enough" [START_REF] Barral | Large deviations for the local fluctuations of random walks[END_REF]. We proceed as follows : At a scale τ , the trace is chopped into k τ intervals {I j,τ = [(j -1)τ, jτ [, j = 1, . . . , k τ } of length τ and we have (almost-surely), for all α ∈ R :

f τ (α, ǫ τ ) = 1 τ log # j : Ij,τ Φ(X s )ds ∈ [α -ǫ τ , α + ǫ τ ] k τ and lim τ →∞ f τ (α, ǫ τ ) = f (α). (16)
In practice, for the empirical estimation of the large deviations spectrum, we use a similar estimator as the one derived in [START_REF] Barral | On the estimation of the large deviations spectrum[END_REF] and also used in [START_REF] Loiseau | Modeling TCP throughput : an elaborated large-deviations-based model and its empirical validation[END_REF]. At scale τ , we compute for each q ∈ R the values of Λ ′ τ (q) and Λ ′′ τ (q), where Λ τ (q) = τ -1 log k -1 τ kt j=1 exp q Ij,τ Φ(X s )ds . Then, for each value of τ , we count the number of intervals I j,τ verifying the condition in expression [START_REF] Barrat | Dynamical Processes on Complex Networks[END_REF] and estimate the scale-dependant empirical log-pdf f τ (α, ǫ τ ), with the adaptive choices derived in [START_REF] Barral | On the estimation of the large deviations spectrum[END_REF] :

α τ = Λ ′ τ (q) and ǫ τ = -Λ ′′ τ (q) τ . (17)
Let us now illustrate the LDP in the context of the specific VoD use case, where X would correspond to (i, r), the bi-variate Markov process. Φ(X) is i, the observable and τ 0 Φ(X S ) ds = i τ corresponds to the average number of users with a period τ .

Numerical Interpretations

For ease of computation we estimate the Large Deviation Spectrum for cases where I max = 30, R max = 60. We also choose the parameters accordingly (so that it does not saturate with the maximum value) for buzz and buzz-free scenarios. For the first case β 1 = 0.1, β 2 = 0.8, γ = 0.7, µ = 0.3, l = 1.0, a 1 = 0.006 and a 2 = 0.6. For the buzz-free case : β 1 = β 2 = β = 0.1, γ = 0.7, µ = 0.3, l = 1.0. Intrinsically, Large Deviation Principle naturally embeds the time scale notion into the statistical description of the aggregated observable at different time resolutions. As expected, the theoretical LD spectra displayed in Figure 12(a) reach their maximum for the same mean number of users. This apex is the almost sure value as described in Section ??. As the name suggests almost sure workload (α a.s ) corresponds to the mean value that we almost surely observe on the trace. More interestingly though, the LD spectrum corresponding to the buzz case, spans over a much larger interval of observable mean workloads than that of the buzz-free case. This remarkable support widening of the theoretical spectrum shows that LDP can accurately quantify the occurrence of extreme, yet rare events. Plots (b)-(c) of Figure 12 compare theoretical and empirical large deviation spectra obtained for the two traces. For each given scale (τ ) the empirical estimation procedure yields one LD estimate. These empirical estimates at different scales superimpose for a given range of α. This is reminiscent of the scale invariant property underlying the large deviation principle. If we focus on the supports of the different estimated spectra, the larger the time scale τ is, the smaller becomes the interval of observable value of α. This is coherent with the fact that for a finite trace-length the probability to observe a number of current viewers, that in average, deviates from the nominal value (α a.s ) during a period of time (τ ) decreases exponentially fast with τ . To fix the ideas, the estimates of plot (c), indicate that for a time scale τ = 400 sec., the maximum observable mean number of users is around 5 with probability 2 400•(-0.02) ≈ 35.10 -5 (point A), while it increases up to 9 with the same probability (2 100•(-0.08) ) for τ = 100 sec. (point B).

f (α) f (α) f (α) α = i τ α = i τ α = i τ

Resource management policies

Retuning to our VoD use case, we now sketch two possible schemes for exploiting the Large Deviation description of the system to dynamically provision the allocated resources :

-Identification of the reactive time scale for reconfiguration : Find a relevant time scale that realizes a good trade-off between the expectable level of overflow associated to this scale and a sustainable opex cost induced by the resources reconfiguration needed to cope with the corresponding flash crowd. -Link capacity dimensioning : Considering a maximum admissible loss probability, find the safety margin that it is necessary to provision on the link capacity, to guarantee the corresponding QoS.

Identification of the reactive time scale for reconfiguration

We consider the case of a VoD service provider who wants to determine the reactivity scale at which it needs to reconfigure its resource allocation. This quantity should clearly derive from a good compromise between the level of congestion (or losses) it is ready to undergo, i.e. a tolerable performance degradation, and the price it is willing to pay for a frequent reconfiguration of its infrastructure. Let us then assume that the VoD provider has fixed admissible bounds for these two competing factors, having determined the following quantities :

α * > α a.s. : the deviation threshold beyond which it becomes worth (or mandatory) considering to reconfigure the resource allocation. This choice is uniquely determined by a capex performance concern. σ * : an acceptable probability of occurrence of these overflows. This choice is essentially guided by the corresponding opex cost. Let us moreover suppose, that the LD spectrum f (α) of the workload process was previously estimated, either by identifying the parameters of the Markov model used to describe the application, or empirically from collected traces. Then, recalling the probabilistic interpretation we surmised in relation [START_REF] Gusella | Characterizing the variability of arrival processes with indexes of dispersion[END_REF], the minimum reconfiguration time scale τ * for dynamic resource allocation, that verifies the sought compromise, is simply the solution of the following inequality :

τ * = max{τ : P{ i τ ≥ α * } = ∞ α * e τ fτ (α) dα ≥ σ * }, (18)
with f τ (α) as defined in expression [START_REF] Barrat | Dynamical Processes on Complex Networks[END_REF]. From a more general perspective though, we can see this problem as an underdetermined system involving 3 unknowns (α * ,τ * and σ * ) and only one relation [START_REF] Jammalamadaka | Testing exponentiality by comparing the empirical distribution function of the normalized spacings with that of the original data[END_REF]. Therefore, and depending on the sought objectives, we can imagine to fix any other two of these variables and to determine the resulting third so that it abides with the same inequality as in expression (18).

Link capacity dimensioning

We now consider an architecture dimensioning problem from the infrastructure provider perspective. Let us assume that the infrastructure and the service providers have come to a Service Level Agreement (SLA), which among other things, fixes a tolerable level of losses due to link congestion. We start considering the case of a single VoD server and address the following question : What is the minimum link capacity C that has to be provisioned such that we meet the negotiated QoS in terms of loss probability ? Like in the previous case, we assume that the estimated LD spectrum f (α) characterizing the application has been priorly identified. A rudimentary SLA would be to guarantee a loss free transmission for the normal traffic load only : this loose Inria QoS would simply amount to fix C to the almost sure workload α a.s. . Naturally then, any load overflow beyond this value will result in goodput limitation (or losses, if there is no buffer to smooth out exceeding loads). For a more demanding QoS, we are led to determine the necessary safety margin C 0 > 0 one has to provision above α a.s. to absorb the exact amount of overruns corresponding to the loss probability p loss that was negotiated in the SLA. From the interpretation of the large deviation spectrum provided in Section ??, this margin C 0 is determined by the resolution of the following inequality :

C 0 : ∞ αa.s.+C0 τmax τmin e τ •f (α) dτ dα ≤ p loss : ∞ αa.s.+C0 e τmax•f (α) -e τmin•f (α) f (α) dα ≤ p loss (19)
In this expression, τ min is typically determined by the size Q of the buffers that is usually provisioned to dampen the traffic volatility. In that case,

τ min = Q α -(α a.s. + C 0 ) , (20)
corresponds to the maximum burst duration that can be buffered without causing any loss at rate α > C = α a.s. + C 0 . As for τ max , it relates to the maximum period of reservation dedicated to the application. Most often though, the characteristic time scale of the application exceeds the dynamic scale of flash crowds by several orders of magnitude, and τ max can then simply be set to infinity. With these particular integration bounds, Equation ( 19) simplifies to

C 0 = C -α a.s. : ∞ C -1 f (α) e Q α-C f (α) dα ≤ p loss , (21)
a decreasing function of C, which can be solved using a simple bisection technique.

As long as the server workload remains below C, this resource dimensioning guarantees that no loss occurs. All overrun above this value will produce losses, but we ensure that the frequency (probability) and duration of these overruns are such that the loss rate remains conformed to the SLA. The proposed approach clearly contrasts with resource over-provisioning that does not seek at optimizing the capex to comply with the loss probability tolerated in the SLA. The same provisioning scheme can straightforwardly be generalized to the case of several applications sharing a common set of resources. To fix the idea, let us consider an infrastructure provider that wants to host K VoD servers over the same shared link. A corollary question is then to determine how many servers K can the fixed link capacity C support, while guaranteeing a prescribed level of losses. If the servers are independent, the probability for two of them to undergo a flash crowd simultaneously is negligible. For ease and without loss of generality, we moreover suppose that they are identically distributed and modeled by the same LD spectrum f (k) (α) = f (α) with the same nominal workload α (k) a.s. = α a.s. , k = 1, . . . K. Then, following the same reasoning as in the previous case of a single server, the maximum number K of servers reads :

K = arg max K (C -K • α a.s. ) ≤ C 0 , (22)
where the safety margin C 0 is defined as in expression [START_REF] Varadhan | Large deviations[END_REF]. Then, depending on the agreed Service Level Agreements, the infrastructure provider can easily offer different levels of probability losses (QoS) to its VoD clients, and adapt the number of hosted servers, accordingly.

Conclusion

Many applications deployed on a cloud infrastructure, such as a Video on Demand service, are well known for undergoing highly volatile demand, making their workload hard to qualitatively and quantitatively characterize. Adopting a constructive approach to capture the VoD users' behavior, in this report we proposed a simple, concise and versatile model for generating the workload variations in such context. We also devised an heuristic identification procedure that aims at estimating the parameters values of the model from a single collected trace. First, we numerically evaluated the accuracy of this procedure using several synthetic traces. Our experiments show that the procedure introduces little bias and typically recovers the actual parameters value with a relative error of about 10%. Second, we apply this same procedure against two real video server workload traces. Obtained results demonstrate that, once the model has been calibrated, it succeeds to reproduce the statistical behavior of the real trace (in terms of both the steady-state probabilities and the autocorrelations for the workload time series). Moreover, owing to the constructive nature of our model, the estimated values of the parameters provide valuable insight on the application that it would be difficult, or even impossible, to infer from the raw traces. The captured information may answer questions of practical interest to cloud oriented service providers, like : is the application workload mostly driven by spontaneous behaviors, or is it rather subject to a gossip phenomenon ? Furthermore, a key-point of this model is that it permits to reproduce the workload time series with a Markovian process, which is known to verify a Large Deviation Principle (LDP). This particularly interesting property yields a large deviation spectrum whose interpretation enriches the information conveyed by the standard steady state distribution : For a given observation (workload trace), LDP allows to infer (theoretically and empirically) the probability that the time average workload, calculated at an arbitrary aggregation scale, deviates from its nominal value (i.e. almost sure value).

We leveraged this multiresolution probabilistic description to conceptualize two different management schemes for dynamic resource provisioning. As explained, the rationale is to use large deviation information to help network and service providers together to agree on the best capexopex trade-off. Two major stakes of this negotiation are : (i) to determine the largest reconfiguration time scale adapted to the workload elasticity and (ii) to dimension VoD server so as Inria to guarantee with upmost probability the Quality of Service imposed by the negotiated Service Level Agreement. More generally though, the same LDP based concepts can benefit any other "Service on Demand" scenarii to be deployed on dynamic cloud environments.

Figure 2 -

 2 Figure 2 -Basic schematics of a VoD system with transcoding/caching servers reserve a sufficient amount of server outgoing bandwidth to sustain continuous media delivery.When multiple VoD providers (such as Netflix) are on board to use cloud services from cloud providers, there will be a market between VoD providers and cloud providers, and commodities to be traded in such a market consist of bandwidth reservations, so that VoD streaming performance can be guaranteed.As a buyer in such a market, each VoD provider can periodically make reservations for bandwidth capacity to satisfy its random future demand. A simple way to achieve this is to estimate expectation and variance of its future demand using historical demand information, which can easily be obtained from cloud monitoring services. As an example, Amazon Cloud-Watch provides a free resource monitoring service to Amazon Web Service customers for a

2 Figure 3 -

 23 Figure 3 -Markov chain diagram representing the evolution of the Current viewers (i) and Past Viewers (r) populations with a Hidden Markov Model.

Figure 4 -

 4 Figure 4 -Illustration of our model ability at generating different dynamics of workload I(t). See

Figure 5 -

 5 Figure 5 -Evolution of the exponential test statistics (10) applied to the traces of Figure 4. Dotted vertical lines locate the actual value of µ for each case ; dot markers on each curve indicate the estimated value µ corresponding to the minimum point of the statistical test Tµ.

Figure 6 -

 6 Figure 6 -Evolution of the number of active past viewers. Comparison of the actual (non observable)

Figure 7 -

 7 Figure 7 -Weighted linear regression of Ω -1 vs x corresponding to the three traces of Figure 4. Superimposed are the linear trends fitted on the respective data.

Inria

Figure 9 -

 9 Figure9-Evolution of the Mean Square Error versus the data length N in a log-log plot. For the sake of conciseness, we only show here the results corresponding to the case (b) of Table1.

Figure 10 -

 10 Figure 10 -Real workload time series corresponding to a VoD server demand from [20]. Initial trace was scaled up by a factor of 10 to increase the mean workload. Trace is chopped into two separate processes (Trace I and II) corresponding to different activity levels.

Figure 11 -

 11 Figure 11 -Comparison of the empirical steady-state distribution and of the autocorrelation function of the real (blue curves) and the fitted (red curves) traces. Top two plots correspond to trace I and bottom plots correspond two trace II.

Inria

Figure 12 -

 12 Figure 12 -Large Deviations spectra corresponding to the traces of Figure 4. (a) Theoretical spectra for the buzz free (blue) and for the buzz (red) scenarii. (b) & (c) Empirical estimations of f (α) at different scales from the buzz free and the buzz traces, respectively.

Figure 13 -

 13 Figure 13 -Dimensioning K, the number of hosted servers sharing a fixed capacity link C. The safety margin C0 is determined according to the probabilistic loss rate negotiated in the Service Level Agreement between the infrastructure provider and the VoD service provider.

Table 1 -

 1 Parameters value used in the workload model to generate the three traces plotted in Fig.4. The last two rows correspond to the theoretical mean workload of Eq. (4) and to the sample mean value estimated from the traces.

		case (a)	case (b)	case (c)
	β 1	4.762 × 10 -4 3.225 × 10 -5 2.439 × 10 -5
	β 2	0.0032	0.0032	0.0032
	γ	0.0111	0.0020	0.0011
	µ	5 × 10 -4	3.289 × 10 -5	2.5 × 10 -5
	l	10 -4	10 -4	10 -4
	a 1	10 -7	10 -7	10 -7
	a 2	0.0667	0.0667	0.0667
	E(i)	1.92	15.68	44.72
	Emp. mean i	1.74	16.72	45.23

Table 2 -

 2 Estimated Parameters from traces I and II separately.

		β 1	β 2	γ	µ	l	a 1	a 2
	I	0.0013	0.0084	0.0039	0.0028	0.0032	3.13.10 -4	0.022
	II	0.0049	0.0183	0.0118	0.0095	0.0005	1.32.10 -5	0.041
								Inria

RR n°8072

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Aurélie Villard

Luc Petit

Alan Lelah

Daniel Brissaud

Eco-design for microelectronic products: measurement and integration of environmental criteria

Keywords: eco-design, microelectronic product, life cycle assessment, chips packaging

Introduction

In a context where more and more micro-products are embedded in day-to-day applications, designers of microelectronic components are dealing with higher and higher technologic stakes linked to circuits' miniaturization and energetic efficiency research. That trend justifies that an effort be made in the microelectronic industry to make chips cleaner. Before envisaging designing chips taking into account environmental considerations, the designer must have indications on the environmental performances of the different manufacturing processes that could be used. This paper, illustrated by a case of study focused on the chip packaging, describes how to characterize environmental issues so as to influence a decision to be taken during the development of the manufacturing process involved. Chips are usually encapsulated in a package prior to being embedded into electronic system. This package performs four key functions:

1. Provide a sturdy set of leads or balls that allow a chip to be connected to the system in which it will operate 2. Provide physical protection for the chip against breakage or scratching 3. Provide environmental protection against damage from chemicals, moisture or gases 4. Dissipate the heat generated by the chip during operation. The package is manufactured with a process flow made of about fifteen steps. It begins when the finished silicon wafer is split up into individual die. The die is then assembled into the package which can be handled in the final application. Awareness on environmental performances of packaging-technology is a key issue for microelectronic industry for three main reasons. First, although the technological complexity lies mainly in silicon treatment, package design has a significant impact on the performance and reliability of the chip(s) it house and its business cost can be a large fraction of the total manufacturing cost [START_REF] Andrae | Life cycle assessment of integrated circuit packaging technologies[END_REF]]. The technological roadmaps aim to reduce the mass and volume and add the possibility to encapsulate several chips on a same package. These roadmaps drive designers to always imagine innovative solutions to ensure a solid and minimalist protection for more and more complex and fragile chips. Secondly, with two European directives WEEE and RoHS, legislative pressure is focused on package. Indeed the package represents 95% of the mass of the chip. In the last decade, compliance to RoHS and so the search for health neutrality is leading to high technical challenges for package designers. In a study of the shift from Tin/Lead to lead-free solders paste, [START_REF] Ecoinvent | Attributional and Consequential Environmental Assessment of the Shift to Lead-Free Solders[END_REF]] contributed to demonstrate that RoHS compliance process is not likely to result on environmental improvement and that LCA helps to complete the knowledge on the consequences of such a shift through a life cycle approach. Finally, the environmental impact of a chip is far to from being negligible compared to the semiconductor component. According [START_REF] Higgs | Developing an overall CO2 footprint for semiconductor products[END_REF]], the package substrate is the top supply chain material impacting on carbon emissions. The package manufacturing part is such a big consumer of metals, plastics and organic compounds that its impact on the environment is absolutely not negligible.

Eco-design tool description

Usually, concerns of package designer are driven by different constraints and must reach a compromise. The reasons for change usually come from standard industrial values: cost, quality and efficiency. At designer-level, they are materialized by electrical performances, volume, weight, testability, reliability, choice of materials (avoid banned substances) and cost. The environmental criterion is not measured and so cannot weigh in decision-making processes. Eco-design for package is a conceptual approach in which the designer can consider the environmental impacts of the package he designs. The constraints linked to environmental impacts become a decisive factor to be systematically integrated in the earlier stages of packaging technology development.

Properties of the eco-design tool

Methodology for the development There have been several previous approaches [START_REF] Krishnan | Quantifying the Environmental Footprint of Semiconductor Equipment Using the Environmental Value Systems Analysis (EnV-S)[END_REF]], [START_REF] Griese | Sustainable development of microelectronic technology processes integration of ecodesign[END_REF]], [START_REF] Somani | Environmental Impact Evaluation Methodology for Emerging Silicon-Based Technologies[END_REF]], [START_REF] Harland | Design for Environment in New Semiconductor Manufacturing Technologies[END_REF]] to develop environmental design and evaluation tools for the microelectronic industry. The direct duplication of such methods for a specific industrial context poses several problems because of their theoretical complexity, the nature of the necessary information or the in-house organization of work in the company. So the deployment of an original tool for eco-design of packaging technologies has been preferred: it has been built directly from the collaboration with concerned design teams. A user-centred methodology was applied regarding designers' specificities, and mainly the data they can provide in the earlier stages of development. Several working groups have been organized to understand designers' expectations. An opinion survey helped to define the outline of the results they judge to be useful and easy-to-use. There was so a deliberate choice to exclude any other considerations than environment (economic, technical, or quality) for a better legibility of the results for this pilot project. It gives a tool admittedly using non-original concepts but really in line with STMicroelectronics designer's expectations and adapted to packages specificities. The next section describes the tool which is an interface for simulating design options during development of a packaging technology.

Functionalities

Eco-design tool for package designers is applied on new packages under development that are not yet real materialized. It tends to consider the impact generated by a package through its life cycle by: 1. Measuring it: estimate if there is benefit or loss when a change occurs in the process flow that means weighing the environmental footprint of a design decision 2. Controlling it: ensure that the impact is not getting worse compared to previous generations 3. Improving the impact: underline weak points. If there are several likely options, they are compared and the best option available is promoted. In order to fulfil the 3 goals, the first stage consists of acquiring a better knowledge of the effect on the environment of the process change. Life cycle assessment methodology (LCA) was selected as the support for marking design option regarding a large spectrum of environmental impacts [Bhander 2004]. LCA is acknowledged worldwide as a powerful technique which aims at presenting a global view of the potential impacts derived from the product under consideration, over its full life cycle. Since the last decade, environmental analysis of chips [START_REF] Williams | Environmental impacts of microchip manufacture[END_REF]], [START_REF] Boyd | Life-cycle assessment of computational logic produced from 1995 through 2010[END_REF]] has been progressing toward life cycle LCA. However, it is important to keep in mind that main challenges still remain in developing LCA tools to improve data modelling, keep pace with technology changes and handle uncertainty [START_REF] Plepys | The environmental impacts of electronics. Going beyond the walls of semiconductor fabs[END_REF]]. Seven indicators, called µ-KEPIs for "Key Environmental Performances Indicators for microelectronic industry", related to the main concerns of microelectronic industrial sector have been recently developed by STMicroelectronics and are listed in table 1. The main contributions to impact categories (emissions, raw materials, chemical substances, water or electricity) are computable. In the results display screen, Pareto graphs are available for each µ-KEPI: the ten most damaging areas, the ten most damaging elementary steps, the sources which together count for more than 90% of the damage. When a new process step is inspired from an existing one, the designer can also compare to this parent step.

LCA implementation

The environmental evaluation is applied separately to the process steps used in package manufacturing. The impact of the package is the sum of the impact of each step. However, to tackle the difficulty to model a complex system such as a microelectronic product, a specific scope for the LCA is required. [START_REF] Williams | The 1.7 Kilogram Microchip: Energy and Material Use in the Production of Semiconductor Devices[END_REF]] developed the concept of the secondary materialization. He proposed that "increasingly complex products require additional secondary materials and energy to realize their lower entropy form". Contrarily to conventional products, the model of the chip must focus on the secondary raw materials used by the manufacturing processes rather than on its material structure. For package manufacturing, the two levels, direct and indirect, are likely to be modelled with a bottom-up approach. Practically, material and energetic flows, in-coming or out-coming, were inventoried and the environmental analysis is assessed with SimaPro, a LCA-dedicated software.

Chip life cycle starts when raw materials and energetic resources are extracted from the soil and ends with the disposal of the chip. As illustrated in the figure 1, the scope of LCA is limited to manufacturing including supply chain and the elimination of the substances which remain on the chip. A basic average end of life scenario has been created for chips eliminated in 2011 with different destinations (incineration, landfill or recycling channels). The end of life phase of the package is considered in a fragmented way by considering individually each material which composes the chip. So effects of end of life are allocated to the process step where the material is added. For instance, effect of gold recycling are attributed to wire bonding and effects of polymers incineration are attributed to moulding process.

Environmental evaluation of a manufacturing process by R&D process engineering

LCA of a new process

The case presented in this section studies a modification of the manufacturing conditions in wire bonding that seeks to largely decrease the cost of the operation: gold wire is replaced by similar copper-made wire (length and thickness). The benefits or losses in terms of environmental performances are not a priori obvious. Due to heavier extraction processes, the environmental impact of gold is largely higher than copper one: the environmental footprint of the two metals are extracted from Eco-Invent database [START_REF] Ecoinvent | Attributional and Consequential Environmental Assessment of the Shift to Lead-Free Solders[END_REF] and summarized in Table 2. Water eco-toxicity kg TEG eq 2.34E+02

1.75E-02 Nevertheless, copper has different physical properties which have led to intrinsic process modifications: to avoid oxidation of the copper wire, a tool, called "copper kit" is stacked on the equipment to create a neutral ambiance. This kit induces extra consumption of electricity and N2/H2 gaseous mix is injected in the chamber during the process. The throughput in equipment is decreased, meaning extra equipment occupation and so an augmentation of electricity consumption. In this case, the relative comparison with the seven indicators highlighted that the increase of electricity and chemicals consumptions is offset by the selection of a greener metal. The Figure 2 presents on the 7 µ-KEPIs the comparative results of the two techniques of bonding.

Figure 2. Comparison of the environmental footprints of the two wire bonding techniques

Four environmental aspects are improved against only one which is getting worse: copper wire is preferable than a gold one. In selecting a more eco-friendly material, changing the package composition and modifying its processes, LCA stressed that the environmental pressure of a device is reduced. This result is all the more true that yield of copper process is going to be improved. In some cases, LCA results are equivocal regarding the seven indicators: the designer, who is not an environmental expert, is unable to choose the best option available. In that case the analysis is coupled with specific data treatment to compare it with in-house environmental tolerance thresholds. These thresholds, fixed in function of in-house strategic issues, allow normalization of the scores and then highlighting the most important aspects.

Quick-LCA for processes

In some cases, the recipe has already been elaborated but needs some updates to match the new package specificities. To do this, an option for assessing "Quick-LCA" has been specifically developed for each process. The µ-KEPIs can be deduced from referent package by integrating parametric models on process units [START_REF] Murphy | Development of Parametric Material, Energy, and Emission Inventories for Wafer Fabrication in the Semiconductor Industry[END_REF]]. The key parameters of the process step are identified and the functions which make vary the flows are defined. The engineer must just specify the modifications compared to a referent package.

Figure 3. Quick-LCA interface for wire bonding

In the wire bonding case (figure 3), if the parameter "wire length" is modified, in-coming and outcoming flows can be deduced from a similar recipe at the condition that other parameters are unchanged. Indeed, the length of wire to be bonded is directly proportional to both the mass of metal consumed and the duration of the process: the quantities needed for LCA (mass of metal, volume of gas N2/H2 and electricity) are easily calculated.

Experimentation of the eco-design tool

The functionalities of the tool are tested by a R&D design team during a design review. A design review is a punctual meeting for evaluating the progress of a project. The team is constituted with 8 process engineers and the project leader. The package ζ is under development: no prototype has been yet tested. The experimentation aims at measuring, by means of the eco-design tool, the future environmental impact of the package ζ once it will be manufactured and integrated in its final electronic application.

Preliminary studies to prepare the design review

The package ζ is designed in order to answer new technical specifications in terms of functionalities, and dimensions. It is requested for innovative complex integrated-circuits needing an extra-thin protection. Five packages (α, β, γ, δ and ε) have been identified by the product manager as potential former or parent technologies. The five packages have the same body size (12x12 mm) but different thickness, number of connections, balls diameter and connecting options (wire bonding or flip chip pads). The package ζ has mixed characteristics coming from these five pre-existing packages but its process assembly is more complex because of a thicker substrate, a more important number of electric connections and so a high density of electrical connections per square millimetre.

During the design review

Once the LCAs of the five packages and theirs recipes are recorded, the environmental evaluation of the package ζ can be easily assessed with R&D team help. First, the manufacturing process flow of the package is captured. The figure 4 describes the theoretical path of the questions asked in turn to the process engineer during the design review. The table 3 summarizes the data collected during the meeting according a theoretical path of questioning: the more likely flow, the maturity status of each process step (already existing, modified from a pre-existing recipe, new) and in the two first cases which are the parent packages. It appears that only three packages (α , β and δ) are indeed necessary to conduct the analysis on the package ζ. For the processes which have never been tested before, in this case laser opening of interconnections through organic resin, the process engineer in charge of the development is asked to estimate the parameters of the future recipe. All inputs and by-products in terms of fluids, materials, energy and emitted substances must be evaluated. An uncertainty warning is associated to this step and the scores are noted as temporary and will be re-evaluated once the recipe parameters will be better-known. For already existing processes, µ-KEPIs are either directly downloaded from parent package or "Quick-LCA" interface is used if there are different.

After the design review

Thanks to the tool, projected environmental footprint of a package was known several months before industrialization. Results are presented to the design team: the environmental footprint of the package

Conclusion and perspectives

The environmental evaluation of a packaging technology under-development has been successfully assessed thanks to the collaboration with a design team. It was part of the elaboration and test of a prototype tool presented in this paper. It is expected that this prototype tool will incline designers to systematically include environmental criteria into their decision-making process. This experiment was a first step of eco-design in the company and gave several perspectives before being systematically used by package designers. First of all, it has proven that environmental assessment of new package is possible. Mainly, the use of a simplifying hypothesis that design is mostly made by analogy has been confirmed.

In such an eco-design process, the collaboration and active participation of all the design team is required. To ensure viability and potential duplication at company-scale, the time allocated to this new task should be minor. That is why designers' constraints are carefully evaluated to be in line with their day-to-day tasks. First, it is advisable to further simplify process modelling in the tool: it will be possible both by asking for less data for the life cycle inventory and by reducing the number of key parameters in Quick-LCA. Secondly, the data to be collected for LCA of process recipes, except input quantities of raw materials and fluids, are not easily available for a R&D engineer because they have rather never been considered so far. To tackle this difficult, a full collaboration between R&D centres and production areas is required and should be established.

Even if the results have not been directly used and integrated for the design, for the first time, the designers were aware about environmental impact of their activities and that any design decisions is likely to modify the environmental footprint of the future chip. The interpretation of the results allows the whole design team to define trends and find axes for improvement. In a process of continuous improvement, each time a LCA of a package or a process step is assessed, the seven indicators will be recorded in a database. The more this database will grow, the more the range of processes used in the company will be covered and so faster will be the LCA of new package. This experience gave encouraging results and has validated a method able to arouse designers' interest about environmental aspects without any unreasonable constraints.

Figure 1 .

 1 Figure 1. Scope for package LCA Data related to recipes used to run the process are collected either by the R&D Process Engineer who has created the recipe or by the Production Process Engineer who manages industrialization of this recipe. Raw materials properties (density and composition) are collected with suppliers' help via the Material Declaration process or the Material Safety Datasheet. Some suppliers were able to give the impacts of their activity per kilogram or per production unit. Every raw material used during manufacturing is modelled into SimaPro.

Figure 4 .

 4 Figure 4. Theoretical path of the questions used to fill the interface

ζ

  compared to the previous generations and a Pareto with the relative environmental footprint of each manufacturing process step of the package ζ. Figures 5 and 6 show the type of results obtained in the case of the global warming indicator and with a relative scale taking package ζ as reference.

Figure 5 .

 5 Figure 5. Global Warming Potential for package ζ and its parent packages

Table 1 . Indicators considered in the tool -µ-KEPIs

 1

	Indicators	Unit	Sources
	Water consumption during package manufacturing (tape and pure water) cubic meters	x
	Electricity consumption by the equipments during package manufacturing mega Joules	x
	Abiotic depletion	kg Sb eq.
	Eutrophication Greenhouse effect	Kg PO 4 eq kg CO 2 eq.	CML
	Photochemical oxidation	kg C 2 H 2 eq
	Water eco-toxicity	kg TEG eq IMPACT 2002+

Table 2 . Environmental footprint of 1 meter-long wires made of gold and copper

 2

	Indicators	Unit	Gold wire Copper wire
	Abiotic depletion	kg Sb eq.	9.82E-04	4.81E-08
	Eutrophication	Kg PO 4 eq	1.12E-04	1.31E-08
	Greenhouse effect	kg CO 2 eq.	8.00E-02	6.29E-06
	Photochemical oxidation kg C 2 H 2 eq	2.08E-05	1.65E-08

Table 3 . Data collected during design review (reproduction of the interface)

 3

	More likely manufacturing flow ζ	Status	Interface	Parent package(s)
	Substrate manufacturing	Modified	Quick-LCA	β and δ
	Wafer bumping	New	New LCA	x
	Wafer back grinding	Modified	Quick-LCA	β or δ
	Wafer sawing	Modified	Quick-LCA	α or β or δ
	Top ball attach	Modified	Quick-LCA	α or β or γ or δ
	Cleaning	Already existing	Download	β or δ
	Substrate baking	Already existing	Download	α or β or γ or δ
	Flip chip and reflow	Modified	Quick-LCA	β or δ
	Substrate baking	Already existing	Download	α or β or γ or δ
	Plasma	Already existing	Download	β or δ
	Underfill dispensing and curing	Modified	Quick-LCA	β or δ
	Plasma	Already existing	Download	α or β or γ or δ
	Injection moulding and curing	Modified	Quick-LCA	α
	Laser treatment opening	New	New LCA	x
	Ball attach	Modified	Quick-LCA	δ
	Laser marking	Already existing	Download	α or β or γ or δ or ε
	Package singulation	Already existing	Download	α or β or γ or δ or ε

Acknowledgments

This research is funded by STMicroelectronics. The authors wish to thank all the members of the pilot design team, CPA Grenoble of STMicroelectronics, who takes an active part in the eco-design project for package designers.