Christian Ghiaus
email: christian.ghiaus@insa-lyon.fr

Florent Alzetto

Design of experiments for Quick U-building (QUB) method for building energy performance measurement

Keywords: QUB method, model identification, modelling and measurement error, global thermal conductance, exponential matrix, eigenvalues and eigen vectors, time constants

Quick U-building (QUB) is a method for short time measurement of energy performance of buildings, typically one night. It uses the indoor air temperature response to power delivered to the indoor air by electric heaters. This paper introduces a method for estimating the expected measurement error as a function of the amplitude and the time duration of the input signal based on the decomposition of the time response of a state-space model into a sum of exponentials by using the eigenvalues of the state matrix. It is shown that the buildings have a group of dominant time constants, which gives an exponential response, and many very short and very large time constants, which have a small influence on the response. The analysis of the eigenvalues demonstrates that the QUB experiment may be done in a rather short time as compared with the largest time constant of the building.

Δ𝑇 𝑐

-difference between the indoor and outdoor temperatures during cooling [°C] 𝛥𝑇 0 𝑐 -difference between the indoor and outdoor temperatures during cooling at time 𝑡 0 [°C] Δ𝑇 ℎ

-difference between the indoor and outdoor temperatures during heating [°C] 𝛥𝑇 0 ℎ -difference between the indoor and outdoor temperatures during heating at time 𝑡 0 [°C] Vectors and matrices 𝐟 -vector of sources of heat flow rate 𝐛 -vector of sources of temperature 𝐮 -vector of inputs 𝐱 -vector of states

Introduction

The areal heat transfer coefficient is an intrinsic indicator of the energy performance of buildings in steady-state. Therefore, the measurement of the heat transfer coefficient of residential buildings, especially in a short period of time (less than a day), is of particular interest.

The thermal transmittance of a flat wall (or the U-value) is defined as the "heat flow rate in the steady state, 𝑄 ̇, divided by area, 𝐴, and by the temperature difference between the surroundings of both sides of a flat uniform system, Δ𝑇 (ISO 7345, 2018):

𝑈 = 𝑄 Ȧ Δ𝑇 (1)
Heat transfer coefficients are defined for the heat transmission to external environment, to adjacent spaces, and to the ground (ISO 13790, 2008;ISO 52016-1, 2017;ISO 13789, 2017).

In steady-state, the U-value of the building elements can be obtained by dividing the heat flow rate over the temperature difference on the two sides of a wall (ISO 9869-1, 2014;ASTM C 1046ASTM C -95, 2001;;ASTM C 1155ASTM C -95, 2001;;Janssens, 2016). This method is valid if:

-the thermal properties of the materials and the heat transfer coefficients do not change during the measurement; -the variation of the heat stored in the wall is negligible as compared with the heat transferred through the wall [START_REF] Thébault | Refinement of the ISABELE method regarding uncertainty quantification and thermal dynamics modelling[END_REF].

Similarly, the areal heat transfer coefficient, 𝐻′, is defined as "the heat flow rate between the building (i.e. the internal environment) and the external environment, 𝑄 ̇, divided by the area, 𝐴, and the difference of temperature between the internal and external environment, Δ𝑇" (ISO 7345, 2018):

𝐹 𝑆 ≡ 𝐻 ′ = 𝑄 Ȧ Δ𝑇 (2)
Therefore, the heat transfer coefficient of the building is:

𝐻 = 𝑄 Δ𝑇 (3)
Since steady-state conditions are never achieved for whole buildings, the methods used for measuring the overall thermal transmittance are based on one or a combination of the following procedures [START_REF] Deconinck | Comparison of characterisation methods determining the thermal resistance of building components from onsite measurements[END_REF]:

-imposing steady-state, e.g. hot-box method for a whole wall (ISO 8990, 1994); -considering mean values over a long period of time, longer than three days and multiple of 24h (ISO 9869-1, 2014); -correcting the steady state method for storage effects (ISO 9869-2, canceled); -identifying model parameters by using long-term experiments. [START_REF] Fels | PRISM: An Introduction[END_REF] [START_REF] Fels | PRISM: An Introduction[END_REF] Steady-state or dynamic models are used for parameter identification in long term experiments, typically a season or a whole year. In the first approach, the areal heat transfer coefficient is determined from the energy signature which represents the relation between the consumption and the temperature difference between indoors and outdoors. Basically, the areal heat transfer coefficients are determined by regression [START_REF] Fels | PRISM: An Introduction[END_REF][START_REF] Hitchin | Daily energy consumption signatures and control charts for airconditioned buildings[END_REF][START_REF] Vesterberg | A single-variate building energy signature approach for periods with substantial solar gain[END_REF]Danov, Carbonell, Cipriano, & Marti-Herrero, 2013;Ghiaus, 2006) or robust regression (Ghiaus, 2006). In the second approach, the areal heat transfer coefficient is estimated from the coefficients of dynamic models [START_REF] Jimenez | Models for describing the thermal characteristics of building components[END_REF][START_REF] Roels | On site characterisation of the overall heat loss coefficient: Comparison of different assessment methods by a blind validation exercise on a round robin test box[END_REF][START_REF] Naveros | Setting up and validating a complex model for a simple homogeneous wall[END_REF][START_REF] Raillon | An efficient Bayesian experimental calibration of dynamic thermal models[END_REF]Janssens, 2016;[START_REF] Naveros | Physical parameters identification of walls using ARX models obtained by deduction[END_REF].

Since measuring the performance in quasi steady-state requires typically a few weeks, shorter-time methods for unoccupied houses were investigated. In these methods, the measurements are done when the electrical appliances and other internal heat gains are zero, with the exception of the electrical energy used for the experiment, which is counted. A widely used method for measuring the overall heat transfer coefficient of a house is co-heating. The heat transfer coefficient of the whole building is determined, including the air renewal rate, by steady-state measurement and regression analysis (Bauwens & Roels, 2014). The monitoring is done typically for 5 -10 days with measurements at a time steps of 40 -60 min. Variants include measurement of losses through air infiltration, solar aperture, and ground floor.

PSTAR is a unified method for building simulation and short-time tests [START_REF] Subbarao | Short-Term Energy Monitoring (STEM)[END_REF]. It uses energy balance equations, one-time measurements, and hourly measured data in order to obtain a description of the building at macrolevel. The main idea is to obtain a quasi-steady state by maintaining the temperature constant for several hours through the night in order to reduce at insignificancy the effect of thermal storage. The tests are done for about 36 hours during two nights, one night for coheating and another one for cooling down. The thermal conductance is measured during the last two hours. The method uses the measurement of the air-infiltration rate. The main issue with this approach is that it assumes that steady-state is obtained in a rather short time (one night).

Another class of methods is based on finding the parameters of the models from the response to excitation pulses. ISABELE method uses least squares identification of the parameters of a first order dynamic model obtained from the response to a step cooling (Bouchié R. , et al., 2015;Bouchié R. , Alzetto, Brun, Boisson, & Thebault, 2014;[START_REF] Thébault | Refinement of the ISABELE method regarding uncertainty quantification and thermal dynamics modelling[END_REF]. The experimental protocol has three steps:

-no heating power injected in order to assess the initial thermal energy stored in the building; if the building has quasi-constant indoor temperature for several weeks, this step is not required; -heat the building for at least two days to homogenous internal temperature having values at least 10 °C larger than the temperature of the first step; -turn off the heater for at least two days and measure the free-running temperature variation in time.

The response factors method measures the global resistance of a building element by identifying the time series parameters of the response to a rectangular or triangular pulse [START_REF] Rasooli | A response factor-based method for the rapid insitu determination of wall's thermal resistance in existing buildings[END_REF]Brisken & Mitalas, 1956). The triangular shape of the excitation in heating was obtained by heating up to 70 -90 °C with a radiative heater for 15 min; in cooling, the radiative heater was removed, and a fan and an ice bag were used. A main requirement is that the temperature on the outdoor side of the wall remains constant; therefore, the outdoor wall is insulated from the exterior conditions by using a box in which the temperature is controlled [START_REF] Rasooli | A response factor-based method for the rapid insitu determination of wall's thermal resistance in existing buildings[END_REF].

QUB method consists on a positive step input followed by a zero step, i.e. the excitation is rectangular. The outdoor temperature is considered constant during the experiment. The parameters of a first order dynamic model are deduced from these measurements. This method measures the heat transfer coefficient, 𝐻, of a building in a rather short time, usually one night. QUB method was tested on many experimental and real houses. Some of these results have been already published:

• on a house constructed inside an environmentally controlled chamber at University of Salford [START_REF] Meulemans | QUB/e: A Novel Transient Experimental Method for in situ Measurements of the Thermal Performance of Building Fabrics[END_REF]Alzetto, Pandraud, Fitton, & Heusler, 2018;Alzetto, Farmer, Fitton, Huges, & Swan, 2018); • on one of the twin houses of Fraunhofer Institute for Building Physics (Alzetto, Pandraud, Fitton, & Heusler, 2018); • on a house of 280 m 2 in southern France [START_REF] Mangematin | Quick measurements of energy efficiency of buildings[END_REF];

• on an apartment in a high-raise building in Stockholm [START_REF] Meulemans | An assessement of QUB/e method for fast in situ measurements of thermal performance of building fabrics in cold climates[END_REF].

The response factors, ISABELE and QUB methods have similarities in the experimental procedure: short time experiment, rectangular or triangular excitation, constant outdoor temperature. However, there are important differences. In the measurement procedure, response factors are designed for a building element (wall) and considers the heat flow rate and surface temperatures while ISABELE and QUB methods are designed for the whole building and consider the equivalent indoor air temperature. In the method used for interpreting the results, response factors method uses experimental values to find out the parameters of the time series, ISABELE identifies the parameters of a model with five resistances and one capacitance, while QUB method uses slopes to find the parameters of an exponential response.

A common issue not treated in these methods is the design of experiments: what is the amplitude and the duration of the excitation to be used in order to obtain minimum experimental errors. This paper responds to this question for QUB method. A second issue is the justification of relatively short time experiment as compared with the largest time constants. This paper treats this aspect by analysing of the eigenvalues of a detailed simulation model.

Principle of QUB method

The principle of QUB method is to apply two power pulses to the building and to determine the global conductivity and effective capacity from these measurements by finding the parameters of a first order dynamic model [START_REF] Mangematin | Quick measurements of energy efficiency of buildings[END_REF] (Figure 1). The model on which the measurements are fitted is: The two experiments conducted in QUB method are: 1) Heating: all the rooms are heated at roughly the same temperature with controlled electrical heaters; the total injected power is 𝑃 = 𝑃 ℎ .

2) Free-cooling: the building is left to cool. The power during the cooling period, 𝑃 = 𝑃 𝑐 , is from interruptible sources, such as internet connection and measurement devices (Alzetto, Pandraud, Fitton, & Heusler, 2018).

In QUB method, the two parameters, 𝐻 𝑄𝑈𝐵 and 𝐶, of the model (4) are obtained from the two slopes, 𝛼 𝑐 and 𝛼 ℎ given by eq. ( 62) and eq. ( 55), respectively:

𝐻 𝑄𝑈𝐵 = 𝑃 ℎ 𝛼 𝑐 -𝑃 𝑐 𝛼 ℎ 𝛥𝑇 0 ℎ 𝛼 𝑐 -𝛥𝑇 0 𝑐 𝛼 ℎ (5)
and

𝐶 = 𝑃 ℎ 𝛥𝑇 0 𝑐 -𝑃 𝑐 𝛥𝑇 0 ℎ 𝛼 ℎ 𝛥𝑇 0 𝑐 -𝛼 𝑐 𝛥𝑇 0 ℎ (6)
When the slopes are taken after time 𝑡 0 , their expressions are given by (see Annex 1)

𝛼 ℎ = 𝑒 -𝑡 0 𝜏 (𝑃/𝐶 -Δ𝑇 0 ℎ 𝐻 𝑄𝑈𝐵 /𝐶) (7)
and

α c = 𝑒 -𝑡 0 𝐻 𝑄𝑈𝐵 𝐶 ( 𝑃 𝑐 𝐶 - Δ𝑇 0 𝑐 𝐻 𝑄𝑈𝐵 𝐶 ) (8)
Therefore

𝐻 𝑄𝑈𝐵 * = 𝑃 ℎ 𝛼 𝑐 -𝑃 𝑐 𝛼 ℎ 𝛥𝑇 0 ℎ 𝛼 𝑐 -𝛥𝑇 0 𝑐 𝛼 ℎ = 𝐻 𝑄𝑈𝐵 (9)
and

𝐶 * = exp (-𝑡 0 /𝜏) 𝑃 ℎ 𝛥𝑇 0 𝑐 -𝑃 𝑐 𝛥𝑇 0 ℎ 𝛼 ℎ 𝛥𝑇 0 𝑐 -𝛼 𝑐 𝛥𝑇 0 ℎ = exp (-𝑡 0 /𝜏)𝐶 (10)
Note that 𝐻 𝑄𝑈𝐵 has the same expression regardless of the time at which the slope is considered, in origin or at time 𝑡 0 > 0, while the value of 𝐶 depends on the time at which the slope is calculated. In order to find the heat capacity 𝐶 * when slopes at time 𝑡 0 > 0 are considered, the transcendental equation (10) needs to be solved. Equation ( 9) is important because it shows that the overall heat transfer coefficient, 𝐻 𝑄𝑈𝐵 , can be estimated based on the slopes measured at moments different from the initial time. As shown hereafter, the buildings have a dominant time constant, which gives an exponential response, and many short and large time constants, which have a small influence on the response. If the measurement is done after a time which corresponds to the steady-state for the short time constants, the time response is an exponential. Therefore, QUB method allows us to measure the global conductance by measuring the slopes of this exponential response at moments different of the initial time, long enough to settle the response due to short time constants and short enough to not be influenced by the very large time constants.

Optimal design of experiments for QUB method

Knowing the outdoor temperature 𝑇 𝑜 and the power in cooling phase 𝑃 𝑐 , the aim of the optimal design of experiments for QUB method is to find the heating power 𝑃 ℎ and the heating time 𝑡 𝑄𝑈𝐵 which minimize the total error of the measured heat transfer coefficient, 𝐻 𝑄𝑈𝐵 . The total error is composed by the intrinsic error of the method and the measurement error. The proposed procedure is:

-From a thermal network of the building obtain its state-space model.

-From the state-space model, calculate the overall conductance 𝐻 ≡ 𝑅 -1 (an exact value) that will be considered as reference.

-Introduce the boundary conditions (temperatures and heat flows) and initial conditions and calculate 𝐻 𝑄𝑈𝐵 by simulating the QUB experiment on the detailed model.

-Draw a contour map of the error 𝜀 𝑄𝑈𝐵 = 𝐻 𝑄𝑈𝐵 -𝐻 𝐻 as a function of the heating power 𝑃 ℎ , the time duration 𝑡 𝑄𝑈𝐵 , and the measurement errors. -Chose an experiment which minimizes the error and respects technological constraints like maximum available power and admissible temperature. -Compare the experimental response of the building with the response of the model.

Total error

The total error depends on the intrinsic error of the method and on the error coming from measurements. The intrinsic error coming from QUB method is

𝜀 𝑄𝑈𝐵 = 𝐻 𝑄𝑈𝐵 -𝐻 (11)
where 𝐻 𝑄𝑈𝐵 is the overall heat transfer coefficient calculated by QUB experiment (Figure 1, equation ( 5)), and 𝐻 is the overall heat transfer coefficient between outdoor and indoor air, calculated from a model of the building.

The measurements needed for 𝐻 𝑄𝑈𝐵 are affected by errors. In order to estimate the error from measurements, let us consider a function of several variables 𝑓(𝑥 1 , … , 𝑥 𝑛 ). We want to evaluate the error of the function when the approximate values 𝑎 1 , … , 𝑎 𝑛 are known and the exact values 𝑥 1 , … , 𝑥 𝑛 are unknown.

If:

-the approximate measured values 𝑎 1 , … , 𝑎 𝑛 are statistically independent; -the approximate values have the errors normally distributed; -the errors 𝜎 1 = 𝑎 1 -𝑥 1 , … , 𝜎 𝑛 = 𝑎 𝑛 -𝑥 𝑛 are the standard deviations; -the errors are small; then the error (i.e. the standard deviation) of the function is [START_REF] Taylor | An Introduction to Error Analysis[END_REF][START_REF]Evaluation of measurement data -Guide to the expression of uncertainty in measurement[END_REF]:

𝜎 𝑓 = 𝑓(𝑎 1 , … , 𝑎 𝑛 ) -𝑓(𝑥 1 , … , 𝑥 𝑛 ) = √ (𝜎 1 𝜕𝑓 𝜕𝑥 1 ) 2 + ⋯ + (𝜎 𝑛 𝜕𝑓 𝜕𝑥 𝑛 ) 2 (12)
The estimation of the error propagation in equation ( 5), needs estimation of measurement errors. The errors 𝜀 𝑃ℎ = 𝑃 ℎ -𝑃 ̅ ℎ ; 𝜀 𝑃𝑐 = 𝑃 𝑐 -𝑃 ̅ 𝑐 , 𝜀 𝛥𝑇ℎ = 𝛥𝑇 ℎ -𝛥𝑇 ̅ ℎ , 𝜀 𝛥𝑇𝑐 = 𝛥𝑇 𝑐 -𝛥𝑇 ̅ 𝑐 are obtained from the characteristics of the power devices and temperature sensors. The errors on slope, 𝜀 𝑎ℎ = 𝛼 ℎ -𝛼 ̅ ℎ ; 𝜀 𝑎𝑐 = 𝛼 𝑐 -𝛼 ̅ 𝑐 , are estimated from the correlation coefficient obtained when the linear regression is done.

Let us consider that 𝐻 is the exact value and 𝐻 𝑄𝑈𝐵 is the measured value, which depends on errors on 𝛼 ℎ , 𝛼 𝑐 , 𝑃 ℎ , 𝑃 𝑐 , 𝛥𝑇 ℎ and 𝛥𝑇 𝑐 (Figure 1).

We will note: The error coming from QUB method, equation ( 11), and the error coming from measurement, equation ( 13), are statistically independent. Then, the total error is:

𝜀 𝐻 = √ (𝐻 𝑄𝑈𝐵 -𝐻 ) 2 + 𝜀 𝐻𝑚 2 (14)
and the total relative error is

𝜀 𝐻% = √ (𝐻 𝑄𝑈𝐵 -𝐻 ) 2 + 𝜀 𝐻𝑚 2 𝐻 (15)
A reasonable estimation of the error coming from measurement is about 0.5 -1 °C [START_REF] Strachan | IEA Annex 58: Full-scale Empirical Validation of Detailed Thermal Simulation Programs[END_REF]. Therefore, modelling errors of 1 °C may be considered in the range of the measurement error, while modelling errors of 2 -3 °C are usual [START_REF] Strachan | Whole model empirical validation on a full-scale building[END_REF].

Global conductance of a building

The aim of QUB method is to measure the global heat conductance (or resistance) of a house. While the thermal transmittance of a building element is well defined (ISO 7345, 2018;ISO 13789, 2017), the areal heat transfer coefficient of the building, 𝐻′, given by equation ( 2), requires "a conventional definition of internal temperature, external temperature, area and the different contributions resulting in the heat flow rate. The heat flow rate may optionally include the contributions of heat transmission through building envelope, thermal bridges, ventilation, solar radiation, etc. The area may optionally be the envelope area, the floor area, etc. The temperature difference may optionally be a weighted temperature difference." (ISO 7345, 2018). Proposals and justifications for these conventional definitions are discussed hereafter.

10/26

One boundary temperature

The global conductance is well defined in steady-state if there is only one outdoor temperature, 𝑇 𝑜 , for the building (Figure 2.a):

𝐻 ≡ 𝑃 𝜃 𝑖 -𝑇 𝑜 (16)
This is a common interpretation of the overall conductance. Equation ( 16) is valid for steadystate. Since steady-state can never be fully achieved for a building, methods such as degreedays or degree-hours consider that the conductance, 𝐻 ≡ 𝑈𝐴, is constant and estimate the global transmittance by the integral in time [START_REF] Ghiaus | Potential for free-cooling by ventilation[END_REF]:

𝐻 = ∫ 𝑃𝑑𝑡 𝑡 𝑓𝑖𝑛𝑎𝑙 0 ∫ (𝜃 𝑖 -𝑇 𝑜 )𝑑𝑡 𝑡 𝑓𝑖𝑛𝑎𝑙 0 (17)
In dynamic model, the conductance is the reciprocal of the static gain, which is the total resistance (Figure 2c).

A simple thermal circuit shows that the outdoor temperature, 𝑇 𝑜 , and the power delivered to the building, 𝑃, are sources (Figure 2b). In a bloc diagram, the sources are inputs (Figure 2c, d &e). For a linear model, the blocs are transfer functions which can be represented generally by 𝐾𝐺(𝑠), where 𝐾 is the static gain and 𝐺(𝑠) is the dynamic gain. For strictly causal systems, the dynamic gain 𝐺(𝑠) tends to zero when time tends to infinity, 𝑠 → 0 ⟺ 𝑡 → ∞. For a building, the static gain of the outdoor temperature is 𝐾 𝑇 = 1 (from physical reasons, when 𝑃 = 0, if 𝑡 → ∞, then 𝜃 𝑖 → 𝑇 𝑜 ). The static gain of the power is the global resistance, 𝐾 𝑃 = 𝑅 ≡ 1/𝐻 (Figure 2d).

More boundary conditions

When there are more boundary conditions (Figure 2e), the definition is not so straightforward.

The output in steady-state is:

𝜃 𝑖 = ∑ 𝐾 𝑖 𝑖 𝑇 𝑖 + 𝐾 𝑝 𝑃 (18)
where 𝐾 𝑖 are the static gains of the inputs that are temperatures and 𝐾 𝑃 is the static gain of the power, 𝑃. This is true, whether all outdoor temperatures are the same, ∀𝑖, 𝑇 𝑖 = 𝑇 𝑜 , or not. Note that ∑ 𝐾 𝑖 𝑖 = 1 [START_REF] Ghiaus | Experimental estimation of building energy performance by robust regression[END_REF].

Therefore, by considering that the global conductance is a function of the supplied power and the temperature difference between indoor and outdoor, equation ( 16) will give errors in steady-state when there are multiple exterior boundary temperatures.

Multi-zone buildings

In the case of a multi-zone building, i.e. more indoor temperatures, the principle is to find the mean or equivalent thermal conductance. For a multi-zone building, we need the equivalent (mean) temperature, 𝜃 ̅ , and the equivalent (global) conductance, 𝐻 ≡ 𝑈𝐴, which produce the same heat loss over the whole surface area 𝐴:

𝑃 = 𝑈𝐴(𝜃 ̅ -𝑇 𝑜 ) (19)
as the temperature differences (𝜃 𝑖 -𝑇 𝑖 ) over the surface areas 𝐴 𝑖 :

𝑃 = ∑ 𝑈 𝑖 𝐴 𝑖 (𝜃 𝑖 -𝑇 𝑖 ) 𝑖 (20)
Therefore, for determining 𝑈, the mean temperature 𝜃 ̅ and the transmittances of the elements 𝑈 𝑖 𝐴 𝑖 are needed.

Mean temperature

In order to determine the total conductance, 𝐻 ≡ 𝑈𝐴, from equation ( 19), the equivalent mean temperature, 𝜃 ̅ , is needed. It is the temperature which gives the same accumulation of energy in the zone as if there were distributed temperatures:

∑ 𝑚 𝑖 𝑐𝜃 ̅ 𝑖 = ∑ 𝑚 𝑖 𝑐𝜃 𝑖 𝑖 (

Total conductance obtained from zone temperatures and conductance of elements

The total conductance 𝐻 ≡ 𝑈𝐴 can be defined in function of the zone temperatures and the conductance of the elements. From ( 19) and ( 20) it results

𝑈𝑆(𝜃 ̅ -𝑇 𝑜 ) = ∑ 𝑈 𝑖 𝐴 𝑖 (𝜃 𝑖 -𝑇 𝑜 ) 𝑖 (24)
where

𝐴 = ∑ 𝐴 𝑖 𝑖 ( 25
)
which has two unknowns: 𝜃 ̅ and 𝑈𝐴 ≡ 𝐾. Without reducing the generality, let us consider two zones for which,

𝑈 1 𝐴 1 (𝜃 1 -𝑇 𝑜 ) + 𝑈 2 𝐴 2 (𝜃 2 -𝑇 𝑜 ) = 𝑈(𝐴 1 + 𝐴 2 )(𝜃 ̅ -𝑇 𝑜 ) (26)
which yields the global conductance

𝑈𝜃 ̅ = 𝑈 1 𝐴 1 𝜃 1 + 𝑈 2 𝐴 2 𝜃 2 𝐴 1 + 𝐴 2 (27)
Therefore, from ( 16) and ( 27), the equivalent conductance for multi-zone is

𝐻 ≡ 𝑈𝐴 = ∑ 𝑈 𝑖 𝐴 𝑖 𝜃 𝑖 𝑖 • ∑ 𝑚 𝑖 𝑖 ∑ 𝑚 𝑖 𝜃 𝑖 𝑖 (28)
There are problems with equation ( 28):

• the conductances 𝑈 𝑖 𝐴 𝑖 are difficult to obtain in the case of coupled conductiveradiative heat transfer; • the temperatures 𝜃 𝑖 are needed.

Total conductance obtained from state-space model

Reformulated in terms of thermal circuits, the problem of finding the equivalent total conductance, eqs. ( 19) and (20), becomes: transform the complete thermal circuit characterized in steady-state by [START_REF] Ghiaus | Experimental estimation of building energy performance by robust regression[END_REF][START_REF] Strang | Computational Science and Engineering[END_REF]:

𝟎 = -𝐀 𝑇 𝐆𝐀𝛉 + 𝐀 𝐓 𝐆𝐛 + 𝐟 (29)
into a circuit with only one equivalent conductance (Figure 2b): 0 = -𝐺𝜃 ̅ + 𝐺𝑇 𝑜 + 𝑓 (30) so that the load is the same 𝑓 = ∑ 𝑃 𝑖 𝑖 = ∑ 𝐟.

13/26

If the temperature sources are zero, or if they are the same for all boundaries, i.e. 𝐛 = 𝟎, then the thermal circuit (29) can be written in matrix form as

𝐊𝛉 = 𝐟 (31)
where 𝐊 ≡ 𝐀 𝑇 𝐆𝐀. Without reducing the generality, let us consider a two-zone building (i.e. the size of 𝛉 is two). If there is only one outdoor temperature, then it becomes the reference temperature, 𝑇 𝑜 = 0. Therefore, the problem of finding the equivalent total conductance becomes: The advantage of using the state-space model for calculating the global conductance is evident when the heat transfer phenomena are coupled and the coefficients 𝐾 𝑖𝑗 from equation ( 32) do not have simple expressions. For a state space model,

transform
{ 𝐱̇= 𝐀𝐱 + 𝐁𝐮 𝐲 = 𝐂𝐱 + 𝐃𝐮 (38)
the static gain 𝐾 𝑖𝑗 of the transfer function between input 𝑗 and output 𝑖 can be obtained as the steady state response to a unitary step [START_REF] Rowell | Analysis and Design of Feedback Control Systems[END_REF])

𝐲 𝑠𝑠 = [-𝐂𝐀 -𝟏 𝐁 + 𝐃]𝐮 𝑆𝑆 (39)
where 𝐮 𝑆𝑆 is the vector of unitary step inputs 𝑢 1 (𝑡), … , 𝑢 𝑛 (𝑡):

𝐮 𝑆𝑆 = [ 𝑢 1 (𝑡) … 𝑢 𝑛 (𝑡) ] (40)
The static gains 𝐾 𝑖𝑗 between all outputs 𝑖 and the input 𝑗 is the steady-state value of the transfer function for input 𝑗 when 𝑢 𝑗 = 1 and 𝑢 𝑖≠𝑗 = 0. By using the solution (39) of the steady-state model and equation ( 22) for the mean indoor temperatures, the building heat transfer coefficient or the overall thermal conductance, 𝐻, is then

𝐻 = ∑ 𝑃 𝑖 𝒊 ∑ 𝑚 𝑖 𝜃 𝑖 𝒊 ∑ 𝑚 𝑖 𝒊 -𝑇 𝑜 (41)
where -𝑃 𝑖 is the power delivered to the air in zone 𝑖; -𝑚 𝑖 -mass of air in zone 𝑖; -𝜃 𝑖 -temperature of zone 𝑖 obtained from the output vector 𝐲 𝑠𝑠 given by equation ( 39); -𝑇 𝑜 -outdoor temperature.

Design of experiments for QUB method

The design of experiment is based on the calculation of the error of the QUB method, 𝜀 𝑄𝑈𝐵% = (𝐻 𝑄𝑈𝐵 -𝐻)/𝐻, as a function of the power 𝑃 ℎ and the time duration of the experiment, 𝑡 𝑄𝑈𝐵 (Figure 1). Examples are given for two buildings: 1. A test building with one-thermal zone (bungalow) considering only the temperatures of the outdoor air and the ground. 2. A multi-zone building with multiple exterior conditions.

Obtaining the response of the model to QUB test

The values used in the definition of 𝐻 𝑄𝑈𝐵 given by equations ( 5) or ( 9) are obtained from the step response of the state-space model (38) with initial conditions 𝐱(0) [START_REF] Rowell | Analysis and Design of Feedback Control Systems[END_REF]:

{ 𝐱(𝑡) = 𝑒 𝐀𝑡 𝐱(0) + 𝐀 -1 (𝑒 𝐀𝑡 -𝐈)𝐁𝐮 𝐲(𝑡) = 𝐂𝐱 + 𝐃𝐮 (42)
where the matrix exponential of square matrix 𝐀 is the state transition matrix 𝑒 𝐀𝑡 ≡ 𝚽(𝑡) = 𝐕𝑒 𝚲𝒕 𝐕 -1 , with 𝐕 the matrix of eigenvectors and 𝚲 the matrix of eigenvalues of matrix 𝐀. The initial state in equation ( 42) is obtained from the boundary conditions, 𝐮(0):

𝐱(0) = -𝐀 -1 𝐁𝐮(0) (43)
The value of the indoor temperature as a response to a step input is calculated from (42) -at time 𝑡 𝑄𝑈𝐵 by considering the initial state 𝐱(0) and -at time 2𝑡 𝑄𝑈𝐵 by considering the initial state as 𝐱(𝑡 ℎ ) ≡ 𝐱 𝑠1 .

Obtaining the response as a sum of exponentials

The step response given by equation ( 42) can be written as a function of eigenvectors 𝐕 and eigenvalues 𝚲 as

𝐲(𝑡) = 𝐂𝐌𝑒 𝚲𝑡 𝐌 -1 𝐱(0) + (𝐂𝐀 -1 𝐌)𝑒 𝚲𝑡 (𝐌 -1 𝐁𝐮) -𝐂𝐀 -1 𝐁𝐮 + 𝐃𝐮 (44)
Let us emphasize the coefficients of the exponentials by rewriting (44) as:

𝐲(𝑡) = 𝐂𝐌 𝐝𝐢𝐚𝐠(𝐌 -1 𝐱(0))𝑒 𝚲𝑡 + 𝐂𝐀 -1 𝐌 𝐝𝐢𝐚𝐠(𝐌 -1 𝐁𝐮)𝑒 𝚲𝑡 -(𝐂𝐀 -1 𝐁 + 𝐃)𝐮 (45)
The right-hand of equation ( 45) has three terms:

• the 1 st term is the exponential response influenced by the initial conditions;

• the 2 nd term is the exponential response influenced by the step input;

• the 3 rd term is the steady-state value (which does not depend on the initial conditions).

One-thermal zone house (bungalow)

Let us consider a single zone parallelepipedal bungalow (Figure 3). The walls are of insulated galvanized steel (surface of ceiling and ground floor 13.5 m 2 , surface of exterior walls 37.6 m 2 ). The windows are double glazed (surface 3.88 m 2 ). The model takes into account the infiltration rate (0.5 vol/h) and the radiative exchange between the walls (Figure 4). The spatial discretization of the walls can be done arbitrarily small; the one-slice discretization shown in Figure 4 is used only for convenience of representation.

6.

3.1 Only one boundary temperature (i.e. outdoor temperature)

The total resistance (𝑅 ≡ 1/𝐻) of the model of the bungalow is determined by

𝐲 𝑠𝑠 = [-𝐂𝐀 -𝟏 𝐁 + 𝐃]𝐮 = 𝑅 = 0.0191 °C/W (46)
for the unitary step input 𝐮 = [𝑇 𝑂 𝑇 𝑂 𝑇 𝑂 𝑇 𝑂 𝑇 𝑂 𝑃] 𝑇 = [0 0 0 0 0 1] 𝑇 . For our example, the largest time constant is 𝑇 = 26.4 • 10 3 s = 7.34 h (Figure 5d).

The results are given in Figure 5. Power larger than 943 W is needed for heating. Figure 5a shows the intrinsic error of the QUB method as a function of heating power, 𝑃 ℎ , and time, 𝑡 ℎ .

It can be seen that for a heating time larger than 4 h, the intrinsic error of the method is zero for any power (Figure 5a). However, measurements are affected by errors. Figure 5b shows the estimated error obtained by adding the measurement error given by equation ( 74). It can be noticed that important errors result by diminishing the heating time and the heating power. In order to obtain robust estimations, the heating time needs to be larger than 3 hours and the power larger than about 2000 W. The value of the heat transfer coefficient found by the QUB experiment is 𝐻 𝑄𝑈𝐵 = 61.2 W/K which can be compared with the initial supposed value 𝐻 𝑄𝑈𝐵 = 52.4 W/K obtained from the model used for the design of the experiment (equation ( 46)). The values of the heat transfer coefficient found by QUB experiment was compared with the value found by co-heating; as in other experiments, the differences between QUB and co-heating were about 5 -10 % [START_REF] Mangematin | Quick measurements of energy efficiency of buildings[END_REF][START_REF] Meulemans | An assessement of QUB/e method for fast in situ measurements of thermal performance of building fabrics in cold climates[END_REF]Alzetto, Farmer, Fitton, Huges, & Swan, 2018;Alzetto, Pandraud, Fitton, & Heusler, 2018;Alzetto, Gossard, & Pandraud, 2014) The fact that small measurement errors can be obtained for relatively short heating duration can be explained by the time constants and the coefficients of the exponential response of the model given by equation ( 45) and shown in Figure 5d. For the model of the bungalow, the time constants can be grouped in five classes: a) Very short time constants (< 0.1 h) with significant coefficients in equation ( 45). Their response is so quick that they act as a step response at the time scale of the heating and cooling. b) Small time constants (0.5 h … 1 h) but with very small coefficients in equation ( 45).

Due to the value of their coefficients, these exponentials do not have a significant contribution to the response. c) Medium time constants (1 … 2 h) with significant coefficients in equation (45). These are the time constants that have a significant contribution to the exponential response in the QUB method.

d) Medium time constants (2 … 3 h) but small coefficients in equation ( 45). Due to the values of their coefficients, these time constants have small influence on the response. e) Large time constants (> 6 h) with large coefficients in equation ( 45). These time constants have little influence on the response due to the relatively short heating and cooling time used in the experiment.

Therefore, the response of the system may be considered a pure exponential plus a constant when the durations of heating and cooling are larger than 3 h (4 times the time constant of 0.8 h). This observation shows an important feature of the QUB method: it allows measuring the heat transfer coefficient in a short time but this time needs to be long enough in order to obtain the steady state for the small time-constants which have significant coefficients and short enough to have small values of the exponential corresponding to large time constants.

Figure 5d shows another interesting feature of the thermal models of buildings. If the model is obtained by a fine discretisation, the time constants are grouped: many small-time constants and a few, rather clearly separated, larger time constants. Figure 5c shows the time response for the conditions of the experiment. Two features can be noticed: the fast response due to the short time constants and the exponential response due to the medium time constants associated with significant coefficients.

Figure 6 Error evaluation of heat transfer coefficient for the bungalow when: a) the global conductance, 𝐾 𝑃 , is considered; b) the experimentally measurable global conductance, 𝐻, is considered.

More boundary temperatures for the bungalow

Usually, the residential buildings have a non-adiabatic ground floor. Inherently, the QUB method, as any other measurement method, will do an error if the ground temperature is not taken into account. Figure 6a shows the intrinsic error of the QUB method when the value of the conductance is the correct one, i.e. representing the static gain of power calculated from the physical values; in this case, the value of 𝐻 does not depend on the boundary temperatures. Figure 6b shows the errors when the total conductance is obtained with equations ( 16) and/or (17) (i.e. by considering only the outdoor air temperature), which is a common interpretation of the total conductance; this time, the value of 𝐻 depends of the boundary temperatures. Figure 6 shows that, by increasing the power, the error is reduced, i.e. the influence of the ground floor in the response is reduced.

Case of a typical house

Finally, let us consider a typical two-floor house with the ground floor area of 93.3 m 2 , the total floor area 186.6 m 2 and total volume of 502 m 3 (Figure 7). There are triple glazed windows with aluminum frames; external walls are made of silicate bricks of 18 cm with external mineral wool insulation. The ceiling has an external insulation made of 54 cm of mineral wool. There is no crawl space below the ground floor; there are 31 cm of expanded polystyrene and 25 cm of reinforced concrete. The value of the air leakage rate measured by blower door is 𝑛 50 = 0.34 vol/h at 50 Pa.

It can be seen from Figure 8a that, if the outdoor temperature is 19.5°C, the optimal design of the experiment is for a power of about 2500 W and for a heating time of 6 h. The error obtained in this case is about 7 % (coming from the ground floor) and the equivalent indoor temperature is 20.3°C. Figure 8a shows the valley in which the optimum is located. The experiment is very sensitive to power variation and less sensitive to time duration. Figure 8b shows that the response is quasi-instantaneous at the beginning (corresponding to the small time-constants), followed by an exponential form (corresponding to the medium time constants). This is explained by the distribution of the 149 time-constants and the amplitude of their corresponding coefficients obtained by using equation (45) (Figure 8c). Similarly to the case of the bungalow, the time constants can be grouped in five classes: a) Very short time constants (0.5 • 10 -1 h … 2 • 10 -1 h) with significant coefficients in equation ( 45). Their response is so quick that they act almost as a step response at the time scale of the heating and cooling. b) Small time constants (2 • 10 -1 h … 10 0 h) but with very small coefficients in equation ( 45).

Due to the value of the eigenvectors, they do not have a significant contribution to the response. c) Medium time constants (around 10 0 h) with significant coefficients in equation ( 45). These are the time constants that have a significant contribution to the exponential response in the QUB method. d) Medium time constants (around 10 1 h) but small coefficients in equation ( 45). Due to the values of their coefficients, these time constants have small influence on the response. e) Large time constants (> 10 1 h) and large coefficients in equation ( 45). These time constants have little influence on the response due to the relatively short heating and cooling time.

Therefore, at the time scale of the QUB experiment, the system with these time constants gives an exponential response (Figure 8b). The distribution of time constants and of their corresponding coefficients justifies the use in QUB experiment of a short time as compared with the largest time constants of the building.

Conclusions

A method to design QUB experiments was developed and experimentally tested on two real buildings: a bungalow and a typical house. The method gives the error of the overall thermal transmittance, 𝐻, as a function of the heating power and the heating duration. Initial conditions, outdoor temperature and power during cooling may be introduced in the design of the experiment. The method is valid if the outdoor conditions change insignificantly during the experiment. Therefore, the experiment is conducted during the night time.

The relatively small-time duration of QUB experiment is explained by the distribution of the time constants. There are a lot of time constants that have small values, a few that have medium values and significant coefficients, and a few with large values and significant coefficients. The QUB method takes advantage of this distribution of the time constants by designing an experiment at the time scale of the medium time constants (about a few hours). Therefore, a good design of the experiment of QUB method allows to obtain measurement errors less than 10 % for real in-situ houses with experimental durations shorter than 12 h.

This method allows the experimenter to select the best input parameters (heating power and duration) based on a priori knowledge of the building envelope to be measured and on the weather forecast. It could also be used as a decision tool to make or not a measurement depending on the predicted error that would be committed as a function of building constructive characteristics and weather. This method could help improving the data analysis by selecting best period for fitting the thermal model on the measured data.

Annexes

Annex 1: Slope of the response curve of QUB test in heating and cooling In heating, it is supposed that the space is submitted to two inputs: the outdoor temperature, 𝑇 𝑜 , and the power, 𝑃 ℎ . In order to emphasize this, by considering Δ𝑇 ≡ 𝜃 𝑖 -𝑇 𝑜 , let's write equation ( 4) as 𝐶𝜃 𝑖 ̇= -𝐻 𝑄𝑈𝐵 𝜃 𝑖 + 𝑃 ℎ + 𝐻 𝑄𝑈𝐵 𝑇 𝑜 .

(47)

Equation ( 47) is linear; it can be solved by applying the superposition of two equations:

𝐶𝜃 ̇𝑖 = -𝐻 𝑄𝑈𝐵 𝜃 𝑖 + 𝑃 ℎ (48)
and 𝐶𝜃 ̇𝑖 = -𝐻 𝑄𝑈𝐵 𝜃 𝑖 + 𝐻 𝑄𝑈𝐵 𝑇 0 ℎ .

(49)

The general solution of equation ( 48) is: is obtained. The solution of equation ( 48) is:

𝜃 𝑖 =
𝜃 𝑖 = 𝑃 ℎ 𝐻 𝑄𝑈𝐵 (1 -exp (- 𝑡 𝜏 )) (51)
The solution of equation ( 49) is:

𝐶𝜃 𝑖 ̇= -𝐻 𝑄𝑈𝐵 ′ (𝜃 𝑖 -𝑇 0 ℎ ) ; 𝛥𝑇 ≡ 𝜃 𝑖 -𝑇 0 ℎ , 𝛥𝑇 = 𝑐 1 exp(-𝑡/𝜏) ; 𝛥𝑇(0) = 𝛥𝑇 0 , 𝛥𝑇 = 𝛥𝑇 0 exp(-𝑡/𝜏) . which is the slope in origin of the sum of two exponential responses given by equation ( 53).

In QUB method, the tangent to the exponential response at a time 𝑡 = 𝑡 0 > 0 is used. The equation of the tangent at 𝑡 = 𝑡 0 is: with the slope 𝛼 ℎ = 𝑒 -𝑡 0 𝜏 (𝑃/𝐶 -Δ𝑇 0 ℎ 𝐻 𝑄𝑈𝐵 /𝐶)

Δ𝑇 ℎ = (1 -𝑒
Assuming that the outdoor temperature, 𝑇 𝑜 , is constant and defining: 𝛥𝑇 𝑐 = 𝜃 𝑖 -𝑇 𝑜 the difference between the indoor and outdoor temperatures; The tangent in the origin of the exponential given by eq. ( 60) has the equation

Δ𝑇

  outside temperature difference, Δ𝑇 ≡ 𝜃 𝑖 -𝑇 𝑜 ; 𝑃 -power delivered by electric heaters to the indoor air; 𝐶 -effective heat capacity; 𝐻 𝑄𝑈𝐵 -heat transfer coefficient.

Figure 1

 1 Figure 1 Principle of QUB method (Alzetto, Gossard, & Pandraud, 2014)

Figure 2

 2 Figure 2 Definition of the total conductance: a) on a resistance; b) on a 1 st order thermal circuit model; c) on a bloc diagram of a 1 st order model (which corresponds to the thermal circuit from b); d) on a bloc diagram of a model with only one temperature as an input; e) on a bloc diagram of a model with more temperatures as inputs.

  Simplifying by the specific mass and by the height of the zones,

Figure 3

 3 Figure 3 General view of the bungalow

Figure 5

 5 Figure 5 Design of experiment for a bungalow: a) relative intrinsic error 𝜀 𝑄𝑈𝐵% as a function of power and experiment duration; b) sum of intrinsic and measurement error c) time response for 𝑡 𝑄𝑈𝐵 = 3 h and 𝑃 ℎ = 2000 W; d) Time constants (stars) and amplitudes of the corresponding exponentials (circles) of the model.

Figure 7 :

 7 Figure 7: Thermal zones of the house considered in the experiments: (a) ground floor (b) first floor.

Figure 8

 8 Figure 8 Design of experiment for a house: a) relative error 𝜀 𝑄𝑈𝐵% as a function of power and duration of experiment; b) time response for 𝑃 ℎ = 2500 W and 𝑡 𝑄𝑈𝐵 = 6 h; c) Time constants (stars) and amplitudes of the corresponding exponentials (circles) of the model.

  𝑖 | 𝑡=0 -𝑇 𝑜 the difference at initial time, 𝑡 = 0, equation (4) becomes 𝐶 𝑑 𝑑𝑡 Δ𝑇 𝑐 = -𝐻 𝑄𝑈𝐵 Δ𝑇 c + 𝑃 𝑐 (58) with initial condition 𝛥𝑇 𝑐 = 𝛥𝑇 0 𝑐 . The general solution is Δ𝑇 𝑐 = 𝑐 1 exp (

  𝛼 ℎ , 𝛼 𝑐 , 𝑃 ℎ , 𝑃 𝑐 , 𝛥𝑇 ℎ and 𝛥𝑇 𝑐 the measured values that are affected by absolute errors which are relatively small as compared to the true values.-𝛼 ̅ ℎ , 𝛼 ̅ 𝑐 , 𝑃 ̅ ℎ , 𝑃 ̅ 𝑐 , 𝛥𝑇 ̅ ℎ and 𝛥𝑇 ̅ 𝑐 the true (exact) values.-𝜀 𝑎ℎ = 𝛼 ℎ -𝛼 ̅ ℎ ; 𝜀 𝑎𝑐 = 𝛼 𝑐 -𝛼 ̅ 𝑐 ; 𝜀 𝑃ℎ = 𝑃 ℎ -𝑃 ̅ ℎ ; 𝜀 𝑃𝑐 = 𝑃 𝑐 -𝑃 ̅ 𝑐 , 𝜀 𝛥𝑇ℎ = 𝛥𝑇 ℎ -𝛥𝑇 ̅ ℎ , 𝜀 𝛥𝑇𝑐 = 𝛥𝑇 𝑐 -𝛥𝑇 ̅ 𝑐 the absolute errors (i.e. the measured minus the exact value) of these variables, respectively. If similar devices and instruments are used, it is reasonable to consider that 𝜀 𝑎ℎ = 𝜀 𝑎𝑐 = 𝜀 𝑎 , that 𝜀 𝑃ℎ = 𝜀 𝑃𝑐 = 𝜀 𝑃 and that 𝜀 𝛥𝑇ℎ = 𝜀 𝛥𝑇𝑐 = 𝜀 𝛥𝑇 .The probable error of 𝐻 𝑄𝑈𝐵 obtained from measurements is (see Annex 2):

	-𝜀 𝐻𝑚 = [(𝜀 𝛼	𝜕𝐻 𝑄𝑈𝐵 𝜕𝛼 ℎ + (𝜀 𝛥𝑇 ) 2 + (𝜀 𝛼 𝜕𝐻 𝑄𝑈𝐵 𝜕𝐻 𝑄𝑈𝐵 𝜕𝛼 𝑐 𝜕𝛥𝑇 ℎ ) 2 + (𝜀 𝛥𝑇 ) 2 + (𝜀 𝑃 𝜕𝛥𝑇 𝑐 𝜕𝐻 𝑄𝑈𝐵 𝜕𝑃 ℎ ) ] ) 𝜕𝐻 𝑄𝑈𝐵 2 1/2 2	+ (𝜀 𝑃	𝜕𝐻 𝑄𝑈𝐵 𝜕𝑃 𝑐	)	2	(13)
							9/26

  Annex 2: Probable error from measurementsThe true value of 𝐻 𝑄𝑈𝐵 is𝐻 𝑄𝑈𝐵 (𝛼 ̅ ℎ , 𝛼 ̅ 𝑐 , 𝑃 ̅ ℎ , 𝑃 ̅ 𝑐 , 𝛥𝑇 ̅ ℎ , 𝛥𝑇 ̅ 𝑐 ) = 𝐻 𝑄𝑈𝐵 (𝛼 ℎ -𝜀 𝛼 , 𝛼 𝑐 -𝜀 𝛼 , 𝑃 ℎ -𝜀 𝑃 , 𝑃 𝑐 -𝜀 𝑃 , 𝛥𝑇 ℎ -𝜀 𝛥𝑇 , 𝛥𝑇 𝑐 -𝜀 𝛥𝑇 )By developing the right-hand side in a series expression, we obtain where 𝐻 ̅ 𝑄𝑈𝐵 ≡ 𝐻 𝑄𝑈𝐵 (𝛼 ̅ ℎ , 𝛼 ̅ 𝑐 , 𝑃 ̅ ℎ , 𝑃 ̅ 𝑐 , 𝛥𝑇 ̅ ℎ , 𝛥𝑇 ̅ 𝑐 ) and 𝐻 𝑄𝑈𝐵 ≡ 𝐻 𝑄𝑈𝐵 (𝛼 ℎ , 𝛼 𝑐 , 𝑃 ℎ , 𝑃 𝑐 , 𝛥𝑇 ℎ , 𝛥𝑇 𝑐 )

	with the slope
	α c =	𝑃 𝑐 𝐻 𝑄𝑈𝐵	-Δ𝑇 0 𝑐 /𝜏 =	𝑃 𝑐 𝐻 𝑄𝑈𝐵	-	Δ𝑇 0 𝑐 𝐻 𝑄𝑈𝐵 𝐶	(62)
	The equation of the tangent at 𝑡 = 𝑡 0 is
	Δ𝑇 𝑐 = (1 -𝑒 -𝑡 0 𝜏 )	𝑃 𝑐 𝐻 𝑄𝑈𝐵	+ 𝑒 -𝑡 0 𝜏 [Δ𝑇 0 𝑐 + (	𝑃 𝑐 𝐶	-	𝑐 Δ𝑇 0 𝜏	) (𝑡 -𝑡 0 )]	(63)
	with the slope
	α c = 𝑒 -𝑡 0	𝐻 𝑄𝑈𝐵 𝐶 ( 𝑃 𝑐 𝐶	-	Δ𝑇 0 𝑐 𝐻 𝑄𝑈𝐵 𝐶	)	(64)
														(65)
	𝐻 ̅ 𝑄𝑈𝐵 =
	𝐻 𝑄𝑈𝐵 -𝜀 𝛼	𝜕𝐻 𝑄𝑈𝐵 𝜕𝛼 ℎ	-𝜀 𝛼	𝜕𝐻 𝑄𝑈𝐵 𝜕𝛼 𝑐	-𝜀 𝑃	𝜕𝐻 𝑄𝑈𝐵 𝜕𝑃 ℎ	-𝜀 𝑃	𝜕𝐻 𝑄𝑈𝐵 𝜕𝑃 𝑐	-𝜀 𝛥𝑇	𝜕𝐻 𝑄𝑈𝐵 𝜕𝛥𝑇 ℎ	-𝜀 𝛥𝑇	𝜕𝐻 𝑄𝑈𝐵 𝜕𝛥𝑇 𝑐	(66)
	and the partial derivatives are
	𝜕𝐻 𝑄𝑈𝐵 𝜕𝛼 ℎ	=	𝛥𝑇 𝑐 (𝑃 ℎ 𝛼 𝑐 -𝑃 𝑐 𝛼 ℎ ) 𝜎 2	-	𝑃 𝑐 𝜎	(67)
	𝜕𝐻 𝑄𝑈𝐵 𝜕𝛼 𝑐	= -	𝛥𝑇 ℎ (𝑃 ℎ 𝛼 𝑐 -𝑃 𝑐 𝛼 ℎ ) 𝜎 2	+	𝑃 ℎ 𝜎	(68)
	𝜕𝐻 𝑄𝑈𝐵 𝜕𝑃 ℎ	=	𝜎 𝛼 𝑐								(69)
	𝑐 = Δ𝑇 0 𝑐 + ( 𝜕𝐻 𝑄𝑈𝐵 𝜕𝑃 𝑐 = -𝛼 ℎ 𝜎	𝑃 𝑐 𝐻 𝑄𝑈𝐵	-	Δ𝑇 0 𝑐 𝜏	) 𝑡	(61) (70)

 Bouchié, R., Alzetto, F., Brun, A., Weeks, C., Preece, M., Muhammad, A., & Sisinni, M. (2015).

Methodologies for the Assessment of Intrinsic Energy Performance of Buildings Envelope. PERFORMER Project. Brisken, W., & Mitalas, G. (1956). Heat Load Calcula-tions by Thermal Respons. ASHRAE Transaction, 62, 391. Danov, S., Carbonell, J., Cipriano, J., & Marti-Herrero, J. (2013). Approaches to evaluate building energy performance from daily consumption data considering dynamic and solar gain effects. Energy and Buildings,57,[110][111][112][113][114][115][116][117][118]


Romuald Coupan

Frédéric Plantier

Jean-Philippe Torré

Christophe Dicharry

Pascale Senechal

Fabrice Guerton

Peter Moonen

Abdel Khoukh

Sid Ahmed Kessas

Mehrdji Hemati

Creating innovative composite materials to enhance the kinetics of CO 2 capture by hydroquinone clathrates

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

with a low molecular weight, such as carbon dioxide (CO 2 ), methane (CH 4 ) and nitrogen (N 2 ), can be trapped as guest mole cules in a clathrate structure. A variety of host molecules, such as water, phenol and hydroquinone (HQ), have been studied in recent years and are described in literature [2]. Scientists have been fasci nated by gas clathrates for more than a century, but these solid compounds are now the focus of even greater attention owing to their possible practical applications such as CO 2 capture and sequestration. This new interest arises from their capacity to store large quantities of gases, more or less selectively depending on the clathrates considered and the gas mixture used.

To begin with, the most studied clathrate compounds are ''gas hydrates" [3], which generally form at high pressure (a few MPa) and low temperature (a few K above 273 K) with water and gas. For the past 15 years or more, they have been the focus of many studies [4], which included adding chemical promoters or using different multiphase contactors [5 7] to overcome the many criti cal limitations that prevent hydrate based gas separation pro cesses from being deployed, from pilot to industrial scale. However, the use of gas hydrates to separate CO 2 from gas mix tures such as CO 2 /CH 4 and CO 2 /N 2 is still hindered by the following main obstacles: (i) the overall cost of the process is not competitive compared to others such as reactive absorption [START_REF] Herri | CO2 capture by using hydrates: From experimental evidences to economic impossibility as a post combustion process[END_REF], except when the raw gas is already pressurized (as the gas compression costs are prohibitive in most cases); (ii) the difficulties in handling, con tacting and transporting a multiphase mixture composed of at least three phases that have to coexist (water or water/oil mixture, gases and solids) for the reaction to take place; (iii) the necessity to maintain a relatively low temperature in the apparatus for the gas hydrates to form (generally slightly above 273 K); and (iv) the low rate of hydrate formation and insufficient selectivity of gas hydrates toward CO 2 in certain gas mixtures, particularly CO 2 / CH 4 [START_REF] Ricaurte | CO 2 Removal from a CO 2 -CH 4 gas mixture by clathrate hydrate formation using THF and SDS as water-soluble hydrate promoters[END_REF][START_REF] Ricaurte | Combination of surfactants and organic compounds for boosting CO 2 separation from natural gas by clathrate hydrate formation[END_REF].

Secondly, HQ is a solid hydroxy aromatic compound recently identified in literature as a potential innovative clathrate candidate for gas storage and separation [11 13]. Interestingly, for different gas mixtures such as CO 2 /N 2 , CO 2 /H 2 and CO 2 /CH 4 [START_REF] Lee | Preferential occupation of CO 2 molecules in hydroquinone clathrates formed from CO 2 /N 2 gas mixtures[END_REF][START_REF] Lee | Selective CO 2 trapping in guestfree hydroquinone clathrate prepared by gas phase synthesis[END_REF], the enclathration reaction with HQ seems to be highly selective to CO 2 . The thermodynamic pressure temperature conditions (i.e. the clathrate phase equilibria) in which some HQ clathrates form and dissociate are mentioned in literature [16 19] and the equilib rium data and detailed characterization of CO 2 HQ clathrates were recently published [20 21]. But although HQ clathrates seem very interesting for practical applications such as CO 2 separation by direct gas/solid reaction (i.e. without using any solvent), the gas clathrate formation kinetics need to be improved before these compounds can be used on an industrial scale. Previous studies involving tests carried out using powdered HQ have already attempted to do so [START_REF] Lee | Preferential occupation of CO 2 molecules in hydroquinone clathrates formed from CO 2 /N 2 gas mixtures[END_REF][START_REF] Lee | Selective CO 2 trapping in guestfree hydroquinone clathrate prepared by gas phase synthesis[END_REF][START_REF] Allison | Clathration by phenol and quinol. Part. II Kinetics[END_REF]. But very fine HQ particles would be difficult to handle in an industrial gas solid contactor (such as a fixed bed reactor) for several reasons: (i) the gas circulating through the HQ powder causes a significant pressure drop that in some cases could be excessive, (ii) the very fine particles can be blown into the process equipment, and (iii) operators must apply stringent safety measures when handling large quantities of a fine pulverulent material such as HQ during loading and unloading of the reactor. Indeed, reactivity is limited by the contact area between the gas and the solid HQ, and depends on the specific area of the reactive medium used. Consequently, in an attempt to enhance the gas solid enclathration reaction, we developed a com posite material made of porous particles impregnated with HQ that can be used as a reactive medium in a clathrate based gas separa tion process.

There are two procedures for impregnation: wet impregnation (abbreviated as WI below) and dry impregnation (abbreviated as DI below) in a fluidized bed. Creation of a composite material by WI consists of a dip coating method involving several steps: impregnation, filtration and drying [START_REF] Jiang | A review of wet impregnation -an alternative method for the fabrication of high performance and nano-structured electrodes of solid oxide fuel cells[END_REF]. DI in a fluidized bed was previously developed in the Laboratoire de Génie Chimique (chemi cal engineering laboratory) of Toulouse and patented by the INPT and the CNRS in 2001 [START_REF] Hemati | Fabrication de catalyseur d'un nouveau type par imprégnation en lit fluidise[END_REF]. Unlike WI, creation of a composite material by DI in a fluidized bed is a ''one pot process" in which the stages of impregnation, filtration and drying are carried out in just one apparatus. This technique consists in spraying a solu tion containing the organic precursor onto a warm fluidized bed of porous particles chosen as the support material for the precur sor. This condensed technique was initially used for depositing metal nanoparticles on porous supports (such as silica gel, alumina and activated carbon) and has shown that in operating conditions: (i) the deposition efficiency is close to 100%, (ii) the deposit is quasi uniform, and (iii) the size of the supported metal nanoparticles is controlled by the pore diameter [ 25 29]. In the present work, the composite materials synthesized by the two procedures (i.e. WI and DI) are evaluated as reactive media for CO 2 capture. Both the gas storage capacity of the reactive media and data on the CO 2 cap ture kinetics were measured and compared to the results obtained with pure powdered HQ. All experiments were performed at 3.0 MPa and 323 K by means of a gravimetric method using a mag netic suspension balance. These experimental conditions were chosen for two main reasons: (i) CO 2 HQ clathrates can form under relatively high temperature (e.g. 323 K) and pressure [START_REF] Coupan | Experimental determination of phase equilibria and occupancies for CO 2 , CH 4 , and N 2 hydroquinone clathrates[END_REF] condi tions, which could be economically very advantageous when the raw gas is already hot and pressurized, the case for example of a production gas containing CO 2 /CH 4 ; (ii) CO 2 HQ clathrates may reach their maximum capacity when the formation pressure is close to 3.0 MPa [START_REF] Peyronel | On some new clathrates of hydroquinone[END_REF], justifying the operating pressure chosen for these tests.

Experimental section

Materials

The organic precursor used in this study is HQ (purity of 99.5 mol%) provided by Acros Organics. The solvent used in WI and DI is absolute ethanol (purity greater than 99 mol%). This choice was motivated by the high solubility of HQ in this solvent (e.g. 0.2063 mol% at 308.2 K) [START_REF] Li | Solubility of hydroquinone in different solvents from 276[END_REF] and its low boiling temperature of 351 ± 2 K (calculated mean values obtained from the NIST) that facilitates the drying step. The N 2 , CO 2 and helium (He) used for the experiments (minimum mole fraction purity of 99.995%) are pur chased from Linde Gas SA. The porous supports used are silica par ticles provided by SiliCycle. The experiments are conducted with two types of silica particles: high grade irregular silica particles (SiliaFlash Ò , referred to as SF below) and high purity analytical grade spherical silica particles (SiliaSphere Ò , referred to as SS below) with tight particle size distribution. Information concerning native SF and SS particles is given in the Supplementary Material.

Experimental set up and methods

Prior to impregnation, a saturated solution of HQ in ethanol was prepared at 308 K based on the solubility data and a correlation proposed by Li and coworkers [START_REF] Li | Solubility of hydroquinone in different solvents from 276[END_REF]. Before each experiment, the solid particles were dried for 24 h at 423 K in a muffle furnace (from Nabertherm).

Wet impregnation

WI is performed in hermetic glass vessels. An oven (model 800, from Memmert) is used to maintain the temperature at a constant value with an accuracy of ±0.5 K. A vacuum filtration system com posed of a Buchner funnel connected to a vacuum pump is used for same way as for the global capture capacity and reaction kinetics by using this type of composite material.

Acknowledgments

We would like to acknowledge the entire work group involved in the ORCHIDS project and the Carnot Institute ISIFoR (Institute for the Sustainable engIneering of Fossil Resources). We extend our thanks to the staff of the Service Analyse et Procédés of the University of Toulouse for their help. We also wish to thank Mathieu Molinier, Victor Cabrol and Sofiane Bekhti for their tech nical assistance and the Total E&P, ''Gas Solutions" department for its financial backing. We would also like to thank Total for pro viding UMS 3360 DMEX with the Zeiss Xradia Versa 510 used to carry out the tomographic acquisitions reported in this article.

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cej.2017.05.038.

SUPPLEMENTARY MATERIAL 1. INFORMATION CONCERNING THE NATIVE SUPPORT MATERIALS

The native support materials -SiliaFlash® (SF) and SiliaSphere® (SS) silica particles -used in this study underwent particle size characterization and gas porosimetry in order to experimentally determine their textural parameters. Our experimental data and those provided by the supplier are given in Table SI1 and Table SI2 respectively. As the batch of SF particles presented a large size dispersion (i.e. from 200 to 1,000 µm), it had to be separated into two fractions to allow homogenous fluidization: a low fraction of particles from 200 to 500 µm (called L-SF below) and a high fraction of SF particles from 500 to 1,000 µm (referred to as H-SF below).

Note that the SS particles were not evaluated by gas porosimetry, as the PSD was out of the measurement range of the equipment. For the native SF particles, the PSD analysis revealed mesopores of about 38 ± 5 nm, in good agreement with supplier data (pore diameter of 33 nm provided). The operating conditions applied for developing the composite materials by Dry Impregnation (DI), particularly the minimum fluidization velocities of the particles obtained at room temperature, are listed in Table SI3.

Table SI3. Operating conditions applied for the DI of SiliaFlash® (SF) and SiliaSphere® (SS). detector from Oxford Instruments.

Support material H-SF L-SF SS

TECHNICAL DETAILS OF ANALYTICAL TECHNIQUES

The mean particle diameter and particle size distribution were evaluated by sieving and by using a laser diffraction particle size analyzer from Malvern Instruments (model Mastersizer 2000). The laser diffraction particle size analyzer works for particles of about 0.02 µm to 2,000 µm.

The morphology of the particles was obtained using a field emission gun scanning electron microscope (SEM-FEG) from JEOL (model JSM 7100F TTLS). This SEM is equipped with an energy dispersive X-ray (EDX) spectrometer. The EDX technique detects X-rays emitted by the sample while it is being bombarded by an electron beam and is able to characterize the elemental composition of the surface analyzed.

The precursor content, defined as the mass of HQ per unit mass of support, was determined by thermogravimetric analysis (TGA) by measuring the mass loss rate of the sample subjected to a temperature ramp from ambient temperature to 773 K in an oxidizing atmosphere (air).

The thermal signature of the impregnated compound was verified by differential scanning calorimetry (DSC). A Model Q600 analyzer from TA Instruments was used to perform both TGA and DSC.

The attrition resistance was determined by a simple attrition test. A sample of the composite material was sieved for 3 hours and its precursor content was measured by TGA. The value obtained was compared with the precursor content of the original sample.

The textural parameters -specific area, total pore volume, and pore size distribution (PSD)were estimated from N 2 adsorption-desorption isotherms at 77 K using an automatic Micromeritics TriStar II 3020 system. This system can be used to evaluate mesoporous materials (i.e. materials with a pore size of 2 to 50 nm) with a specific area reaching several thousands of m 2 /g. Before the measurements were taken, samples were purified at 308 K and 10 Pa for 24 hours. Specific areas were determined based on the Brunauer-Emmett-Teller (BET) method [Brunauer, S., Emmett, P.H., Teller, E., 1938. Adsorption of gases in multimolecular layers. J. Am. Chem. Soc. 60, 309-319]. Both the total pore volume and the pore size distribution were calculated using the Barrett-Joyner-Halenda (BJH) method [Barrett, E. P., Joyner, L.G., Halenda, P.P., 1951. The determination of pore volume and area distributions in porous substances. I. Computations from nitrogen isotherms. J. Am. Chem.

Soc. 73, 373-380] applied to the desorption branch of the N 2 isotherm.

Structural data were obtained by cross-polarization / magic angle spinning (CP/MAS) solidstate 13 C NMR experiments using a Bruker Avance 400 spectrometer. The sample was placed in a zirconia rotor of 7 mm and the spectra were collected at ambient temperature at a resonance frequency of 100.61 MHz and a magic angle spinning rate of 7 kHz. The contact time and pulse delay were 3.8 ms and 5 s respectively. Chemical shifts were referenced to tetramethylsilane at 0 ppm and carbonyl signal of glycine at 176 ppm as external reference.

The three-dimensional microstructure was obtained by X-ray tomography using a Zeiss Xradia Versa 510. Samples were subjected to the ambient test conditions (28°C) for at least


José Horta
email: jose.horta@telecom-paristech.fr

Daniel Kofman
email: daniel.kofman@telecom-paristech.fr

Eitan Altman
email: eitan.altman@sophia.inria.fr

David Menga
email: david.menga@edf.fr

Mathieu Caujolle
email: mathieu.caujolle@edf.fr

Real-time enforcement of local energy market transactions respecting distribution grid constraints

Future electricity distribution grids will host a considerable share of the renewable energy sources needed for enforcing the energy transition. Demand side management mechanisms play a key role in the integration of such renewable energy resources by exploiting the flexibility of elastic loads, generation or electricity storage technologies. In particular, local energy markets enable households to exchange energy with each other while increasing the amount of renewable energy that is consumed locally. Nevertheless, as most ex-ante mechanisms, local market schedules rely on hour-ahead forecasts whose accuracy may be low. In this paper we cope with forecast errors by proposing a game theory approach to model the interactions among prosumers and distribution system operators for the control of electricity flows in real-time. The presented game has an aggregative equilibrium which can be attained in a semi-distributed manner, driving prosumers towards a final exchange of energy with the grid that benefits both households and operators, favoring the enforcement of prosumers' local market commitments while respecting the constraints defined by the operator. The proposed mechanism requires only oneto-all broadcast of price signals, which do not depend either on the amount of players or their local objective function and constraints, making the approach highly scalable. Its impact on distribution grid quality of supply was evaluated through load flow analysis and realistic load profiles, demonstrating the capacity of the mechanism ensure that voltage deviation and thermal limit constraints are respected.

I. INTRODUCTION

We consider future distribution grids in a context of massive deployment of Renewable Energy Sources (RES) and storage technologies. Maintaining quality of supply in such a context will require infrastructure reinforcements as usual, but also an active participation of Distribution System Operators (DSO) into the energy management of flexible Distributed Energy Resources 1 (DER). This can be achieved through Demand Side Management, in particular, by gaining access to flexibility from residential prosumers. A myriad of demand side management mechanisms have been proposed and applied for the scheduling of such flexible resources, from which the most relevant ones are those that provide incentives to 1 Flexible loads, controllable generation and storage resources. prosumers for their active participation in distribution grid energy management as selfish agents [START_REF] Molderink | Management and control of domestic smart grid technology[END_REF]- [START_REF] Horta | Novel market approach for locally balancing renewable energy production and flexible demand[END_REF]. In particular, we will consider local renewable energy markets that enable households to agree with each other one hour-ahead of time on the exchange of renewable energy. Such markets, as well as all energy management mechanisms that rely on ex-ante agreements (day-ahead, hour-ahead, etc.), imply agents subscribing commitments on future energy flows, as a result of market transactions. Nevertheless, such flows cannot be forecasted with high precision due to their dependency on human behavior and on exogenous variables controlled by nature (temperature, sunshine and wind). Thus, in order to cope with the uncertainty of demand and local renewable production, a real-time control of energy resources must provide the means to drive households towards a final exchange with the grid that benefits both prosumers and DSOs.

In this work we address the control of electricity flows on a residential low voltage distribution grid in real-time. We propose a game theory model that enables households to decide which electricity flows to exchange with the distribution grid in order to minimize their costs. Such costs depend on the average strategy of all households and also on a penalty signal associated to a coupling constraint. This aggregative game enables households to optimize their own costs while reducing voltage rises/drops and current flows on the lines.

The main contributions of this work are the following:

• We rely on game theory to model the interactions between a DSO and prosumers for the control of their electricity flows in real-time. • The proposed game has a competitive aggregative equilibrium that can be attained through dynamics that only require one-to-all broadcast of price signals that do not depend on the amount of players. • To the best of our knowledge our work is the first to propose a real-time mechanism for the control of residential electricity flows in a context of massive deployment of renewable energy sources and storage, capable of enforcing both voltage deviations and thermal limit constraints.

The article is structured as follows. We start by presenting related work in Section II. Then we introduce the system under study in Section III. The model and the aggregative game are described in Sections IV and V respectively. In Section VI we present the simulations and results. Finally, in Section VII we conclude the article and provide some perspectives.

II. RELATED WORK

We focus on the literature related to real-time control of electricity flows on distribution grids. Several centralized approaches have been proposed for real-time energy management in distribution grid systems [START_REF] Cecati | Combined operations of renewable energy systems and responsive demand in a smart grid[END_REF]- [START_REF] Borghetti | Short-term scheduling and control of active distribution systems with high penetration of renewable resources[END_REF]. The difficulty in such approaches is that flexible DER are owned and controlled by households, which hold the local information needed for a centralized control and will release part of this information only if adequate incentives are provided. In [START_REF] Vytelingum | Agent-based micro-storage management for the smart grid[END_REF] and [START_REF] Ramchurn | Agent-based control for decentralised demand side management in the smart grid[END_REF] authors propose decentralized mechanisms with the objective to flatten the aggregated demand of a large set of households. In [START_REF] Li | Optimal demand response based on utility maximization in power networks[END_REF] they consider the goal of the electricity supplier is to maximize social welfare, which is achieved in a distributed fashion by households optimizing their own benefits.

Game theory has been applied mainly to day-ahead energy scheduling rather than to real-time control, and particularly to modeling the interactions between the electricity supplier and its clients, with the goal of minimizing energy costs rather than enforcing DSO constraints [START_REF] Mohsenian-Rad | Autonomous demand-side management based on game-theoretic energy consumption scheduling for the future smart grid[END_REF], [START_REF] Li | Optimal demand response based on utility maximization in power networks[END_REF], [START_REF] Samadi | Advanced demand side management for the future smart grid using mechanism design[END_REF]- [START_REF] Deng | Residential energy consumption scheduling: A coupled-constraint game approach[END_REF].

In [START_REF] Mohsenian-Rad | Autonomous demand-side management based on game-theoretic energy consumption scheduling for the future smart grid[END_REF] they propose a Stackelberg game model that allows electricity suppliers to define prices leading to an equilibrium that minimizes Peak to Average Ratio. In [START_REF] Soliman | Game-theoretic demand-side management with storage devices for the future smart grid[END_REF] they follow the same approach from [START_REF] Mohsenian-Rad | Autonomous demand-side management based on game-theoretic energy consumption scheduling for the future smart grid[END_REF], but with a strictly convex cost function. In [START_REF] Deng | Residential energy consumption scheduling: A coupled-constraint game approach[END_REF], an online version of a scheduling mechanism for flexible appliances is proposed, which copes with price prediction errors. The literature related to real-time control of electricity flows on distribution grids does not consider RES and storage resources or does not take particular care of voltage deviations or thermal constraints.

III. SYSTEM DESCRIPTION

The low voltage distribution grid under study has one or more three-phase four-conductors feeders to which households (with different distributions of resources: PV, batteries and flexible loads) are connected, all equipped with smart meters and with the possibility of controlling the flows they exchange with the electricity grid thanks to smart inverters [START_REF] Pignier | Energy management device and its associated method[END_REF]. Households agree on the exchange of energy with each other through a local energy market described in Section III-A. Then, a Real-Time Control mechanism (RTC) described in Section III-B copes with forecast errors by driving households towards a final exchange with the grid that benefits the prosumer and respects the DSO's quality of supply requirements. The RTC mechanism is combined with Dynamic Phase Switching (DPS) as described in Section III-C for balancing load across phases.

A. Local renewable energy market

We consider the scheduling mechanism proposed in [START_REF] Horta | Novel market approach for locally balancing renewable energy production and flexible demand[END_REF], where houses agree one hour-ahead of time, through a local market implemented over a blockchain-based transactive platform, on the amount of energy they exchange with each other on a 10 minutes basis. In addition to the exchange with neighbors, any remaining energy flow is contracted with the supplier. One hour ahead of every time slot, each house fixes the schedule of battery usage until the end of the day (moving window), imposing cycling and depth of discharge constraints aimed to extend the battery lifespan. The DSO provides agents with rewards for locally trading renewable energy when is best for the grid, which are attributed ex post to incentive agents to enforce their contracted quantities in real-time by absorbing all or part of their forecast errors with their flexibility budget2 .

B. Real-time control mechanism

Just before each time slot begins, prosumers update their forecast for the next 10 minutes and decide up to which point the error will be absorbed by their flexibility budget or by the electricity supplier. Deviating from the optimal flexibility schedule has an associated cost due to battery (or comfort) degradation, while deviating from the hour-ahead market committed quantities implies loosing part of the reward as well as the beneficial price. The costs of relying on extra flexibility, the final electricity price and the allocation of hourahead market rewards are described in Sections IV-C, IV-D and IV-E respectively. In addition, during the time slots of high excess of production or demand, households will be charged a penalty fixed by the DSO through the transactive platform, associated to a constraint on the average flow that neighbors exchange with the grid. Such a price signal allows the DSO to influence the decision of households towards a final exchange with the grid that reduces voltage rises/drops and current flows.

C. Dynamic phase allocation

We need to ensure the flows across phases are as evenly distributed as possible in order to improve RTC mechanism performance, which is not aimed at balancing load across phases. We assume the operator is capable of deploying solid state switches in a subset of households and of controlling dynamically their allocation to phases as described in [START_REF] Horta | Augmenting DER hosting capacity of distribution grids through local energy markets and dynamic phase switching[END_REF]. This functionality is complementary to the RTC mechanism.

IV. MODEL

We consider a set of competitive agents H = {1, 2, . . . , H} representing prosumers and a set N = {1, 2, . . . , N } of logical nodes on the feeder where the connections of several households are aggregated3 , referred to as points of common coupling (PCC). A household is connected to only one PCC. Households control the amount of energy exchanged with the grid by adapting the demand from batteries or flexible loads 4 . The goal of the operator is to influence households' decisions in order to reduce voltage deviations and current flows.

A. State of agents

The state of agents is composed by demand and production forecasts, by the quantities committed on the hour-ahead local market, and by the state of the battery at the beginning of the current time slot, determined by the State of Charge (SoC) and by the accumulated charge/discharge quantities. For each agent i ∈ H, the ex-ante mechanism scheduled the usage of the battery for the current time slot to be ŝi , in order to exchange xi = xi m + xi u with the grid, where xi m and xi u are the amounts traded with neighbors and with the supplier, at prices P m and P u (P f in the case of an injection) respectively, under a forecasted gap ri between production ĝi and consumption li according to the following balance equation

xi m + xi u + ŝi = li -ĝi = ri . (1)
A closer to real-time forecast or a Very Short Term Load Forecast [START_REF] Hsiao | Household electricity demand forecast based on context information and user daily schedule analysis from meter data[END_REF] enables each house to update its gap between consumption and production as r i = ri + ri , where ri = li -g i is the error on hour-ahead production and load forecasts.

B. Agents' strategies

The real-time control of each household aims to decide at which extent the error ri will be absorbed by the grid or by demand flexibility. Thus, each agent i decides on a strategy (x i , s i ), where s i ∈ R N represents the amount of energy charged (negative) or discharged (positive) to/from the batteries and x i ∈ R N is the total amount of energy to be exchanged with the grid on the current time slot 5 (includes the energy exchanged with neighbors and with the electricity supplier). The relationship between the battery flow s i and the grid flow x i is x i + s i = l ig i = r i , the same as in [START_REF] Molderink | Management and control of domestic smart grid technology[END_REF].

C. Batteries

The amount of energy stored in the battery of household i at the end of current time slot is represented by e i ∈ R N . Its evolution in time is given by e i = e i 0s i , where e i 0 represents the energy on the battery at the beginning of current time slot.

Flexibility operations must respect the following constraint:

max[e i 0 -E i , s i , e i ] ≤ s i ≤ min[e i 0 -E i , s i , e i ] (2)
bounding battery usage to the intersection of three segments: a) Energy capacity: the battery cannot be discharged or charged beyond its minimum or maximum energy levels, E i and E

i , which usually correspond to a SoC around 10-20% and 80-90% respectively, in order to preserve its lifetime. The scheduling mechanism limits the energy capacity available for hour-ahead trading to a range of SoC between 25% and 75%, reserving a share of the capacity for the RTC mechanism, which can use a depth of discharge between 10% and 90%. b) Power capacity: maximum energy that can be discharged/charged on a single time slot, s i and s i respectively. 5 All but one of the vector components are equal to zero. In future work this will not be the case, due to electric vehicles that can connect to any PCC. c) Cycling capacity: e i and e i are the charging and discharging cycling constraints. For instance, the hour-ahead scheduling mechanism limits battery usage to one daily cycle, while up to an extra cycle per day can be exploited by the RTC mechanism. This extra capacity is allocated to individual time slots in proportion to the corresponding day-ahead forecasted production, in order to adapt the flexibility available during sun hours when forecast error variability is higher.

For batteries being operated myopically in real-time under such constraints, the degradation costs can be expressed as

C i e (s i ) = a i (s i -ŝi ) T (s i -ŝi ) + b i T (s i -ŝi ) (3)
where s i and C i e represent the use of demand flexibility and its corresponding costs, regardless of the source of flexibility being the battery or the elasticity of demand. In the case of batteries, a i = a i s and b i = b i s are parameters that must be chosen adequately to ensure that the costs for an extra use of the battery are compensated by the economic incentives, depending on battery characteristics [START_REF] Ma | A distributed charging coordination for large-scale plug-in electric vehicles considering battery degradation cost[END_REF] and on the flexibility budget reserved for the RTC mechanism. For the case of flexible loads, a i = a i l and b i = 0 represent a purely quadratic disutility function as was already proposed in [START_REF] Samadi | Optimal real-time pricing algorithm based on utility maximization for smart grid[END_REF], [START_REF] Gkatzikis | The role of aggregators in smart grid demand response markets[END_REF].

D. Final electricity prices

With respect to prices, we need to take into account the two components of the hour-ahead commitment xi , which are the quantity xi m traded with neighbors at price P i m and the quantity xi u contracted with the supplier at price P u (or at Feed-In Tariff (FIT) P f if the prosumer is selling). The price to pay for the effectively exchanged flow x i will be the following:

C i p (x i ) = (x i -xi m )p + xi m P i m (4)
where

p = P u , P f , if x i >= xi m , else.
Such a price allocation models the fact that a household i that consumes or injects more than what was agreed on the market is considered to have enforced the market transaction; the excess of consumption or injection will be charged or payed at the price P u or P f respectively. Otherwise, a house that consumes 6 less than agreed on the market will pay x i m P i m , as if she had consumed the committed quantity, but will receive (x i mx i )P f , which corresponds to automatically selling the difference to the supplier at a lower price (P f < P m < P u ).

E. Reward allocation

The reward attribution depends on the deviation from the hour-ahead committed quantities as follows:

C i r (x i ) = a i r (x i -xi ) T (x i -xi ) -R (5)
where the value R is defined one hour-ahead by the DSO, for instance as a percentage of the transacted quantity. The reward gets reduced quadratically as the agent i deviates from its hour-ahead committed quantity xi . Note that if the deviation is too big the reward becomes a penalty.

F. Aggregative Constraint

The strategies of households are coupled by an aggregative constraint that enables DSO to indirectly enforce voltage deviations and thermal constraints. The constraint over the average strategy is as follows:

c ≤ 1 H H i=1 x i ≤ c (6)
where c and c ∈ R N are the maximum permitted aggregated flows on each PCC in order to conservatively enforce voltage deviations and thermal constraints. We denote as C ⊂ R N , the set of 1 H H i=1 x i such that ( 6) is satisfied. The constraints are obtained by the DSO through load flow sensibility analysis [START_REF] Borghetti | Short-term scheduling and control of active distribution systems with high penetration of renewable resources[END_REF]. For instance, for a period of high injection we progressively increase the aggregate injection and we allocate the increase to individual household in proportion to their injection on the previous time slot. When a constraint is detected, the injections on the previous iteration are considered as the maximum power transit supported by the grid for the time slot. This procedure requires specific knowledge about the grid infrastructure (nodes, lines, impedances, etc.), which is not necessarily available for distribution grids. Such information can be estimated by using different measurements of power injection, voltage and/or current variations [START_REF] Zhou | Simplified calculation of voltage and loss sensitivity factors in distribution networks[END_REF]- [START_REF]Learning topology of the power distribution grid with and without missing data[END_REF]. This would require short measurement campaigns after which a model of the grid can be obtained. Machine learning techniques can be applied to detect when a new measurement campaign for a model update is necessary.

G. Problem definition

Each agent i's goal is to find a strategy (x i , s i ) such that

x i ∈ arg min y J i (y, u) (7)
s.t.

(1), (2),

y ∈ [-x i , x i ] (8)
where, J i (y, u) is the local cost function that depends on the individual flow y and on u( 1 H H j=1 x j , λ), a control (price 7 ) signal broadcasted (on the transactive platform) by the DSO, which depends on the average strategy flow of houses and on λ, a penalty associated to the coupling constraint in [START_REF] Cecati | Combined operations of renewable energy systems and responsive demand in a smart grid[END_REF]. The parameter x i in constraint (8) corresponds to the maximum power capacity contracted with the supplier by household i. We denote with X i the set of feasible values of x i .

The cost function J i : R N × R N → R ∪ ∞ is defined as:

J i (x i , u) = f i (x i ) + u T x i (9)
7 negative in case of peak of production and positive in a peak of demand.

where the function

f i (x i ) = C i e (x i ) + C i r (x i ) + C i p (x i
) reflects the costs that depend only on the individual strategy flow, which are the sum of battery/comfort degradation costs, reward allocation and final electricity pricing8 .

V. COMPETITIVE AGGREGATIVE GAME

The problem described above forms a competitive aggregative game, as the optimal response of an agent depends on the aggregate response of the rest of players and they all share a common penalty associated to a coupling constraint. Grammatico proved in [START_REF] Grammatico | Dynamic control of agents playing aggregative games with coupling constraints[END_REF] that such a game has an aggregative equilibrium under assumptions on functions {f i } H i=1 being l-strongly convex and on compactness, convexity and Slater's qualification [START_REF] Boyd | Convex optimization[END_REF] of the sets {X i } H i=1 and C. Note that we need C ⊆ 1 H H i=1 X i to avoid the operator fixing constraints that are not attainable with the feasible responses of households. For this we assume that the DSO can estimate the set 1 H H i=1 X i or that, for instance, households could be asked to communicate their flexibility budget for every time slot in which the real-time control mechanism is activated. Aggregative equilibrium is defined as follows.

Definition 1:: Aggregative equilibrium. A tuple ((x i * ) H i=1 , λ * ) is an Aggregative equilibrium, for the game described in [START_REF] Shi | Real-time energy management in microgrids[END_REF] with the coupling constraint in [START_REF] Cecati | Combined operations of renewable energy systems and responsive demand in a smart grid[END_REF] if

1 H H i=1 x i * ∈ C,
and for all i ∈ H,

x i * ∈ arg min y∈X i f i (y) + D 1 H y + H j =i x j * + λ * T y. (10)
A. Iterative Process Such an equilibrium can be attained by following the semidistributed approach proposed by Grammatico [START_REF] Grammatico | Dynamic control of agents playing aggregative games with coupling constraints[END_REF]Section III]. The mechanism relies on households responding optimally to the incentive signal u. The iterative process to attain an aggregative equilibrium starts at the beginning of the 10 minutes time slot with a period of around 10 seconds 9 , which should enable households to adapt to penalty changes, while providing enough iterations for the semi-distributed algorithm to converge to the equilibrium. At the beginning of k-th iteration, the operator will estimate the aggregated optimal response A(u (k-1) ) on the previous iteration by measuring the flows on each PCC of the corresponding feeder. Then it will update the signal u (k) = κ(t, u (k-1) ), using the dynamic control law κ proposed in [START_REF] Grammatico | Dynamic control of agents playing aggregative games with coupling constraints[END_REF]Section III]. The signal will then reach the prosumers through the next published block on the transactive platform. At the detection of a change in signal u, each agent will proceed to update its optimal response x i * (k) , by updating its battery optimal flow s i * (k) from (1). In each iteration, prosumers decide on the energy flow to exchange with the grid on the whole 10 minutes, while relying on forecast done at the beginning of the time slot and on measurements of energy flows and battery usage up to the current iteration. The relationships on the optimal flows for the entire time slot and the flows up to iteration k are as follows

x i * (k) = x i 0k + x i kT (11a) s i * (k) = s i 0k + s i kT (11b)
where x i 0k corresponds to the flow exchanged with the grid up to the k-th iteration and x i kT to the future energy flow up to the end of the time slot. The notation is the same for battery flows. Note that we assume households have access to the cumulated energy exchanged with the grid x i 0k , through the smart meter, and to the corresponding battery usage s i 0k , through the smart inverter/battery controller.

The model we apply is deterministic, as it does not consider that the state of agents can change during the negotiation process towards the equilibrium. It relies then on the forecast done by households at the beginning of each time slot to be perfect. An analysis of the sensibility to very short term forecast errors and the possibility of applying stochastic game theory are subjects of future research.

Under such a deterministic scenario, each household will define the battery setting for the rest of the time slot as follows

s i kT = r i -x i * (k) -s i 0k (12)
where ( 1) and (11b) were used. Here we assume that households set their batteries to the value s i kT so that if the flows are maintained up to the end of the time slot the total flow exchanged with the grid during the time slot would be

x i * (k) .
The iteration procedure towards an aggregative equilibrium is summarized as follows Algorithm 1 Dynamic control of competitive optimal responses [START_REF] Grammatico | Dynamic control of agents playing aggregative games with coupling constraints[END_REF] Initialization: t ← 0;

• The DSO chooses u (0) ;

Iterate until convergence:

• DSO broadcasts u (0) to all agents

• Each agent i ∈ H computes x i * (u (t) ),
• and define battery setting for the rest of the time slot s i tT .

• DSO measures average best response A(u (t) ),

• obtains u (t+1) = κ(t, u (t) )

t ← t + 1

VI. SIMULATIONS AND RESULTS

A. Procedure and system scenario

The simulation of the real-time control mechanism relies on local energy market results and on an allocation of households to phases updated by a DPS mechanism. Those two additional mechanisms are simulated following the procedures described in [START_REF] Horta | Novel market approach for locally balancing renewable energy production and flexible demand[END_REF] and [START_REF] Horta | Augmenting DER hosting capacity of distribution grids through local energy markets and dynamic phase switching[END_REF] respectively.

We consider the forecast used by households for the hourahead market to be a day-ahead persistence forecast, and the updated forecast just before the time slot to be perfect, which allows us to apply a deterministic model.

With respect to the monitoring period, the RTC mechanism is only launched if during the previous time slot we detected a voltage deviation of 9% or a current flow above 70% of the thermal limit. To simulate such monitoring we realize load flow studies for every time slot in high sun hours (from 10 am to 2 pm), and when a risk of constraint is detected the RTC mechanism is launched.

We consider a distribution grid composed of one feeder with 50 households, 60% of which are equipped with a battery and 80% with a PV panel. All households participate on the market and are eligible for dynamic phase switching. With respect to the level of aggregation of households on each logical PCC, we consider one PCC per phase of the feeder. This means that all the houses on the same phase are considered to share the same logical PCC. This keeps the model simple and independent of the physical structure of the feeder and the distribution of households. Furthermore, the measurements needed by the DSO to implement the real-time control mechanism are currently available at this level of aggregation without the need of deploying further measurement devices.

B. Simulation tools and parameters

To implement the mechanism described in Algorithm 1, we use MATLAB [START_REF]MATLAB, version 8.6.0[END_REF] together with the Gurobi [START_REF]Gurobi optimizer reference manual[END_REF] optimization suite. For the power flow analysis, we rely on the Distribution Network Simulation Platform (DisNetSimPl) developed by EDF R&D. We rely on an electricity network model conformed by a 20 kV/410 V transformer of 160 kVA rated power and a feeders with 70 mm 2 aluminum power lines.

The parameters used for the simulations were the following: Load profiles -: We use realistic synthetic consumption data obtained from the Multi-agent Simulator of Human Behavior SMACH [START_REF] Amouroux | Simulating human activities to investigate household energy consumption[END_REF] as input for the household optimization problem. The load curves correspond to 7 winter days of consumption on a 1-minute basis from 50 households of mixed profiles.

Production profiles -: We consider the same synthetic production curve for all the PV panels (all located on the same area), but the profile varies on a daily basis.

Batteries -: We consider ideal batteries of 9 kWh total capacity. Up to 6 kWh (10% to 77%) are allowed to be used for local market exchanges in order to preserve the battery lifespan while reserving a share (1.2 kWh) of the capacity to the RTC mechanism, which can use up to 7.2 kWh (10% to 90%) of the capacity and perform an additional cycle 10 for absorbing forecast errors. For the case of flexible appliances, whose elasticity is entirely reserved for the RTC mechanism, we consider up to 1 kWh of energy is available for coping with forecast errors 11 .

Electricity prices -: We consider a common supplier offering a Time Of Use pricing with two levels: 15 c/kWh from 12 am to 4 pm and 20 c/kWh from 5 pm to 11 pm. While the FIT is considered to be 10 c/kWh.

C. Results and discussion

We analyze the performance of the RTC mechanism during the time slots with high renewable energy production (between 10 am and 2 pm). In particular, we are interested in cases where we observe that thermal or voltage limits are breached, or close to the limits, for the scenario without RTC. From the 7 days analyzed on the simulations we focus only on the most critical one. First, we analyze voltage deviations with respect to the nominal value U n and then the reduction of current intensity over the lines. We compare three cases, one without any control from the DSO in real-time, one with DPS only and one applying the RTC mechanism coupled with DPS. For the first two cases households optimize their electricity bill with the signal u being 0 for all time slots.

1) Reduction of Voltage deviations: We start by analyzing the impact of the RTC mechanism in voltage deviations along the feeder. In Figure 1, we see four surfaces, the red (planar) surface shows the voltage limit during the period in which RTC was active (from time slot 64 up to 74), the dark blue one shows the voltage deviations observed without any DSO control, the light blue one with only DPS being applied and the green one applying the RTC mechanism coupled with DPS. While for the cases without RTC we can observe the voltage limits being violated (when the corresponding surfaces go above the limit), the application of RTC combined with DPS achieves consistent reductions of (maximum) voltage deviations. This can be observed as the green surface (RTC + DPS) goes below the others starting from time slot 64 and up to time slot 74. These results demonstrate the capacity of RTC combined with DPS for avoiding voltage limit breaches. 11 A higher amount could be considered taking into account that only a water heater can consume more than 5 kWh a day 2) Reduction of current flows through the lines: With respect to the current flows through the lines, in Figure 2 we show the current intensity along the Feeder during the monitoring period. The red (planar) surface shows the thermal limit during the period in which RTC was active. During this period we can clearly appreciate considerable reductions of current intensity all along the feeder avoiding the violation of thermal constraints that are observed for the two cases without RTC 12 . Such reductions are enough to avoid or postpone the replacement of entire line segments, which directly translates into a considerable reduction of infrastructure investment. 3) Dynamic control mechanism performance: With respect to the performance of the dynamic control mechanism to attain the aggregate equilibrium, we observed an improvement of the convergence when combining RTC with DPS. Without the DPS mechanism the mean amount of iterations was over 60 (10 minutes), which would not be applicable to our setting without further tuning. The DPS mechanism brought the mean convergence time to 5 minutes, with a maximum of 39 iterations, which is adequate for our setting. This is probably due to the fact that our level of logical PCC aggregation is the phase, and if the loads and their DER are better distributed among phases the mechanism converges faster. For instance, if one of the phases has little or no flexibility budget, then the convergence will be slower. After convergence, each household just needs to put the battery to track the optimal flow to exchange with the grid during the rest of the time slot.

VII. CONCLUSIONS AND PERSPECTIVES

Local energy markets enable households to exchange energy with each other while increasing the amount of renewable energy that is consumed locally. Such markets increase the capacity of distribution grids for hosting renewable energies, but, as all ex-ante mechanisms, local market schedules rely on hour-ahead forecasts whose accuracy may be low. We propose a game-based real-time control mechanism to cope with forecast errors by driving households to a final exchange with the grid that benefits both the prosumer and the DSO. The proposed game has an aggregative equilibrium which can be attained in a semi-distributed manner with a number of iterations independent of the amount of households.

The performance of the mechanism is evaluated through load flow analysis and realistic load curves for a scenario with 50 households, where 80% of them are equipped with PV panels and 60% with storage. The simulations show that without the control mechanism the grid would not be capable of hosting such level of penetration of renewable energies. However, when the real-time control mechanism is applied in combination with dynamic phase allocation, maximum voltage deviations and current intensities are considerably reduced, avoiding the violations of voltage deviation and thermal constraints. We show that the proposed real-time mechanism is capable of adjusting the flows issued from local market commitments towards an equilibrium that optimizes the prosumers' electricity bill while increasing the capability of the distribution grid to support the energy transition.

The performance of the proposed real-time market mechanism can be further improved by better selecting the parameters of the dynamic control scheme and of the local cost functions. The economic performance for the DSO and the performance of the equilibrium with respect to social welfare require further analysis. With respect to the combination with the dynamic phase switching approach, the selection of households and the switching mechanism could be specifically adapted to further reduce voltage deviations. We are currently working on a sensibility analysis to very short term load/production forecast errors that could provide valuable insights for the development of stochastic game models.

Fig. 1 .

 1 Fig. 1. Maximum voltage deviations along the feeder.

Fig. 2 .

 2 Fig. 2. Current intensity level along the feeder.

Flexibility from batteries/loads reserved for absorbing forecast errors.

e.g.: start of a branch or where points of measurement are deployed.

We model the flexibility of loads as a small battery. A comprehensive model of demand elasticity is out of the scope of this work.

The logic is similar if an agent agrees to inject energy for their neighbors.

Note that C i p (x i ) are convex piecewise linear functions, while functions f i are strictly convex and l-strongly convex for any l ∈ 0, 2(a i + a i r ) .

While the time between blocks on the blockchain-based transactive platform described in[START_REF] Horta | Novel market approach for locally balancing renewable energy production and flexible demand[END_REF] is of around 5 seconds.

This would be similar to allowing batteries to cycle up to two times in a day. Nevertheless, we preallocate the flexibility budget to individual time slots proportionally to the forecasted energy production, without the possibility of accumulating the budget that is not used.

The rebound effect after the end of the RTC period, as in most demand response applications, happens when current and voltage levels are safer.