Jacek Kabziński
email: jacek.kabzinski@p.lodz.pl

Is Extreme Learning Machine Effective for Multisource Friction Modeling?

Keywords: extreme learning machine, friction modeling, adaptive control 1

The aim of this contribution is to discuss suitability of extreme learning machine (ELM) approach for modeling multisource friction for motion control purposes. The specific features of multisource friction in mechatronic systems are defined, the main aspects of friction modeling by a standard ELM are investigated and some modifications are proposed to make it more suitable for specific demands of the discussed task. This allows to formulate some general remarks concerning properties of ELM for function approximation.

Introduction

Extreme Learning Machine (ELM) is relatively new neural modeling technique. It is intensively promoted by numerous researchers and attracts increasing attention for its fast and simple training algorithm. Several review papers [START_REF] Huang | Trends in extreme learning machines: A review[END_REF][START_REF] Liu | Is Extreme Learning Machine Feasible? A Theoretical Assessment (Part I)[END_REF][START_REF] Lin | Is Extreme Learning Machine Feasible? A Theoretical Assessment (Part II), Neural Networks and Learning Systems[END_REF] describe the basic form of ELM and discus characteristic features of the method. The range of possible applications of ELM covers as well regression and classification problems and learning schemes include batch processing and sequential on-line training. In this contribution we concentrate on function approximation problems with batch data processing approach. The authors of ELM stress that the network structure is simple, hidden neurons parameters are selected on random "without seeing the data" and remain constant, output weights are calculated once for a batch of data, with the use of Moore-Penrose generalized inverse. Hence, without iterative tuning of parameters by back-propagation, ELM achieves very short learning times and, with sufficient number of neurons, sufficient accuracy in terms or modeling RMSE. On the other hand, several publications report that ELM may: create ill-condition numerical problems [START_REF] Zhixiang | A modified extreme learning machine with sigmoidal activation functions[END_REF], introduce overfitting, require too big number of neurons [START_REF] Liu | Is Extreme Learning Machine Feasible? A Theoretical Assessment (Part I)[END_REF], that the influence of number of neurons and parameters of random weights selection on the resulting performance is unclear and requires analysis [START_REF] Parviainen | Connection between Extreme Learning Machine and Neural Network Kernel[END_REF]. Having in mind the "no free lunch" rule, i.e. remembering that it is impossible to find "the best" solution for all possible problems, we investigate the effectiveness of ELM modeling for a special problem -friction (i.e. friction force or torque) modeling for motion control in mechatronic systems. This problem is important for applications, and it constitutes an interesting benchmark for modeling techniques. Although several memoryless models of friction as a function of relative speed of moving surfaces are available [START_REF] Wojewoda | Hysteretic effects of dry friction: modelling and experimental studies[END_REF], parametrical identification of these models is difficult because of high nonlinearity. The problem becomes more complex if several sources of friction are present in the plant, acting together against the motion. In this case the use of AI models is acceptable, and the control algorithms benefit from linear parameterization of such models [START_REF] Kabziński | One-Dimensional Linear Local Prototypes for Effective Selection of Neuro-Fuzzy Sugeno Model Initial Structure[END_REF]. In this contribution the main features of friction modeling by a standard ELM are investigated and some modifications are proposed to make it more suitable for specific demands of the discussed task.

Multisource Friction Model

The friction modeling and parametric identification is extremely complex if the drive system has to match friction forces from several sources. A simple example is a mover acting between two different surfaces. As a matter of fact the presence of several sources of friction is a typical situation. Usually friction is caused by a load machine, bearings, a motor, etc. Of course it is possible to use several different models of friction -one for each recognized source, but this will increase the number of parameters drastically and will make the whole model too complicated. As all friction forces accumulate, working against the motion, it will be impossible to distinguish among the models and to identify all parameters. The multisource friction model may be defined as

= + - , + (1)
where S is the number of sources, v is the velocity between two surfaces in contact. For each source of friction is the Coulomb friction, the Stribeck effect is described by and the Stribeck velocity is , while denotes a coefficient of viscous friction [START_REF] Wojewoda | Hysteretic effects of dry friction: modelling and experimental studies[END_REF]. The function , influences the shape of the characteristics. Several possibilities were reported in the literature [START_REF] Wojewoda | Hysteretic effects of dry friction: modelling and experimental studies[END_REF], for example

, = , , = | | , , = ! | | " . (2)
The exact form of , is not justified theoretically and must be adjusted to the considered application. Generally, it represents the transient between the boundary lubrication and the hydrodynamic lubrication zones, it starts with 1 for = 0 and disappears for big . The data collected to identify the friction model are obtained during some experiments, from available measurements, or from a proper observer output [START_REF] Kabzinski | Fuzzy Friction Modeling for Adaptive Control of Mechatronic Systems[END_REF]. Let us denote the measured or estimated values of friction by %

. The common features of the obtained data are:

• Changing parameters , , results in various shapes of the characteristics , but several values of parameters may lead to the same curve , so the identification of all parameters in 1 by any curve fitting method is impossible.

• The data are subject to measurement or observation inaccuracies, that may be represented as a uniformly distributed noise.

• It is possible to identify the resulting viscous friction coefficient = ∑ from the steady-state data, collected for sufficiently big velocities, as ∆ ≈ ∆ , although it is impossible to distinguish among .

• The number of the data is limited. Measurement/estimation generates an uniform random noise in [-0.5, 0.5] N. The data are normalized to the range [0 1] by scaling factors L MNO , L PQ for the output and input respectively. As = ∑ is identified using linear regression for the data obtained from big speed motion and , ≈ 0 for > 3 FGH , the final model, obtained from the data S T P , U P = VL PQ P , L MNO W % P -P XY , Z = 1, … , \] will be

F = ^LMNO _`a L PQ Z + bc | | < 3 FGH L MNO _`a 3L PQ FGH Z + bc | | > 3 FGH . ( 3
)
where _`a • stands for the Artificial Intelligence Model, that will be provided by EML or other methods. It is assumed that the Stribeck curve is symmetrical and modeling for positive velocity allows to obtain a complete model, although generalization to nonsymmetrical characteristics is straightforward. \ =200 points are selected on random from 0, 3 FGH and are used for training. The point W0, 0 X, where 0 is calculated as the mean of measurements for four smallest values of velocity, and W3 FGH , 3 FGH X are added to the learning dataset to obtain the model covering the velocity interval 0, 3 FGH . The friction curves are plotted in fig. 1

Fuzzy Friction Modeling

The fuzzy friction modeling technique was discussed in details in [START_REF] Kabzinski | Fuzzy Friction Modeling for Adaptive Control of Mechatronic Systems[END_REF]. It is known that the careful initial model selection [START_REF] Kabziński | One-Dimensional Linear Local Prototypes for Effective Selection of Neuro-Fuzzy Sugeno Model Initial Structure[END_REF] enables to obtain accurate and relatively simple Takagi-Sugeno-Kang models, and the use of nonlinear consequences helps to improve the model performance [START_REF] Kabziński | TSK Fuzzy Modeling with Nonlinear Consequences[END_REF]. Here fuzzy modeling is used only for comparison with EML models. The standard initialization technique implemented as genfis1 in Matlab fuzzy toolbox is used together with AFNFIS approach to model training [START_REF] Jang | ANFIS: Adaptive-network-based fuzzy inference system[END_REF]. The results obtained for different number of bell-shaped membership functions (MFs) are presented in fig. 2,3. Although the training error is decreasing with number of membership functions, the test error is increasing, starting from very moderate number of MFs. The overfitting is visible in fig. 2. Hence the lack of generalization is the most important drawback of fuzzy modeling. The choice of 4 membership functions seems reasonable, it makes the model to be linear subject to 8 parameters and nonlinear subject to 12 more, so 20 parameters are optimized. The performance of 4 membership functions model is presented in fig. 5

Standard EML Model

ELM was originally proposed by Huang et al. [START_REF] Huang | Trends in extreme learning machines: A review[END_REF] as a learning algorithm for the singlehidden-layer neural network. The main concepts behind the ELM are that: (i) the weights and biases of the hidden nodes are generated randomly and are not adjusted and (ii) the output weights are determined analytically. The output function of ELM is of the form (for a one output case, T is the input variable, i.e. the normalized velocity in friction modeling problem):

opq T = r ℎ T q ( 7
)
where M is a number of hidden neurons, ℎ T = ℎ T , … , ℎ F T k is the output vector of the hidden layer with respect to the input v, and r = r , … , r q k is the vector of the output weights (between the hidden layer and the output node). Different activation functions (AFs) may be used in the hidden neurons. The most popular are: sigmoid, Gaussian, multiquadric, step, triangular and sine AFs. The sigmoid AF is:

ℎ T = 1 1 + Tt -u T + v (8)
where u , I = 1, . . a is the input weight of the k-th hidden neuron, v is its bias. The parameters of hidden nodes u and v , are randomly generated according to any continuous probability distribution instead of being explicitly trained. That process is independent of the training data. If S T P , U P = VL PQ P , L MNO W % P -P XY ,

Z = 1, … , \] (9)
compose a set of learning data, the N-dimensional vector of model outputs F is compared with the target vector F, where

F = opq T , … , opq T j k , = U , … , U j k (10)
As

F = w ℎ ⋯ ℎ q ⋮ ⋱ ⋮ ℎ j ⋯ ℎ q j { w r ⋮ r q { = |r (11)
the output weights r are found by minimizing the approximation error:

e = ‖|r -‖ i (12)
The optimal solution is:

r M~O = | ! (13)
| ! is the Moore-Penrose generalized inverse of matrix | . The standard configuration of ELM assumes sigmoid AFs with random parameters uniformly distributed in [-1,1]. This was used for the benchmark modeling. The model performance is demonstrated in fig. 9. Mean values from 20 runs are presented. As it is plotted in fig. 10, the output weights and the condition coefficient of matrix || k became very big numbers. Therefore using more than few hidden neurons does not make any sense. The results obtained for 6 hidden neurons are presented in fig. 11. The achieved modeling errors are e OfGPQ , e Og O , e kOg O = 0.0232, 0.0233, 0.0072 and are comparable with the scores obtained by the fuzzy model. The generalization properties of standard ELM are much better -no overfitting with increasing number of neurons is observed. The noise corrupting the data set does not allow to decrease the modelling errors below the limits obtained with several hidden neurons, and from numerical point of view the model becomes ill-conditioned for bigger number of neurons.

Improving ELM

The ELM drawbacks described in the previous section motivate the for modifications that will improve the modeling performance, i.e. the model accuracy, ability of generalization and the numerical stability.

Activation functions

The change of the activation function type from sigmoid to any other infinitely differentiable function, with the same weights and biases selection procedure, will not improve the model. The errors behave very similarly for different smooth AFs (fig. 13) and also •b € | k | and •r M~O • remain at comparable levels. Of course the application of AFs with a compact support, like membership functions used in fuzzy modeling will change the situation, making ELM to behave like a fuzzy model.

Regularization

The well-known remedy for numerical problems in ELM, caused by ill-conditioned set of linear equation, is regularization of the least-square problem [START_REF] Huang | Trends in extreme learning machines: A review[END_REF]. Instead of minimizing (12) the weighted problem is considered, with the performance index:

e ' = ‖r‖ i + ƒ‖|r -‖ i ( 14
)
where C is a design parameter. The optimal solution is now [START_REF] Huang | Trends in extreme learning machines: A review[END_REF]:

r 'M~O = 1 ƒ `+ | k | | k (15)
and decreasing C allows to improve the condition number of ' `+ | k |. Inevitably this will degrade the modelling accuracy, so a compromise is required. The proper choice of C depend on the problem structure and it is difficult to formulate general rules. In 15). Moreover the model instability (growing errors) is observed for N>50. Smaller values of C will stabilize the model but the modeling errors will be bigger.

Fig. 15

Average training and test errors versus no of hidden neurons for C=10 5 and sigmoid AFs.

Fig.16

Average training and test errors versus no of hidden neurons for enlarged variability of sigmoid AFs.

Activation function variability

The ELM is not able to learn the features from the data, as a fully trained neural network does. It projects the input data into the feature space the randomly weights happen to specify, and learns a linear mapping in that space. Therefore this nonlinear mapping of the data into a feature space should be able to extract the features sufficient for predicting the target variable of a regression task. The importance of sufficient AFs variability was noticed previously [START_REF] Parviainen | Connection between Extreme Learning Machine and Neural Network Kernel[END_REF]. The point is, to find the way to correct the random mechanism of weights and biases creation to increase the variability without losing too much from the simplicity of random selection. The first step to enlarge the variability of sigmoid AFs is to increase the range of uniform random selection of input weights. It must be suitably fixed to guarantee that the sigmoid operation neither remains linear nor too strongly saturates in the input interval.

In the presented application the range u ∈ -10,10 seams suitable. Next the biases may be selected to ensure that the range of sigmoid function is sufficiently large. For a positive input weight u > 0 sigmoid function is strictly increasing in [0,1]. Assuming that the minimal value of the sigmoid AF shell be smaller than the given c , and the maximal value should be bigger than c i , 0 < c < c i < 1 yields

ℎ 0 = 1 1 + Tt -v < c , ℎ 1 = 1 1 + Tt -u -v > c i (
… † = -u + ˆ 1 c i -1 , … ‡ = -ˆ 1 c -1 (17)
If the selected weight is negative u < 0, the bias v shall be selected on random, uni- Even better results are obtained applying the variability enlargement and the regularization together. The results for 10 neurons and C=10 5 are plotted in fig. 19. The achieved modeling errors are e OfGPQ , e Og O , e kOg O = 0.0231, 0.0234, 0.0074 . As it is visible from fig. 20 the generalization ability is preserved for increasing number of neurons -the training and test errors remain at the same level. No overfitting was observed for number of neurons<50. It is remarkable that the norm of output weights remains small: ‖r‖ < 6.

formly from v ∈ … ‡, … † ,

Direct minimization of the condition coefficient

Good conditioning of the linear equation ( 11) is vital to derive accurate ELM with moderate output weights. Therefore some hidden neurons parameters may be not selected randomly, but may be used to minimize •b € | k | directly. So, the hidden neurons will be divided among two subsets, say RAND and OPT. First, the AF parameters of neurons in RAND will be selected on random, next the parameters of neurons in OPT will be calculated to minimize Of course in this case we cannot claim that the weights and biases selection process is so fast and simple as in the basic ELM algorithm. The choice of the optimization algorithm is also an open problem. It is difficult to assume anything about the shape of the minimized function, as it depends on the training data and randomly selected parameters of neurons in RAND. The exemplary plot of the condition number as a function of the weight and the bias, when first ten sigmoid neurons belong to RAND and the neuron number 11 is the only one in OPT, is presented in fig. 21. As it is visible, the minimized function is not smooth and multimodal. Therefore gradient type algorithms are not applicable, and good heuristic algorithm, that is able to avoid the local minima is recommended. The optimization result will be influenced by the initial condition, as well as by the range of the decision variables, i.e. the weights and biases of the neurons in OPT. All these makes the whole process quite complicated and time consuming. The result highly depends on the effectiveness of the applied optimization algorithm. Selecting a half of neurons as RAND with random selection with the enlarging the AFs variability approach, and the other half in OPT with minimization performed by Particle Swarm Optimization algorithm, provides results presented in fig. 22. The achieved modeling errors are e OfGPQ , e Og O , e kOg O = 0.0231, 0.0235, 0.0076 . The mean value of the output weights norm is smaller than 17 for any number of neurons. So the results are comparable with these provided by application of AF variability enlargement and regularization together.

•b € | k | .

ELM Models in Adaptive Control

The final test for the ELM model of friction is its behaviour in a control system. Here we investigate the applicability of the proposed models in a simple adaptive backstepping, tracking position control. The motion dynamics of the considered system is given by

v dt dx = friction e F F dt dv m - = . ( 18
)
where m is a forcer mass, Fe is a thrust force, and Ffriction represents all types of friction forces. The mover speed is v and the position is x. Although equation ( 18) is written according to linear motion convention, it may be also used for rotational movement description, if we read m as a moment of inertia and consider torques instead of forces. It is assumed that the thrust force is proportional to the control variable (the motor current)

) ( ) ( t i t F e ⋅ = ϕ , (2)
and coefficient ϕ is not known exactly. The desired smooth position trajectory will be denoted by d

where 1 1 0 > k is a design parameter, we will be able to describe the tracking error dynamics as

+ ⋅ - = + ⋅ - - = & & & , ( 4
) v e k x v v e d d - ⋅ + = - = 1 1 2 & (5) ϕ ϕ m m F i v m v m v m e m o friction d o o d o o = - - ⋅ = ⋅ - ⋅ = ⋅ , 1 ( ) 2 1 1 1 e e k k x v d d + ⋅ - ⋅ + = & & & . ( 7
)
The control variable i has to compensate the function

friction d o F v m D ϕ 1 - ⋅ = & , (8)
and to assure fast tracking. We will use a model D ˆ for D, incorporating the fuzzy friction model (17,18). The general structure of D ˆ will be given by

ξ A T D = , (9)
where T A ˆ is a vector of adaptive parameters and ξ is known. The structure with two adaptive parameters:

[ ] k m o T , = A         = ) ( 1 v F v Fm N d ϕ & ξ , (10)
⋅ - + = ⋅ ξ ε A & . ( 29
)
The Lyapunov function ( )

A Γ A 2 1 1 2 2 2 1 - + ⋅ + = T o e m e V . ( 11
)
with positive definite symmetric matrix Γ is selected to investigate the system stability. Taking any of the adaptation laws ,e e are uniformly ultimately bounded [START_REF] Khalil | Nonlinear Systems[END_REF]. For example with adaptation performed according to (31a) we get

ξ Γ = 2 ˆe A & , or ΓA A δ ξ - Γ = 2 ˆe & , or ΓA A 2 2 2 1 2 ˆe e e + - Γ = δ ξ & . ( 12a
2 2 2 2 2 1 1 2 2 2 2 2 1 1 2 1 2 1 ε ε +       - - - ≤ ⋅ + ⋅ - + ⋅ - = e k e k e e k e k V & (13)
and is negative outside

2 max 2 2 2 1 1 ε k e e > + ,       - = 2 1 , min 2 1 k k k (14)
Finally the control algorithm with the adaptation law (31c) was selected. Cooperation with several friction models was tested (fuzzy model with 4 membership functions and 3 different ELM models). Although all models may be used for adaptive control, it must be stressed that the more accurate fuzzy model results in smaller control signals. Taking this aspect into account, the ELM model with 10-neurons, enlarged variability of AFs and regularization (presented in fig. 19) was the best one. The application of models with big values of output weights is sensitive to numerical aspects of control implementation -the discretization time must be shortened. The proposed modeling technique was also tested on the real life data collected from a linear permanent magnet motor connected with a pneumatic linear motor. So at least three sources of friction are present: 1) the friction in liner bearings of the electric motor, 2) the friction in linear bearings of the pneumatic motor trolley, 3) the friction between the pneumatic motor cylinder and piston. The results confirm conclusions from simulations presented here.

Conclusions

The conclusions from the presented model comparison are collected in table 2.

The presented grade concerns the friction modeling problem, but some conclusions seem to be general for curve approximation by ELM.

• The increase of number of hidden neurons is often limited by numerical instability or brings insufficient benefits because of noisy data. Engineering applications with moderate number of neurons will provide appropriate accuracy of models and facilitate practical implementations. • The regularization itself is often not sufficient to improve the conditioning and strongly spoils the accuracy.

• The sufficient variability of the applied activation functions is very important.

It does pay to modify the random choice of parameters to improve the variability and to make the feature mapping more effective. The presented research demonstrates that extreme learning machine provides an interesting alternative to fuzzy modeling, although the application is not so universal and straightforward as it is recommended by many publications.

Fig. 2 .

 2 Fig.1 3-source friction. Fig. 2. The data collected for training and testing.

Fig. 3 Fig. 4

 34 Fig. 3 Training and test errors versus no of bell-shaped membership functions

, 6 , 7 .

 67 The obtained values of train, test and theoretical test errors are e kOg O , e Og O , e kOg O = 0.0224, 0.0241, 0.0075 . Similar results were obtained with other types of membership functions. For example the errors achieved with Gauss-type MFs are plotted in fig. 8.

Fig. 5

 5 Fig. 5 Fuzzy model with 4 MFs and the trading data.

Fig. 6

 6 Fig.6 Training and test errors versus no ofGaussian-type membership functions.

Fig. 7

 7 Fig. 7 Bell-shaped MFs of the trained fuzzy model

Fig. 9

 9 Fig.[START_REF] Kabziński | TSK Fuzzy Modeling with Nonlinear Consequences[END_REF] Average training and test errors versus no of hidden neurons.

Fig. 10

 10 Fig.10 Average norm of output weights and average •b € | k |

Fig. 11

 11 Fig. 11 Modeling with 6-neurons ELM Fig.12 Activation functions of 6-neurons ELM

Fig. 13

 13 Fig. 13Average training and test errors versus no of hidden neurons for sigmoid, sinus and radial basis AFs.

Fig. 14

 14 Fig.14 •b € V ' `+ | k |Y as a function of number of neurons and regularization coefficient C

  the considered case (200 samples, sigmoid AFs) an exemplary plot of •b € V ' `+ | k |Y is presented in fig. 14. Unfortunately, if C provides acceptable condition coefficient, the accuracy of the model is too poor (fig.

Fig. 18

 18 Fig.18 Activation functions of 10-neurons ELM and enlarged variability of AFs.

Fig. 19 5 .Fig. 20

 19520 Fig. 19 Modeling with 10-neurons ELM, enlarged variability of AFs and regularization with C=10 5 . Fig. 20 Average training and test errors versus no of hidden neurons for enlarged variability of sigmoid AFs and regularization with C=10 5 .

Fig. 21

 21 Fig. 21 An exemplary plot of log •b € | k | as a function of the 11 th neuron parameters.

Fig. 22

 22 Fig. 22 Average training and test errors versus no of hidden neurons for optimization of conditioning by the half of neurons.

x , actual position by x and the tracking error by x x e d - = 1 .

 1 The velocity will be 'virtual control' for position tracking. If we choose the desired velocity d

  was selected, where F is the friction model (3) (with the derived ELM model inside), o m ˆ is supposed to adapt the changing inertia and k ˆ corrects

1 .

 1 It is assumed that "the best" adaptive parameters T * A exist, such that the model εmax < ∞ . The difference between actual and "the best" parameters is denoted as

  ,b,c) with small positive δ , allows to prove that the Lyapunov function derivative (32) is negative outside a certain, bounded set, and so 2 1

•

  It is possible to find reasonable bound for the biggest Stribeck velocity maxC D < FGH . The benchmark problem considered here is obtained from the model (1) with parameters given in table 1.

Table 1 .

 1 Parameters of 3-source friction model (SI units).

	I = 1 I = 2 I = 3	4 6 5	8 16 20	0.01 0.005 0.02	2.5 50 5

  and the normalized data T P , U P are presented in fig.2. The testing dataset consists of another 200 points selected on random from 0, 3 FGH RMSE is used to measure the training error e OfGPQ , or the test error e Og O

	Also so called "theoretical testing error", compering the model with the exact function
	used to generate he data, is considered in simulations. It is calculated using \t=500 values distributed equidistantly in 0, 3 FGH
	e kOg O = h	1 \U	jO P	W_`a T P -L MNO	P -P X	i
		e OfGPQ = h	1 \	P j	_`a T P -U P	i	(4)

8

  Training and test errors versus no of Gaussian-type membership functions

		1													0.0234
												Training Data
												ANFIS Output
		0.9
															0.0232
		0.8
															0.023
	Ffriction [normalized]	0.5 0.6 0.7												training error	0.0228
															0.0226
		0.4
															0.0224
		0.3
		0.2	0	0.1	0.2	0.3	0.4	0.5 v [normalized]	0.6	0.7	0.8	0.9	1	0.0222	0	10	20	30	40	50 iterations	60	70	80	90	100
															0.03
		in1mf1				in1mf2			in1mf3			in1mf4											train err.
		1																						test err.
																									teor.test err.
															0.025
		0.8
	Degree of membership	0.4 0.6											FIS errors	0.015 0.02
		0.2
															0.01
		0
															0.005	2	3	4	5	6	7	8	9	10	11	12
			0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1								No of MFs
									input1

  where … ‡, … † are given in(17). The application of the proposed method to the considered friction modeling problem results in the performance presented in fig.16. The achieved modeling errors are e OfGPQ , e Og O , e kOg O = 0.0210, 0.0234, 0.0070 . The results for 10 neurons are presented in fig.17,18. Comparing figures 18 and 12 we notice that the variability of AFs was noticeably increased. The output weights are smaller and ill-conditioning is eliminated.

	Fig. 17 Modeling with 10-neurons ELM and
	enlarged variability of AFs.

Table 2 .

 2 Suitability of various models for multisource friction modeling

			Accuracy	Generaliza-	Training time Numerical
				tion		conditioning
						and stability
	Fuzzy model	sufficient	poor	moderate	good
	Basic ELM	sufficient	moderate	short	bad
	ELM	with	poor	poor	short	good -if we
	regularization				compromise
						with the accu-
						racy
	ELM with en-	good	good	short	good
	larged varia-
	bility of AFs
	ELM with en-	good	the best	short	the best
	larged varia-
	bility of AFs +
	regularization
	ELM enlarged	good	good	long	good
	variability of
	AFs and con-
	ditioning opti-
	mization

& & & & ,(6)


Mehmet Dilan Mert
email: dilanmert@gmail.com

Ulvi Şimşek

Suat Özdemir
email: suatozdemir@gazi.edu.tr

Privacy Preserving Metering Protocol in Smart Grids

Keywords: Smart Grid, Privacy Preserving, Time Perturbation, Secret Sharing Scheme

Due to dramatically increasing energy need traditional grid becomes inadequate to manage the energy generation and distribution system. Smart Grid is an efficient, user-friendly solution to today's energy management and generation problem. However, Smart Grid introduces many privacy issues that do not exist in traditional grids. For example, with Smart Grid it is possible to obtain consumers' daily life patterns from their energy usage data. In this paper, we tackle with this problem and propose a novel privacy-preserving communication protocol. The proposed protocol is based on time perturbation and Shamir's Secret Sharing (SSS) scheme. In the proposed protocol, consumption reports are generated by smart meters and sent to a data collection center using Laplace Distribution based time perturbation. In addition, SSS is employed to prevent attacks during multi-hop data transmission. The time perturbation prevents malicious users to see the actual data while it is stored in the data center. Security analysis shows that the proposed protocol preserves privacy of consumer data during both data transmission and data storage.

Introduction

Traditional power grid is considered inefficient and unreliable since it just distributes the energy to consumers but does not provide any control over power consumption.

In order to satisfy the rapidly increasing energy demand, the energy suppliers not only increase energy generation capacity but also allow consumers to control the way that they consume energy. Smart Grid (SG) integrates traditional power grid and information technologies to provide efficient monitoring and control of energy generation and consumption. SG usually consists of innovative sensing and control systems that are able to perform real time monitoring of power generation, transmission and usage. SG enables consumers and energy providers to analyze consumption data and forecast the power utilization [START_REF] Rottondi | Distributed Privacy-Preserving Aggregation of Metering Data in Smart Grids[END_REF][START_REF] Erkin | Privacy-Preserving Data Aggregation in Smart Metering Systems[END_REF].

A typical SG has three main segments: power generation unit, transmissiondistribution network, and smart meters (as shown in Fig. 1). In this paper, we focus on smart meters that present a number of challenges in sensing, analyzing, and communication. Smart meters offer several new opportunities that do not exist in traditional meters. For example, a smart meter can be remotely read and provide shorter finegrained energy reading intervals composed to traditional ones. This allows users to achieve more efficient energy usage and hence balance overall network status.

Fig. 1. A Typical Smart Grid

Despite its advantages, SG introduces many risks and challenges regarding security and privacy of both consumers and energy providers. Consumers' privacy takes place on the top of these concerns. Smart meters directly cause many of these challenges. For example, disclosing energy consumption data measured by smart meters creates a privacy threat for consumers. Using this data the daily life patterns like presence/absence, number of people or usage pattern of home appliances can be obtained [START_REF] Rottondi | Privacy-Preserving smart metering with multiple data consumers[END_REF].

In multi energy provider SG's consumers' data are usually stored in a common data storage environment (such as cloud) where each provider should only access its own consumers' data. However, cloud based common storage systems pose another privacy problem. In such an environment, an honest-but-curious (HBC) attacker can perform a time series analysis on consumer data and obtain his/her daily life pattern.

To address above issues, in this paper, we propose a protocol that protects the consumers' privacy during both data transmission and storage. The proposed protocol employs SSS scheme that enables smart meters partition the measurement reports and send each part of the report to be transmitted over a different path. Hence, during data transmission eavesdroppers cannot obtain the consumption or patterns unless they predict and listen all the paths. In addition, smart meters do not instantly send their measurements but send them in batches. Before transmitting these measurement batches smart meters use a time based perturbation scheme to measured values to prevent time series analysis of data at the data collection server. To the best of our knowledge, this is the first work that employs SSS scheme and time based perturbation in SG to ensure consumer privacy during both data transmission and storage. In the literature, there are studies that use these techniques to ensure privacy during data transmission, however our work differs from them by considering not only data transmission but also data storage and ensuring the privacy of data both during data transmission and data storage. Our work shows that combination of SSS scheme and temporal perturbation algorithms for SG is enough to ensure consumer privacy during both data transmission and storage for HBC and eavesdroppers. In addition, the proposed protocol is intended to provide privacy at the data collection server where the measurement data of multiple consumers is stored.

The rest of this paper is structured as follows. In Section II, security requirements for SG are explained in detail. In Section III, we present system and threat model along with assumptions. The proposed protocol is explained in Section IV. The performance evaluation is presented in Section V using theoretical analysis and simulation results. Finally, the paper is concluded in Section VI.

Security Requirements of Smart Grids

Compared to traditional grid structure SG requires higher degrees of network connectivity to support the new sophisticated features that bring new vulnerabilities [START_REF] Özdemir | Security Technology for Smart Grid Networks[END_REF]. SG incorporates smart metering equipment and traditional grid to obtain information about electricity consumption of consumers. For this reason, user privacy is an important issue. Since electricity consumption patterns not only show how much energy consumed but also give what appliances are used, when the consumer is at home, at work, or travelling [START_REF] Hadley | Smart-Grid Security Issues[END_REF]. It is shown that this kind of data is valuable and important for financial or political analysis [START_REF] Lu | EPPA: An Efficient and Privacy-Preserving Aggregation Scheme for Secure Smart Grid Communications[END_REF].

To prevent malicious actions that can threat consumer privacy, the following security requirements must be provided in SG communications [START_REF] Lu | EPPA: An Efficient and Privacy-Preserving Aggregation Scheme for Secure Smart Grid Communications[END_REF].

 Confidentiality: This service protects individual consumer's reports from attackers. It prevents unauthorized disclosure of information that is not open to public and individuals [START_REF] Lu | EPPA: An Efficient and Privacy-Preserving Aggregation Scheme for Secure Smart Grid Communications[END_REF][START_REF] Komninos | Survey in Smart Grid and Smart Home Security: Issues, Challenges and Countermeasures[END_REF].  Data Integrity: The accuracy and consistency of data must be provided during data transmission and storage [START_REF] Yan | A Survey on Cyber Security for Smart Grid Communications[END_REF][START_REF] Komninos | Survey in Smart Grid and Smart Home Security: Issues, Challenges and Countermeasures[END_REF].  Authorization: It refers to prevent people or systems to interact with SG without permission. It provides a way to distinguish between legitimate and illegitimate users [START_REF] Yan | A Survey on Cyber Security for Smart Grid Communications[END_REF].  Non-repudiation: This service ensures that users and energy service provider cannot deny a transmission of a message that has real time pricing and metering data [START_REF] Yan | A Survey on Cyber Security for Smart Grid Communications[END_REF][START_REF] Wang | A survey on the communication architectures in smart grid[END_REF].

3

System Models and Assumptions

System Model

The proposed system model consists of three components (as shown in Fig. 2). First, we assume that a smart meter is available in every house or industrial building that measures electricity usage and report the measurement to the data center over other smart meters via power line and/or wireless communication. A smart meter generates reports in every time slot and sends this information to data collection center in batches. Before sending of the report, a smart meter buffers the report for random time to perturb the measurement time of the report. Then, the smart meter divides the report to several shares using SSS scheme. Each share of the report is transmitted over a different path so that the privacy of the report is ensured. Consequently, to obtain the whole report, an attacker must listen over each of the distinct paths that are used to forward the shares. How the reports are routed to the data collection center is out of scope of this paper. Our second assumption is that the data collection center is directly connected to SG via at least n smart meter. The task of the center is to collect the shares of each report and recover the original report for the users such as billing companies, distribution / transmission systems or third-party services request detailed information for billing or load balancing or another purpose. Finally, there are billing companies, distribution / transmission systems or third-party services that request data collection center to get its own consumers' data.

Fig. 2. System Model

Threat Model

In this paper, we consider eavesdroppers and HBC type adversaries. An eavesdropper listen communication channels to obtain information about consumers. HBC adversaries are generally located on data centers. They follow the procedures correctly and do not provide any false information on purpose. However, a typical HBC attacker is curious to probe the detailed privacy of the consumers. For example, he/she can be an insider operator who secretly sells out consumer profiles. He is interested in the indi-viduals' time-series measurements, and attempts to infer the detailed information about participants' behaviors and activities.

The Proposed Protocol

In this section, we propose our privacy-preserving smart metering protocol for SG communications. The proposed protocol ensures consumers' privacy both during data transmission and storage. In order to achieve this, the proposed protocol employs SSS scheme that enables smart meters partition reports and send each part of the report to be transmitted over a different path. Hence, during data transmission eavesdroppers cannot obtain the consumption or patterns unless they predict and listen all the paths. In addition, smart meters do not instantly send their measurements but send them in batches. Before transmitting these measurement batches smart meters use a time based perturbation scheme to measured values to prevent time series analysis of data at the data storage environment. In what follows we explain the techniques used in the proposed protocol in detail.

Time Perturbation

Existing work about perturbation is divided into two parts, data-oriented and contextoriented. In data-oriented approaches, the measurement data is protected from external attacks like eavesdropping communication or internal attackers like compromised nodes. In context-oriented approaches, timing or location information is protected [START_REF] Li | Privacy preservation in wireless sensor networks: A state-of-the-art survey[END_REF]. In real-time applications, consumer's privacy reveals with obtaining time-series measurement data by attackers. The outcome is inferring consumer's behavioral pattern to access time information of measurement data. Temporal perturbation is based on hiding the time information of smart meter messages to preserve electricity usage pattern from attackers. In the literature, temporal perturbation is widely used to provide privacy. For example, in [START_REF] Kamat | Temporal Privacy in Wireless Sensor Networks[END_REF], the authors propose a buffering scheme that uses exponential distribution to create obfuscation on temporal information of sensing data in a WSN. A package-buffering scheme is proposed to obfuscate temporal information for the problem of temporal privacy.

In our protocol, the temporal perturbation method is basically works as explained.

Smart meters read measurement at each time slot and buffer them for a random delay to send the data collection center. The delay must be an advanced time that the data is buffered for sending to other nodes at the perturbed time. With this perturbation, the data cannot be changed and sent at reading time directly. The resulting time-series of data do not represent the original sequence. For this reason, with perturbing time information of messages, attackers cannot directly retrieve to data and cannot infer daily life pattern of consumers'.

Shamir's Secret Sharing Scheme

Shamir's Secret Sharing (SSS) scheme was proposed in [START_REF] Shamir | How to share a secret[END_REF] for sharing secrets in privacy. According to the scheme, a secret is divided in w parts (shares) and w parts are sent to the target. In order to recover the secret at least t shares are required where 𝑡 ≤ 𝑤 [START_REF] Rottondi | Distributed Privacy-Preserving Aggregation of Metering Data in Smart Grids[END_REF].

The SSS scheme works as explained below [START_REF] Rottondi | Privacy-Preserving smart metering with multiple data consumers[END_REF]: In each smart meter, after measuring the electricity consumption, the time of measurement is perturbed based on the proposed protocol. Then, prepared report is divided into w shares and then all generated w shares sent to the data collection center over other smart meters.

 𝑚 ∈ 𝑍

Proposed Protocol

We propose a privacy preserving communication protocol, which consists of four parts; generate report, perturb temporal information and generate-send report shares and recover-store reports. The first three-phases are performed at every round in every smart meter. The last phase is performed at the data collection center that manages storing reports and presenting them to the third-party organizations.

Generate Report

To provide a real-time measurement; at every time slot, a consumer's electricity consumption measured by the smart meter. The smart meter ID and a time stamp are appended to this measurement. The time stamp is used to prevent replay attacks and to realize the time of measurement.

𝑴

To ensure data integrity of the generated report, the seed information is generated with a hash function and added to report as follows.

𝒔𝒆𝒆𝒅 = 𝐻(𝑚𝑒𝑡𝑒𝑟 𝑖𝑑 * 𝑇 𝑗 * 𝐾 𝑖 ) (2)
𝑴 𝒊 𝒋′ = 𝑚𝑒𝑡𝑒𝑟 𝑖𝑑 + 𝑡𝑠(𝑇 𝑗 ) + 𝑚𝑒𝑎𝑠𝑢𝑟𝑒𝑚𝑒𝑛𝑡 + 𝑠𝑒𝑒𝑑

Perturb temporal information

To disorientate attackers who eavesdrops communication channels between smart meters and data collection center, the time information of report is perturbed before report being sent. After generating report and the seed, to perturb time information, below steps are performed.

number N is produced according to Laplace (0, b) distribution with zero mean and the variance of 2b 2 . If N is greater than or equal to 0, below steps are followed. But if N is less than 0, continue to produce N until greater than 0.

The time information is updated with N like follows;

𝑇 𝑗 ′ = 𝑇 𝑗 + 𝑁 (4) 𝑡𝑠(𝑇 𝑗 ′ ) = 𝑡𝑠(𝑇 𝑗 + 𝑁) (5)
To indicate the real time slot to the center, N value is sent secretly with the new message:

𝑁 ′ = 𝐸(𝑁) (6)
The message is constructed by new time information;

𝑀 𝑖 𝑗′′ = 𝑚𝑒𝑡𝑒𝑟 𝑖𝑑 + 𝑡𝑠(𝑇 𝑗 ′ ) + 𝑚𝑒𝑎𝑠𝑢𝑟𝑒𝑚𝑒𝑛𝑡 + 𝑁 ′ + 𝑠𝑒𝑒𝑑 (7)
After generating new report for current time slot, the report is buffered until the sending time is expired.

Generate shares and send report.

When the time is expired for buffered reports, the reports are divided according to SSS scheme as explained previous section.

𝑀 𝑖 𝑗′′ / 𝑤 (8)
Generated w shares are sent through other smart meters to the data collection center. With t shares, the data collection center could recover the report. Because of possible packet loss, all created w shares are sent to data collection center to be sure to obtain reports for each time slot. Choice of smart meters sent through depends on the routing protocol. In this paper, we assume that any multi hop routing protocol can be used.

Data Storage

The data collection center handles collecting reports from smart meters and storing them privately. Privacy preserving during data transmission is provided with above processes but at data center side it is still a problem. The center is responsible to keep reports in secure and response third party organizations for billing or load balancing.

To provide the security of data storage, the center stores shares in database without recovering the report and performing any other operations, when shares, which are sent from smart meter, reached the center. Hence, the center concatenates the report from related shares when the request message is arrived with the third party organization's predetermined private key. It provides that the private data are not directly accessible to other party organizations. The request responses are unique that it is represented measurement of one smart meter and time slot. Furthermore, the center concatenates shares to obtain the report. After report generation, the center gets time stamp of sending time, secret time slot information and measurement. Furthermore, timestamp represents sending time of the report that is perturbed in smart meter. To obtain the measurement time information, the delayed time represented N is used for the measurement time information.

𝑅𝑒𝑝𝑜𝑟𝑡: 𝐶𝑜𝑛𝑐𝑎𝑡𝑎𝑛𝑒𝑡𝑒(𝑡 𝑘 (𝑀 𝑖 𝑗′′ )) , k = 0 to t

Privacy Analysis and Simulation Results

In this section, we analyze the performance of the proposed protocol that provides privacy preserving communication in SGs. We first evaluate the protocol in terms of privacy, and then present the simulation results.

In this protocol, we provide privacy for consumers in two phases including packet encryption. First phase is temporal perturbation which prevents obtaining customers' daily life pattern from time-series data. In this phase, the generated measurement reports are sent out in batches with random delay to perturb the time. Randomization process is based on Laplace Distribution, so an attacker must repeat the same Laplace Distribution in order to eliminate time perturbation and obtain the ordered measurement reports. The security of Laplace Distribution process is show in [START_REF] Dwork | The Algorithmic Foundations of Differential Privacy[END_REF], hence the data stored in the data center is guaranteed to be secure against HBC type adversaries who can access the data.

Second phase is secret sharing during data communication. In order to cope with eavesdroppers who listen communication channels to obtain data of customers, generated reports are divided into parts according to SSS scheme and sent over different paths. The attacker must have at least t shares and selected prime numbers to recover original report. The security of SSS scheme is shown below. Similar proofs can be found in [START_REF] Bogdanov | How to securely perform computations on secret-shared data[END_REF][START_REF] Rosulek | Shamir Secret Sharing Alternative Proof of Security[END_REF][START_REF] Lipmaa | Lecture 9: Secret Sharing, Threshold Cryptography, MPC[END_REF] as well.

There are two distributions that are identical. There are specified 𝑡 constraints. The interpolation theorem tells us that there is only one polynomial 𝑓 satisfying all these constraints. Sharing process chooses a polynomial 𝑓 with 𝑓(0) = 𝑚 and selects the other 𝑘 -1 coefficients uniformly from 𝑍𝑝.

From the construction of polynomials, for all 𝑀 𝑖 𝑗′ , 𝑀 𝑖 𝑗′′ ∈ 𝑍𝑝, probabilities 𝑃𝑟[𝑚 =

𝑚'] are equal to 1/ 𝑝 𝑘-1 as it is the same distribution for all 𝑚. Hence, SSS scheme is perfectly secure and does not depend on the computational power of any party.

Simulation Results

In this section, we present simulation results and assess the performance of our proposed protocol. In our simulation model, we assume there are 100 smart meters, and 10 percent of smart meters are malicious nodes. We also assume the measurement reports are divided into 8 shares.

Simulation results for malicious node actions are analyzed by increasing threshold value, which is the number of shares to recover a report. The results show that as the threshold value increases the attackers' efficiency decreases, i.e., the number of obtained reports is diminished (shown as Fig 3). Therefore, in order to increase the privacy protection, threshold value t must be chosen as a large number. However, increasing t also increases the communication overhead of SG, therefore tradeoff between these two must be taken into account.

Fig. 3. Obtained reports possibility depend on threshold

We also evaluated the effect of the number of malicious nodes on the proposed protocol. Fig 4 . shows the obtained number of reports as the number of malicious nodes in the network increases. In order to avoid this situation, threshold value must be increased.

Fig. 4. Obtained reports depend on count of malicious nodes

It must be noted that the simulation results consider only SSS scheme and do not consider the time perturbation. Since our protocol employs time perturbation, even if an attacker obtains a complete report he/she will not be able to extract a meaningful data from that report as long as he/she does not have the secret information of the perturbation scheme.

Related Work

There are many studies that focus on the privacy issues of SG. In [START_REF] Rottondi | Privacy-Friendly Appliance Load Scheduling in Smart Grids[END_REF], SSS scheme is used for domestic appliances service request. In [START_REF] Rottondi | Privacy-Preserving smart metering with multiple data consumers[END_REF], customer data are collected via SSS scheme that encrypted data are aggregated with homomorphic properties. The communication channel is described as secure and the external attackers are not assumed to eavesdrop the communication channel. In [START_REF] Li | Secure Information Aggregation for Smart Grids Using Homomorphic Encryption[END_REF] and [START_REF] Li | Secure Information Aggregation for Smart Grids Using Homomorphic Encryption[END_REF] secure aggregation is used with homomorphic encryption that protects user privacy via aggregation tree.

Temporal privacy is also investigated by several papers in the literature. In [START_REF] Li | Privacy preservation in wireless sensor networks: A state-of-the-art survey[END_REF], dataoriented and context-oriented perturbation methods are presented. In [START_REF] Kamat | Temporal Privacy in Wireless Sensor Networks[END_REF], authors proposed a buffering and delaying protocol in WSNs to obfuscate the events that detected by sensors. This protocol preserves temporal privacy but could not be employed in SG in terms of importance of time information.

Conclusion

This paper proposes a privacy-preserving communication protocol for SG based on time perturbation and SSS scheme. The protocol provides privacy to participants and provides protection against attackers such as eavesdroppers, HBC or compromised intermediate smart meters. The proposed protocol ensures that intermediate smart meters and eavesdroppers could not infer any information from shares without least t shares of report. Even if a packet obtained by an attacker, he/she could not obtain measurement information that is measured by smart meter. With time perturbation operation, attackers could not access to daily life pattern through time-series measurement data because the sequence of data is changed for every measurement time.

To provide privacy of consumers to third-party organizations, the data center could only recover data to respond monthly bill requests or instant load balancing requests. Our work demonstrates that SSS scheme and time perturbation are the effective way for both consumer and data center privacy.

  𝑞 is the secret, 𝑞 is a prime number which is greater than possible secrets.  To divide the secret into 𝑤 shares, choose 𝑡 -1 integer random numbers 𝑝 1 , 𝑝 2 , … , 𝑝 𝑡-1 which are in range of [0, 𝑞 -1].  Compute the 𝑗 𝑡ℎ share(𝑥 𝑗 , 𝑦 𝑗 ), 1 ≤ 𝑗 ≤ 𝑤. 𝑥 𝑗 are integer distinct numbers like 1, 2, 3, ... and 𝑦 ≡ 𝑚 + 𝑝 1 𝑥 𝑠 + 𝑝 2 𝑥 𝑠 2 + ⋯ + 𝑝 𝑡-1 𝑥 𝑠 𝑡-1 mod 𝑞.  After the equations above, the secret can be recovered with t shares.


Zorana Bankovi¢
email: zorana.bankovic@imdea.org

Umer Liqat
email: umer.liqat@imdea.org

Pedro López-García

A Practical Approach for Energy Ecient Scheduling in Multicore Environments by combining Evolutionary and YDS Algorithms with Faster Energy Estimation

Keywords: Scheduling, energy eciency, multicore systems, evolutionary algorithms, YDS

Energy ecient scheduling and allocation in multicore environments is a well-known N P -hard problem. Nevertheless approximated solutions can be eciently found by heuristic algorithms, such as evolutionary algorithms (EAs). However, these algorithms have some drawbacks that hinder their applicability: typically they are very slow, and if the space of the feasible solutions is too restricted, they often fail to provide a viable solution. In this paper we propose an approach that overcomes these issues. The approach is based on a custom EA that is fed with predicted information provided by an existing static analysis about the energy consumed by tasks. This solves the time ineciency problem. In addition, when this algorithm fails to produce a feasible solution, we resort to a modication of the well-known YDS algorithm that we have performed, well adapted to the multicore environment and to the situations when the static power becomes the predominant part. This way, we propose a combined approach that produces an energy ecient scheduling in reasonable time, and always nds a viable solution. The approach has been tested on multicore XMOS chips, but it can easily be adapted to other multicore environments as well. In the tested scenarios the modied YDS can improve the original one up to 20%, while our EA can save 55 -90% more energy on average than the modied YDS.

Introduction

Energy ecient scheduling has gained a lot of interest in the recent past. A great number of publications, e.g., [START_REF] Marco | Analytic clock frequency selection for global DVFS[END_REF], try to present it as a mixed integer linear optimisation problem, which can be solved using mixed integer linear programming, or using a heuristic approach. However, these algorithms become impractical or fail to deliver a solution as the problem size grows. There is a signicant group of publications on using EAs for the problem of optimal scheduling and allocation in multiprocessor systems that allow Dynamic Voltage and Frequency Scaling (DVFS), e.g., the approach presented in [START_REF] Mezmaz | A bi-objective hybrid genetic algorithm to minimize energy consumption and makespan for precedence-constrained applications using dynamic voltage scaling[END_REF] aims to minimize both energy and makespan as a bi-objective problem.

Energy ecient scheduling and allocation in multicore environments with enabled DVFS is a well-known N P -hard problem. Nevertheless approximated solutions can be eciently found by heuristic algorithms, such as evolutionary algorithms (EAs). An example is our previous EA-based algorithm [START_REF] Bankovi¢ | Stochastic vs. Deterministic Evolutionary Algorithm-based Allocation and Scheduling for XMOS Chips[END_REF]. In our setting, we want to solve the general scheduling problem where the tasks have arbitrary release times and deadlines, and where preemption and migration of tasks are allowed, but the problem still remains N P -hard for arbitrary release times and deadlines of the tasks which are not agreeable, 1 as it was proven in [START_REF] Albers | Speed scaling on parallel processors[END_REF]. Our algorithm has been adapted for its application to multicore XMOS chips, but it could easily be adapted to any multicore environment, ranging from small scale embedded systems up to large scale systems, such as data centers. Its rst practical implementation relied on an existing analytical model for calculating the energy consumption for programs running on these chips [START_REF] Kerrison | Measuring and modelling the energy consumption of multithreaded, multi-core embedded software[END_REF]. The energy model is limited to the cases when all the cores belong to the same chip, thus they must run at the same voltage and frequency level at each moment. For this reason, in this work we solve the problem of the so-called global DVFS, when all the cores always have the same voltage and frequency.

In our rst EA implementation, each individual in each generation is evaluated by using the above mentioned program energy model, which requires the execution traces of the programs. Given that the traces can be huge, even for small programs, such evaluation introduces a huge amount of overhead. In order to overcome this issue, in this paper we use an existing static analysis which, at compile time, without the need of executing the programs, and in a few seconds, gives a safe estimation of the energy consumed by programs. As energy consumption often depends on (the size of) input data, which is not known at compile time, the static analysis provides the energy as a function of the input parameters, which is calculated once such input values are known at runtime. The energy consumption estimated by using the static analysis for a given scheduling is computed as the sum of the energies of the tasks running on dierent cores. This gives a safe upper bound on the total energy consumption, although it may be less precise than the estimations computed with the program energy model mentioned previously. This may reduce possible energy savings, nevertheless, the information it provides is still good enough to decide which scheduling is better, and the gain in speed of the algorithm is huge: the simulation time is reduced from a few hours to a few minutes.

However, the EAs can have trouble in nding a viable solution, in the sense that not all task deadlines are met, if the task deadlines are too tight. In order to overcome this problem, we have adapted the standard YDS 2 algorithm [START_REF] Yao | A scheduling model for reduced cpu energy[END_REF] (explained in Section 2) to multicore environments. As we will see later, our 1 Two tasks are agreeable if the task with later release time also has a later deadline. 2 The name is created using the rst letter of the authors' last names. experimental results show that if the EA nds a viable solution, it is better than the one obtained by our modied YDS algorithm in terms of energy savings.

For these reasons, the approach we propose (depicted in Figure 1) consists of the following steps:

1. Perform the static analysis of the input tasks to estimate the energy consumed by each of them. 2. Execute the EA using such estimations.

If the EA provides a viable solution, i.e., all the task deadlines are met, this is the nal solution.

Otherwise, execute our modied YDS algorithm and take its output as the nal solution.

We can distinguish two energy models in Figure 1:

Instruction-level energy model: gives an estimation of the energy consumed by the execution of a single instruction, which in general depends on its inputs, context, etc. However, for simplicity, the model assigns a constant value to each instruction. It is used by an abstract-interpretation based static analysis to infer the energy consumed by a program.

Program-level energy model: a formula that gives the energy consumption of the entire program, as presented in [START_REF] Kerrison | Measuring and modelling the energy consumption of multithreaded, multi-core embedded software[END_REF], which is used by our modied YDS algorithm.

In summary, in this paper we propose a time ecient scheduling approach, which always provides a viable energy ecient solution, and scales well as the input size grows. The main original contributions of our approach are:

1. The combination of an EA algorithm that resorts to a modied YDS algorithm, which always provides a viable solution. 2. Use of static analysis for energy estimation at compile time to guide the EA process, which results in signicant speed-up in solving the scheduling problem, and hence the practicability of our approach, while still providing a solution with signicant energy savings.

3. An improvement of the YDS algorithm, which eciently solves the static power issue in the situations the chip cannot be switched o (explained in Section 2).

The rest of the paper is organised as follows. Section 2 gives more details about our proposed approach. Section 3 presents an experimental evaluation of the approach. Finally, some conclusions are drawn in Section 4.

Our Proposed Approach

Evolutionary Algorithm The approach we propose is based on the NSGA-II [START_REF] Deb | A fast elitist multi-objective genetic algorithm: NSGA-II[END_REF] multiobjective evolutionary algorithm with two objectives: the execution time and the total energy consumption, where both should be minimised. The objectives are clearly in conict, since the application of DVFS reduces energy, but increases execution time. This justies the usage of a multiobjective algorithm. NSGA-II has been proven in the literature to perform good when the number of objectives is small [START_REF] Bankovi¢ | Stochastic vs. Deterministic Evolutionary Algorithm-based Allocation and Scheduling for XMOS Chips[END_REF], and in our case we have only two.

Since the output of the multiobjective approach is a set of solutions which form the (approximated) Pareto front, we can choose the solution that meets some given energy and/or time requirements. Usually we pick the solution with the minimal energy consumption among those that meet the given time bound (if applicable).

Individual Representation The problem that we are solving is the optimal (in terms of energy or time, possibly under some requirements involving them) allocation and scheduling of a set of tasks, where each task is dened by its: Unique ID. Release time, i.e., the moment when the task becomes available. Deadline, i.e., the latest moment when the task has to nish. Number of clock cycles, as a good approximation of the execution time.

Thus, the solution to this problem has to contain the following information:

The core(s) where each task will be executed. Since we allow task migration, a task can be allocated to more than one core. The current voltage and clock frequency (V, f ) state to exploit DVFS. The time periods when the tasks are executed. The number of clock cycles each task will execute in the dierent periods marked by the preemption and migration of that task. This allows to express the number of cycles a task will execute before it is preempted, as well as the number of cycles it will execute after it is resumed in the same or in a dierent core, etc.

Having in mind these requirements, we have designed the solution representation as shown in Figure 2, which does not introduce signicant overhead when executing the EA. Any given task has a positive (unique) number as its ID.

Each gene representing a task ID is followed by a gene representing the number of cycles of the task that will be executed without any preemption. The order of task IDs represents the order of their temporal execution. We also use negative two digit numbers to encode the spatial allocation of the tasks. The rst digit represents the core where the tasks are being executed and the second one an encoding of the (V, f ) state of that core. As it will be explained in Section 3, Table 1, the number of dierent cores and states is nite, as well as the number of their combinations. The tasks following the allocation code are executed on that coded location. For instance, on Figure 2 we read: on core 1 in state 2, 48 cycles of task 1 will be executed, and 77 cycles of task 5, in this exact order, etc.

Our approach allows a random order of allocation codes, in order to solve the most general problem. However, if two consecutive allocation codes have dierent (V, f ) states, this means that the tasks allocated to the cores they represent will not be executed in parallel, since all of the cores have to run at the same (V, f ) at any moment, and thus the (V, f ) state has to be changed before the second group of tasks is executed. For example, in Figure 2 the allocation code following -12 is -24, which means that the chip will be rst in the (V, f ) state 2 and all tasks allocated to core 1 will be executed sequentially on that core. After they nish their execution, the core will change its (V, f ) state from 2 to 4, and the tasks allocated to core 2 will be executed sequentially on core 2. If the second allocation code were -22 instead of -24, then the (V, f ) state would not change, and the tasks allocated on cores 1 and 2 would be executed in parallel. Population Initialisation Individuals in the initial population are created by randomly assigning tasks to random cores in random (V, f ) settings with equal probability. However, in order to provide a load balanced solution (as much as possible), the probability of choosing a core decreases as its load increases, which is given by the following formula:

P rob = 1 N umberOf Cores - CurrentCoreLoad T otalLoad (1)
where CurrentCoreLoad stands for the current load of the core expressed as the number of cycles, while T otalLoad stands for the total number of cycles of all the tasks on all cores. According to the formula, at the beginning of the initialisation process, all cores have the same probability of being chosen, while this probability decreases as the core becomes loaded, and is close to 0 when the load reaches the state where it is equally distributed in all cores. If during the initialisation process a newly calculated probability value of a core is below 0, the value is rounded to 0, and no new load will be assigned to that core. These random solutions do not always have to provide a viable solution, i.e., some of the tasks might miss their deadlines. For this reason, the mutation operator and the objectives are designed to deal with this problem.

Solution Perturbations Given the unique nature of the individual representation, we have designed new crossover and mutation operators. The individuals that participate in the crossover are selected by using the standard tournament selection process.

The Crossover Operator Since our solution allows task migration, a given task ID can appear more than once, so we cannot apply any of the existing permutation-based crossover operators. Thus, we have designed our own operator, where each child will preserves the order of genes representing the task allocation and scheduling from one parent, and only the genes representing the number of cycles can be taken from the other parent. In this way, the produced ospring is a combination of both parents, and is at the same time a viable solution to the problem. The process is depicted in Figure 3 for the most simple case of 2 cores, 2 tasks and 2 states, with one possible output. We can observe that the rst child, C1 takes the scheduling and allocation from the rst parent, P 1, and the cycle distribution from the second parent, P 2, while the second, C2, takes the scheduling and allocation from P 2 and the cycle distribution from P 1. The Mutation Operator The mutation operator can perform dierent actions involving one or two tasks. Consider two tasks i and t. When choosing the rst one we give higher probability to the tasks which miss their deadlines, in order to achieve a viable solution as soon as possible. In each generation we perform either one of the following operations with the same probability (depicted in Fig. 4):

Swapping: i and t, together with their corresponding number of cycles, exchange their positions in the solution. However, in order to avoid creating solutions which are not viable, i and t have to belong to the cores that are executed in parallel (dened with consecutive allocation codes that are in the same (V, f ) state). In Fig. 4 we can observe that tasks 1 and 2 are swapped between cores 1 and 2, and both cores are in state 1.

Moving : move i to a random position j. For the same reason as before, the position j has to belong to a core being executed in parallel as i's. In Fig. 4 we can observe that the rst part of task 1 (40 cycles) is moved to core 2, before task 2.

Changing the number of cycles: assigns a dierent number of cycles to all the appearances of task i, in a way the total number of cycles of that task remains the same. In Fig. 4 we can observe that task 1 after the change executes 25 cycles on core 1 in state 1 and 45 cycles on core 2 in state 2.

Objective Functions

Execution Time One objective of our optimisation problem is to minimize the total execution time of the schedule, which is the time spent since the rst task starts its execution until the last task nishes its execution. However, since the initial population is randomly created, it is possible that some of the tasks miss their deadlines, making the solution unviable. Assuming that these solutions can provide some quality genetic material, we do not want to discard them completely, but we penalize them by adding the amount of time the tasks have missed their deadlines to the objective function. Thus, the time objective function for n cores and k dierent tasks is the following:

T = T + 1≤i≤n ( 1≤j≤k x i,j • y j • (s i,j + τ i,j -deadline j )) (2)
where T is the total execution time, given by:

T = max 1≤i≤n 1≤j≤k (x i,j • (s i,j + τ i,j )) -min 1≤i≤n 1≤j≤k (x i,j • s i,j ) (3)
where s i,j ≥ 0 is the moment when task j is scheduled on core i (s i,j = 0 if the task j is not scheduled on core i), τ i,j is the execution time of task j on core i, x i,j is a binary value, that represents whether the task j is executed on core i (x i,j = 1) or not (x i,j = 0). The second part of formula (2) represents the penalisation, where y j is another binary value that expresses whether the task j has missed its deadline, deadline j , (y j = 1) or not (y j = 0).

Energy Consumption This objective represents the total energy consumption of the given schedule. In the most general case it is given by the following formula:

E = 1≤i≤n (P st,i • T + 1≤j≤k (x i,j • p i,j • τ i,j )) (4)
where P st,i is the static power of core i, T is the same as in formula (3), p i,j is the dynamic power of task j when executed on core i, and x i,j and τ i,j are the same as in formula [START_REF] Bankovi¢ | Energy Ecient Allocation and Scheduling for DVFS-enabled Multicore Environments using a Multiobjective Evolutionary Algorithm[END_REF]. In this work we use static analysis to estimate the energy of single tasks, which will be explained in more detail in Section 2, while the energy is the sum of the energies of all the tasks, as given in formula (4).

The Modied YDS Algorithm YDS [START_REF] Yao | A scheduling model for reduced cpu energy[END_REF] is a well known algorithm for energyecient scheduling for single core DVFS-enabled environments. We have chosen it because it always nds a feasible (and optimal) solution that minimises the total energy consumption, it is simple and fast. However, it does not take into account the static power, which nowadays forms an important part of the total power. YDS reduces the frequency and voltage in order to minimise the dynamic power in a way that the execution time of tasks are extended to their deadlines. However, this also results in an increase of the static energy. Thus, there is a critical point from which further reduction of voltage and frequency actually starts increasing the energy consumption. Attempts to reducing the static energy when applying DVFS have mainly tried to group the inactive periods by moving task executions towards their deadline or its release point, and turning o the chip during such periods [START_REF] Jejurikar | Leakage aware dynamic voltage scaling for real-time embedded systems[END_REF]. However, this is not always possible since chip wake up can take more time than available. Our alternative proposal does not turn o the chip, but instead, nds the critical (V, f ) point below which further decrease is not benecial. Thus, our modication of YDS consists of the following steps: 1. In order to decide if it is benecial to further decrease the frequency, and in this way avoid the problem when the increase in energy consumed by the static power leads to the total energy increase, we use the simple slope-based method presented in [START_REF] Miyoshi | Critical power slope: Understanding the runtime eects of frequency scaling[END_REF]. 2. In order to support the multicore system, we propose two dierent heuristics for allocating the tasks to dierent cores: the load balanced solution and the solution where a task allocation leads to the minimal frequency increase.

After the allocation, the YDS algorithm is applied to each core.

3. In order to adapt the algorithm so that it assigns only the frequencies supported by the system, we propose to divide the computational load into two parts and execute them on two supported frequencies in a way the total execution time remains (almost) the same.

A Solution to the Static Power Issue of YDS As mentioned before, we use the simple slope-based method presented in [START_REF] Miyoshi | Critical power slope: Understanding the runtime eects of frequency scaling[END_REF]. The only requirement for its application is the availability of a power model where the static and dynamic power are separated. The main idea of the method is the following. If we x the voltage, the power is a linear function of the frequency, and after applying some simple numeric transformations (explained in [START_REF] Miyoshi | Critical power slope: Understanding the runtime eects of frequency scaling[END_REF] in detail), it can be expressed in this way:

P f = P fmin + m • (f -f min )
(5) where P f denotes that the power depends on f , and f min is the minimal possible frequency, assumed to be the one which permits the execution of tasks to nish at their deadlines, and m is called the slope of the power function. If we can compare energies at dierent frequencies, we will know if it is energy ecient to decrease the frequency. Since power is a function of m, the same applies to energy, and thus it determines the decision of decreasing the frequency. In theory, there should exist a slope at which energy is equal for all frequencies. This slope is called the critical power slope :

m critical = P fmin -P idle f min (6)
If the actual slope m (calculated from Eq. 5) is greater than the critical one then we can decrease the frequency in order to save energy. However, if the slope is lower than the critical one, then the frequency should be increased in order to save the energy. Since the voltage can also change, the slope should be calculated for each (V x , f x ) point for each frequency f x :

m fx critical = P fx -P idle f x (7)
This value is then compared with the actual slope m fx at each (V x , f x ) point with frequency f x . Again, if m fx > m fx critical , we should decrease the frequency in order to save the energy. However, if m fx < m fx critical , the frequency should be increased in order to save the energy.

Optimal Task-Core Allocation Other aspects of adapting YDS to a multicore environment consist of nding an optimal number of cores, and allocation, i.e., assignments of tasks to cores. In this work we do not deal with the rst part of the problem, we just show that an optimal number of cores exists and it is not necessarily equal to the maximal possible number of cores. Regarding the second part, we have tested two possibilities:

1. Assign tasks to the cores so that the load is equally distributed between them.

2. Assign tasks in the way its addition assumes minimal change in frequency: a task t is assigned to the core with minimal activity, measured as the number of clock cycles in its active period (t_release_time, t_deadline).

Assigning Frequencies Supported by the System If the frequency f calculated by YDS is not supported by the system, the total number of cycles ω i is divided in two parts, ω i1 and ω i2 , which are executed on two frequencies f 1 and

f 2 (f 1 ≤ f ≤ f 2 )
supported by the underlying system. The values of ω i1 and ω i2 are calculated by solving the following system of equations:

ω i f ≈ ω i1 f 1 + ω i2 f 2 ω i = ω i1 + ω i12 (8)
Energy Static Analysis as Input In order to estimate the energy consumed by programs without actually running them we use an existing static analysis. It is a specialisation of the generic resource analysis presented in [START_REF] Serrano | Resource Usage Analysis of Logic Programs via Abstract Interpretation Using Sized Types[END_REF] for programs written in a high-level C-based programming language, XC [START_REF] Watt | Programming XC on XMOS Devices[END_REF], running on the XMOS XS1-L architecture, that uses the instruction-level energy cost models described in [START_REF] Liqat | Energy Consumption Analysis of Programs based on XMOS ISA-level Models[END_REF]. The analysis is general enough to be applied to other programming languages and architectures (see [START_REF] Liqat | Energy Consumption Analysis of Programs based on XMOS ISA-level Models[END_REF][START_REF] Liqat | Inferring Energy Consumption at Dierent Software Levels: ISA vs. LLVM IR[END_REF] for details). It enables a programmer to symbolically bound the energy consumption of a program P on input data x without actually running P (x). It is based on setting up a system of recursive cost equations over a program P that capture its cost (energy consumption) as a function of the sizes of its input arguments x. Consider for example the following program written in XC: i n t f a c t ( i n t N) { i f (N <= 0) r e t u r n 1 ; r e t u r n N * f a c t (N -1 ) ; } The transformation based analysis framework of [START_REF] Liqat | Energy Consumption Analysis of Programs based on XMOS ISA-level Models[END_REF][START_REF] Liqat | Inferring Energy Consumption at Dierent Software Levels: ISA vs. LLVM IR[END_REF] would transform the assembly (or LLVM IR) representation of the program into an intermediate semantic program representation (HC IR), that the analysis operates on, which is a series of connected code blocks, represented as Horn Clauses. The analyser deals with this HC IR always in the same way, independent of where it originates from, setting up cost equations for all code blocks (predicates).

f act e (N ) = f act_if e (0 ≤ N, N ) + c entsp + c stw + c ldw + c ldc + c lss + c bf f act_if e (B, N ) =    f act e (N -1) + c bu + 2 c ldw + c sub + + c bl + c mul + c retsp if B is true c mkmsk + c retsp if B is false
The cost of the function f act is captured by the equation f act e which in turn depends on the equation f act_if e , that captures the cost of the two clauses representing the two branches of the if statement, and a sequence of low-level If we assume (for simplicity of exposition) that each instruction has unitary cost in terms of energy consumption, i.e., c i = 1 for all i, we obtain the energy consumed by fact as a function of its input data size (N): f act e (N ) = 13 N + 8.

The functions inferred by the static analysis are arithmetic functions (polynomial, exponential, logarithmic, etc.) that depend on input data sizes (natural numbers). We use them in our scheduling and allocation algorithm to estimate the energy consumed by the dierent tasks involved. Such estimation can be computed very eciently once the input data sizes of the tasks are known, since all the basic arithmetic functions involved can be evaluated in little bounded time.

Experimental Evaluation

XMOS Chips In this work we target the XS1-L architecture of the XMOS chips as a proof of concept. Although these chips are multicore and multithreaded, in this work we assume a single core architecture with 8 threads, which is the architecture for which we have an available energy model. In this case, we can use the algorithm and representation of individuals described in Section 2 by considering that a thread in our experiments is conceptually equivalent to a core executing tasks sequentially, as described previously. We refer the reader to [START_REF] Bankovi¢ | Stochastic vs. Deterministic Evolutionary Algorithm-based Allocation and Scheduling for XMOS Chips[END_REF] for a description of a representation of individuals whose allocation codes include three digits, representing: a core, a thread running in parallel on that core, and a (V, f ) state.

In the XS1-L architecture, the threads enter a 4-stage pipeline, meaning that only one instruction from a dierent thread is executed at each pipeline stage. If the pipeline is not full, the empty stages are lled with N OP s (no operation). Eectively, this means that we can assume that the threads are running in parallel, with frequency F/N , where F is the frequency of the chip, and N = max(4, numberOf T hreads). DVFS is implemented at the chip level, which means that all the cores have the same voltage and frequency at the same time. In order to apply DVFS, we need a list of Voltage-Frequency (V, f ) pairs or ranges that provide a correct chip functioning. We have experimentally concluded that the XMOS chips can function properly with the voltage and frequency levels given in Table 1.

Task Set In order to test our proposed approach, we use two dierent groups of task sets. The rst group is made up of small tasks, where the EA training with the program-level energy model takes around one day to complete. This group is used to show the dierence between the results obtained with the EA trained with the program-level energy model and the EA trained with the energy estimations obtained by the static analysis. In this group, we use four dierent arithmetic programs: fact(N), for calculating the factorial of N, fibonacci(N) for calculating the Nth Fibonacci number, sqr(N) for computing N 2 and power_of_two(N) for computing 2 N . In total, we have created a set of 22 tasks to be scheduled, corresponding to the execution of the previous programs with dierent inputs N .

In the second group we use real world programs, where the EA training based on the program-level energy model is not practical: fir(N), i.e., Finite Impulse Response (FIR) lter, which in essence computes the inner-product of two vectors of dimension N, a vector of input samples, and a vector of coecients, and biquad(N), which is a part of an equaliser implementation, based on a cascade of Biquad lters, whose consumed energy depends on the number of banks N. We have used four dierent FIR implementations, with dierent number of coecients: 85, 97, 109 and 121. Furthermore, we have used four implementations of the biquad benchmark, with dierent number of banks: 5, 7, 10 and 14. We have tested our approach in scenarios of 16 and 32 tasks, each one corresponding to such implementations. The tasks corresponding to the same implementation have dierent release times and deadlines.

The energy consumed by the programs is inferred at compile time by the static analysis described in Section 2. Such energy is expressed as a function of a parameter N , the size of the input, which is only known at runtime. Such functions are given in Table 2 for 3 of the 6 dierent voltage and frequency levels used in this work (for conciseness, as the functions for each program have the same complexity order, but dierent coecients). The static analysis assumes that a single program (task) is running on one thread on the XMOS chip, while all other threads are inactive. In this implementation, the EA algorithm approximates the total energy of a schedule by adding the energies of all the tasks. Although in this way we loose precision, the estimation still provides precise enough information for the EA to decide which schedule is better. Testing Scenarios In our current implementation, we assume no dependency between the tasks since it is not supported by the available energy models. The release times and deadlines of the dierent tasks are set in dierent scenarios in order to experimentally show the benets of DVFS and optimal scheduling, where all the tasks have dierent release times and deadlines, with tighter deadlines; and that it is important to take into account the static power, especially in the case of loose deadlines. Scenario 1: Tasks with Loose Deadlines In this scenario the release time of a task k, denoted T k rel is a random moment between 0 and the total execution time at the maximal frequency of all the tasks executed sequentially on a single core. Also, the deadline of a task is a random moment between T k rel + 10 × T k maxf and T k rel + 20 × T k maxf , where T k maxf denotes the execution time of the task at maximum frequency. This way we achieve a scenario with loose deadlines even at a smaller frequency. Scenario 2: Tasks with Tight Deadlines Here, the release time is the same as in Scenario 1. However, the deadline of a task is a random moment between

T k rel + 5 × T k maxf and T k rel + 7 × T k
maxf . This way we get tighter deadlines, but also provide a set of tasks which are schedulable on the given platform. Note that the deadlines become even tighter as the frequency decreases. Results: The Improved YDS Table 3 shows the savings of our improved YDS algorithm presented in Section 2 compared to the original YDS, for dierent number of cores and for two dierent ways of task allocation: Alloc. 1, where the load is evenly distributed between the cores, and Alloc. 2, where the addition of a task implies a minimal increase in the frequency. The energy saving resulting from a particular scheduling is calculated using the following formula:

Y DS_original -Y DS_modif ied Y DS_original • 100 (9)
In Table 3, energy savings are achieved in all, but in two cases with tight deadlines. A possible reason could be the fact that all the threads need to have the same frequency always, which means that the maximal necessary frequency is assigned, which is not necessarily be optimal for all the threads. However, in the case of loose deadlines, the savings are much more signicant, since the static power plays a more important role.

Results: EA vs. improved YDS Both EA and YDS are implemented in C++. EA extends the MOGAlib library [START_REF] Gaál | Multi-Level Genetic Algorithms and Expert System for Health Promotion[END_REF] for multiobjective genetic algorithms. In the EA, the population of 200 individuals is evolved for 150 generations. The probability of both crossover and mutation is 0.9. The mutation is assigned higher probability than usual due to its important role for reaching a viable solution. Since the result of the optimisation process is a set of possible solutions which form the approximated Pareto front, we take the solution with minimal energy consumption where all task deadlines are met. Table 4 presents results comparing the EA trained with the energy estimations provided by static analysis, versus the improved YDS algorithm presented in Section 2. In the rst column, the energy of the nal solution calculated by using the program-level energy model is given (EA s ). The second column gives the energy of the nal scheduling obtained by the modied YDS algorithm (referred as Y DS m ) using the program-level energy model, while the third column gives the energy saving of the EA trained with static analysis compared to YDS. Finally, the last column shows the energy saving obtained with the EA trained with the program-level energy model(EA m ) [START_REF] Bankovi¢ | Energy Ecient Allocation and Scheduling for DVFS-enabled Multicore Environments using a Multiobjective Evolutionary Algorithm[END_REF], which is only applicable in the scenarios with a small number of numeric tasks. Each row shows statistics for each scenario taken from 10-20 runs of the algorithm for the same scenario, where CI0.01 and CI0.05 represent 99% and 95% condence intervals, meaning that we can claim with 99(95)% certainty that the nal result will fall in these intervals.

In order to perform the comparison between the EA and Y DS m , in the case of EA and tight deadlines, we present the results when the EA can nd a viable solution. However, the EA does not always provide a viable solution in all the scenarios with tight deadlines created as explained previously. As we can see in Table 4, if the EA nds a viable solution, it always performs better than the Y DS m . We can also observe that the EA trained with the programlevel energy model achieves better results. However, the EA trained with the energy estimations by static analysis still achieves very good results, but with the training process that lasts around 10 minutes, compared to around 24 hours of training the EA with the program-level energy model, which makes it much more practical.

Conclusions and Future Work

In this work we propose a holistic approach for optimal scheduling, allocation and voltage(V ) and frequency(f) assignment in multicore environments, adapted for multicore XMOS chips. The main part of our approach is based on our custom developed EA-algorithm, which relies on static analysis to eciently estimate the energy of input tasks. The use of such static analysis improves signicantly the speed of the EA training process, thus allowing its real world applicability. Furthermore, since in the case of very tight task deadlines the EA can fail in providing a feasible solution, we add one more stage based on the YDS algorithm, which has been adapted for multicore environments and improved in a way it takes into account the static power. In this way, we have developed an ecient approach which is capable of providing energy savings in each possible scenario.

However, although the use of static analysis based estimations still provides energy savings, we have seen that better results can be achieved with more precise energy estimations. For this reason, we plan to use a static analysis of the energy consumed by concurrent programs, which is expected to provide additional savings. As a future work, we also plan to study the eect of dierent versions of the crossover and mutation operators in dierent situations, which could enable adding a heuristic for choosing the optimal version for each possible scenario.

Fig. 1 .

 1 Fig. 1. Overview of our scheduling approach.

Fig. 2 .

 2 Fig. 2. An example of (part of) a solution (i.e., individual) representation.

Fig. 3 .

 3 Fig. 3. An example of a crossover operation.

Fig. 4 .

 4 Fig. 4. Examples of mutation operations.

Table 1 .

 1 Viable (V, f ) pairs for XMOS chips. The cost of low-level instructions, which constitute an energy cost model, is represented by c i where i ∈ {entsp, stw, ldw, ...} is an assembly instruction. Such costs are supplied by means of assertions that associate basic cost functions with elementary operations.

	V oltage(V )	0.95	0.87	0.8	0.8	0.75	0.7
	f requency(M Hz)	500	400	300	150	100	50
	instructions.

Table 2 .

 2 Energy functions inferred by static analysis for 3 dierent pairs of voltage (V)/frequency (MHz).

		V = 0.70	V = 0.75	V = 0.80
		F = 50	F = 100	F = 150
	f act(N )	60.5 N + 46	35 N + 26.7	27 N + 20.5
	f ib(N )	87.19 × 1.62 N +	50.32 × 1.62 N +	38.68 × 1.62 N +
		26.7 × (-0.62) N -74.7 15.44 × (-0.62) N -43 11.85 × (-0.62) N -33.2
	sqr(N )	21.3 N 2 + 121 N	12.3 N 2 + 69.8 N	9.48 N 2 + 53.7 N
		+39.1	+22.5	+17.3
	powerOf 2(N )	55.1 × 2 N -39	63.7 × 2 N -39	24.49 × 2 N -30
	f ir(N )	74.93 N + 124.5	43.36 N + 71.9	33.41 N + 55.2
	biquad(N )	386 N + 128	223.6 N + 74.2	172.5 N + 57.2

Table 3 .

 3 Energy savings obtained by the modied YDS vs. the original YDS (%).

		Tight deadlines Loose deadlines
	#Cores Alloc. 1 Alloc. 2 Alloc. 1 Alloc. 2
	1	4.18	4.18	6.21	6.21
	2	1.5	4.26	14.67	14.67
	3	-5.26	3.17	14.67	14.67
	4	2.22	2.77	8.8	8.8
	5	-3.28	3.47	11.18	11.18
	6	0.95	4.34	11.82	11.82
	7	4.8	3.03	10.9	10.9
	8	19.36	5.61	10.56	10.56

Table 4 .

 4 EA vs. improved YDS in dierent scenarios.

		EAs(µJ) Y DSm(µJ) (Y DSm-EAs) Y DSm	(%)	(Y DSm-EAm) Y DSm	(%)
		A scenario with 22 small numeric tasks and loose deadlines
	Mean	14.3	33.1	56.8		76.57
	CI 0.01 11.6 -17	NA	48.64 -64.95	67.87-85.27
	CI 0.05 12.2 -16.4	NA	50.45 -63.14	70.05-83.09
		A scenario with 22 small numeric tasks and tight deadlines
	Mean	14.6	34.8	60.92		69.83
	CI 0.01 11.5-17.7	NA	49.14 -66.95	57.18-57.18
	CI 0.05 12.2 -17	NA	51.15 -64.94	60.34-60.34
	A scenario with 16 tasks made of Biquad and FIR lters and loose deadlines
	Mean	4.38	35.3	87.59		NA
	CI 0.01 3.4 -5.3	NA	85 -90.37	NA
	CI 0.05 3.7 -5.1	NA	85.55 -89.52	NA
	A scenario with 16 tasks made of Biquad and FIR and tight deadlines
	Mean	14.5	35.4	59.04		NA
	CI 0.01 9.4-19.6	NA	44.63 -73.45	NA
	CI 0.05 10.6 -18.4	NA	48.02 -70.06	NA
	A scenario with 32 tasks made of Biquad and FIR lters and loose deadlines
	Mean	17.85	68.16	73.81		NA
	CI 0.01 10.8 -25	NA	63.32 -84.15	NA
	CI 0.05 12.5 -23.3	NA	65.82 -81.66	NA
	A scenario with 32 tasks made of Biquad and FIR lters and tight deadlines
	Mean	29.43	68.16	56.82		NA
	CI 0.01 0.72 -51.6	NA	24.3 -89.44	NA
	CI 0.05 12.5 -46.3	NA	32.07 -81.66	NA

Acknowledgements. The research leading to these results has received funding from the European Union 7th Framework Programme under grant agreement 318337, ENTRA -Whole-Systems Energy Transparency, Spanish MINECO TIN'12-39391 StrongSoft and TIN'08-05624 DOVES projects, and the Madrid M141047003 N-GREENS project.