Pietro Aricò
email: pietro.arico85@gmail.com

Gianluca Borghini

Dr Ilenia Graziani

Jean-Paul Imbert

Géraud Granger

Raïlane Benhacene

Dr Simone Pozzi

Dr Linda Napoletano

Gianluca Di Flumeri

Prof Alfredo Colosimo

Ing Gianluca Borghini

Ing Géraud Granger

Ing Railane Benhacene

Ing Gianluca

Di Flumeri

Prof Fabio Babiloni

Ing Pietro Aricò

Air-traffic-controllers (ATCO): neurophysiological analysis of training and workload

Keywords: EEG, ECG, EOG, cognitive learning, training assessment, mental workload, ATM, ATCO

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

INTRODUCTION

Controlled airspace is divided into sectors. An en route sector is a region of airspace that is typically situated at least 50 km from an airport for which an associated ATCO has responsibility. ATCOs have to accept aircraft into their sector; check aircraft, issue instructions, clearances, and advice to pilots and hand aircraft off to adjacent sectors or to airports. When the aircraft leaves the airspace assigned to the ATCO, control of the aircraft passes onto ATCO controlling the next sector (or to the tower ATCO). As is typical in many real-world complex systems, this environment imposes multiple concurrent demands on the operator, in fact in the en route air traffic control environment, the system that confronts the air traffic controller comprises a large number of aircrafts coming from a variety of directions, at diverse speeds and altitudes, heading to different destinations [START_REF] Gronlund | Role of memory in air traffic control[END_REF]. ATCOs have two main goals. The primary goal is to ensure that aircraft under jurisdiction adhere to International Civil Aviation Organization (ICAO) mandated separation standards. For example, one of the most common separation standards requires that aircraft under radar control be separated by at least 1,000 feet vertically and 5 nautical miles horizontally. The secondary goal is to ensure that aircraft reach their destinations in an orderly and expeditious manner. These goals require the ATCO to perform a variety of tasks, including monitoring air traffic, anticipating loss of separation (i.e., conflicts) between aircraft, and intervening to resolve conflicts and minimize disruption to flow. (For an extensive compilation of the tasks and goals of en route control, see [START_REF] Rodgers | Conversion of the CTA, Inc., en route operations concept database into a formal sentence outline job task taxonomy[END_REF]). Total world airline scheduled passenger traffic in terms of passenger-kilometers is projected to grow at an annual average rate of 4.4% over the period 2002 to 2015, according to forecasts prepared by the ICAO (2004). To accommodate predicted traffic growth there is a need to increase en route airspace capacity through the introduction of new air traffic management systems, controller tools and procedures. But, the consensus among research and operational communities is that it is really important to understand the factors that drive mental workload if they are to improve airspace capacity [START_REF] Christien | Air traffic complexity indicators ATC sectors classification[END_REF], [START_REF] Majumdar | The Factors Affecting Airspace Capacity in Europe: A Cross-Sectional Time-Series Analysis Using Simulated Controller Workload Data[END_REF]. Most research has focused on identifying characteristics of the air traffic picture that create task demand for ATCOs (e.g., [START_REF] Grossberg | Quarterly Rep. of the FAA Office of Air Traffic Evaluations and Analysis)[END_REF]- [START_REF] Manning | Investigating the Validity of Performance and Objective Workload Evaluation Research[END_REF]). Others argue that there is no simple linear relationship between task demand and workload (e.g., [START_REF] Athènes | Complexity and controller workload: Trying to bridge the gap[END_REF], [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF]). Several current research groups agree with Sperandio's [START_REF] Sperandio | Variation of Operator's Strategies and Regulating Effects on Workload[END_REF] view that a relationship between task demand and workload can be better understood by considering how ATCOs use strategies to manage their resources and regulate their workload [START_REF] Majumdar | The Factors Affecting Airspace Capacity in Europe: A Cross-Sectional Time-Series Analysis Using Simulated Controller Workload Data[END_REF], [START_REF] Athènes | Complexity and controller workload: Trying to bridge the gap[END_REF], [START_REF] Averty | Evaluating A New Index Of Mental Workload In Real ATC Situation Using Psychophysiological Measures[END_REF]- [START_REF] Histon | The Impact of Structure on Cognitive Complexity in Air Traffic Control[END_REF]. Factors such as skills, training, experience, fatigue and other "stressors" all mediate the relationship between task demands, safety and performance of the ATCO. Hence, it is easy to understand how quantitative information about skills level and mental states could help to evaluate the ATCO's workload level and to decide if they might need more training before working into real environments. Several studies described the correlation of spectral power of the EEG bands with the complexity of the task that the subjects are performing [START_REF] Borghini | Measuring neurophysiological signals in aircraft pilots and car drivers for the assessment of mental workload, fatigue and drowsiness[END_REF]. In fact, an increase of electroencephalographic power spectral density especially over the frontal cortex in the theta band (4 -7 Hz) and an EEG PSD decrease in the alpha band (8 -12 Hz) over the parietal cortex have been observed when the required mental workload, the task's complexity and the amount of information to be processed increase. Furthermore, it has been suggested that also an increased Heart Rate (HR) could be related with an increased mental workload and engagement. On the contrary the eyesblink duration and rate are inversely correlated with the increase of the mental workload and attention levels [START_REF] Borghini | Measuring neurophysiological signals in aircraft pilots and car drivers for the assessment of mental workload, fatigue and drowsiness[END_REF]. The hypotheses of the study are that i) as the EEG theta PSD increases and the EEG alpha PSD decreases with these cognitive phenomena , at the end of the training period such PSDs increment and decrement should be different from the beginning of the training period , therefore such trends could be taken as indexes of the correct acquisition of procedural skills and of less request of cognitive resources for the correct execution of the task, ii) the combined use of EEG features and HR can improve the reliability of the measure with respect of using the single information. Also, these biosignals can be used for the real-time evaluation of the operator's workload level, in a real Air Traffic Management (ATM) scenarios. Such hypotheses have been tested on a group of 6 subjects who succeeded in the 5-days-training-period and who were asked to execute the experimental task for two more weeks in order to evaluate their mental workload under three different difficulty levels.

MATERIALS AND METHODS

Experimental subjects and ATM simulation task

A group of six healthy volunteers has been selected in terms of age (21±4 years) and previous computer game skills. The subjects have been asked to learn how to execute correctly an ATM task (labyrinth, LABY), that never did before, under easy (E), medium (M) and hard (H) conditions, randomly selected and proposed to them in order to avoid any habituation and expected effects. The LABY microworld is a functional simulation of Air Traffic Control (ATC), provided by the software engineers of the French ENAC, that captures the underlying processes involved in electronic air traffic management with a simplified version of the operational human-machine interface. Microworlds are computerbased human-in-the-loop simulation environments that offer testing, behavioral/physiological measurement, and training capabilities, with the flexibility to build various scenarios [START_REF] Brehmer | Experiments with computer-simulated microworlds: Escaping both the narrow straits of the laboratory and the deep blue sea of the field study[END_REF], [START_REF] Gonzalez | The Use of Microworlds to Study Dynamic Decision Making[END_REF]. LABY is a dynamic environment whereby a controller must issue directional commands to guide aircrafts along a predetermined route, whilst avoiding potential conflicts and dealing concurrently with other incoming information. The LABY microworld is based upon the main task of guiding N plane(s) around a predetermined route, indicated by a green path (Fig. 1). Participants must input numerical values such as heading, flight level, speed, etc., in order to direct flight around the trajectory and to avoid any conflicts or obstacles which may occur during the flight-route. Penalties are applied if the aircrafts deviate off the route or if other constraints are not met. The difficulty of the task can be altered according to how many aircrafts the participant have to control, the number and type of clearances required over the time and the number/trajectory of other interfering flights. In the first week, the subjects trained for 5 consecutive days (SESSIONS T1÷T5) and their neurophysiological signals have been recorded in the first (T1), in the third (T3) and in the fifth (T5) session. The behavioral (task performance) and subjective workload perception data have been collected every day. After the training period, the subjects were asked to execute the LABY on two consecutive days in the second week and on a day after a week since the last experimental session. At the end of each experimental condition, the subjects filled the NASA-Task Load indeX (TLX, [START_REF] Hart | Development of NASA-TLX (Task Load Index): Results of Empirical and Theoretical Research[END_REF]) questionnaire for the evaluation of the perceived workload of the proposed task.

Acquisition of the brain activity and of the physiological signals: Electroencephalogram (EEG) and physiological signals, including vertical electrooculogram (EOG) and electrocardiogram (ECG)

, have been recorded by the digital monitoring BEmicro system (EBNeuro system). The sixteen EEG channels (FPz, F3, Fz, F4, AF3, AF4, C3, Cz, C4, P3, Pz, P4, POz, O1, Oz and O2), the ECG and the EOG channels have been collected simultaneously with a sampling frequency of 256 (Hz). All the EEG electrodes have been referenced to both earlobes, and the impedances of the electrodes were kept below 10 (kΩ). The bipolar electrodes for the heart activity have been placed on the Erb's point, while the bipolar electrodes for the EOG have been positioned vertically on the left eye.

EEG analysis:

The acquired EEG signals have been digitally band-pass filtered by a 4 th order Butterworth filter (low-pass filter cut-off frequency: 30 (Hz), high-pass filter cut-off frequency: 1 (Hz)) and then segmented in epochs of 2 seconds, 0.125 secondsshifted. The EOG signal has been used to remove eyes-blink artefacts from the EEG data by using the Gratton method [START_REF] Gratton | A new method for off-line removal of ocular artifact[END_REF]. The EEG PSD has then been estimated by using the Fast Fourier Transform (FFT) in the EEG frequency bands defined for each subject by the estimation of the Individual Alpha Frequency (IAF) value [START_REF] Klimesch | EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis[END_REF]. The PSDs in the theta and alpha bands have then been analyzed by estimating the Coefficient of Determination (r 2 ), or r-square, between the considered experimental condition and the reference condition. As 0 < r 2 < 1 by definition, a signed r 2 has been derived by multiplying the coefficient of determination by the sign of the slope of the corresponding linear model of the regression analysis. In this way, it has been possible to obtain information not only about if the two datasets were different, but also about the direction of such difference. A Stepwise Linear Discriminant Analysis (SWLDA, [START_REF] Aricò | A brain computer interface system for the online evaluation of ATCs' workload[END_REF]) has been used to select the most relevant spectral features to discriminate the mental workload levels. In particular, the classifier was trained using data from one triplet (Easy, Medium and Hard) and the extracted parameters were tested over the other remaining triplets within the same session (INTRA cross-validations) or the other sessions (INTER cross validations). Several moving average samples (NMA) have been applied to the output of the classifiers (WEEG,): NMA(1) = 0.125 (sec), NMA(8) = 1 (sec), NMA(16) = 2 (sec), NMA(32) = 4 (sec), NMA(64) = 8 (sec). The moving average was expected to increase the stability and the accuracy of the index with the drawback of introducing delays in the workload estimation, inducing a decrease of the workload refresh rate.

ECG and EOG analysis:

As well as for the EEG, the ECG and the EOG signals have been band -pass filtered, respectively 1-8 (Hz) and 8-16 (Hz), and then segmented in epochs of 8 seconds, 0.125 secondsoverlapped. The HR and the EBR have been estimated by calculating the distance between consecutive peaks occurring in the ECG and in the EOG signals. In particular the R-peaks and the eyeblinks peaks have been used, and then they have been normalized by the z-score transformation with respect to the reference condition, in which the subjects watched the task interface without responding to them. As for the EEG, also for the HR parameter a workload index (WHR) has been calculated by using the SWLDA at different output rates: NMA(1) = 0.125 (sec), NMA(8) = 1 (sec), NMA(16) = 2 (sec), NMA(32) = 4 (sec), NMA(64) = 8 (sec).) Fusion workload index: A Fusion workload index has been calculated as a combination of the WEEG and the WHR indexes. In particular, the two SWLDA classifiers outputs have been synchronized, because their different delays (EEG: 2 (sec) overlapped of 125 (msec); HR: 8 (sec) overlapped of 125 (msec)), and then a new workload index (Fusion based workload index, WFusion) has been computed as a linear combination of the WEEG and the WHR score (Equation 1).

𝑊 𝐹𝑢𝑠𝑖𝑜𝑛 = 𝑎𝑊 𝐸𝐸𝐺 + 𝑏𝑊 𝐻𝑅 (1)
The coefficients a and b of the linear combination have been estimated for each subject by means of a simple LDA performed considering the EEG and the HR score distributions (WEEG and WHR) calculated over the cross validations for the three different difficulty levels (Figure 2).

Classifier performance analysis:

The dataset deriving from the three days of workload evaluation have been re -organized in 15 triplets (5 triplets per session) of runs (Easy, Medium and Hard). All the possible cross-validations have been considered, training the classifier with one triplet and testing the extracted features over the remaining triplets. The values of the Area Under Curve (AUC) of the Receiver Operating Characteristic (ROC, [START_REF] Bamber | The area above the ordinal dominance graph and the area below the receiver operating characteristic graph[END_REF]), describing the accuracy of the system, have been calculated from the outputs of the classifier (for each of the different output rates).

Workload score distributions analyses: The workload score distributions of the single conditions (Easy, Medium and Hard) has been calculated using the same approach of the AUC evaluation, thus

Statistical analyses:

The results derived from the different analyses have been then validated by the statistical analyses performed by using the STATISTICA software (Statsoft). For the Training Protocol, the one-way repeated measures ANOVA (Confidence Interval, CI = .95) was used for all the neurophysiological data (dependent variables) with the SESSIONS (3 levels) as independent variable. Such factor has three levels, one for each day of the week in which the EEG recording was made (T1, T3 and T5). For the Workload Protocol, statistical analyses over the i) classifier performances, ii) workload scores distribution and iii) NASA-TLX scores have been performed.

i)

Three repeated measures ANOVA (CI = .95) have been performed, one for each classifier (EEG, HR and Fusion based), using the different comparison couples of difficulty levels (Easy vs Hard, Easy vs Medium and Medium vs Hard, 3 levels), and two of the "moving average lengths" (NMA(x), x={1, 64}, 3 levels) as within factors and the related AUC values as dependent variable, for all the subjects. We selected only two moving average values because the number of subjects. Also, a repeated measures ANOVA has been performed by considering the three classifiers (EEG, HR and Fusion based, 3 levels) and the two moving average lengths (NMA(x), x={1, 64}, 3 levels) as within factors, and the AUC values averaged for the three couples of conditions (Easy vs Hard, Easy vs Medium and Medium vs Hard) as dependent variables for all the subjects. In addition, a Duncan post-hoc test has been performed in order to test the effects between all the factors.

ii) Three repeated measures ANOVA (CI = .95) have been performed, one for each classifier (EEG, HR and Fusion based), using the difficulty conditions (Easy, Medium and Hard, 3 levels) and Cross-validation type (INTRA and INTER, 2 levels), for each subject, as independent variables and the related workload index distributions (WEEG, WHR and WFusion) as dependent variables..

iii)

A one-way ANOVA (CI=.95) was performed on the NASA-TLX scores (dependent variable) with the difficulty condition (Easy, Medium and Hard) as independent variable.

RESULTS

Training improvement assessment

LABY performance analysis: Throughout the training sessions, the performance of the subjects increased continuously in terms of mean performance level and accuracy. Figure 3 shows the performance's index adopted across the different training days. By the inspection of Fig. 3 it is easy to note the simultaneous increase of the performances level and the decrease of the amplitude of the standard deviations in the learning curve. On the second day of training, all the subjects reached at a good level of performance (almost the 90%) and since the third day (T3), they could reach performance level higher than 95%. The one-way ANOVA performed on the global LABY score showed significant differences across the sessions (F(4, 20) = 17.74 with a p < .00001, η 2 p=.78). The Duncan post-hoc test showed that the first two sessions (T1 and T2) were statistically different from all the others (p < .0001) while the last three ones (T3, T4 and T5) were not statistically different to each other. Frontal PSD theta: The ANOVA results reported in Figure 4 show a statistical significant modulation of the of EEG PSD in theta band over the frontal areas (EEG channels: AF3, AF4, F3, Fz, and F4) across the different training sessions (F(2, 10) = 4.18, p < .048, η 2 p=.45) and the Duncan's post-hoc test confirmed these differences p < .03. It is evident that in the central session (T3), when the subjects have been supposed to have learnt how to execute correctly the task and focused the cognitive resources for improve their performances, the frontal PSD theta reached the highest increment respect all the other sessions. Heart Rate and Eyeblinks Rate analysis: Figure 6 and7 show the results of the ANOVA analysis of the autonomic parameters of HR and of EBR, respectively. The HR shows that the subjects were emotively engaged in correspondence of the central training session (T3), as the HR in T3 was the highest one, and that at the end of the training period they were more relaxed with the experimental task, as both the HR and the EBR decreased and increased, respectively. In fact, the Duncan's post-hoc tests reported significant (p<.01) differences between the HR and EBR values of the first (T1) and last (T5) training session. In addition, the EBR z-score shows how the subjects kept to pay attention to the task, as it was negative even at the end of the training.

Perceived mental workload: NASA-TLX analysis:

The one-way ANOVA for the NASA-TLX data (Figure 8) shows significant differences among the training sessions (F(4, 20)=7.67 and p< .00065, η 2 p=.61). The post-hoc test allowed to found out that the average scores of the NASA-TLX were statistically different up to the fourth session (T4), whereas the T4 and T5 sessions were perceived similar in terms of perceived workload.

Workload evaluation and classification

Classifier performance analyses: The ANOVA analyses (Figure 9) revealed a significant increment in the classifier' performance at higher moving average values for both the EEG (F(1, 5)=372.14, p=0.00001, η 2 p=.99) and the Fusion based classifiers (F(1, 5)=18.13, p=0.008, η 2 p=.78), while no differences have been highlighted for the HR based classifier (F(1, 5)=1.34, p=.29, η 2 p=.21). In addition, the post-hoc test showed that AUC values calculated using all the three classifiers in the "Medium vs Hard" couple were significantly lower (all p<.001) than the other two ones. Finally, the ANOVAs highlighted a significant interaction between the Classifiers and the Moving averages (F(2, 10)=25.71, p=0.0001, η 2 p=.84). In particular, post-hoc test showed that the fusion based classifier (WFusion) showed higher AUC values than the EEG based classifier (WEEG) at short refresh times (0.125s, p<.001), and higher AUC values than both the EEG and the HR classifiers at long refresh times (8s, p<.05). (Figure 9).

EEG-based HR-based Fusion

-based 0,125 1 2 4 8 Refresh time (seconds) 0,5 0,6 0,7 0,8 0,9 1,0
Area Under Curve (AUC)

.

NASA-TLX analyses:

The ANOVA results revealed a main effect of the difficulty levels (F(1,10)=8.38, p=.016, η 2 p=.46). Subjective perception of the workload increased as the difficulty of the task increased. This result is consistent with the score distribution analyses (Figure 10), showing a high reliability of the estimated mental workload index.

DISCUSSION

The neurophysiological signals, the task performance scores and the experienced workload describe a story in which the subjects found their own strategies for the correct execution of the proposed task (LABY) and then got completely confident with the execution of it. At the central part of the training period (T3) the cognitive and emotive engagement reached the highest levels, as the frontal PSD theta and the HR showed the highest increment. The trends of the parietal PSD alpha and of the EBR showed how the subjects kept to pay attention to the execution of the task. In fact, the parietal PSD alpha and the EBR continued to decrease and to increase, respectively, up to the last training session of the first week (T5). From a perception point of view, the NASA-TLX scores demonstrated that the subjects experienced less workload throughout the sessions, especially at the end of the training period (T5) respect to the beginning of it (T1).

Once the subjects became trained with the LABY task, an algorithm able to estimate in real-time the mental workload of the user has been tested on the next two weeks, by using the combination of EEG and ECG signals. There are several studies in literature in which these biosignals are used for assessing the mental workload, but normally they are used separately. In this study, it has been demonstrated that the combination of EEG and HR allows to differentiate significantly the workload level over three different difficulty levels, showing a high discrimination accuracy (AUC > .7). Furthermore, the fusion of the EEG and HR information allows to significantly increase the reliability of the algorithm with respect of using only EEG or HR alone. Also, the subjective features used for the evaluation of the mental workload remained stable over a week, between the "week 2" and "week 3" of the experimental protocol., so, even after a week, it might not be necessary to recalibrate the algorithm with new EEG data. The aspects related to the classifier stability and accuracy are highly important for the usability of the system. In fact, to use such system in real environments, it could be enough to calibrate the algorithm with the specific parameters of the operator only once and then just use them without further adjustments maintaining a high reliability and stability over, at least, a period of one week. Finally, the calculated workload index showed the same trend of the NASA-TLX workload assessment.

CONCLUSIONS

Two protocols have been presented in this study, the training and the workload evaluation of ATCOs by means of neurophysiological signals. The integration of information derived by the brain activity, through the EEG, and the physiological signals of ECG and of EOG with the supervision of Experts could be used as possible innovative "cognitive metric" for evaluating the degree of the learning and the training progress throughout their periods of formation. Also, this method could be applied when the comparison between subjects is required in terms of required cognitive resources for the execution of specific tasks. In fact, after a fixed period of training it could be possible i) to quantify how well the subjects can complete a task, in terms of cognitive resources necessary to the correct execution, and ii) to compare subject's cognitive performances by estimating the neuro-physiological EEG, HR and EBR parameters presented in this study. In addition, an algorithm able to estimate the mental workload of an operator by using the combination of EEG rhythms and HR signals has been proposed. It has been demonstrated that i) the system is able to significantly differentiate three workload levels related to the three difficulty level of the task employed; ii) the subjective features used for the evaluation of the mental workload remain stable over a week; iii) the combination between the information derived from the EEG and the HR signals allows to significantly increasing the reliability of the system. Finally, iv) the subjective evaluation of the workload shows the same trend of the physiological workload indexes (WEEG, WHR, WFusion). The evaluation of the mental workload by using the information derived by biosignals, allows to have an objective and more reliable measure than using the subjective questionnaires, such as the NASA-TLX. Also, another advantage with respect to the subjective measurements (e.g. NASA-TLX) is the assessment of the real-time variations of the workload within the same task, e.g. each 125 msec. Further experiments will be performed to even further test and extend the long term use of the algorithm. The present study, carried out in a laboratory environment, should be replicated on a larger sample size (more than six subjects) and in a more realistic scenario involving professional operators.

Figure 1 .

 1 Figure 1.

Figure 2 .

 2 Figure 2.

Figure 3 .

 3 Figure 3.

Figure 4 .

 4 Figure 4. Parietal PSD alpha: Figure. 5 shows the trend of the parietal EEG PSD in alpha band over the EEG channels P3, Pz and P4, represented as variation of signed r-square. Repeated measures ANOVA showed significant differences of the parietal PSD alpha (F(2, 10)=9.95 with an associated p = .0042, η 2 p=.67) and a decreasing trend of the spectral PSD from T1 to T5 across the training sessions.

Figure 5 .

 5 Figure 5.

Figure 6 .

 6 Figure 6.

Figure 7 .

 7 Figure 7.

Figure 8 .

 8 Figure 8.

Figure 9 .

 9 Figure 9. Workload score distribution: The ANOVA analyses (Figure 10) revealed that the score distributions related to the different difficulty conditions (Easy, Medium and Hard) for all the three classifiers were significantly separated (EEG-based: F(2,10)=8.90, p=.006, η 2 p=.64; HR-based: F(2,10)=4.93, p=.032, η 2 p=.50, Fusion-based: F(2,10)=10.05, p=.004, η 2 p=.67). Furthermore, no significant differences were found between the workload scores related to the INTER and the INTRA crossvalidations (EEG-based: F(1,5)=.22, p=.66, η 2 p=.04; HR-based: F(1,5)=1.40, p=.29, η 2 p=.22, Fusionbased: F(1,5)=1.78, p=.24, η 2 p=.26).

Figure 1 .

 1 Figure 1. The LABY is a dynamic environment whereby an ATC must issue directional commands

Figure 2 .

 2 Figure 2. Fusion based workload index assessment (WFusion). The Fusion workload index (WFusion)

Figure 3 .

 3 Figure 3. The trend of the global LABY score across the five different training sessions (T1-T5). The

Figure 4 .

 4 Figure 4. Mean EEG PSD (r-square) in theta band over the frontal EEG channels AF3, AF4, F3, Fz

Figure 5 .

 5 Figure 5. Parietal EEG PSD in alpha frequency band during the training period. The graph reports the

Figure 6 .

 6 Figure 6. Heart Rate (z-score) values across the training sessions. The trend shows how in the central

Figure 7 .

 7 Figure 7. Eyesblink rate (z-score) values across the training sessions. The values are all negative

Figure 8 .

 8 Figure 8. Average NASA -TLX scores of the training sessions. After each training session the

Figure 9 .

 9 Figure 9. Mean values and related standard errors (CI = .95) of the AUC values achieved using the

Figure 10 .

 10 Figure 10. Mean values and related standard errors (CI = .95) of the distributions of the workload

This work is co-financed by EUROCONTROL on behalf of the SESAR Joint Undertaking in the context of SESAR Work Package E -NINA research project. This work is also supported by BrainSigns srl, with a grant from Regione Lazio, FILAS, named "BrainTrained". CUP F87I12002500007.


Giovanni Neglia
email: giovanni.neglia@inria.fr

Giuseppe Di Bella
email: giuseppe.dibella@telecomitalia.it.

Laura Giarré

Ilenia Tinnirello

Scalable and Privacy-Preserving Admission Control for Smart Grids

Keywords: smart grids, direct load control, admission control, privacy, large deviation

Energy demand and production need to be constantly matched in the power grid. The traditional paradigm to continuously adapt the production to the demand is challenged by the increasing penetration of more variable and less predictable energy sources, like solar photovoltaics and wind power. An alternative approach is the so called direct control of some inherently flexible electric loads to shape the demand. Direct control of deferrable loads presents analogies with flow admission control in telecommunication networks: a request for network resources (bandwidth or energy) can be delayed on the basis of the current network status in order to guarantee some performance metrics. In this report we go beyond such an analogy, showing that usual teletraffic tools can be effectively used to control energy loads. In particular we propose a family of control schemes which can be easily tuned to achieve the desired trade-off among resources' usage, control overhead and privacy leakage.

Contrôle d'admission passant à l'échelle et préservant la vie privée pour les réseaux intelligents

Résumé : Demande et production d'énergie doivent être constamment équilibrées dans le réseau électrique. Le paradigme traditionnel qui consiste à adapter dynamiquement la production à la demande est contesté par la montée en puissance de sources d'énergie plus variables et moins prévisibles, comme l'énergie solaire photovoltaïque et l'énergie éolienne. Une approche alternative est le contrôle direct de certaines charges électriques intrinsèquement flexibles dans le but de modeler la demande. Le contrôle direct des charges reportables présente des analogies avec le contrôle d'admission des flux de données dans les réseaux de télécommunication : une demande pour les ressources du réseau (bande passante ou énergie) peut être différée selon l'état actuel du réseau afin de garantir certaines mesures de performance. Dans ce rapport, nous allons au-delà d'une telle analogie et montrons que les outils de télétrafic habituels peuvent être efficacement utilisés pour contrôler les charges électriques. En particulier, nous proposons une famille de schémas de contrôle qui peuvent être facilement configurés pour atteindre le compromis souhaité entre utilisation de ressources, coûts du contrôle et révélation d'informations privées.

Mots-clés : Réseaux intelligents, contrôle direct des charges électriques, contrôle d'admission, grandes déviations

Introduction

Direct Load Control (DLC) allows energy utilities to control electric loads at the customers' premises. In the past DLC was used in critical situations to prevent blackouts by shutting down these loads. More recently, an extensive use of DLC has been advocated as a way to shape energy demand peaks or provide other ancillary services [START_REF] Mathieu | State estimation and control of electric loads to manage real-time energy imbalance[END_REF][START_REF] Hao | Ancillary service to the grid through control of fans in commercial building hvac systems[END_REF]3].

In this report we consider a scenario where DLC functionalities are deployed at a large set of small deferrable energy loads, like appliances at residential users. As we are going to discuss, the required additional intelligence and communication capabilities may be introduced through smart plugs, without the need to replace older "dumb" appliances. For such a scenario, we propose a simple control mechanism that guarantees that the instantaneous power demand exceeds a given bound with probability smaller than . Only a stochastic characterization of the power demand from each class of appliances is required. This mechanism combines two different operation paradigms. In the first one, appliances need to ask a controller the permission to start, and the controller will limit the number of simultaneously active appliances to n(t). In the second one, an activation probability function-p(t)-is broadcast periodically to all the appliances; appliances do not notify the controller but they start with probability p(t) and postpone their decision to the time t + T with probability 1 -p(t). The first operation paradigm requires more communication exchanges between the appliances and the controller and reveals more information about the customers' habits. The second paradigm works in an open-loop fashion on the basis of historical aggregated data and, then, does not disclose any private information. At the same time the lack of an exact knowledge of the current number of active appliances generates a lower average utilization of the resources in order to satisfy the constraint. We simply combine the two paradigms by means of a probability q: when an appliance wants to start operating, it will ask the permission to the controller with probability q, and it will decide autonomously using the function p(t) with probability 1 -q. The parameter q can then be chosen in order to achieve the wished trade-off among resources' usage, control overhead and privacy leakage. In fact, by increasing q we gradually i) reduce privacy by exposing more the energy profile of each user, ii) increase control overhead because the controller needs to directly interact with a large number of appliances, iii) increase the efficiency. Fig. 1 qualitatively depicts the effect of the parameter q. The detailed description of our system is in Sec. 3, and our analysis (Sec. 4) and experiments (Sec. 6) allow to quantify these trade-offs.

In our mechanism, the control policy is determined by the two functions p(t) and n(t). In Sec. 4, we show how techniques developed for flow admission control in packet networks can be used to determine a stationary control policy (i.e. p(t) = p and n(t) = n), when the appliances' activation rate is assumed to be time-invariant. Our report contributes then to show how teletraffic engineering tools can be advantageously used also in the context of future smart power grids. In Sec. 5, we derive time-variant control policies for the more realistic case when the appliance activation rate i) is time-varying and ii) needs to be estimated. Our main goal in the present report is not to design a full-fledged solution for DLC, but rather to show how teletraffic tools (such as those used for Call Admission Control in telecommunication networks) can be effectively used to achieve a large scale deployment of the DLC approach to shape energy peaks.

Related Work

Direct load control is a specific mechanism of demand side management, that allows electric utilities to turn specific users' appliances on and off during peak demand periods and critical events. Most of the current large scale DLC programs work on thermostatic loads [START_REF] Mathieu | State estimation and control of electric loads to manage real-time energy imbalance[END_REF][START_REF] Hao | Ancillary service to the grid through control of fans in commercial building hvac systems[END_REF], such as air conditioners and heating systems, because they allow a fine-tuning regulation of power demand. Alternative approaches consider electric vehicles or other battery-empowered appliances [START_REF] Kempton | Electric vehicles as a new power source for electric utilities[END_REF] which, beside acting as adaptive loads, can also reinject energy in the grid. as adaptive loads to be used for regulation. Multiple load typologies, including interruptile or deferrable loads such as the pool pumps, have been also considered for responding to different frequency components of the regulation signal [START_REF] Meyn | Ancillary service to the grid from deferrable loads: The case for intelligent pool pumps in florida[END_REF]. Our focus in this report is on deferrable loads and simple control policies requiring limited intelligence in the network. Direct load control programs have been activated in pilot projects for both residential and industrial customers. For residential customers, some utilities also offer commercial services [3] based on the deployment of radio-controlled switches on air-conditioning units or electric water heaters at the customers' premises. For industrial customers, the control mechanisms are very simple and often based on voice dispatch (i.e. a telephone network) [START_REF] Callaway | Achieving controllability of electric loads[END_REF]. Alternative communication infrastructures are currently considered for future large scale deployments [7]. While power grid operators envision the utilization of advanced metering infrastructure, mainly based on power line communications (PLC) and GPRS technologies, for transporting DLC signals, Internet connections available in most households may enable more flexible programs managed by energy suppliers or third parties [START_REF] Ancillotti | The role of communication systems in smart grids: Architectures, technical solutions and research challenges[END_REF]. Other solutions are considering the use of PLC from energy meters to zone data concentrators, which can access the Internet by means of a public telecommunication network [START_REF]Demand-response technology roadmap[END_REF]. Our solution can benefit from data concentrators' broadcast functionalities to control a large number of appliances, but it is otherwise essentially architecture-agnostic.

The specific DLC control policy can be designed on the basis of dynamic programming optimization [START_REF] Hsu | Dispatch of direct load control using dynamic programming[END_REF], fuzzy logic-based decisions [START_REF] Bhattacharyya | A fuzzy logic based approach to direct load control[END_REF], or other profit maximization schemes [START_REF] Sheble | Direct load control-a profit-based load management using linear programming[END_REF]. In this report we consider simpler approaches rooted in teletraffic engineering which intrinsically work better for large scale systems. Indeed, loads on the electrical grid are multiplexed at different aggregation levels (distribution transformers, primary station, prize zone), similarly to traffic from data sources multiplexed at a router of a different hierarchical level. The use of teletraffic tools is not completely novel. Network calculus has been recently exploited to size energy batteries [START_REF] Ghiassi-Farrokhfal | Toward a realistic performance analysis of storage systems in smart grids[END_REF] and transformers [START_REF] Ardakanian | On the use of teletraffic theory in power distribution systems[END_REF]. Queuing theory is also used in [START_REF] Di Bella | Modeling Energy Demand Aggregator for Residential Users[END_REF] for sizing the population of customers subscribing a DLC program under a given maximum tolerable delay for activating the controlled appliances. In this report we consider a different application of teletraffic tools to smart grids, highlighting the analogy between DLC and call admission control.

Inria

System Description

We consider the problem of peak shaving for an energy supplier, which wants to reduce its customers' consumption during the time of the day when energy costs are higher. To achieve these goals, the energy suppliers may interact with the distribution system operators (DSO), and/or with novel intermediate figures called load aggregators [START_REF] Ruiz | A Direct Load Control Model for Virtual Power Plant Management[END_REF] by means of a communication infrastructure transporting the control messages. Both possibilities are illustrated in Fig. 2 where dashed lines indicate control messages flows. The solution based on the DSO can take advantage of broadcast functionalities deployed at the data concentrators. In fact broadcasting of low-rate control signals can be very efficiently performed by PLCs in the low-voltage (LV) distribution grid [START_REF]An overview to new plc applications and new challenges -maxium integrated products[END_REF]. The aggregator can exploit standard Internet connectivity to interact with appliances at customers' premises. The two approaches can coexist. Our solution can operate in both the scenarios, so we will talk generically about the load controller to denote the entity which drives the appliances. The energy supplier specifies an high-level command for the controller in terms of maximum tolerable probability to exceed a given power demand in the controlled area. The load aggregator is responsible to meet this requirement by deciding about the activation of deferrable loads, such as dish-washers and laundry-machines, in the controlled area. On the basis of the demand prediction, the requirements are mapped into a control signal that is broadcast to all the controlled households for deciding about the admission or the deferral of load activation requests. Each household is equipped with a gateway able to receive the control signal from the controller and to interact with the domestic appliances by means of local area technologies (such as ZigBee or PLC). Smart appliances can be natively equipped with a programmable interface able to communicate with the gateway, while dumb appliances can be controlled by means of smart plugs [START_REF] Ganu | nplug: A smart plug for alleviating peak loads[END_REF].

In our control system, the signal broadcast by the load controller represents the probability p(t) to accept or defer a load activation request during the intervals of the day. At each household, the gateway can autonomously process a novel activation request on the basis of this function: the activation request is accepted with probability p(t) and deferred of a fixed time interval T with probability 1 -p(t). This operation applies to a fraction 1 -q of all the requests. The others are directly forwarded to the load controller and the gateway waits for an ACK from the load controller in order to accept the request. The controller maintains a cap n(t) on the number of appliances active at a given time instant (among those it is aware of, i.e. those whose activation request has been sent to the controller). The ACK can then be delayed until some previously activated appliance does not stop operating. This mechanism allows the load controller to actually estimate the time-varying power demand, but also to have a tighter control on the aggregate power demand and then to achieve a more effective resource utilization in the controlled area.

Fig. 3 summarizes the actuator model at the household: different control operations can be programmed by the load controller by specifying q and p(t). For example, when q is set to 1, all the activation requests are forwarded to the load controller by means of unicast transmissions. This implies a better control on the aggregated power demand, because the number of active appliances is known and the only source of randomness is due to the appliance consumptions. On the contrary, when q is very small, most of communications are unidirectional (from the load controller to the households), while decisions can be taken locally with minimal delays. However, the aggregate power consumption is affected by two sources of randomness (the one related to the number of appliance active, and the other one related to the appliance consumption) which require a lower admission rate for not exceeding the power constraint.

Control in a Stationary Setting

Our goal is to design a direct load control system able to impose a certain threshold on the overall power absorbed from a set of appliances with known statistical properties. In the following, we will consider each appliance power consumption as a stochastic process, characterized by a random variable X with probability density function equal to f X (x).

Let m(t) be the number of active appliances at time t, then the total absorbed power is

P (t) = m(t) i1
X i . First we consider that the energy utility wants to guarantee that P (t) does Inria not exceed the threshold P with probability larger than , i.e.:

Pr   m(t) i=1 X i > P   < . (1)
Then, in Sec. 4.4, we show how similar results can be derived for a different type of constraints on P (t).

Classic large deviation results for q = 1

If q = 1, i.e. if all the queries are forwarded to the controller, our problem is to determine the maximum value n, such that

Pr n i=1 X i > P < . (2)
The same problem has been considered in the context of Call Admission Control (CAC) in telecommunication networks. The interpretation of the quantities is different: the purpose is to determine the maximum number n of homogeneous data flows, each with instantaneous rate X i in order to guarantee that the aggregate rate on a link exceeds the value P (e.g. the link capacity) with probability at most . Different approaches have been proposed to determine n. Here, we introduce a simple one based on basic large deviation results (see e.g. [START_REF] Kelly | Stochastic Networks[END_REF]Ch. 6]), but more sophisticated ones could also be applied [START_REF] Chang | Performance Guarantees in Communication Networks, ser. Performance Guarantees in Communication Networks[END_REF]. Our main goal is hereafter to show how teletraffic tools can in general be reused in this different context, rather than develop the solution in all its details.

Let M X (s) = ln E[e sX ] be the cumulant generating function. From the Chernoff bound and basic properties of the cumulant generating function, it follows easily that:

ln Pr(X 1 + ... + X n > nc) ≤ n inf s≥0 [M X (s) -sc] (3)
Moreover, Cramer's theorem proves that this bound is asymptotically tight:

lim n→∞ 1 n ln Pr(X 1 + ... + X n > nc) = inf s≥0 [M X (s) -sc] (4)
Practically speaking, when n is large, Cramer's theorem can be used to approximate the probability that the sum of n independent random variables exceeds a bound ( P in our case) as follows:

ln Pr(X 1 + ... + X n > P ) ≈ inf s≥0 [nM X (s) -s P ] ,
and then n can be determined as the largest integer such that

inf s≥0 [nM X (s) -s P ] < ln .
The set of values {0, 1, . . . , n} is called the acceptance region for the admission controller.

The approach can be easily generalized to a finite number of appliances' classes and the concept of effective bandwidth can be defined for each class in order to easily derive a subset of the acceptance region [START_REF] Kelly | Stochastic Networks[END_REF]Ch. 6].

Extension to q < 1

In this section we show how the results above can be extended to the scenario where a subset of the appliances do not perform any query to the controller, but start autonomously their operation. For the moment we assume that such appliances turn on according to a homogeneous Poisson process with rate λ c . In the following Sec. 5 we explain how this rate can be determined taking into account the retrial mechanism described in Sec. 3 and how one can deal with the actual time-variant activation process. For the moment, we just anticipate that λ c is a linear function of 1 -q and of the instantaneous probability of activation p: λ c ∝ (1 -q)p. In order to keep the equations simple, we assume that all the appliances have the same activation time equal to D, but the results can be easily extended to the case when D is a random variable.

Let N be the random number of appliances starting autonomously. Under the assumptions indicated above, N is distributed as a Poisson random variable with parameter E

[N ] = Λ c = λ c D. The condition in Eq. (2), then becomes Pr   n i=1 X i + n+N j=n+1 X j > P   < , (5)
where (X i ) i∈N is a sequence of independent and identically distributed random variables. Large deviation results hold for large systems, for example when both the number of random variables (n) and the threshold to be exceeded (nc) diverge as in Eq. ( 4). In Eq. ( 5) we need to let both the first and the second addend scale in the same way, otherwise one of them would become negligible in comparison to the other. We assume then that the mean Λ c of the number of appliances starting autonomously scales linearly with n:

Λ c = nΛ 0 c = n(1 -q)pDλ 0 c
, where we also took into account λ c 's dependence on p and 1 -q. It is then possible to prove the following result: Theorem 1. Let (X i ) i∈N be a sequence of independent and identically distributed random variables with cumulant generating function M X (s), and N be a Poisson random variable with mean nΛ 0 c . We assume that M X (s) is defined for any s and Pr{X > c} > 0. It holds:

lim n→∞ 1 n ln Pr   n i=1 X i + n+N j=n+1 X j > nc   = inf s≥0 M X (s) + Λ 0 c e M X (s) -1 -sc .
Proof. Let us define {Y i,k , i = 1, . . . n, k ∈ N} to be a set of independent random variables distributed as X and {N i , i = 1, . . . n} to be a set of independent Poisson random variables with mean Λ 0 c . The sum of the variables N i has the same distribution of N and then:

n i=1 Ni k=1 Y i,k d = n+N j=n+1 X j . Let Z i X i + Ni k=1 Y i,k , for i = 1, . . . n. It follows that n i=1 Z i d = n i=1 X i + n+N j=n+1 X j .
All the variables Z i are independent and identically distributed with cumulant generating function M Z (s) = M X (s)+Λ 0 c e M X (s) -1 . The thesis follows by applying directly Cramer's theorem to In a finite-size system, this result is used to approximate the probability that the total power consumption exceeds the bound P as follows:

ln Pr   n i=1 X i + n+N j=n+1 X j > P   ≈ inf s≥0 nM X (s) -λ c D e M X (s) -1 -s P .
and then the acceptance region is made by all the pairs of values (p, n) such that:

inf s≥0 nM X (s) + λ c D e M X (s) -1 -s P ≤ ln . (6)
In the next subsection, we characterize such region for a simple case.

A numerical example

In order to illustrate the admission control rule determined by the results above, we consider a toy example where there is a single class of appliances, which would like to activate according to a Poisson process with rate λ. Then, there will be a Poisson process with rate qλ to the query-response system, and a spontaneous activation process of appliances with rate p(1 -q)λ (we ignore retrials for the moment). We also assume that the appliance's consumption profile can be modeled as the two-state Markov process in Fig. 4. The appliance stays in the high consumption state on average 10 minutes and consumes 1.5kW, it stays in the low consumption state on average 20 minutes and consumes 500W. Literature suggests that Markov processes can indeed accurately model appliances' consumption [START_REF] Widén | A high-resolution stochastic model of domestic activity patterns and electricity demand[END_REF]. For this appliance, the cumulant generating function can be easily calculated as M X (s) = ln 1/3e s1.5 + 2/3e s0.5 . Then inequality (6) can be used to determine the acceptable pairs (p, n) with λ c = p(1 -q)λ:

inf s≥0 nM X (s) + p(1 -q)λD e M X (s) -1 -s P ≤ ln . (7)
Fig. 5 shows the frontiers of 10 different acceptance regions calculated for λ = 12, D = 90 minutes, q = 0.1, 0.2, . . . , 1, = 0.1 and P = 0.8λD E[X] (the 80% of the average absorbed power in absence of the control). Each curve2 is made by the pairs (p, n) for which the constraint in ( 6) is satisfied with an equality. The acceptance region is then made by all the points that are below this curve. This is evident analytically because the left term in (7) is decreasing in n and p, but it is also intuitively clear because of the physics of the system: for a given acceptable pair (p, n), if we limit more p or n (p < p or n < n), then less appliances will be admitted in the system and the constraint will be even more satisfied (both (p , n) and (p, n ) are acceptable). We observe the expected trade-off between n and p: on the frontier as p increases, n needs to decrease. The more appliances we decide to admit through the query-response system, the more we have to limit through p the expected number of appliances which will activate autonomously.

Another interesting remark about the frontiers in Fig. 5 is that they look almost linear. This can clearly significantly simplify the determination of the acceptance region and leads to important conclusions about the efficiency of resources' usage. For this reason, we investigate deeper such behavior. The term e M X (s) grows very fast for s ≥ 0 and because of this the point of minimum for the left-hand side of Eq. ( 7) is close to 0, specially for larger p(1 -q). It is then possible to consider a first order Taylor approximation (e M X (s) -1) ≈ M X (s) and the frontier can be approximated as follows:

inf s≥0 (n + (1 -q)pDλ) M X (s) -s P = ln ,
For a given value of q, the frontier is made by pairs (p, n) such that:

n + (1 -q)pDλ = constant. ( 8
)
To completely characterize this linear relation, it is enough to observe that for p = 0 the maximum value for n (let us denote it as n * p=0 ) does not depend on q, i.e. all the frontiers pass by the point (0, n * p=0 ) and they are described by the following linear equation parameterized in q:

n + (1 -q)pDλ = n * p=0 . (9)
If there are n appliances activated through the query-response mechanism, the expected number of appliances in the system is n + (1 -q)pDλ and then the expected power consumption is E

[P ] = (n + (1 -q)pDλ)E [X]
. This means that for a given q value all the points in the frontier described by Eq. ( 9) have the same resources' usage E [P ] = n * p=0 E [X] while satisfying the constraint. Moreover, for any value of q we can achieve the same efficiency by selecting (n, p) on the corresponding frontier. Although Eq. ( 9) stems from a linear approximation, it suggests that the efficiency (in terms of average utilization) of our control mechanism is not very sensitive neither to the pair (n, p) (as long as it is a point of the frontier), nor to the parameter q. This is confirmed by our numerical analysis: for example we observed less than 5% reduction of the maximum expected load changing q from q = 1 to q = 0 despite the fact that for q = 0 the control needs to deal also with the variability of the number of appliances admitted. However, the importance of direct queries in practical applications is not only due to efficiency reasons, but also to the possibility to accurately estimate the appliance activation rate as discussed in Sec. 5. Appendix B further discusses the relation between n, p, q and average utilization.

A different constraint

Before describing how the above control policy can actually be used in a case when the activation pattern is time-varying and unknown, we want to briefly show that teletraffic tools can also be used to deal with different constraints. In particular, the constraint Pr(P (t) > P ) ≤ ignores the time duration during which the constraint P is violated. It may be significant for the energy utility to target the following different type of constraint.

Let us define an energy bucket as follows: the bucket can store an energy capacity equal to B, it is filled with constant rate P and emptied with the variable rate P (t). Let L(t) denote the energy stored in the bucket at time t. Fig. 6 describes on the left side the energy-buffer. While inspired by the token-bucket used to shape data traffic, this energy-bucket is only used to express a constraint. If one requires that L(t) > 0, ∀t then P (t) can occasionally exceed the constraint P by δP , but only for the time interval [t, t + L(t)/δP ]. The ratio behind such bound is that smaller violations would be tolerated for a longer time. Given the incertitude on the energy demand, the control policy cannot guarantee the deterministic bound described above, but rather its probabilistic version:

Pr(L(t) = 0) < , ∀t. (10)
The possibility to use teletraffic tools (and in particular large deviation results) also to deal with the constraint in [START_REF] Hsu | Dispatch of direct load control using dynamic programming[END_REF] derives from the following equivalence: the energy-bucket on the left of Fig. 6 is equivalent to the queue in the same figure on the right, where the input and the output process have been reversed, i.e. a traffic with rate P (t) is offered to a queue with service rate P . Specifically, if L(0) = B and Q(0) = 0, then L(t) = Q(t) for every t ≥ 0. This result is quite evident, but we have not found it in the literature, so we prove it in Appendix C. The important consequence for our purposes is that controlling the activation of electric appliances so that the aggregate consumption rate satisfies constraint ( 10) is equivalent to admit data flows to a queue so that the probability to lose packets because of buffer overflow is smaller than . In fact (under large deviation approximation):

Pr(Q(t) > B) ≈ Pr(Q(t) = B) = Pr(L(t) = 0).

Similar results for the acceptance region hold [19, Ch. 6]:

Pr(L(t) = 0) ≈ exp sup t≥0 inf s≥0 stnα(s, t) -s(B + P t) , where α(s, t) = 1

st log E e sX(0,t) , with X(0, t) being all the random energy consumption of one appliance during [0, t].

This example has further shown the potential relevance of applying teletraffic tools to load control in smart grids.

Control in a Time-Variant Setting

In this section, we show how the analytical results derived in the previous section can be practically used when the appliances' activation process is not stationary, its rate is unknown, and probabilistically controlled appliances retry to activate some time later.

Indeed, the usage of electric appliances exhibits a strong time-of-the-day effect. For example Fig. 7 shows the activation rate of washing machines over 15 minutes intervals, as derived from data in [START_REF] Miceli | Sustainable development and energy saving laboratory[END_REF]. We need then to take into account the effect of a time variant λ(t), which is generally unknown even if historical data may be available. Moreover, because of the retrial mechanism from probabilistically controlled appliances, the actual activation rate results time-varying even when the spontaneous activation rate λ is constant.

Inria

Because of the time variability of the process to control, the control policy needs to be timevarying too, i.e. in general we will have n(t) and p(t). There are different possible choices on how to jointly adapt the two control actions, which lead to different performance in terms of resources' usage (see Appendix B for the effect of n and p on the efficiency), communication requirements (e.g. if p is constant, it does not need to be transmitted periodically) and fairness between the two groups of appliances-those controlled probabilistically and those controlled through query-response-(e.g. in terms of delay before the activation). For simplicity, in this paper, we decided to consider n(t) = const = n constant and then to compensate for process changes by dynamically tuning p(t). The control is time-slotted with time intervals of length T c . Without loss of generality we consider that the control starts at time t = 0, and we denote by p k the value of the control action during the k-th time slot, i.e. p(t) = p k for t ∈ [(k -1)T c , kT c ).

For the sake of simplicity we will also assume that the retrial delay T is equal to T c and the activation time of the appliance is a multiple of T c (D = dT c ), even if the three parameters are in general independent.

Fig. 8 shows the model of the whole system we are going to describe in the following sections.

The plant

As it is usual in control theory, we call plant the combination of the process under control and the control actuator. The input is the spontaneous activation process with rate λ k during the k-th time slot, i.e. the process of the activation instants in absence of any form of control (we omit time dependence in the figure). The control system assumes λ(t) to be constant during a control slot, while this is not necessarily the case. In Sec. 6 we evaluate the effects of such an approximation. We keep assuming that the point process of all the activation time instants can be correctly modeled by a (non-homogeneous) Poisson process with rate λ(t), because it originates from the superposition of many independent individual choices (users deciding to turn on their appliances). The initial requests are randomly split in two independent Poisson processes with rate respectively qλ(t) and (1 -q)λ(t). The rate of appliances that will consider to activate autonomously in slot k is λ eq,k . It holds:

λ eq,k = (1 -q)λ k + λ eq,k-1 (1 -p k-1 ), (11)
where the first term is due to all the appliances that are considering to activate themselves for the first time during the k-th time slot, and the second is due to those that have already considered this decision in the (k -1)-th slot and have postponed it to the current one. The actual activation rate of probabilistically controlled appliances3 is

λ c,k = p k λ eq,k .
Finally, the actual number of active appliances m pc,k admitted through the probabilistic control and active at time kT c is equal to those arrived in the interval [kT c -D, kT c ]. This is a Poisson random variable with expected value

E [m pc,k ] = k h=k-d+1 λ c,h .

The controller

The number of appliances active at a given time does not exceed n. Let m qr,t denote the number of appliances active at time t.

Estimators. The controller does not know the state of the plant (e.g. how many appliances are taking the decision to activate autonomously). It needs then to estimate the rates λ and λ eq . An estimate for quantity x is denoted as x. In particular two different set of estimates will be useful: λp k and λp eq,k will estimate the sequences until the current slot k, while λf k+1 and λf eq,k+1

will be used as a prediction for the slot k + 1 in order to determine p k+1 . In this paper we consider simple estimators for these quantities, but we show that they work reasonably well in Sec. 6. Given N r,k the number of queries received during the k-th slot, the Maximum Likelihood Estimator for λ k is simply λp k = N r,k /(qT c ). Clearly, this estimate could be improved if some a-priori is available (e.g. from historical data). The prediction for slot k + 1 is λf k+1 = λp k . The estimates λp eq,k and λf eq,k+1 are simply obtained as a function respectively of ( λp k , λp eq,k-1 ) and of ( λf k+1 , λp eq,k ) using the same relation that links λ eq,k to λ k , i.e.

λp eq,k = (1 -q) λp k + λp eq,k-1 (1 -p k-1 ).
λf eq,k+1 = (1 -q) λf k+1 + λp eq,k (1 -p k ). Control logic. The controller determines p k+1 on the basis of the acceptance region derived as described in Sec. 4. For simplicity we consider a linearized frontier, but the following operation could be adapted to the case when the frontier is described by a more generic function. Moreover, in Sec. 4.3 we have shown that the linear approximation works very well. Let n (q) (p) be the frontier parametric equation for a given q value. We denote the minimum value of n on the frontier by n (q) min (it is 0 if the frontier intersects the segment [0, 1]). Moreover let p (q) max max{p|n (q) (p) ≥ n (q) min }, i.e. the largest probability value in the frontier. For the sake of fairness Inria between the probabilistically controlled appliances and the querying ones, we would like the system to work in a point with coordinates (p, n) approximatively in the middle of the frontier. Under the linear approximation, and considering the maximum arrival rate in the control period, the corresponding value of n is n (q) min +λ max (1-q)p (q) max /2. It has to be remarked that the expected number of appliances at time t that have been activated through the query-response mechanism can be at most t t-D λ(τ )q dτ (see discussion in Appendix B). Then we have set n as follows n = min 0.9 max

t∈[Tsc,Te] t t-D λ(τ )q dτ , n (q) min + λ max (1 -q) p (q) max 2 . ( 12
)
We then select the control probability p in order to keep the total expected number of appliances constant, i.e.

Λ k+1 + n = const, ( 13
)
where Λ k+1 is the expected number of active probabilistically controlled appliances at the end of the k + 1-th slot, and the constant is equal to n

(q) min + λ max (1 -q)p (q)
max . 4 In the stationary case without retrials, it was simply Λ k+1 = λ c D = λ(1 -q)pD. Here, we can express Λ k+1 as the sum of two terms, one (denoted as Λ hist,k+1 ) due to all the probabilistically controlled appliances already in the system at the begin of slot k + 1 (or equivalently at the end of slot k), the other due to the estimated number of appliances which will activate during slot k + 1. We can estimate Λ hist,k+1 as follows

Λhist,k+1 = d-2 h=0 p k-h λp eq,k-h T c ,
and then Λk+1 = Λhist,k+1 + λf eq,k+1 T c p k+1 .

Finally, p k+1 can be iteratively derived from Eq. ( 13) and Eq. ( 14) as:

p k+1 = n * q -n -Λhist,k+1 λf eq,k+1 T c .
This may lead to a too prudential strategy, because the number m qr,k+1 of active appliances in the query-response queue can be significantly smaller than n, specially at the begin of the control period (the queue fills initially with rate λq) and then the configuration above would lead to a severe underutilization. A solution is to predict the number of active appliances in the query-response queue during the residual control period (i.e. mqr,h for h > k), on the basis of the queue status at time k and expected activation requests. Indeed, the activity state of the enqueued appliances can be deterministically calculated from their activation/de-activation time instants, while the number of future arrivals can be predicted from the current estimate of the arrival rate λf k+1 . In particular let a qr,h be the total number of number of appliances in the query-response system (active or waiting) at slot h. This can be estimated as sum of two terms. The first term is the number of appliances currently (at slot k) in the query-response system which will be in the system at slot h. It can be immediately calculated from the current status of the queue (number of active appliances, the activation time instants, and number of appliances in the queue), we denote it as a p qr,h . The other term (â f qr,h ) requires the prediction of 4 The reader may wonder why we have introduced n (q) min and p (q) max and in particular why the linearized frontier is not expressed simply by Λ k+1 + n = n * p=0 . The answer to these questions is in Appendix B. future arrivals from now h and can be estimated as âf

RR n°8769

qr,h = λ f k (h -k)T c .
The estimate for m qr,h is then mqr,h = min(n, a p qr,h + âf qr,h ). Then, p k+1 can be calculated replacing n in Eq. ( 13) with the most pessimistic forecast until the end of the control period, i.e. with max{ mqr,h , h > k}. This is the approach we adopted.

Numerical Results

In this section we show the performance of the system described in Sec. 5 in a realistic setting. In particular we consider 30000 washing machines under control. The instantaneous power consumption of a washing machine is assumed to follow the simple model in Fig. 4 and the activation time is D = 90minutes. The spontaneous activation rate λ(t) is the one in Fig. 7, derived from experimental data in [START_REF] Miceli | Sustainable development and energy saving laboratory[END_REF]. The largest expected power demand P max is at around time 11.00am. We assume the energy supplier sets the constraint as Pr(P > P = 0.8P max ) < 0.1 in the interval [T s , T e ] = [10.00am, 11.30am]. In order to be sure to satisfy the constraint at t = 10.00am, the control needs to start at T sc = 8.30am.

Fig. 9 plots the evolution of the power demand with and without control for q = 0.5 and T c = 15minutes together with the probability signal p(t). We observe that the controller does not actually affect the system (p(t) = 1) until t = 10.00am and that it actually manages to maintain the absorbed power below P for the whole duration of the control interval. Observe also how power consumption significantly increases after T e . This is due to the fact that a severe constraint has been imposed for a long time interval. Moreover, the power increase can be made smooth by gradually increasing P after T e .

The following figures 10 and 11 show the utilization-evaluated as ratio between the average power demand and P -and the probability to exceed the bound at time t = 10.26am estimated over 1000 simulations for different values of q ∈ [0, 1] and for T c = 1, 5, 15, 30 minutes. The time instant of observation falls in the interval where there is the largest activation rate after a period when the rate has been almost constantly increasing. It is then a particularly critical instant for The upper plots of figures 10 and 11 correspond to the ideal case when the controller has perfect estimation of the average request rate in the next timeslot, i.e. λf k+1 = λ k+1 . We observe that, as we expected (see Appendix B), utilization increases as q increases, but no more than 3%. The probability to exceed the bound roughly reflects the same trends, with higher probabilities corresponding to higher utilization. In any case the probability values are well below = 0.1. The curves almost overlap for all the values of T c but T c = 30minutes. This is due to the fact that the actual activation request rate is constant over 15 minutes time intervals, then for T c = 1, 5, 15 minutes, the knowledge of the average rate in the next slot corresponds to the knowledge of the actual rate. For T c = 15 minutes, there is some mismatch for some estimates, but it does not affect much the operation, while for T c = 30 minutes, the average arrival rate is a bad predictor for the actual arrival rate.

The lower plots in figures 10 and 11 show the same metrics when the simple estimators described in Sec. 5 are used. In this case, we expect an increase of utilization, due to the fact that the controller will usually underestimate λ k+1 for t < 10.00am (because the arrival rate keeps increasing), so it will select a too high probability p k+1 allowing the activation of a number of appliances larger than the correct value. This error has larger consequences for small q, when a larger percentage of appliances is activated through the probabilistic control. For T c different from 30 minutes, there is only a slight increase in the utilization, but still it has a significant effect on the overload probability. In particular for q < 0.1 the constraint is no more satisfied. For T c = 30 minutes, the controller uses the average rate measured in [9.30am, 10.00am] to estimate the arrival rate during the interval [10.00am, 10.30am] with about a 20% of relative error. This justifies the bad performance achieved with this setting. We observe that, even without using better estimators, one could simply counteract the estimation errors, for example by reducing p k+1 by a given factor corresponding to the maximum variability of the arrival rate from a control slot to the following one.

We also show the curves for the time-average utilization (Fig. 12 and Fig. 13), and overload probability (Fig. 14 and Fig. 15) over the whole control interval (for both the cases of exact and estimated λ k+1 ). A final aspect to evaluate is how the communication overhead at the controller depends on the request arrival rate, the control interval [T sc , T e ] and the parameters q and T c . Here we refer to the scenario when the controller is deployed at the DSO and broadcast functionalities are available at data concentrators (as discussed in Sec. 3). Let η be the number of data concentrators. Each query forwarded to the controller during the interval [T sc , T e ] generates an ACK message and during this interval the controller needs to send one message every T c to each of the η concentrators. Then, the expected number of messages at the controller is:

2 Te Tsc qλ(t)dt + T e -T sc T c η (15)
An alternative solution is to rely only on queries. In this case the expected number of queries is 2 Te Tsc λ(t)dt + η, because all the queries will be forwarded to the controller and one message needs to be sent to all the concentrators (and then broadcast to all the appliances) to specify the interval of control.5 Fig. 16 plots the expected number of messages obtained from Eq. ( 15) for different values of the number of concentrators η and T c = 15minutes, as well as for the solution relying only on queries (the horizontal line). Different values of T c would simply shift up or down the costs without changing the slope. The intersection between one of the oblique lines and the horizontal one determines the q value for which the two approaches have the same communication overhead. For smaller values of q our solution is more advantageous in this regard, and obviously in terms of privacy.

Conclusion

In this report we propose a DLC scheme for smart grids that can work with a large number of dumb appliances, and can be configured for providing the desired trade-off between resource utilization, communication overhead and privacy. The control system has been designed by using well-established teletraffic tools for admission control, which work well for large scale systems, in order to provide probabilistic guarantees on the power demand in the controlled zone. The main idea is controlling a class of electric appliances by combining a centralized queryresponse system with a probabilisitic system (periodically programmed by the central controller) able to take local decisions. The highest is the ratio between the activation requests that are processed locally and the ones that are forwarded to the controller, the lowest are the information leakage on user habits and the communication overheads, at the expense of control efficiency (in terms of resource utilization and probability to exceed the required power demand). [7] H. Farhangi, "The path of the smart grid," Power and Energy Magazine, IEEE, vol.

A Alternative Proof of Theorem 1

Let Y n = n i=1 X i and W n = n+N j=n+1 X j where N is a Poisson r.v. with parameter nΛ 0 (and then W n is a compound Poisson r.v.). It holds M Yn (s) = nM X (s) and M Wn (s) = nΛ 0 c e M X (s) -1 . Define also Z n = Y n + W n and let f 1 (s) = M X (s) and f 2 (s) = Λ 0 c (e M X (s) -1)
. We apply Chernoff bound to the variable Z n -nc, we obtain:

ln Pr(Z n -nc > 0) ≤ inf s≥0 ln E e s(Zn-nc) = inf s≥0 [ln E e sYn + ln E e sWn -nsc] = inf s≥0 [nM X (s) + nΛ 0 c e M X (s) -1 -nsc] = n inf s≥0 [f 1 (s) + f 2 (s) -sc] .
M X (s) is defined for any s and then derivable (see [START_REF] Ganesh | Big queues, ser. Lecture notes in mathematics[END_REF]). Given that M Zn (s) = n M X (s) + Λ 0 c e M X (s) -1 , also M Zn (s) is defined and derivable for any

s. M Zn (0) = E [Z n ] = n E [X] (1 + Λ 0 c ) < nc then M Zn (s) -c is decreasing for s = 0. Moreover M Zn (s) ≥ n Pr(X > c)e sc
and then M Zn (s) -nsc diverges when s diverges. It follows that it exists a finite s * > 0 at which M Zn (s) -nsc is minimized (M Zn (s) is also convex, then s * is unique) and M Zn (s * ) -nc = 0, i.e. f 1 (s * ) + f 2 (s * ) = c.

We can then consider c 1 = f 1 (s * ) and c 2 = f 2 (s * ) and then s * is a point of minimum for the functions f 1 (s) -sc 1 and f 2 (s) -sc A similar large deviation result for a compound Poisson process is in [START_REF] Tao | Large deviation estimate in ruin theory[END_REF] and leads to:

B On the Frontiers

If we choose a control setting (p, n) and there are actually n appliances activated in the system through the query-response mechanism, then the expected power consumption is We say that a control setting uses the resources more efficiently (shortly, it is more efficient) than another, if it can sustain a larger expected power demand without violating the bound. The linear approximation for evaluating the frontiers of the acceptance regions given in Eq. ( 9) leads to interesting consequences about efficiency, because it implies that all the pairs (p, n) on the frontiers have the same efficiency independently from q. This result is counter-intuitive: we could expect that the query-response control offers a better handle on the total power consumption, because the exact number of working appliances is known, and only the randomness of individual power absorption determines the variability of power consumption. On the other hand, the probabilistic control has another source of variability, i.e. the number of appliances actually admitted (that is a Poisson random variable). In this section we are going to show that this is the case, and the independency is simply and artefact of the linear approximation. The first observation is that the frontiers are not segments but in any case they are convex curves. Indeed, it can be shown from Eq. (7) that the acceptance region is concave because its complement is the intersection of a family of half-planes and then convex. Moreover, for p = 0 the maximum number of appliances that can be admitted does not depend on q, then all the frontiers intersect the y-axis in the same point (0, n * p=0 ). Fig. 17 shows then a qualitative representation of what the frontiers look like. The loci with constant expected power consumption are characterized by n + p(1 -q)Dλ = const.

E [P ] = (n + p(1 -q)Dλ)E [X] . (16)

Frontiers

Constant load

We compare the slope of such segments (-(1 -q)Dλ) with the slopes of the tangents to the frontiers. Given a pair (p, n) on a frontier, let s * ≥ 0 be a value at which the infimum of nM X (s) + p(1 -q)λD e M X (s) -1 is attained. The tangent to the frontier in (p, n) has slope:

-(1 -q)λD e M X (s * ) -1 M X (s * ) ≤ -(1 -q)λD,

Figure 1 :

 1 Figure 1: Different trade-offs among resources' usage, control overhead and privacy leakage achievable by tuning the parameter q.

Figure 2 :

 2 Figure 2: Reference scenario: communication infrastructure and control system for DLC managed by a DSO and/or a load aggregator.

Figure 3 :

 3 Figure 3: Household control model in terms of activation probability p(t) and ratio of forwarded queries q.

Figure 4 :

 4 Figure 4: Two-state Markov process as a simplified model of the appliance power consumption.

Figure 5 :

 5 Figure 5: Acceptance regions of the pairs (p, n) for different values of the querying probability q.

Figure 6 :

 6 Figure 6: Analogy between the energy bucket on the left and a work-conserving queuing system on the right.

Figure 7 :

 7 Figure 7: Aggregate activation rate for 30000 washing machines.

Figure 8 :

 8 Figure 8: System block diagram.

Figure 9 :

 9 Figure 9: Instantaneous absorbed power (with and without the applied control), and the computed activation probability p(t).

Figure 10 :

 10 Figure 10: Resources utilization at 10.26am.

Figure 11 :

 11 Figure 11: Probability to exceed P at 10.26am.

Figure 12 :

 12 Figure 12: Average resources utilization over the entire control interval [T s , T e ] with known λ.

Figure 13 :

 13 Figure 13: Average resources utilization over the entire control interval [T s , T e ] with the estimation of λ.

Figure 14 :

 14 Figure 14: Average probability to exceed P over the entire control interval [T s , T e ] with known λ.

Figure 15 :

 15 Figure 15: Average probability to exceed P over the entire control interval [T s , T e ] with the estimation of λ.

2 Figure 16 : 1 n

 2161 Figure 16: Expected number of exchanged messages for different numbers of data concentrators.

lim n→∞ 1 n 1 n

 11 ln Pr (W n > nc 2 ) = inf s≥0 Λ 0 c e M X (s) -1 -sc 2 = inf s≥0 [f 2 (s) -sc 2 ] .Both the upper and the lower bound for 1/n ln Pr (Z n > nc) converge to inf s≥0 [f 1 (s)+f 2 (s)sc], then it follows:lim n→∞ ln Pr (Z n > nc) = inf s≥0 M X (s) + Λ 0 c e M X (s) -1 .

Figure 17 :

 17 Figure 17: Qualitative behavior of the frontiers. The segments are the geometric loci with constant expected power consumption.

Inria 5 23 C

 523 Control in a Time-Variant Setting 12 5.1 The plant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 5.2 The controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Equivalence between the Energy-Bucket and the Queue 25 Inria RESEARCH CENTRE SOPHIA ANTIPOLIS -MÉDITERRANÉE 2004 route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

  8, no. 1, pp. 18-28, January 2010.

		0.07
	Average Overload Probability	0.01 0.02 0.03 0.04 0.05 0.06
		0
		0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
						Probability q

RR n°8769

An anonymous reviewer from IEEE CDC

suggested this simpler proof. Our original proof is in Appendix A.Inria

Obviously n can only assume integer values, but here we are considering a continuous relaxation. In practice, the values can then be lower-rounded.RR n°8769

It is possible to show that also the point process of such activations is a Poisson process[START_REF] Neglia | Unidirectional Direct Load Control through Smart Plugs[END_REF].RR n°8769

If the interval is not specified the appliances would need to query the controller all the day along, with an expected number of messages 2 one day λ(t)dt + η.Inria

because for any x it holds e x ≥ x + 1. Taking into account these considerations, we have plotted in Fig. 17 two loci with constant power demand (dashed lines). The figure confirms then our intuition: for a given frontier (a given q value), the larger n, the higher the efficiency. Despite all these considerations, as Fig. 17 shows and Eq. ( 9) justifies analytically, this effect is negligible.

There is another aspect to be mentioned. We have assumed above that by selecting a pair (p, n) we can actually have n appliances activated through the query-response system, but this may not be the case if q is small. In particular if qλD << n, the number of appliances activated through the query-response system will be smaller than n most of the time. Observe in fact that qλD is the expected number of requests to the controller during an interval n. If such number is much smaller than n, the cap will be ineffective most of the time. Practically speaking the maximum value of n to be considered for a given q is qλD. This remark is reflected in our way to set n, see Eq. [START_REF] Sheble | Direct load control-a profit-based load management using linear programming[END_REF].

We conclude this appendix with some additional consequences of Eq. ( 17) and the convex shape of the frontiers on the use of the linear approximation. If we consider a segment with slope -(1 -q)Dλ passing by a point (p (q) max , 0), the segment is guaranteed to lie inside the acceptable region. On the contrary, a segment with the same slope passing by (0, n * p=0 ), would lie outside the acceptable region (even if very close to it). This is the reason why in Sec. 5 we have introduced p (q) max to configure the system instead to rely on Eq. ( 9).

C Equivalence between the Energy-Bucket and the Queue

We show the equivalence for generic time-variant input/output rates. In particular we replace P with T (t).

We write the equations for the two dynamics:

If L(0) = L 0 and B(0) = B -L 0 , it can be easily checked that it will be dQ(t) dt = -dL(t) dt for any t ≥ 0. Then if L(t) is the solution of [START_REF] Ganu | nplug: A smart plug for alleviating peak loads[END_REF] with