Interfacial Instability during Granular Erosion

Gautier Lefebvre, Aymeric Merceron, and Pierre Jop * Surface du Verre et Interfaces, UMR 125 CNRS/Saint-Gobain, 39, quai Lucien Lefranc, F-93303 Aubervilliers, Cedex, France (Dated: September 24, 2015)

The complex interplay between the topography and the erosion and deposition phenomena is a key feature to model granular flows such as landslides. Here, we investigated the instability that develops during the erosion of a wet granular pile by a dry dense granular flow. The morphology and the propagation of the generated steps are analyzed in relation to the specific erosion mechanism. The selected flowing angle of the confined flow on a dry heap appears to play an important role both in the final state of the experiment, and for the shape of the structures. We show that the development of the instability is governed by the inertia of the flow through the Froude number. We model this instability and predict growth rates that are in agreement with the experiment results.

Understanding sediment transport in Nature is essential to model landscape evolution, such as the transport processes in rivers and their formation [START_REF] Montgomery | [END_REF]2]. Patterns are known to spontaneously develop at the bottom of the river bed depending on the coupling between the flow and the bed geometry through the erosion/transport laws [3,4]. Granular flows such as landslides can also alter the underlying ground. However, contrary to the fluid case, the erosion mechanisms are less clear when the granular material is dense. Thus recent studies have focused on the evolution of the interface between an erodible layer and a granular flow [5][6][7], especially to model the entrainment rates, the velocities and the runout distances. Understanding such processes can also give insights on past and present climates for example [8][9][10]. Moreover, although the instabilities of riverbeds have been extensively characterized [11], the studies of instabilities in granular flows were generally focused on the flowing layer, thus describing avalanche fronts, roll waves, upward traveling waves [12][13][14][15], providing a better understanding of granular flows. As a result, little attention has been paid to the evolution of a granular bed under a solid mechanical load [16,17]. Finally, while previous studies on the erosion rate of a cohesive media by a granular flow focused on the physical properties of the cohesive materials (liquid bridges, tensile strength, elastic modulus) [7,18,19], the coupling between the erosion mechanism and the interface morphology with the granular flow is still an open question.

In this Letter, we explore the coupling between a flow of dry granular material and a cohesive granular bed. Our experiment exhibits an instability, as a train of steps appears along the initially flat interface. Although this instability shares similarities with the erosion of a dense cohesive bed by fluid flow [3], we show that this system presents unique features linked to the specific properties of granular matter. Using a granular rheology and an erosion law, we can model the phenomenon. The threshold and growth rates of the initial instabilities are in good agreement with the theory, providing insights into the erosion mechanism itself.

We performed experiments in a quasi-two-dimensional cell made of parallel vertical glass plates, whose gap width W was set at 6 and 12 mm. The same polydisperse glass beads of diameter 200-400 µm were used for both the cohesive granular bed and the flow. The cohesive material was obtained by mixing a small quantity of water with the beads to obtain a low water content of 1% wt, which corresponds to the formation of individual capillary bridges [20]. This material was introduced and packed into the cell to form a heap with a planar free-surface of length ranging from 30 to 35 cm. We injected the dry granular material at the top left of the cell so a flowing layer would from above the wet heap (Fig. 1a). The two remaining parameters are the initial inclination of the heap θ i and the constant injected flow rate Q. Both parameters influence the flow properties, mainly the stresses at the interface τ b and the mean velocity u.

As shown by the time evolution of the heap in figure 1a, erosion occurs at the interface between the cohesive pile and the granular flow, without deposition. Wet beads are slowly extracted by the action of the dry flow. Thus, the volume of the cohesive heap is gradually reduced during the experiment until it reaches a steady state in which the interface is flat and the dry granular material flows without eroding the wet bed. During this erosion process, two different interface morphologies can be observed. Either the eroded interface remains initially flat and parallel to the initial slope, or we observe the creation of step-shaped structures for rather high inclinations and high flow rates. These steps grow and propagate upstream, as shown in figures 1a,b, then disappear one after the other at the top of the channel. Once they all have disappeared, the erosion stops. In both cases, the final interface is flat and forms an angle θ f with the horizontal with θ f < θ i .

Figure 1b shows the evolution of the heap free surface with time in an unstable case. The first minutes of the development of the instability see the creation of many structures of small sizes, rapidly merged by a coarseninglike process. The surviving steps go through a growth phase lasting roughly 30 minutes, until they reach a stable shape (height and length): The back of the steps is not eroded, and appears as straight parallel lines while the front propagates upwards at a constant speed. We note that, close to the entrance and exit of the channel, the steps are influenced either by the injection point or by the accelerating flow. A single experiment may produce steps with different sizes, even once they are well established. However all the steps display a similar steep front, and a flat back with a well defined slope. Figure 2a shows the shapes of four different steps obtained for the same flow rate and different initial slopes. The front is well fitted by a parabola, whereas the back forms an angle θ f with the horizontal, which is equal to the final erosion angle. This means that the final equilibrium state is quickly reached locally, i.e. at the back of the steps. Moreover, the quantitative fitting of the parabola corresponds to a free fall under the gravity with an initial velocity of 22 cm.s -1 (solid red lines) which is comparable to the mean velocity u of the flow on the final slope (11.6 cm.s -1 ). This collapse of different step profiles also shows that the shape of the higher fronts are a simple continuation of the smaller ones (fig 2a). We interpret this similar shape as a consequence of the take-off of the granular flow occurring on the front of the steps. As the grains are almost free-falling after this point, there is no feedback from the lower part of the chute toward the upper part of the front of the steps through the dilute flow. This decoupling explains why the shape is kept, regardless of its size. It is worth noting that the erosion mechanism by this dilute flow is also probably different from the one initially occurring with the dense flow. This morphology and this propagation share similarities with the head-cut [21] or step and pool formation in riverbeds [3].

In the following, we focus on the final slope observed for both types of erosion morphology in this system. First let us comment on the stress evolution at the interface. For experiments driven at constant flow rate as these ones, one easily shows that the shear stress is increasing when the slope decreases using granular scaling laws (a result valid along an inclined plane [22], in presence of sidewall friction [23,24] or captured by empirical laws close to the flow threshold [START_REF] Deboeuf | Andreotti[END_REF]). Moreover, recent experiments have shown that the erosion of a wet granular surface results from a stochastic phenomenon [7]. These remarks rule out any intrinsic erosion threshold due to the wet pile. Therefore we measured the final angle θ f for different flow rates, in channels of different widths, using stainless steel and glass beads, while keeping the same erodible heap (wet glass beads). Concurrently to each experiment, we directly measured the flow slope in the very same conditions, injecting the same flow rate in an initially empty cell (without the wet pile). In this well-known granular-flow geometry on a dry heap (called heap flow), the flowing layer stabilizes over a static dry heap at a given angle depending on the flow rate [23,24]. Figure 2b shows that this heap-flow angle θ eq is equal to the final erosion angle, which proves that the erosion limit depends only on the flow properties. If the interface angle became smaller than the heap-flow angle, then a static dry layer would form between the wet grains and the flow preventing any further erosion process to occur.

This criterion thus explains the angle of the back of the steps. One last consequence of the localization of the erosion on the step fronts is the reduction of the average erosion rate up to a factor 2 (see Suppl. Mat. [26]). As already mentioned, the erosion of the cohesive heap does not always lead to the generation of steps. We explored under which conditions the interface remains flat. Due to the analogy with the shape and the propagation of the structures led by the pool and step instabilities in riverbeds [3], we expect to observe an inertial destabilization if the physical origin is identical. The underlying mechanism would be the phase lag between the topography and the local erosion rate due to inertia. Such phase lag exists for anti-dune or step-pool formations between the transport coefficient or the erosion and the shape of the interface [3,4] [START_REF]erosion (or transport) process is driven by the inertia of the flow: The shear stress is the link between erosion and inertia in water flows[END_REF]. We expect therefore a stability criterion based on the Froude number Fr = u/ √ gh cos θ i , which compares the inertia with the gravity effects, where g cos θ i is the projection of the gravity perpendicular to the slope. We first check the relevance of a critical Froude number by plotting the different states observed in the two channels in the Fr-θ i diagram. As shown by the gray strip in the figure 3a, our measurements show a critical Froude number, around 0.75, above which the instability develops. The mechanism is thus the delay of the local erosion rate with respect to the topography induced by the flow inertia.

Consequently, in the parameter space (θ i , flow rate) plotted in figure 3b, we can now draw the limits of the stable domain using iso-Fr lines. To compute the coordinates of these lines, we use the now established rheology for such dense dry granular flows µ(I) [START_REF] Midi | [END_REF], which links the internal friction coefficient to the inertial number. More precisely, assuming the flat interface to be rigid at the time scale of the flow, we follow the approach proposed by Jop et al. [24] to compute the vertical velocity profile: we integrate numerically the velocity profile of a steady 2D granular flow on a rigid rough plane with a noslip boundary condition, and including friction at sidewalls [26]. After adjusting the parameters of the model on other dry heap experiments, we are able to draw the limits in the phase diagram (Fig. 3b). The solid line corresponds to the equilibrium angle of the flow over a dry heap. The dashed line is an iso-Froude curve (Fr = 0.87, see later) which delimit the stable region.

To further understand the mechanism of the instability, we model the flow using a Saint-Venant approach in a long-wave approximation [29](see Suppl. Mat. [26]). In the following, the influence of sidewalls will only be taken into account in the shape of the vertical velocity profile in the basic state. We discard this frictional term in the equations for sake of simplicity and because h/W < 0.48 in the experiments [24]. The reference framework is the flat initial interface (the x-z axes in fig. 1a) eroded at constant rate < dη/dt >= -E 0 , which represents the normal velocity of the mean interface. The evolution of the thickness h, the elevation of the interface η and the mean velocity u are then described by three equations: the mass conservation inside the flow of constant volume fraction and the momentum equation, both integrated over the depth z, and the erosion law. We verified the scale separation between the characteristic times of the slow erosion process (few cm/hr) and of the flow (10 cm/s). This quasistatic evolution allows to neglect the mass flux into the flow due to erosion and the time derivative in the momentum equation:

ρ ∂(αhu 2 ) ∂x = ρgh cos θ tan θ -µ b (u, h) - ∂h ∂x , (1)
where α takes into account the vertical velocity profile (e.g. 5/4 for a Bagnold velocity profile or 4/3 for a linear velocity profile as in our case). The left term is the inertia, and in the right hand side, from left to right, the terms are the weight of the layer, the basal friction and the derivative of the pressure. The basal friction is modeled by the µ(I) rheology [12,[START_REF] Midi | [END_REF]. Following recent results on granular erosion [7], the erosion rate ∂η ∂t is proportional to -

Fcap F exp(- Fcap F )
, where F cap is the capillary force responsible for the cohesion and F is the mean force acting on grains. We used forces instead of stresses because the model was developed at the grain scale. In addition, we focus here only on the physical mechanism of the instability not on the influence of the cohesion on erosion [7]. As shown by the Fig. 3a,b, the erosion instability is linked to the inertia of the flowing grains, so F must depend on their inertia. A recent numerical study [19] proposes that the average force of impacting grains is linked to their kinetic energy. We thus assume that F is proportional to (ρu 2 ) β/2 with β = 1.2 [19]. The parameters are made dimensionless using the steady state erosion rate of the flat interface E 0 , the thickness h 0 , the flow rate per unit width q 0 , and the mean force F 0 : h = h/h 0 , x = x/h 0 , η = η/h 0 , ũ = uh 0 /q 0 , t = tE 0 /h 0 . After linearization around the initial steady state on a flat interface of slope θ i we decompose the flow into normal modes ∝ exp -i( kx + ωt ) and find the growth rate:

σ = (ω) = c k2 (αFr 2 -1) (a -b) 2 + k2 (αFr 2 -1) 2 , ( 2
)
where a and b are ∂µ b /∂u and ∂µ b /∂h [12][26], and c = β(F cap /F 0 -1). A critical Froude number is found above which the interface is unstable: Fr c = α -1/2 = 0.87 for α = 4/3. We use this value to plot the iso-Froude curve in figure 3b showing a good agreement with the experimental results. Moreover, we can compare the predicted initial growth rate to the experimental ones, measuring the height of the steps in the early exponential regime. Figure 3c shows the measured growth rates of the initial steps rising from the flat interface for different experiments. The most unstable mode is obtained The curves correspond to a prediction of the theoretical growth rate (eq. ( 2)) for the most unstable mode k = ∞ (dashed line) and for a wave number in agreement with the experiment h0k = 0.1 (solid line) at an inclination of θi = 30 o . Inset: Direct comparison of growth rates using the measured wave numbers k as a parameter.

with k = ∞ whose growth rate diverges at the threshold (dashed line). However, the model predictions for kh 0 = 0.1 and θ i = 30 o are in good agreement with the experimental data considering the error bars, noting that the experimental k is not strictly constant but close to 0.1. The inset of Fig. 3c shows the comparison when using the experimental k. We checked that choosing a different shape factor α for the steady-state velocity profile does not change qualitatively the graph. Our study indicate that inertia drives this instability. Whenever the granular flow is supercritical (Fr > Fr c ), the interface does not remain flat. This result is important to model correctly the evolution of the interface: if the mean erosion-force acting on the interface was linked to the weight of the layer or to the shear stress, the instability would not be predicted. Our findings strongly support the role of the impact of individual grains [19] in granular erosion mechanisms.

In conclusion, we have studied a new erosion instability occurring with granular flows over a cohesive granular bed. Above a critical Froude number, the interface is unstable and the shape of the steps is governed by the properties of the granular flows. We model with success the mechanisms in the framework of a Saint-Venant approach, extending the relevance of the depth-averaged methods to model complex evolutions of granular flows [29,30]. However, this model cannot predict any selection of a finite wavelength, additional longitudinal dissipation term in the momentum equation may be able to reduce the growth rate of too short wavelengths [29,30]. We finally identify that the inertia of grains governs the erosion mechanisms in dense granular flows. Further work is required to model the dilute erosion mechanisms in the fully developed state. We believe that our results could lead to a better understanding of geophysical flows.

FIG. 1 .

 1 FIG. 1. (a) Successive lateral views of a typical experiment. The shallow flowing granular layer appears brighter above the wet heap. The first frame is the initial state, and the others are snapshots taken during the erosion. The black arrow follows one step. The flow rate is 3.3 g/s and the initial inclination 35.8 o . The width of one picture corresponds to 27 cm. (b) Successive profiles of the interface, from dark blue to light red, every 15 minutes in the 6 mm channel, showing the steady propagation of steps over an initial slope angle of 39.4 o .

FIG. 2 .

 2 FIG. 2. (a)Superimposition of steps from a given experiment (black and blue dashed lines, solid green), and of a bigger step (dotted black), built from a higher initial slope angle: 43.7 o . The shifted red lines are a fitting curve composed of a straight line of slope 31 o followed by a parabolic curve corresponding to a free fall. (b) Angle of the heap at the end of the experiment θ f , versus the steady angle of flow of the same granular material down a dry heap θeq. Different flow rates are used to modify θeq, and the cohesive heap is made of glass beads.

FIG. 3 .

 3 FIG. 3. (a) Evolution of the interface with the initial Fr. The symbols depict the evolution of the interface. The borders of the gray strip are 0.7 and 0.8. (b) (θi-flow rate) phase diagram in the 12 mm-wide channel. The solid line is the heap flow angle θeq, under which erosion is not possible. The dashed line is the curve Fr = 0.87. (c) Dimensionless growth rates of the most unstable mode as a function of the Froude number in experiments for glass beads in the 6 mm channels (•) and 12 mm channel (•) and steel beads in 12 mm channel ( ). The curves correspond to a prediction of the theoretical growth rate (eq. (2)) for the most unstable mode k = ∞ (dashed line) and for a wave number in agreement with the experiment h0k = 0.1 (solid line) at an inclination of θi = 30 o . Inset: Direct comparison of growth rates using the measured wave numbers k as a parameter.

The authors would like to acknowledge B. Andreotti who pointed out the cyclic step analogy and to thank E. Dressaire for her help with proofreading.


Máximo Pérez-Romero

Enrique Romero-Cadaval

Adolfo Lozano-Tello

João Martins
email: jf.martins@fct.unl.pt

Rui Lopes
email: rm.lopes@campus.fct.unl.pt

An Innovator Nonintrusive Method for Disaggregating and Identifying Two Simultaneous Household Loads

Keywords: Monitoring, nonintrusive, disaggregation, identifying, simultaneously, S-Transform

At the present time the monitoring systems are important in some areas, such as electric power supply industry and household environment because they provide useful information to energy storage and management tasks. A new nonintrusive monitoring method is proposed in this paper and it is able to disaggregate and identify two loads working simultaneously using a single measuring sensor and a least squares regression algorithm based on discrete form of the S-Transform.

Introduction

The energy efficiency status report [1] showed by Joint Research Centre, or JRC, points out electricity consumption and efficiency trends in the EU-27 in 2012. According to the consumption statistics, final energy consumption in the EU-27 residential sector increased by 1.69%, furthermore the lowest consumption level of the last 20 years was reached in 2007. In terms of final residential energy consumption means 26,65% of total final energy consumption in the year 2010 showed in the breakdown of the Fig. 1. This sector was the second most consuming part of total consumption, so the residential field is important in energy efficiency programmes and policies.

Smart grid concept is a communication infrastructure to connect different distributed energy resources with three main targets, that is, by adapting the demand side to the distributed energy sources, transforming the distributed energy resources into systems providing active energy and supporting services, such as reactive power supply, to stabilize the distributed network and using storage devices in demand side [START_REF] Gregor | What the term Agent stands for in the Smart Grid Definition of Agents and Multi-Agent Systems from an Engineer's Perspective[END_REF]. Smart grids enable a two-way exchange of information and power between producers and consumers, and this leads to increase transparency and promote responsible energy saving measures on the consumer's side. It is important that consumers are involved to obtain effective smart electricity systems through trust, understanding and clear tangible benefits and consumers become proactive consumers or prosumers, being necessary to know individual loads on demand side. Load monitoring systems are important to manage smart grids because they provide accurate power and energy usage information to power usage prediction and management. There are two usual methods in literature to monitor appliances: intrusive and nonintrusive methods.

Intrusive methods use sensors in each appliance to measure electric signals, such as voltage and current waveform. This measuring manner is expensive because installation needs time and one sensor for each appliance. For overcoming this drawback researchers have studied other procedure to extract the individual energy information of each device. It is called Nonintrusive Load Monitoring (NILM) or Nonintrusive Appliance Load Monitoring (NIALM) and it consist of a single sensor placed in the main breaker level to measure current and voltage and an algorithm for disaggregate and identify each one of appliance. This question was firstly studied by George W. Hart [START_REF] Hart | Nonintrusive Appliance Load Monitoring[END_REF] of Massachusetts Institute of Technology (MIT) and published in the 90´s.

The methods for disaggregating different appliances can be divided into two groups: the steady state analysis and the transient state analysis. Regarding first group, it is based on the change in steady state power and requires low sampling rates. Hart´s work was already mentioned at the end of the previous paragraph. Next, Cole y Albicki [START_REF] Cole | Algorithm for nonintrusive identification of residential appliances[END_REF] have proposed an extension to the Hart method. It uses edges and slopes as appliance features together with steady-state power draw. Finally, Baranski and Voss have developed an approach [START_REF] Baranski | Non-Intrusive Appliance Load Monitoring Based on an Optical Sensor[END_REF], [START_REF] Baranski | Genetic Algorithm for Pattern Detection in NIALM Systems[END_REF], [START_REF] Baranski | Detecting Patterns of Appliances from Total Load Data Using a Dynamic Programming Approach[END_REF] for their system consisting of a histogram made up by historical data. As far as transient state analysis is concerned, transient characteristics can be used to disaggregate different household devices because they are tied to type of a given appliance. One of the most important researchs belongs to Leeb [START_REF] Leeb | Develop ment and Validation of a Transient Event Detector[END_REF], [START_REF] Leeb | Transient Event Detection in Spectral Envelope Estimates[END_REF]. An extension of work mentioned above consists of investigating the line voltage transients appeared when appliances are switched on or off (Cox and Leeb,[START_REF] Cox | Transient Event Detection for Nonintrusive Load Monitoring and Demand Side Manage ment Using Voltage Distortion[END_REF]). Finally Chang proposed a transient energy detector [START_REF] Chang | Load Identification in Neural Networks for a Nonintrusive Monitoring of Industrial Electrical Loads[END_REF].

A new method for disaggregating and identifying different household loads based on the S-transform [START_REF] Stockwell | Localization of the complex spectrum: the S transform[END_REF] is proposed in this paper (other S-transform method was used by [START_REF] Martins | A Novel Nonintrusive Load Monitoring System Based on the S-Transform[END_REF]). This work is part of [START_REF] Pérez-Romero | Local Energy Management Unit for Residential Applications[END_REF]. The novelty is the use of mathematical regression along with the S-transform and the development of algorithm to be implemented with the software LabVIEW 2012. In this way, this paper is based on such topic by disaggregating and identifying loads working individually or simultaneously and by sending such kind information in real-time to a decision central system.

S-Transform

The S-transform is a time-frequency analysis and representation technique proposed by Stockwell, Manshinha and Lowe in 1996. It is consists of a moving and scalable localizing Gaussian window leading to a time-frequency representation. Phase information is absolute because the origin of the time axis is taken as the fixed reference point.

The analysis begins with the target signal h(t). W(τ,f) is the Continuous Wavelet transform (CWT) of the signal using a specific mother wavelet and its equation appears in (1).

( ) ( ) ( , ) , . W f h t w t f dt τ τ ∞ -∞ = - ∫ ( 1
)
where t is the time, f is the frequency and w(t-τ,f) is the specific mother wavelet. Then

( ) ( ) 2 2
( , ) ( , ) , .

i f i f S f W f e h t w t f e dt π τ π τ τ τ τ ∞ - - -∞ = = - ∫ (2)
The mother wavelet is the first expression of (3).

( ) ( )

2 2 2 2 1 ( , ) ; ; ( ) 0. 2 t i ft w f g t e g t e f T f π σ λ τ σ λ σ π - - = = = = > (3)
The formula g(t) is a Gaussian window and the parameter σ means the Gaussian window width ,where λ sets frequency resolution. In this paper the parameter λ has the value 1. As the window width is the inverse of frequency, S-Transform results in a multiresolution analysis of the signal. Therefore, a high value of frequency gives high resolution of time and a low value of frequency gives high resolution of frequency.

Finally, by combining the expressions ( 4) and (2) the complete definition of the continuous S-Transform is given by [START_REF] Baranski | Non-Intrusive Appliance Load Monitoring Based on an Optical Sensor[END_REF].

( ) ( ) 2 2 2 2 ( , ) 2
t f i f t f w t f e e τ π τ τ π -- - - - = (4)
( )

( ) ( ) 2 2 2 2 2 ( , ) . 2 t f i f t i f f S f h t e e e dt τ π τ π τ τ π -- ∞ - - - -∞ = ∫ (5)
To efficiently implement on a computer the S-Transform it is necessary to develop the discrete form of the S-Transform. Stockwell et al. rewrote the continuous S-Transform by using Fourier transform H(α) of time series h(t) as the equation [START_REF] Baranski | Genetic Algorithm for Pattern Detection in NIALM Systems[END_REF] shows to get the discrete form. So, S-Transform can be run as computer algorithm by using the computational efficiency of Fast Fourier Transform and Inverse Fourier Transform in [START_REF] Baranski | Genetic Algorithm for Pattern Detection in NIALM Systems[END_REF].

From the discrete time series h[kT] corresponding to h(t) with a sample time T the equation of the discrete Fourier transform is given by [START_REF] Baranski | Detecting Patterns of Appliances from Total Load Data Using a Dynamic Programming Approach[END_REF].

( ) ( ) 2 2 2 2 2 2 ( , ) 0. i f S f H f e e d f π λ α πατ τ α α - ∞ -∞ = + ≠ ∫ (6)
[ ]

1 2 0 1 0,1,..., 1. N i nk N k n H h kT e n N NT N π - - =   = = -     ∑ (7)
where k and n = 0, 1, …, N -1.

Therefore, the discrete S-Transform is gotten by setting f → n/(NT) and τ → jT resulting in ( 8) and ( 9), with λ equal to 1.

N m i mj n N m n m n S jT H e e n NT NT π π - - = +     = ≠         ∑ (8)
[ ]

1 0 1 ,0 0. N m m S jT h n N NT - =   = =     ∑ (9)
where j, m and n = 0, 1, …, N -1.

The S-Transform consists of a complex matrix (called S-matrix) with dimension N by M where rows correspond to frequency and columns correspond to time.

Case Study

The proposed system consists of parts showed in Fig. 2: a monitor, a light bulb, a single current sensor of the measuring board, a point of common coupling (PCC), a control platform sbRIO 9631 of National Instruments consisting of a FPGA (hardware part) and a real-time processor (software part) and a notebook being only used to program the control platform. Overall system communicates with the central server called CEMIS.

The current sensor measures the current demanded by complete system and its information is sent to analog input of the control platform for disaggregating and identifying what loads are switched on or switched off. The two loads are connected in parallel to the PCC and one or both of them can be working simultaneously. The programming methodology of disaggregation and identification tasks, main topic of this paper, is based on codesign hardware-software and developed in LabVIEW 2012. The FPGA will get data from current sensor because the analog input is connected to it directly. It will also perform preprocessing tasks such as scale current data, convert from data type of fixed-point to unsigned word of 32 bits and write to a memory buffer. The sample frequency is 10 kHz and the number of points N is 201. The resultant matrix "S" is a transposed complex matrix with N rows and M = N/2+1 columns, that is, 201x102. The task of identification consists of the least squares regression [START_REF] Cox | Transient Event Detection for Nonintrusive Load Monitoring and Demand Side Manage ment Using Voltage Distortion[END_REF].

1 0 -1 1 0 -1 ( ) ( ) 2 1 2 1
.

N i i i i z a x a y e = - + = ∑ ( 10
)
where z i is the S-Transform of the overall measured current signal, x i is the S-Transform of the pattern signal of linear load (light bulb), Fig. 4(a), y i is the S-Transform of the pattern signal of electronic load (monitor), Fig. 4(b), a 1 and a 2 are the amplitudes of such pattern signals, respectively, and e is the least squares error. The coefficients a 1 and a 2 are obtained in ( 12) from ( 11) and they can be used in the time domain because the S-Transform meets the linearity criteria.

( ) .

1 2 1 1 2 0; N i i i i i e z a x a y x a = ∂ = - - = ∂ ∑ ( ) 1
N N N N N i i i i i i i i i i i i i BD CE AE CD a a A x B y C x y D z x E z y AB C AB C = = = = = - - = = = = = = = - - ∑ ∑ ∑ ∑ ∑ (12)

Experimental Results

On one side, for the case where a single light bulb is switched on the Fig. 5(a) shows signal fitted and computed by the system algorithm (white colour) and measured current signal of light bulb (yellow colour) and Fig. 5(b) shows the S-Transform of measured signal. The coefficients a 1 and a 2 have values 0,368337 and 0,0163964 ( Table 1), respectively, and identify a linear load as white colour signal of the Fig. 5(a) illustrates. Table 1 (this table illustrates the coefficients and error for each appliance and the combination of both.) and its waveform appears in Fig. 6(a) whereas its S-Transform is shown by Fig. 6(b). Finally, the main case consists of one light bulb and monitor working simultaneously. The system calculates the values of coefficients a 1 and a 2 and, according to the gotten results, it can be said that two loads have been disaggregated and identified. Fig. 7

Conclusion and Further Work

The main feature of this new nonintrusive method is to disaggregate and identify two loads, a linear device and nonlinear device, by using S-Transform and least squares regression with one cycle of the time series signal. It is an advantage because the S-Transform uses computational efficiency of FFT and IFFT and other parameters are not necessary, such as harmonic analysis and reactive power, reducing computations. Furthermore, the system needs a single current sensor reducing cost. In the future, more loads will be connected to the system and the least squares regression will extend to detect such devices although the complexity will increase. There is also an open work to identify devices switched on by using neural networks.

Fig. 1 .

 1 Fig. 1. Final energy consumption breakdown into sectors in the EU-27, 2010 (Eurostat).

Fig. 2 .

 2 Fig. 2. System diagram.

Fig. 3 (

 3 Fig. 3(a) shows the FPGA algorithm flowchart. The software part will read the data from buffer, convert from unsigned word of 32 bits to data type of double precision, extract one signal cycle, execute the S-transform and disaggregate and identify the loads by means of least squares regression as it is shown in the flowchart of Fig. 3(b).

Fig. 3 .

 3 Flowcharts: (a) FPGA algorithm, (b) Processor algorithm.

Fig. 4 .

 4 Fig. 4. (a) Current signal of light bulb. (b) Current signal of monitor.

Fig. 5 .

 5 (a) Measured and fitted current signal of light bulb (front panel of LabVIEW 2012). (b) S-Transform of measured current signal.On the other hand, when a single electronic load or monitor is turned on the system discovers a nonlinear load currently working. The coefficients are shown in

  (a) shows time series signals and Fig. 7(b) presents the S-Transform.

Fig. 6 .Fig. 7 .

 67 (a) Measured and fitted current signal of monitor (front panel of LabVIEW 2012). (b) S-Transform of measured current signal. (a) Measured and fitted current signal of light bulb and monitor (front panel of LabVIEW 2012). (b) S-Transform of measured current signals.

2 Relationship to Collective Awareness Systems The

  Collective Awareness Platforms are ICT systems taking advantage of distributed knowledge and data from real environments, that is, Internet of things, to generate open databases. These types of platforms are scalable resulting in a big data source for different applications, such as real-time household power consumption data for managing energy demand of consumers, becoming collective tools.

Table 1 .

 1 Values of coefficients.

	LOAD	a 1	a 2	Error
	Light bulb	0,368337	0,0163964	0,003888227
	Monitor	0,0450099	0,309192	0,220845
	Light bulb and monitor	0,409668	0,321245	0,181646

2 2

Acknowledgements. This work has been developed under support of Telefónica Chair of University of Extremadura.


F A Barata
email: fbarata@deea.isel.ipl.pt

Rui Neves-Silva

Distributed MPC for Thermal Comfort in Buildings with Dynamically Coupled Zones and Limited Energy Resources

Keywords: Multi-zone thermal comfort, electro-thermal analogy, DMPC, limited energy resource

This paper presents a distributed predictive control methodology for indoor thermal comfort that optimizes the consumption of a limited energy resource using a demand-side management approach. The building divisions are modeled using an electro-thermal modular scheme. For control purposes, this modular scheme allows an easy modeling of buildings with different plans where adjacent areas can thermally interact. The control objective of each subsystem is to minimize the energy cost while maintaining the indoor temperature in the selected comfort bounds. In a distributed coordinated environment, the control uses multiple dynamically coupled agents (one for each subsystem/zone) aiming to achieve satisfaction of available energy coupling constraints. The system is simulated with two zones in a distributed environment.

Introduction

Buildings are responsible for a large share of our global energy use. Energy use is in fact the main determinant of a buildings global environmental footprint, considering its total life span.

Reducing energy consumption in the building stock is a trend in the world today, partly because of economic, partly because of environmental reasons. Energy dissipation depends on the construction of house, the materials from which the house is built, insulation during the year, outside temperature, and additional sources of energy. Consumption of energy is predominantly determined with a selection of materials and architectural solutions, and it can be further reduced with efficient management of heating or cooling. An effective heating/cooling management is provided in the framework of predictive control, particularly Model Predictive Control (MPC), has been granted to reduce and optimize the energy consumption in the residential sector namely to deal with temperature set points regulations [START_REF] Moroşan | Building temperature regulation using a distributed model predictive control[END_REF][START_REF] Ma | Predictive Control for Energy Efficient Buildings with Thermal Storage[END_REF][START_REF] Freire | Predictive controllers for thermal comfort optimization and energy savings[END_REF] and load management [START_REF] Giorgio | A Model predictive Control Approach to the Load Shifting Problem in a household Equipped with an energy Storage Unit[END_REF][START_REF] Zong | Application of predictive control for active load management in a distributed power system with high wind penetration[END_REF].

The distinct advantages of MPC based control solutions compared to classical controllers are: using also relevant future information in making control decisions using predicted profiles (e.g. ambient temperature and solar irradiation on outside walls), routine handling of multi-input-multi-output (MIMO) systems; routine respecting of system constraints (e.g. finite amount of heating/cooling power in room possible, desired room temperature spans); and explicit orientation of control actions towards the goal which can be using least energy possible, or spending least money possible, or causing least CO 2 emissions possible, or their combination.

Optimal control for indoor environment requires preservation of comfort conditions for buildings occupants and minimization of energy consumption and cost [START_REF] Barata | Model Predictive Control for Thermal House Comfort with Limited Energy Resources[END_REF]. Basically, MPC makes a tradeoff between energy savings and thermal comfort. The MPC has also advantage when controlling distributed systems [START_REF] Negenborn | Multi-Agent Model Predictive Control with Applications to Power Networks[END_REF], [START_REF] Scattolini | Architectures for distributed and hierarchical Model Predictive Control -A review[END_REF]. Distributed Model Predictive Control (DMPC) algorithms are the state of the art in complex control problems with many interconnected subsystems. DMPC allows the distribution of decision-making while handling constraints in a systematic way. DMPC strategies can be characterized by the type of couplings or interactions assumed between constituent subsystems [START_REF] Trodden | Distributed model predictive control of linear systems with persistent disturbances[END_REF].

The method of subsystems sharing coupled constraints can be seen in [START_REF] Trodden | Distributed model predictive control of linear systems with persistent disturbances[END_REF], [START_REF] Keviczky | Decentralized Receding Horizon Control for Large Scale Dynamically Decoupled Systems[END_REF] being the strategy here presented a DMPC with coupled constraints (renewable energy must be shared by all divisions) and dynamically coupled zones. Thus, in a distributed coordinated environment, the control uses multiple dynamically coupled agents (one for each subsystem/division) aiming to achieve satisfaction of coupling constraints.

The desired approach here presented intends to take advantage from the innovative technology characteristics provided by future Smart Grids (SGs) [START_REF] Siano | Demand response and smart grids-A survey[END_REF]. In the smart world, simple household appliances, like dishwashers, clothes dryers, heaters, air conditioners will be fully controllable in order to achieve the network maximum efficiency.

Compared with the aforementioned literature, the novel contributions of this work are related with the existence of a system with coupled constraints and dynamically coupled zones in a cost function with distinct objectives allowing thermal comfort with a consumption to weatherize the divisions inside the available power constraints.

The remaining of this paper is organized as follows: Section 2 presents the technological contribution of this paper and Section 3 presents the system architecture. In Section 4 is introduced the house dynamic model and the model predictive controller is presented in Section 5. Some results and analysis are shown in Section 6 and conclusions are drawn in Section 7.

Relationship to Collective Awareness Systems

Collective awareness is an important issue to the use of groupware systems and virtual collaboration. From the Collective Awareness Systems (CAS) perspective, distributed networks with multi-agents may represent systems that are able to collectively contribute to a same objective. The analysis in this paper is on a line of technological contributions related with the CAS characteristics. Future appliances will be linked in the grid and will be fully controllable, monitored and regulated in real time in order to collaborate to achieve energy efficiency, balance between demand and supply, intelligent load control consumer comfort and CO 2 emissions reduction.

System Architecture

In this paper the large scenario considers a distributed network that involves a residential community, with electricity power source generated by their own renewable energy park. Hence, two kinds of energy are considered, the green that is from the renewable source and the red energy, from the grid. The energy from the grid is always available, although at a higher price, and it is only consumed when the green energy is not enough to satisfy the demand. The houses/divisions may have different plans, thermal loads, thermal characteristics, occupancy and comfort temperature bounds, and consequently with different energy needs for heating/cooling the spaces. The idea is to apply a predictive control law to maintain the temperature and power consumption inside their bounds. The green energy is limited and predictable and must be shared by all houses/divisions. The red resource consumed for comfort implies a penalty in the final cost function [START_REF] Barata | Model Predictive Control for Thermal House Comfort with Limited Energy Resources[END_REF] due to the soft constraint violation. This penalization means that the maximum available green resource was exceeded. It is considered that the outside temperature, disturbances and daily comfort temperature bounds are known by each system inside the predictive horizon (N). The green resource that is not consumed at a certain instant is stored in batteries or delivered to the grid.

Thermal Model of the House

The idea here presented is to apply the principle of analogy between two different physical domains that can be described by the same mathematical equations. Thus, a linear electrical circuit represents the building and the state-space equations are obtained by solving that circuit. Here, the temperature is equivalent to voltage, the heat flux to current, the heat transmission resistance is represented by electrical resistance and the thermal capacity by electrical capacity. The equivalent circuit of the building is obtained by assembling models of the walls, windows, internal mass, etc. In the case of single-zone buildings, interior walls are being part of the internal thermal mass while exterior walls are forming the building envelope. Several approaches can be seen in [START_REF] Hazyuk | Optimal temperature control of intermittently heated buildings using Model Predictive Control: Part I e Building modelling[END_REF][START_REF] Zong | Application of predictive control for active load management in a distributed power system with high wind penetration[END_REF] where is shown that building models can be simpler or more complex depending on the objective. Fig. 2 shows the used thermal-electrical modular approach. The model presented in (1-3) is a low order model describing the dominant dynamics for division i [START_REF] Ma | Predictive Control for Energy Efficient Buildings with Thermal Storage[END_REF] with adjacent area, which can be considered suitable for control proposes.

( )

i i i Pd losses heat i i Q Q Q C dt dT + - = 1 , (1)
∑ = - + - = n j ij i j eq i out losses R T T R T T Q i i 1 , (2)
i th wall window i roof eq R R R R R i i i + = // // , with ∑ = materials i window window R R , ∑ = materials i roof roof R R and ∑ = materials i wall wall R R (3)
For several adjacent zones (1) can be generically written as follows, The plant model representation (1) for one division with an adjacent zone can be approximated by a discrete model using Euler discretization [START_REF] Luyben | Process Modeling, Simulation and control for Chemicals Engineers, 2ºedition[END_REF] with a sampling time of ∆t.

                                +                     + + + +                                 - - - - - - =               n n n n n eqn out eq out eq out n n n n n n eqn n n n n n n eq n n eq n u u u C C C C Pd C R T C Pd C R T C Pd C R T T T T C R C R C R C R C R C R C R C R C R C R C R C R C R C R C R T T T M L O M M M L L M M L M L M M O L & M & & 2 1 2 1 2 2 2 2 1 1 1 2 1 2 1 2 1 2 2 2 2 2 21 2 2 2 21 1 1 1 12 1 1 1 12 1 2 . 1 1 0 0 0 1 0 0 0 1 1 1 ... 1 1 1 1 1 ... 1 1 1 1 1 1 ... 1 1 1 1 (4) Tout T1, C1, Qheat1, Req1 Qlosses1, QPd1 R12 T2, C2, Qheat2, Req2 Qlosses2, QPd2 T3
) ( ) ( ) ( ) ( ) 1 ( k d k T D k u B k T A k T i j i i i i i i + + + = + , ( 5
)
where

i eq oa i d i ij i j i i i ij i eq ij eq i C R t T C t P d t R C T D C t B t R C R R R A i i i i ∆ + ∆ = ∆ = ∆ =         ∆ + - = , , , 1
, u 1 (k) is the necessary heat/cooling power, T i (k) is the indoor temperature, d i (k) is a disturbance signal resulting from P di the external disturbances (kW) (e.g. load generated by occupants, direct sunlight, electrical devices or doors and windows aperture to recycle the indoor air), and T oa , the temperature of outside air (ºC).

Model Predictive Control Cost Function

Fig. 3 shows the implemented MPC scheme. The controllers from the areas that are thermally coupled interchange information about their state prediction as can also be seen in [START_REF] Negenborn | Multi-Agent Model Predictive Control with Applications to Power Networks[END_REF]. At each time step, each one of the agents must solve his MPC problem. The cost function objectives are: minimize the energy consumption to heating and cooling; minimize the peak power consumption; maintain the zones within a desired temperature range and maintain the used power within the green available bounds. The generic problem to be solved by each agent, assumes the following form:

{ } ( ) ( ) ∑ ∑ ∑ = + + = + + - + - = + + + + + + ∆ N k t k t t k t N k t k t t k t t N t t t N k t k t u u t u γ , γ , ε , ε U, 1 2 | 2 | 1 2 | 2 | 2 | 1 2 | 1 0 2 | ,..., max min γ γ ψ ε ε ρ φ (6)
subject to the following constraints,

t k t t k t j t k t t k t t k t d DT Bu AT T | | | | | 1 + + + + + + + + + = , (7)
t k t t k t T T T | | + + + ≤ ≤ - ε ε , ( 8
)
t k t A t k t A i i U U U | | + + + ≤ ≤ - γ γ , (9) 0 , , , | | |
| ≥ + + + + t k t t k t t k t t k t ε ε γ γ . ( 10
)
In ( 6), ݑ represents the power control inputs, φ is the penalty on peak power consumption, ρ is the penalty on the comfort constraint violation, ψ the penalty on the power constraint violation and N is the length of the prediction horizon. In [START_REF] Scattolini | Architectures for distributed and hierarchical Model Predictive Control -A review[END_REF], ߝ and ߝ are the vectors of temperature violations that are above and below the desired comfort zone defined by T and T . In [START_REF] Trodden | Distributed model predictive control of linear systems with persistent disturbances[END_REF], ߛ and ߛ are the power violations that are above or lower the maximum,

Results

The presented results were obtained with an optimization Matlab routine. As a first approach towards developing a control structure it is considered an individual house with two divisions thermally coupled. The used parameters that characterize the house are showed in Table 1. It is considered that the divisions have the same thermal characteristics and penalties in the cost function. The planned thermal perturbations are, however, different for the two divisions (Table 1 and Fig. 4(b)). The division that uses the available green energy first is Division 1 (D1) and Division 2 (D2) uses only the remainder. By this reason, the maximum available energy to D1 is always the maximum green available stock, and for D2 is given by (11): The outdoor temperature forecasting (T oa ) has the profile present in [START_REF] Barata | Distributed Model Predictive Control for Thermal House Comfort with Auction of Available Energy[END_REF]. The comfort limits and available renewable resource vary during the 24 period and it can be seen that both indoor temperature and consumed power are always maintained inside the constrained bounds (Fig. 5 and6). Taking advantage of the predictive knowledge of the thermal disturbance and making use of the space thermal storage, it can also be seen that in both divisions the MPC changes the indoor temperature in anticipation to the thermal disturbance (Fig. 5 For each one of the divisions it can be seen in Fig. 7 that the "Real Cost" is much lower than the cost of only consuming the red resource "Red Cost" at a higher fixed price. To illustrate the benefit of the power constraint penalization in (6), Fig. 7(b) shows that maintaining all the other features and changing only the penalty value D1 had a higher cost compared with Fig. 7(a).

) : ( ) : ( ) : ( 2 1 2 N k k u N k k U N k k U A A + - + = + . (11)

Conclusions

In this paper, a distributed MPC control technique was presented along with a thermal-electrical modular scheme in order to provide thermal house comfort. The solution obtained solves the problem of control of multiple subsystems dynamically coupled subject to a coupled constraint. Each subsystem solves its own problem by involving its own and adjacent rooms state predictions and also the shared constraints. Changing the penalty values, the consumer can choose in each division between indoor comfort and lower costs. It could be observed through the simulations and results analysis that suitable dynamic performances were obtained.

Fig. 1 .

 1 Fig. 1. Global scenario.

Fig. 2 .

 2 Fig. 2. Generic schematic representation of thermal-electrical modular analogy for several divisions.

Q

  cooling losses (kW), i T the inside temperature (ºC), i C the equivalent thermal capacitance (kJ/ºC), and the external thermal disturbances (kW) (e.g. load generated by occupants, direct sunlight, electrical devices or doors and windows aperture to recycle the indoor air). In (2) out T is the outdoor temperature (ºC), ij R the thermal resistance between divisions,

Fig. 3 .

 3 Fig. 3. Block diagram of the implemented MPC system in each division.

Fig. 4 .

 4 Fig. 4. (a) outdoor temperature forecasting (T oa ); (b) thermal disturbance profile (P d ).

Fig. 5 .

 5 Fig. 5. Division 1 (a) power profile; (b) indoor temperature profile.

Fig. 6 .

 6 Fig. 6. Division 2 (a) power profile; (b) indoor temperature profile.

  (b) and 6(b)).

Fig. 7 .

 7 Fig. 7. (a) Cost profile with ψ=500 (b) Cost profile with ψ=1.

Table 1 .

 1 Scenario parameters.

C(

kJ/ºC) Rroof (ºC/W) Rwalls(ºC/W) Rth(ºC/W) Rwindows(ºC/W) Rtotal(ºC/W) R12(ºC/W)

	9.2×10 3	0.192	0.031	0.038	0.023	0.050	0.030
	Divisions	ρ ρ ρ ρ	ψ ψ ψ ψ	Φ	∆ ∆ ∆ ∆t	N	T(0) (ºC)
	D1	100	500	2	1	24	22
	D2						23