MSc Eric Vorger
email: eric.vorger@mines-paristech.fr

PhD Patrick Schalbart

PhD Bruno Peuportier

Integration of a Comprehensive Stochastic Model of Occupancy in Building Simulation to Study how Inhabitants Influence Energy Performance

Building energy simulation (BES) is currently used to design comfortable and energy efficient buildings, e.g. by comparing architectural alternatives. However, usual assumptions on occupancy are too simplistic and do not correspond to real situations. According to performance monitoring experiments, measured buildings consumptions are generally higher than predicted ones. Among the causes of these differences, the role of occupant's behaviour is identified as particularly important. In the tools used by professionals, occupancy is modelled by conventional ratios (e.g. number of persons per m²) and profiles. Besides being inaccurate this representation leaves no place for diversity. On the contrary, the stochastic model presented here takes into account households' and inhabitants' variability in terms of socio-demographic characteristics, schedules, use of electrical appliances, and adaptive behaviour. Instead of one simulation with conventional scenario leading to a unique energy consumption value, a series of simulations is conducted and yields a statistical distribution. For each simulation, virtual households are created through a probabilistic procedure according to dwellings' properties and each occupant is defined by a set of characteristics (age, sex, employment status, etc.). These characteristics condition households' electrical appliances ownership and occupants' activity scenarios, generated through a stochastic model calibrated on data from a French Time Use Survey (TUS). The activity scenarios are used as inputs to simulate the use of electrical appliances, and predict adaptive actions (depending on external and ambient conditions). A case study on a residential building located in Lyon (France) illustrates how the energy consumption probability distribution obtained after a thousand simulations can be used in a process of energy performance guarantee (EPG).

INTRODUCTION

Many monitoring experiments performed in the past years on new or retrofit buildings revealed that buildings energy consumptions exceeded the previsions of the simulation tools. This gap can partly be explained by differences between the real building envelope and the simulated one but it seems that the effective conditions of occupancy play an even more important role. [START_REF] Andersen | The influence of occupants' behaviour on energy consumption investigated in 290 identical dwellings and in 35 apartments. Healthy Buildings[END_REF], for example, highlighted this point by observing very dispersed annual heating consumptions for similar dwellings under the same climate. Moreover, the role of occupant's behaviour is particularly important in efficient buildings whose performance is very sensitive to solar and internal heat gains as well as to heat losses due to ventilation. In the simulation tools, occupancy is modelled by conventional ratios and profiles. A new approach is necessary to increase the relevance of simulation results (in terms of energy consumption and thermal comfort), avoid design errors and represent behavioural diversity. Controlling the uncertainties related to occupancy is also essential to develop a process of EPG, currently a major issue for the building sector. An original model of occupancy that tackles this issue in the case of residential buildings is described in the second part. The third one presents its application to the issue of EPG. Suggested improvements for the model and possible applications are discussed in the conclusion.

DESCRIPTION OF THE INTEGRATED BEHAVIOURAL MODEL

The proposed integrated behavioural model is composed of several submodels in order to: -create realistic virtual inhabitants defined by a set of socio-demographic characteristics (no equivalent model has been found in the literature),

-generate inhabitants' activity scenarios according to their characteristics, using the detailed model developped by [START_REF] Wilke | A bottom-up stochastic model to predict building occupants' time-dependent activities[END_REF], which was calibrated with data from a TUS, -equip households, and describe appliances' duty cycles and standby powerdata from surveys were used to define ownership probabilities depending on households' characteristics and measurements were used to generate a diversity in appliances' characteristics, while reference models [START_REF] Richardson | Domestic electricity use, A highresolution energy demand model[END_REF][START_REF] Tanimoto | Validation of methodology for utility demand prediction considering actual variations in inhabitant behaviour schedules[END_REF][START_REF] Widén | A high-resolution stochastic model of domestic activity patterns and electricity demand[END_REF] consider the average national ownership rate and a unique power use for each type of appliance, -simulate the use of appliancesdata from measurement campaigns made possible a calibration for each type of appliance, while previous models used a residual electrical power to comply with aggregated dwelling's electricity consumptions, -simulate actions of occupant's on windows using Haldi and Robinson's model. The way the different parts of the model are articulated and coupled to the dynamic BES software Pleiades+Comfie is summarised in Figure 1.

Figure 1

Overall scheme of the integrated behavioural model coupled to the BES.

By repeating the simulation according to the Monte Carlo method, the model generates a distribution of the energy consumption relative to occupancy instead of the usual unique value. Temperature setpoints control is not integrated (indoor temperatures were measured in the building considered for the case study) nor actions on shading devices (which are difficult to model in residential buildings since the driving factors are numerous and vary from case to case : the role of temperature and privacy can be important, whereas the link with visual comfort is strongly prominent in office buildings).

Stochastic creation of virtual households

The first submodel creates virtual inhabitants for every housing at the beginning of a simulation. Since series of simulations are conducted, the same housing is stochastically inhabited by different virtual households at each simulation (e.g. a five-persons family, or a single occupant, etc.). Each inhabitant is described by a set of socio-demographic characteristics that are consistent:

-for a given individual (e.g. the retired status must be in agreement with the age), -between individuals of the same households (e.g. if the first individual is a married man, father of one child and having a high income level, his wife and child must be generated too, their ages must be consistent and they must share his living standard), -between inhabitants and the housing, for which the simulation tool user is always able to provide the location (urban/rural), the number of rooms and the type (detached house/apartment/other).

Individuals' characteristics are strongly intercorrelated, therefore they are determined sequentially. From the number of rooms, the type of housing and the location there are probabilities for different number of members of the household; knowing the size of the household there are probabilities for different types of households; knowing the type of household there are probabilities for the age and gender of his reference member, etc. The whole procedure draws on large datasets from the last population census in France and associated surveys.

Generation of activity scenarios for each occupant

Presence and activities modeling allows localising the inhabitants in the different rooms and assigning corresponding metabolic heat gains. But more importantly, it is a necessary input to simulate the adaptive actions and the use of appliances and lighting. Wilke et al.'s model (2013), based on Markov processes and multinomial logit models, was adapted to provide weekly scenarios (then reproduced identically for every week of the simulation) for every inhabitant. In this model, the probability to start an activity and the associated distribution of duration both depend on the day of the week, the time of the day, and the socio-demographic characteristics of the individuals. To illustrate its outputs, Figures 2 and 3 respectively present the average daily profile of the activities for a population of 15,441 people (corresponding to the individuals of the TUS) and a daily scenario for a random occupant.

Figure 2

Average daily profile of the activities for a population of 15,441 individuals. Example of a daily scenario for a random individual (the colors legends are the same as in Figure 2).

Before linking activity scenarios to appliances operation, households have to be equipped and their equipments described. For all these steps, the example of the dishwasher is used as an illustration.

Households' electrical appliances ownership

This submodel equips households with appliances, according to probabilities relative to their characteristics. More than thirty appliances are modeled, including several instances of the same equipment (e.g. principal and secondary television or computer have different characteristics). According to the type of appliance, probability of ownership corresponds to the average national rate, or depends on the age of the reference person 1 , the type of household, its income or a combination of these factors. In the case of dishwashers, probabilities of ownership are determined as functions of the type of households and its income (c.f. Table 1) according to results from (INSEE, 2012).

Description of appliances' duty cycles

This submodel assigns a duty cycle (constituted of a square signal or a succession of square signals) and a standby power to each appliance. Appliances of a same type can have various operating powers according to their technology, dimension and preferred mode of use. Data are mostly provided by several comprehensive measurement campaigns realised by the engineering office ENERTECH in partnership with EDF and ADEME 2 . For instance, 65 dishwashers' consumptions were measured with a 10 min resolution during one year (ENERTECH, 2008). Dishwashers' duty cycles can be divided into four phases. Water is heated during phases 1 and 3 while cold water is used in phases 2 and 4. A fifth phase of active drying is optional. Each appliance is associated with a nominal power ranging from 1,500 to 2,100 W 3 . Phases correspond to a percentage of this power and a duration (both potentially variable). In agreement with observations of the campaign, an appliance in a specific household is associated to a unique duty cycle. Table 2 summarises the description of the operating phases. A seasonal coefficient is affected to the nominal power to reflect the cold water temperature variation. The three values of the coefficient, corresponding to winter, summer and midseason (respectively 1.09, 0.82 and 0.98), were obtained by dividing the measured average daily consumption of the period by the measured annual average daily consumption. Figure 4 shows a good agreement between the modeled cycles consumptions and the measured ones.

Figure 4

Comparison between modeled and measured cycles consumptions.

Simulation of the use of appliances -Starting probabilities calibration

The use of appliances is linked to occupants' activities through pragmatic hypotheses. In the case of dishwashers, it is assumed that they operate just after the end of a meal. Therefore, every time a period of "Meals and snacks" ends, a random number is sampled on [0, 1] and compared to the dishwasher starting probability to determine wether a cycle begins. For a few appliances, among which dishwashers, time delay functions are integrated to allow shifting the starting time. In the case of dishwashers, delays mostly concern after-dinner cycles, when people do not need to use the appliance before the next morning. Shifting is then used to benefit from lower electricity prices of the night. To have a number of cycles during the night in coherence with measurements it is supposed that 25 % of the cycles between 6 p.m. and 11 p.m. are shifted after 11 p.m., with an average delay of 4 hours. A main objective of the model is to accurately reproduce the load variations during a day. To do so, the starting probabilities are calibrated as functions of the hour, so that the average daily load curve obtained through a great number of simulations fits the average daily load curve taken from the measurements. Figure 5 illustrates the obtained similarity of the two curves (the root mean square error per hour is 2.3 Wh) using the calibrated starting probabilities presented in Table 3.

Figure 5

Comparison of the modeled and measured average daily load curves.

Simulation of the use of appliances -Evaluation of the ability of the model to produce diversity

By choosing this way of calibration, it is ensured that the model reproduces an average number of cycles and an average annual consumption that are in good agreements with measurements. It is also necessary to evaluate its ability to produce a realistic diversity. The proper behaviour of the model on that point is illustrated by Figure 6 which shows, for example, that 10 % of the dishwashers have a consumption higher than 500 kWh/year.

Figure 6

Comparison of the modeled and measured annual consumptions.

Application of the electrical consumption model including artificial lighting modeling

Artificial lighting use is modeled similarly to appliances, except that switch on probabilities depend on the hours of sunrise and sunset. The calibration data are provided by annual measurements of every light bulbs in 100 housing units in France (ENERTECH, 2004).

As an illustration of the whole electrical consumption model outputs, Figure 7 presents a load curve of a random four room apartment for the first week of the year (Monday 1 th .january to Sunday 7 th ). Appliances are regrouped by type (e.g. dishwashers are part of the "washing" group with washing machines, tumble dryers, clothes irons and vacuum-cleaners).

Figure 7

Example of a weekly electrical load curve for a four room apartment.

One can observe for instance the importance of standby consumption of the audiovisual part for this housing and the concentration of peak loads around 8 a.m., noon and 7 p.m. The derived internal heat gains scenario presents the same trends with lower peaks because only part of the electric energy consumption is converted into heat gains (e.g. for dishwashers 40 % of the energy is lost as greywater).

Occupants' actions on windows

Whereas the previous submodels were part of a preprocessing that may be seen as a "scenario generator", windows' openings and closings interact with the thermal calculation at each time step since such actions depend on inside and outside temperatures, while air flows influence the heat gains or losses. Haldi and Robinson's model (2009) is used for this part. Its predictions address the binary action "open/close" but gives no indication on the opening rate. Therefore, when an opening happens, an entering airflow is sampled between 5 and 15 vol/h and assumed constant during the opening duration.

APPLICATION TO EPG OF A RESIDENTIAL BUILDING

One objective of the integrated behavioural model is to assess the uncertainties related to occupancy in order to develop a process of EPG. In this context, we studied a 4 storeys residential building with a 1,048 m² area located in the suburb of Lyon. It includes 16 apartments, modelled by a thermal zone each, plus a central circulation corresponding to a seventeenth zone (Figure 8 a). The exercise carried out in this case study was the following one: this building was retrofitted including an insulation of the outside walls; heating consumption, air exchange, and indoor temperature were measured before the renovation; the goal was to estimate the heating consumption after the renovation with a confidence interval allowing the project manager to make a commitment on a guaranteed performance. Then, this performance was compared to the actual consumption, which was measured during the year following the renovation. The uncertainty due to occupants' behaviour was calculated with the Monte Carlo method (1,000 simulations including the integrated behavioural model were conducted to draw a smooth gaussian distribution of the heating loads but fewer are required to reach the convergence on the average value and the standard deviation). The thermal characteristics of the building after renovation are given in Table 4. Results are summarised in Figure 8 b. We first consider that the heating temperature was maintained to its value before renovation (approached by a 20.2 C constant scenario). Using this heating scenario and scenarios from the French regulatory method for presence and internal gains4 , the annual heating loads were 38 kWh/m², whereas the Monte-Carlo method leads to an average value of 45 kWh/m² with a standard deviation of 3.5 kWh/m² (blue curve). Accepting a risk of exceeding of 2.5 %, the guaranteed performance would be 52 kWh/m². The average value for annual heating loads and the guaranteed performance threshold were respectively 18 % and 37 % higher than the value obtained through the regulatory calculation method.

The actual loads (62 kWh/m²) were superior to the guaranteed performance value, essentially because of a "rebound effect" on the heating setpoint : average temperatures measured in winter after the renovation approached 23 °C (compared to 20.2 °C before). Guarantee obviously cannot longer apply if the requirements of the inhabitants increase after the renovation. For such a heating setpoint (red curve), the guaranteed performance would have been 75 kWh/m². A comprehensive stochastic model of occupants' behaviour in residential buildings is proposed. It integrates an original model for the creation of virtual individuals described by a set of sociodemographic parametres. This allows a high degree of refinement in the generation of schedules and in the attribution of equipment to households according to statistical data. The use of appliances and lighting is modeled on the basis of inhabitants' activities with a higher accuracy than existing models from the literature, through data from several large measurement campaigns. A reference model for interactions of occupants with windows was adapted. The whole model is coupled to a dynamic BES tool with no more necessary input than the building description (but any available information on inhabitants' characteristics or equipment can be filled by the user). The distribution of the simulation outputs is than obtained using the Monte-Carlo method. In the case study, it was observed that the monitoring of temperature setpoints in winter is essential in a process of EPG (knowing that it can be easily and reliably implemented at low cost). If not, considering the possibilities of evolutions of the temperature setpoints leads to a wide distribution, unexploitable in this context. However, this aspect being under control, a commitment on a guaranteed value becomes possible. In this case the guaranted value would be 38 % higher than the value calculated with conventionnal deterministic scenarios.

In further studies the whole coupled model will be encapsulated within the statistical software R. The influence of different parametres of the integrated model will be analysed to identify the aspects on which attention should be paid in terms of improvement and data collection. The outputs of interest will include heating loads but also domestic electricity consumptions and comfort indicators.

Figure 3

 3 Figure 3

Figure 8

 8 Figure 8 (a) Thermal zones of the building (design interface of the simulation tool Pleiades+Comfie). (b) Heating loads distributions for two temperature setpoints corresponding to measuremements before and after renovation.

Table 1 . Dishwashers' Probabilities of Ownership

 1

	Household type\Income quartile	1 st	2 nd	3 rd	4 th
	One person household	25 %	30 %	35 %	40 %
	Couple alone	60 %	65 %	70 %	80 %
	Couple with children	65 %	80 %	85 %	90 %
	Single-parent family	35 %	45 %	45 %	55 %
	Other type of household	45 %	50 %	50 %	70 %

Table 2 . Dishwashers' Duty Cycles

 2

	Phases	N°1	N°2	N°3	N° 4	N° 5
	% of the nominal power	100	10 -30	40 -80	5	40 -80
	Duration (min)	20	10 -60	20	10	5

1

The "reference person" is a demographic concept which indicates a kind of "head of the household". 2 EDF : Eléctricté De France ; ADEME : Agence De l'Environnement et de la Maîtrise de l'Energie 3 Maximal power can be higher but peak loads occur on periods much shorter than the 10 min time step of the model. A superior temporal resolution isn't necessary to evaluate heat gains in a building or even to model the electricity demand of several housing units because of thermal mass and averaging effects. 30th INTERNATIONAL PLEA CONFERENCE 16-18 December 2014, CEPT University, Ahmedabad

Table 3 . Dishwashers' starting probabilities at the end of a meal as a function of time

 3

	0h -8h30	8h30 -13h	13h -14h	14h -19h	19h -21h	21h -24h
	0.10	0.16	0.22	0.17	0.25	0.15

Table 4 . Thermal characteristics of the building

 4

	Facades	U = 0.22 W/(m².K)	Thermal bridges	Psi = 356 W/K
	Loggias	U = 0.23 W/(m².K)	Mechanical ventilation flow	0.3 vol/h
	Low floor	U = 0.49 W/(m².K)	Air infiltrations flow	0.15 vol/h
	Roof terrace	U = 0.29 W/(m².K)	T° setpoint before works	20.2 °C
	Windows	U = 1.9 W/(m².K)	T° setpoint after works	22.8 °C
	CONCLUSION

16-18 December 2014, CEPT University, Ahmedabad

30th INTERNATIONAL PLEA CONFERENCE 16-18 December 2014, CEPT University, Ahmedabad

Presence and internal gains during the night and the week-end. Using also the conventionnal scenario for heating (19 °C during presence), the calculated annual heating loads were 29 kWh/m².


M Robillart

P Schalbart

B Peuportier

Extraction of heating control rules from the dynamic programming method for load shifting in energy-efficient building

Keywords: Rule extraction, optimal control, load shifting

Extraction of heating control rules from the dynamic programming method for load shifting in energy-efficient building

M Robillart, P Schalbart, B Peuportier

INTRODUCTION

Existing control schemes

In modern construction of buildings, the main objectives for the control systems are to save energy [START_REF] Ferguson | Predictive thermal control of building systems[END_REF], to increase comfort [START_REF] Mathews | Developing cost efficient control strategies to ensure optimal energy use and sufficient indoor comfort[END_REF] and to reduce peak electricity demand [START_REF] Greensfelder | An investigation of optimal control of passive building thermal storage with real time pricing[END_REF]. To meet such objectives, control systems have to be able to anticipate the weather, the occupancy, and the solar and internal gains. Dounis and Caraiscos (2009) reviewed many advanced control systems meeting such objectives. For instance, during a summer period, control systems are used to maintain comfort using passive cooling (Braun et al, 2001), to reduce energy consumption of air conditioning [START_REF] Chahwane | Advanced control systems engineering for energy and comfort management in a building environment-A review[END_REF], or to control solar protections [START_REF] Nielsen | Quantifying the potential of automated dynamic solar shading in office buildings through integrated simulations of energy and daylight[END_REF]. During a winter period, control systems are used to decrease the energy consumption of the heating system [START_REF] Le | Gestion optimale des consommations d'énergie dans les bâtiments[END_REF] or to reduce peak demand [START_REF] Malisani | Investigating the ability of various buildings in handling load shiftings[END_REF].

Load shifting

Recently, numerous efforts have been made to reduce electricity peak-demand. In Europe, these peaks mostly appear during winter periods and are due to heating systems. For example in France, the building sector represents 68 % of the final electricity consumption (ADEME, 2012). To guarantee the grid stability, some studies have been done on electrical load shifting.

Thanks to electricity demand-side response (DSR), the consumer demand for energy can be modified through various methods such as financial incentives or education. Many economical models are used by the demand side response programs. Two categories may be distinguished: time based programs and incentive based programs [START_REF] Marwan | Demand Side Response to Mitigate Electrical Peak Demand in Eastern and Southern Australia[END_REF]Federal Energy Regulatory Commission, 2006). Examples of application of time based programs are Time Of Use (with fixed electricity prices for off-peak and peak hours), or Real Time Pricing (with variable electricity tariffs). For incentive based programs, an example is the Direct Load Control, which allows to turn specific appliances on and off during peak demand periods.

At the level of the individual houses, the electricity peak reduction can be achieved thanks to a careful architectural design to efficiently manage solar gains [START_REF] Ferguson | Predictive thermal control of building systems[END_REF]). An advanced control system can also be used to reduce heating consumption. Such control can be based on power tariff [START_REF] Hämäläinen | Cooperative consumers in a deregulated electricity market -dynamic consumption strategies and price coordination[END_REF]Pineau et Hämäläinen, 2000) or the use of the thermal mass of the building to shift part of electricity consumption [START_REF] Wyse | Investigation into the Time-Shifting of Domestic Heat Loads[END_REF][START_REF] Hong | The influence of thermal storage on microgeneration flexibility[END_REF]. For instance, [START_REF] Favre | Application of dynamic programming to study load shifting in buildings[END_REF] used the dynamic programming method to shift the building consumption. The proposed method consisted in over-heating the building in the hours before the peak based on weather forecast, and occupancy and internal gains schedules for the next 7 days. However, this exact optimisation method is time-consuming and can hardly be applied to real-time control.

Rule extraction

In developing an operational strategies framework, exact optimisation results can be used to extract simplified control rules that are implementable in real-time.

This approach was first applied in water resource management. The application was to developed simplified control rules for reservoir management based on the results of offline model predictive control (MPC) [START_REF] Wei | Optimal tree-based release rules for real-time flood control operations on a multipurpose multireservoir system[END_REF]. The approach has recently been used in the building context. For instance, [START_REF] May-Ostendorp | Extraction of supervisory building control rules from model predictive control of windows in a mixed mode building[END_REF] used many data mining techniques (generalised linear models, classification and regression trees and adaptive boosting) to extract rules from offline MPC results for a mixed mode building operated during the cooling season. To our knowledge, this approach was never applied to shift the heating consumption in building.

The present study is based on the results of [START_REF] Favre | Application of dynamic programming to study load shifting in buildings[END_REF]. and its objective is to develop operational strategies to shift the heating load in building. A new methodology is proposed to extract decision models from dynamic programming results and then compare them.

MODELS

Thermal model of the building

The building is modelled considering spatial zones of homogenous temperature. For each zone, each wall is meshed according to the finite volume technique with a uniform temperature and thermal capacity. Another mesh is added for the zone's air and furniture. Energy conservation equations are written on each mesh within the building and form a system of equations:

(1) with the thermal capacity of the node , -the temperature of the node , -the solar and internal gains (due to heating, occupancy and other appliances), -the heat losses by conduction, convection and radiation.

Repeating energy conservation equations for each mesh leads to a linear time-invariant system [START_REF] Peuportier | Simulation tool with its expert interface for the thermal design of multizone buildings[END_REF], temporal variation terms being added in the simulation:

(2) with the node temperature vector, -the diagonal thermal capacity matrix, -the driving forces (climate parameters, heating, etc.), -the output vector (indoor temperatures accounting for air and wall surfaces), -, , , the state, input, output and feedforward matrices, respectively.

In order to perform simulation, it is important to know the occupancy of the building which defines the emission of heat by the inhabitants and appliances, and the thermostat setpoint influencing the heating equipment. Another important aspect is the weather model influencing heat losses and solar gains. The data regarding house occupancy and weather models were included in the driving forces vector .

The high order linear model (2) needed to be reduced because its state dimension was too large to allow a fast convergence of the optimisation algorithm. A reduction method (modal reduction) was thus applied to lower the state dimension. In this work, the building energy simulation tool COMFIE was used [START_REF] Peuportier | Simulation tool with its expert interface for the thermal design of multizone buildings[END_REF].

Optimisation algorithm

The dynamic programming method was developed by Bellman (1957). It is a sequential optimisation method which examines all possible ways to solve an optimisation problem and provides, given a discretisation, an optimal set of commands over a period.

To apply dynamic programming, a state variable describing the system is used and discretised temporally:

∈ , ⊂ (3)
where is the set of possible states and the dimension of . The control vector can be chosen in a set ⊂ (the set of possible controls) where is the dimension of the control vector:

∈ , ⊂ (4)
One can act on the system state through the control variable . The state space equation of the dynamical system . is thus:

, 1 , (5)
A value function is defined, which is the cost to go from to 1 :

,

Under these assumptions, a finite-horizon decision problem takes the following form:

max , (7)
subject to the constraints (3) and ( 4) and the state space equation ( 5). denotes the optimal value that can be obtained by maximising the objective function subject to the assumed constraints. The dynamic programming method is then applied to break this decision problem into smaller sub-problems. Bellman's principle is thus used: "An optimal policy has the property that whatever the initial state and initial decision are, the remaining decisions must constitute an optimal policy with regard to the state resulting from the first decision". Equation ( 7) becomes:

max , , (8)
Figure 1 shows how the dynamic programming operates: In Figure 1, the state variable is discretised in five states. Optimisation aims at minimising the state variable. In this example, thanks to Bellman's principle, four solutions can be discarded either because they trespass the constraints or because they reach the same states as solutions with lower costs.

For the application of dynamic programming in building context, the chosen state variable and the cost function are defined in §5.2.1.

Rule extraction: Generalised linear model

The generalised linear model (GLM) framework was used to derive simplified decision models from the dynamic programming results allowing a small computational expense adapted to real time control. GLM models the relationship between regressors (explanatory variables) and response . It consists of three elements:

-a random component (the response is assumed to be generated from a particular probability distribution), -a deterministic component (a linear combination of explanatory variables ), -a link function (that provides the relationship between the linear combination of explanatory variables and the mean of the distribution function).

We have thus to estimate the following model:

(9) with - . the link function, -
the expected value of .

The unknown parameters and are typically estimated with maximum likelihood [START_REF] Gill | Generalized Linear Models[END_REF].

METHODOLOGY

The following section describes techniques employed to extract decision models from dynamic programming results. Dynamic programming was used to generate training data and validation data to identify a GLM's parameters and to evaluate its performance, respectively. It was done in two stages.

As a first step, Test Reference Year-type (TRY) weather data were used to perform optimisation using dynamic programming and to elaborate an optimal strategy. This optimal strategy (training data) was then used to identify the GLM's parameters with the same weather data.

As a second step, the predictive capacity of the model was measured on locally recorded weather data. We compared the GLM's results with the optimal strategy calculated by dynamic programming (validation data).

Model identification

Model identification was done in a four-step process (Figure 2). First, all data used by dynamic programming were collected (Test Reference Year-type weather data, electricity tariff, occupancy). Secondly, the optimal strategy was elaborated using dynamic programming (training data). Thirdly, we identified the GLM's parameters thanks to optimal strategy. Fourthly, the resulting control models were implemented within the simulation platform (COMFIE).

Model comparison

Model comparison was done in a three-step process (Figure 3). Firstly, all data used by the optimisation method and the GLM were collected (local weather data, electricity tariff, occupancy). Secondly, we performed optimisation with dynamic programming and GLM to determine optimal strategy and operational strategies respectively. Finally, we compared performances of operational strategies against the optimal strategy.

Building description

The building under study is a single-family house based on an actual experimental passive house being part of the INCAS platform built in Le Bourget du Lac, France, by the National Solar Energy Institute (INES). The studied house has two floors and a total living floor area of 89 m² (Figure 4). The North facade has only two small windows whereas 34 % of the South facade is glazed. The building's façades include double ( 1.1 . . , 0.6) and triple on the North ( 0.7 . . , 0.45) glazing windows of various dimension. The south façade also includes solar protection for the summer period. The house is highly insulated with a high thermal mass as shown in Table 1. The main thermal bridges as well as the building air tightness have been carefully designed and implemented. The house is heated by an electrical resistance integrated in an efficient heat recovery ventilation system. According to thermal simulation results using the thermal model described in §3.1, the annual heating load is 14 kWh.m -2 .

Optimisation parameters

Dynamic programming parameters

The chosen state variable for dynamic programming is the total energy stored in the building, calculated as follows [START_REF] Favre | Application of dynamic programming to study load shifting in buildings[END_REF]) :

(10) with the reference temperature chosen at 0°C, -the number of nodes An upper and lower bound of this state variable was defined according to its initial value.

Then it was discretised in 800 nodes.

To ensure thermal comfort in the building, indoor temperature had to be maintained between 19°C ( ) and 26°C ( ). We considered a typical four people family occupancy: the building was non-occupied only during the working days from 8 a.m. to 5 p.m.. Each occupant emitted 80 W due to their metabolism, and internal gains from appliances were also considered during occupied hours. The heating power was in the range of 0 W ( ) to 5000 W (

).

The model of the building was mono-zonal and the optimisation was done over 34 days (which was the maximal decision problem's horizon solved by dynamic programming), with one hour time step, to generate training and validation data. The goal of dynamic programming was to minimise the heating cost of the building by determining a set of commands (heating power ) with constraints on thermal comfort and heating power. Thus, the finite-horizon decision problem took the following form:

Electricity tariff

To shift electricity demand, a time-of-use pricing was considered (Table 2):

Weather data

Meteonorm data from Chambery (to generate training data) and local weather conditions data (to generate validation data) measured at the Chambéry airport which is 300 meters away from the building, were used to perform simulations. Test Reference Year-type (TRY) were used to develop GLM because these weather data represent the typical long-term weather patterns. Thus, GLM's results were adjusted with the long-term average climatic conditions. Then, we used local weather conditions data to assess GLM's behaviour in real conditions.

Meteorological features are summarised in Table 3.

Skill evaluation

Objective criteria for evaluating the predictive quality of the model were required. Therefore, the following indicators were used to assess its performance:

-the mean absolute error (MAE), between heating powers calculated by dynamic programming and GLM, -the average heat power, -the cumulative cost, -the percentage of high peak hours which are load shifted, -the percentage of peak hours which are load shifted, -the thermal discomfort rate representing the number of hours when the indoor temperature falls below 19°C (in %), -the thermal discomfort rate representing the number of hours when the indoor temperature rises above 26°C (in %).

RESULT ANALYSIS

Explanatory variables

Explanatory variables that can be measured in building were used to develop GLM. Thus, to determine the heating power at time step Δ , we used explanatory variables at time step Δ : outdoor temperature , global horizontal irradiance and electricity tariff ( . Explanatory variables at time step were also considered: indoor temperature , and heating power.

Models developed

To apply generalised linear model (GLM), we had to define the link function. That is why we changed the response variable as a proportion of maximum heating power (5000 W). For example, a heating power at time step Δ of 2500 W, corresponds to a predicted variable by GLM of 50 %. The statistical model used by GLM is thus a multiple logistic regression and the link function is the logit function ln . This model was used to relate the proportion of maximum heating power to predictor variables through the logistic link function:

ln 1 (14)
Five models were developed, each one using all or some of the training data (Table 4).

Training data were divided into three groups: off-peak hours training data (TD OPH ), peak hours training data (TD PH ) and high peak hours training data (TD HPH ).

In the implementation, the heating power at time step Δ was set at 0 W in the following cases:

-during high peak hours for GLM_2, GLM_3, GLM_4 and GLM_5 models, -during peak hours for GLM_3 and GLM_5 models.

These choices were done in order to ensure load shifting during peak and high peak hours. However, in the implementation, the heating power at time step Δ was set at 0 W during high peak hours et peak hours for both models. The same logic was applied for GLM_2 and GLM_4 models.

Results

Each GLM model was implemented in the building energy simulation tool COMFIE. Table 5 summarises GLM models' results obtained on validation data. The dynamic programming reference results are described in the DP column.

The resulting model predictions of GLM_1 and GLM_2 and the original optimised sequence are presented in Figure 5. We can clearly observe that GLM_1 and GLM_2 did not follow the dynamic programming's behaviour. For example, we can see that the GLM_2 model performed significantly worse than dynamic programming, with a very high thermal discomfort rate (93 %) and an indoor temperature exceeding 30°C (30.6°C). GLM_1 had a similar behaviour with a significant cumulative cost (137 €) and a high mean absolute error (111 %).

Figure 6 shows that predictions of GLM_4 and GLM_5 are also different from the optimised results. For instance, we can see that GLM_4 and GLM_5 had a significant cumulative cost (105 € and 92 € respectively) and a high average power (1347 W and 1309 W respectively). Moreover, GLM_4 and GLM_5 had a relatively large mean absolute error (88 % and 71 % respectively).

However, Figure 7 illustrates the interesting behaviour of GLM_3. Firstly, due to its design, no electricity was consumed during high peak hours and peak hours. Secondly, it had a cumulative cost and an average heat power close to dynamic programming (72 € and 1023 W compared to 68 € and 936 W for DP). Thirdly, its mean absolute error (40 %) and mean relative error (9 %) were reasonable. Finally, GLM_3's computational time was 150 times smaller than the dynamic programming method, using a desktop computer. The GLM_3 model presented a satisfactory behaviour and seemed a possible candidate to be used as simplified control system. However, on some occasions, it did not respect the thermal comfort constraints. Therefore, an improved controller was considered that switched heating on as soon as the indoor temperature was below 19°C.

Application controller

An ideal on-off controller was considered. It was applied during peak and high peak hours as GLM_3 did not work during these periods. Its control law switched between the minimum heating power (0 W) and the maximum heating power (5000 W). The ideal on-off controller was switched on when the indoor temperature fell below 19°C -ɛ and was switched off when the indoor temperature rose above 19°C + ɛ (in order to respect the 19°C set point temperature). Assuming that ɛ tended toward 0, the deadband of the on-off controller ( ε tended toward 0. The use of this ideal on-off controller aimed at assessing maximum performance of GLM_3 + controller.

Electricity tariff

Dynamic programming GLM_3

The control law was the following:

-During off peak hours

P t Δt GLM_3 T t T 0 T t (15)
-During peak and high peak hours

P t Δt Controller on T t 19° Controller off T t 19°C ε ε → 0 (16)
The obtained results are shown in Table 6. We can notice the interesting behaviour of GLM_3 + controller. Firstly, thanks to the on-off controller, GLM_3 + controller respected temperature constraints (the lowest temperature reached was 19°C). Then, we can see a slight deterioration of peak hours and high peak hours shifted. For example, GLM_3 + controller had 92 % of high peak hours which were load shifted in comparison with the 100% of GLM_3 (and the 99 % of dynamic programming). Similarly, GLM_3 + controller had 95 % of peak hours which were load shifted. It was less efficiency than GLM_3 (100 %) but it was better than dynamic programming (88 %). Finally, GLM_3 + controller had a cumulative cost (72.9 €) and an average heat power (1029 W) close to GLM_3 and dynamic programming. Consequently, adding an on-off controller with GLM_3 enabled to improve the GLM_3's behaviour and to respect temperature constraints. Figure 8 shows the GLM_3's behaviour both with and without the on-off controller. To plot GLM_3+controller's graph, the heating power was averaged over one hour time periods.

CONCLUSION

Dynamic programming method has been used to study load shifting of heating systems in an energy-efficient building. Due to its computational expense, a statistical technique (generalised linear model) has been introduced that allow for the extraction of logistic decision models from the dynamic programming results. This method models the relationship between explanatory variables and a response variable. The results showed that generalised linear models were able to imitate the general characteristics of the dynamic programming results, with a much smaller computational expense and limited overshooting of the setpoint.

To improve the GLM's behaviour, an on-off controller was added that switched heating on as soon as the indoor temperature did not respect temperature constraints. The results showed that the GLM+on-off controller respected temperature constraints and that there were a slight deterioration of peak hours and high peak hours shifted. Therefore, rule extraction (generalised linear model) is a promising technique for developing operational control strategies. Given their simple mathematical formulation, generalised linear models coud be implemented in real time building systems control.

NOMENCLATURE

Latin

Figure 1 :

 1 Figure 1: Dynamic programming description

Figure

  Figure 2: Model identification

Figure 4 :

 4 Figure 4: view of the house (west and south façades)

  at time step , -_ the electricity cost at time step t. -the duration of the optimisation period

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Heating power calculated by dynamic programming, GLM_1 and GLM_2 (Third week)

Figure 8 :

 8 Figure 8: Heating power calculated by dynamic programming, GLM_3 and GLM_3 + controller (Third week)

  This work has been supported by the French Research National Agency (ANR) through the "Villes et Bâtiments Durables (VDB) 2012" program (project ANR-12-VBDU-0006).

  . « Chiffres clés du bâtiment ».Bellman, Richard. 1957. « E. 1957. Dynamic Programming ». Princeton UniversityPress. BellmanDynamic programming1957. Braun, James E., Kent W. Montgomery, et Nitin Chaturvedi. 2001. « Evaluating the Performance of Building Thermal Mass Control Strategies ». HVAC&R Research 7 (4): 403-428. doi:10.1080/10789669.2001.10391283.

Table 1 :

 1 Building description

		External wall	Ground	Intermediate floors	Ceiling	Interior partition
				16 cm
	Composition	15 cm thick 20 cm of extruded polystyrene	20 cm insulation concrete slab 20 cm external	concrete concrete slab screeds and girders 12 cm	40 cm of glass wool	4 cm of glass wool
				floor
	U (W.m -2 .K -1 )	0.15	0.15	2.2	0.09	0.96

Table 2

 2

		: Electricity prices
		Off-peak hours	Peak hours	High peak hours
	Hours	12 a.m. to 9 a.m.	9 a.m. to 5 p.m. 10 p.m. to 12 a.m.	5 p.m. to 10 p.m.
	Cost per kWh (€)	0.0864	0.1275	0.255

Table 3

 3

	: Weather data
	Training data	Validation data
	Minimum temperature (°C)	-9.10	-14
	Average temperature (°C)	1.44	-0.22
	Maximum temperature (°C) Average global horizontal irradiance (W.m -2 ) Maximum global horizontal irradiance (W.m -2 )	11.50 58 486	11.33 60 569

Table 4 :

 4 Training data GLM_3 and GLM_5 models were different because they did not have the same training data. Indeed, GLM_3 was trained only on off-peak hours training data (TD OPH ) whereas GLM_5 was trained on complete training data (TD OPH , TD PH , TD HPH ).

		Off-peak hours	Peak hours	High peak hours
	GLM_1	TD OPH	TD PH	TD HPH
	GLM_2	TD OPH	TD PH	-
	GLM_3	TD OPH	-	-
	GLM_4	TD OPH	TD PH	TD HPH
	GLM_5	TD OPH	TD PH	TD HPH
	As a more specific example,

Table 5 :

 5

			GLM models' results
		GLM_1	GLM_2	GLM_3	GLM_4 GLM_5	DP
	Average heat power (W)	1353	1811	1023	1347	1309	936
	Cumulative cost (€)	137	158	72	105	92	68
	High peak
	hours load	0	100	100	100	100	99
	shifted (%)
	Peak hours
	load shifted	0	0	100	0	100	88
	(%)
	TI min (%)	0	0	8	0	0	0
	TI max (%)	3	93	0	2	0	0
	T min / T max	21.2 / 26.15	21.4 / 30.6	18.4 / 23.8	20.7 / 26.1	19.6 / 26 19 / 23.4
	MAE (%)	111	153	40	88	71	-

Table 6

 6

		: GLM_3 + controller results
		GLM_3	GLM_3 + Controller	DP
	Average heat power (W)	1023	1029	936
	Cumulative cost (€)	72.2	72.9	68
	High peak
	hours load	100	92	99
	shifted (%)
	Peak hours
	load shifted	100	95	88
	(%)
	TI min (%)	8	0	0
	TI max (%)	0	0	0
	T min / T max	18.4 / 23.8	19 / 23.8	19 / 23.4
	MAE (%)	40	41	-

9th International Conference on System Simulation in Buildings, Liege, December 10-12, 2014


H Apriyanto
email: haris.apriyanto@enseeiht.fr

G Ravet

O D Bernal

M Cattoen

F Lizion

H C Seat

V Chavagnac

A Multimode Fiber Refractive Index Sensor

Keywords: Refractive index sensor, Refractometer, Multimode fiber sensor, Evanescent wave, Stripped cladding I

We present a simple refractive index sensor based on multimode fibers realized by stripping three different lengths of cladding. We theoretically explain and experimentally validate three mechanisms occurring in the multimode fiber for sensing surrounding refractive index. This sensor has been demonstrated at 1550 nm for a wide range of refractive index variation from 1.3164 to 1.608. Our sensor is very sensitive for measuring refractive index values lower than but close to the refractive index value of the fiber core.

INTRODUCTION

The refractive index (RI) sensor, also commonly known as the refractometer, is widely used in the food, chemical, biomedical industries, as well as in environmental monitoring, etc. Due to their good sensitivity, RI sensors have been used for measuring methane concentration by using cryptophane based traps in polymeric PDMS to absorb the methane molecules leading to refractive index variation from 1.412 to 1.413 with a sensitivity 5.5 RIU/nM at 830 nm of wavelength [START_REF] Boulart | A novel, low-cost, high performance dissolved methane sensor for aqueous environments[END_REF].

The principal techniques employed for measuring RI, an optical phenomenon, consist of prisms [START_REF] Astrua | Prism refractive index measurement at INRiM[END_REF][START_REF] Zhou | Intensity modulated refractive index sensor based on optical fiber Michelson interferometer[END_REF], fiber-based surface plasmon resonance (SPR) [START_REF] Pathak | SPR Based Fiber Sensor to Measure Refractive Index of Glycerol and Acetone[END_REF], microstructured fiberbased SPR [START_REF] Santos | SPR microstructured D-type optical fiber sensor configuration for refractive index measurement[END_REF], Fresnel reflection at an optical fiber interface [START_REF] Xu | Simple fiber-optic refractive index sensor based on fresnel reflection and optical switch[END_REF][START_REF] Zhao | High-resolution and temperature-insensitive fiber optic refractive index sensor based on fresnel reflection modulated by Fabry-Perot interference[END_REF][START_REF] Zhao | Simple fresnel reflection-based optical fiber sensor for multipoint refractive index measurement using an AWG[END_REF], fiber Bragg grating [START_REF] Iadicicco | Thinned fiber Bragg gratings as refractive index sensors[END_REF][START_REF] Iadicicco | Thinned fiber Bragg gratings as high sensitivity refractive index sensor[END_REF], fiber ring laser [START_REF] Zhang | Refractive index sensor based on fiber ring laser[END_REF], nanometric plasmonic [START_REF] Jin | A novel nanometeric plasmonic refractive index sensor[END_REF], multimode fiber [START_REF] Fukano | High-sensitivity fiber-optic refractive index sensor based on multimode interference using smallcore single-mode fiber for biosensing[END_REF][START_REF] Chen | Multimode fiber Mach-Zehnder interferometer for measurement of refractive index[END_REF][START_REF] Fiber | High refractive index liquid level measurement via coreless multimode fiber[END_REF], and tapered multimode optical fiber [START_REF] Kim | High resolution refractive index sensing with cladded multimode tapered optical fibre[END_REF].

In this paper, we present a simple yet elegant multimode fiber sensor for measuring RI over a relatively wide dynamic range where the sensing area is a section of stripped fiber cladding. The impact of the length of stripped cladding section is also investigated to estimate the sensor's sensitivity over the useful RI range. The advantages of this sensor are its simplicity, cost-effectiveness, wide range and high sensitivity for RI values lower than but near to the RI of the fiber core.

II. PRINCIPLE OF OPERATION

In wave optics, the propagation of the injected laser beam in the optical fiber can be represented in terms of modes. To define the number of modes (M), the V-parameter can be used and is a function of the relative RI difference between the fiber core (𝑛 𝑐𝑜 ) and cladding (𝑛 𝑐𝑙 ), operational wavelength (𝜆), and core radius (a), governed by the relationship [START_REF] Saleh | Fundamental of Photonics[END_REF]

 𝑀 = 4 𝜋² 𝑉² 
where V is given by 𝑉 = 2𝜋𝑎 𝜆 𝑛 𝑐𝑜 ² -𝑛 𝑐𝑙 ² .

In a step-index fiber, each mode that propagates in the optical fiber has a unique propagation constant which is proportional to the effective guide index (n eff ) and wavenumber (k) [START_REF] Wartak | Computational Photonics: An Introduction with MATLAB[END_REF] given by

𝛽 = 𝑘𝑛 𝑒𝑓𝑓 (2)
However, through ray optics analysis, this mode can also be represented as a ray propagating in the core region at a certain angle (θ a ), as shown in Fig. 1, with where the possible angles within which the ray can propagate in the core should respect θ c ≤ θ ≤ 90°. The guided modes are thus within 𝑛 𝑐𝑜 𝑠𝑖𝑛 𝜃 𝑐 ≤ 𝑛 𝑐𝑜 𝑠𝑖𝑛 𝜃 ≤ 𝑛 𝑐𝑜 𝑠𝑖𝑛 90, corresponding to possible effective indices of 𝑛 𝑐𝑙 ≤ 𝑛 𝑒𝑓𝑓 ≤ 𝑛 𝑐𝑜 . The number of guided modes, and hence the transmitted optical power, in the fiber are thus highly dependent on RI variation of the outer cladding layer, potentially rendering the MMF a highly sensitive refractometer. For RI measurements, a section of the MMF outer cladding layer is removed and the exposed fiber core serving as the sensing region is arranged to be in contact with the external medium subject to RI variation.

𝑛 𝑒𝑓𝑓 = 𝑛 𝑐𝑜 𝑠𝑖𝑛𝜃 (3)
The length of the stripped cladding section is here investigated as a function of the MMF sensitivity to RI variations and is correlated to 3 possible operating conditions occurring in the sensing region with respect to the RI of the external sensing medium (n ms ), namely: a) n ms ≤ n cl , b) n cl ≤ n ms ≤ n co , and c) n co ≤ n ms .

For n ms ≤ n cl (or Zone I), the RI of the medium is detected via the phenomenon of evanescent wave absorption that exists during propagation when the modes are guided by total internal reflection (TIR). The evanescent field is generated at the core-cladding interface, and penetrates into the cladding (Fig. 2) to a certain depth (d p ) [START_REF] John | Evanescent Wave Fibre Optic Sensors : Design, Fabrication and Characterization[END_REF] given by θ c θ

𝑑 𝑝 = 𝜆 0 2𝜋 𝑠𝑖𝑛 2 𝜃- 𝑛 𝑐𝑙 𝑛 𝑐𝑜 2 (4)
Further, the evanescent wave is attenuated in the cladding. The cladding attenuation coefficient γ is proportional to the total attenuation coefficient α and to the fraction of the guided power r, 𝛾 = 𝑟𝛼 , with 𝑟 = 𝑃 𝑐𝑙𝑎𝑑𝑑 𝑃 𝑡𝑜𝑡𝑎𝑙 . Since we removed the cladding and replaced it by a medium for sensing, the power transmitted at the sensing length (L) is given by For the second condition, illustrated in Fig. 3, where n cl ≤n ms ≤ n co (Zone II), the RI variation is due to 2 phenomena: a) modification to the critical angle since 𝜃 𝑐 = 𝑎𝑟𝑐 sin(𝑛 𝑐𝑙 / 𝑛 𝑐𝑜 ), leading to the variation of the number of guided modes, and b) the evanescent wave absorption by the external medium since TIR still occurs under these conditions. Hence, there should be a stronger variation of the transmitted power than Zone I over the same variation in RI. Finally, when n co ≤ n ms (Zone III), the laser beam in the MMF core is propagated by external reflection. This phenomenon occurs since the core index is lower than that of the cladding at an incident angle close to 90°. This implies that the external reflection only occurs for lower order modes. In this condition, no evanescent wave interacts with the sensing medium as illustrated in Fig. 4. Since the sensing medium has higher RI than the core, a portion of the ray at the boundary will be transmitted to the exterior while the rest will be reflected back into the core. Using Snell's law, the reflected angle (θ r ) which is equal to the incident angle (θ t ), and the transmitted angle (θ t ) are related by The resulting total reflectivity R is thus given by

𝑃 𝐿 = 𝑃(0)exp (-𝛾𝐿) (5)
𝑅 = 1 2 (𝑟 𝑝 ² + 𝑟 𝑠 ²) (7)

III. EXPERIMENT SETUP

The sensing fiber employed is a multimode plastic-cladsilica (PCS) fiber with a numerical aperture (NA) of 0.48, and core and cladding diameters of 200 µm and 230 µm, respectively. The cladding layer is removed to obtain three different sensing lengths of 1 cm, 2.5 cm and 4 cm. Further, in order to excite all possible modes in this MMF, a microscope objective (MO) with an NA of 0.65 is employed.

To detect the RI of the medium, the variation of the optical power transmitted in the fiber is measured. In addition, to compensate ambient temperature perturbations, a differential set-up is exploited, with one MMF serving as the reference fiber (P ref ) and the second as the sensing fiber (P sens ), as shown in Fig. 5. The RI variation of the medium is obtained from a combination of glycerol and water for RI from 1.3164 to 1.4569 (at a wavelength of 1550 nm) [START_REF] Saunders | The refractive index of common solvents and solutions at 1550 nm[END_REF], and calibrated oils for RI beyond 1.4569.

IV. RESULTS AND DISCUSSION

The RI measurements are carried out by normalizing the power transmitted by the sensing MMF with respect to the reference MMF over the RI range 1.3164 -1.608. From the available NA and fiber core data (n co ~1.44402), the fiber cladding index, n cl , is estimated to be 1.362.

The experimental results are plotted in Figs. 6 and7 for three different lengths of sensing region. To facilitate the analysis of the results, both curves are divided into the three conditions (or zones) as previously explained, i.e. a) n ms ≤1.362, b) 1.362≤n ms ≤ 1.44402, and c) n ms ≥ 1.44402. Subsequently, to estimate the sensitivities of this sensor, we derived the function of the curves in Fig. 6 for each of the sensing lengths, the results of which are plotted in Fig. 7. According to Fig. 7, for Zone I, a better sensitivity is achieved by the 4 cm stripped cladding, compared with the 1 cm and 2.5 cm stripped claddings. This confirms the exponentially decreasing transmitted power due to evanescent wave absorption as a function of the stripped cladding length described by [START_REF] Santos | SPR microstructured D-type optical fiber sensor configuration for refractive index measurement[END_REF]. The highest sensitivity is achieved at the end of Zone II with a -32.93 RIU -1 , -28 RIU -1 and -22.98 RIU -1 for the 1 cm, 2.5 cm and 4 cm stripped cladding, respectively. This sharp decrease is due to the larger losses of the lower order modes. However, at the beginning of Zone II (1.365-1.417), which can be applied for methane sensing, the longer cladding section (4 cm) has a better sensitivity than the others. For Zone III, since external reflection occurs in the sensing region, most of the optical power is transmitted to the exterior of fiber, thus significantly reducing the guided power in the core. According to ( 6) and ( 7), the guided optical power increases for increasing RI. Due to a lower number of reflections (N) at the coremedium boundary, the guided power is thus higher for the shorter stripped length. By adapting [START_REF] Zhao | High-resolution and temperature-insensitive fiber optic refractive index sensor based on fresnel reflection modulated by Fabry-Perot interference[END_REF], the guided power can be calculated by 𝑃(𝑁) = 𝑃 0 (𝑅) 𝑁 , with P 0 the incident power. For a fiber with a core diameter a and a sensing length L, N≈L/(a tan θ i ). Since the value of R in [START_REF] Zhao | High-resolution and temperature-insensitive fiber optic refractive index sensor based on fresnel reflection modulated by Fabry-Perot interference[END_REF] is always less than 1, more reflections lead to a reduced guided power. Figs. 6 and7 show that in the case of external reflection, the shorter stripped clad fiber has a better sensitivity.

V. CONCLUSION

A refractive index sensor based on the multimode fiber was conceived to measure refractive indices from 1.3164 to 1.608 at a wavelength of 1550 nm. We described three sensitivity zones with respect to the RI of the cladding and the core of the fiber. In the first zone, a stripped clad fiber of 4 cm has a better sensitivity than both the 1 cm and 2.5 cm devices. The sensitivity is highest at the end of the second zone, and is better for the 1 cm stripped cladding. Finally, for the third zone, less power is guided but increases for increasing refractive index, with a better sensitivity obtained for a shorter length of the stripped clad.

Fig. 1 .

 1 Fig. 1. Principle of propagation in step-index fiber.

Fig. 2 .

 2 Fig. 2. Evanescent wave in the cladding region (Zone I).

Fig. 3 .

 3 Fig. 3. Contributions from critical angle variation, and evanescent wave absorption for ncl≤nms≤ nco (Zone II).

Fig. 4 .

 4 Fig. 4. External reflection in MMF (Zone III).

  coefficients for the P-and S-polarizations are given by Fresnel relations[START_REF] Taylor | Fresnel Equations[END_REF]:𝑟 𝑝 = 𝑡𝑎𝑛 (𝜃 𝑖 -𝜃 𝑡 ) 𝑡𝑎𝑛 (𝜃 𝑖 +𝜃 𝑡 ), 𝑟 𝑠 = 𝑠𝑖𝑛 (𝜃 𝑖 -𝜃 𝑡 ) 𝑠𝑖𝑛 (𝜃 𝑖 +𝜃 𝑡 )[START_REF] Xu | Simple fiber-optic refractive index sensor based on fresnel reflection and optical switch[END_REF]

Fig. 5 .

 5 Fig. 5. Schematic diagram of the experimental set-up.

Fig. 6 .

 6 Fig. 6. Optical power ratio (Psens/ Pref) versus RI variation.

Fig. 7 .

 7 Fig. 7. Sensitivity of the sensor versus RI variation (P=Psens/Pref).

Sensing region (stripped cladding section)

ACKNOWLEDGMENT

Haris Apriyanto acknowledges the Indonesia Endowment Fund for Education (LPDP) for a PhD scholarship.