Sylvain Deville

Guillaume Bernard-Granger

Influence of surface tension, osmotic pressure and pores morphology on the densification of ice-templated ceramics

Keywords: Porosity, Sintering, Ice-templating

We attempt here to clarify the influence of the nature of the solvent on the ice-templating process, based on literature experimental results of porosity and suspension concentration. In particular, we compare the relative influence of surface tension, osmotic pressure and the pore size and radius of curvature on the densification of the green bodies. Variations of the osmotic pressure when changing the solvent can yield significantly different behaviour during the freezing stage and therefore significantly affect the densification behaviour, an effect not apparent when using the Shanti model. In addition, and equally important, a modification of the macropores morphology is induced by the evolution of the solvent crystals morphology, and resulting in smaller and more intricate crystals and macropores which facilitates the densification. The radius of curvature is probably also an essential parameter controlling the densification of ice-templated ceramics.

Introduction

The ice-templating process has attracted an increasing interest over the last few years, in particular in the field of porous ceramics [START_REF] Deville | Freeze-Casting of Porous Ceramics: A Review of Current Achievements and Issues[END_REF] . The porosity in these materials is templated by the solvent crystals, which are subsequently removed by sublimation to obtain a green body, which in turn can be densified by conventional sintering methods. The process share many similarities with the existing ceramic processing routes -such as slip casting or pressing-from the point of view of particles packing in the green body. Nevertheless, some specificities of the process should be investigated with greater attention. In particular, the choice of the solvent, which can be made for pragmatic reasons such as a higher solidification temperature or the absence of specific equipment -e.g., a freeze-dryer-for solvent removal, can also have a great influence over the characteristics of the porous bodies obtained. The solvent play a key role on the redistribution and packing of particles which occur during the solidification step. Here, we attempt to clarify the influence of the nature of the solvent on the process, based on literature data. In particular, we compare the relative influence of surface tension and the pore size and radius of curvature on the densification of the green bodies.

Results and Discussion

The nature of the solvent is one of the many parameters which can be selected or adjusted in the ice-templating process. Water is the most commonly used, for obvious practical reasons, but also for the unique directional lamellar morphologies it usually provides. Camphene has also attracted a lot interest, since the solidification temperature is much higher and avoids the use of low temperatures, and sublimation of the camphene crystals can occur at room temperature and pressure. Other solvents or solvent mix have been considered, such as glycerol, dioxane or camphor-naphtalene. The relationships between the slurry concentration (in volume) and the total porosity for these solvents and many different materials are plotted in figure 1. Of course, a great dispersion of the values is observed, due variations of the processing conditions, slurries formulation, nature and characteristics of the powders, and densification conditions. We can nevertheless already notice that the porosity seems to be systematically much lower when camphene is used, compared to samples processed with water. Interestingly, the total porosity is even lower with the glycerol or dioxane water mix, although only a few reports can be found and data scattering -as observed for water and campheneis unknown.

A first parameter to take into account is the volume change of the solvent during the solidification: the solidification of water leads to a slight volume increase, while camphene slightly shrinks. This effect alone can nevertheless not account for the magnitude of differences observed experimentally (complete densification in some cases).

The local density in the green body will depend on the particles redistribution and concentration behaviour during the solidification step. Only a few parameters should account for the characteristics of particle packing in the ice-templated porous bodies.

Here we use the model developed by Shanti et al. [START_REF] Shanti | Particle Redistribution During Dendritic Solidification of Particle Suspensions[END_REF] , which have proved to be reliable so far. Particles are rejected by the growing crystals and concentrate in the intercrystals space, until the osmotic pressure exceeds the capillary pressure. At this point, the solidliquid interface invades the interparticles space. The concentration at which this bal-ance is lost is referred to as the breakthrough concentration. The breakthrough concentration therefore defines, after the sublimation step, the relative proportion of macropores, corresponding to pure ice crystals, and micropores, corresponding to the interparticle spaces in the concentrated particles regions between the pure ice crystals. If the breakthrough concentration is lower, the fraction of macropores is lower and the fraction is micropores is greater. This balance can be of great importance for densification, since micropores are easily removed during the sintering step. Let us first assume that the observed differences in densification behaviours originate from variations of the breakthrough concentration.

The breakthrough concentration φb can be expressed as

3 1 2 4 ÷ ÷ ø ö ç ç è ae - = g p f f R kT m b (eq. 1)
or

W m b - = f f (eq.2) with 3 1 2 4 ÷ ÷ ø ö ç ç è ae = g pR kT W (eq.3)
φm being the volume fraction of particles at maximum packing, T the temperature, R the particle radius (assuming particles are spherical) and γ the surface tension of the solvent. The only parameter of the solvent playing a role here is therefore the surface tension.

When plotting W vs. surface tension (fig. 2), water and camphene exhibits a very similar behaviour for usual particle sizes (0.2-1 μm). Variations of W become important for low surface tension and low particle size. This is more clearly observed when plotting W vs. particle size, for water and camphene (fig. 3). The critical value of particle size, at which W becomes non-negligible (>1%), is around 300 nm for camphene and only 50 nm for water. Although the typical average particle size in commonly used ceramics powders is around 0.2-1 μm, the presence of a population of small particles (<100 nm) is almost systematic. It is also known that such a population of small particles, even if very limited, can have a tremendous influence over the osmotic pressure of colloidal suspensions [START_REF] Rami | Colloidal Stability for Concentrated Zirconia Aqueous Suspensions[END_REF] . In our case, an increase of the osmotic pressure induces an important shift of the breakthrough concentration. We can therefore assume that in all cases, the breakthrough concentration will be lowered when the surface tension is lowered. The fraction of micropores in the green bodies processed with camphene should be greater than that of those processed with water. If we assume that this has a first order influence over the densification mechanism, since micropores are easily removed, a greater fraction of micropores should make densification easier, so that we should expect lower porosity content with camphene or other low surface tension solvents. The experimental results are in good agreement with this theoretical result:

porosity is systematically lower with camphene than with water. The differences in porosity between water-based and camphene-based materials are nevertheless too large to be accounted for only by the previously exposed reason. A variation of the surface tension, with the consequences previously exposed, should results in variation of the porosity of a few percent (typically 1 to 5%), no more, as shown in figure 2 and3. The experimental data show that the variations in porosity can be as high a 20% with camphene, and even up to 40% with camphor-naphtalene. The balance between micropores and macropores is therefore probably not the most important parameter controlling the densification, according to the model of Shanti et al.

The Shanti paper, in deriving their expression for φb, actually neglects the effect of the solvent on the osmotic pressure P of the colloidal particles. That solvent effects cannot be neglected is illustrated by measurements of the osmotic pressure of silica particles (R=0.25μm) dispersed in water (.01M NH4Cl) and in toluene [START_REF] Guo | Aggregation Effects on the Compressive Flow Properties and Drying Behavior of Colloidal Silica Suspensions[END_REF] . In toluene the osmotic pressure, at a given particle volume fraction, is found to be several orders of magnitude greater than in water. Thus changing the solvent can have a huge effect on the measured osmotic pressure. The analysis given in the Shanti paper cannot account for this effect as they use an expression for the osmotic pressure that is independent of the solvent, and is based on a hard-sphere equation of state. The solid-liquid surface energy of water varies from 0.07 to 0.03 N/m 5 , depending on its purity, while that for toluene is approximately 0.01 N/m [START_REF] Hatase | Discovery of homogeneous-nucleation-based crystallization in simple glass-forming liquid of toluene below its glass-transition temperature[END_REF] . Applying these values to the equation 1 yields (assuming φm = 0.64):

φb (water) = φm -.006 = 0.634 φb (toluene) = φm -.008 = 0.632 Thus, according to the Shanti model, it would seem changing the solvent has little effect on the breakthrough concentration. However, as noted above the Shanti derivation neglects the effect of the solvent on the osmotic pressure. Instead of the Shanti equation one should use the more general equation of the osmotic pressure derived by White [START_REF] White | Capillary rise in powders[END_REF] (assuming spherical particles) to determine φb:

( )R b b f g f - = P 1 3 (eq.4)
where P is the measured osmotic pressure. To determine φb from this equation one chooses φb until the left hand and right hand sides match. If we apply this equation to the experimental data of Guo [START_REF] Guo | Aggregation Effects on the Compressive Flow Properties and Drying Behavior of Colloidal Silica Suspensions[END_REF] we obtain: φb (water) = 0.63 φb (toluene) = 0.4

Thus changing the solvent would significantly reduce the breakthrough concentration. This seems much closer to the values required in the camphene case to account for the variations of the experimental observations. One way to test this more general relation would be to determine the osmotic pressure of the colloidal particles when dispersed in both water and camphene, and the other solvents such as dioxane and glycerol. If they differ by several orders of magnitude, it is likely the more general White equation will yield significantly different breakthrough concentrations.

Finally, we propose here that an additional parameter is controlling the densification of the green bodies, that is, the local curvature of the solid interface. It is well established that the diffusion processes, associated to grain growth and densification during solid-state sintering of ceramic powders, are governed by the curvature of the interfaces (solid / vapour for densification, solid / solid for grain growth). The driving force is the pressure difference across the interface, DP, having the well known expression:

R DP g 2 = (eq. 5)
with g the surface energy of interest and R the radius of curvature of the interface. The lower the radius of curvature of the interfaces, the higher the driving force for grain growth and/or densification.

In the case of ice-templated ceramics, very different morphologies are obtained when different solvent are used [START_REF] Deville | Freeze-Casting of Porous Ceramics: A Review of Current Achievements and Issues[END_REF] . In particular, in the case of camphene, highly dendritic structures are obtained, therefore with a low radius of curvature of the interface. When using water, lamellar morphologies of greater pores dimensions and with a less intricate structure are obtained. In the case of glycerol and dioxane/water suspensions, the structures are even more intricate, with smaller and more dendritic pores. The experimental results summarized in figure 1 are in good agreement with this hypothesis. Green bodies with structures exhibiting a low local radius of curvature will facilitate the densification and yield porous structures with a lower porosity content, which is schematically represented in figure 4. For an equivalent initial slurry concentration, the final total porosity decreases when the low local radius of curvature decreases. For green bodies where the pores are very small and intricate [START_REF] Araki | Room-Temperature Freeze Casting for Ceramics with Nonaqueous Sublimable Vehicles in the Naphthalene-Camphor Eutectic System[END_REF] , the porosity can be completely removed by densification.

Conclusions

Based on a review of the literature data and an analysis of the model developed by Shanti et al., we propose here that the densification of ice templated ceramics is governed both by the influence of the osmotic pressure during the solidification stage and by the morphology of the macropores, and more specifically the radius of curvature of the solid/vapour interface. The choice of a solvent for the process, such as water, camphene, dioxane or glycerol, will have a partial influence through the variation of the surface tension, which affects local particles packing and therefore the balance between macroporosity and microporosity, but this alone cannot explain the magnitude of porosity variations observed experimentally. The Shanti derivation neglects the effect of the solvent on the osmotic pressure. Instead of the Shanti derivation one should, for instance, use the more general equation of the osmotic pressure derived by White.

Variations of the osmotic pressure when changing the solvent can yield significantly different behaviour during the freezing stage and therefore significantly affect the densification behaviour, an effect not apparent when using the Shanti model. In addition, and equally important, a modification of the pores morphology is induced by the evolution of the solvent crystals morphology when changing the solvent. Highly dendritic and small crystals will facilitate the densification of the green body and yield porous materials with a lower total porosity. Determination of the osmotic pressure in the typical ice-templating suspensions are required to validate these propositions.

Fig 1 :

 1 Fig 1: Total porosity vs. slurry concentration, literature data 2,8-31 . The samples made with water, and camphene are respectively in white and dark. Other solvent are also represented, although fewer results can be found. The porosity of samples made with camphene is almost systematically lower than when made with water.

Fig 2 :

 2 Fig 2: W vs. surface tension, for various particle sizes, and for a fixed temperature of 273K. The surface tension of water is very sensitive to the presence of impurities and usually lower than the value for pure water.

Fig 3 :Fig. 4 :

 34 Fig 3: W vs. particle size, for camphene at 319K and water at 263 K

Table 1 :

 1 Surface tension of commonly used solvents[START_REF] Bayander | Solid-liquid interfacial energy of camphene[END_REF]

	Solvent	Surface tension (N/m)
	Water	0.07 to 0.03
	Glycerol	0.065
	Dioxane	0.03
	Camphene	0.0044

Corresponding author, E-mail : sylvain.deville@saint-gobain.com


Peter W Glynn

Marvin K Nakayama

Bruno Tuffin

Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets

We develop simulation estimators of measures associated with the tail distribution of the hitting time to a rarely visited set of states of a regenerative process. In various settings, the distribution of the hitting time divided by its expectation converges weakly to an exponential as the rare set becomes rarer. This motivates approximating the hitting-time distribution by an exponential whose mean is the expected hitting time. As the mean is unknown, we estimate it via simulation. We then obtain estimators of a quantile and conditional tail expectation of the hitting time by computing these values for the exponential approximation calibrated with the estimated mean. Similarly, the distribution of the sum of lengths of cycles before the one hitting the rare set is often well-approximated by an exponential, and we analogously exploit this to estimate tail measures of the hitting time. Numerical results demonstrate the effectiveness of our estimators.

INTRODUCTION

Many stochastic processes possess a regenerative structure, so the process "probabilistically restarts" at an increasing sequence of regeneration times; e.g., see [START_REF] Shedler | Regenerative Stochastic Simulation[END_REF] and [START_REF] Kalashnikov | Topics on Regenerative Processes[END_REF]. Suppose the process rarely visits some set A of states, and we are interested in estimating (performance or risk) measures associated with the tail distribution of the hitting time T to A . For example, in a stable GI/G/1 queue, the set A may correspond to a large number of customers in the system (e.g., buffer overflow), so A is rarely hit. In a highly reliable Markovian system consisting of a collection of components that fail and get repaired, system failures occur when certain combinations of components are down; in this case, the set A corresponds to the failed states, which are rarely visited.

Under various assumptions and asymptotic regimes, as visits to A become rarer, the distribution of the ratio of the hitting time T to A divided by its expectation µ converges weakly to an exponential; see Chapter 3 of [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF]. These results generalize Rényi's theorem (Proposition 1.1.2 of [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF]), which establishes that as p → 0, the product of p times the sum of a geometrically distributed number (with parameter p) of independent and identically distributed (i.i.d.) nonnegative random variables with finite mean converges weakly to an exponential. For our regenerative setting, the weak convergence motivates approximating the distribution F of T by an exponential with mean µ. As µ is unknown, we estimate it via simulation to calibrate the approximation. We then obtain estimators of the q-quantile (for a fixed 0 < q < 1) and the conditional tail expectation (CTE) of T by computing these values for the calibrated exponential approximation, where the CTE is the conditional expectation of T given that it exceeds its q-quantile.

We extend the idea by exploiting similar exponential approximations to the distribution G of the sum S of the lengths of the cycles before the one in which the rare set A is hit. The exponential approximation depends on the unknown mean η of S, and we use simulation to estimate η, which provides us with an estimator for the distribution G of S. We further simulate to estimate the distribution H of the time V to hit A in the cycle in which A is visited. We can then express the hitting time T to A as the sum of S and V . The regenerative property guarantees that S ∼ G and V ∼ H are independent, so the distribution F of T is the convolution of G and H. Taking the convolution of our simulation estimators of G and H thus leads to an estimator of F, and we then compute the q-quantile and CTE of the estimated F. We present numerical results showing the effectiveness of our methods. The rest of the paper unfolds as follows. Section 2 describes the problem mathematically and develops the notation. Section 3 explains the asymptotic regimes under which hitting the set A is a rare event, and discusses the weak convergence of T /µ and S/η to exponentials. Section 4 (resp., 5) exploits the resulting exponential approximation to T (resp., S) to develop our estimators of the q-quantile and the CTE of the hitting time. We give numerical results in Section 6, and concluding remarks appear in Section 7.

PROBLEM DESCRIPTION AND NOTATION

Consider a continuous-time stochastic process X = [X(t) : t ≥ 0] evolving on a state space S . For A ⊂ S a subset of states (e.g., "failed states"), define T = inf{t ≥ 0 : X(t) ∈ A } as the hitting time (or first passage time) to A . Let F be the cumulative distribution function (CDF) of T . For fixed 0 < q < 1, our goal is to estimate the q-quantile

ξ = F -1 (q) ≡ inf{t : F(t) ≥ q} (1)
of F and the conditional tail expectation (CTE)

γ = E[T | T > ξ ]. (2)
In the finance context (e.g., Section 2.2 of [START_REF] Mcneil | Quantitative Risk Management: Concepts, Techniques, Tools[END_REF] and [START_REF] Hong | Monte Carlo methods for value-at-risk and conditional value-at-risk: A review[END_REF]), a quantile is often called a value-at-risk (VaR), and the CTE is also known as the expected shortfall or the conditional value-at-risk (CVaR). We assume that X is (classically) regenerative, with 0 = Γ 0 < Γ 1 < Γ 2 < • • • as the sequence of regeneration times of X, so the process "probabilistically restarts" at each Γ i ; see p. 19 of [START_REF] Kalashnikov | Topics on Regenerative Processes[END_REF]. For example, an irreducible continuous-time Markov chain (CTMC) on a finite state space is regenerative, with successive hits to a fixed state forming a sequence of regeneration times. For i ≥ 1, let

τ i = Γ i -Γ i-1 , and the process [X(Γ i-1 + s) : 0 ≤ s < τ i ] is called the ith (regenerative) cycle of X, which has length τ i . As X is regenerative, (τ i , [X(Γ i-1 + s) : 0 ≤ s < τ i ]), i ≥ 1,
is a sequence of i.i.d. pairs of cycle lengths and cycles. Let τ be a generic copy of τ i . For i ≥ 1, let T i = inf{t ≥ 0 : X(Γ i-1 + t) ∈ A } be the time elapsing after Γ i-1 until the next hit to A . For x, y ∈ ℜ, let x ∧ y = min(x, y) and x ∨ y = max(x, y). Let I (•) be the indicator function, which takes value 1 (resp., 0) when its argument is true (resp., false). The regenerative property implies that (τ i , T i ∧ τ i , I (T i < τ i )), i ≥ 1, is an i.i.d. sequence of triplets. Let N(0) = 0, and for j ≥ 1, let N( j) = inf{i > N( j -1) : T i < τ i } be the index i of the cycle corresponding to the jth cycle in which A is hit. For j ≥ 1, let M( j) = N( j) -N( j -1) -1, which is the number of cycles that do not hit A between cycles N( j -1) and N( j). We can then express the hitting time to A as

T = M(1) ∑ i=1 τ i + T M(1)+1 , (3)
where M(1) may depend on τ 1 , τ 2 , . . . , τ M(1) and T M(1)+1 .

We next give a stochastically equivalent representation of T in (3) in terms of independent random variables. Let W be a random variable having CDF G W with

G W (x) = P(τ ≤ x | τ < T ). (4)
Also, let V be a random variable with CDF H, where

H(y) = P(T ≤ y | T < τ). (5)
Let M be a geometric random variable with P(M = k) = p(1p) k for each k ≥ 0, where

p = P(T < τ). (6)
Let W 1 ,W 2 , . . . be i.i.d. copies of W , which are independent of V and M, where V and M are also independent. Define S = ∑ M i=1 W i . Let G be the CDF of S, and let η = E[S]. The regenerative property of X ensures that

T D = S +V, with S ∼ G independent of V ∼ H, (7)
where D = denotes equality in distribution. Define µ = E[T ], which is the expected hitting time to the set A . As is well known (e.g., see [START_REF] Goyal | A unified framework for simulating Markovian models of highly dependable systems[END_REF] and [START_REF] Glynn | On the estimation of the mean time to failure by simulation[END_REF]), the regenerative structure of X allows us to express µ as a ratio

µ = E[T ∧ τ] p ≡ ζ p (8)
with p from (6), and both the numerator and denominator in [START_REF] Ecuyer | Approximating zero-variance importance sampling in a reliability setting[END_REF] are expectations of cycle-based quantities.

ASYMPTOTIC REGIMES

To develop estimators of the q-quantile ξ = F -1 (q) of the CDF F of T D = ∑ M i=1 W i +V and the CTE γ, we consider some approximations that require p in (6) to be small. For a theoretical framework to accommodate this, we parameterize the problem by introducing a rarity parameter ε > 0 and examine the behavior of F ≡ F ε as ε → 0, where we assume that

p ≡ p ε → 0 as ε → 0. (9)
We now provide examples that use such parameterizations. In the first example the rarity comes from a receding set A ≡ A ε of failed states, with step-wise probability distributions independent of the parameterization. In the second example, it is the opposite: the transitions of the discrete-event system depend on the parameterization, but the set A of failed states does not.

Example 1 For a stable GI/G/1 queue with first-in, first-out discipline, let X(t) denote the number of customers in the system at time t ≥ 0, where the first customer arrives at time t = 0 to an empty system. The process X with the state space S = {0, 1, 2, . . .} is regenerative with the beginnings of busy periods as regeneration times; e.g., see p. 16 of [START_REF] Kalashnikov | Topics on Regenerative Processes[END_REF]. We are interested in the distribution of the time when X first hits a high level b ε ≡ 1/ε , where the interarrival-and service-time distributions do not vary with ε. Thus, we let the set

A ≡ A ε = {b ε , b ε + 1, b ε + 2, . . .}, and T ε = inf{t ≥ 0 : X(t) ∈ A ε }
represents the first time that the queue length hits b ε -1. Theorem 1 of [START_REF] Sadowsky | Large deviations theory and efficient simulation of excessive backlogs in a gi/gi/m queue[END_REF] shows that (9) holds.

Example 2 Consider a highly reliable Markovian system (HRMS), as studied in [START_REF] Shahabuddin | Importance sampling for highly reliable Markovian systems[END_REF], [START_REF] Nicola | Techniques for fast simulation of models of highly dependable systems[END_REF], and [START_REF] Rubino | Markovian models for dependability analysis[END_REF]. The system comprises a finite set of components, each of which has exponentially distributed lifetimes and repair times. The components may be of different types, which may require different classes of repairmen. The system is failed when certain combinations of components are failed. At time t = 0, all components are operational. We can model the evolution of the system as a CTMC X = [X(t) : t ≥ 0] on a finite state space S , where each state in S specifies the number of failed components of each type, along with any necessary information about the queueing of failed components for each repairman class. Thus, X is regenerative, with hits to the state with all components operational as regeneration times. Assume that the system is highly reliable in the sense that the components' failure rates are much smaller than their repair rates. We then parameterize component failure rates as positive powers of the rarity parameter ε, and assume that repair rates and the set A of failed states do not change as ε → 0. Under appropriate assumptions (see [START_REF] Shahabuddin | Importance sampling for highly reliable Markovian systems[END_REF]), we have that (9) holds. In this setting, µ ε = E ε [T ε ] is often called the mean time to failure (MTTF).

Because we now actually are considering a family of models indexed by the rarity parameter ε > 0, we should thus write

T = T ε , W = W ε , M = M ε , V = V ε , p = p ε , F = F ε , G = G ε , G W = G W,ε , H = H ε , P = P ε , E = E ε , µ = µ ε , etc.
However, we often omit the subscript ε to simplify notation.

Under a variety of different sets of assumptions under which ( 9) is true, the scaled hitting time T ε /µ ε converges weakly to an exponential: for each x ≥ 0,

P ε (T ε /µ ε ≤ x) → 1 -e -x as ε → 0, ( 10
)
where we recall that

µ ε = E ε [T ε ].
To handle settings as in Example 1 of the GI/G/1 queue with receding sets A ε , we can apply Theorem 3.4.1 of [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF], which provides conditions to ensure [START_REF] Nakayama | Quick simulation techniques for estimating the unreliability of large regenerative models of highly reliable systems[END_REF] holds for a Harris-recurrent Markov chain, where the set A ε varies with ε so that (9) remains valid. For the HRMS in Example 2, the system dynamics (transition probabilities and holding-time distributions) change with ε, which requires a "triangular array" formulation; Theorem 3.2.5 of [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF] gives general conditions (not only for HRMSs) guaranteeing the validity of (10) when p ε , G W,ε , and H ε are allowed to depend on ε (see also [START_REF] Nakayama | Quick simulation techniques for estimating the unreliability of large regenerative models of highly reliable systems[END_REF] for the specific HRMS context).

It is often the case that when ( 10) is true, we further have that the sum

S ε = ∑ M ε i=1 W ε,i satisfies P ε (S ε /η ε ≤ y) → 1 -e -y as ε → 0 (11)
for each y ≥ 0, where we recall that

η ε = E ε [S ε ].
For example, if we assume that V = V ε ≡ 0 in [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF], the conditions in Theorem 3.2.5 of [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF] also ensure that [START_REF] Nicola | Techniques for fast simulation of models of highly dependable systems[END_REF] holds.

Throughout the rest of the paper, we apply non-simulation approximations based on the asymptotic results in ( 10) and [START_REF] Nicola | Techniques for fast simulation of models of highly dependable systems[END_REF], and then calibrate the approximations by estimating the unknown parameters µ ε and η ε via simulation. To distinguish the estimators and approximations for the different methods, we adopt the following notational convention. For an unknown quantity α, such as a CDF or parameter, we let α denote a non-simulation approximation to α. Also, we let α denote an estimator of α constructed from simulation-generated data.

APPROXIMATING THE CDF F OF T BY AN EXPONENTIAL

The limiting result [START_REF] Nakayama | Quick simulation techniques for estimating the unreliability of large regenerative models of highly reliable systems[END_REF] suggests that for small ε (which we now drop to simplify the notation),

F(t) = P(T ≤ t) = P(T /µ ≤ t/µ) ≈ 1 -e -t/µ ≡ F exp (t) (12)
for each t ≥ 0, where we recall that µ = E[T ]. We will next use the approximation F exp to F to obtain approximations to the q-quantile ξ = F -1 (q) in (1) for a fixed 0 < q < 1 and the CTE γ in (2). The exponential approximation in [START_REF] Rubino | Markovian models for dependability analysis[END_REF] motivates approximating ξ by

ξ exp = F -1 exp (q) = -µ ln(1 -q). ( 13
)
For the CTE, if T has exactly CDF F exp , then its CTE is

γ exp = ξ exp + µ = µ[1 -ln(1 -q)] (14)
by the memoryless property of F exp . But in ( 12), [START_REF] Sadowsky | Large deviations theory and efficient simulation of excessive backlogs in a gi/gi/m queue[END_REF], and ( 14), the parameter µ is unknown, so we next calibrate our approximations by estimating µ via simulation.

Simulation Estimator of µ

As we saw in [START_REF] Ecuyer | Approximating zero-variance importance sampling in a reliability setting[END_REF], the expected hitting time can be represented as a ratio µ = ζ /p. Because the asymptotic result ( 10) needs (9) to hold, the denominator p = P(T < τ) is a rare-event probability, so we will estimate it with importance sampling (IS); see Chapters V and VI of [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF] for an overview of this variance-reduction technique. But the numerator ζ = E[T ∧ τ] can be more efficiently handled by crude simulation (i.e., without IS), so we will independently estimate ζ and p. [START_REF] Goyal | Measure specific dynamic importance sampling for availability simulations[END_REF] call this approach measure-specific importance sampling, which we implement as follows.

To estimate the numerator ζ = E[T ∧ τ] in (8), we generate T i ∧ τ i , i = 1, 2, . . . , s, as s i.i.d. copies of T ∧ τ sampled using crude simulation. An estimator of ζ is then

ζ ≡ 1 s s ∑ i=1 T i ∧ τ i . (15)
Independently of the simulation runs employed to construct ζ in [START_REF] Shedler | Regenerative Stochastic Simulation[END_REF], we use IS to estimate the denominator p = P(T < τ) in [START_REF] Ecuyer | Approximating zero-variance importance sampling in a reliability setting[END_REF] as follows. Applying a change of measure, write

p = E[I (T < τ)] = E [I (T < τ)L], (16)
where E denotes expectation under the IS measure, and L is the corresponding likelihood ratio. The representation (16) motivates the following approach to estimate p. Let (I (T i < τ i ), T i ∧ τ i , L i ), i = 1, 2, . . . , r, be i.i.d. copies of (I (T < τ), T ∧ τ, L) generated via IS. We then estimate p by

p = 1 r r ∑ i=1 I (T i < τ i )L i . (17)
We combine the two estimators ζ from (15) and p from (17) to obtain

µ = ζ p (18)
as the simulation estimator of µ in (8).

Simulation Estimators of F, ξ , and γ

We next employ the simulation estimator µ from (18) to calibrate the approximate CDF F exp in (12) of T . Specifically, we replace µ in ( 12) by µ to obtain

F exp (t) = 1 -e -t/ µ (19)
as a simulation estimator of F(t) for each t ≥ 0.

Similarly, for the q-quantile ξ = F -1 (q) of F, we approximated it in ( 13) by ξ exp = -µ ln(1q). Replacing µ by its simulation estimator µ from (18) leads to the estimator

ξ exp = F -1 exp (q) = -µ ln(1 -q) (20)
of ξ . Finally, for the CTE, we replace ξ and µ in ( 14) by their simulation estimators ξ exp and µ to obtain

γ exp = ξ exp + µ = µ[1 -ln(1 -q)] (21)
as a simulation estimator of the CTE γ.

The simulation estimators F exp in (19) of F, ξ exp in (20) of ξ , and γ exp in (21) of γ are based on the approximation [START_REF] Rubino | Markovian models for dependability analysis[END_REF], which becomes more accurate as the rarity parameter ε → 0 in [START_REF] Nakayama | Quick simulation techniques for estimating the unreliability of large regenerative models of highly reliable systems[END_REF]. But for an actual physical system, we have a fixed value of ε > 0, so the exponential approximation F exp (t) in (12) typically will not exactly equal F(t). Thus, the simulation estimators F exp , ξ exp , and γ exp will often be biased.

APPROXIMATING THE CDF G OF S BY AN EXPONENTIAL

Recall that T D = S +V , where S ∼ G is independent of V ∼ H by [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF]. We next devise methods that separately estimate G and H to estimate the CDF F of T , its q-quantile ξ , and the CTE γ. To do this, we will approximate G by an exponential CDF, which is motivated by the asymptotic result in [START_REF] Nicola | Techniques for fast simulation of models of highly dependable systems[END_REF].

By [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF], we can write the CDF F of T as a convolution

F(t) = G H(t) = H(t -x) dG(x), (22)
where denotes the convolution operator. The convergence in [START_REF] Nicola | Techniques for fast simulation of models of highly dependable systems[END_REF] suggests that for small ε (which is now dropped to simplify the notation), we have the parametric approximation

G(x) = P(S ≤ x) = P(S/η ≤ x/η) ≈ 1 -e -x/η ≡ G exp (x). (23)
As η = E[S] is unknown, we will use simulation to estimate it, as will be discussed in Section 5.1, along with the estimation of the CDF H of V . Let ξ = F -1 (q) be the q-quantile of F. We next obtain another representation for the CTE γ, which only relies on the regenerative property of X and does not require the limiting result [START_REF] Nicola | Techniques for fast simulation of models of highly dependable systems[END_REF] nor the approximation in (23).

Theorem 1 If F is continuous at ξ , then the CTE satisfies

γ = 1 1 -q x [1 -H(ξ -x)] dG(x) + y [1 -G(ξ -y)] dH(y) . (24)
Proof. First write the CTE as

γ = E[T | T > ξ ] = E[T I (T > ξ )] P(T > ξ ) = E[T I (T > ξ )] 1 -q , ( 25
)
where the last step follows from the continuity of F at ξ . Express T D = S +V using [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF], so the numerator in the right side of (25) then satisfies

E[T I (T > ξ )] = E[S I (S +V > ξ )] + E[V I (S +V > ξ )] = E E[S I (V > ξ -S) | S] + E E[V I (S > ξ -V ) | V ] = E[S P(V > ξ -S | S)] + E[V P(S > ξ -V | V )] = E S [1 -H(ξ -S)] + E V [1 -G(ξ -V )]
because of the independence of S and V by [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF]. Thus, (24) immediately follows.

Simulation Estimators of η, G, and H

As the approximate CDF G exp in (23) and the CTE representation in (24) depend on the unknown η, G, H, and ξ , we next describe simulation estimators for them. We first explain how to handle η

= E[S] = E[∑ M i=1 W i ]. Writing S = ∑ ∞ i=1 W i I (M ≥ i), we see that η = ∞ ∑ i=1 E[W i I (M ≥ i)] = E[W ] ∞ ∑ i=1 E[I (M ≥ i)] = E[W ] ∞ ∑ i=1 P(M ≥ i) = 1 -p p ν (26)
by the independence of M and W i , i ≥ 1, which are i.i.d. with mean ν = E[W ]. Thus, (26) shows that η can be expressed as a function of expectations of cycle-based quantities, where we recall that W ∼ G W , with G W the conditional CDF of τ given that τ < T , as defined in (4). Hence,

ν = E[τ | τ < T ] = E[τI (τ < T )]/(1 -p), and we estimate E[τI (τ < T )] by (1/s) ∑ s i=1 τ i I (τ i < T i )
, where (τ i , I (τ i < T i )), i = 1, 2, . . . , s, are from the same i.i.d. crude-simulation data used to construct ζ in [START_REF] Shedler | Regenerative Stochastic Simulation[END_REF]. Then we estimate ν via

ν = 1 (1 -p)s s ∑ i=1 τ i I (τ i < T i ), ( 27
)
where p is from (17). Therefore, we arrive at our estimator of η based on (26) as

η = 1 -p p ν. ( 28
)
We then replace η in (23) by its estimator η from (28) to obtain

G exp (x) = 1 -e -x/ η (29)
as a parametric estimator of G exp (x) in (23).

We next discuss how to use importance sampling to estimate the CDF H of V , where we recall that H is the conditional CDF of T given that T < τ, as in [START_REF] Hong | Monte Carlo methods for value-at-risk and conditional value-at-risk: A review[END_REF]. Note that

H(y) = P(V ≤ y) = P(min(T, τ) ≤ y | T < τ) = P(min(T, τ) ≤ y, T < τ) p = E[I (min(T, τ) ≤ y, T < τ)] p = E [I (min(T, τ) ≤ y, T < τ) L] p , (30)
where, as before, E is the expectation operator under importance sampling, as in ( 16), and L is the corresponding likelihood ratio. We previously obtained the IS estimator p in (17) to handle the denominator p of (30). For the numerator, let (I (T i < τ i ), T i ∧ τ i , L i ), i = 1, 2, . . . , r, be the same i.i.d. copies of (I (T < τ), T ∧ τ, L) obtained through IS that we also employed to construct the estimator p. We then build a nonparametric estimator H of H as

H(y) = 1 p • r r ∑ i=1 I (T i ∧ τ i ≤ y, T i < τ i ) L i . ( 31
)
While other nonparametric estimators of H can also be constructed, e.g., by interpolating H in (31), or by using kernel methods, etc., we will only work with H in (31).

Simulation Estimators of F, ξ , and γ

Now that ( 29) and (31) provide simulation estimators for the CDFs G and H of S and V , respectively, we use them to build an estimator for the CDF F of T . The representation of F in ( 22) as a convolution of G and H suggests estimating F(t) by

F (t) ≡ G exp H(t) = H(t -x) d G exp (x). ( 32
)
The following result works out an expression for F (t). Proposition 1 The estimator F in (32) of the CDF F of T satisfies

F (t) = 1 - 1 p • r r ∑ i=1 I (T i < τ i ) L i e -(t-A i ) + / η , (33)
where A i = T i ∧ τ i , η is defined in (28), and x + = max(x, 0).

Proof. Put (29) and ( 31) into (32) and use the fact that G exp (x) = 0 for x < 0 to get

F (t) = 1 p • r r ∑ i=1 I (T i < τ i ) L i t 0 I (A i ≤ t -x) d G exp (x) = 1 p • r r ∑ i=1 I (T i < τ i ) L i (t-A i ) + 0 d G exp (x) = 1 p • r r ∑ i=1 I (T i < τ i ) L i G exp ((t -A i ) + ),
which equals (33) by ( 29) and ( 17).

The corresponding estimator of the q-quantile ξ

= F -1 (q) is ξ = F -1 (q). ( 34
)
Computing ξ may require applying a numerical root-finding method, such as the bisection method or Newton's method, which can be computationally costly.

We next give a simulation estimator for the CTE γ.

Proposition 2 When ξ , G, and H in (24) are replaced by their respective estimators ξ in (34), G exp in (29), and H in (31), the resulting CTE estimator is

γ = 1 1 -q 1 p • r r ∑ i=1 L i I (T i < τ i ) ( ξ ∨ A i ) + η e -( ξ -A i ) + / η , (35)
where

A i = T i ∧ τ i .
Proof. The first term in the outer square brackets in (24) is expressed in terms of 1 -H(•), and note that

1 -H(y) = p p - 1 p • r r ∑ i=1 I (A i ≤ y) I (T i < τ i ) L i = 1 p • r r ∑ i=1 [1 -I (A i ≤ y)] I (T i < τ i ) L i = 1 p • r r ∑ i=1 I (A i > y, T i < τ i )L i . (36)
where the third equality holds by (17). Thus, replacing ξ , G, and 1 -H in the first term inside the outer square brackets of (24) by their respective estimators leads to

x [1 -H( ξ -x)] d G exp (x) = 1 p • r r ∑ i=1 L i I (T i < τ i ) ∞ x=0 x I (A i > ξ -x) d G exp (x) = 1 p • r r ∑ i=1 L i I (T i < τ i ) ∞ x=( ξ -A i ) + x d G exp (x) = 1 p • r r ∑ i=1 L i I (T i < τ i ) ( ξ -A i ) + + η e -( ξ -A i ) + / η , (37)
where we recall η is defined in (28).

The second term inside the outer square brackets of (24) becomes

y [1 -G exp ( ξ -y)] d H(y) = 1 p • r r ∑ i=1 L i I (T i < τ i ) ∞ y=0 y [1 -G exp ( ξ -y)] I (A i ∈ dy) = 1 p • r r ∑ i=1 L i I (T i < τ i ) A i [1 -G exp (( ξ -A i ) + )] = 1 p • r r ∑ i=1 L i I (T i < τ i ) A i e -( ξ -A i ) + / η . (38)
By replacing the two terms inside the outer square brackets of (24) by their estimators (37) and (38), we obtain the estimator of (24) as

γ = 1 1 -q 1 p • r r ∑ i=1 L i I (T i < τ i ) ( ξ -A i ) + + η + A i e -( ξ -A i ) + / η , which equals (35) because ( ξ -A i ) + + A i = ( ξ -A i + A i ) ∨ (0 + A i ) = ξ ∨ A i .

NUMERICAL RESULTS

We now present numerical results for our estimators, from Sections 4.2 and 5.2, of the CDF F of the hitting time T , its qquantile ξ , and the CTE γ. Due to space limitations, we focus on only one simple model of an HRMS, as in Example 2 of Section 3. Specifically, the HRMS has three component types, with five components of each type. There are 15 repairmen, so failed components never queue for repair. The system is up whenever at least two components of each type work, so the failure set A comprises states having at least four components down of one type. Each component has failure rate ε and repair rate 1. We consider two versions of the model: one with ε = 10 -2 , and the other has ε = 10 -4 .

To implement our methods, we need to specify the IS distribution employed to construct the estimators p in (17) and H in (31), which subsequently are used in the estimators of F, ξ , and γ. [START_REF] Rubino | Markovian models for dependability analysis[END_REF] provide an overview of IS techniques designed to simulate HRMSs, where the basic idea is to increase the probability of failure transitions. In our experiments, we applied the IS approach known as Zero-Variance Approximation (ZVA) of [START_REF] Ecuyer | Approximating zero-variance importance sampling in a reliability setting[END_REF]. ZVA produces estimators of the MTTF µ in [START_REF] Ecuyer | Approximating zero-variance importance sampling in a reliability setting[END_REF] 32) for ε = 10 -2 (left) and ε = 10 -4 (right).

For the methods in Section 4, computing F exp , ξ exp , and γ exp in ( 19)-( 21) requires simulation to estimate only µ = ζ /p. To do this, we simulated for each ε a total of 10 4 cycles, of which we allocated s to construct the estimator ζ in (15) using crude simulation, and we sampled the remaining r cycles with IS to build the estimator p in (17). [START_REF] Goyal | A unified framework for simulating Markovian models of highly dependable systems[END_REF] derive the optimal allocation of s and r for a fixed total budget to minimize the work-normalized asymptotic variance of the MTTF ratio estimator µ in (18). We ran pilot simulations with 10 3 cycles without IS (resp., with IS) for the numerator (resp., denominator) to estimate the parameters determining the optimal allocation. Adding constraints that at least 10% of the 10 4 second-stage cycles are for simulating the numerator and the same for the denominator, we got the optimal r = 9000 and s = 1000 for both values of ε in our example. To build F in (33) of Section 5, we used the same s + r = 10 4 simulated second-stage cycles to compute the estimators ν in (27) and H in (31).

In the discussion below, we refer to F exp (resp., F ) and its corresponding estimators of ξ and γ in (20) and (21) (resp., (34) and ( 35)) as exponential (resp., convolution) estimators. For comparison, we also built an empirical distribution F from n = 10 4 i.i.d. observations T1 , T2 , . . . , Tn , of T for ε = 10 -2 generated via crude simulation, where F(t) = (1/n) ∑ n i=1 I ( Ti ≤ t). (We could not obtain an empirical distribution for ε = 10 -4 , as we will explain later.) To define the corresponding estimators of ξ and γ, let T1:n ≤ T2:n ≤ • • • ≤ Tn:n be the sorted Ti values. Then the empirical estimators of the q-quantile and CTE are ξ = F-1 (q) = T nq :n and γ = [1/((1q)n)] ∑ n i= nq Ti:n , respectively, where • is the ceiling function.

The left side of Figure 1 plots F(t), F exp (t), and F (t) for ε = 10 -2 , which shows that the three curves closely align. But the CPU times are not of the same order of magnitude, as seen in Table 1 for the estimation of quantiles (for each particular method, computing the CDF and quantile estimators requires roughly the same time). The right side of Figure 1 plots the exponential and convolution estimators of F(t) for ε = 10 -4 , but we were not able to obtain the empirical distribution. Indeed, it would have required a CPU time of more than one year to sample 10 4 observations of T for ε -4 using crude simulation (on average, a single run of T takes more than one hour). But as ε shrinks, the approximations in [START_REF] Rubino | Markovian models for dependability analysis[END_REF] and (23) become more accurate by virtue of the limiting results in [START_REF] Nakayama | Quick simulation techniques for estimating the unreliability of large regenerative models of highly reliable systems[END_REF] and [START_REF] Nicola | Techniques for fast simulation of models of highly dependable systems[END_REF], so the CDF estimators F exp and F should be close to the true F for ε = 10 -4 . Moreover, Table 1 shows that for this model, the CPU times for the exponential and convolution estimators do not increase when ε decreases.

Table 1 contains results for the empirical, exponential, the convolution estimators of the q-quantile, for three values of q. We used the bisection method to numerically compute the inverse in (34) for the convolution estimator. For ε = 0.01, the exponential and convolution estimators of ξ are close to the empirical estimator, but with much less computational effort expended. For the exponential estimator, we include a biased 95% confidence interval (CI) based on the CI for µ; see (20). (For ε fixed, the approximation in [START_REF] Rubino | Markovian models for dependability analysis[END_REF] leads to the estimators F exp and ξ exp having bias, which do not go away as the sample sizes grow large. In contrast, the bias of the ratio estimator (18) of µ vanishes as sample sizes increase.) For ε = 10 -4 , we are not able to provide an empirical estimator of ξ , but the estimator (18) of µ is so good that the (biased) CI of ξ based on the exponential estimator has a relative width of only about 0.1%. However, the exponential and convolution quantile estimators differ by about 0.3%, so the (biased) exponential CIs do not include ξ . This may indicate that ξ and ξ exp have different levels of bias.

Table 2 provides results for the CTE estimators. These numbers exhibit similar behavior to what we saw with the quantile estimators: much smaller CPU times for the exponential and convolution estimators than for the empirical, and very narrow (biased) exponential CIs (see ( 21)), with the convolution estimators for ε = 10 -4 just above the exponential CIs.

CONCLUDING REMARKS

We used simulation to calibrate approximations to the distribution F of the hitting time T to a rarely visited set A of states for a regenerative process. Section 4 (resp., 5) approximated the CDF F of T (resp., CDF G of S) by an exponential (12) (resp., (23)) based on the asymptotic result (10) (resp., [START_REF] Nicola | Techniques for fast simulation of models of highly dependable systems[END_REF]), which requires the rarity parameter ε → 0. But for an actual physical system, we have a fixed ε > 0, which introduces bias in both exponential approximations. Chapter 3 of [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF] provides upper and lower bounds to the true CDFs F and G, and we are investigating employing the bounds to obtain upper and lower bounds for the q-quantile ξ and the CTE of the true CDF F.

  having the (desirable) Bounded Relative Error (BRE) property, with the (even better) Vanishing Relative Error (VRE) holding under certain conditions. Plots of the empirical estimator F(t) of F(t), the exponential estimator F exp (t) from (19), and the convolution estimator F (t) from (

		1					1
		0.8					0.8
	Estimated F(t)	0.4 0.6				Estimated F(t)	0.4 0.6
		0.2			Empirical		0.2
					Exponential				Exponential
		0			Convolution	0			Convolution
		0	0.5	1	1.5	2	0	2	4	6
				t	•10 7		t		•10 14
	Figure 1:

Table 1 :

 1 Quantile estimators.

	ε	q	Empirical 95% CI	CPU	Expon. Est.	Expon. 95% CI	CPU	Convol. Est.	CPU
	0.01	0.1	(1.701e+05, 1.971e+05)	890 sec	1.830e+05	(1.764e+05, 1.896e+05)	0.3 sec	1.865e+05	0.4 sec
	0.01	0.5	(1.206e+06, 1.271e+06)	890 sec	1.204e+06	(1.161e+06, 1.247e+06)	0.3 sec	1.227e+06	0.4 sec
	0.01	0.9	(3.958e+06, 4.135e+06)	890 sec	4.000e+06	(3.856e+06, 4.143e+06)	0.3 sec	4.075e+06	0.4 sec
	10 -4	0.1	N/A	N/A	1.757e+13	(1.756e+13, 1.758e+13)	0.3 sec	1.762e+13	0.4 sec
	10 -4	0.5	N/A	N/A	1.155e+14	(1.154e+14, 1.157e+14)	0.3 sec	1.159e+14	0.4 sec
	10 -4	0.9	N/A	N/A	3.840e+14	(3.838e+14, 3.842e+14)	0.3 sec	3.850e+14	0.4 sec

Table 2 :

 2 CTE estimators

	ε	q	Empir. Est.	CPU	Expon. Est.	Expon. 95% CI	CPU	Convol. Est.	CPU
	0.01	0.1	1.964e+06	890 sec	1.920e+06	(1.851e+06, 1.989e+06) 0.3 sec	1.956e+06	0.4 sec
	0.01	0.5	3.011e+06	890 sec	2.941e+06	(2.836e+06, 3.046e+06) 0.3 sec	2.996e+06	0.4 sec
	0.01	0.9	5.915e+06	890 sec	5.737e+06	(5.531e+06, 5.942e+06) 0.3 sec	5.844e+06	0.4 sec
	10 -4	0.1	N/A	N/A	1.839e+14	(1.834e+14, 1.845e+14) 0.3 sec	1.848e+14	0.4 sec
	10 -4	0.5	N/A	N/A	2.817e+14	(2.809e+14, 2.826e+14) 0.3 sec	2.831e+14	0.4 sec
	10 -4	0.9	N/A	N/A	5.495e+14	(5.479e+14, 5.512e+14) 0.3 sec	5.523e+14	0.4 sec

ACKNOWLEDGMENTS

This work has been supported in part by the National Science Foundation under Grant No. CMMI-1537322. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.


H Touati

H Castel-Taleb

B Jouaber

S Akbarzadeh

A Khlass

Modeling and performance evaluation of the eICIC/ABS in H-CRAN

In this paper, we propose mathematical models to evaluate the performance of the interference remediation technique eICIC/ABS (enhanced Inter-Cell Interference Coordination / Almost Blank Sub-frame) in the context of Heterogeneous Cloud based Radio Access Networks (H-CRAN) architecture and 5G networks. The objective is to propose a dynamic resource management tool to ease decisions on the activation/deactivation of micro-cells as well as on the distributions of subframes among macro and micro cells. First, we propose a Markov chain based model that fits the behavior of the considered scheme and allows the analysis of the cell throughput according to traffic load, radio conditions and the distribution of available resources among macro and micro cells. Then, we propose an approximation model with a closed form formula. The two models are validated and evaluated in terms of accuracy and computation time. Numerical results are compared to matlab simulations that reproduce realistic radio conditions. Results show that both models are accurate. However, the closed form approximation is less complex and provides faster results.

Introduction

The H-CRAN architecture proposed for 5G cellular networks combines the Heterogeneous Networks concept (HetNets) with the Cloud Radio Access Network (C-RAN) architecture. The HetNet concept allows the use of small cells within macro cells in order to enhance edge user's radio conditions. In return, the small cells use the same radio resources and may increase interference levels. For network operators, deciding on the planning, activation and tuning of small cells is not an easy task. Many dynamic parameters should be considered including user's radio conditions, QoS and spectrum and power usage. The C-RAN architecture is also referring to Centralized RAN. It can help in achieving HetNet objectives since decisions can be centralized and based on global knowledge. In order to reduce and/or control the interferences among macro and micro Cells, eICIC and CoMP (Coordinated Multi-Point) mechanisms have been proposed by 3GPP for LTE-Advanced. In the eICIC/ABS (eICIC/Almost Blank Sub-frame), part of the radio resources (subframes) are delegated to the micro cells. During its allocated subframes, the macro cell keeps silent and the micro cell can transmit at reduced interference. With this approach, one of the most important challenges is how to determine the amount of sub-frames to delegate to micro cells. In this paper, we are interested in proposing a dimensioning tool for the eICIC/ABS mechanism. Our objective is to provide an easy to use mathematical model that can help, in real time, network managers for their decision making regarding dynamic activation and dimensioning of the small cells.

System model and assumptions

In this paper, we consider a H-CRAN cluster composed of a central cell and its 6 first tier neighboring cells, numbered from 2 to 7. All these cells are macro cells. The cluster also includes a micro cell deployed within the central macro cell. In the following, we will focus our analysis on users within the central region, covered by both the central macro and the micro cell. We focus our studies on downlink. Users are randomly distributed over the central region. Their serving cells are determined based on the higher received power. The amount of radio resources allocated to each user is determined based on (i) the data size the user is willing to receive, (ii) traffic load within the cell and (iii) the optimal MCS (Modeling and Coding Scheme). This latter is determined according to user's radio conditions subject to interference, attenuation and noise. In practice, the optimal MCS is selected based on the SINR (Signal to Interference and Noise Ratio). The SINR of user u is obtained through the following equation:

SINR(u)[db] = 10log 10 ( PRx(u)[W ] ∑ 7 i=2 PRx i (u)[W ] + No[W ] ) (1)
where PRx(u)[W] (resp. PRx i (u)[W ], for i = 2, . . . , 7) are the received power (in watt) from the serving cell (resp. the interfering cell i ), and No[W ] is the thermal noise power expressed in watt over the entire considered bandwidth. According to [Neta], the thermal noise per Hz is equal to -174. Note that the received power given in equation ( 1) depends on the pathloss PL(u)[db] between the user u and the BS (Base Station) in db. For the macro cell, the pathloss PL(u)[db] is given by:PL(u

)[db] = 128.1 + 37.6 × log 10 (R[km]
) and for the micro cell:

PL(u)[db] = 38 + 30 × log 10 (R[m]) where R[km] (resp. R[m]
) is the distance between the user and the BS in kilometers (resp. meters) [Neta]. Once the SINR of each user is determined, the optimal MCS is selected according to the 3GPP correspondence tables [Netb].

In LTE, radio resources are assigned to users at a frame basis. Each frame is composed of 10 sub-frames.

The BS allocates to each UE a set of GRBs (Group of Resource Blocks). The number of RBs (Resource Blocks) within a GRB depends on the considered bandwidth [Netb]. In the following, we denote by S the total number of available GRBs in a frame and by Nb RB the number of RBs within a GRB. For each user u, with a size (size(u)), and a MCS (MCS(u)), the number of required GRBs/frame is obtained as follows:

Nb GRBs (u) = Nb Subch (u)
Nb RB

where Nb Symb (u) = size(u) MCS(u) determine the number of required OFDM symbols and Nb Subch (u) =

Nb Symb (u) 7×12
determine the number of corresponding sub-channels (i.e. RBs).

Markov chain

We suppose that in the macro-cell (resp. micro-cell), clients arrive according to a Poisson process of rate λ M , (resp. λ m ). We suppose that arrivals are divided into classes corresponding to the batch size (to simplify, we consider that clients of class i have a bach size of i). The batch size of a client corresponds to the number of GRBs required to its transmission within a frame. Let λ M,i and J M (resp. λ m,i , and J m ) arrival rates of class i and the number of classes in the macro cell (resp. micro cell). We have

∑ J M i=1 λ M,i = λ M (resp. ∑ J m
i=1 λ i = λ m ) the total arrival rate in the macro cell and the micro cell respectively. In order to remedy the interference, S servers (which refers to S GRBs/frame) are shared between the macrocell and the micro-cell according to the number of ABS allocated to the micro-cell. During the frame duration, we consider the arrivals of the macro and the micro, and the services of at most S M = S -S × (nb ABS /10) where nb ABS represents the number of ABS assigned to the micro-cell, and S m = S × (nb ABS /10) GRBs. We represent the behavior of eICIC technique as two discrete-time Markov chains one for the macro and the other for the micro cell. We denote by Y M (n) (resp. Y m (n)) the number of GRBs occupying the macro (resp. the micro) at frame n. At each end of the frame time (n), we consider service of clients arrived Formatting a submission for AlgoTel at the precedent frame time n -1, and we add clients arrived in the present frame time n. Let A M,i (n) (resp. A m,i (n)) the number of arrivals at frame n in the macro cell (resp. micro cell). Then ∑

J M i=1 i × A M,i (n) (resp. ∑ J m i=1 i × A m,i (n)
represents the total GRBs necessary to arrivals in the macro cell (resp micro cell). So Y M (n) and Y m (n) for n ≥ 0 represent discrete-time Markov chains, and the evolution equations are given as follows:

Y M (n) = min(S M , [Y M (n -1) -S M ] + + J M ∑ i=1 i × A M,i (n))
(2)

Y m (n) = min(S m , [Y m (n -1) -S m ] + + J m ∑ i=1 i × A m,i (n)) (3)
where

Y M (0) = min(S M , ∑ J M i=1 i × A M,i (0)), and Y m (0) = min(S m , ∑ J m i=1 i × A m,i (0)).
The proposed Markov chain is so slow for a real-time dimensioning tool, due to the fact that we consider batch arrivals that generate very large number of transitions, multiserver (500 for 20MHz, 1000 for 40 Mhz, ...), and it takes 2 hours to simulate a point of the curve (when we consider 20 MHz) to ensure a 95% confidence interval. That's why we had to look for an approximation using a queueing model with a closed form.

Closed form equations

We propose to model the system with a multi-server loss queueing system. As the arrivals are supposed Poisson with different batch sizes (a client has a batch size of j GRBs), then we consider the loss discrete-time M/G/S/S queue [D.82] with batch arrivals and grouped departures to represent the behavior of the macro-cell. In this queueing model, each arrival corresponds to a batch of size j, which occupies simultaneously j servers if there is enough place [D.82], otherwise it is lost. Since in our model, the admission is made by element of the batch (or GRB) , and not by considering the whole batch (or the client batches), then the queueing model will increase the blocking probability. We propose to consider the M/G/S/S queueing model with Poisson simple arrivals, and deterministic service times. Each client occupy one server, but the service time is approximated by the mean service time of all batches sizes of the serving cell (macro or micro). So in this model, clients are not differentied by classes (or different service time) they have the same service time equal to the mean batch size of the clients. So in the exact model, if a client occupy j resources (or GRBs) in a frame time, then in the queueing system the client occupy one resource for j frame times. From M/G/S/S queue, we can use the Erlang loss formula [D.82] in order to derive the blocking probability of clients in the serving cell. From the blocking probabilities we can derive the throughput for each cell, and also the total throughput. Let P M (resp. P m ) the blocking probability of the macro cell (resp. the micro cell). They are computed as follows from Erlang loss formula [D.82]:

P M = ρ S M M S M ! ∑ S M n=0 ρ M n n!
where ρ M = λ M µ M , and

1 µ M = ∑ J M i=1 i × λ i,M
λ M and

P m = ρ Sm m Sm! ∑ Sm n=0 ρm n n!
, where

ρ m = λ m µ m , 1 µ m = ∑ J m i=1 i × λ i,m
λ m . So we can derive the throughput X M for the macro-cell as X M = λ M × (1 -P M ), the throughput X m = λ m × (1 -P m ) for the micro cell, and the total overall throughput is X = X M + X m .

Numerical results

We consider the following parameters: f=2 Ghz, BP = 20 Mhz which corresponds to 500 GRBs/frame and 1GRB= 4 subchannels [Netb], for the macro cell (resp. micro cell) the radius is R macro =3500m (resp. R micro =600m), maximum transmit power is 40dBm (resp. 30dBm), antenna gain= 18 dBi (resp. 6dBi), and client size= 1024 bits. Using Matlab simulator, we generate uniformly in the macro cell the clients whose number follows Poisson distribution. First, the simulator associate each user to his serving cell based on the highest power reception. Then, it calculates for each client the number of GRBs required as explained in the previous section. So clients with the same number of GRBs are regrouped into classes, the corresponding arrival rates and service times will be injected into the proposed Markov chain simulator and into the M/G/S/S formulas so that these two models are based on real scenarios. The figure below compares the performance results (the macro throughput, the micro throughput and the overall throughput) obtained by the simulator of real behavior, the proposed Markov chain and the M/G/S/S. As shown in these figures, in the considered scenarios (macro without/with micro) the throughputs obtained

Conclusion

In this paper, two models are proposed and validated through simulations for the eICIC/ABS mechanisms in the context of H-CRAN. The first is a Markov chain based model. It reflects the behavior of the system with batch arrivals corresponding to users with variable radio conditions and requiring more or less resources to transmit their data. The second model is an approximation using the M/G/S/S queue that provides a closed form result. Both models provide similar results in terms of throughput and system capacity for both the micro and the macro cells. Furthermore, both models results are evaluated by comparing with simulations that reproduce a real behavior of radio conditions. The proposed models can then be used as a dimensioning tool for the considered mechanism. The M/G/S/S approximation can in addition be used in real time to adapt in real time the eICIC/ABS mechanism according to user distribution, to load and to radio conditions.

  Fig. 1: Validation of the Markov chain and the M/G/S/S approximation with the Markov chain and the M/G/S/S formulas match perfectly with the simulation results of the real behavior. Thus, we can validate the proposed Markov chain and the simple M/G/S/S.