Jean-Jerôme J Da Costa
email: jean-jerome.da-costa-soares@ifp.fr

Fabien Chainet

Benoît Celse

Marion Lacoue-Nègre

Cyril Ruckebusch

Noémie Caillol

Didier Espinat

Kriging Modeling to Predict Viscosity Index of Base Oils

Keywords: Kriging, Multilinear Regression, Hydrotreatment, Hydrocracking, Base oil, Viscosity Index, Total liquid effluent

published or not. The documents may come

Introduction

Despite of high costs of investment due to rough operating conditions, hydrocracking has become an attractive option for the upgrading of vacuum gas oil (VGO). The main reason is the increasing demand for middle distillates (e.g. kerosene and diesel cuts) [START_REF] Motaghi | Slurry-phase hydrocracking-possible solution to refining margins[END_REF][START_REF] Butler | Maximize liquid yield from extra heavy oil Next-generation hydrocracking processes increase conversion of residues[END_REF] . VGO cut is a petroleum fraction which contains compounds with boiling points above 370°C [START_REF] Riazi | Characterization and properties of petroleum fractions[END_REF] . It is obtained by distillation of atmospheric residue of crude oil [START_REF] Wauquier | Le raffinage du pétrole: Pétrole brut, produits pétroliers, schémas de fabrication[END_REF] . The relatively high amounts of organic nitrogen and sulfur in VGO makes it not valuable as they are undesirable in the final products (e.g. diesel, base oil). The purpose of hydrocracking units is to convert VGO into more valuable fractions, conforming to ever more stringent product quality specifications [START_REF] Rana | A review of recent advances on process technologies for upgrading of heavy oils and residua[END_REF] . Many types of hydrocracking process are currently used in the petroleum industry. It consists in two successive steps: hydrotreatment step (HDT) which is performed to remove organic nitrogen and sulfur, by hydrodenitrogenation (HDN) and hydrodesulfurization (HDS) reactions, and hydrocracking step (HCK) where large hydrocarbons compounds are broken up and hydroisomerized simultaneously. Hydrocracking unit provides a total effluent which is submitted to atmospheric distillation. Light and middle distillates (naphtha, kerosene, diesel) are then obtained. The unconverted fraction (UCO) is upgraded to base oil using solvent or catalytic dewaxing [START_REF] Daage | Zeolites for cleaner technologies[END_REF] . As products quality depends on hydrocracking experimental conditions, process optimization requires time consuming and costly experiments. High Throughput Experimentation (HTE) units are then increasingly used but provide a too small quantity of total liquid effluent (liqtot) to perform distillation step. Furthermore, the standard measurement methods (ASTM and EN methods) of petroleum cuts properties also require relatively long time and a significant quantity of products. For example, viscosity index of base oils acquisition requires at least 106ml of UCO and 14 hours analysis time. Thus, modeling physicochemical properties of hydrocracking products is a key point. In addition, it allows better understanding of these properties.

VI measures the effect of temperature variation on the oil viscosity [START_REF] Astm D | Standard Practice for Calculating Viscosity Index from Kinematic Viscosity at 40 and 100[END_REF] . This characteristic is essential to ensure engine maintenance under different weather conditions. The measurement of VI is obtained from a logarithmic scale based on two reference oils: the former is a L series lube oil which contains a high amount of naphtenic compounds and represents a high viscosity temperature dependence; the latter is a H series lube oil with a high amount of paraffinic compounds which refers to a low viscosity temperature dependence. They are assigned a VI value of 0 and 100 respectively. The VI is then determined by interpolation or extrapolation (if a given oil sample has a VI above 100) of the tabulated H and L viscosities [START_REF] Astm D | Standard Practice for Calculating Viscosity Index from Kinematic Viscosity at 40 and 100[END_REF] . VI of base oils is particularly difficult to model since it depends on complex molecular parameters of their constituent compounds. First works [START_REF] Guo | Rheological behaviour of a new kind of fuel -residual oil, coal and water triplex synfuel[END_REF][START_REF] Evdokimov | Thermophysical properties and phase-behaviour of asphaltene-containing petroleum fluids[END_REF][START_REF] Martinez-Boza | Steady-state flow behaviour of synthetic binders[END_REF][START_REF] Reddy | A thermodynamic model for predicting n-paraffin crystallization in diesel fuels[END_REF] aimed at finding a relation between VI and temperature through thermodynamic parameters. VI was correlated to the flow of Gibbs activation energy using Arrhenius-type law [START_REF] Andrade | The structure of the atom[END_REF] , but Verdier et al. [START_REF] Verdier | A critical approach to viscosity index[END_REF] have recently shown that these relations are not always adequate for petroleum cuts because they do not include the molecular effects. The VI dependence on molecular composition was first illustrated through the project led by the American Petroleum Institute 14 . Lynch [START_REF] Lynch | Process Chemestry of Lubricant Base Stocks[END_REF] suggested to rank hydrocarbon type according to their VI value. He observed that n-paraffinic compounds have higher VI, followed by monobranched isoparaffins, then multiple branched isoparaffins and mononaphtenes and monoaromatics with long alkyl chains, etc.

Many works were then proposed in the literature to model VI from physicochemical and analytical data of petroleum cuts. Spectroscopic techniques such as [START_REF] Verdier | A critical approach to viscosity index[END_REF] C nuclear magnetic resonance (NMR) [START_REF] Adhvaryu | Quantitative NMR Spectroscopy for the Prediction of Base Oil Properties[END_REF] or near infra-red (NIR) [START_REF] Baptista | Multivariate near infrared spectroscopy models for predicting the iodine value, CFPP, kinematic viscosity at 40 °C and density at 15 °C of biodiesel[END_REF] are increasingly used to characterize petroleum cuts through multivariate analysis. Thus, Kobayashi et al. [START_REF] Kobayashi | Viscosity properties and molecular structure of lube base oil prepared from Fischer-Tropsch waxes[END_REF] correlated the VI of base oils obtained from hydrocracking of Fischer-Tropsch waxes to the ratio /

where ACN and ABN denote the average carbon number and the average branching number respectively and were estimated from [START_REF] Verdier | A critical approach to viscosity index[END_REF] C NMR spectra of base oils. Sharma et al. [START_REF] Sharma | Effects of hydroprocessing on structure and properties of base oils using NMR[END_REF][START_REF] Sharma | Influence of Chemical Structures on Low-Temperature Rheology , Oxidative Stability , and Physical Properties of Group II and III Base Oils[END_REF] Multivariate regression models were also used to predict physicochemical properties of petroleum products.

Sastry et al. [START_REF] Sastry | Determination of Physicochemical Properties and Carbon-Type Analysis of Base Oils Using Mid-IR Spectroscopy and Partial Least-Squares Regression Analysis[END_REF] developed Partial Least of Squares (PLS) models from Attenuated Reflectance Infra-Red (ATR-IR) spectra to analyze lube oils covering a wide range of VI (29 to 151). Values of 0.95 and 1.05 were obtained for and root mean square error (RMSE), respectively. Artificial Neural Network (ANN) were also applied [START_REF] Wu | Artificial neural network model to predict cold filter plugging point of blended diesel fuels[END_REF] .

Shea et al. [START_REF] Shea | Modeling Base Oil Properties using NMR Spectroscopy and Neural Networks Properties using NMR Spectroscopy and[END_REF] used structural molecular parameter estimated from [START_REF] Verdier | A critical approach to viscosity index[END_REF] C NMR data to predict VI of base oils. The model provided a very good of 0.99.

Kriging is mainly used and developed in fields of geophysics [START_REF] Dubrule | Comparing splines and kriging[END_REF][START_REF] Matheron | Traité de géostatistique appliquée[END_REF] , spatial analysis [START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Cressie | The origins of kriging[END_REF] and computer experiments [START_REF] Kleijnen | Kriging metamodeling in simulation: A review; Systemvorraussetzung: Acrobat Reader[END_REF] . To our knowledge, the use of kriging for modeling and predicting physicochemical properties in petroleum field has not been reported in the literature despite the fact that kriging might have many interesting features for this kind of applications. In particular, kriging can deal with both linear and nonlinear data structures.

Analytical data that are currently used to characterize petroleum cuts can be divided into two groups: (i) the basic physicochemical properties (density, refractive index, distillation curves, etc.) 3 and (ii) the spectral data: Whether the use of spectral data combined with multivariate regression models enables to obtain good results for the prediction of VI, basic properties remain much easier to collect and handle, despite the fact that they provide less chemical information. Furthermore, spectroscopic NMR data are not available in refinery. The use of some basic properties to model another one requires to define an adapted correlation to apply classical (non)linear regression methods (which allows direct interpretation of the coefficients). However, the results obtained are strongly affected by the choice of the nature and the number of descriptors used.

This work investigates the potential use of kriging to predict VI of base oils from basic properties of total liquid effluent. The prediction of base oil from properties of an intermediary cut is not common in the literature and constitutes an important challenge of the petroleum industry. The current paper is structured as follows. The first section details databases used in this work (base oils origins and data processing), the second section gives a theoretical approach of kriging method. A third section discusses the choice of descriptors using principal components analysis (PCA) and the quality of the developed models. This paper is focused on the once-through without recycle configuration [START_REF] Ward | Hydrocracking processes and catalysts[END_REF] .

Experimental and methods

In this paragraph data collection is described first. Details about the measurement of petroleum cuts properties are also given. A description of kriging is then proposed using a basic example. Finally, information about computational tools is given.

Hydrocracking process and Data recording

Oil samples were produced from experimental pilot plants (IFPEN, France, Solaize). An experimental scheme for data recording is given in Figure 1. Vacuum Gas Oil (VGO) is analyzed using standard methods (Table 2) before being hydrotreated using HDT reactor under various operating conditions (LHSV, pressure) and typical catalysts. The resulting total liquid effluents (HDT liqtot) are then collected and also analyzed using standard methods. One part of HDT liqtot is distillated to obtain petroleum products [START_REF]Standard Test Method for Distillation of Petroleum Products at Atmospheric Pressure. (32) NF ISO[END_REF] (Figure 1). The second part is hydrocracked using HCK reactor under various operating conditions (LHSV, temperature, pressure, Organic Nitrogen) and typical suitable different catalysts (Figure 1). The resulting total liquid effluents (HCK liqtot) are also collected and analyzed before being distillated [START_REF]Standard Test Method for Distillation of Petroleum Products at Atmospheric Pressure. (32) NF ISO[END_REF] . In both case the distillation step provides unconverted fractions (UCO) corresponding to a boiling range up to a chosen temperature of 370°C. Finally, base oils are obtained by solvent dewaxing of UCO and their VI are measured using standard NF ISO 2909 method 32 to complete databases. As mentionned above, VI of base oils acquisition remains constraining for the refiners since its measurement involves to produce at least 106 ml of UCO and 14 hours analysis time around (distillation of liqtot and solvent dewaxing included). A total of twelve different VGO feed were used to produce all these samples. The properties of these VGO are given in Table 1.

Measurement of physicochemical properties

Various physicochemical properties are commonly measured or estimated to characterize petroleum cuts in refinery. The details of standard methods and well-known correlations used in this work are given in Table 2.

Note that refers to the temperature on simulated distillation following ASTM D2887 method [START_REF] Astm D | -16A. Test Method for Boiling Range Distribution of Petroleum Fractions by Gas Chromatography[END_REF] .

, and , are the part of petroleum cuts that corresponds to a boiling range up to 370°C measured before and after the given process respectively. Thus, is the conversion of 370 + petroleum cut [START_REF] Becker | Hydrotreatment modeling for a variety of VGO feedstocks: A continuous lumping approach[END_REF] . In addition, it should be noted that only VI was measured for base oil samples. All other properties used as potential descriptors were measured on feedstock (VGO), hydrotreated liqtot or hydrocracked liqtot.

Theoretical basis

This paragraph gives a theoretical approach for kriging through a basic example.

Kriging: a spatial interpolation method

Several methods have been developed for spatial interpolation in various disciplines [START_REF] Jin | Spatial Interpolation Methods: A Review for Environmental Scientists[END_REF] . Consider the example illustrated in Figure 2 where a variable (represented by a color scale) was measured at ten different locations distributed in the two-dimensional space ( , ). The problem consists in estimating at a new location ( ) using available observations. Almost all spatial interpolation methods share the formula in Eq.1 which provides the estimation of y at the location as weighted average of the sampled data, ∑

where is the estimated value of , is the measured value at the sampled location and is the weight assigned to this sampled point. all the available data of a region of interest to derive the estimation and thus capture general trends. However, these methods require sampled locations to be homogeneously distributed. By contrast, kriging is a local approach. It operates within a small area around the point being estimated and capture local or short-range variation [START_REF] Burrough | Principles of geographical information systems; Spatial information systems[END_REF] . For kriging, the weights depend on (i) the geometrical configuration of the sampled data (spatially close data contain redundant information), (ii) their distance from the new location and (iii) the structural characteristics of the variable to estimate (anisotropy, regularity, etc.). In practice, the distribution of weights is governed by a mathematical model (we refer to section 2.4.3). Figure 3a and Figure 3b respectively represent a Gaussian and a linear distribution of weights for a two-dimensional space and an isotropic variable . The corresponding formulae are also given for each model. In these formulae, , are the coordinates of any sampled point in the space centered at the new point, , and are internal parameters of the illustrated distributions. When a Gaussian structure is used, it defines a local ellipsoidal neighborhood around the unsampled location (Figure 4a). Only the sampled data that have location inside this neighborhood affect the estimated value (i.e. 0 1,2, 3, 4, 5, 8, 9). In addition, the exponentially decrease with the square of the Euclidian distance from the unsampled location. When a linear structure is used (Figure 4b), the neighborhood is rectangular and the linearly decrease with the distance from unsampled point. In the current example (cf. Figure 4b), only sampled data 4 and 6 have non-zero weights. Thus, it is more obvious to express the kriging estimator as follows:

∑ ∈ (2)
where is the subset of sampled indexes whose weights are non-zero. Note that kriging is an exact method because it generates a prediction that is the same as the observed value at a sampled point [START_REF] Burrough | Principles of geographical information systems; Spatial information systems[END_REF] .

Statistical approach

The goal of classical regression is to approximate with a known mathematical function (linear, polynomial, etc.). Let denotes the study area; can then be expressed as in Eq. 3 at any location ∈ :

(

) 3
where f represents the deterministic part (the global trend) and is the random part (assumed to be normally distributed with zero mean and finite homogeneous variance).

A kriging approach is quite different and can be written as in Eq. 4:

(4)

where is the deterministic part of the expected value of and is a random variable such as the set , ∈ Α describes a multiGaussian process [START_REF] Stein | Gaussian approximations to conditional distributions for multi-Gaussian processes[END_REF] . According to the form of there are three different types of kriging [START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF] :

 The simple kriging (SK) where is a known constant ;

 The ordinary kriging (OK) where is an unknown constant;

 The universal kriging (UK) where ∑ is a linear combination of functions of location .

The simple kriging predictor

In case of simple kriging, , ∈ Α is assumed to describe a two-order stationary process. This stationarity assumption implies that (i) the expected value of is invariant within the study area and (ii) the covariance between and only depends on the separation vector . In what follows will denote the covariance function related to , ∈ Α such that:

a) b)
3 Results and discussions

PCA observation and descriptors research

A principal component analysis (PCA) [START_REF] Jolliffe | Principal component analysis[END_REF] was first carried out to observe HDT and HCK databases and to discuss on potential descriptors (see Table 2) relevance for the prediction of VI. The central idea of the PCA is to reduce the dimensionality of a given data set consisting of a large number of interrelated variables, while retaining as much as possible the variance contained in the data. This is achieved by transforming to a new set of variables called "principal components" (noted PCs), which are uncorrelated (constraints of orthogonality), and which are ordered so that the first few retain most of the variance in all the original data. It should be noted that the variable VI did not contribute to the calculation of principal components.

Figure 5a shows the projection of data whereas Figure 5b illustrates the scores of the original variables on the first principal components PC1 and PC2 that represent respectively 63.4% and 31.6% of the total explained variance for HDT database. The presence of three clusters that refer to different feedstock origins is observed (Figure 5a). A gradient of VI is also clearly observed. This implies that some information related to feedstock origins will be required to well-model the VI and that VI is well explained by PC1 and PC2.

As observed in Figure 5b, PC1 is correlated to the paraffinic carbon content (Cp), the API.Gravity and Kw factor (Table 2), and anticorrelated to the variables density (d) and refractive index (IR) and aromatic carbons contents (Ca). Thus, PC1 refers to information about the density of the HDT total liquid effluent analyzed. PC2 is correlated to volatility properties (TM, MeABP, etc.). Furthermore, the closest variable to the VI is the Watson characterization factor that can be understood as a ratio volatility/density. Thus, the observations above can be interpreted as follows: in case of hydrotreatment, total liquid effluents that provide base oils with highest VI are those which has high volatility and few density (i.e. low aromatic contents (Ca) [START_REF] Riazi | Characterization and properties of petroleum fractions[END_REF] ).

By opposition, total liquid effluent with high paraffinic contents provides base oils with highest VI. These observations are in accordance with the correlations proposed by Sarpal et al. [START_REF] Sarpal | Characterization by 13C n.m.r. spectroscopy of base oils produced by different processes[END_REF] and Sharma et al. [START_REF] Sharma | Effects of hydroprocessing on structure and properties of base oils using NMR[END_REF] . This dependency can be more complex since isoparaffinic and n-paraffinic, compounds that constitute all paraffinic content, can have very different rheological properties according to the branching structure [START_REF] Lynch | Process Chemestry of Lubricant Base Stocks[END_REF] . The aromatic compounds influence that tend to decline the VI of petroleum cuts was also observed by Verdier et al. [START_REF] Verdier | A critical approach to viscosity index[END_REF] and Sastry et al. [START_REF] Sastry | Determination of Physicochemical Properties and Carbon-Type Analysis of Base Oils Using Mid-IR Spectroscopy and Partial Least-Squares Regression Analysis[END_REF] .

The projection of data and the scores of original variables on PC1 (80% of the total variance explained) and PC2 (11% of total variance explained) for HCK database are represented in Figure 6a and Figure 6b respectively.

The results are quite different that in the previous case. Indeed, no general trend can be observed in Figure 6a for the sample distribution regarding to VI. Two data points are relatively far from the others due to their particularly high paraffinic contents. They illustrate the importance to make the distinction between n-paraffinic and isoparaffinic compounds for modelling VI. PC1 is now correlated to density (d), refractive index (IR), naphtenic carbons contents (Cn) and information about volatility (TM, MeABP, etc.) and anticorrelated to Ca, Cp, the VI of feedstock (VI.chg) and the conversion rate ( ) (Figure 6b). Although VI is not as well explained as in hydrotreatment, it is correlated to PC1 and slightly anticorrelated to PC2. These observations can be interpreted as follows: the simultaneous decrease of aromatics and n-paraffins contents in total liquid effluent, due to the conversion of aromatics into naphtenes and the isomerization of n-paraffins during the hydrocracking process, reduces their effect on density and volatility. Consequently, both properties are now related to the naphtenes and isoparaffins contents. The variable the most correlated to VI is now the VI of feedstock (VI.chg).

That implies that feedstock characteristics significantly affect the VI of base oils. However, the role of paraffinic compounds is more difficult to determine for HCK database since the variable Cp does not allow to discriminate normal and iso paraffinic carbons. Overall, the presence of paraffinic carbons in total effluent improve the VI of base oil whereas naphtenic carbons (Cn) tend to decrease it. These observations are in total accordance with Sarpal et al. [START_REF] Sarpal | Characterization by 13C n.m.r. spectroscopy of base oils produced by different processes[END_REF] . 2)

Models comparison

A variable selection was performed to select final descriptors among all variables used in the PCA analysis. This step is necessary as some of these variables are strongly related, that results in redundant information. Thus, it may be not relevant to use all of them in a regression model. Linear and kriging models were developed using same descriptors in each case. Obtained models were then compared by evaluating performances on training and validation data sets respectively.

For HDT database, parity plots obtained by applying a leave-one-out cross validation to the training set are overlapped for linear and kriging model in Figure 7a: the points predicted using MLR are represented by red stars whereas those predicted with kriging are illustrated by black circles. Both models seem provide very close prediction for the leave-one-out step. This observation is supported by RMSECV values of 2.4 and 1.9 for linear and kriging model respectively (Table 4). Similar observations can be done by considering validation set.

Indeed, the parity plots of linear and kriging represented in Figure 7b are still very close. RMSEP values of respectively 2.4 and 2.1 (Table 4) for linear and kriging models were obtained. Furthermore, the percentages of points at a precision level under the confidence interval (IC) of the standard method (±2 points of VI) are equivalent (63.41%) for both models. All these previous observations show that kriging is slightly better than linear model in case of HDT database. Globally, both models are well-performing since their precision are close to the IC of the measurement method. The fact that linear model provides good performances is consistent with the gradient of VI observed in the previous PCA analysis (Figure 5a).

For HCK database, significant differences can be noted between linear and kriging. Firstly, the parity plots obtained by applying leave-one-out to training data set (Figure 8a) and using model on validation data a b

(Figure 8b) seem globally less performing for linear model than for kriging one. Furthermore, a RMSECV value of 4.9 was obtained for RLM whereas kriging provided a value of 3.5 (Table 4). Secondly, statistics on validation set are also clearly better for kriging than for MLR. Indeed, a RMSEP value of 2.6 was obtained for kriging against 3.6 for MLR (Table 4). Kriging also provides a better percentage of points at a precision level under the confidence interval (IC) of the standard method (62.5% against 46.87%). The observations above show that kriging has clearly better performances than linear model in case of HCK database.

Conclusion

Kriging is currently used in geostatistics for modelling field values in localized geographical regions. The dimension is generally less than 3 in this field (spatial coordinates and altitude). The use of kriging in situations where the number of descriptors increases remains problematic as the covariance structure may be too complex to define. This might explain why no example can be found in the literature attempting the modelling of physicochemical data.

In this work, both kriging and linear models were developed to predict VI of base oils produced from hydrotreatment and hydrocracking. In the first case of base oils obtained from hydrotreatment where linear modeling provided good performances in accordance with the reproducibility of the standard measurement method, kriging modeling provided similar performances. For the second case of base oils produced by hydrocracking, MLR performances were insufficient whereas kriging modelling improved performances with identical input variables. The divergence in the MLR model performances according to the considered HDT or HCK database is mainly due to the impact of catalysts which clearly differs according to the referred process.

Indeed, catalysts used for hydrotreatment which are generally made of metals such as Molybdenum or Tungsten related to Cobalt or Nickel (NiMo, NiW, etc.) less affect the structure of molecular compounds than those of hydrocracking (presence of zeolite). Globally, the obtained results are particularly important since the descriptors used are not referring to the final product (base oil) but are measured from both feedstock and on the total liquid effluent before the distillation and dewaxing steps. Thus, kriging allows to access to a very important base oil characteristic without performing distillation and solvent dewaxing steps in case of hydrotreatment as well as in a b

hydrocracking case. To our knowledge, this is the first model for the prediction of VI of base oil using kriging but also for any petroleum products properties and especially from the total effluent global properties. This breakthrough induces an important gain concerning time analysis and sample volume consumption.

Another advantage for kriging model is that it provides a measure of prediction uncertainties as it uses stochastic approach [START_REF] Wang | Kriging regression in digital image correlation for error reduction and uncertainty quantification[END_REF] . Kriging weights estimation are obtained by solving a linear system [START_REF] Dubrule | Comparing splines and kriging[END_REF] , This method would be faced with similar problem in the case of linear model when the number of descriptors will be more significant (as for instance considering spectral data). It should be interesting to adapt this approach to high multivariate situations ( ≫ ).

Figure 1 . 2 Databases

 12 Figure 1. Experimental scheme of HCK process and for data recording

Figure 2 .

 2 Figure 2. Repartition of sampled data for basic example in a two-dimensional case; colored points = sampled data; 0 = new point

Figure 3 .

 3 Figure 3. Example of local structures for weights repartition; a) -Gaussian structure ; b) -linear model; green zone delimited area of influence around any unsampled location

Figure 4 .

 4 Figure 4. a -Weights distribution when using Gaussian structure ; b -Weights distribution using linear structure

Figure 5 .Figure 6 .

 56 Figure 5. a -Projection of HDT data samples on first factorial plane; circles: training samples; squares: validation samples; color scale: from blue (low VI) to red (high VI); b -Correlation circle of input variables in the first factorial plane (see Table2)

Figure 7 .Figure 8 .

 78 Figure 7. Parity plots for database 1 (HDT); a -On training set (leave-one-out); b -On validation set

  and Sarpal et. al[START_REF] Sarpal | Correlation of structure and properties of groups I to III base oils[END_REF] used a linear relationship to model VI of base oils from molecular parameters (different types of carbon contents,[START_REF] Verdier | A critical approach to viscosity index[END_REF] C NMR peaks area corresponding to different branching structures). These works provided important information about the VI dependence on molecular structure of petroleum cuts, but due to the limited number of samples (6 to 8), model robustness can be questioned. Verdier et al.[START_REF] Verdier | A critical approach to viscosity index[END_REF] used a similar approach to model VI of a set of VGOs, base

	oils obtained from hydrotreated VGOs and base oils produced from hydrocracked VGOs. An average absolute
	deviation of 3 VI points and a	of 0.96 were obtained for the 20 oils studied. This study has drawn attention
	because of the variety of samples analyzed, deriving from both hydrotreatment and hydrocracking processes
	covering a wide range of VI (-38 to 146).

Table 1 : Properties of the twelve different VGO used to produce samples

 1

	VGO	Origin	Density at 15°C (g/cm 3 )	Nitrogen (wt ppm)	Sulphur (wt%)	IBF -FBF (°C)	Kinematic viscosity at 70°C	Kinematic viscosity at 100°C
	A	Taching	0,8830	1015,00	0,1255	320,90 -622,60	15,97	7,51
		(chinoise)
		SR
	B	Gansu (Chine)	0,8974	1190,00	1,0409	348,00 -584,00	82,25	12,55
	D	Iranian Saniya	0,9375	1300,00	2,8743	321,90 -615,10	124,74	11,11

Table 2 :

 2 Measured and evaluated properties of feedstock, total liquid effluents and base oil samples

		Property		Standard Methods	References	Range of value
								HDT liqtot	HCK liqtot
		Viscosity Index (VI) Density (d) (in g.cm -3 ) Refractive Index (IR) Simulated distillation	ASTM D2270 or NF ISO 2909 ASTM D7042 ASTM D1218 ASTM D2887	7 / 32 35 36 33	0.85 -0.92 1,44 -1,49	0.78 -0.87 1.41 -1.46
	( ,	5, 10, 15, … , 90, 95
		(in °C) Carbon type		ASTM D3238		37
	characterization (Ca, Cp,
	Cn) (in %) Number of aromatic rings	ASTM D3238		37
	(Ra) Molecular Mass (M) (in g.mol -1 ) Mean Average Boiling	IFP 0010350 API		Based on 3 3	266 -450	156 -366
	Point (MeABP) Watson characterization factor (Kw) Average Temperature (TM)	.	.			3 38	325 -465	203 -385
		(in °C) Average Weighted		2	3	4	38	393 -506	267 -441
	Temperature (TMP) (in °C) Specific gravity (Spgr) API Gravity	.	7 1.001 141.5 131.5	3 3
	Conversion rate	(in %)	,			,
					,

Table 4 :

 4 Summary of statistics for linear and kriging models

	Data HDT	Model Linear	RMSECV 2.4	RMSEP 2.4	τ ±IC (%) 63.4	τ ±2IC (%) 92.4
	Database
		Kriging	1.9	2.1	63.4	95.1
	HCK	Linear	4.9	3.5	46.9	75
	Database
		Kriging	3.6	2.6	62.5	90.5

-6620 -97

, (5) Various types of covariance structure are available in the literature [START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF][START_REF] Isaaks | Applied geostatistics[END_REF] .

Let

denotes a new location and the value to predict. Note the simple kriging estimator of .

is defined as an optimal linear estimator (i.e. unbiased with minimal variance). The following linear system was obtained by combining this definition to the stationary multiGaussian assumptions:

Thus, the simple kriging estimator corresponds to the linear estimator of which has weights solution of the system. Note that optimal weights directly depends on the covariance structure.

This probabilistic approach of kriging is a key point since it allows to consider various approaches in modelling of prediction errors including homogeneous (homoscedasticity) and heterogenous (heteroscdeasticity) situations .

Furthermore, although interpolating functions are usually froced to pass through training data, a "nugget effect" referring to an effective variance related to measured values may be introduced in the covariance structure. samples to the validation set. In all cases, data were divided by performing a space filling algorithm in order to ensure that the most part of predicted points are in situation of interpolation [START_REF] Dicedesign | Two R Packages for Design and Analysis of Computer Experiments[END_REF] .

Assessment of models quality

Model performances are commonly evaluated by estimating well-known statistics on training and validation sets separately. Since kriging is an exact interpolation method, there is no modelling error at training locations. Thus, a leave-one-out (LOO) error estimation method was used to compare models on training set.

The leave-one-out is a special case of cross-validation where the number of folds equals the number of instances in the data set. The learning algorithm is applied once for each instance, using all other instances as training set and using the selected instance as a single-item test set [START_REF] Mang | Encyclopedia of Lubricants and Lubrication[END_REF] . Each predicted value is then compared with the measured one. The statistics used to evaluate models performances on training and validation test are given in Table 3.

Variable selection

A variable selection model algorithm was applied to determine suitable final descriptors for modelling VI among all properties presented in Table 2. These properties were chosen by hydrocracking expertise, then they are all relevant. However, they may contain redundant information due to strong correlations. Thus, it may be not relevant to use all of them to develop a linear model. The variable selection algorithm allows to select final descriptors among all available assuming that the model is linear. The best subset of variables is the one that maximizes a predefined criterion. In this work, we used the adjusted correlation coefficient ( ) which is a penalized release of the classical 47 .

Used software

All models have been implemented using statistical software R (version 3.


Lei Mo
email: lei.mo@inria.fr

Angeliki Kritikakou

Olivier Sentieys
email: olivier.sentieys@irisa.fr

Approximation-aware Task Deployment on Asymmetric Multicore Processors

Asymmetric Multicore Processors (AMP) are a very promising architecture to deal efficiently with the wide diversity of applications. In real-time application domains, in-time approximated results are preferred than accurate -but too late -results. In this work, we propose a deployment approach that exploits the heterogeneity provided by AMP architectures and the approximation tolerance provided by the applications, so as to increase as much as possible the quality of the results under given energy and timing constraints. Initially, an optimal approach is proposed based on problem linearization and decomposition. Then, a heuristic approach is developed based on iteration relaxation of the optimal version. The obtained results show 16.3% reduction in the computation time for the optimal approach compared to the conventional optimal approaches. The proposed heuristic approach is about 100 times faster at the cost of a 29.8% QoS degradation in comparison with the optimal solution.

I. INTRODUCTION

Due to the wide range of nowaday's application domains, such as multimedia, encryption and network, the applications require different system resources, impose different design constraints and have different optimization objectives. Conventional Symmetric Multi-core Processors (SMP) have become a less promising option, as they cannot simultaneously meet theseoften conflicting -requirements and goals. An overprovisioned core wastes area and energy, whereas an underprovisioned core fails to provide the required performance [START_REF] Sentieys | A survey of techniques for architecting and managing asymmetric multicore processors[END_REF]. Thus, Asymmetric Multi-core Processors (AMP) are proposed to deal with the application diversity. AMPs provide the required heterogeneity, since they consist of cores that differ in microarchitectural features, such as pipeline design, issue/fetch width, and cache hierarchy -and not merely in frequency/voltage, as SMPs [START_REF] Sentieys | A survey of techniques for architecting and managing asymmetric multicore processors[END_REF].

In real-time application domains, approximated results obtained in time are preferred than accurate results obtained after the deadline. For instance, in a real-time video application, at each period, an imperfect, but acceptable, quality image is initially produced from the received data. Then, this image can be further refined depending on the available system resources [START_REF] Aydin | Optimal rewardbased scheduling for periodic real-time tasks[END_REF]. Similar applications can be found in many other areas like mobile target tracking, real-time heuristic research and control engineering [START_REF] Yu | Quality-driven dynamic scheduling for real-time adaptive applications on multiprocessor systems[END_REF], [START_REF] Mendez-Diaz | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF]. In order to maximize the quality of the obtained application results and at the same time meet the design constraints, proper deployment approaches have to explore both the architecture heterogeneity and the approximation tolerance of the applications.

The majority of the deployment approaches on multicore processors (SMP/AMP) focuses on precise computation tasks [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF]- [START_REF] Mahmood | Energy-aware real-time task scheduling in multiprocessor systems using a hybrid genetic algorithm[END_REF]. Usually these approaches aim to minimize the energy consumption under real-time constraints. On the other hand, Indep. Alloc.

Migr.

[14] √ √ √ √

the deployment approaches that focus on tasks that tolerate approximation, such as Imprecise Computation (IC) tasks [START_REF] Liu | Imprecise computations[END_REF], aim to maximize the Quality of Service (QoS) under energy and real-time constraints. The state-of-the-art on IC task deployment approaches can be classified based on whether: 1) the tasks are independent or dependent, 2) the platform is an SMP or an AMP, and 3) the solution is optimal or heuristic. Table I positions our work with respect to several representative papers from the literature. Most existing works deal with the independent IC tasks deployment problems on AMP/SMP. Works on deploying dependent IC tasks on AMP/SMP are rare and they solve different problems compared to our work, e.g., task-to-processor allocation is fixed and given upfront [START_REF] Yu | Quality-driven dynamic scheduling for real-time adaptive applications on multiprocessor systems[END_REF], [START_REF] Yu | Dynamic scheduling of imprecise-computation tasks on real-time embedded multiprocessors[END_REF], and no task migration is considered in [START_REF] Mo | Controllable QoS for imprecise computation tasks on DVFS multicores with time and energy constraints[END_REF]. The extension of the optimal deployment methods from dependent to independent tasks and from SMPs to AMPs is not straightforward, as additional nonlinear constraints are usually introduced into the problem. On the other hand, the heuristics of previous works usually adopt multi-step optimizations, e.g., [START_REF] Mendez-Diaz | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF], [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in realtime MPSoC systems[END_REF], [START_REF] Mo | Energy-quality-time optimized task mapping on DVFS-enabled multicores[END_REF] and they require that the variables can be solved separately. Hence, they need to be re-designed when the problem formulation changes.

In contrast, the heuristic proposed in this paper can solve general Mixed-Integer Linear Programming (MILP) problem. The goal of this work is to address the problem of dependent IC tasks deployment on AMP, so as to increase the QoS of the obtained results and at the same time meet the timing and energy constraints. Our first contribution is the formulation of the deployment problem as a Mixed-Integer Non-Linear Programming (MINLP) and its safe transformation into an MILP. Our second contribution is an Optimal Deployment Approach (ODA) based on Benders decomposition [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF] and the algorithm's convergence analysis. Our third contribution is a Heuristic Deployment Approach (HDA), which relaxes the number of required iterations between the ODA subproblems. Finally, the extended experimental results show 16.3% less computation time for the ODA compared to the conventional state-of-the-art optimal approaches. HDA runs ∼100 times faster than ODA with a cost of 29.8% in QoS.

The rest of the paper is organized as follows. Section II introduces the system model. Section III presents the problem formulation. Section IV provides the problem linearization method. Section V and Section VI propose the ODA and HDA algorithms, respectively. Section VII evaluates the performance of the proposed methods. Section VIII concludes this study. II. SYSTEM MODEL 1) Task Model: Definition 2.1: (IC task [START_REF] Liu | Imprecise computations[END_REF]) A task can be logically decomposed into a mandatory subtask and an optional subtask. The mandatory subtask must be completed before the deadline to produce the acceptable result, whereas the optional subtask can be left incomplete at the cost of reduced quality.

We consider a task set T consisting of N dependent and periodic IC tasks {τ 1 , . . . , τ N } released at time 0 and sharing a common hyper-period H. The task set T is modeled by a Directed Acyclic Graph (DAG) G(V, E), where the nodes V correspond to the IC tasks and the edges E indicate the data dependencies between the tasks. Each task τ l is described by a tuple {o l , M l , O l , t s l , d l , g l }. A task τ l is decomposed into a mandatory subtask and an optional subtask characterized in number of execution cycles, i.e. M l and o l , respectively. O l is the maximum optional cycles of task τ l . M l and O l are measured in Worst Case Execution Cycles (WCECs). t s l , d l and g l are the start time, the deadline and the period of task τ l , respectively. We assume that the tasks are non-preemptive and H is the least common multiple of {g 1 , . . . , g N }.

2) Platform Model: We consider an AMP platform with M processors {θ 1 , . . . , θ M }. Each processor θ k is characterized by a given supply voltage and frequency pair (v k , f k ). According to the characteristics of many AMP platforms (e.g., ARM big.LITTLE) [START_REF] Chwa | Energy and feasibility optimal scheduling on big.LITTLE platforms[END_REF], we introduce the concept of cluster. M processors are divided into R clusters {W 1 , . . . , W R }. Each cluster W r consists of a set of symmetric processors that have the same frequency characteristics (e.g., minimal, maximal and operating frequencies). Therefore,

v k = v k and f k = f k , if k, k ∈ W r .
In addition, we define λ l,k ∈ (0, 1] as the execution efficiency of processor θ k , when it executes task τ l [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF]. Correspondingly, the Worst Case Execution Time (WCET) of task τ l , when it is executed on processor θ k , is calculated as M l +o l f k λ l,k . 3) Energy Model: Processors can operate in two modes: idle and active. A processor is said to be in the active mode, if the processor currently executes a task. Otherwise, it is in the idle mode. The power consumption of a processor θ k is expressed as P c k = P s k + P d k , where P s k = C s k v ρ k k is the static power of the processor ready to execute (being either on the active or idle mode),

P d k = C d k f k v 2 k
is the dynamic power of task execution, and C s k , ρ k and C d k are constants depending on processor type [START_REF] Chwa | Energy and feasibility optimal scheduling on big.LITTLE platforms[END_REF]. We assume that when a processor has no task to execute, it goes into idle mode. The transition time and energy overhead is considered very small compared to the one required to execute a task and is assumed to be incorporated into the execution time and energy of the task [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF]. This power consumption model is widely adopted by existing works [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF], [START_REF] Chen | Energy optimization for real-time multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF], [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in realtime MPSoC systems[END_REF]. The system is energy-constrained [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in realtime MPSoC systems[END_REF] in the sense that: 1) the energy budget E s is fixed and cannot be replenished during the hyper-period H, and 2) with E s budget, all mandatory subtasks {M 1 , . . . , M N } are ensured to finish, but not all optional subtasks {O 1 , . . . , O N } can complete their executions.

' 1  ' 2  ' 3  ' 4  ' 5  ' 6  ' 7  ' 8  1  2  3  4 

III. PROBLEM FORMULATION

Let N {1, . . . , N }, N {1, . . . , 2N }, M {1, . . . , M } and R {1, . . . , R}. Denote τ 2l-1 and τ 2l as the mandatory and the optional subtasks of task τ l (∀l ∈ N ), respectively. Based on the task and the platform models mentioned above, we consider 1) task level migration [START_REF] Chwa | Energy and feasibility optimal scheduling on big.LITTLE platforms[END_REF] on the AMP (i.e., task τ l in one cluster can migrate to the other cluster after the execution of the mandatory subtask, and 2) subtask τ 2l starts running only after subtask τ 2l-1 is finished. Note that the subtask set {τ 2l-1 , τ 2l } (∀l ∈ N ) can be rewritten as {τ i } (∀i ∈ N ). The dependency between the subtasks {τ i } in one hyper-period H are described by a symmetric Subtask Dependency Decision (SDD) matrix s [s i,j ]. If s i,j = 1, subtasks τ i and τ j are dependent with each other, otherwise, s i,j = 0. In addition, we introduce an Execution Order Decision (EOD) matrix p [p i,j ]. If p i,j = 1, subtask τ i precedes subtask τ j and τ j is the closest task of τ i , otherwise, p i,j = 0. Fig. 1 shows the DAG of an illustration example with 4 dependent IC tasks. The corresponding SDD matrix s and EOD matrix p are as follows: s =    0 1 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 1 0    , p =    0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0

   .
The problem consists of an objective function to maximize the system QoS subject to a set of real-time and energy constraints. We have to simultaneously cope with the constraints of 1) task level mitigation, 2) subtask dependency, 3) task deadline, and 4) energy supply. Under these constraints, we determine 1) on which processor should the subtask be executed (subtask allocation and mitigation), and 2) when should the subtask start and end (subtask scheduling and optional subtask adjustment). To formulate the task deployment problem, we introduce the following variables: 1) define q [q i,k ] and let the binary variable q i,k = 1, if subtask τ i is assigned to processor θ k , otherwise, q i,k = 0; 2) define u [u i,j ] and let the binary variable u i,j = 1, if subtask τ i precedes subtask τ j (i.e., τ i starts after the end time of τ i ), otherwise, u i,j = 0; 3) define t s [t s i ] and let the continuous variable t s i be the start time of subtask τ i ; 4) define o [o l ] and let the continuous variable o l be the optional subtask of task τ l . It has been shown in [START_REF] Aydin | Optimal rewardbased scheduling for periodic real-time tasks[END_REF], [START_REF] Yu | Quality-driven dynamic scheduling for real-time adaptive applications on multiprocessor systems[END_REF] that in several application domains (e.g., image processing), a linear function models the cases where the system QoS increases uniformly during the optional execution. We introduce a linear function g l (o l ) for task τ l , and our aim is to maximize QoS function l∈N g l (o l ). The Primal Problem (PP) is given by

PP : min q,u,o,t s - l∈N g l (o l ), (1)
s.t. C 1 -C 7 .
The constraint descriptions are as follows:

• Subtask allocation constraint C 1 : k∈M q i,k = 1 (∀i ∈ N ). Each subtask τ i is executed on only one processor.

• Task migration constraint C 2 : q 2l-1,k + q 2l,k ≤ 1 (∀l ∈ N , ∀k, k ∈ W r , k = k , ∀r ∈ R).
The migration of task τ l is performed among different clusters, as the processors in the same cluster are homogeneous. Similar to the previously published work [START_REF] Chwa | Energy and feasibility optimal scheduling on big.LITTLE platforms[END_REF], we consider task migration occurs at no cost or penalty.

• Subtask non-preemptive constraints C 3 : t e i ≤ t s j +(2-q i,k - q j,k )H + (1 -u i,j )H (∀i, j ∈ N , s ij = 0, i = j, ∀k ∈ M) and C 4 : t e j ≤ t s i + (2 -q i,k -q j,k )H + u i,j H (∀i, j ∈ N , s ij = 0, i = j, ∀k ∈ M). In C 3 and C 4 , t e 2l-1 = t s 2l-1 + k∈M q 2l-1,k M l f k λ 2l-1,k and t e 2l = t s 2l + k∈M q 2l,k o l f k λ 2l,k (∀l ∈ N
) are the end time of subtasks τ 2l-1 and τ 2l , respectively. If subtasks τ i and τ j are executed on the same processor θ k (i.e., q i,k = q j,k = 1), C 3 and C 4 are meaningful, else, C 3 and C 4 are always satisfied. With q i,k = q j,k = 1, if subtask τ i precedes subtask τ j (i.e., u i,j = 1), we have t e i ≤ t s j and t e j ≤ t s i + H; if u i,j = 0, we have t e i ≤ t s j + H and t e j ≤ t s i .

• Subtask dependency constraint C 5 : t s j + (1 -p i,j )H ≥ t s i + p i,j (t e i -t s i ) (∀i, j ∈ N , i = j). If p ij = 1
, subtask τ i precedes subtask τ j and τ j is the closest task of τ i . C 5 is meaningful for task start time t s i and t s j (i.e., t s j ≥ t e i ). If

p ij = 0, C 5 can be ignored since the inequality t s j + H ≥ t s i is always true. • Task deadline constraint C 6 : t e 2l = t s 2l + k∈M q 2l,k o l f k λ 2l,k ≤ d l (∀l ∈ N ). Each task τ l should be finished before its prede- fined deadline d l (0 ≤ d l ≤ H). • Energy constraint C 7 : k∈M [t d k P c k + (H -t d k )P s k ] ≤ E s .
The total energy consumed by M processors during the hyper-period H should not exceed the energy supply E s . Let

t d k = l∈N q 2l-1,k M l f k λ 2l-1,k + q 2l,k o l f k λ 2l,k
be the time required to execute all the tasks assigned to processor θ k . Therefore, H -t d k is the time of processor θ k being in the idle state during the hyper-period H. For tractability reasons, we consider o i as continuous variables. When PP is solved, we round the result down. This impact of one cycle is negligible, since the tasks usually execute typically hundreds of thousands of cycles [START_REF] Cortes | Quasi-static assignment of voltages and optional cycles in imprecise-computation systems with energy considerations[END_REF].

IV. MINLP LINEARIZATION

Due to the product q 2l,k o l of the optimization variables in the constraints C 2 -C 7 , PP is an MINLP problem, which is N Phard. Due to page limit, the proofs of the lemmas are omitted.

Lemma 4.1: Given constants s 1 , s 2 > 0 and two constraint spaces

S 1 = {[h, b, x]|h = bx, -s 1 ≤ x ≤ s 2 , b ∈ {0, 1}} and S 2 = {[h, b, x]| -bs 1 ≤ h ≤ bs 2 , h + bs 1 -x -s 1 ≤ 0, h -bs 2 -x + s 2 ≥ 0, b ∈ {0, 1}}, then S 1 S 2 . Since q 2l,k ∈ {0, 1} and 0 ≤ o l ≤ O l , based on Lemma 4.1,
an auxiliary variable (continuous) h 2l,k is introduced into PP to replace the nonlinear item q 2l,k o l . Thus, PP is rewritten as

PP1 : min q,u,o,h,t s - l∈N g l (o l ), (2)
s.t. C 1 -C 10 ,
where

t e 2l = t s 2l + k∈M h 2l,k f k λ 2l,k and t d k = l∈N q 2l-1,k M l f k λ 2l-1,k + h 2l,k f k λ 2l,k . C 8 : h 2l,k ≤ q 2l,k O l , C 9 : h 2l,k ≤ o l and C 10 : h 2l,k + (1 -q 2l,k )O l ≥ o l (∀l ∈ N , ∀k ∈ M)
are the additional constraints introduced by the MINLP linearization. The PP1 is an MILP problem, as the binary variables (q, u) and the continuous variables (o, h, t s ) are coupled with each other linearly, which is much easier to solve than the PP. Remark 4.1: Since 1) the objective functions of PP and PP1 are the same, and 2) the variable replacement does not change the feasible region of the problem (S 1 S 2 ), solving PP1 is equivalent to solving PP, i.e., the optimal values of the objective functions of PP and PP1 are the same.

V. OPTIMAL TASK DEPLOYMENT ALGORITHM

In this section, we propose an optimal algorithm ODA to efficiently solve PP1. Using concise notions, the objective function and constraints of PP1 can be rewritten as Lemma 5.1: Solving the MP (SP), we obtain a lower bound Φ l (k) (upper bound Φ u (k)) of the optimal value of PP2, Φ * .

PP2 : min x,y Φ = f T y (3)
2) MP and SP Iterations: At the initial iteration, i.e., k = 0, we set {x(0)|Ax(0) b 1 }, Φ l (0) = -∞ and Φ u (0) = +∞. Denote ε as a small positive tolerance. The iteration stops when Proof: In Case 2, since 1) Φ(k) < λ(k) T (Cx(k) -b 2 ) due to Lemma 4.1, where Φ(k) is the optimal solution of MP at the k th iteration, and 2) x(k) is not the optimal solution of PP2 due to Φ u (k) -Φ l (k) > ε, the non-optimal solution x(k) is excluded by the constraint Φ ≥ λ(k) T (Cx -b 2 ). In Case 3, the SP has no feasible solution under the given solution x(k). This problem is feasible if the positive variables ξ [ξ i ] are introduced to relax the constraints. We construct a Feasibility Check Problem (FCP) and solve its dual problem (DFCP): where λ

Φ u (k) -Φ l (k) ≤ ε.
FCP : min y,ξ 0 1 T ξ (7)
[ λi ] are the Lagrange multipliers. Let ξ(k) and λ(k) be the optimal solutions of FCP and DFCP at the k th iteration, respectively. If the SP is infeasible, this implies that some constraints cannot be satisfied, and, thus, their related relaxation variables are non-zero. We have 1 T ξ(k) > 0, and further, 1 T ξ(k) = λ(k) T (Cx(k) -b 2 ) > 0 due to the strong duality. Therefore, the infeasible solution x(k) can be excluded by the constraint 0 ≥ λ(k) T (Cx -b 2 ).

Lemma 5.2: The lower bound sequence {Φ l (k)} is increasing and the upper bound sequence {Φ u (k)} is decreasing.

Theorem 5.2: At each iteration with a new FC or IC added into the MP, the solution obtained by ODA converges to the global optimal value Φ * within a finite number of iterations.

Proof: Based on Theorem 5.1 and Lemma 5.2, as well as the fact that the dimension of binary variables x is finite, the solution converges to the global optimal value within a finite number of iterations.

Remark 5.2: The reasons why DSP is solved instead of SP are that 1) the two problems are equivalent due to the strong duality, and 2) the FCs and the ICs are constructed by the solution of the DSP. The reasons why DFCB is solved instead of FCB are that 1) the two problems are equivalent due to the strong duality, and 2) λT (Cx -b 2 ) is a function with respect to variables x, but not 1 T ξ, i.e., 0 ≥ 1 T ξ is an invalid constraint for the MP.

VI. HEURISTIC TASK DEPLOYMENT ALGORITHM Although the solution found by ODA is optimal, this method cannot be used to efficiently solve large problem instances. This is because: 1) the MP is an ILP, which is hard to solve directly compared to the LP-based SP, and 2) the number of MP constraints increases with the number of iterations, due to the addition of FCs and ICs. In order to circumvent these difficulties, we propose a heuristic method. The structure of PP/PP1 shows that 1) the problem has a large number of binary and continuous variables, and 2) these variables are highly coupled with each other. Thus, it is difficult to design a heuristic based on the traditional multi-step optimization methods, such as [START_REF] Mendez-Diaz | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF], [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in realtime MPSoC systems[END_REF]. These methods usually fix one type of variables to determine the other type of variables. Using this multi-step method, the problem in the latter steps maybe unfeasible due to the solution provided by the previous steps. In contrast, inspired by the structure of ODA, we propose a novel heuristic algorithm HDA to efficiently solve PP2.

Theorem 6.1: The FC and IC generated by solving the DSP with x (k) do not exclude the optimal solution (x * , y * ) of PP2, where x (k) is an arbitrary feasible solution of the MP.

Proof: If the DSP has a bounded solution λ (k) with x (k), a new FC is generated, i.e., Φ ≥ λ (k) T (Cx-b 2 ). On the other hand, if the DSP has an unbounded solution with x (k), a new IC is generated, i.e., 0 ≥ λ (k) T (Cx -b 2 ), where λ (k) is the solution of DFCP with x (k). Next, we prove that the solution (x * , y * ) does not violate the above FC and IC. If the DSP has a bounded solution λ (k) with x (k), suppose that x * and Φ * = f T y * violate the FC, i.e., Φ * < λ (k) T (Cx * -b 2 ). This contradicts the fact that Φ * is the optimal value of the objective function of the DSP with x * (i.e., Φ * = max λ 0 λ(k) T (Cx *b 2 ) ≥ λ (k) T (Cx * -b 2 ). Hence, the FC will not exclude the optimal solution x * . If the DSP has an unbounded solution with x (k), x (k) is excluded by the IC. Since x * = x (k), x * does not violate the IC.

Since MP is an ILP, the feasible solution x (k) can be found through heuristics, such as the Feasibility Pump (FP) method [START_REF] Fischetti | The feasibility pump[END_REF]. The HDA structure is similar to the ODA structure except that 1) the optimal MP solution x(k) is replaced by a feasible MP solution x (k) during the MP and SP iteration, and 2) the iteration stops, when the SP has a bounded solution for first time (Fig. 2).

s k = 1.478, C d k = 0.471, ρ k = 0.379 C s k = 1.191, C d k = 0.153, ρ k = 0.757
Compared with the existing heuristics, HDA is able to solve general MILP problem, even if the problem variables cannot be calculated separately. In addition, we know how far the feasible solution found by HDA from the optimal one.

VII. SIMULATION RESULTS

The AMP platform used in the experiments consists of four Cortex-A15 ("big") along with four Cortex-A7 ("LITTLE") cores [START_REF] Chwa | Energy and feasibility optimal scheduling on big.LITTLE platforms[END_REF]. Table II summarizes the AMP parameters obtained from [START_REF] Chwa | Energy and feasibility optimal scheduling on big.LITTLE platforms[END_REF]. The parameters of the experimental set-up are depicted in Table III. The ranges of the mandatory cycles M l and the maximum optional cycles O l are obtained from the work of [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in realtime MPSoC systems[END_REF], [START_REF] Jayaseelan | Temperature aware task sequencing and voltage scaling[END_REF], where analytic experiments took place to measure the range in cycles for MiBench and MediaBench benchmarks. To set the task deadlines we introduce r l and ts l , i.e. the relative deadline and the temporary start time of task τ l , respectively. If task τ i precedes task τ j and τ j is the closest task of τ i , the temporary end time of τ i is set to the temporary start time of τ j (i.e., te i = ts i + r i = ts j ). If τ l is the first task in one hyper-period, ts l = 0. r l and r l are the minimum and the maximum time required to execute the precise version of the tasks with {M l + O l } (∀l ∈ N ) cycles, respectively. E h is the average minimum energy required to execute the precise version of the tasks with {M l + O l } (∀l ∈ N ) cycles. η is an energy efficiency factor. Not all the optional subtasks can be fully executed, i.e., {O l } (∀l ∈ N ), as this violates the real-time and energy constraints. The QoS function of the tasks is adopted from [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in realtime MPSoC systems[END_REF], [START_REF] Mo | Energy-quality-time optimized task mapping on DVFS-enabled multicores[END_REF]. The use of arbitrary task graphs does not affect our problem formulation, as different task graphs only change the values of the problem parameters {A, C, D, f , b 1 , b 2 }. The simulations are performed on a laptop with quad-core 2.5-GHz Intel i7 processor and 16-GB RAM, and the algorithms are implemented in MATLAB 2016a.

d l = ts l + r l H = max l {d l } r l ∈ [r l , r l ] r l = min k M l +O l f k λ l,k r l = max k M l +O l f k λ l,k

Objective function

Energy supply

N l=1 α l o l , 0 ≤ α l ≤ 1, N l=1 α l = 1 Es = ηE h E h = M H -N l=1 min ∀k M l +O l f k λ l,k M k=1 P s k M + N l=1 min ∀k M l +O l f k λ l,k P s k + P d k Fixed parameters Tuned parameters R M Min/Max/Step N η 2 8
10/50/10 0.8/0.9/0.1

Initially, we compare the performance (QoS and computation time) of the ODA and the HDA with: 1) optimal approach -Branch and Bound method (B&B) [START_REF] Boyd | Branch and bound methods[END_REF], known to provide optimal solution for the MILP problem, and 2) stochastic approach -Genetic Algorithm (GA), which is provided by the MATLAB optimization toolbox.

The QoS achieved by ODA with parameters N and η varying is shown in Fig. 3(a). We observe that 1) the QoS increases with η under the given N as more optional cycles can be executed, and 2) the QoS is not always increasing with N under the given η, because when N is changed, a new set of task parameters is generated. We further compare the QoS gain between ODA, HDA, B&B and GA in Fig. 3(b). The QoS gain "ODA vs HDA", "ODA vs B&B", and "ODA vs GA" shows the statistical property of data sets Qo(N,η)-Q h (N,η)

Qo(N,η) , Qo(N,η)-Q b (N,η) Qo(N,η)
, and

Qo(N,η)-Qg(N,η) Qo(N,η)
, where

Q o (N, η), Q h (N, η), Q b (N, η
) and Q g (N, η) are the QoS achieved by ODA, HDA, B&B and GA under the given N and η parameters, respectively. From Fig. 3(b), we observe that 1) the solutions given by B&B and ODA are same, and 2) ODA achieves higher QoS than HDA (29.8% in average) and GA (7.6% in average).

The computation time of ODA and HDA with parameters N and η varying is presented in Fig. 3(c). We observe that 1) η does not change the problem size, and, thus, its influence on computation time is limited, and 2) the influence of parameter N on the computation time of HDA is small, compared with ODA, and 3) the computation time of HDA is also influenced by using FP to solve the MP, since we need to repeatedly run FP until an arbitrary feasible solution is found. Fig. 3(d) shows the computation time gain between ODA, HDA, B&B and GA. Similarly, the computation time gain "ODA vs HDA", "B&B vs ODA", and "GA vs ODA" shows the statistical property of data sets To(N,η)-T h (N,η) , where T o (N, η), T h (N, η), T b (N, η) and T g (N, η) are the computation time of ODA, HDA, B&B and GA under the given N and η parameters, respectively. With N increasing, the computation time of ODA and B&B grows. However, ODA takes a shorter computation time than B&B (16.3% in average). For an optimization problem, its computational complexity increases significantly with the number of variables and constraints. Hence, solving iteratively smaller problems with less variables and constraints (MP and SP) is more efficient than solving a single large problem. This result is in line with the comparison in [START_REF] Randazzo | A comparison of optimal methods for local access uncapacitated network design[END_REF]. Although GA can solve non-linear programming problem, e.g., MINLP, the optimality of the solution is hard to guarantee. Compared with ODA, the GA structure is more complex, and, thus, GA has a longer computation time (32.6% in average). The problem transformation from MINLP-based PP to MILP-based PP1 is necessary, as it can simplify the problem structure, and, thus, the optimal solution is much easier to find.

We further explore the convergence of our approaches by providing the number of iterations required to find the optimal (feasible) solution for ODA (HDA) in Fig. 4(a). For ODA more constraints are added to the problem with N increasing, and, thus, a higher number of iterations is required to find the optimal solution. However, for HDA, the number of required iterations is not always increasing with N . This is because 1) the initial solution x(0) is randomly selected as long as it satisfies the constraint Ax(0) b 1 , and 2) the solution found by HDA, i.e., x (k), is an arbitrary feasible solution.

Finally, we explore the behavior of our method by using ODA to solve the PP1 (QoS-OPT) and the energy-aware task mapping problem (NRG-OPT). The NRG-OPT aims to minimize the objective function k∈M [P c k t d k + (H -t d k )P s k ] (total energy consumption) under the constraints C 1 -C 6 and C 8 -C 10 . The results show that when the NRG-OPT is used for IC-task deployment, the QoS is always equal to 0. This behavior is expected as the smaller the optional subtask, the lower the energy consumed to execute it. Fig. 4(b) shows that QoS-OPT consumes more energy than NRG-OPT, since QoS-OPT maximizes QoS and, thus, executes more optional subtasks than NRG-OPT. However, the energy consumed by QoS-OPT is always less than the supplied energy E s to satisfy the energy supply constraint. Hence, using QoS-OPT to perform IC-task deployment, we can balance QoS-enhancing with energy-usage. This study formulated the approximation-aware task deployment problem on platforms as an MILP, with a goal of maximizing QoS without violating the timing and the energy constraints. The problem is optimally solved by the proposed ODA algorithm. A novel heuristic algorithm HDA is also presented to further reduce the computation time. The numerical results show that ODA is guaranteed to converge to the optimal solution. It can reduce computation time by up to 46.8% and improve QoS by up to 21.1% compared with the existing approaches. HDA can find a feasible solution within a negligible computation time with an average cost of 29.8% in QoS compared with ODA.

Fig. 1 .

 1 Fig. 1. Task graph of illustration example.

5 . 1 :

 51 s.t. Ax b 1 , Cx + Dy b 2 ,where x (y) is the vector of binary (continuous) variables. f is the vector of the objective function coefficient. A, C and D are the constraint coefficient matrices. b 1 (b 2 ) is a u (v)-dimensional vector. Instead of simultaneously solving the binary variables x and the continuous variables y, ODA: 1) decomposes PP2 into two smaller subproblems with less variables: an ILP-based Master Problem (MP) with binary variables x and an LP-based Slave Problem (SP) with continuous variables y, and 2) solves them by using the solution of one to the other. By doing so, the computational complexity is significantly reduced.1) MP and SP formulation: Based on the structure of PP2, at the k th iteration, the MP and the SP are formulated asMP : Φ l (k) = min x, Φ Φ (4) s.t. Ax b 1 , C 11 , C 12 , SP : Φ u (k) = min y 0 f T y (5) s.t. Cx(k) + Dy b 2 ,where x(k) is the optimal solution of the MP at the k th iteration.C 11 is the Feasibility Constraints (FCs), Φ ≥ λ(i) T (Cx -b 2 ) (∀i = 1, . . . , m) and C 12 is the Infeasibility Constraints (ICs), 0 ≥ λ(j) T (Cx -b 2 ) (∀j = 1, . . . , n), where m + n = k.Remark Since the objective function of the PP2 only contains the continuous variables y, and the MP only considers the binary variables x, an auxiliary variable (continuous) Φ is introduced into the MP as the objective function. Φ and Φ have the same physical meaning, as shown in the Proof 5.1. Although the MP includes the continuous variable Φ, this problem can be solved by only considering the binary variables x[START_REF] Mo | Energy-quality-time optimized task mapping on DVFS-enabled multicores[END_REF].

(a) Step 1 :

 1 Solve the Dual Slave Problem (DSP) DSP : max λ 0 λ T (Cx(k) -b 2 ) (6) s.t. f + D λ 0. where λ [λ i ] are the Lagrange multipliers. (b) Step 2: Based on the solution of the DSP, a new constraint is generated. Case 1: If DSP is infeasible, PP2 has no feasible solution. Case 2: If DSP has a bounded solution λ(k), a new FC: Φ ≥ λ(k) T (Cx -b 2 ) is generated. Case 3: If DSP has an unbounded solution, i.e., λ(k) T (Cx(k)b 2 ) = +∞, a new IC: 0 ≥ λ(k) T (Cx -b 2 ) is generated. With a new FC added into C 11 (or a new IC added into C 12 ) at the (k + 1) th iteration, the MP is solved again to obtain a new solution x(k + 1) for the next iteration. 3) Theoretical analysis: Theorem 5.1: The non-optimal and infeasible values of the binary variables x are excluded by the constraints C 11 and C 12 .

1 -

 1 s.t. Cx(k) + Dy b 2 + ξ, λ 0, D T λ 0,

Fig. 2 .

 2 Fig. 2. The structure of ODA and HDA algorithms.

Definition 6 . 1 :

 61 ( -optimal solution) Denote (x, y) as the feasible solution of PP2, and Φ = f T y. If |Φ -Φ * | ≤ , (x, y) is an -optimal solution. Lemma 6.1: Denote Φ l (k) and Φ u (k) as the objective function values of the MP and the SP with (x (k), y (k)), respectively, where (x (k), y (k)) is the solution of PP2 found by HDA. The solution (x (k), y (k)) is an -optimal solution, where = Φ u (k) -Φ l (k).

  To(N,η) , T b (N,η)-To(N,η) T b (N,η) , and Tg(N,η)-To(N,η) Tg(N,η)

Fig. 3 .

 3 Fig. 3. QoS (a)-(b) and computation time (c)-(d) achieved by the different methods with N and η varying.

Fig. 4 .

 4 Fig. 4. Further exploration of the ODA and HDA behavior VIII. CONCLUSIONS

TABLE II FREQUENCY

 II /VOLTAGE LEVELS AND POWER MODEL PARAMETERS OF A BIG CORE AND A LITTLE CORE

					big								LITTLE
	v k (V )	0.93	0.96	1.0	1.04	1.08	1.1	1.15	1.2	1.23	v k (V )	0.9	0.94	1.01	1.09	1.2
	f k (GHz)	0.8	0.9	1.0	1.1	1.2	1.3	1.4	1.5	1.6	f k (GHz)	0.25	0.3	0.4	0.5	0.6
	Power		C

  , O l ∈ [4 × 10 7 , 6 × 10 8 ]

	TABLE III
	EXPERIMENTAL SET-UP
	Task characteristics
	M l

ACKNOWLEDGMENT

This research is funded by INRIA post-doctoral research fellowship program, and is partially funded by ANR ARTE-FACT (AppRoximaTivE Flexible Circuits and Computing for IoT) project (Grant No. ANR-15-CE25-0015).