Eric Allain

N Kanari

I Gaballah

Kinetics of Oxychlorination of Magnesium Oxide

HAL is

I. INTRODUCTION

Magnesium is a major element of the chromite ores. It is present in the spinel structure of chromite having the general formula (Fe 2+ , Mg)(Cr, Al, Fe 3+ ) 2 O 4 and also in the gangue of chromite concentrates that is essentially composed of magnesium silicates such as olivine and serpentine. [1] Recently, chlorination was studied as a potential method to improve the quality of poor chromite concentrates and for the recovery of specific elements (Cr, Fe, ...) from chromite concentrates [1][START_REF] Gaballah | Final Report of European Union Contract No[END_REF][START_REF] Kanari | Proc. of the 126 th TMS Annual Meeting 'EPD Congress 1997[END_REF] . The oxychlorination of chromite concentrates by Cl 2 +O 2 at about 900 °C allowed almost full extraction of iron and partial recovery of chromium as chromium oxychloride (CrO 2 Cl 2 ). It was found [START_REF] Kanari | Light Metals[END_REF] that the kinetics parameters of the chromite oxychlorination such as activation energies and reaction orders changed with the reaction extents.

Magnesium metal is produced mainly either by the electrolysis of purified MgCl 2 from natural brine or by silicothermy. The use of light weight alloys of magnesium by the transport industry allows the production of lighter weight vehicles and consequently reduces the fuel consumption and CO 2 emissions [START_REF] Froes | The Science, Technology, and Applications of Magnesium[END_REF] . Moreover, magnesium oxide is employed in some applications, such as rubber processing, pharmaceuticals, and electrical insulator uses, which require a highpurity MgO [START_REF] Leroy | Production of high purity magnesium oxide by selective chlorination[END_REF] . The common impurities contained in MgO are Fe, Si, Ca ... oxides. For these reasons, kinetic studies of the reactions of MgO, Cr 2 O 3 , Fe 2 O 3 with Cl 2 +O 2 were undertaken. This paper focused on the oxychlorination kinetics of magnesium oxide in isothermal conditions using thermogravimetric analysis 'TGA'.

No studies concerning the oxychlorination of magnesium oxide were found in the literature. Most of investigations were devoted to the thermodynamical and kinetics aspects of the reactions of MgO with Cl 2 , Cl 2 +C, HCl, and HCl+C and they were summarized in Reference [7].

II. MATERIALS AND EXPERIMENTAL PROCEDURE

The sample used in this study was MgO powder having a purity of 98 pct and supplied by PROLABO (Paris, France). Results of Scanning Electron Microscopy (SEM) and X-ray diffraction (XRD) analyses of the magnesium oxide sample are showed in Figure 1. The main identified phase was MgO.

The TGA tests for the oxychlorination of MgO were performed using 40 mg of sample in a vertical experimental set-up schematized previously [START_REF] Kanari | Proc. of the 127 th TMS Annual Meeting 'EPD Congress 1998[END_REF] . Its main unit was a CAHN 1000 microbalance having a sensitivity of 10 µg. The oxychlorination gas mixture composed of Cl 2 +O 2 was purified through H 2 SO 4 and CaCl 2 columns before its introduction in the reactor. The sample was homogeneously distributed in the quartz crucible having a section of about 0.5 cm 2 . This crucible was hooked to the balance using quartz rods. The sample was preheated in nitrogen atmosphere up to the desired temperature before the introduction of the chlorinating gas mixture.

Exhaust gases were purified from chlorine excess by NaOH solution before their release to the atmosphere.

III. RESULTS

A. Effect of Cl 2 /(Cl 2 +O 2 ) ratio on the oxychlorination of MgO.

The chlorination of MgO by Cl 2 can be represented by Eq. [1]. The melting and boiling points of produced magnesium chloride are 714 °C and 1412 °C, respectively [START_REF] Anonymous | Handbook of Chemistry and Physics[END_REF] . As the TGA technique was used to measure the reaction rate, it was mandatory to insure that the rate of MgCl 2 volatilization was higher than that of its formation. One possible way to decrease the reaction rate of MgO with chlorine is to increase the partial pressure of oxygen in the system. Thus, a series of experiments was performed at 950 °C, using a Cl 2 +O 2 gas mixture having a molar ratio of Cl 2 /(Cl 2 +O 2 ) varying from 0.33 to 1.00. Figure 2 shows the evolution of the percentage weight loss (pct WL) of magnesium oxide versus time during its oxychlorination using the conditions mentioned above. Only weight losses were observed when using gas mixtures containing less than 90 pct of Cl 2 . A weight gain followed by a weight loss was recorded when the chlorine content in the gas mixture was higher than 90 pct. This last result indicates that the volatilization rate of magnesium chloride was lower than that of MgO chlorination and consequently the overall reaction rate was probably controlled by the MgCl 2 volatilization rate. For this reason, the kinetics parameters of MgO oxychlorination were studied using the Cl 2 +O 2 gas mixtures having a chlorine content lower or equal to 80 pct. In these conditions the observed weight loss was a direct measure of the extent of the overall oxychlorination reaction.

MgO

+ Cl 2

→ MgCl 2 + 1/2 O 2 [1]

B. Effect of gas flow rate

In order to determine the intrinsic parameters of the gas-solid reactions, it was necessary to minimize the effect of external mass transfer phenomena [START_REF] Szekely | Gas-Solid Reactions[END_REF] . This can be achieved by using a high gas velocity of reactive gases that don't affect the reaction rate. To determine the optimum gas velocity on the oxychlorination rate of MgO, isothermal tests were conducted at 950 °C, using gas mixtures having a linear velocity 'V g ' ranging from 12 to 62 cm/min. Figure 3 gives the evolution of the oxychlorination rate as a function of the gas velocity. The reaction rate increases with the gas velocity up to V g ≈ 40 cm/min. Beyond this value, the rate of oxychlorination was independent of the gas velocity. Therefore, the gas velocity was kept constant and equal to 50 cm/min to study the effect of the other parameters on the oxychlorination rate of MgO.

C. Effect of temperature

The temperature effect on MgO oxychlorination was studied between 850 °C and 1025 °C, using isothermal conditions. An oxychlorinating gas mixture of Cl 2 +O 2 having a molar ratio Cl 2 /O 2 = 4 was used. Figure 4 gives the evolution of pct WL as a function of time for different temperatures. At 850 °C, more than 24 hours were required to reach a reaction extent X * of about 0.8, while the same reaction extent was achieved at 1025 °C for a reaction time of about 1 hour.

* X = weight of the reacted fraction/ initial weight

The data of Figure 4 were used to determine the effect of temperature on the reaction rate for a reaction extent between 0.05 and 0.40. Figure 5 gives the results in an Arrhenius plot. The apparent activation energy 'E a ' of MgO oxychlorination was about 214 kJ/mol. On the other hand, it was observed that the value of the apparent activation energy remained almost constant for reaction extents equal to or lower than 0.70. Such a high value suggests that the overall rate of MgO oxychlorination could be controlled by the rate of the chemical reaction.

Different mathematical equations [START_REF] Szekely | Gas-Solid Reactions[END_REF][START_REF] Gaballah | [END_REF] were tested to fit the experimental data obtained at different temperatures. Best fitting was obtained using Eq. [2] to describe the evolution of X versus time for data of Figure 4. Figure 6 gives curves obtained by fitting the chlorination data (0 ≤ X ≤ 0.7) of MgO using this equation for temperatures between 850 °C and 1025 °C. The correlation coefficient of linearized data was higher than 0.998.

1-(1-X) 1/2 = k t [2]
where k is a constant and t is chlorination time.

D. Effect of partial pressures of the reactive gases

The effects of partial pressure of Cl 2 +O 2 , Cl 2 and O 2 were determined in isothermal conditions in order to evaluate the apparent reaction orders with respect to these gases. The partial orders with respect to the reagents could be equal to the stoichiometric coefficients of the reaction only in the case of elemental reactions. For complex reactions the reaction orders differ from these stoichiometric coefficients. They are obtained experimentally and they may be integer, improper fraction, nil, positive or negative numbers [START_REF] Pannetier | Chimie Générale-Cinétique Chimique[END_REF] .

Effect of partial pressure of Cl 2 +O 2 gas mixture

The effect of partial pressure of Cl 2 +O 2 ranging from 0.25 to 1.00 atm was studied at 950 °C using a Cl 2 +O 2 +N 2 gas mixture having a constant Cl 2 /O 2 ratio equal to 4. This ratio was chosen because it was the optimum value for the selective oxychlorination of the chromite concentrate [1][START_REF] Gaballah | Final Report of European Union Contract No[END_REF][START_REF] Kanari | Proc. of the 126 th TMS Annual Meeting 'EPD Congress 1997[END_REF] .

Its use in this case facilitates the comparison of kinetic data of MgO oxychlorination with those of chromite concentrate. Figure 7 (a) shows the evolution of pct WL as a function of time whilst the relationship between the reaction rate (0.05 ≤ X ≤ 0.40) and P (Cl 2 +O 2 ) is illustrated in Figure 7 (b). An apparent reaction order with respect to Cl 2 +O 2 of about 0.65 was deduced from this plot.

Effect of partial pressure of chlorine

This effect was also studied at 950 °C keeping the oxygen partial pressure constant and equal to 0.33 atm while varying the chlorine partial pressure from 0.20 to 0.67 atm. Nitrogen was used as a dilution gas and the total gas velocity was about 50 cm/min. The obtained isotherms were plotted in Figure 8 (a), showing the evolution of pct WL versus time. Figure 8 (b) indicates that the reaction rate was almost proportional to P Cl 2 since the apparent reaction order with respect to chlorine is close to one. On the other hand, n Cl 2 was higher than n (Cl 2 +O 2 ) (Figure 7 b and 8 b).

Effect of partial pressure of oxygen

Having determined the effects of Cl 2 +O 2 and Cl 2 partial pressures on the reaction rate it was worthwhile to examine the influence of oxygen partial pressure. The tests were carried out at constant chlorine partial pressure (0.33 atm) while the partial pressure of O 2 was varied from 0.20 to 0.67 atm. Figure 9 The apparent reaction order with respect to oxygen was about -0.37 indicating that the oxygen had a negative effect on the chlorination rate of MgO. This result was expected and agrees with Eq. 1.

The effects of (Cl

IV. DISCUSSION

As mentioned previously, this study was performed to explain the kinetics of chromite oxychlorination. More details concerning the reaction kinetics of chromite with Cl 2 +O 2 are available in References 1 through 4. However, some typical curves of chromite oxychlorination were displayed in Figure 11. It was noticed that the chlorination of chromite was characterized by a rapid reaction rate followed by slow one for reaction extents higher than 0.35. As an example at 1000 °C, 35 pct of the sample reacted after a reaction time of 15 minutes, while about 12 hours were necessary to reach a reaction extent of 70 pct. This phenomenon was explained by a different reactivity of the two main components of the chromite, FeCr 2 O 4 and MgCr 2 O 4 , towards oxychlorination gas mixture. The apparent activation energy for reaction extents higher than 0.40 was calculated and corresponds to the oxychlorination of MgCr 2 O 4 [START_REF] Kanari | Light Metals[END_REF] .

Figure 12 shows the Arrhenius plots for the oxychlorination of MgO, Cr To confirm that the kinetics parameters of MgCr 2 O 4 oxychlorination are controlled by the MgO oxychlorination ones, the effects of the partial pressure of the reactive gases on the reaction rate for the oxychlorination of MgO, Cr 2 O 3 [START_REF] Gaballah | Proc. of the 126 th TMS Annual Meeting 'EPD Congress 1997[END_REF] and MgCr 2 O 4 [START_REF] Kanari | Light Metals[END_REF] are compared in Table I [START_REF] Kanari | Light Metals[END_REF] 0.43 0.64 -0.25 0.39

LeRoy [START_REF] Leroy | Production of high purity magnesium oxide by selective chlorination[END_REF] used the low reactivity of magnesium oxide towards chlorine in order to produce high-purity MgO by selective chlorination of natural calcined magnesite at temperatures ranging from 500 °C to 800 °C. The impurities such as CaO and Fe 2 O 3 , contained in the MgO, were converted to chlorides. Ferric chloride was volatilized and calcium chloride was separated by water washing of the residue. The MgO losses ranged from 0.09 to 6.28 pct according to this author. The following paragraphs describe the possibility to remove iron oxides from MgO by oxychlorination.

Figure 13 shows the TGA plots of pct WL versus temperature for the MgO and Fe 2 O 3 during their treatment in Cl 2 +O 2 using non isothermal conditions. This figure clearly shows that the reactivity of Fe 2 O 3 towards an oxychlorinating gas mixture was higher than that of MgO. Almost full chlorination of hematite and volatilization of reaction products were achieved at about 825 °C while the loss of magnesium oxide was negligible.

More consistent information was obtained during the reaction of these solids with Cl 2 +O 2 using isothermal tests. These oxides were chlorinated at 850 °C and 950 °C. Figure 14 illustrates the evolution of pct WL as a function of time. Hematite sample was almost fully chlorinated for a reaction time of about 50 and 10 minutes, at 850 and 950 °C, respectively. Using the same conditions, the weight losses of MgO were less than 5 pct. In our experimental conditions, CaO could be chlorinated even in presence of high oxygen partial pressure in the oxychlorinating gas mixture [START_REF] Kanari | Proc. of the 127 th TMS Annual Meeting 'EPD Congress 1998[END_REF] thus leading to the purification of MgO from its main impurities.

V. CONCLUSIONS

Different kinetics parameters of the MgO oxychlorination were investigated using the TGA technique. Results indicate that the oxychlorination of magnesium oxide between 850 °C and 1025 °C was characterized by an apparent activation energy of about 214 kJ/mol indicating that the overall rate of the process was probably controlled by the rate of the chemical reaction.

The reaction rate of MgO oxychlorination at 950 °C was almost proportional to the partial pressure of chlorine. The reaction orders, at this temperature, with respect to Cl 2 +O 2 and O 2 were 0.65 and -0.37, respectively. Negative value of n O 2 is consistent with the nature of the following reaction :

MgO + Cl 2 → MgCl 2 + 1/2 O 2
Experimental results also demonstrated that the oxychlorination could be used as an efficient method to purify MgO from iron oxides. Almost full extraction of iron can be achieved between 850 °C and 950 °C with MgO losses less than 5 pct. Figure 12 : Arrhenius plots of MgO, Cr 2 O 3 [START_REF] Gaballah | Proc. of the 126 th TMS Annual Meeting 'EPD Congress 1997[END_REF] and MgCr 2 O 4 [START_REF] Kanari | Light Metals[END_REF] oxychlorination. [START_REF] Gaballah | Proc. of the 126 th TMS Annual Meeting 'EPD Congress 1997[END_REF] and MgCr 2 O 4 [START_REF] Kanari | Light Metals[END_REF] oxychlorination.

Kinetics of Oxychlorination of

  (a) traces the evolution of pct WL versus time during the oxychlorination of MgO at 950 °C. This figure clearly shows that the time required to reach a given reaction extent increases with the augmentation of the partial pressure of O 2 in the gas mixture. The evolution of the reaction rate as a function of the oxygen partial pressure is illustrated in Figure 9 (b).

  Magnesium Oxide N. KANARI, E. ALLAIN * , and I. GABALLAH Mineral Processing and Environmental Engineering team, LEM a , CNRS b UMR 7569, ENSG c , INPL d , BP 40, 54501 Vandoeuvre, France * University of Missouri-Rolla, School of Mines and Metallurgy, Center for Pyrometallurgy 210 Fulton Hall, 65401 Rolla, MO, USA

Figure 1 :

 1 Figure 1 : Results of (a) SEM and (b) XRD of the magnesium oxide sample.

Figure 2 :

 2 Figure 2 : Effect of Cl 2 /(Cl 2 +O 2 ) molar ratio on the MgO weight loss.

Figure 3 :

 3 Figure 3 : Effect of gas velocity on the oxychlorination rate of MgO at 950 °C.

Figure 4 :

 4 Figure 4 : Isotherms of oxychlorination of magnesium oxide for (a) T ≤ 925 °C and (b) T ≥ 950 °C.

Figure 5 :

 5 Figure 5 : Arrhenius plot of oxychlorination of MgO by Cl 2 +O 2 .

Figure 6 :

 6 Figure 6 : Fitting of oxychlorination data of MgO using Eq. [2] for (a) T ≤ 925 °C and (b) T ≥ 950 °C.

Figure 7 :

 7 Figure 7 : Isotherms of oxychlorination of MgO (a), using various Cl 2 +O 2 partial pressures and (b) apparent reaction order with respect to Cl 2 +O 2 .

Figure 8 :

 8 Figure 8 : Isotherms of oxychlorination of MgO (a), using various Cl 2 partial pressures and (b) apparent reaction order with respect to Cl 2 .

Figure 9 :

 9 Figure 9 : Isotherms of oxychlorination of MgO (a), using various O 2 partial pressures and (b) apparent reaction order with respect to O 2 .

Figure 10 :

 10 Figure 10 : Effect of gases composition on the reaction rate of oxychlorination of MgO.

Figure 11 :

 11 Figure 11 : Typical isotherm plots for the oxychlorination of chromite.

Figure 13 :

 13 Figure 13 : Non isothermal chlorination of MgO and Fe 2 O 3 using Cl 2 +O 2 .

Figure 14 :Fig. 1 : 4 Fig. 2 :Fig. 3 :Fig. 4 :Fig. 5 :Fig. 6 :Fig. 7 .

 1414234567 Figure 14 : Oxychlorination of MgO and Fe 2 O 3 at 850 °C and 950 °C.

Fig. 8 :Fig. 9 :Fig. 10 :Fig. 11 :Fig. 12 :

 89101112 Fig. 8 : Isotherms of oxychlorination of MgO (a), using various Cl 2 partial pressures and (b) apparent reaction order with respect to Cl 2 . (a) (b)

Fig. 13 : 4 Fig. 14 :

 13414 Fig. 13 : Non isothermal chlorination of MgO and Fe 2 O 3 using Cl 2 +O 2 .

  2 +O 2 ), Cl 2 and O 2 partial pressures on the oxychlorination rate of MgO at 950 °C were compared in Figure 10. The significance of the effects of reactive gases partial pressures, expressed as apparent reaction orders, clearly follows the following sequence : n Cl 2 > n (Cl 2 +O 2 ) > n O 2 . The apparent reaction order of n (Cl 2 +O 2 ) was almost equal to the algebraic sum of those of n Cl 2 and n O 2 .

  . The reaction order with respect of oxygen had a positive value only for the oxychlorination of Cr 2 O 3 because it favors the formation of CrO 2 Cl 2 according to Eq. 3. The apparent reaction orders with respect to Cl 2 +O 2 , Cl 2 and O 2 for the oxychlorination of MgO and MgCr 2 O 4 contained in chromite were closer, compared with those of Cr 2 O 3 . This suggests that the oxychlorination kinetics of MgCr 2 O 4 is mostly controlled by the oxychlorination rate of MgO.

	1/2 Cr 2 O 3	+ Cl 2 + 1/4 O 2	→	CrO 2 Cl 2	[3]
			Table I. Comparison of reaction orders with respect to reactive
			gases of different solids' oxychlorination.
	Solids		n (Cl 2 +O 2 )	n Cl 2	n O 2	n Cl 2 + n O 2
	Cr 2 O 3 [12]	1.29	1.08	0.23	1.31
	MgO		0.65	0.98	-0.37	0.61
	MgCr 2 O 4

ACKNOWLEDGMENTS

This work was performed in the frame of contract N° BRE2-CT92-0173 thanks to the financial support of the European Union (DG-XII). The authors thank Dr. H. L. Schmidt for discussion, suggestions and help. The authors are indebted to Dr. A. Bonazébi, Drs. M. Djona, S. Ivanaj, N. Menad, N. Mirghaffari for discussions and help on different subjects and to Mrs. C. Richard for the technical and administrative support.


K Abdelouahab

M Pelcat

J Serot

F Berry

C Bourrasset

J C Quinton

Hardware Automated Data ow Deployment of CNNs

Deep Convolutional Neural Networks (CNNs) are the state of the art systems for image classi cation and scene understating. However, such techniques are computationally intensive and involve highly regular parallel computation. CNNs can thus bene t from a signi cant acceleration in execution time when running on ne grain programmable logic devices. As a consequence, several studies have proposed FPGAbased accelerators for CNNs. However, because of the huge amount of the required hardware resources, none of these studies directly was based on a direct mapping of the CNN computing elements onto the FPGA physical resources. In this work, we demonstrate the feasibility of this so-called direct hardware mapping approach and discuss several associated implementation issues. As a proof of concept, we introduce the 2 open source tool, that is able to automatically transform a CNN description into a platform independent hardware description for FPGA implementation.

Introduction

Convolutional Neural Networks (CNNs) [START_REF] Lecun | Gradient Based Learning Applied to Document Recognition[END_REF] have become a de-facto standard that increased the robustness and accuracy of machine vision systems. It is possible nowadays to build high performance image classi cation systems by deploying large-scale, pre-trained CNNs models. However, this accuracy comes at the price of a high computational cost as state of the art CNNs may require up to 40 GOPs to classify a single frame [START_REF] Canziani | An Analysis of Deep Neural Network Models for Practical Applications[END_REF]. As a result, implementing CNNs with real-time constraints is challenging task. A possible way to address this challenge is to take advantage of the massive ne grain parallelism o ered by FPGA devices to embody the large amount of intrinsic parallelism exhibited by CNN-based algorithms. In this case, the problem boils down to nd an adequate and e cient mapping between the computation model of the latter and the execution model supported by the former. Based on our previous experience in the implementation of real-time vision applications on FPGA-based platforms [START_REF] Sérot | High-Level Data ow Programming for Recon gurable Computing[END_REF], we advocate the use of a stream-based data ow model to solve this mapping problem. In this approach, a CNN-based algorithm is described as graph of data ow actors exchanging data through unidirectional channels and this graph is statically and physically mapped onto the target FPGA using a library of pre-de ned computing elements to implement actors.

In the sequel, we demonstrate the feasibility of this so-called Direct Hardware Mapping (DHM) approach for implementing realistic CNN-based applications onto Field-Programmable Gate Arrays (FPGAs). Moreover, we introduce 2, a software framework providing a fully automated implementation path for CNNs onto FPGAs using the DHM approach. The 2 tool is compatible with the widely used Ca e deep learning framework [START_REF] Jia | Ca e: Convolutional Architecture for Fast Feature Embedding[END_REF] and generates platform independent synthetizable VHDL code. In other words, we introduce in this work a tool that automatically maps a Ca e pre-trained model onto an FPGA device.

CNNs : Computations and parallelism sources

CNNs are a category of feed forward arti cial neural networks that are bio-inspired by the visual cortex of the brain. The huge improvement of CNN-based algorithms was made possible by two factors: On one hand, the availability of massive-sized annotated image data-sets [START_REF] Deng | Imagenet: A large-scale hierarchical image database[END_REF] allowed to train robust large scale feature extractors and accurate classi ers. On the other hand, the growth of high performance processors and, especially Graphics Processing Units (GPUs), provided the computational power required to train deeper and more complex neural networks [START_REF] Chetlur | cuDNN: E cient Primitives for Deep Learning[END_REF]. A typical CNN structure, as shown in gure 1, will perform a succession of convolutions interspersed with sub-sampling layers. The last stages include typically two or three fully connected neural network for classi cation tasks. The depth (number of layers) of a CNN ensures better accuracy and less over-tting and recent networks are usually very deep with 8 to 19 layers [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]. two subsampling layers and one fully connected stage (FC).

Convolution layers

Convolutional layers are the most computationally intensive and are responsible -in a typical implementation -for more than 90% of the CNN execution time [START_REF] Cong | Minimizing computation in convolutional neural networks[END_REF]. Each layer (l) extracts N feature maps from C input channels by performing N convolutions of size K × K on each input. This ltering is followed by the application of a non-linear activation function act and a bias term b n to each set of features. As shown in equation 1, N × C convolutions are required to process a given layer.

f (l ) [n, i, j] = b (l ) [n] + C c=1 K p=1 K q=1 Φ (l ) [c, i + p, j + q].w (l ) [n, c, p, q] (1)
where • f (l ) is a tensor of output feature maps of layer (l)

• b (l ) [n] is the bias term applied to applied to feature n

• Φ (l ) is a tensor of input feature maps of layer (l)

• w (l ) is tensor of pre-learned lters

As already pointed out in [START_REF] Motamedi | Design space exploration of FPGA-based Deep CNNs[END_REF], the computations described in equations 1 exhibit a large amount of potential parallelism:

• Inter Layer parallelism: CNNs have a feed-forward hierarchical structure consisting of a succession of data-dependent layers. Layers can therefore that can be executed in a pipelined fashion where the execution of layer (l) can start before the execution of layer (l -1) ends.

• Inter neuron parallelism: Each neuron of a layer is independent when processing features. Thereby, a full data-parallelism can be exploited when computing concurrently each of the N (l ) element of equation 1

• Inter convolution parallelism: All of the convolutions performed by a single neuron can also be evaluated simultaneously by computing concurrently the C (l ) convolutions of equation 1.

• Intra convolution parallelism: 2D image convolution can be implemented in a pipelined fashion [START_REF] Richard | Parameterized convolution ltering in a eld programmable gate array[END_REF] allowing the K × K multiplications to be computed concurrently in equation 1

Subsampling layers

A common operation when conceiving CNNs is to periodically insert subsampling (or pooling) layers inbetween successive convolutional layers. These downsample the inputs by selecting the average, or, more commonly, the maximum of a given neighborhood of each pixel as described in equation 2

f (l ) [n, i, j] = max p,q ∈[1:K ] Φ (l ) [n, i + p, j + q] (2)
Pooling layers reduce the amount of parameters required to process the next stages of the network, which controls over tting in one hand and decrease the computation load on the other.

Fully connected layers

A Fully Connected (FC) neural network -with usually 3 or 4 hidden layers-terminates CNNs and acts as a classi er. In this case, no parameters are shared across the feature-maps (feature maps and learned parameters have the same dimension). In this case, FC layer activations are computed with the inner product operation followed by a bias o set as detailed in equation 3

∀l = 1 : L (Number of FC layers) ∀n = 1 : N (Number of output feature maps) f (l ) [n] = act       b (l ) [n] + C (l ) c=1 < ϕ (l ) [c], w (l ) [n, c] >       (3)
where <, > denotes the the inner product operator.

3 Direct Hardware Mapping of CNN entities

Data ow processing of CNNs

The foundations of data ow Models of Computation (MoC) were formalized by [START_REF] Jack | A Preliminary Architecture for a Basic Data-ow Processor[END_REF] in order to create an architecture where multiple fragments of instructions can process simultaneously a stream of data. Programs respecting data ow semantics are described as a network (graph) of fundamental processing units commonly called actors and communicating abstract data messages called tokens on unidirectional First-In First-Out (FIFO) channels.

In terms of architecture-application matching, the CNN's layout ts naturally with a stream-based model of computation. In other words, the operations involved in feed forward propagation of a CNN -described in the latter section-can be executed following the stream-based data ow MoC. In fact, CNNbased algorithms can be modeled as modeled as data ow process networks (DPNs) where nodes correspond to processing actors and edges correspond to communication channels. Each actor follows a purely datadriven execution model where execution ( ring) is triggered only by the availability of input operands.

The DHM approach consists of physically mapping entirely graph of actors onto the target device. Each actor becomes a computing unit with its speci c instance on the FPGA and each edge is mapped to a signal.

DHM of Convolution layers

As stated in section 2.1, convolutional layers are the most computation intensive tasks in a given network. However, DHM approach fully exploits all the parallelism sources of theses layers. All neurons of a layer are mapped on the device to take advantage of intra-neuron parallelism (Fig 2-a). In neurons, each convolution is mapped separately (Fig 2-b) and nally, within a convolution engine, each multiplier is instantiated separately (Fig 2-c). As an example, gure 3 illustrates how a convolution layer C1 (C = 3, N = 5, K = 3) extracts 5 features from a 3-channel input pixel ow. In this example, 15 convolution and 5 activation blocks are mapped onto the FPGA as a result of the layer graph transformation, which corresponds to 135 multiplications, 20 summations and 5 activations. Direct Hardware Mapping of CNNs completely removes the need for an external memory to store intermediate results or parameters. Moreover, thanks to the fully pipelined execution model, the global throughput is only limited by the maximum clock frequency. However, these advantages come at the cost of a high resource consumption since the whole graph has to mapped onto the physical resources of the FPGA. In certain cases, this could limit the complexity of the CNNs that can be handled by the DHM approach. It is crucial, therefore, to ensure that the core operations involved in CNN actors can be translated e ciently in hardware. The most important issues, by far, are those related to on-chip memory requirements on one hand, and the implementation of arithmetic operators on the other hand.

ϕ 0 ϕ 1 . . . ϕ C η 0 η 1 . . . η N f 0 f 1 f N (a) ϕ C . . . ϕ 1 ϕ 0 b 0 conv 0C . . .
Σ 4 Σ 3 Σ 2 Σ 1 Σ 0 act 4 act 3 act 2 act 1 act 0 f (C1) 4 f (C1) 3 f (C1) 2 f (C1) 1 f (C1) 0

Neighborhood extraction

Data ow-based processing of convolutions -such in [START_REF] Richard | Parameterized convolution ltering in a eld programmable gate array[END_REF]-can be divided into 2 parts: neighborhood extraction (NE) and Multiply-ACCumulation (MAC). Neighborhood Extraction (NE) relies on bu ers to grant a full access to the K (l ) × K (l ) neighbors of each pixel (as shown in gure 4). Such an architecture is advantageous since it can directly extract the neighborhood of streams of pixels each clock-cycle.

Multiply Accumulate (MAC) performs a multiplication of neighborhood pixels with pre-learned kernels then accumulates the result to output feature maps. As long as the access to full neighborhood pixels is guaranteed, each of the multiplications of can be performed in a parallel way using K (l ) × K (l ) multipliers (as shown in Fig 2-c). Combining NE and parallel MAC strategy fully exploits the intra Kernel parallelism of CNNs which grants high acceleration to convolutions and, consequently, the feature extraction process.

Neighborhood Extraction Factorization (NEF)

When adopting a DHM approach, it is possible to factorize the neighborhood extraction process in order to optimize the memory print of convolutional layers. In this case, it is possible to rely only on on-chip memory bu ers to process a hole convolutional layer. Thus, since multiple neurons in a given layer have same input features to process (only the convolution kernels change), the neighborhood extraction entity can be factorized for each input feature map. This will divide the memory requirements of each layer by a factor N (l ) (cf gure 5). For instance, the rst layer of the AlexNet CNN (N=96,C=3,K=11) would require 96×3×11×11 = 34KB of bu er memory to be processed, while a factorization of neighborhood extractors leads to 96 times less memory requirements (0.3KB). Full results of NEF on Alexnet layers are detailed in gure 6. Several studies [START_REF] Gupta | Deep Learning with Limited Numerical Precision[END_REF][START_REF] Gysel | Hardware-oriented Approximation of Convolutional Neural Networks[END_REF] have demonstrated that CNNs, and more generally deep learning applications, usually tolerate approximate computations with short xed-point arithmetic. Frameworks such as Ristretto [START_REF] Gysel | Hardware-oriented Approximation of Convolutional Neural Networks[END_REF], for example, can perform ne-tuning of data representation in order to support xed-point numerical representations with variable data lengths. In particular, an 8-bit (resp. 2-bit) precision is su cient to infer the AlexNet [START_REF] Krizhevsky | ImageNet Classi cation with Deep CNNs[END_REF] (resp. LeNet [START_REF] Lecun | Gradient Based Learning Applied to Document Recognition[END_REF]) CNNs with little to no degradation in classi cation accuracy. The DHM approach advocated in this paper can indeed take advantage of this to signi cantly reduce the amount of required hardware resources by rst inferring the minimal required precision and then deriving the size of the hardware resources to exactly match this precision.

Multiplications with Logic Elements

Convolutions require many multiplications. If these multiplications are implemented using hardwired Digital Signal Processing (DSP) blocks within the target FPGA, this dramatically limits the complexity of the CNN that can be implemented. For instance, the second layer of the LeNet5 network (C = 6, N = 16, K = 5) requires 2400 multipliers. This number largely exceeds the number of DSP blocks provided by many FPGAs and, especially by embedded devices. We overcome this problem by systematically forcing the synthesis tool to implement multiplications with logical elements instead of DSP blocks, leading the resulting implementations to rely on AND gates and trees of half-adders [START_REF] Altera | Implementing Multipliers in FPGA Devices[END_REF].

In addition, we take advantage of the fact that in the case of CNNs the convolution kernels -and hence the second operand of most of multiplications -are actually constants (derived from the o ine training stage). It is therefore possible to use a specialized version for those multiplier instances. While this approach limits the exibility of the system -it requires to re-compile and re-synthesise the VHDL design whenever parameters values are changed -, it delegates to the synthesis tool the task to perform low-level area and performance optimizations. More particularly, multiplications by 0 (resp 1) are removed (resp. replaced by a simple signal connection) and multiplications by a power of 2 are implemented using shift registers. 5 The H

Multiplicand

utility

The H 2 framework is set of tools built upon the principles and optimization techniques described in the previous section. It is capable of automatically generating a platform independent hardware description of a CNN from a Ca e model [START_REF] Jia | Ca e: Convolutional Architecture for Fast Feature Embedding[END_REF]. First, layer speci cations (Layer type, Number of input channels C, Number of output features N , kernel size K) are extracted from the Ca e model and the learned parameters are read, rounded to a xed-point representation format and written as generic parameters in a con guration le. Second, a top-level VHDL le is created by transforming the data ow graph described in Ca e. The top-level instantiates a set of generic layers parametrized according to the Ca e model specications. These layers are described using a small number of basic prede ned actors. These actors, written in a structural VHDL, follow the data ow execution semantics discussed in the latter sections. The output is a platform independent VHDL code that can be implemented on the FPGA device using the adequate synthesis tool. The H 2 framework and the library of CNN actors supporting the DHM approach are open-source and available. online 1

_ WIDTH => CONV2 _ IMAGE _ WIDTH, KERNEL _ SIZE => CONV2 _ KERNEL _ SIZE, NB _ IN _ FLOWS => CONV2 _ IN _ SIZE, NB _ OUT _ FLOWS => CONV2 _ OUT _ SIZE, KERNEL _ VALUE => CONV2 _ KERNEL _ VALUE, KERNEL _ NORM => CONV2 _ KERNEL _ NORM, BIAS _ VALUE => CONV2 _ BIAS _ VALUE ) port map( clk => clk, reset _ n => reset _ n, enable => enable, in _ data => pool1 _ data, in _ dv => pool1 _ dv, in _ fv => pool1 _ fv, out _ data => conv2 _ data, out _ dv => conv2 _ dv, out _ fv => conv2 _ fv );
...

Experimental Results with Haddoc2

As a proof of concept, we have implemented, using the H 2 framework, FPGA-based accelerators for three CNN-based applications, listed in Table 2. The rst one is the Ca e version of the LeNet5 [START_REF] Lecun | Gradient Based Learning Applied to Document Recognition[END_REF] CNN that requires 20.78 MOPs to process a frame of size 28x28. The second application is the face detector used in [START_REF] Farabet | An FPGA-Based Stream Processor for Embedded Real-Time Vision with CNNs[END_REF] which requires 622.08 MOPs to process a 320x240 frame. The last one is introduced in [START_REF] Heikki Huttunen | Car type recognition with Deep Neural Networks[END_REF] to perform car type classi cation and requires 268.28 MOPs to process 96x96 frames. The two rst CNNs have been trained using Ca e while the third model has been directly downloaded as a Ca e pre-trained model. Table 2 gives parameter values for each CNN convolutional layer. LeNet5 and CarType CNNs have 2 convolutional layers while FaceDetect has 3. The corresponding hardware descriptions of each network have been automatically generated using Haddoc2 on an Intel i7-4770 CPU and were synthesised on two FPGA devices using respectively Intel Quartus 16.1 and Xilinx Vivaldo 2016.4. Table 3 reports post-tting results of the LeNet-5 accelerator on an embedded Intel Cyclone V 5CGXFC9E7 device using 3 implementation strategies. In the rst case, only DSP blocks are used to map the CNN multiplications. The resulting hardware requires 72× the available resource of the device. The second case features an implementation of multiplication based on logic elements and requires 3.8× the available logic. Using tailored multipliers reduces resources by a factor of 8.6×, tting the CNN accelerator onto an Intel Cyclone V device. 4 details post tting results on two embedded FPGA platforms: the Intel Cyclone V 5CGXFC9E7 and the Xilinx Kintex7 XC7Z045FBG. To the best of our knowledge, these numbers are the rst to demonstrate the applicability of a DHM-based approach for the implementation of CNNs on embedded FPGAs. The three hardware accelerators t onto the embedded devices with no o -chip memory requirement. The memory footprint shown in post tting reports corresponds to line bu ers used by the data ow-based convolution engine and both synthesis tools instantiate LUT-based memory blocks to implement these bu ers. As expected when using DHM, the logic utilization in the FPGA grows with the the topology of the CNN. However, in all the studied cases, the resources are su cient to support direct hardware mapping. Finally, the same table reports timing analysis results of the three generated hardware accelerators. With a peak frequency of 62.3 MHz for the CarType CNN, DHM grants a maximum computation throughput of 1813 GOPs/s. For the face detection neural network, the presence of a third convolutional layer in the pipeline drops the maximum frequency to 56.7 MHz (i.e 357 GOPs/s) in the Cyclone device, which corresponds to 164 classi cations/sec on 512x512 images with a 3-multiscale pyramid. Several studies leverage on FPGA computational power and hardware exibility to implement the feedforward propagation of CNNs. A non exhaustive review of these can be found in [START_REF] Gri N Lacey | Deep Learning on FPGAs: Past, Present, and Future[END_REF]. In most of approaches, acceleration of CNN-based applications is provided by mapping a limited subset of processing elements onto the target device. This is the case for example in [START_REF] Ovtcharov | Accelerating Deep Convolutional Neural Networks Using Specialized Hardware[END_REF] where authors describe an accelerator for the AlexNet CNN [START_REF] Krizhevsky | ImageNet Classi cation with Deep CNNs[END_REF] implemented on a large Stratix V FPGA which, to the best of our knowledge, outperforms most state-of-the-art implementations in terms of computational and outperformed most of state-of-the-art implementations such [START_REF] Chakradhar | A Dynamically Con gurable Coprocessor for Convolutional Neural Networks[END_REF][START_REF] Farabet | CNP: An FPGA-based processor for Convolutional Networks[END_REF][START_REF] Peemen | Memory-centric accelerator design for Convolutional Neural Networks[END_REF]. Most of these designs are FPGA based accelerators for convolution with a relatively similar architecture of parallel processing elements associated with embedded hardcore processors running a software layer. Other approaches like [START_REF] Zhang | Optimizing FPGAbased Accelerator Design for Deep Convolutional Neural Networks[END_REF] relies on analytical design scheme using the roo ine model and loop tiling to propose an inference engine where the attainable computation roof of the FPGA is reached. This loop tilling optimization is performed on a C code then implemented in oating point on a Virtex 7 485T using Vivaldo HLS Tool.

As it has been seen in the latter sections, feed forward propagation is an algorithm that intrinsically suits to data ow processing. Thus, dedicated stream processors for CNNs have been proposed. The most notable contribution was neuFlow [START_REF] Farabet | NeuFlow: A runtime recon gurable data ow processor for vision[END_REF]: A runtime recon gurable processor for real-time image classi cation. In this work, Farabet and al. introduced a grid of processing tiles that were con gured on runtime to build a data ow graph for CNN applications. It was associated to "luaFlow": a data ow compiler that transforms a high-level ow-graph representation of an algorithm (in a Torch environment [START_REF] Collobert | Torch. NIPS Workshop on Machine Learning Open Source Software[END_REF]) into machine code for neuFlow. Such architecture was implemented on a Virtex 6 VLX240T and provided a 12 fps categorization for 512x375 images. Thus, NeuFlow transformed a CNN graph into a set of data ow instructions, where each instruction is described as an hardware con guration of 2D-processing elements called Processing tiles (PTs). Execution of the graph is carried out by sequencing the instructions on the target FPGA. This approach requires an external memory to store intermediate results, which in turn, even with the help of a DMA, limits the nal speedup.

By contrast, the DHM approach and H 2 tool introduced in the present work performs all processing on the y and does not require an external memory to store intermediate results. Throughput is therefore not limited by o -chip memory bandwidth. Previous works in [START_REF] Abdelouahab | A Holistic Approach for Optimizing DSP Block Utilization of a CNN Implementation on FPGA[END_REF] describe a rst version of Haddoc that relied on the Caph [START_REF] Sérot | High-Level Data ow Programming for Recon gurable Computing[END_REF] , a High-Level Synthesis (HLS) tool to provide data ow-based hardware accelerators for CNNs on FPGAs. While this implementation operated at very high frame-rates (800 classi cations/sec on 256 × 256 images), the over-head that comes with the HLS heavily restrained the size of CNNs to be implemented.

Figure 1 :

 1 Figure 1: An example of a CNN topology with 3 convolutional layers (C1,C2,C3)two subsampling layers and one fully connected stage (FC).

∀l = 1 :

 1 L (Number of con layers) ∀n = 1 : N (Number of output feature maps) ∀i =1 : Ix (feature map rows) ∀j = 1 : I ( feature map columns)

∀l = 1 :

 1 L (Number of pool layers) ∀n = 1 : N (Number of output feature maps) ∀i =1 : Ix (feature map rows) ∀j = 1 : I ( feature map columns)

Figure 2 :

 2 Figure 2: The 3 levels of DHM implementation of CNN entities: (a) in convolution layers, (b) in neurons, (c) in convolution engines

Figure 3 :

 3 Figure 3: Applying the 3 levels of DHM ( g 2) to a dummy convolutional layer C1 (N=5, C=3, K=3): 15 separate convolution engines (135 Multipliers and 15 adders) plus 5 adders and 5 activation blocks are required to process the layer in a full parallel fashion. (bias omitted)

Figure 4 :

 4 Figure 4: Architecture of a 3 × 3 neighborhood extractor : 2 Bu ers with image length size are required to perform a 3 × 3 convolution on streams of pixels p i j

Figure 5 : 3 conv1

 53 Figure 5: Data-path of a convolutional layer (bias omitted): The factorization of neighborhood extraction process reduces the memory bu ers by a factor of 5 when compared to gure 3

Figure 6 :

 6 Figure 6: Ratio of memory requirements between architectures w/ and wo/ NEF for Alexnet convolutional layers: 390% less memory is required when factorizing the neighborhood extractors

Figure 7 :

 7 Figure 7: Hardware generation: the CNN layer arrangement is described in the top-level les while kernel parameter values and layer speci cation are written on the con guration le.

Table 1 :

 1 Resource utilization of a random 3 × 3 convolution engine on an Altera Cyclone V device with di erent implementations.

		Variable		Constant
	LE Based	ALM: 380 (0.67 %) ALM : 121 (0.21 %) DSP : 0 (0 %) DSP : 0 (0 %)
	DSP Based	ALM : 71 (0.12 %) DSP : 10 (6.41 %)	ALM : 70 (0.12 %) DSP : 7 (4.48 %)
	.prototxt			toplevel.vhd
		H	2
	.ca emodel			params.vhd
	Ca e			Hardware

  .

	Listing 1: Ca e description of a con layer	Listing 2: Generated VHDL code of the layer
	name: "LeNet"	...
		architecture RTL of lenet is
	...
		...
	layer {
	name: "conv2"	conv2: convLayer
	type: "Convolution"	generic map(
	bottom: "pool1"	PIXEL _ SIZE	=> PIXEL _ SIZE,
	top: "conv2"	IMAGE
	param {
	lr _ mult: 1
	}
	param {
	lr _ mult: 2
	}
	convolution _ param {
	num _ output: 50
	kernel _ size: 3
	stride: 1
	weight _ filler {
	type: "xavier"
	}
	bias _ filler {
	type: "constant"
	}
	}
	}
	...

Table 2 :

 2 Topology of the convolutional layers of studied CNNs.

		LeNet5 [19] FaceDetect [12] CarType [15]
	Input size	28 x 28	320 x 240	96 x 96 x3
	Layer parameters N C K N C	K	N C	K
	conv1+maxpool	20 1 5 6	1	7	32 3	5
	conv2+maxpool	50 20 5 10 6	7	32 32	5
	conv3	---30 10	3	--	-
	Kops/Pixel	26.5	6.3		29.1

Table 3 :

 3 Resource utilization by a DHM LeNet5 CNN with di erent implementations strategies for multipliers.

		DSP-based	LE-based	LE-based + const.
	Logic Usage (ALM)	NA	433500 (381%)	50452 (44%)
	DSP Block usage	24480 (7159 %)	0 (0%)	0 (0%)
	Table

Table 4 :

 4 Resource Utilization of the Haddoc2-generated convolutional layers of studied CNNs with 5-bit representation on: a-an Intel Cyclone V FPGA, b-a Xilinx Kintex 7 FPGA.

			LeNet5 [19] FaceDetect [12] CarType [15]
		Logic Elements (ALMs) 50452 (44%)	6158 (5%)	48243 (42%)
		DSP Blocks 1	0 (0 %)	0 (0%)	0 (0%)
	a	Block Memory Bits	2752 (1%)	41408 (1%)	28320 (1%)
		Frequency	69.14 MHz	56.7 MHz	66.0 MHz
		Processing capabilities 1832 GOPs/s	357 GOPs/s	1920 GOPs/s
		Slices	48114 (88%)	6221 (11%)	49082 (89%)
		DSP Blocks 1	0 (0%)	0 (0%)	0 (0%)
	b	LUTs as Memory	420 (1%)	1458 (2%)	1154 (1%)
		Frequency	62.13 MHz	44.41 MHz	62.3 MHz
		Processing capabilities 1646 GOPs/s	279 GOPs/s	1813 GOPs/s
	7 Related work

https://github.com/KamelAbdelouahab/haddoc2


Katia Bachi
email: bachi.katia13@gmail.com

Cédric Chauvière
email: cedric.chauviere@uca.fr

Hacène Djellout
email: hacene.djellout@uca.fr

Karim Abbas

PROPAGATION OF EPISTEMIC UNCERTAINTY IN QUEUEING MODELS WITH UNRELIABLE SERVER USING CHAOS EXPANSIONS

Keywords: Unreliable queueing model, Epistemic uncertainty, Chaos expansions, Orthogonal polynomial, Sobol' indices AMS 2010 subject classifications. 90B22, 60K25, 65C05, 65C50

In this paper, we develop a numerical approach based on Chaos expansions to analyze the sensitivity and the propagation of epistemic uncertainty through a queueing systems with breakdowns. Here, the quantity of interest is the stationary distribution of the model, which is a function of uncertain parameters. Polynomial chaos provide an efficient alternative to more traditional Monte Carlo simulations for modelling the propagation of uncertainty arising from those parameters. Furthermore, Polynomial chaos expansion affords a natural framework for computing Sobol' indices. Such indices give reliable information on the relative importance of each uncertain entry parameters. Numerical results show the benefit of using Polynomial Chaos over standard Monte-Carlo simulations, when considering statistical moments and Sobol' indices as output quantities.

Introduction

Queueing systems with breakdowns are widely used to model problems occuring in computer, manufacturing systems, and communication networks. Typically, in queueing models the performances measurements are assessed for deterministic parameters values. However, in practice, the exact values of these parameters are not well known (uncertain), and they are generally estimated empirically from few experimental observations. This lack of precise information propagates uncertainties in the output measure and such study is the purpose of this work.

Since the pioneering work of Thiruvengadam [START_REF] Thiruvengadam | Queuing with breakdowns[END_REF] and Avi-Itzhak and Naor [START_REF] Avi-Itzhak | Some queuing problems with the service station subject to breakdown[END_REF], many authors have been investigating the queueing systems with server breakdowns, see for example [START_REF] Cao | Analysis of an M/G/1 queueing system with repairable service station[END_REF][START_REF] Li | Reliability analysis of M/G/1 queueing systems with server breakdowns and vacations[END_REF][START_REF] Wang | Reliability analysis of the retrial queue with server breakdowns and repairs[END_REF] and references therein. In this paper, we choose a functional approach to the analysis of the dependence of performance measures of certain queues, such as the M/M/1 and the M/G/1 queue with breakdowns, with respect to some input parameters. More specifically, denoting the probability of a server breakdown by θ, we seek to compute the stationary distribution of the queue-length process, denoted by π θ . For a fixed value of θ and a finite waiting capacity, π θ can be found numerically by solving π θ Ξ θ = π θ and π θ = 1, where Ξ θ denotes the transition probability matrix of an embedded jump chain. The definition of the embedded jump chain will depend on the type of the queue: for the M/M/1 queue we use the sample-chain, embedded at appropriate Poisson-times, and for the M/G/1 queue we will embed the chain at departure and repair moments; details will be provided later in the paper. In case of a large or infinite waiting capacity, π θ can be obtained via a Laplace-Stieltjes transform [START_REF] Abate | An introduction to numerical transform inversion and its application to probability models[END_REF]. Solving for π θ involves numerical inversion, which can be computationally demanding [START_REF] Shortle | Waiting-time distribution of M/DN /1 queues through numerical Laplace inversion[END_REF]. The computation of π θ is a challenging problem and a variety of approaches have been proposed in the literature for approximately or indirectly solving the stationary distribution. The predominant approach is to obtain either the generating function of π θ or an analytical expression for π θ containing a Laplace-Stieltjes transform, see, for example, [START_REF] Abate | An introduction to numerical transform inversion and its application to probability models[END_REF][START_REF] Baccelli | Queueing systems with breakdowns in data base modeling[END_REF]. Also numerical solutions by means of the matrix geometric method [START_REF] Neuts | Matrix-geometric solutions in stochastic models[END_REF] are available, see [START_REF] Mitrany | A Many-Server Queue with Service Interruptions[END_REF][START_REF] Neuts | A Markovian queue with N servers subject to breakdowns and repairs[END_REF] for details.

In performance analysis, one is not only interested in evaluating the system for specific set of parameters but also in the sensitivity of the performance with respect to these parameters. For example, in a queueing model with breakdowns, the breakdown probability is a parameter of key interest and, in this paper, we will analyze the dependence of π θ on θ, which is significantly more challenging than evaluating π θ for a fixed θ. Most of the time, it is assumed that these stochastic models are solved for fixed parameters values. However, the parameters of the model are determined through insufficient statistical data (a limited number of observations), leading to uncertainty in the assessment of their values. This parametric uncertainty, induced from the incomplete information of the parameter, is called epistemic uncertainty [START_REF] Limbourg | Dependability Modelling under Uncertainty: An Imprecise Probabilistic Approach[END_REF][START_REF] Mishra | An unobtrusive method for uncertainty propagation in stochastic dependability models[END_REF][START_REF] Winkler | Uncertainty in probabilistic risk assessment[END_REF].

In order to estimate the uncertainty of the parameters in the performance measures of the model, two complementary approaches may be used: uncertainty analysis and sensitivity analysis. Uncertainty analysis consists in modelling input parameters as random variables, see, for example, [START_REF] Cacuci | Sensitivity and uncertainty analysis[END_REF][START_REF] Helton | Treatment of uncertainty in performance assessments for complex systems[END_REF][START_REF] Helton | Uncertainty and sensitivity analysis in the presence of stochastic and subjective uncertainty[END_REF]. Then, the sensitivity analysis aims at determining the relative contribution of individual parameters. Several approaches for the propagation of uncertainty have been developed, including interval arithmetic [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF][START_REF] Rocco | Distribution systems reliability uncertainty evaluation using an interval arithmetic approach[END_REF], Taylor series expansion [START_REF] Dhople | A parametric uncertainty analysis method for Markov reliability and reward models[END_REF][START_REF] Granger | Comparative analysis of uncertainty propagation methods for robust engineering design. Uncertainty: A guide to dealing with uncertainty in quantitative risk and policy analysis[END_REF][START_REF] Ouazine | A functional approximation for retrial queues with two way communication[END_REF][START_REF] Shooman | Probabilistic Reliability: An Engineering Approach[END_REF][START_REF] Takhedmit | A parametric uncertainty analysis method for queues with vacations[END_REF], moments [START_REF] Granger | Comparative analysis of uncertainty propagation methods for robust engineering design. Uncertainty: A guide to dealing with uncertainty in quantitative risk and policy analysis[END_REF][START_REF] Mattia | Comparative analysis of uncertainty propagation methods for robust engineering design[END_REF][START_REF] Shooman | Probabilistic Reliability: An Engineering Approach[END_REF], Monte Carlo analysis [START_REF] Granger | Comparative analysis of uncertainty propagation methods for robust engineering design. Uncertainty: A guide to dealing with uncertainty in quantitative risk and policy analysis[END_REF][START_REF] Mattia | Comparative analysis of uncertainty propagation methods for robust engineering design[END_REF][START_REF] Shooman | Probabilistic Reliability: An Engineering Approach[END_REF] . Overviews of these approaches are available in several reviews [START_REF] Winkler | A Comparative Review of Sensitivity and Uncertainty Analysis of Large-Scale Systems -II: Statistical Methods[END_REF][START_REF] Helton | Uncertainty and sensitivity analysis techniques for use in performance assessment for radioactive waste disposal[END_REF][START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF][START_REF] Iman | An investigation of uncertainty and sensitivity analysis techniques for computer models[END_REF][START_REF] Ionescu-Bujor | A comparative review of sensitivity and uncertainty analysis of large-scale system, I: Deterministic methods[END_REF][START_REF] Ronen | Uncertainty Analysis[END_REF].

Since the seminal work of Ghanem an Spanos in the 90's [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF] and later generalized by Xiu and Karniadakis [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF], Polynomial Chaos has spread over a broad scientific community, see [START_REF] Sepahvand | Uncertainty quantification in stochastic systems using polynomial chaos expansion[END_REF], [START_REF] Chauvière | Computational modeling of uncertainty in time-domain electromagnetics[END_REF]. Then, due to the specific form of these polynomials, it quickly became apparent that this representation was suitable to compute Sobol' indices [START_REF] Sudret | Global sensitivity analysis using polynomial Chaos expansion[END_REF]. To the best of our knowledge, this is the first atempt to use chaos expansion to investigate the propagation of the epistemic uncertainty through the performance measure of the queueing systems.

The main objective of this work is to develop a numerical procedure to investigate the sensitivity and the propagation of epistemic uncertainty in the input parameters of a queueing systems with unreliable server. As a typical example, a customer arrives at the queueing system according to a Poisson process with rate λ, which we consider as an input parameter of the model. We assume that this parameter is not precisely known and therefore it can be model as a random variable of known probability density function. Clearly, the output measure is a function of this input random variable. The epistemic uncertainty is propagated through a functional relationship of the type Y (ω) = g(λ(ω)) that links the input parameter to the output measures. In this work, special attention is given to the output stationary distribution. From a statistical point of view the quantity of interest may be the moments of the functional stationary distribution. Such quantities can easily be computed from the coefficient of the polynomial Chaos representation of g(λ(ω)). Over the past few years, several techniques have emerged to compute the coefficients of polynomial Chaoses. They can be classified in two main families: intrusives methods and non-intrusive methods [START_REF] Pettersson | Polynomial Chaos methods for hyperbolic partial differential equations[END_REF]. In our context, Y (ω) can be explicitely expressed as a functional relationship and therefore the two approaches are of equal complexity. In this paper, projection method based on the Gauss quadrature rules are used to obtain the chaos coefficients. Then statistical moments and Sobol' sensitivity indices can easily be computed from polynomials chaos expansion.

The paper is organized as follows. The next section, is devoted to restate some basic notions about orthogonal polynomial and Chaos expansion. Then, the Sobol' indices which are an important tool in sensitivity analysis are introduced. Epistemic uncertainty in input parameters of two queueing systems with unreliable server is considered in Section 4. In the first queue, a single input parameter is considered as random. On the other hand, the second queueing system considers four input random variables and a sensitivity analysis is performed. For the two queues, numerical results for the computation of statistical moments are obtained using both chaos expansion and Monte-Carlo simulation. Finally, Section 5 draws some conclusions.

Polynomial Chaos Expansion

Polynomial Chaos expansion shares many elements with power series involving orthogonal polynomials, which are commonly used in the spectral community. This is why the following subsection is devoted to restate some results concerning orthogonal polynomials.

Orthogonal Polynomials: univariate case.

Let V be the real vector space of all polynomials in one variable with real coefficients and let be the positive inner product on V defined as

< u, v >= I u(x)v(x)f (x)dx ∀u, v ∈ V
where f : I ⊂ R → R + is a nonnegative integrable function of x hereinafter referred to as a 'weight function'. The set of polynomials {Ψ n } n≥0 are said to be orthogonal with respect to the weight function f (x), if Ψ n is a polynomials of degree n and

< Ψ n , Ψ m >= I Ψ n (x)Ψ m (x)f (x)dx = h 2 n δ n,m , n, m ∈ N, (2.1)
where δ n,m is the Kronecker's delta function, and h n are non-zero constants. We recall that, for orthogonal polynomials of degree d = 0, Ψ 0 is always equal to one (Ψ 0 = 1). Furthermore, the system (2.1) is called orthonormal if h n = 1.

The most general way to build such polynomials rely on the three-term recurrence relation

Ψ n+1 (x) = (x -a n )Ψ n (x) -b n Ψ n-1 (x) Ψ 0 (x) = 1, Ψ -1 (x) = 0 (2.2) with          a n = < xΨ n , Ψ n > < Ψ n , Ψ n > , n ∈ N b n = < Ψ n , Ψ n > < Ψ n-1 , Ψ n-1 > , n ∈ N * . ( 2.3)
(see Gautschi [START_REF] Gautschi | Orthogonal polynomials: computation and approximation[END_REF] for more details). Furthermore, according to Favard's theorem, for a given weight function f , corresponds a unique set of coefficients (a n , b n ) n∈N . From a numerical point of view, the integrals appearing into (2.3) can be evaluated using Fejér quadrature rule, see [START_REF] Rahman | Extended Polynomial Dimensional Decomposition for Arbitrary Probability Distributions[END_REF], for example. Then, taking an orthogonal polynomial basis of degree n, {Ψ 0 , Ψ 1 , . . . , Ψ n }, any sufficiently regular function u : I ⊂ R → R may be represented by its projection Π n u on such basis, i.e.

Π n u(x) = n i=0 u i Ψ i (x), (2.4)
where the coefficients of the projection can be computed by evaluating u i = <u,Ψ i > <Ψ i ,Ψ i > . Depending on the regularity of the function u and the choice of the polynomial basis, upper bounds of the truncation error u -Π n u can be derived [START_REF] Gottlieb | Numerical analysis of spectral methods: theory and applications[END_REF]. In particular, when u is C ∞ , one can observe the so-called spectral convergence (the truncation error decays exponentially fast with respect to n). The construction of multivariate orthogonal polynomial basis rely on univariate polynomial basis as will be described in the following section.

Orthogonal Polynomials: multivariate case.

The setting is now the real vector space V d of all polynomials in d variables with real coefficients, together with the positive inner product on V d defined as

< u, v >= I u(x)v(x)f (x)dx ∀u, v ∈ V d , (2.5)
where, now, f :

I ⊂ R d → R + is a nonnegative integrable function of x =(x 1, x 2, . . . , x d ).
We further assume that the weight function f can be written as the product of univariate weight functions, i.e. f (x 1, x 2, ..,

x d ) = f 1 (x 1 )f 2 (x 2 ) . . . f d (x d )
, where f i :

I i ⊂ R → R + are nonnegative integrable function.
It would be tempting to construct the elements of the multivariate polynomial basis in a similar fashion i.e. as a tensor product of monovariate elements. However, such construction leads to uncomplete basis (all monomials up to a certain degrees are not represented). To ensure that {Ψ 0 (x), Ψ 1 (x), .., Ψ P p d (x)} is a multivariate polynomial basis of degree d with (P p d + 1) elements, the procedure described in the sequel is usually adopted. We first define Ψ α 1 ...α d (x) as the tensor product of the elements of a univariate polynomial basis, of degree n i.e.

Ψ α 1 ...α d (x) = d i=1 Ψ α i (x i ), (2.6)
where α i ∈ {0, 1, . . . , n}. However, not all the elements of the form (2.6) are retained when constructing the multivariate polynomial basis (that would lead to an uncomplete basis of degree nd with (n + 1) d elements). Instead, for a given degree p, only the elements that satisfy

d i=0 α i ≤ p in (2.6
) are kept and a one to one correspondence between the multi-index (α 0 , . . . , α d ) and the i th element Ψ i (x) of the multivariate basis is set. Proceeding that way, it can be shown that the number of elements of a complete multivariate polynomial basis of degree p is

P p d + 1 = p + d d = (p + d)! d!p! . (2.7)
Similarly to (2.4), any function u : R d → R may be represented by its projection Π p u on such basis, i.e.

Π p u(x) = P p d i=0 u i Ψ i (x). (2.8)
Having set the basic framework of orthogonal polynomials, we now look into the way they can be efficiently used in the field of probability, where they are usualy refered as Polynomial Chaos (PC).

Polynomial Chaos expansion.

Denote by (Ω, A, P) the probability space, where as usual Ω is the set of all possible outcomes, A is a σ-algebra over Ω, and P is a function A → [0, 1] that gives a probability measure on A. Consider an R d -valued independent random vector

X = (X 1 , • • • , X d ) that describes input uncertainties.
The probability law of X may be defined by the probability density function

f X (x) = d i=1 f i (x i ), (2.9)
where f i (x i ) is the marginal probability density of X i defined on (Ω i , A i , P i ). Let us now denote L 2 (Ω i , A i , P i ) the space of real random variables with finite second order moments, i.e. such that

E(X 2 i ) = x 2 i f i (x i )dx i < ∞, (2.10)
where E stands for the mathematical expectation. L 2 (Ω i , A i , P i ) is a Hilbert space which can be provided with a set of complete orthogonal basis {Ψ i j (x)} j≥0 that are consistent with the density of X i . For example, Legendre polynomials are associated with uniform distributions whereas Hermite polynomials are associated with Gaussian distributions. Similarly, L 2 (Ω, A, P) is provided with a set of complete multivariate orthogonal basis {Ψ j (x)} j≥0 which, in turn, is consistent with the density of X. A,P) for i = 1, . . . , ℓ be a mathematical model. For the sake of simplicity, we will consider only one component of this model, denoted by Y , since the same procedure apply identically to all the other components. Since Y is assumed to belong to L 2 (Ω, A, P), it can be represented as [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF][START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]

Let Y = (Y 1 , . . . , Y ℓ ) : Ω → R ℓ such that Y i ∈ L 2 (Ω,
Y (X 1 , .., X d ) = z 0 Ψ 0 + ∞ i 1 =1 z i 1 Ψ 1 (X i 1 ) + ∞ i 1 =1 i 1 i 2 =1 z i 1 i 2 Ψ 2 (X i 1 , X i 2 ) + + ∞ i 1 =1 i 1 i 2 =1 i 2 i 3 =1 z i 1 i 2 i 3 Ψ 3 (X i 1 , X i 2 , X i 3 ) + • • • , (2.11)
which, after some rearranging, can be rewritten in a more convenient way as

Y (X) = ∞ j=0 y j Ψ j (X) (2.12)
Then, similarly to equation (2.8), this serie is truncated by keeping terms up to a degree p

Y (X) ≈ Π p Y (X) = P p d j=0 y j Ψ j (X). (2.13) Example 2.1. Let X = (X i 1 , X i 2
) be a random normal vector. The functional approximation of the response Y = Y (X) may be approximated with Hermite-Chaos expansions according to (2.13). If we fix the degree of the Chaos expansions to be d = 2, we obtain the series expansion with P 2 2 = 6 terms as follows:

Y ≈ 5 i=0 y i Ψ i (X 1 , X 2 ) (2.14) = y 0 Ψ 00 + y 1 Ψ 10 + y 2 Ψ 01 + y 3 Ψ 11 + y 4 Ψ 20 + y 5 Ψ 02 ,
where

Ψ i 1 i 2 (X 1 , X 2 ) = Ψ i 1 (X 1 )Ψ i 2 (X 2 )
is the product of univariate Hermite polynomials of degree i 1 and i 2 satisfying i 1 +i 2 ≤ 2. Such polynomials can easily be computed according to the recurrence formula (2.2) together with (2.3) with a n = 0 and b n = n. Should we have used Legendre polynomials instead of Hermite's ones, they would be built by taking a n = 0 and b n = n 2 (4n 2 -1) . Table 1 shows the bivariate Hermite polynomial basis constructed from univariate basis {1, X 1 , X 2 1 -1} and {1, X 2 , X 2 2 -1} and the link between the main index j and the multi-index (i 1 , i 2 ). For completness, the same Table also provides bivariate Legendre polynomial basis constructed from univariate basis {1, X 1 , X 2 1 -1/3} and {1, X 2 , X 2 2 -1/3}. [START_REF] Pettersson | Polynomial Chaos methods for hyperbolic partial differential equations[END_REF]. Here we have used the first one which consists in premultiplying (2.13) by Ψ j (X) and by taking the expectation of the resulting product. Using the orthogonality of the PC basis, most of the terms cancel and we end up with

j (i 1 , i 2 ) Ψ i 1 i 2 (Hermite) Ψ i 1 i 2 (Legendre) 0 (0, 0) 1 1 1 (1, 0) X 1 X 1 2 (0, 1) X 2 X 2 3 (1, 1) X 1 X 2 X 1 X 2 4 (2, 0) X 2 1 -1 X 2 1 -1/3 5 (0, 2) X 2 2 -1 X 2 2 -1/3
y j = E(Y (X)Ψ j (X)) = I⊂R d Y (x)Ψ j (x)f (x)dx for j = 0, 1.., P p d .
(2.15)

The above integral can be evaluated through different techniques: going from rough Monte-Carlo sampling simulation to Gaussian quadrature rule or sparse quadrature rules when the dimension d of the input random is high. Here we evaluate such integrals using Gaussian quadrature rules which take the form

I⊂R d Y (x)f (x)dx ≈ Ng 1 i 1 =1 Ng 2 i 2 =1
...

Ng d i d =1 ω i 1 ω i 2 ..ω i d Y ( x i 1 , x i 2 .., x i d ), (2.16
y 0 = E(Y (X)), (2.17)
i.e. the first coefficient of the PC expansion is the expectation of the random response of the system. Similarly, by considering the approximation of

Y 2 Y (X) 2 ≈ P p d i=0 P p d j=0 y i y j Ψ i (X)Ψ j (X), (2.18)
and taking the expectation on each side, the orthogonality of the PC basis leads to

E(Y (X) 2 ) = P p d i=0 y 2 i , (2.19)
from which the variance of the random response of the system can easily be deduced. Furthermore, the PC decomposition provides a convenient way of computing Sobol' indices, as explained in the next section.

Sensitivity Analysis

The purpose of sensitivity analysis, is to investigate the influence of each input parameter and their possible interactions onto the output measures. They can be casted into two main families: local analysis based on a local perturbation around an average value and global analysis that consider input parameters as random variables and decompose the output variance into several components. The Sobol' indices belong this last type of family.

Sobol' indices.

As in the previous section, we consider the mathematical model Y = Y (X), where the input parameter X = (X 1 , . . . , X d ) are d independant random variables belonging to L 2 (Ω i , A i , P i ) for i = 1, .., d and similarly, Y is assumed to belong to L 2 (Ω, A, P). In 1993, Sobol [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] proposed an indicator of the influence of the input parameter X i defined by

S i = V (E(Y /X i )) V (Y ) = V i V , (3.1)
commonly termed 'first order Sobol' indices'.

V i = V (E(Y /X i )) is the conditionnal variance of Y with respect to X i and V = V (Y )
is the total variance of Y . Similarly, sensitivity indices of higher order can be defined by first introducing the following decomposition of the total variance, which is valid for any output Y belonging to L 2 (Ω, A, P)

V = d i 1 =1 V i + 1≤i 1 <i 2 ≤d V i 1 i 2 + 1≤i 1 <i 2 <i 3 ≤d V i 1 i 2 i 3 + • • • + V i 1 ...i d , (3.2)
where

V i 1 = V (E(Y /X i 1 )) V i 1 i 2 = V (E(Y /X i 1 , X i 2 )) -V i 1 -V i 2 V i 1 i 2 i 3 = V (E(Y /X i 1 , X i 2 , X i 3 )) -V i 1 i 2 -V i 1 i 3 -V i 2 i 3 -V i 1 -V i 2 -V i 3 • • • (3.3) V i 1 ...i d = V - d i=1 V i - 1≤i 1 <i 2 ≤d V i 1 i 2 -• • • - 1≤i 1 <i 2 ...<i d-1 ≤d V i 1 ...i d-1
Then, Sobol' indices of order k are given by

S i 1 ...i k = V i 1 ...i k V . (3.4)
Although Sobol' indices could be computed by estimating integrals apearing into equations (3.3), such a procedure would be both computationaly expensive and hardly tractable. Instead, computing those indices from the PC representation of the random output Y turns out to be an efficient alternative and has been the method of choice for many years since the seminal work of Sudret [START_REF] Sudret | Global sensitivity analysis using polynomial Chaos expansion[END_REF]. For that, we define by

I i 1 ,i 2 ,...,is (s ≤ d) the set of d-dimensional vectors α = (α 1 , . . . , α d )
with α = 0 and 0 ≤ α 1 ≤ . . . ≤ α d ≤ d that selects elements of the PC basis Ψ α 1 ...α d defined by (2.6) containing solely the variables X i 1 , X i 2 , . . . , X is . This way, the multi-indices defined by I i will select elements of the PC basis depending only on the variable X i . Similarly, I i,j will select elements of the PC basis that depend on X i and X j , to the exclusion of any other variable. With this notation, Sobol's indices of order k are simply given as a function of the PC coefficients as follows:

S i 1 ...i k = 1 V (α 1 ,...,α d )∈I i 1 ,...,i k y 2 α 1 ,...,α d . (3.5)
For a problem with d input random parameters, it can be shown that (2 d -1) Sobol' indices can be computed for each output random quantity of interest. When the number of input r.v. is large, the number of Sobol' indices grows exponentially and it becomes difficult to draw information from these statistics. This is why, in 1996, Homma and Saltelli [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF] introduced the total sensitivity indices S T i (i = 1, . . . , d) which measures the total effect of the i th random input parameter. It is defined as the sum of all sensitivity indices S i 1 ...i k (k = 1, . . . , d) for which, one of the indices i 1 , i 2 , . . . , i k is equal to i, i.e.

S T i = d k=1 (i 1 ,...,i k )∈J k i S i 1 ,...,i k ,
where

J k i is the set of k-dimensionel vectors (i 1 , . . . , i k ) with 1 ≤ i 1 < . . . < i k ≤ d
, such that one of its components is equal to i. We can also compute these indices from the PC coefficients as follows:

S T i = 1 - 1 V (α 1 ,...,α d )∈I i y 2 α 1 ,...,α d ,
where I i is the complementary set of I i .

Example 3.1.

The computation of the sensitivity indices for the example 2.1 where, p = 2 and d = 2 gives

S 1 = 1 V (y 2 1,0 + y 2 2,0 ), S 2 = 1 V (y 2 0,1 + y 2 0,2 ), (3.6)
S 1,2 = 1 V y 2 1,1 , with V = y 2 1,0 + y 2 0,1 + y 2 1,1 + y 2 2,0 + y 2 0,2 (3.7)
In accordance with (3.2), we see that the variance V of the output can be divided into three parts:

V = S 1 + S 2 + S 1,2 .
The first part S 1 represents the influence of the first input random variable alone; the second S 2 represents the effect of the second one alone whereas S 1,2 accounts for the combined effect of the two input random variables.

Similarly, according to the definition of the total sensitivity indices, we have

S T 1 = S 1 + S 1,2 , and
S T 2 = S 2 + S 1,2 .

Uncertainty analysis in unreliable queueing models

In this section, we consider the PC expansion for propagating the uncertainty in performance measures of queueing models with breakdowns and repairs, due to epistemic uncertainties in the model input parameters. Section 4.1 studies the M/G/1/N queue with breakdowns, when the perturbation of a single parameter is introduced. Then, section 4.2 addresses the case of the M/M/1/N queue with breakdowns, for multiple random input parameters.

4.1.

The M/G/1/N queue with breakdowns and repairs. Consider a finite capacity M/G/1/N queue with server subject to breakdowns and repairs. Assume that the customers arrive at the system according to a Poisson stream with rate λ, the service times are general and independent identically distributed with mean 1/µ and we denote the service time distribution by F (x). There can be at most N customers in the queue (including the one being served), and customers attempting to enter the queue when there are already N customers present are lost. The service discipline is assumed to be FCFS. At the beginning of each service, there is a probability θ that the server breaks down (and the customer is sent back to the queue) and enters a repair state. The time of repair is exponentially distributed with rate r. The only instant of time when a server can breaks down is right at the beginning of a service. This system is modeled by a Markov chain X = {X n : n ∈ N}, embedded at service completions and completion of a repair and has state space S = {0, 1, . . . , N -1}. Let Ξ = (a i,j ) i,j∈S be the transition probability matrix of the Markov chain X. The matrix Ξ has the following form:

Ξ =                     a 0 a 1 a 2 a 3 • • • a N -2 1 - N -2 k=0 a k a 0 a 1 a 2 a 3 • • • a N -2 1 - N -2 k=0 a k 0 a 0 a 1 a 2 • • • a N -3 1 - N -3 k=0 a k 0 0 a 0 a 1 • • • a N -4 1 - N -4 k=0 a k . . . . . . . . . . . . 0 0 0 0 • • • a 0 1 -a 0                     , (4.1)
where

a k = θ ∞ 0 e -λx (λx) k k! dF (x) + (1 -θ) r r + λ λ λ + r k , k = 0, . . . , N -2. (4.2)
Note that the Markov chain X is unichain and therefore the stationary distribution exists. Let π denote the stationary distribution of the queue-length process embedded at service completions and completion of a repair in the M/G/1/N queue with breakdowns and repairs (see [START_REF] Abbas | A critical account of perturbation analysis of Markov chains. Markov Process[END_REF] for details). In the sequel, we consider π as a function of the breakdown probability θ, denoted by π(θ). Assume that the probability θ is estimated from insufficient statistical data, and hence has uncertainty associated with it. In the next section, we will discuss a functional approach based on PC expansion for computing the epistemic uncertainty in stationary distribution π(θ), due to epistemic uncertainties in θ.

Epistemic uncertainty in queueing system.

In this model, we will assume that the uncertainty only affects one parameter. To analyze the propagation of epistemic uncertainty, the mean and the variance of the responses will be computed using the PC expansion. Let us consider π(θ) = (π 0 (θ), • • • , π N -1 (θ)) be the random response or the model output of the queueing system. We consider the following perturbation model on the uncertain parameter:

θ(ω) = θ + σ θ ε(ω), (4.3)
where θ(ω) is a random variable with uniform density on the interval [0, 1]. The parameters θ and σ θ represent the mean and the standard deviation of θ and they can be estimated by statistical methods. ε(ω) can be considered as a random noise inflicted on θ and it modelizes epistemic uncertainty. We assume that ε(ω) is a random variable with uniform distribution on [-1, 1]. Then, the approximation of any component π ℓ (θ) of the stationnary distribution by the PC expansions writes

π ℓ (θ) ≈ P p d i=0 y i Ψ i (θ),
where {Ψ i (θ)} 0≤i≤P p d form an orthonormal Legendre polynomial basis, and y i are the coefficients of the approximating series expansion, see (2.13) .

Numerical results.

In this section, the parameters N , λ and r are deterministic and their values are set to N = 7, λ = 1, r = 0.4, respectively. The only input random parameter is θ and it is given by θ = 0.5 + 0.28 ε, ε ∼ U (-1, 1).

For numerical computations, we propose two types of distribution for the service time, namely, Erlang (E 2 ) and Hyperexponential of second order (H 2 ). Those two distributions are described in the following subsections. 4.1.3. Epistemic uncertainty in the M/E 2 /1 queue with server breakdowns and repairs. In (4.2), the density F ′ (x) = f (x) of the service time is assumed to be Erlang of second order, i.e.

f (x) = µ 1 µ 2 µ 2 -µ 1 (e -µ 1 x -e -µ 2 x )1 [0,+∞) (x),
where the rates µ 1 = 4 and µ 2 = 2.

We determine the stationary distribution vector with respect to the random variable θ(ω) by solving the system π θ Ξ θ = π θ and π θ = 1, where Ξ θ denotes the transition probability given by (4.1). Therefore, the random outputs of interest are the stationary distributions π i , i ∈ {0, . . . , 6}. Their projection is performed on a monovariate PC basis of degree n = 4 and the coefficients are computed using a 6-points Gaussian quadrature rule as explained in section 2.4. Numerical results for the expectations and the variances using a Monte-Carlo simulations of a sample size N M C = 1000 and N M C = 100000 are given in Table 2 and Table 3, respectively. Results are compared with those of the PC expansion. We note that the Monte Carlo simulations (MC in short) converge to the PC results; however the later are obtained at a fraction of the cost of the MC simulations. Similarly moments of higher order (skewness and kurtosis) are given in Tables 4 and5, and the same conclusions can be drawn. Figure 1 show the density of the output stationary distributions for a random uniform input of the form (4.3). The shape of the densities show that they are far from being uniform. Furthermore, the plots are coherent with Tables 4 and5. Indeed, a positive skewness reflects a random variable with higher probabilities in the left part of its support, and vice versa. In this section, the density F ′ (x) = f (x) of the service time is now assumed to be Hyperexponential of second order i.e.

E(π i ) PC MC (N M C = 1000) MC (N M C = 100000) π 0 0.
V (π i ) PC×10 -4 MC×10 -4 (N M C = 1000) MC×10 -4 (N M C = 100000) π 0 1.
f (x) = (γµ 1 e -µ 1 x + (1 -γ)µ 2 e -µ 2 x )1 [0,∞) (x),
where the rates are µ 1 = 3/2, µ 2 = 3; and γ = 0.3. All the other parameters are set as in the previous section, including the random input. Similarly, Tables 6,7, 8 and 9 show the moments of the stationary distributions when the Hyperexponential law is used instead of the Erlang one for the service time. Here also, we can observe the good convergence properties of the PC, showing the robustness of the method independently of the law of the service time.

4.2.

The M/M/1 queue with server breakdowns and threshold-based recovery policy.

Consider a finite capacity M/M/1 queue with server breakdowns and threshold-based recovery policy. If a server provides service for a customer, the server may experience a breakdown, and the repair begins when the number of customers present in the system exceeds some prespecified threshold level q (1 ≤ q ≤ N ). The stream customers arrive at the queue according to a Poisson law of parameter λ. Arriving customers form a single waiting line based on the order of their arrivals Mariginal probability density f πi of of the steady-state vector in M/H 2 /1/7 queue and the server can serve only one customer at a time. There can be at most N customers present at the queue (including the one in service), and customers attempting to enter the queue when N customers are already present are lost. A single server with exponential distributed service times with rate µ is considered. We assume that the server can break down only if the system is not empty. The lengths of breakdowns are identically distributed and follow an exponential distribution with rate α. The duration of reparation of the server is assumed to follow an exponential law with parameter β. After the server being repaired, it switches to working state and continues to provide service until the system becomes empty. The arrival flow customers, service times, breakdown times and repair times are assumed to be mutually independent input random variables. The service discipline is FCFS. A typical state of this system may be denoted by X(t) = {(N (t), Y (t)); t ≥ 0} , where N (t) is the number of customers in the queue at time t, and Y (t) is a random variable representing the server state at time t. If at time t, the server is experiencing a breakdown, then Y (t) = 1; otherwise, the server is in the working state and Y (t) = 0. The stochastic process X(t) is a continuous time Markov chain whose state space S = {(n, 0) : n = 0, . . . , N } ∪ {(n, 1) : n = 1, . . . , N }. The infinitesimal generator matrix Q of the continuous time Markov chain X(t) has the following block-tridiagonal form:

(θ) in M/E 2 /1/7 queue E(π i ) PC MC (N M C = 1000) MC (N M C = 100000) π 0 0.
V (π i ) PC ×10 -4 MC×10 -4 (N M C = 1000) MC×10 -4 (N M C = 100000) π 0 4.
Skew(π i ) PC MC (N M C = 1000) MC (N M C = 100000) π 0 0.
Q =                A 0 B 0 0 0 0 0 0 0 0 . . . D 1 B 1 F 1 0 0 0 0 0 0 . . . 0 D 2 B 2 F 2 0 0 0 0 0 . . . 0 0 D 3 B 3 F 3 0 0 0 0 .
0 0 0 0 0 0 0 . . . D N E N                , (4.4)
where

B 0 = λ, A 0 = -λ, D 1 = µ 0 , B i = -(λ + µ + α) α 0 -λ , i = 1, • • • q -1, F j = λ 0 0 λ , j = 1, • • • , N -1, D k = µ 0 0 0 , k = 2, • • • N M s = -(λ + α + µ) α β -(β + λ) , s = q, • • • N -1, E N = -(α + µ) α β -β .
The Markov transition diagram of the process is given below can be performed via a recursive scheme [START_REF] Yang | Cost analysis of a finite capacity queue with server breakdowns and threshold-based recovery policy[END_REF] , or by using matrix analytic methods [START_REF] Neuts | Matrix-geometric solutions in stochastic models[END_REF]. In the sequel, we follow a different method, and we will provide a functional approach based on the PC expansion for obtaining the stationary distribution π m,n in terms of some parameters, which are obtained under the epistemic uncertainties.

4.2.1. Epistemic uncertainty in queueing system.

In this model, we consider the epistemic uncertainty in all the input parameters of the queueing system, λ, µ, α, β. However, now π m,n (θ) is assumed to be the random output response of the model of the queueing system, where θ = [θ 1 = λ, θ 2 = µ, θ 3 = α, θ 4 = β] is the vector that describes the model parameters. We consider an epistemic uncertainty in the input parameter θ of known probability density function. Due to the uncertainty in input parameters, the responses π m,n (θ(ω)) are considered as random variables. The uncertainty analysis of the functional π m,n may be quantified by the computation of its mean and variance and moments of higher order.

Furthermore, we proceed to a sensitivity analysis in order to determine the most sensitive parameters with respect to the performance measures. A convenient way to compute those quantities is to approximate the functional π m,n (θ(ω)) in the form of the PC expansion. The sensitivity analysis introduced in section 3 help us to point out which parameters generate a higher sensitivity with respect to the stationary distribution. We consider the following perturbation models for all the uncertain parameters:

θ j (ω) = θj + σ j ε j (ω), j = 1, • • • , 4
where θ j (ω) is random variable, θj , σ j , are respectively, estimated mean and the standard deviation of the random variable θ j obtained by statistical methods. ε j (ω) is a random noise inflicted on θ i which models epistemic uncertainty. Here, ε j (ω) is assumed to follow a standard normal distribution

. Let θ = [θ 1 (ω) = λ, θ 2 (ω) = µ, θ 3 (ω) = α, θ 4 (ω) = β], π m,n (θ(ω)) and ε = [ǫ 1 , . . . , ǫ 4 ]
be random fields.

The approximation of the measure π m,n (θ(ω)) in the form of multivariate PC expansions is given by:

π m,n (θ) ≈ P p 4 i=0 y i Ψ i (θ),
where Ψ i (θ) is the multivariate orthonormal Hermite polynomial and y i are the coefficients of the approximating series expansion, (2.13).

Numerical application.

As an illustrative example, the queueing system described in section 4.2 is studied. The capacity is set to N = 5 and the threshold is equal to q = 3. In order to study the sensitivity of the queue, we assume that the exact value of all the input parameters, α, β,λ, µ, are not precisely known and we propose the following perturbation on these parameters:

α = α + σ α ε 1 , ε 1 ∼ N (0, 1) β = β + σ β ε 2 , ε 2 ∼ N (0, 1) λ = λ + σ λ ε 3 , ε 3 ∼ N (0, 1) µ = μ + σ µ ε 4 , ε 4 ∼ N (0, 1)
so α, β, λ, µ, are Gaussian random variables. The mean and the standard deviation of these random variables are set respectively to μ = 7.3, λ = 2, α = 3, β = 4, σ α = 0.04, σ β = 0.02, σ λ = 0.04 and σ µ = 0.02.

Table 10 gives the first order Sobol' indices for all the components of the stationary distribution vector. From this table, we see that for all but π 0,1 and π 1,3 , the most influencial input random parameter is λ. efffect of two input random parameter is negligible in comparison with the effect of a single parameter. The same is also true for higher order indices since the total Sobol' indices are almost equal to first order ones, as shown by Table 12. Note that for the components π 0,1 and π 1,3 , it is the parameter α which is most influencial.

When considering π 0,1 and π 1,3 , we can set α as being the only random parameters and when considering all the other components of the stationary distribution, we can consider the parameter λ as sole input random variable. This is what we do as a second numerical experiments. Then, Tables 13 to 16 compare the statistical moments of the stationary distribution vector when the four input parameters are random (first column) to the case where only the most influencial parameter is random (second column). As expected, the two columns of the tables match rather well.

Conclusion

In this paper, we have developped a numerical approach based on polynomial chaos expansion, to study the sensitivity and the propagation of the epistemic uncertainty in queueing models that occurs with unreliable servers. models of queueing systems have been investigated. In the first model (M/G/1/N queue with breakdowns and repairs), the epistemic uncertainty only affects one input parameter, whereas in the second model (M/M/1/N queue with server breakdowns and threshold-based recovery policy), it affects four parameters. In the latter case, a sensitivity analysis using Sobol' indices is performed. When considering the stationary distribution as output quantity of interest, it was shown that the parameter λ of the Poisson law modelling the customers arrival at the queue, was the most influencial factor. This finding was then confirmed by considering λ as the only random input parameter and by setting the three remaining ones to their average values. In that case, it was shown that the statistical moments of the output measure were relatively insensitive to the other parameters. Finally, comparisons with Monte-Carlo simulations showed the good convergence properties of the chaos expansion.

7 Figure 1 .

 71 Figure 1. Mariginal probability density f πi of the steady-state vector π(θ) in M/E 2 /1/7 queue

7 Figure 2 .

 72 Figure 2. Mariginal probability density f πi of of the steady-state vector in M/H 2 /1/7 queue

Figure 3 .

 3 Figure 3. Markov transition diagram of the M/M/1/N queue with server breakdowns and threshold-based recovery policy

Table 1 .

 1 Bivariate Hermite and Legendre Polynomial Basis 2.4. Computing PC coefficients. The way to compute the PC coefficients {y j } 0≤j≤P p d appearing into equation (2.13) can be casted into two different families: projection methods and regression methods

Table 2 .

 2 Expected value of the steady-state vector in M/E 2 /1/7 queue

		0237	0.0230	0.0237
	π 1	0.0413	0.0422	0.0413
	π 2	0.0660	0.0667	0.0660
	π 3	0.1018	0.1026	0.1018
	π 4	0.1557	0.1563	0.1558
	π 5	0.2393	0.2401	0.2392
	π 6	0.3722	0.3719	0.3722

Table 3 .

 3 Variance of the steady-state vector in M/E 2 /1/7 queue

		7753		1.6686	1.7809
	π 1	3.6413		3.6789	3.6423
	π 2	5.2349		5.1333	5.2175
	π 3	5.0800		4.9922	5.0508
	π 4	1.8737		1.8585	1.8732
	π 5	2.1440		2.0101	2.1534
	π 6	58.7699		61.0693	58.7199
		Skew(π i )	PC	MC (N M C = 1000) MC (N M C = 100000)
		π 0	0.5053	0.5429	0.5064
		π 1	0.3565	0.2911	0.3550
		π 2	0.1730	0.1435	0.1701
		π 3	-0.0722	-0.0927	-0.0730
		π 4	-0.5848	-0.6616	-0.5922
		π 5	-0.7140	-0.7968	-0.7036
		π 6	0.0543	0.0692	0.0540

Table 4 .

 4 Skewness coefficient of the steady-state vector in M/E 2 /1/7 queue

	Kurt(π i )	PC	MC (N M C = 1000) MC (N M C = 100000)
	π 0	2.0384	2.1425	2.0341
	π 1	1.8968	1.9122	1.8953
	π 2	1.7864	1.7682	1.7900
	π 3	1.7463	1.7513	1.7495
	π 4	2.0301	2.0290	2.0385
	π 5	2.1921	2.3603	2.1972
	π 6	1.7759	1.8007	1.7724

Table 5 .

 5 Kurtosis coefficient of the steady-state vector in M/E

2 /1/7 queue 4.1.4. Epistemic uncertainty in the M/H 2 /1 queue with server breakdowns and repairs.

Table 6 .

 6 Expected value of the steady-state vector in M/H 2 /1/7 queue

		0397	0.0407	0.0396
	π 1	0.0602	0.0588	0.0602
	π 2	0.0853	0.0847	0.0853
	π 3	0.1179	0.1179	0.1179
	π 4	0.1619	0.1615	0.1619
	π 5	0.2233	0.2238	0.2233
	π 6	0.3118	0.3103	0.3117

Table 7 .

 7 Variance of the steady-state vector in M/H 2 /1/7 queue

		7367	4.5451	4.7370
	π 1	6.7973	6.5517	6.8120
	π 2	7.1124	7.1342	7.1237
	π 3	4.7642	4.5830	4.7726
	π 4	0.7665	0.8177	0.7653
	π 5	5.4708	5.2979	5.4699
	π 6	64.0300	62.3300	64.0396

Table 8 .

 8 Skewness coefficient of the steady-state vector in M/H 2 /1/7 queue

		3481	0.2807	0.3531
	π 1	0.1478	0.2379	0.1474
	π 2	-0.0881	-0.0925	-0.0896
	π 3	-0.4210	-0.4300	-0.4230
	π 4	-1.2389	-1.1806	-1.2381
	π 5	-0.3119	-0.3561	-0.3101
	π 6	0.2751	0.3242	0.2747
	Kurt(π i )	PC	MC (N M C = 1000) MC (N M C = 100000)
	π 0	1.8643	1.9138	1.8624
	π 1	1.7524	1.7200	1.7568
	π 2	1.7297	1.6881	1.7266
	π 3	1.8885	1.8168	1.8935
	π 4	3.3357	3.7223	3.3412
	π 5	1.7598	1.7821	1.7586
	π 6	1.8392	1.8446	1.8372

Table 9 .

 9 Kurtosis coefficient of the steady-state vector in M/H 2 /1/7 queue

Table 11 ,

 11 which gives the second order Sobol' indices, shows that the combined

	π m,n \S θ i	S α	S µ	S β	S λ
	π 0,0	0.3799 0.0366 0.0046 0.5789
	π 0,1	0.7942 0.0047 0.0096 0.1913
	π 1,1	0.0144 0.0038 0.0078 0.9740
	π 0,2	0.0632 0.0160 0.0080 0.9126
	π 1,2	0.1650 0.0003 0.0109 0.8236
	π 0,3	0.0445 0.0473 0.0073 0.9008
	π 1,3	0.7812 0.0151 0.0418 0.1619
	π 0,4	0.0502 0.0264 0.0006 0.9229
	π 1,4	0.3120 0.0124 0.0411 0.6345
	π 0,5	0.0528 0.0198 0.0044 0.9229
	π 1,5	0.1512 0.0085 0.0396 0.8005

Table 10 .

 10 First order Sobol' indices for stationary distribution vector in M/M/1/5 queue with server breakdowns and threshold-based recovery policy π m,n \S θ i ,θ j S α,µ × 10 -8 S α,β × 10 -8 S α,λ × 10 -6 S µ,β × 10 -9 S µ,λ × 10 -8 S β,λ × 10 -7

	π 0,0	0.6972	1.3636	0.2051	1.8411	4.9130	1.1427
	π 0,1	0.0144	2.8434	0.0127	0.4145	0.0016	0.6002
	π 1,1	0.0361	0.0175	0.5309	0.3338	0.0569	1.9217
	π 0,2	0.0516	0.7513	0.0183	0.0181	0.0249	0.8644
	π 1,2	0.0263	0.0365	0.0143	0.0109	0.0013	7.3456
	π 0,3	0.8840	0.0198	0.3924	0.0283	0.4031	0.8201
	π 1,3	9.8593	1.6277	0.3311	0.1568	0.0997	0.0147
	π 0,4	0.1210	7.1050	4.0803	4.9641	0.0012	8.0038
	π 1,4	0.0137	0.0368	0.4395	0.9306	0.6662	8.6512
	π 0,5	0.0113	0.2362	0.5077	0.6212	0.0029	0.1801
	π 1,5	0.0250	0.0959	0.0187	0.0433	0.0010	0.0427

Table 11 .

 11 Second order Sobol' indices for stationary distribution vector in M/M/1/5

	queue with server breakdowns and threshold-based recovery policy
	π m,n \S T θ i	S T α	S T µ	S T β	S T λ
	π 0,0	0.3799 0.0366 0.0046 0.5789
	π 0,1	0.7943 0.0047 0.0096 0.1914
	π 1,1	0.0144 0.0038 0.0078 0.9741
	π 0,2	0.0634 0.0160 0.0080 0.9128
	π 1,2	0.1652 0.0002 0.0109 0.8237
	π 0,3	0.0445 0.0473 0.0074 0.9009
	π 1,3	0.7812 0.0151 0.0419 0.1619
	π 0,4	0.0502 0.0264 0.0005 0.9229
	π 1,4	0.3121 0.0124 0.0411 0.6345
	π 0,5	0.0528 0.0198 0.0044 0.9230
	π 1,5	0.1514 0.0085 0.0396 0.8007

Table 12 .

 12 Total order Sobol' indices for steady-state vector in M/M/1/5 queue with server breakdowns and threshold-based recovery policy

Table 13 .

 13 Expected Value of the stationary distribution in M/M/1/5 queue with server breakdowns and threshold-based recovery policy V (π m,n ) PC ×10 -6 (1 r.v.) PC×10 -6 (4 r.v.)

	E(π m,n ) PC (1 r.v.) PC (4 r.v.)
	π 0,0	0.2639	0.2639
	π 0,1	0.0723	0.0723
	π 1,1	0.1084	0.1084
	π 0,2	0.0495	0.0495
	π 1,2	0.1827	0.1827
	π 0,3	0.0636	0.0636
	π 1,3	0.0927	0.0927
	π 0,4	0.0428	0.0428
	π 1,4	0.0523	0.0523
	π 0,5	0.0261	0.0261
	π 1,5	0.0457	0.0457
	π 0,0	14.711	25.338
	π 0,1	0.7178	0.8879
	π 1,1	2.4845	2.5455
	π 0,2	0.4445	0.4905
	π 1,2	4.2497	5.0989
	π 0,3	0.7925	0.8793
	π 1,3	0.6677	0.8450
	π 0,4	1.4357	1.5550
	π 1,4	0.9757	1.5432
	π 0,5	1.2546	1.3592
	π 1,5	3.4476	4.3111

Table 14 .

 14 Variance of the stationary distribution in M/M/1/5 queue with server breakdowns and threshold-based recovery policy

Table 15 .

 15 To illustrate the applicability of the proposed approach, two Skew(π m,n ) PC (1 r.v.) PC (4 r.v.) Skewness coefficient of the steady-state vector in M/M/1/5 queue with server breakdowns and threshold-based recovery policy Kurt(π m,n ) PC (1 r.v.) PC (4 r.v.)

	π 0,0	0.0350	0.0395
	π 0,1	0.0443	0.0415
	π 1,1	0.0350	0.0371
	π 0,2	-0.0722	-0.0411
	π 1,2	0.0200	0.0362
	π 0,3	-0.0536	-0.0523
	π 1,3	-0.0316	-0.0556
	π 0,4	-0.0011	-0.0042
	π 1,4	-0.0622	-0.0116
	π 0,5	0.0488	0.0576
	π 1,5	0.0374	0.0668
	π 0,0	3.0016	3.0020
	π 0,1	3.0037	3.0046
	π 1,1	3.0016	3.0024
	π 0,2	3.0071	3.0064
	π 1,2	2.9999	3.0023
	π 0,3	3.0041	3.0035
	π 1,3	3.0017	3.0060
	π 0,4	2.9983	2.9979
	π 1,4	3.0044	2.9988
	π 0,5	3.0013	3.0028
	π 1,5	2.9995	3.0052

Table 16 .

 16 Kurtosis coefficient of the steady-state vector in M/M/1/5 queue with server breakdowns and threshold-based recovery policy