Addendum-Erratum to Dissipative Systems Analysis and Control, Springer-Verlag London, 2nd Edition, 2007 (the first figure indicates the page number)

• Notation section: λ(A): an eigenvalue of A ∈ IR n×n .

• Notation section: σ(A): the set of eigenvalues of A ∈ IR n×n (i.e. the spectrum of A).

• 17, 18, 28: line 2 after (2.22), and also in (2.24), (2.29), (2.30), line 3 in the equation above (2.75): change u T (t) to u t (τ ), and u T (jω) to u t (jω).

• 17: line 8, the integration is done over τ ∈ [0, t]

• 32: obviously the expression for P (ω e ) under (2.95) is meaningless. Please contact Prof. O. Egeland for details.

• 59, lemma 2.55: rational

• 81, Lemma 3.11: µ min (A) denotes the minimum eigenvalue of the symmetrized matrix 1 2 (A + A T ). • 81: σ min (B) > 0 is equivalent to Ker(B) = {0}, and to m < n and B of full rank m.

• 81, last line: it is H(jωµ), not T (jωµ)

• 89 line 6: it is H(jω), not T (jω)

• 88, line -3: it is H(jω), not T (jω)

• 101: in the framed set of implications/equivalences, the last implication has to be reversed and in the last line this is P = P T ≥ 0 (hence the last condition is really sufficient for PRness as indicated in the paragraph in between the two framed sets).

Reference [444] is also worth looking at.

• 119: in (3.134) it is review, vol.49, no 4, pp.545-592, 2007, and in particular section 6 therein. Their Theorem 6.2 is close to Theorem 2.49 page 57 in the book.

∂V f ∂x (x)[Ax + Bu] + w(x, u) ≥ 0 • 121,
• 147, example 3.75: it is Ax = 0 if x < 0. The notation x + means max(0, x).

• 149, the dissipativity of linear complementarity systems is throughly investigated in Camlibel, Iannelli, Vasca, "Passivity and complementarity", Mathematical Programming A, 2013, DOI: 10.1007/s10107-013-0678-4 . Closely related results concerning Lur'e set-valued systems are in Brogliato and Goeleven, "Existence, uniqueness of solutions and stability of nonsmooth multivalued Lur'e dynamical systems", Journal of Convex Analysis, vol.20, no 3, 2013, and "Well-posedness, stability and invariance results for a class of multivalued Lur'e dynamical systems", Nonlinear Analysis: Theory, Methods and Applications, vol.74, pp.195-212, 2011.

• 167, in Remark 3.96 this is [0, ∞) (see Theorem 3.91)

• 175, preservation of dissipativity after time-discretization has recently been tackled in S. Greenhalg, V. Acary, B. Brogliato, "Preservation of the dissipativity properties of a class of nonsmooth dynamical systems with the (θ, γ)-algorithm", Numerische Mathematik, 125(4), 2013. See several references that concern this issue in this article. The issue that is tackled concerns the preservation of the strorage function, the supply rate, and the dissipation function, and therefore yields much more stringent conditions than the mere preservation of dissipativity (with possible different supply rate, dissipation function after discretization).

• 170, about the Popov's line. In the book by Aiserman and Gantmacher from 1965 it is pointed outy that Popov criterion also holds in the case of a negative slope 1 r . • 183, Fact 8: the inequality is strict in (1 ≤ p < +∞), because obviously this does not work for the infinity norm (take f (•) as a constant non zero function).

• 183: Fact 7: in the proof, we compute the smallest lower-bound.

• 192, in (4.18) all T are t

• 197: second line of definition 4.26: u(t) ∈ U.

• 197-198: the role of the additive constants β which basically accounts for initial stored energy, has also been investigated in D.H. Hill: 'Dissipative nonlinear systems: basic properties and stability analysis", Proc. of IEEE CDC, pp.3259-3264, 16-18 December 1992.

• 204, Theorem 4.43: It is supposed that there exists a storage function V (•) such that V (x ⋆ ) = 0, so that V r (x ⋆ ) = 0. From the definition of the required supply in Definition 4.36, this also assumes that the system is reachable from x ⋆ . This is Theorem 2 in [510].

• 202: more on reversibility and its relationship with reciprocity, for LTI systems, is in [511,sections 8 and 9]. In particular Theorem 8 of [51] provides a way to check reversibility.

• 209 and 210: in (4.45) the last term is y T (s)Su(s), and same mistake in line 4 after Theorem 4.56.

• • 230,Lemma 4.91: in (4.94) this is V (t, x).

• p.247, line 10: alll should be all.

• 248 and 249, replace all (4.108) by (4.163)

• 251, Theorem 4.111, line 2: IR m×m

• Section 5.1: Most of the results presented in this section have been stated in two articles by D.J. Hill and P. Moylan: "Stability results for nonlinear feedback systems", Automatica, vol.13, pp.377-382, 1977, and "General instability results for interconnected systems", SIAM J. Control and Optimization, vol.21, no 2, pp.256-279, 1983. The errata for these papers include the following. In the SIAM paper, the first condition that appears in equation ( 12), Theorem 7, is redundant. In the Automatica article, condition (iii) of Theorem 4 can be removed, while condition (ii) should be stated with T Ā(•) a nonnegative real valued function.

• 254, 255: in (4.185) replace ∂u ∂x in the first term of the left-hand side by ∂u ∂t . Line after lemma 4.115: do the same for the equality after "One has", and delete the "is". Also the power 2 is missing in the integrand. In (4.188) multiply the first term in the left-hand side by 1 2 . These two Lemmas are taken from H. Brézis, Analyse Fonctionnelle, Théorie et Applications, Masson, Paris, 1983, sections X.1 and X.3.

• 267: from Lemma 5.13 and Corollary 3.4 it may be deduced that the solutions P = P T of the KYP Lemma LMI, are > 0 if (C, A) is observable.

• 298: in Theorem 5.68, this is "...with input u(•)..."

• 305: line before section 5.9.4: Bounded

• 306: Theorem 5.71 is taken from M.R. James et al, SIAM J. Control Optimization, vol.43, no 5, pp.1535Optimization, vol.43, no 5, pp. -1582Optimization, vol.43, no 5, pp. , 2005.

• 339: In definition 6.37: a symmetric positive definite matrix R(x); in equation (6.95), this is (J(x) -R(x)) ∂H 0 ∂x + .... • 366: generally speaking BV functions do not satisfy that for any t, there exists a σ such the function is continuous on [t, t + σ). This is because at t a BV function may possess an accumulation of jumps on the left, and an accumulation of jumps on the right. But in mechanical systems with frictionless unilateral constraints and piecewise analytic data, accumulations of impacts on the right do not exist (result of Ballard in 2000). So in this particular case the generalized velocity has this property.

• 368: line -8, it is -∂ψ V (x).

• 382: ξ in the line after (7.10) refers to the ξ in (6.184)

• section 3.9.4: Notice that for a maximal monotone operator F : IR n → IR n , x → y = F (x), with F = ∂f where f (•) is convex, proper, lower semi continuous, then the "input-output" product satisfies x T y = f (x) + f * (y), where f * (•) is the dual function of f (•).

• 409: CTCE means Cross Terms Cancellation Equality

• 431 line 2 section 7.8.3: (6.162) is (6.156), and λ z 1 is as in (6.169) • 533, Lemma A.65: the last condition is S(I -R † R) = 0 or equivalently (I -RR † )S T = 0.

• 536, line 11 in the proof of Lemma A.71: realization

• 564, reference [470]: 1992

  line 5 after remark 3.38: Thus H(-jω, jω) is the condition.... • 123, left hand side of (3.148): P is G • 124, line 2: P is G • 125: line 5 after Theorem 3.44: numerically • Section 3.9 (Lur'e problem and absolute stability): it's also worth reading the survey by R. Shorten et al, "Stability Criteria for Switched and Hybrid Systems", SIAM

  230 and 116: Lemma 3.36 page 116, and Lemma 4.91 page 230, are taken from the book referenced 145 by Faurre, Clerget and Germain. In fact Lemma 3.36 is presented in the book by Faurre et al as a corollary of Lemma 4.91.


Sophie Ponthieux
email: sophie.ponthieux@insee.fr

Marie Cottrell
email: cottrell@univ-paris1.fr

Living conditions: classification of households using the Kohonen algorithm

Living conditions: classification of households using the Kohonen algorithm

Sophie Ponthieux and Marie Cottrell

Abstract:

In the analysis of poverty and social exclusion, indicators of living conditions are some interesting non-monetary complements to the usual measurements in terms of current or annual income. Living conditions depend in fact on longer term factors than income, and provide further information on households' actual resources that allow to compare more accurately between living standards. But in counterpart, a difficulty comes from the qualitative nature of the information, and the large number of dimensions and items that may be taken into account; in other words, living conditions are difficult to "measure". A consequence is that very often, the information is either used only partly, or reduced into a global score of (bad) living conditions, that results from counting "negative" items, and the qualitative dimension is lost. In this paper, we propose to use the Kohonen algorithm first to describe how the elements of living conditions are combined, and secondly to classify households according to their living conditions.

The main interest of a classification is to make appear not only quantitative differences in the "levels" of living conditions, but also qualitative differences within similar "levels".

Introduction

Since the 1970s in United-Kingdom, more recently in France, the use of non-monetary indicators in the analysis of poverty is developing, and poverty or exclusion are studied both in terms of income and in terms of living conditions [START_REF] Townsend | Poverty in the United-Kingdom[END_REF][START_REF] Nolan | Resources, deprivation and poverty[END_REF]. Living conditions include a great number of domains. [START_REF] Dickes | Ressources financières, bien-être subjectif et conditions d'existence[END_REF] lists ten of them: dwelling, durables, food, clothing, financial resources, health, social relations, leisure, education and work. Not all the existing studies include this complete set of domains. The choice of including or not one of those domains may be based on two arguments: in the first one, the main hypothesis is that the subjects are rational in their behavior, which leads to select only the domains where privations are supposed to decrease or disappear when the financial resources increase [START_REF] Mack | Poor Britain[END_REF]; the second one is based on the notion of "standard" [START_REF] Townsend | Deprivation[END_REF][START_REF] Dickes | Ressources financières, bien-être subjectif et conditions d'existence[END_REF], and leads to consider any domain as soon as all the subjects are -at least potentiallyinvolved. Each domain is described by a list of items, that can be, for example, a characteristic of the dwelling ("is there central heating/bath or shower/etc. in the house"?), a particular consumer good ("do you have a color TV/a car/a washing machine/etc."?), the ability to afford ("new clothes/a week's holiday away/etc."?); for each item, the respondent for the household declares to have or have not, to be able to afford or not.

Whatever the option in choosing the domains, the main difficulty is to deal with a great quantity of information (depending on the number of dimensions retained and the number of items within each dimension) that is mainly qualitative. One possibility is to construct a "score" of (bad) living conditions, which corresponds to the total number of "negative" items (for example, [START_REF] Lollivier | Pauvreté d'existence, monétaire ou subjective sont distinctes[END_REF]. The interest of this approach is that high scores will always characterize households who deal with cumulative difficulties; nevertheless it will not allow to distinguish between households who have smaller scores, but difficulties of completely different nature; and since there is no method for weighting the items, not having a TV set has exactly the same value as not being able to afford buying meat/chicken/fish every second day, even though it is clearly not the same nature of hardship.

An alternative to using a score is to classify households into groups that would be consistent both in terms of living conditions "level" and in terms of (bad) living conditions nature. For this purpose, several classification methods can be used; in what follows, we use a classification algorithm relating to self-organizing maps, the Kohonen algorithm. The data are from the third wave of the French part of the European community households panel.

In a first step, we classify only the modalities, in order to obtain a good description of living conditions characteristics: how they are combined, i.e. what are the most frequent associations between the items modalities. The result show some neat oppositions not only between "positive" and "negative" modalities but also between groups of items, that could suggest some kind of qualitative gradation in living conditions hardship. In a second step, we classify the observations, searching for a consistent grouping of households described only by their living conditions. The classifications obtained tend to confirm that beyond differences in the "level" of (bad) living conditions, there are significant within level differences in the nature of the difficulties.

Data, variables of the analysis and method

The data source is the French part of the European Community Households Panel, here in its third wave (year 1996). It provides the required detailed information about material living conditions (dwelling, environment, durable goods, deprivations); we also know the household's income, and whether the household's respondent considers that the monetary resources allow to live from "very comfortably" to "with great difficulty" ("subjective" living conditions after). Finally, the source covers three dimensions of households' poverty: in terms of monetary income, in terms of material living conditions, and in terms of "subjective" living conditions; this will allow to characterize the classes of living conditions from the two other points of view.

The observations are households. For the classifications, we use only their characteristics relating to material living conditions. Living conditions are described by 10 items about the dwelling (5 about convenience and 5 about problems), 4 items about environmental topics, 6 items for the durables and 6 items about deprivations, a total of 26 dummy variables, that is 52 modalities (detailed in Appendix, table A1). For each item, the "negative" modality (having a problem, not having an item, not being able to afford) is always coded "1" vs. "0" for the "neutral" modality.

The observations are also described by a set of general characteristics (Appendix, table A2): type of household, average age of the adults (persons aged 17 years and over), number of children under 17 years, current income per consumption unit, and type and location of the dwelling. In addition, we also characterize households by an indicator of monetary poverty, an indicator of subjective living conditions, and scores of material living conditions (total and partialby domain -). These general characteristics are used only to compare between classes, not to construct them.

Only the observations with no missing variable for all these descriptors are kept for the analysis, that is 6458 households.

For the classifications, the information in input is not always under the same form: we use a response table to classify the characteristics and successively the partial scores (score by domain) then the coordinates (after a multiple correspondence analysis, MCA) to classify the observations. The classifications are done using the Kohonen algorithm (for an introduction to the algorithm and its applications to data analysis, see [START_REF] Kohonen | Self-organization and Associative Memory[END_REF][START_REF] Kohonen | Self-organization and Associative Memory[END_REF][START_REF] Kohonen | Self-Organizing Maps[END_REF][START_REF] Kaski | Data Exploration Using Self-Organizing Maps[END_REF][START_REF] Cottrell | The Kohonen algorithm : a powerful tool for analysing and representing multidimensional quantitative et qualitative data[END_REF]. The main property of the Kohonen algorithm is its property to preserve the topology of the data: after convergence, similar data are grouped into the same class or into neighbor classes. This feature allows to represent the proximity between data, as in a projection, in the Kohonen map. As a further treatment, the Kohonen classes can be clustered into a reduced number of macroclasses (which only contain neighbor Kohonen classes) by using a classical hierarchical classification.

Classif ication of the modalities

In this first step, the objective is to obtain a representation of the combination between the modalities of the whole set of items. We have tested successively a one-dimension and a two-dimensions classifications using a Kohonen algorithm inspired by the MCA (Cottrel, Letremy & Roy, 1994;[START_REF] Cottrell | Multiple correspondence analysis of a crosstabulation matrix using the Kohonen algorithm[END_REF].

The results from the first one are represented along a ten classes string (Kohonen map of one dimension). We obtain seven classes of "negative" modalities, and three classes of "neutral" modalities. On the "negative" side, the first class (reading Figure 1 from the left) associates low standard dwelling (no bath or shower, no hot running water and no indoor toilet) and the absence of a telephone and TV set, which are very common items (owned by about 99% of the households). The next class groups together not having a car and not being able to afford one meal of meat/fish/chicken every second day, that is a very serious deprivation. The third class associates several deprivations and the absence of some "modern" durable goods (no micro-wave oven, no VCR -but it is to be related to the absence of a TV set -, and no dish-washer). The three next classes are characterized by the association of problems relating to the quality of he dwelling (dampness, shortage of space) and environmental problems (pollution, noise), the last class of this group adding the inability to afford one week's holiday away from home. The last of the "negative" classes contains only one characteristic which is the inability to replace broken or worn furniture. The interest of this classification is that it indicates a qualitative gradation in the seriousness of the conditions; a, study of the code-vectors profiles shows also a neat gradation of the negative modalities, which is an interesting result in that it could be usable as a basis for a "weighting" of the items.

In a second classification, the position of the modalities are represented on a 10 x 10 grid (Figure 2.a). At first glance, the map shows two regions, with the « neutral » modalities grouped at the top and on the right of the grid, and the « negative » ones at the bottom and on the left. This division is confirmed when we look at the representation of the distances between classes (Figure 2.b). In the bottom left cell we find again the characteristics of a poor standard of the dwelling, and generally, the bottom line corresponds to attributes of rather bad living conditions. The seriousness of the characteristics tends to decrease when we go towards the top of the map, but it is also interesting to notice that there are many modalities that have no immediate neighbor, suggesting that the hypothesis of cumulative difficulties could be somewhat reductive. If we group now the classes according to their closeness on the map, we obtain 3 groups of "negative" modalities: a first one grouping the characteristics of very serious living conditions (low standard dwelling, absence of very common durable goods, and privations in elementary consumptions), a second one corresponding mainly to other problems relating to the dwelling and environmental disadvantages, and a last group, at the "frontier" between "negative" and "neutral" conditions, characterized by the inability to afford one week's holidays or replacing worn out furniture that suggest a particular status for these items. These results are consistent with those obtained with a MCA (Figure 3).

Classes of households

We try now to classify the households according to their living conditions characteristics.

First, in order to use the full qualitative resource of the initial information, we use in input the coordinates of the observations after a MCA (this corresponds to a transformation of the responses into quantitative values). Then we use only the partial scores calculated by domain of the living conditions, each household being then characterized by 5 scores.

Classification using the whole set of items

The inputs are now the coordinates of the observations, resulting from a "traditional" MCA. The Kohonen algorithm is used to classify the observations in a 8 X 8 grid (since we have about 6000 observations, it could give about 100 observations by class). The 64 classes obtained are then grouped into 5 super-classes (SC1 to SC5), using a hierarchical classification. Figure 4 shows the Kohonen map obtained, and indicates the superclasses. For the analysis, another clustering of 3 super-classes -which results in keeping SC1 as one group A, adding SC2 and SC3 into another group B, and SC4 and SC5 into a third group C -is also used.

The analysis of the classes characteristics (table 1) shows some neat differences between group A and the two others: as for the characteristics of material living conditions only, group A appears to be in the best situation; this is also the case when we look at their monetary poverty rate, which is lower than on average, and their "subjective" living conditions, that appear easier than on average. But there are also significant differences between groups B and C: the households in group B are mainly disadvantaged in the domains of durables and privations, while those in group C are mostly concerned with low standard dwellings. The households in group B are also those facing the highest proportion of environmental problems, which is consistent with the proportion living in large structures and big cities; it is also the group where we find the highest proportion of lone parents households. If we look at their other characteristics, the monetary poverty rate is higher in group B than in group C, and "subjective" living conditions more often declared "very difficult". Going back to the five super-classes for more detail, we note that in fact, there are strong differences within groups B and C. Within group B, SC2 appears to suffer from very serious deprivations (in food and clothes), while group SC3 is mostly concerned by a lack of durables. Within group C, we have two sub-groups living in low standard dwellings, but there is a neat difference in the nature of missing items: mostly heating system and absence of a separate kitchen in SC4, and mostly basic comfort (bath or shower, hot running water and indoor toilets) in SC5; SC5 is also the class where we find the smallest average income and the highest proportion of poor. SC5, which counts for 2% of the households, is also characterized by a high proportion of households counting only one person, older than on average and living in rural areas.

These results tend to confirm that differentiations are not only in the level, but also in the nature of the disadvantages.

Classification using the partial scores

Our purpose now is to use a quantitative measurement of living conditions, but also to introduce some of the qualitative dimension of the information. We start from the simplest way to count (bad)living conditions "level", which consists in calculating a score of "bad points" for the whole set of items and setting one or several thresholds. The problem is then to determine where to put the thresholds to obtain consistent classes. This question is discussed in [START_REF] Lollivier | Pauvreté d'existence, monétaire ou subjective sont distinctes[END_REF]; they propose to use the rate of monetary poverty to set a "poverty line" for living conditions: the score that separates "good" and "bad" living conditions is the one for which the cumulative percent of the distribution of the score is equal (or close) to the monetary poverty rate. Applied to our data, the monetary poverty rate is 10,7% and at this cumulative percent of the score's distribution, the score equals 9 (table 2); so we'll say that living conditions are "good" while the score is under 9, and "bad" from a score of 9. One problem with this method is that it cuts the population into one group having "good" living conditions and one group having "bad" living conditions whatever the items; for example, the total value of living in a dwelling without hot running water and not being able to afford a meal of meat/fish/chicken every second day and not being able to afford buying new clothes is the same as having noisy neighbors and not having a VCR and not having a micro-wave oven. So we have tried in what follows to obtain classes defined using the qualitative dimension of the information and independently from an exogenous threshold.

For this, we have classified the observations according to their partial scores (one score by domain) of living conditions. We have used the Kohonen algorithm to classify the households along a one-dimension map of five classes (Figure 5). The "progression" appears neatly going from best to worse living conditions, but not only in terms of level: each of the 5 domains of living conditions appear to contribute more or less to the households positions (table 3). For example, class 3 is mostly characterized by environmental problems, this being also what makes the difference between classes 1 and 2. Then in classes 4 and 5, the partial scores are above average in almost all the domains of living conditions, except durables in class 4 and environmental problems in class 5.

If we look closer at the other characteristics of the households, we find that classes 4 and 5 have the highest proportions of persons living alone or lone parents, that classes 3 and 4 are more often living in large structures than on average, which is consistent with the proportion of environmental problems. More in the detail of the items, it appears clearly that the households of class 5 combine the highest proportions of difficulties in all the items of almost all the domains; the main difference with the households of class 4 being in the items relating to the standard of the dwelling. Class 3 is neatly disadvantaged on all the items relating to environmental quality, as for class 2, but with a smaller intensity and with almost no other material difficulty. Class 1 is clearly the one benefiting for the "best" material living conditions.

If we compare now the result from using the total score and the threshold set at 9 items and the result of our classification, a first main difference is in the proportion of households who will be said having bad living conditions: respectively 10,8% and (at least) 15%, if we consider in this case only the class 5, and up to almost 25% if we add classes 4 and 5. A second difference is that having more than two classes (one "good" and one "bad") allows to better differentiate between households who may have similar total score, but are different by the domains of their difficulties.

Conclusion

This paper was a first attempt to use the Kohonen algorithm to classify households' living conditions. It has proved a useful tool, in that it allows to construct classes that are neatly distinct and rather easily interpretable. In a first classification using the full qualitative information, we obtain classes of households which separates one class having "satisfying" material living conditions, and several classes of households who face more or less difficulties in almost all the domains considered; this classification makes appear neatly a small group of households characterized firstly by the lack of a basic comfort in their dwelling and the absence of some very common durable goods (telephone and television), as opposed to others groups characterized rather by serious deprivations in food and clothing. As for the results of the second classification based on the partial scores suggest strongly that poverty in terms of living conditions is not only a question of "level", but also a question of what contributes to the level.

The classifications show also that the most serious lacks or privations are often associated with a rather high score of (bad) living conditions, and a low income.

Nevertheless, because they use the qualitative dimension of the information, they show some associations of characteristics that will define particular groups of households, that do not appear if we "measure" living conditions with a score.

Figure 1 -Figure 4 -

 14 Figure 1 -Classification of the modalities in a one-dimension Kohonen map CLB1 CLW1 CLE1 PLT1 FMO1 PHO1 TEV1 CAR1 NCHF1 NVIP1 CLC1 CLR1 PLF1 PLM1 PLE1 PLS1

Figure 5 -

 5 Figure 5 -Means of 5 Kohonen classes of households based on the partial scores

Table 1 a -Characteristics of the classes of households based on the full set of items

 1

		A		B			C		All
		1	2	3	2+3	4	5	4+5
	Proportion of each class	70,6	8,7	5,5	14,2 13,3	1,9	15,2 100,0
	Average total score	2,2	9,3	5,4	7,8	4,7	9,1	5,3	3,5
	Average partial scores:
	dwelling 1	0,1	0,3	0,3	0,3	1,0	2,8	1,2	0,3
	dwelling 2	0,4	0,7	0,8	0,7	0,8	1,2	0,9	0,5
	environment	0,6	0,9	0,9	0,9	0,7	0,6	0,7	0,7
	durables	0,1	0,7	2,1	1,2	0,6	0,7	0,6	0,4
	deprivations	0,6	3,7	1,7	2,9	1,1	2,3	1,2	1,0
	Average income per C.U.	8378 5103 5457 5240 6797 4565 6518 7650
	Proportion of poor	6,1	27,6 19,8 24,6 17,0 33,3 19,0	10,7
	"Subjective" living conditions:
	with great difficulty	2,0	25,7 14,1 21,2	7,2	13,9	8,0	5,7
	with difficulty	9,6	29,8 20,3 26,1 12,4 10,7 12,2	12,3
	with some difficulty	27,6 34,8 40,1 36,8 31,1 41,0 32,3	29,6
	fairly easily	44,6	8,6	22,0 13,8 38,1 27,1 36,8	39,0
	easily and very easily	16,1	1,3	3,4	2,1	11,2	7,4	10,7	13,3
	Type of dwelling
	House, isolated	43,5 28,5 27,4 28,1 33,2 41,5 34,3	39,9
	House, in a neighborhood	21,5 20,3 22,0 21,0 18,8 36,6 21,0	21,4
	Structure <10 units	11,8 15,7 15,8 15,7 19,2 13,8 18,5	13,4
	Structure >=10 units	22,4 34,6 34,5 34,5 26,7	4,1	23,9	24,3
	Other	0,8	0,9	0,3	0,7	2,1	4,1	2,3	1,0
	Type of location
	Rural town	27,5 24,4 26,8 25,4 28,8 45,5 30,9	27,7
	City <10000 inh	11,2	8,7	10,7	9,5	11,4	6,5	10,8	10,9
	10000 to <100000 inh	19,0 20,9 21,2 21,0 20,9 17,1 20,4	19,5
	100000 to <2000000 inh	28,1 32,6 30,2 31,7 27,1 26,0 26,9	28,5
	Paris area	14,2 13,4 11,0 12,5 11,9	4,9	11,0	13,5
	Type of household
	one person household	21,7 33,0 27,1 30,7 33,8 48,8 35,7	25,1
	couple without child	28,2 17,3 21,8 19,0 25,4 15,5 24,2	26,3
	couple with child(ren)	40,1 30,3 35,6 32,4 31,9 19,5 30,4	37,5
	lone parent family	6,3	16,0 10,7 14,0	5,1	8,9	5,6	7,3
	other type	3,8	3,4	4,8	3,9	3,7	7,3	4,2	3,8
	Total number of persons in the
	household	2,6	2,5	2,7	2,6	2,4	2,0	2,3	2,6
	Age of the adults	47,9 46,5 41,9 44,7 40,6 58,2 42,8	46,7
	Number of persons aged <17	0,6	0,6	0,7	0,6	0,6	0,3	0,5	0,6
	Note: the characteristics that are significantly over-represented are in bold

Detailed proportions for the full set of items

		A		B			C		All
		1	2	3	2+3	4	5	4+5
	no bath or shower	2,7	0,0	0,0	0,0	0,7	71,5	9,6	3,4
	no hot running water	0,0	0,0	0,0	0,0	0,2 100,0	12,7	1,9
	not indoor toilet	1,8	0,9	0,6	0,8	2,3	49,6	8,2	2,6
	inappropriate heating system	1,9	17,5	16,7	17,2	45,5	30,1	43,6	9,5
	no separate kitchen	1,0	10,9	8,2	9,8	52,0	26,0	48,8 10,4
	rot in window frames	6,6	15,3	13,6	14,6	16,8	30,1	18,5	9,6
	damp walls, floor, etc.	11,5	19,8	21,2	20,3	22,2	31,7	23,4 14,6
	leaky roof	3,6	6,4	6,8	6,6	7,9	17,9	9,2	4,8
	shortage of space	9,2	19,1	20,6	19,7	18,6	15,5	18,2 12,0
	too dark, not enough light	6,7	13,6	13,0	13,3	16,3	21,1	16,9	9,2
	noise from outside	16,0	26,4	26,3	26,3	23,3	17,1	22,6 18,5
	noise from neighbors	8,3	13,2	15,5	14,1	14,2	9,8	13,6	9,9
	pollution	13,7	19,3	19,5	19,3	14,9	19,5	15,5 14,8
	vandalism, lack of safety	20,7	31,4	27,1	29,7	21,5	16,3	20,8 22,0
	no telephone	0,0	0,0	0,0	0,0	9,1	4,1	8,4	1,3
	no color TV	0,0	0,0	0,6	0,2	8,5	4,9	8,0	1,3
	no car	3,1	14,6	18,9	16,3	10,0	10,6	10,1	6,0
	no microwave	0,2	10,3 100,0	45,0	7,2	14,6	8,1	7,8
	no VCR	4,1	20,7	39,8	28,1	15,6	17,1	15,8	9,3
	no dishwasher	5,4	20,5	48,9	31,5	11,6	15,5	12,1 10,1
	Cannot afford
	meat/fish/chicken every second
	day if wanted	0,2	46,9	0,0	28,7	3,6	16,3	5,2	5,0
	new clothes	0,4	76,1	15,0	52,5	6,6	25,2	8,9	9,1
	keeping the dwelling warm
	enough	3,1	28,0	11,9	21,8	9,9	20,3	11,2	6,9
	replacing worn out furniture	26,2	92,0	68,6	83,0	43,9	61,0	46,0 37,3
	a week's holiday away once a
	year	24,5	79,0	57,9	70,8	36,7	69,1	40,8 33,5
	having friends/family for a drink
	or meal at least once a month	5,8	48,5	20,1	37,5	7,3	33,3	10,6 11,0

Table 2 -Setting a poverty line of living conditions

 2

			Cumulative
	SCORE Percent	Percent
	23	0,0	100,	0,0
	22	0,0	100	0,0
	21	0,0	100	0,0
	20	0,0	100	0,0
	19	0,1	99,9	0,1
	18	0,1	99,9	0,2
	17	0,1	99,8	0,3
	16	0,2	99,6	0,5
	15	0,6	99,4	1,1
	14	0,6	98,8	1,7
	13	0,9	98,2	2,6
	12	1,4	97,3	4,0
	11	1,6	95,9	5,6
	10	2,3	94,4	7,9
	9	2,9	92,1	10,8
	8	3,5	89,2	14,3
	7	3,9	85,7	18,2
	6	5,7	81,8	23,9
	5	6,5	76,2	30,4
	4	9,0	69,6	39,4
	3	8,6	60,6	48,0
	2	14,5	52,0	62,5
	1	14,3	37,5	76,8
	0	23,2	23,2	100

Table 3 -5 Kohonen classes of households based on partial scores

 3

		Basic score		Kohonen classes		All
		0	1	1	2	3	4	5
	Proportion of each class	89,2	10,8	48,7	12,5	14,0	9,9	15,0 100,0
	Average partial scores:
	dwelling 1	0,2	0,9	0,1	0,1	0,1	0,5	1,0	0,3
	dwelling 2	0,4	1,4	0,2	0,0	0,6	1,5	1,2	0,5
	environment	0,6	1,4	0,0	1,0	2,1	1,5	0,6	0,7
	durables	0,3	1,1	0,1	0,1	0,1	0,5	1,6	0,4
	deprivations	0,7	3,8	0,6	0,5	0,5	1,8	3,0	1,0
	Average total score	2,5	11,3	1,5	2,0	3,7	7,0	8,8	3,5
	Average income per C.U.	8021 4591 8251 9092 8671 6026 4613	7650
	Proportion of poor	7,8	34,7	6,3	3,9	4,4	15,3	33,4	10,7
	"Subjective" living conditions:
	with great difficulty	2,6	30,9	1,5	1,3	2,7	9,4	23,4	5,7
	with difficulty	10,5	28,0	8,1	7,9	9,5	20,9	26,7	12,3
	with some difficulty	29,3	32,7	27,3	25,8	28,6	39,5	34,9	29,6
	fairly easily	42,8	7,5	46,4	46,8	43,8	25,5	13,1	39,0
	easily and very easily	14,8	1,0	16,7	18,3	15,4	4,7	1,9	13,3
	Details
	no bath or shower	2,1	13,8	1,4	0,3	0,2	2,3	15,9	3,4
	no hot running water	1,0	9,9	0,4	0,5	0,2	1,7	10,0	1,9
	not indoor toilet	1,6	11,2	0,7	0,4	0,3	3,3	12,4	2,6
	no separate kitchen	7,3	36,5	3,0	4,4	5,0	24,2	35,6	9,5
	inappropriate heating system	8,5	18,2	7,0	3,9	4,1	14,1	24,6 10,4
	rot in window frames	6,9	31,5	2,1	0,0	12,1	31,9	24,6	9,6
	damp walls,floor,etc.	11,5	40,4	5,9	0,0	17,6	42,3	33,8 14,6
	leaky roof	3,4	17,1	1,6	0,0	3,9	16,1	12,8	4,8
	shortage of space	10,1	28,2	5,5	0,0	17,9	33,0	23,8 12,0
	too dark, not enough light	7,0	27,4	3,3	0,0	10,4	27,7	22,7	9,2
	noise from outside	15,8	40,4	0,0	21,2	63,6	45,5	16,2 18,5
	noise from neighbors	7,8	27,4	0,0	10,4	31,5	24,2	12,2	9,9
	pollution	13,0	29,7	0,0	18,1	53,7	31,4	12,6 14,8
	vandalism, lack of safety	19,3	43,8	0,0	50,2	58,6	45,9	19,8 22,0
	no telephone	0,7	6,2	0,1	0,1	0,0	1,4	7,4	1,3
	no color TV	0,9	4,9	0,1	0,0	0,1	0,8	7,6	1,3
	no car	4,2	20,8	1,4	1,7	2,0	9,1	26,3	6,0
	no microwave	5,8	24,4	2,3	1,6	1,7	9,7	35,3	7,8
	no VCR	7,0	27,9	3,4	2,0	3,2	9,7	39,9	9,3
	no dishwasher	7,7	29,7	3,8	3,2	3,7	14,4	39,8 10,1
	Cannot afford
	meat/fish/chicken every second
	day if wanted	1,4	35,2	1,0	1,0	0,8	8,0	23,4	5,0
	new clothes	3,7	53,6	2,9	2,0	2,9	17,7	35,5	9,1
	keeping the dwelling warm
	enough	3,0	39,1	1,6	2,4	2,1	12,8	28,7	6,9
	replacing worn out furniture	30,3	95,3	25,0	20,1	22,0	65,0	87,5 37,3
	a week's holiday away once a
	year	4,8	61,8	23,5	17,3	15,7	54,4	82,5 33,5
	having friends/family for a drink
	or meal at least once a month	26,5	92,0	3,7	3,0	3,2	20,9	42,0 11,0
	APPENDIX

Table A -Detailed frequencies of the items APPENDIX Table B -General descriptors

 A Poverty threshold at 50 % of the median income per C.U.

			Variable		Frequencies (%)	Means
	Type of	House, isolated	LOGT	1	39.9
	dwelling	House, in a neighborhood		2	21.4
		Structure <10 units		3	13.4
		Structure >=10 units		4	24.3
		Other		5	1.0
	Location	Rural town	TUR	0	27.7
		City <10000 inh		1	10.9
		10000 to <100000 inh		2	19.5
		100000 to <2000000 inh		3	28.5
		Paris area		4	13.5
	Type of	one person household	TYM	0	25.1
	household	couple without child		1	26.3
	(children taken	couple with child(ren)		2	37.5
	into account if	lone parent family		3	7.3
	<25 years old)	other type		4	3.8
	Total number of individuals	NBTOT			2.6
	Number of children <17 years old	NB17			0.6
	Mean age of persons at least 17 years old	AGEM			46.7
	"Subjective" living very difficult	SLS	1	5.7
	conditions	difficult		2	12.3
		rather difficult		3	29.6
		rather comfortable		4	39.0
		comfortable and very c.		5	13.3
	Monthly income (FRF) per C.U.(a)	REVUC			7650
	Monetary poverty	poor(b)	POOR	1	10.7
		non poor		2	89.3
	Score for the	Total (all domains)	TOT4			3.5
	material living	Dwelling, convenience	CLOGT			0.3
	conditions	Dwelling, problems	PLOGT			0.5
		Environment	ENVIR			0.7
		Durables	DURAB			0.4
		Deprivations	PRIV			1.0
	Source : Insee, European community households' panel, wave 3 (1996)
	(a) Consumption Unit, using the following equivalence scale : 1 -0.5 -0.3
	(b)

Acknowledgement

The SAS programs used for Kohonen classifications are due to Patrick LETREMY (SAMOS/MATISSE, Université de Paris 1)


Maria Lopez Abelairas

Luis Pereira De Oliveira

Jan J Verstraete
email: jan.verstraete@ifpen.fr

Application of Monte Carlo techniques to LCO gas oil hydrotreating: Molecular reconstruction and kinetic modelling

Keywords: Composition modelling, molecular reconstruction, kinetic modelling, kinetic Monte Carlo, stochastic simulation algorithm, LCO gas oil hydrotreating

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The pressure on the quality and the maximum impurities content of refinery products has been increased in the last decades due to the environmental concerns associated to their use. The suitable performance of the processes aimed to remove these impurities is essential to satisfy the specifications fixed for oil products. Among these processes, hydrotreating (HDT) is one of the most mature technologies. The use of hydrogen at high temperature over a catalyst bed allows to remove sulphur, nitrogen and other undesirable elements from the petroleum distillates such as naphthas or gas oils.

The correctness in the prediction of HDT process performance directly depends on the reliability of the used kinetic model. A lumping strategy, in which molecular components are grouped according to their global properties, has usually been applied in the kinetics models of complex hydrocarbons [START_REF] Quann | Structure-Oriented Lumping: Describing the Chemistry of Complex Hydrocarbon Mixtures[END_REF][START_REF] Quann | Building useful models of complex reaction systems in petroleum refining[END_REF][START_REF] Vynckier | Lumping Multicomponent Mix[END_REF] . However, this approach is no longer manageable due to the elevated number of generated groups and reaction pathways. The limited applicability of lumped models has resulted in the development of molecular-based kinetic models [START_REF] Quann | Building useful models of complex reaction systems in petroleum refining[END_REF][START_REF] Klein | Molecular modeling in heavy hydrocarbon conversions[END_REF][START_REF] Weekman | Lumps, models, and kinetics in practice[END_REF][START_REF] Sanchez | Kinetic Model for Moderate Hydrocracking of Heavy Oils[END_REF][START_REF] Verstraete | Modeling fixed-bed residue hydrotreating processes[END_REF][START_REF] García | Indepth modeling of gas oil hydrotreating: From feedstock reconstruction to reactor stability analysis[END_REF][START_REF] Pereira De Oliveira | A Monte Carlo modeling methodology for the simulation of hydrotreating processes[END_REF] . In this latter approach, reactant species and products are described by a selection of molecules or small groups of them which react following networks formed by reactions at the molecular level and distinguished by a global kinetic constant. The elemental steps are regrouped within the global reactions scheme, thereby reducing the number of reactive species and the size of reaction network. Following this approach, Lopez Garcia et al. [START_REF] García | Indepth modeling of gas oil hydrotreating: From feedstock reconstruction to reactor stability analysis[END_REF] have simulated the hydrotreating of LCO gas oils and, more recently, Pereira et al. [START_REF] Pereira De Oliveira | A Monte Carlo modeling methodology for the simulation of hydrotreating processes[END_REF][START_REF] Oliveira | Molecule-based kinetic modeling by Monte Carlo methods for heavy petroleum conversion[END_REF][START_REF] Oliveira | Simulating vacuum residue hydroconversion by means of Monte-Carlo techniques[END_REF] have proposed a kinetic methodology for its application on conversion of petroleum cuts based on the Stochastic Simulation Algorithm (SSA) developed by Gillespie [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF] .

These kinetics models required a molecular description of the feedstock. However, the high complexity of the petroleum cuts hinders their molecular characterisation even by using the most advanced analytical techniques. To surpass this drawback, a molecular representation of the feedstock can be determined by the application of a molecular reconstruction algorithm. This kind of algorithms provides a set of molecules from petroleum analyses and chemical knowledge. Several approaches have been performed in this sense, from the first model proposed by Liguras and Allen [START_REF] Liguras | Structural models for catalytic cracking. 2. Reactions of simulated oil mixtures[END_REF] to the stochastic reconstruction (SR) developed by Neurock et al. [START_REF] Neurock | Molecular representation of complex hydrocarbon feedstocks through efficient characterization and stochastic algorithms[END_REF] in which a Monte Carlo method is used to sample the objective function.

Hudebine [START_REF] Hudebine | Reconstruction moléculaire de coupes pétrolières[END_REF][START_REF] Hudebine | Reconstruction of Petroleum Feedstocks by Entropy Maximization. Application to FCC Gasolines[END_REF] has characterised different petroleum cuts by the application of a two-steps algorithm based on SR algorithm and on reconstruction entropy maximisation (REM). The combination of both methods (SR-REM) overcomes the drawbacks of its separated application, such as noise of the objective function and computational effort. Posteriorly, Pereira et al. [START_REF] Pereira De Oliveira | Développement d'une méthodologie de modélisation cinétique de procédés de raffinage traitant des charges lourdes[END_REF][START_REF] Oliveira | Molecular Reconstruction of Petroleum Fractions : Application to Vacuum Residues from Different Origins[END_REF] have modified Hudebine's method to introduce a genetic algorithm which performs the minimisation of the objective function.

The current study is focused on the validation of the SR-SEM and SSA algorithms developed by Pereira et al. on the hydrotreating of Light Cycle Oil (LCO) gas oils. A brief description of the stochastic methodology applied to both molecular reconstruction and reaction simulation has been developed. Then, the molecular reconstruction of eleven LCO gas oil has been carried out and the hydrotreating of three of them has been simulated at different operational conditions.

Methodology

The stochastic methodology applied in the current work consists of two steps. In the first one, a molecular reconstruction algorithm is used to model the feedstock by means of a set of molecules which mixture properties are similar to those of the feedstock. This algorithm comprises two coupled methods: stochastic reconstruction (SR) and reconstruction by entropy maximization (REM). The application of the SR method enables to obtain a set of molecules typical for a given type of oil fraction. Then, the molar fraction of molecules is adjusted using the REM method in order to improve the predicted mixture properties.

Finally, the generated set of molecules is introduced as input in the second step concerning the molecule-based kinetic modelling of the conversion process. A kinetic Monte Carlo (kMC) method is applied to model the chemical reaction system at a molecular level. This approach is an alternative to traditional deterministic kinetic methods in which a serial of differential equations (ODE) are involved to obtain a temporal evolution of the species in the system.

Molecular reconstruction

The SR method, from which a representative set of molecules is generated, uses probability distributions functions (PDF's) of structural attributes, such as number of rings or number of side chains, to characterise the oil fraction to be reconstructed [START_REF] Neurock | Molecular representation of complex hydrocarbon feedstocks through efficient characterization and stochastic algorithms[END_REF][START_REF] Trauth | Representation of the Molecular Structure of Petroleum Resid through Characterization and Monte Carlo Modeling[END_REF] . Feedstock analyses, which provide information about atomic abundance (elemental analysis) or structural classes (GC-MS), and the previous knowledge about petroleum cuts should be taken into account to select the proper structural attributes. The selected PDF's are sampled N times using a Monte Carlo method in order to generate an initial set of N molecules (Figure 1). To carry out the transition from PDF's to the set of molecules, the application of a building diagram and a serial of chemical rules is required to avoid the creation of impossible and unlikely molecules. The building diagram is a decision tree that defines the sequence in which PDF's should be sampled to proceed to the posterior assembly of the different attributes, which provides the structure of the molecule. During the assembly process, the application of chemical rules discards the molecules which do not satisfy the thermodynamic or likelihood criteria.

The application of the Monte Carlo sampling procedure led to the creation of a mixture of N molecules, each one with a molar fraction equal to 1/N. The average properties of this mixture are calculated and compared with the experimental data by means of an objective function.

An optimisation method is used to minimise this objective function though the modification of PDF's parameters. The simulated annealing [START_REF] Hudebine | Reconstruction moléculaire de coupes pétrolières[END_REF][START_REF] Trauth | Representation of the Molecular Structure of Petroleum Resid through Characterization and Monte Carlo Modeling[END_REF][START_REF] Trauth | Structure of complex mixtures through characterization, reaction and modelling[END_REF][START_REF] Hudebine | Molecular reconstruction of LCO gasoils from overall petroleum analyses[END_REF] and the genetic algorithms [START_REF] Verstraete | Molecular reconstruction of vacuum gasoils[END_REF][START_REF] Verstraete | Molecular reconstruction of heavy petroleum residue fractions[END_REF][START_REF] Schnongs | Reconstruction moléculaire de coupes pétrolières lourdes[END_REF] are both usual methods applied to the optimisation of the objective function in the field of the stochastic reconstruction. However, Schnongs [START_REF] Schnongs | Reconstruction moléculaire de coupes pétrolières lourdes[END_REF] has demonstrated that the genetic algorithms are more adapted to the stochastic reconstruction than the simulated annealing since the genetic algorithms are less sensitive to the initial value of the distribution parameters and it has also showed a lower oscillation of the value of the objective function.

Once the initial set of molecules is obtained, the second step based on the reconstruction by entropy maximisation (REM) is applied. During this step, the molar fraction of the molecules in the mixture is modified to obtain a mixture with properties closer to those experimentally determined [START_REF] Hudebine | Reconstruction moléculaire de coupes pétrolières[END_REF] .To carry out the adjustment of the parameters, the REM method, based on Shannon's information theory [START_REF] Shannon | A Mathematical Theory of Communication[END_REF] , uses the maximisation of an information entropy criterion (E(x i )) given by the Eq. (1):

𝐸(𝑥 𝑖 ) = -∑ 𝑥 𝑖 • ln (𝑥 𝑖 ) 𝑁 𝑖=1 (1)
with

∑ 𝑥 𝑖 = 1 𝑁 𝑖=1 (2)
where x i corresponds to the molar fraction of molecule i and N is the number of molecules present in the initial mixture determined by applying the SR method. This criterion guaranties that, without constrains or available information, a given molecule cannot be preferred to others, so that the distribution of the predefined set of molecules remains uniform. The introduction of a matrix of constrains from analytical data, mixing rules and the initial set of molecules distorts this uniformity in order to match this information. These constrains are introduced in the entropy maximisation criterion by means of the Lagrange multiplier method.

In the case of linear constrains were used, the problem is reduced to an optimisation of J parameters of a no-linear equation that can be performed by the conjugate gradient method.

Once the Lagrange parameters are determined, the molar fraction of the molecules is calculated in order to obtain the representative mixture of molecules.

Reaction modelling

The goal of the kMC method is to describe the chemical reaction system at the molecular level by following the transformations of a discrete population of molecules. Contrary to deterministic methods, the temporal evolution of the system is not described by a set of differential equations, but by a single probabilistic function termed the chemical master equation (CME) [START_REF] Gillespie | A rigorous derivation of the chemical master equation[END_REF][START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF] . This equation determines the different status of the system in a sequence of discrete-time steps. The evolution of the system is a result from chemical reactions, those have associated a probability function, which depends on the reactant molecule and the rate constant.

However, CME cannot usually be solved by analytical or numerical methods due to its high complexity. In this case, a Monte Carlo procedure developed by Gillespie [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF] and termed Stochastic Simulated Algorithm (SSA) is used to solve CME. This method is based on an event-space approach, i.e., the evolution of the system is followed event by event, resulting in disparate time intervals between two consecutive events or, in this particular case, reactions.

The scheme followed in the application of the kMC method has been illustrated in Figure 2.

The initial step of the procedure consists in the "molecular discretisation" of the set of molecules generated by applying the SR-REM method. During this step, the replication of the molecules in the mixture is performed according to their molar fraction and it is based on the same principle as the procedure developed by Hudebine and Verstraete [START_REF] Hudebine | Molecular reconstruction of LCO gasoils from overall petroleum analyses[END_REF] for the creation of a representative set of molecules for gasoline fractions. A replication factor (F rep ), defined as the maximum number of molecules in the discrete mixture, is multiplied by the mole fraction of each molecule in order to calculate the number of replications. F rep also defines the lowest mole fraction to be retained since the molecules with a mole fraction lower than 0.5/F rep are discarded. The factor 0.5 arises from the rounding off to the closest integer.

In the next step of the kMC procedure, all potential reaction events should be identified and listed taking into account a set of reaction rules that generates the reactions by inspecting the structure of the molecules. Once all potential reaction events are identified, the normalised probability for each reaction event is calculated as the ratio of its stochastic rate (R v ) to the sum of the rates of all reaction events (M), as shown in Eq. (3):

𝑃 𝑣 = 𝑅 𝑣 ∑ 𝑅 𝑖 𝑀 𝑖 (3) with 𝑅 𝑣 = ℎ 𝑣 • 𝑐 𝑣 (4)
where h v is the number of distinguishable combinations of the reactant molecule and c v is the stochastic rate parameter. For monomolecular reactions, the number of distinguishable combinations of the reactant molecules h v equals 1, while for bimolecular reactions, it is equal to the number of available molecules of the second reactant. The stochastic rate parameter c v is closely linked to the deterministic rate parameter (k v ) of each reaction. In the case of monomolecular reactions, k v and c v are equal, while for bimolecular reactions, c v is equal to the ratio of k v by reaction volume V. [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF] The effect of the reaction temperature on the stochastic rate is taking into account by means of the Arrhenius expression reflected in Eq. ( 5):

𝑐 𝑣 = 𝑐 𝑣 (𝑇 𝑟𝑒𝑓 ) • exp (- 𝐸 𝑎 𝑅 • ( 1 𝑇 - 1 𝑇 𝑟𝑒𝑓 )) (5)
where c v (T ref ) is the stochastic rate constant of the reaction type v at the reference temperature T ref , E a is the activation energy, R is the ideal gas constant, T is the system temperature and

T ref is the reference temperature.
The subsequent step of the algorithm is aimed to the determination of the cumulative reaction probability distribution (D R ), which grouped all M reactions that can take place in the mixture at a given time t.

A first random number (RN 1 ) is obtained to determine the reaction time step (∆t) according to the Eq. ( 6), as proposed Gillespie [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF] :

∆𝑡 = - ln (𝑅𝑁 1 ) ∑ 𝑅 𝑖 𝑀 𝑖 (6)
The distribution D R is randomly sampled in order to select the next reaction (µ) to be executed in the next reaction time step. For this purpose, a second random number between 0 and 1 is drawn (RN 2 ) and used to select the next reaction from the cumulative distribution D R , as shown in Eq. ( 7):

𝐷 𝑅 (𝜇 -1) < 𝑅𝑁 2 ≤ 𝐷 𝑅 (𝜇) (7)
Once the reaction and the time step have been defined, the system is updated by the execution of the selected reaction, which implies the replace of the reactant(s) by the product(s) and the increase of simulation time in ∆t. The simulation continues reaction by reaction until the final simulation time is reached and a set of product molecules is obtained.

Due to its stochastic nature, the algorithm should be executed several times in order to provide several sets of molecules which average can be a proper representation of the molecules in the reaction system.

Application to the molecular reconstruction of LCO gas oils

LCO gas oils are gas oil fractions obtained from the catalytic cracking process and composed by a mixture of saturates (paraffins and naphthenes), olefins and aromatic hydrocarbons [START_REF] García | Indepth modeling of gas oil hydrotreating: From feedstock reconstruction to reactor stability analysis[END_REF] .

These mixtures also show a marginal presence of nitrogen but an important content in sulphur, which is present in heterocyclic structures, such as thiophenes and benzothiophenes [START_REF] García | Analysis of aromatic sulfur compounds in gas oils using GC with sulfur chemiluminescence detection and high-resolution MS[END_REF] . Their boiling point is in the range of 120°C to 450°C, which corresponds to the molecules with a number of carbons between C 8 and C 30

29

. According to this data and the mass spectrum, LCO gas oils are formed by molecules with none or one poly-aromatic structure with a number of rings which usually ranges from 0 to 4 [START_REF] Cooper | Aromatic saturation of distillates : an overview[END_REF][START_REF] Castex | Analyse des kérosènes et des gas oils moyens par spectrométrie de masse à moyenne résolution[END_REF][START_REF] Fafet | New Developments in Mass Spectrometry for Group-Type Analysis of Petroleum Cuts[END_REF] .

In the current work, it has been performed the molecular reconstruction of eleven LCO gas oils using the SR and REM algorithms. The PDF's of molecular attributes in Table 1 is used to generated the set of molecules that would represent the LCO gas oils. Two kind of PDF's (histograms and gamma functions) were used to defined molecular attributes. The use of histograms is limited to the attributes with a narrow range of possible values (less than three) while the gamma functions describe the attributes with a higher number of possible values. These PDF's have been defined taking into account the previous knowledge about LCO gas oils composition. Firstly, it is considered that the molecules in mixture are only constituted by carbon , hydrogen and sulphur due to the negligible presence of nitrogen species. Then, nitrogen content is fully assigned to carbon fraction, the most abundant element in mixture.

Another hypothesis is that paraffins, which would also include the marginal olefins fraction, are formed by alkyl chains with at least 8 carbons. It is also considered that, according to experimental data, each molecule is formed by a maximum of 4 cycles, in which linear alkyl chains can be the unique substituents, and that the sulphur containing compounds present in the mixture only belong to the thiophenes family and its derivatives.

In order to define the relationships between distributions and also identify the sampling steps, a building diagram is needed for the reconstruction of the molecules in LCO gas oils (Figure 3). First, the type of molecule to be constructed (saturate or aromatic) is determined by distribution 1. In the case of saturate molecules, distribution 3 allows to determine the number of cyclohexane rings (ranged from 0 to 4). When the number of rings present is null, the molecule to be constructed is a paraffin, which chain length is determined by the application of distribution 8. Contrary, in the case of a naphthenic units with between 1 and 4 cyclohexane rings, it is necessary to determine the number (distribution 6) and type (distribution 7) of the side chains. In the case of aromatic molecules, the number of benzenes, cyclohexanes and thiophenes is determined by distributions 2, 3 and 4, respectively. These three distributions are sampled until the total number of rings is inferior or equal to 4. The number of benzenes is always 1 or higher to avoid the creation of naphthenic molecules in this route. Besides, according to experimental data, most of the thiophene rings do not contain cyclohexane rings, so that distribution 3 is not sampled when thiophene rings are present.

Once the number and type of cycles is defined, the distributions 5 and 6 are used to find the number of side chains in naphthenic and aromatic rings, respectively, and the distribution 7 is applied to determine the type of chain.

Elemental analysis (carbon, hydrogen, sulphur and nitrogen), average molecular weight, specific gravity, simulated distillation, mass spectrometry (MS), [START_REF] Quann | Structure-Oriented Lumping: Describing the Chemistry of Complex Hydrocarbon Mixtures[END_REF] H NMR and [START_REF] Liguras | Structural models for catalytic cracking. 2. Reactions of simulated oil mixtures[END_REF] C NMR were used to characterise the LCO gas oils analytically. Elemental analysis is used to remove the molecules which are composed by an element with a null experimental value. For example, if sulphur presence has not been detected, all sulphureted molecules (benzothiophenes and dibenzothiophenes) are eliminated. In a similar way, all molecules belonging to a spectrometric family which has not been found in mass spectrometry analyses are suppressed.

Spectrometric families are determined according to the number of aromatic and thiophene rings in compounds. Thus the following families can be distinguished: saturates (paraffins and cycloparaffins), monoaromatics (alkylbenzenes, tetralins…), diaromatics (naphthalenes, diphenyls…), triaromatics and compounds with more than three aromatic rings (anthracenes, phenanthrenes…), benzothiophenes and dibenzothiophenes. Meanwhile, simulated distillation enables to discard the molecules with a boiling temperature 5°C under the initial or 5°C over the final boiling point. This correction is introduced due to the bias between the real boiling point and the value obtained from simulated distillation curve 33 . The "experimental" molecular weight was calculated using an API correlation based on specific gravity and simulated distillation 34 .

Sensitivity of the stochastic method

The average properties of the simulated mixture of molecules are calculated by the application of a number of linear mixing rules. An objective function (OF) is used to perform the comparison between calculated and measured properties. This function is the sum of the relative deviations between calculated and experimental data and it is given by Eq. ( 8):

𝑂𝐹 = 1 𝑁 𝑝 • ∑ 𝛿 𝑖 𝑁 𝑝 𝑖=1 (8)
where N p is the number of properties present in objective function and δ i is the relative deviation between the calculated and experimental values of property i which is calculated following Eq. ( 9):

𝛿 𝑖 = 1 𝑁 𝑀,𝑖 • ∑ �𝑋 𝑗,𝑖 𝑒𝑥𝑝 -𝑋 𝑗,𝑖 𝑐𝑎𝑙𝑐 � 𝑋 𝑗,𝑖 𝑒𝑥𝑝 𝑁 𝑀,𝑗 𝑗=1 (9)
where N M,i is the number of measurements of property i (for example, the number of points in the simulated distillation) and X j,i exp and X j,i calc are, respectively, the experimental and calculated values of measurement j of property i.

An elitist genetic algorithm developed by Pereira et al. [START_REF] Oliveira | Molecular Reconstruction of Petroleum Fractions : Application to Vacuum Residues from Different Origins[END_REF] was applied in order to minimise this objective function by the modification of the parameters of the PDF's for the structural attributes. Genetic algorithms [START_REF] Verstraete | Molecular reconstruction of vacuum gasoils[END_REF][START_REF] Verstraete | Molecular reconstruction of heavy petroleum residue fractions[END_REF][START_REF] Schnongs | Reconstruction moléculaire de coupes pétrolières lourdes[END_REF] , as well as simulated annealing [START_REF] Hudebine | Reconstruction moléculaire de coupes pétrolières[END_REF][START_REF] Trauth | Representation of the Molecular Structure of Petroleum Resid through Characterization and Monte Carlo Modeling[END_REF][START_REF] Trauth | Structure of complex mixtures through characterization, reaction and modelling[END_REF][START_REF] Hudebine | Molecular reconstruction of LCO gasoils from overall petroleum analyses[END_REF] and particle swarm optimization [START_REF] Kennedy | Swarm intelligence[END_REF] , are robust global optimizers, a characteristic required for their application in stochastic models. Besides, according to its elitist nature, the applied methodology allows to a certain percentage of individuals to move to the following generation (in this case, a percentage of 50% was considered), with which the best individuals (i.e. those with the lowest OF values) remain in the next generation.

The number of molecules is an important parameter in molecular reconstruction since an excessively low number would lead to an inadequate representation of the attributes of the represented mixture. Otherwise, the use of a large number of molecules would imply a significant increase in calculation time. In order to determine the optimal value for this parameter, the sensibility of the stochastic method has been tested by the generation of 100 mixtures of N molecules from a set of parameter distributions determined for a LCO gas oil. Then, a set of 5000 molecules was considered as optimal to perform the stochastic reconstruction of LCO gas oils.

Construction of the molecular representation of LCO gas oils

The analyses used in this case for the molecular reconstruction were the elemental analysis, MS and simulated distillation. The elemental analysis provides information on the number of thiophenes, while the information on the chemical structure is provided by MS and by simulated distillation in the case of the average molecule size and size distribution. The experimental and simulated properties were compared in Table 2 for a LCO gas oil.

Properties that were not used in reconstruction step (specific gravity, molecular weight, 1 H NMR and [START_REF] Liguras | Structural models for catalytic cracking. 2. Reactions of simulated oil mixtures[END_REF] C NMR) were calculated from the predicted molecular mixture and also compared with the corresponding experimental values.

The properties of the mixture obtained after the application of the SR step are mostly in good agreement with the analytical results. Some deviations are nevertheless observed in the simulated distillation curve. These differences may be explained by the fact that the distillation curve is calculated under the assumption that all molecules are perfectly separated by the use of increasing boiling points. However, this is not strictly true (especially in the region near the initial and final distillation points), since the chromatographic columns are not completely apolar and separation is influenced by the polarity of the molecules [START_REF] Hudebine | Molecular reconstruction of LCO gas oils from overall petroleum analyses[END_REF] .

On the other hand, the coupled application of SR and SEM schemes led to a significant improvement in the agreement between predicted and experimental data for most of the properties (Table 2). The most meaningful effects were observed on more specific and detailed analyses, such as NMR and MS, while only small improvements could be detected in the more global measurements, such as in the case of elemental analysis.

On the following step, a total of eleven LCO gas oils have been reconstructed, and the properties of the predicted mixtures have been compared with the corresponding experimental data in order to perform the validation of the current methodology (Figure 5). It was found a good agreement with the analytical results for most of the evaluated properties. However, a deviation was detected for the predicted average molecular weight, which was mostly placed below the value calculated from experimental data. This may be explained by the alterations in the determination of the distillation curve, which imply a reduction in the boiling point of the molecules with respect to their normal value. It leads to the prediction of molecules smaller than their counterpart and, consequently, to reduce the average molecular weight of the mixture.

These results show that the predicted set of molecules are proper molecular representations of the corresponding LCO gas oils. Some of these simulated mixtures have been used as inputs for the kinetic model of the hydrotreating of LCO gas oil.

Application to the hydrotreating of LCO gas oils

Once the representative set of molecules has been generated for each LCO gas oil, the hydrotreating reactions are simulated by applying the kMC method. The experimental data used in this step has been collected from previous studies developed at IFP Energies nouvelles about the hydrotreating of these oil fractions [START_REF] Lopez Garcia | Analyse de la réactive des composés soufrés dans les coupes pétrolières: cinétique et modélisation de l'hydrotraitement[END_REF] .

Three main type of reactions are involved in the hydrotreating process: hydrodesulphurisation (HDS), hydrodenitrogenation (HDN) and hydrogenation of benzene and cyclohexane rings and of olefins. However, the reactions of hydrodenitrogenation and hydrogenation of olefins will be neglected due to the lack of nitrogen and olefins in the mixture. Similarly, the HDS reactions are limited to heterocyclic structures since it has been considered that sulphur is only present in the mixture as thiophenes and their derivatives.

Hydrogenation of aromatic rings (HDA), an exothermic and reversible reaction, is thermodynamically limited at typical hydrotreating conditions, so that the reverse reaction of dehydrogenation should also be considered. Korre et al. [START_REF] Korre | Hydrogenation of polynuclear aromatic hydrocarbons.2. Quantitative Structure/Reactivity Correlations[END_REF] have indicated that no partially hydrogenated ring compounds, such as dihydronaphthalene and hexahydrophenanthrene, were found in the product stream, which shows that hydrogenation proceeds in a ring-by-ring manner. Besides it was detected that the higher number of aromatic rings, the higher reactivity of the molecule. This fact implies that polyaromatic species, such as pyrene, would be the most reactive whereas that the HDA of monoaromatics, such as alkyl benzenes and tetralines, would be the slowest. Other qualitative trend observed during the process was that, for the molecules with more than two aromatic rings, the hydrogenation of the rings in the middle of the molecule goes faster than in the case of the external rings [START_REF] Stanislaus | Aromatic hydrogenation catalysis: a review[END_REF][START_REF] Korret | Polynuclear Aromatic Hydrocarbons Hydrogenation . 1 . Experimental Reaction Pathways and Kinetics[END_REF] . To this, it should be added the fact that the present of alkyl substituents and of cyclohexane rings increases the hydrogenation rate [START_REF] Korre | Hydrogenation of polynuclear aromatic hydrocarbons.2. Quantitative Structure/Reactivity Correlations[END_REF] .

HDS is an exothermic and basically irreversible reaction at the hydrotreating conditions [START_REF] Girgis | Reaction Networks, and Kinetics in High-Pressure Catalytic Hydroprocessing[END_REF] through which sulphur is removed from the petroleum cuts in form of hydrogen sulphide. The removal of this element is one of the main aims of hydrotreating process due to its catalyst poisoning and pollution potential [START_REF] Oliveira | Simulating vacuum residue hydroconversion by means of Monte-Carlo techniques[END_REF] . The HDS reactivity of sulphureted species closely depends on their structure. According with this, non-heterocyclic compounds present a very higher reactivity at hydrotreating conditions than heterocyclic species, such as thiophenes family and derivatives. Different reactivities are also found within this last compounds family, in which thiophenes (Ts) and benzothiophenes (BTs) show a higher reactivity than dibenzothiophenes (DBTs) [START_REF] Whitehurst | Present State of the Art and Future Challenges in the Hydrodesulfurization of Polyaromatic Sulfur Compounds[END_REF] . The hydrotreating rates of the heterocyclic species are also influenced by the presence of substituent groups in the positions adjacent to the sulphur atom [START_REF] Pereira De Oliveira | A Monte Carlo modeling methodology for the simulation of hydrotreating processes[END_REF] . The desulphurisation process in these species goes through two different pathways:

(1) the direct removal of the sulphur atom without ring hydrogenation (hydrogenolysis pathway) or (2) the saturation of the thiophene and the fused benzene rings prior to the removal of the sulphur atom from the thiophene ring (hydrogenation pathway). The HDS of Ts, BTs and unsubstituted DBTs mostly follows the hydrogenolysis pathway whereas hydrogenation becomes the preferred pathway for DBTs at the same time that the degree of substitution increases [START_REF] Girgis | Reaction Networks, and Kinetics in High-Pressure Catalytic Hydroprocessing[END_REF][START_REF] Whitehurst | Present State of the Art and Future Challenges in the Hydrodesulfurization of Polyaromatic Sulfur Compounds[END_REF] .

Taking into account the above information, the reaction types reflected in Table 3 has been considered during the stochastic modelling of LCO hydrotreating. A set of reaction rules is also established in order to identify the reaction events from the structure of the molecule.

Hydrogenation of aromatic rings and dehydrogenation of saturated rings have been included

in the current model as ring-by-ring reactions while it have been considered that sulphur is removed from thiophene rings as H 2 S. Besides, the rate constants of the HDT reactions must be specified.

In the case of hydrogenation and dehydrogenation reactions, their reaction rates can be calculated from Eq. ( 10) and (11), respectively:

𝑅 𝐻𝑦𝑑𝑟𝑜 = 𝑐 𝐻𝑦𝑑𝑟𝑜 • 𝑝 𝐻 2 (10)
𝑅 𝐷𝑒ℎ𝑦𝑑𝑟𝑜 = 𝑐 𝐷𝑒ℎ𝑦𝑑𝑟𝑜 = 𝑐 𝐻𝑦𝑑𝑟𝑜 𝐾 𝑒𝑞 (11) where c Hydro and c Dehydro are the stochastic reaction constants for hydrogenation and dehydrogenation (s -1 ), K eq is the equilibrium constant (mol•(mol•atm) -1 ), and p H2 is the partial pressure of hydrogen (atm).

The expressions for the calculation of hydrogenation rate constants and the equilibrium constants have been developed in previous studies [START_REF] Pereira De Oliveira | Développement d'une méthodologie de modélisation cinétique de procédés de raffinage traitant des charges lourdes[END_REF] , starting from the quantitative structure/reactivity correlations (QS/RCs) proposed by Korre et al. [START_REF] Korre | Hydrogenation of polynuclear aromatic hydrocarbons.2. Quantitative Structure/Reactivity Correlations[END_REF] . They are given in Eq. ( 12) and ( 13), respectively:

ln�𝑐 𝐻𝑦𝑑𝑟𝑜 � = -4.65 -7.06 • 𝑛 𝐻 1.18 -1.08 • � ∆𝐻 𝑅 0 𝑅 � +0.616 • 𝑁 𝐴𝑅 + 0.330 • 𝑁 𝑆𝑅 -2.80 • 𝑁 𝑇𝑅 - 𝐸 𝑎 𝑅 • � 1 𝑇 - 1 𝑇 𝑅𝑒𝑓 � ( 12
)
ln�𝐾 𝑒𝑞 � = 2.952 -13.215 • 𝑛 𝐻 + 5.196 • 10 -3 • � ∆𝐻 𝑅 0 𝑅 � -0.784 • 𝑁 𝑆𝑅 - ∆𝐻 𝑅 𝑅 • � 1 𝑇 - 1 𝑇 𝑅𝑒𝑓 � (13)
where n H is the stoichiometric coefficient for molecular hydrogen, ∆H R 0 is the heat of reaction at 25°C (J•mol -1 ), N AR is the number of aromatic rings in the molecule, N SR is the number of saturated rings, N TR is the number of thiophenes fused to the aromatic ring, E a is the activation energy (J•mol -1 ), R is the ideal gas constant (J•(K•mol) -1 ), T is the temperature of the system (K), and T ref is the reference temperature (K).

The values of the coefficients in the K eq expression proposed by Korre et al. [START_REF] Korre | Hydrogenation of polynuclear aromatic hydrocarbons.2. Quantitative Structure/Reactivity Correlations[END_REF] have been modified in order to obtain a better agreement between predicted and experimental equilibrium data. The experimental data of a set of 21 compounds, including anthracene, phenanthrene and benzene, have been used to determine the new constants. For the hydrogenation rate constant c Hydro , the N TR term was added to reflect the effect of the presence of thiophene rings on the hydrogenation of aromatic rings, as observed in hydrodesulphurisation kinetics [START_REF] Girgis | Reaction Networks, and Kinetics in High-Pressure Catalytic Hydroprocessing[END_REF] . An exponential coefficient was also added to the n H term to improve the prediction of the kinetics of benzoic compounds. The hydrogen exponential coefficient, the coefficient for the thiophene correction and the independent term in the expression for k SR were determined by minimising the difference between the predicted and experimental data for the gas oil LCO1 (Table 4) during its hydrotreating at 300 °C (Figure 6.a).

A very good agreement between experimental and predicted data was obtained for monoaromatic compounds. Besides, the fast hydrogenation of the aromatic compounds with more than two aromatic rings was also properly reflected but some differences could still be observed for saturates and diaromatic compounds. Despite they followed the same tendency than the predicted data, the experimental data showed a larger removal of diaromatic compounds and also a higher concentration of saturated compounds. The introduction of hydrogen exponential coefficient may be the cause for this deviation since, although it provides a better prediction of benzoic compounds kinetics, it worsens the results for other molecules, such as naphthalenes.

HDS reactions have been classified in five families according to the structure of the reacting molecule and a rate parameter was assigned to each one (Table 4). The reaction rates are the component directly given by the corresponding rate constant following Eq. ( 15):

𝑅 𝐻𝐷𝑆 = 𝑐 𝐻𝐷𝑆 (15)
The rate constants c HDS were obtained from the fit to the experimental data of the LCO gas oil which was previously used to determine the coefficients in Eq. ( 15) (Figure 6.b). In this case, the predicted curves perfectly fit to experimental data of both BT and DBT compounds.

In order to validate the current modelling methodology, the hydrotreating of the three LCO gas oils in Table 4 has been simulated at different conditions of temperature and partial H 2 pressure. These LCO gas oils present different composition profiles. LCO1 and LCO2 both present a significant content in sulphur aromatics (mainly DBT in LCO1 and BT in LCO2) while, in the case of LCO3, the content in these compounds is very lower (less than 1.5% wt) and the addition of sulphides such as dimethyl disulphide (DMDS) is required in order to reach an adequate partial pressure of H 2 S over the process. Besides, LCO3 is mainly composed by saturated compounds (almost 60% wt) unlike the other two gas oils, in which diaromatics are the most abundant compounds.

The predicted results after 1 h of operation have been compared with analytical data in the parity plots depicted in Figure 7. A good agreement with experimental data was found for saturated and monoaromatic compounds with an average relative error lower than 14%. The prediction of the species with three or more aromatic rings was accurate in most of the experimental points and a very good fit was also found for sulphur aromatic compounds.

Performance variations due to modifications in operational conditions, i.e. temperature and hydrogen partial pressure, were contemplated and well-reflected on the current model. On the other hand, as in the case of LCO1, some differences were also detected for diaromatic compounds in LCO3. High complexity of mixtures, which leads to the use of correlations to determine the reactivity of each molecule in hydrogenation/dehydrogenation processes, implies a greater difficulty to accurately predict the content in all aromatic species, although a good global agreement between predicted and experimental data was reached.

Conclusions

A stochastic two-steps procedure has been described and validated on the hydrotreating of three LCO gas oils at different operational conditions. In the first step, a molecular reconstruction algorithm is used to represent the feedstock by means of a set of molecules which mixture properties are similar to those of feedstock. This algorithm is formed by the combination of two methods: stochastic reconstruction (SR) and reconstruction by entropy maximization (REM). The SR method leads to obtain a set of molecules typical of a given feedstock while the molar fraction of molecules is adjusted using the REM method in order to improve the predicted mixture properties. Then, the generated set of molecules is used as input in the second step concerning the molecule-based kinetic modelling of the hydrotreating process and in which a kinetic Monte Carlo (kMC) method is applied. The expressions for the calculation of adsorption, equilibrium and hydrogenation/dehydrogenation rate constants have been obtained by applying the quantitative structure/reactivity correlations (QS/RCs).

The current stochastic methodology has demonstrated to be a valid tool for the reconstruction of LCO gas oils by combining their available analytical data. This procedure can be used not only to get the molecular information that is not provided by analytical data but also to obtain an accurate molecular representation which can be used as input in kinetic models. Monte Carlo techniques have also given good results in their application on the kinetic simulation of the hydrotreating of LCO gas oils. Moreover, the simulation of the reactions does not require a pre-defined kinetic network due to it is generated as the reactions proceed. Besides, the use of the corrected QS/RCs led to a good prediction of molecular reactivities in complex mixtures of hydrocarbons such as LCO gas oils.

  represented. It was observed that the increase of the size of the mixture leads to a drop of the objective function until attaining an asymptotic value at N = 5000 (Figure 4.a). The increase of the number of molecules also implies a rapprochement of the minimum and maximum values and a reduction in the standard deviation of the objective function. However, the drop of the standard deviation is very mild at N > 5000, while the calculation time continues to rise linearly with the size of the mixture (Figure 4.b).

Figure captions Figure 1 .Figure 2 .Figure 3 .

 captions123 Figure captions

Figure 4 .Figure 5 .Figure 6 .Figure 7 .

 4567 Figure 4. Evolution as a function of the number of molecules in the predicted mixture of (a) the average value of the objective function (white diamonds) and the maximum (black triangles) and minimum (black squares) values (b) the CPU time required (white diamonds) and the standard deviation of the objective function (black crosses)

Figure 1

 1 Figure 1.

Figure 2

 2 Figure 2.

Table 1 .

 1 Structural attributes of LCO gas oil molecules

	Structural attribute	Values	Distribution	Parameter names
	Type of molecule	0 or 1	Histogram	0
	Number of benzene rings per core	1, …, 4	Gamma	1 and 2
	Number of naphthenic rings per core	0, …, 4	Gamma	3 and 4
	Number of thiophenes per core	0 or 1	Histogram	5
	Acceptance probability for aromatic carbon CH	0 or 1	Histogram	6
	Acceptance probability for naphthenic carbon CH 2	0 or 1	Histogram	7
	Type of the side chains	0, 1 or 2	Histogram	8 and 9
	Length of the paraffinic chains	> 0	Gamma	10 and 11

Table 2 .

 2 Comparison between experimental and calculated properties after the SR and REM steps

	Experimental	Simulated (after SR)	Simulated (after SR-REM)

Table 3 .

 3 Reaction types and their parameters (c)

	Type of reaction					Reaction	c (s -1 )
	Hydrogenation		+ H 2 2	Eq. (10)
	Dehydrogenation					+ H 2 2	Eq. (11)
	Hydrogenolysis of thiophene and their derivatives
			+	4	H 2	CH 3	+ S H 2	138.90
		S
	Hydrogenolysis of BT and their derivatives
			+	3	H 2	CH 3 + S H 2	55.60
		S
	Hydrogenolysis of DBT and their derivatives
	Hydrodesulphurisation	S				+ H 2 2	+ S H 2	19.40
	Hydrogenolysis of 4-alkyl-DBT and their derivatives
						+ H 2 2	+ S H 2	0.56
	R	S				R
	Hydrogenolysis of 4,6-dialkyl-DBT and their
	derivatives
						2 + H 2	+ S H 2	0.06
	R	S	R			R	R

Table 4 .

 4 Hydrotreating conditions (partial H 2 pressure and temperature) and properties from SR-REM algorithm for the three LCO gas oils used as input in the kMC method

		LCO1	LCO2	LCO3
	Hydrotreating conditions
	pH 2 (atm)	69.1	95.1	46.1
	T (°C)	300/330	300/330 300/310/320
	Composition
	C (%wt)	88.01	88.57	87.23
	H (%wt)	10.51	9.74	12.57
	S (%wt)	1.48	1.69	1.45
	Saturates (%wt)	28.58	19.85	59.84
	Monoaromatics (%wt)	16.14	22.23	9.55
	Diaromatics (%wt)	34.24	38.96	26.20
	Triaromatics and more (%wt)	11.49	8.95	3.19
	BT (%wt)	2.70	6.39	0.67
	DBT (%wt)	6.86	3.61	0.56