Julien Eynard
email: julien.eynard@univ-perp.fr

Stéphane Grieu
email: grieu@univ-perp.fr

Monique Polit
email: polit@univ-perp.fr

Wavelet-based multi

Keywords: wavelet-based multi-resolution analysis, artificial neural networks, outdoor temperature, thermal power consumption, hot water distribution network, district boiler

come

Introduction

The actual European energy context reveals that the building sector is one of the largest sectors of energy consumption. In France, about 25% of Greenhouse Gases (GHG) emissions and 45% of energy consumption are due to buildings [START_REF]et de la Maîtrise de l'Energie, Les chiffres clés du bâtiment en 2006[END_REF]. Consequently, the adopted "Energy Performance of Buildings Directive" (Official Journal of the European Communities, 2002) [START_REF][END_REF], focusing on energy use in buildings, requires all the European Union (EU) members to enhance their building regulations and to improve energy efficiency. With the aim of surmounting the actual energy crisis (mainly caused by both the rarefaction of fossil fuels and excessive energy consumption), managing energy demand, promoting renewable energy and finding ways for energy savings are worldwide concerns [START_REF]Europe's energy position. Present & future[END_REF].

The present paper, the first of a series of three dealing with the modeling and the optimization of a district boiler as part of the OptiEnR project, focuses on the development of a forecasting tool using a wavelet-based multi-resolution analysis and artificial neural networks. The OptiEnR project began in late 2008 and will proceed, at least, until late 2010. It includes researchers from the ELIAUS laboratory of the University of Perpignan Via Domitia (south of France) and engineers from two French companies, Cofely GDF-SUEZ [4] and Weiss France [START_REF] France | [END_REF]. Cofely GDF-SUEZ is currently the leading European brand for environmental and energy efficiency services while Weiss France designs, manufactures and installs automatic boiler rooms for all types of wood, biomass, wastes and fluids. The considered district boiler, which will be described later in the paper, is situated at La Rochelle, west coast of France, and managed by Cofely GDF-SUEZ. It supplies domestic hot water and heats residential and public buildings, using mainly wood and sometimes fuel or gas if necessary. The OptiEnR project focuses on optimizing the performance of the boiler, adding to the plant a thermal storage unit and using a model-based predictive optimal controller. Its main objective is to minimize the use of fossil energy, stocking renewable energy during low-demand periods and using it during high-demand periods. The project consists of the following successive tasks (Figure 1): the first task is to forecast outdoor temperature and thermal power consumption (of the hot water distribution network); the second task focuses on modeling the district boiler; the third task deals with studying the feasibility of integrating to the plant a thermal storage unit and finding the most adequate storage material; the fourth task is to model the thermal storage while the fifth task focuses on optimizing the boiler functioning using both a Model Predictive Controller (MPC) [START_REF] Garcia | Model predictive control: Theory and practive. A survey[END_REF][START_REF] Qin | A survey of industrial model predictive control technology[END_REF] and the forecasted parameters (first task). As abovementioned, the present paper deals with the first task, the next papers of the series dealing with the tasks 2 (paper 2) and 3 to 5 (paper 3). Schemes for hourly temperature forecasting have been mainly developed in the context of short or long-term load forecasting and power utilities management. The complex and nonlinear nature of temperature variations as well as the abundance of historical data suggest that computational intelligence data-based modeling techniques would be good candidates to forecast temperature. In 1996, Liu et al. [START_REF] Liu | Comparison of very short-term load forecasting techniques[END_REF] evaluated fuzzy logic, artificial neural networks and auto-regressive models for very short-term load forecasting. They concluded that the just-mentioned tools can be good candidates for this application. In 1996 and 1997, Khotanzad et al. [START_REF] Khotanzad | An artificial neural network hourly temperature forecaster with applications in load forecasting[END_REF][START_REF] Khotanzad | ANNSTLF: a neural-network-based electric load forecasting system[END_REF] also used backpropagation neural networks to forecast hourly temperature for the next seven days, using daily high and low temperatures and an adaptive daily update of the weights. The main drawback of this approach is that the large size of the neural networks involved leads to a large number of weights to be optimized. In 1998, Hwang et al. [START_REF] Hwang | Power load forecasting by neural network with a learning process for considering overtraining problem[END_REF] also used artificial neural networks and proposed a new and effective learning algorithm allowing solving both the problems of underfitting and overfitting. Yoo and Pimmel [START_REF] Yoo | Short term load forecasting using a self-supervised adaptive neural network[END_REF] developed in 1999 a self-supervised adaptive neural network to perform short term load forecasts. They used this kind of neural networks to extract correlational features from temperature and load data. The levels of error compared favorably with those given by other techniques. In 2000, Sharif and Taylor [START_REF] Sharif | Real-time load forecasting by artificial neural networks[END_REF] used separate artificial neural networks for load forecasting of one hour to four hours ahead. We can also highlight that Gonzàlez and Zamarreño [START_REF] Gonzàlez | A short-term temperature forecaster based on a novel radial basis functions neural network[END_REF] used in 2001 a radial basis functions neural network to develop a short-term temperature forecaster, without using the daily high and low temperature values. They obtained good results. One year later (2002), the same authors, Gonzàlez and Zamarreño [START_REF] Gonzàlez | A short-term temperature forecaster based on a state space neural network[END_REF], proposed a short-term hourly environmental temperature forecaster based on a state space neural network. Its training was based on a random optimization method. Because of the non-stationary characteristic of temperature, training was executed daily with the aim of updating the network weights. The outdoor temperature dynamics was satisfactorily captured. However, and, as it was highlighted by Hippert et al. in 2000 [16], a large input dimensionality compared with the number of training records can lead to unstable networks and/or overfitting. As a consequence, one can remark that Tassadducq et al. [START_REF] Tassadduq | Application of neural networks for the prediction of hourly mean surface temperatures in Saudi Arabia[END_REF] used in 2002 a backpropagation neural network to forecast the temperature at a given hour of the next day, only using the temperature at the same hour of the present day. In 2008, Chaabene and Ben Ammar [START_REF] Chaabene | Neuro-fuzzy dynamic model with Kalman filter to forecast irradiance and temperature for solar energy systems[END_REF] introduced a new methodology for dynamic forecasting of meteorological parameters. An ANFIS [START_REF] Jang | ANFIS: Adaptative-Network-based Fuzzy Inference System[END_REF][START_REF] Jang | Neuro-fuzzy modeling and control[END_REF][START_REF] Lin | Neural fuzzy systems: a neuro-fuzzy synergism to intelligent systems[END_REF] neuro-fuzzy system offered daily time distribution of irradiance and ambient temperature relying on the meteorological behavior during the days before. First, the obtained medium term forecasting was then modeled using an ARMA model [START_REF] Makridakis | Forecasting: methods and applications[END_REF][START_REF] Pandit | Time series and system analysis with applications[END_REF]. Next, a Kalman filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF][START_REF] Cathin | Estimation, control and the discrete Kalman filter[END_REF][START_REF] Hasltine | Critical evaluation of extended Kalman filtering and movinghorizon estimation[END_REF] provided short-term forecasting. The developed models were tested and validated, leading again to very good results.

Because artificial neural networks also suffer from the difficulty of finding both an optimal topology and adequate training parameters [START_REF] Grieu | Prediction of parameters characterizing the state of a pollution removal biologic process[END_REF], some different approaches were proposed. In 1957, Allen [START_REF] Allen | A simplified equation for minimum temperature prediction[END_REF] tested a simple parametric model for minimum temperature prediction. In 1985, Schneider et al. [START_REF] Scheider | 25-hour electric utility load forecasting[END_REF] fitted a two-harmonics Fourier model to the temperature data of the past 21 days to produce a temperature day profile. Hourly forecasts were then obtained by stretching or contracting this profile. Fan and McDonald [30] proposed in 1994 a load forecasting model consisting of time series, nonlinear load-weather functions and a residual function represented by an ARMA model [START_REF] Makridakis | Forecasting: methods and applications[END_REF][START_REF] Pandit | Time series and system analysis with applications[END_REF]. In 1996, Chow and Leung [START_REF] Chow | Neural netwoprk based short-term load forecasting using weather compensation[END_REF] presented a novel technique for electric load forecasting based on neural weather compensation. They proposed a nonlinear generalization of the Box and Jenkins approach [START_REF] Box | Time series analysis: forecasting and control[END_REF][START_REF] Pankratz | Forecasting with univariate Box-Jenkins models: concepts and cases[END_REF] for nonstationary time series forecasting and obtained very good results. In 2000, Chen and Hwang [START_REF] Chen | Temperature prediction using fuzzy time series[END_REF] used fuzzy time series [START_REF] Song | Fuzzy time series and its models[END_REF][START_REF] Song | Forecasting enrollments with fuzzy time series -Part I[END_REF][START_REF] Song | Some properties of defuzzification neural networks[END_REF][START_REF] Song | Forecasting enrollments with fuzzy time series -Part II[END_REF] to deal with temperature forecasting problems and overcome the drawback related to historical data represented by linguistic values. They proposed a new fuzzy time series model, called two-factor time-variant fuzzy time series model. Based on the proposed model, they developed two algorithms for temperature prediction and obtained good forecasting results. Using the same model, Lee et al. [START_REF] Lee | Handling forecasting problems based on twofactors high-order fuzzy time series[END_REF] established in 2006 two-factor high-order fuzzy logical relationships based on historical data. Again, the proposed method got a higher forecasting accuracy rate than the existing methods. In 2004, Abdel-Aal [START_REF] Abdel-Aal | Hourly temperature forecasting using abductive networks[END_REF] used abductive networks to forecast hourly temperature. Abductory inductive mechanism is a supervised inductive machine-learning tool for automatically synthesizing abductive network models from a database of inputs and outputs representing a training set of solved examples. This tool can automatically synthesize adequate models that embody the inherent structure of complex and highly nonlinear systems. Abdel-Aal developed next-day and next-hour models. Performance compares favorably with neural networks models developed using the same data, and with more complex neural networks that requires daily training. Performance is significantly superior to naïve forecasts based on persistence and climatology. In 2008, Dong and Pedrycz [START_REF] Dong | A granular time series approach to long-term forecasting and trend forecasting[END_REF] introduced the concept of granular time series. It can be mainly applied to longterm forecasting and trend forecasting to overcome the curse of dimensionality (which plagues most predictors when carrying out long-term forecasts) and cope with uncertainly present in many times series. Dong and Pedrycz used a technique of fuzzy clustering to construct information granules on a basis of available numeric data present in the original time series. In the sequel, they developed a temperature forecasting model which captures the essential relationships between such information granulates and in this manner constructs a fundamental forecasting mechanism. Dong and Pedrycz demonstrated that the model they proposed comes with a number of advantages which manifest when processing with a large number of data. The same year (2008), Lee et al. [START_REF] Lee | Temperature prediction and TAIFEX forecasting based on high-order fuzzy logical relationships and genetic simulated annealing techniques[END_REF] presented a new method for temperature prediction, based on high-order fuzzy logical relationships and genetic simulated annealing techniques [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF][START_REF] Gen | Genetic algorithms and engineering design[END_REF][START_REF] Goldberg | Genetic algorithm in search, optimization, and machine learning[END_REF][START_REF] Goldberg | Messy genetic algorithms: motivation, analysis, and first results[END_REF][START_REF] Holland | Adaptation in natural and artificial systems[END_REF]. Simulated annealing techniques were used to adjust the length of each interval in the universes of discourse with the aim of increasing the forecasting accuracy rate. The proposed method provided very good accuracy.

One can also highlight some interesting and recent works, although not directly related to temperature forecasting but having many similarities with the proposed methodology. In 2007, González-Romera et al. [START_REF] González-Romera | Forecasting of the electric energy demand trend nad monthly fluctuation with neural networks[END_REF] used artificial neural networks to extract the trend component from monthly electric demand data and then performed separate predictions of both tendency and fluctuation, which were summed up to obtain the series forecasting. A mean absolute percentage error of about 2% was obtained. The same authors, González-Romera et al., proposed a novel hybrid approach to investigate the periodic behavior of the Spanish monthly electric demand series [START_REF] González-Romera | Monthly electric energy demand forecasting with neural networks and Fourier series[END_REF]: this behavior is forecasted with a Fourier series [START_REF] Mitra | Handbook for digital signal processing[END_REF] while the trend is predicted using an artificial neural network. Satisfactory results were obtained. These results improve those reached when only neural networks or ARIMA models [START_REF] Abdel-Aal | Forecasting monthly electric energy consumption in Eastern Saudi Arabia using univariate time-series analysis[END_REF][START_REF] Saab | Univariate modeling and forecasting of energy consumption: the case of electricity in Lebanon[END_REF] were used for the same purpose. As highlighted by the authors, a correct separation of trend and fluctuation and the optimization of the forecasting tools used to carry out the predictions are the key issues in the success of this technique. Furthermore, the combined use of a neural network and a Fourier series provides a simpler structure than that with only one neural network to carry out the whole prediction.

The proposed methodology deals with the concept of time series, even if not only past values are considered for estimating future values, and uses a wavelet-based Multi-Resolution Analysis and multi-layer Artificial Neural Networks. One could speak of "MRA-ANN" methodology. The discrete wavelet transform allows decomposing sequences of past data in subsequences (named coefficients) according to different frequency domains [START_REF] Gao | Neural-wavelet methodology for load forecasting[END_REF]. From these coefficients, Multi-Layer Perceptrons (MLP) are used to estimate future subsequences of 4 hours and 30 minutes. Future values of outdoor temperature and thermal power consumption are then obtained by simply summing up the estimated coefficients. Substituting the prediction task of an original time series of high variability by the prediction, using MLP neural networks (many other neural networks have been tried but no significant improvement of the accuracy was observed), of its wavelet coefficients on different levels of lower variability, is the main idea of the present work. Then, the reconstruction of future values is performed by simply summing up the estimated coefficients. In addition, the sequences of past data are completed, for each of their components, by both the minute of the day and the day of the year to place the developed model in time. First, the paper focuses (section 2) on the considered district boiler and its working, and highlights the required upgrades. Section 3 describes both the data used to develop and validate the forecasting tool and the proposed MRA-ANN methodology. The wavelet-based multi-resolution analysis and the multilayer Perceptron are also described. Sections 4 (outdoor temperature) and 5 (thermal power consumption of the hot water distribution network) present the forecasting results. Section 4 deals also with the impact on forecast accuracy of parameters such as both wavelet order and decomposition level, and with the topology of the ANN used. The final part of the paper concludes the present work and remembers the aim of the OptiEnR project.

The district boiler of La Rochelle

The district boiler of La Rochelle, whose synopsis is shown in Figure 2, is composed of a breaking pressure bottle, a cogeneration plant and two thermal boilers. The first one, a large 4.5 MW wood boiler, uses renewable energy (it is fed with woodchip). The second one, a 7 MW gas-fuel oil boiler, uses fossil energy. Both boilers have smokes and air-to-water heat exchangers. They supply hot water to the collecting hydraulic circuit according to a temperature set-point which is defined from outdoor temperature. During the cold season (from October to May), the wood boiler is continuously running while its heating power is adapted to the demand. The gas-fuel oil boiler functions during very cold periods only, when the wood boiler fails to respond to the demand. The primary hydraulic circuit (3000 m 3 ), or "distribution network", supplies hot water to heat residential and public buildings (for example schools), for a total of 2700 accommodations. Domestic hot water is also produced, for a total of 3500 accommodations. The cogeneration plant, connected to the "return" part of the primary hydraulic circuit, produces electricity using gas and warms up the cold water before it goes back to the collecting hydraulic circuit. The breaking pressure bottle pulls apart the two hydraulic circuits, because of the difference between their respective flows. Table 1 presents the coverage rate of each of the three heat generators. Let us remember that the main goal of the OptiEnR project deals with the minimization of the gas-fuel oil boiler coverage rate, leading to the reduction of the fossil energy consumption. As a consequence, the wood boiler coverage rate will be maximized. To reach this objective, we proposed to add a thermal storage unit to the plant (Figure 2) with the aim of stocking hot water when the demand is low and using it when the demand can't be met by the wood boiler. This is a classical solution to optimize the functioning of boilers used to heat buildings [START_REF] Tanton | Improving the effectiveness of a domestic centralheating boiler by the use of heat storage[END_REF][START_REF] Stritih | Optimization of a thermal storage unit combined with a biomass boiler for heating buildings[END_REF][START_REF] Cao | Optimal design of thermal-energy stores for boiler plants[END_REF]. Thus, the gas-fuel oil boiler will be only used when both the hot water demand is very high and the thermal storage unit is empty. The flow of the water passing through this unit (

) will be adjusted thanks to the control of the thermal storage feed pump (

). With the aim of optimizing the use of this thermal storage unit and reducing the fossil boiler coverage rate, a model predictive controller will define over the next 4 hours and 30 minutes the optimal sequence of and the wood boiler set-point temperature (

), taking into consideration some parameters measured at the district boiler as well as both the forecasted outdoor temperature (

) and thermal power consumption (of the distribution network) (

). The wood boiler ( ) and gas-fuel oil boiler ( ) feed pumps are controlled using a standard on/off controller.

is the distribution network feed pump allowing controlling the water differential pressure according to outdoor temperature variations (Figure 2).

Materials and methods

Database: outdoor temperature and thermal power consumption

Outdoor temperature ( ) is measured at the district boiler and used for computing both the boilers and the primary circuit temperature set-points and the differential pressure set-point. That is why forecasting accurately this parameter, also used by the Model Predictive Controller (MPC) for managing the thermal storage unit, is of paramount interest. When taking a look at the data, one can highlight some interesting characteristics: first, outdoor temperature increases during the considered period (one speaks of "increase trend"); secondly, a 24-hour pseudo-period has been detected; finally, climatic phenomena impact on the parameter variability throughout the day. Of course, we want the proposed forecast methodology, based on a multi-resolution analysis and artificial neural networks, to be able of detecting these characteristics and taking advantage of them to make accurate forecasts. Let us also note that a meteorological station located in La Rochelle measures outdoor temperature and provides to the district boiler operators forecasted temperatures for five future days (from D+3 to D+7), with a sampling time set to 3 hours. With the aim of comparing these measured and forecasted temperatures with the temperatures measured at the district boiler, the considered vectors are upsampled from 3 hours to 5 minutes. Indeed, the sampling time of all the parameters measured and used to develop the boiler model is set to 5 minutes. Table 2 allows comparing both the forecasted and measured outdoor temperatures with the measurements done at the district boiler of La Rochelle.

Table 2. Comparison of both the forecasted and measured outdoor temperatures provided by the meteorological station with the measurements done at the district boiler. FIT: similarity criterion;

MRE: Mean Relative Error.

Outdoor temperature FIT (%) MRE (%)

Measured Table 2 highlights first that both the outdoor temperatures measured at the meteorological station and at the district boiler are quite similar. As a result, one can suppose that the temperature forecasted at the meteorological station is usable for managing the thermal storage unit using a MPC. Indeed, both the temperatures forecasted at the station and measured at the district boiler are also quite similar (Figure 3), with a FIT and a MRE ranging between 17.1% ( ) and 43.1% (

) and between 6.14% ( ) and 9.01% ( ) respectively. Unfortunately, Figure 3 highlights that the peaks of temperature are significantly under-estimated. As a consequence, the MPC could give the order to start the gas-fuel oil boiler unnecessary. So, instead of improving the district boiler performance and decreasing the use of fossil energy, more gas and fuel will be consumed. That is why, finally, the outdoor temperature forecasted at the meteorological station will not be used by the MPC, a forecast module being developed and presented in this paper. ) during the first months of 2009 (from early January to late March). This parameter is not measured in situ and has an exogenous influence on the district boiler functioning because it impacts on the temperature of the water back from the distribution network.

is calculated using the difference between (the temperature of the hot water leaving the breaking pressure bottle and entering the distribution network) and

(the temperature of the cold water coming back from the distribution network) as well as the flow of the water passing through this network (

). is estimated using the opening percentage of the distribution network feed pump valve ( ), related to the water differential pressure ( ), its set-point ( ) and outdoor temperature (Figure 2). As depicted by Figure 4, the thermal power consumption and the outdoor temperature are inversely proportional.

Figure 4. Calculated thermal power consumption of the hot water distribution network (

).

Wavelet-based multi-resolution analysis 3.2.1. Continuous Wavelet Transform (CWT)

The continuous wavelet transform is used with the aim of decomposing a signal into wavelets i.e. into highly localized small oscillations [START_REF] Daubechies | Ten lectures on wavelets[END_REF][START_REF] Heil | Continuous and Discrete Wavelet Transforms[END_REF]. Whereas the Fourier transform [START_REF] Cody | The Fast Wavelet Transform -Beyond Fourier transforms[END_REF][START_REF] Perrier | Wavelet spectra compared to Fourier spectra[END_REF] decomposes a given signal into infinite length sines and cosines, losing all time-localization information, the CWT's basis functions are scaled and shifted versions of the time-localized mother wavelet. So, this mother wavelet is a continuous in both time and frequency function that can be considered as a source function from which scaled and translated basis functions are generated. The CWT is used to obtain a time-frequency representation of a signal that offers very good both time and frequency localization. Considering a mother wavelet , one can obtain the following wavelet series , where is the scale factor and is the translation factor according to the following expression (equation 1) [START_REF] Gubner | Wavelet transforms for discrete-time periodic signals[END_REF]:

For any function , its continuous wavelet transform can be expressed as follows (equation 2): where is the complex conjugate of . is called the wavelet coefficient. The continuous wavelet coefficients can be used with the aim of reconstructing the function in the following way (equation 3): . Basically, the discrete wavelet multi-resolution analysis, commonly based on Daubechies orthogonal wavelet basis [START_REF] Daubechies | Ten lectures on wavelets[END_REF], allows decomposing a signal into approximations (i.e. low frequency coefficients) and details (i.e. high frequency coefficients) using a filter bank composed of both low-pass ( ) and high-pass ( ) filters [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. This process can be repeated times, producing levels of decomposition, but decomposing (downsampling) the approximations

Discrete Wavelet Transform (DWT) and bank of filters

(3) (2) (1)
(the low frequency coefficients) only. Indeed, the high frequency coefficients are neglected. So, a signal can be first decomposed into an approximation and a detail (that is the level 1 of the decomposition). Then can be decomposed into an approximation and a detail (that is the level 2 of the decomposition) and so on. Considering levels of decomposition, the reconstruction process allows recovering the initial signal, summing the details and the approximation of level . As an example, Figure 5 depicts the 3-level decomposition of a signal and its reconstruction [START_REF] Mallat | A wavelet tour of signal processing[END_REF].

The Daubechies wavelets family

Different families of wavelets whose qualities vary according to several criteria can be used for analyzing sequences of data points [START_REF] Percival | Wavelet methods for time series analysis[END_REF]. The main criteria are (1) the speed of convergence to 0 of these functions when the time or the frequency goes to infinity, which quantifies both time and frequency localizations, (2) the symmetry, (3) the number of vanishing moments of and (4) the regularity, which is useful for getting nice features, like smoothness of the reconstructed signal. The most commonly used wavelets are the orthogonal ones (Daubechies, Symlet or Coiflet wavelets) [START_REF] Addison | The illustrated wavelet transform handbook[END_REF]. Because the Daubechies wavelets [START_REF] Daubechies | Ten lectures on wavelets[END_REF] (Figure 6) have the highest number of vanishing moments, this family has been chosen for carrying out the wavelet-based multiresolution analysis of the considered sequences of data points.

The multi-layer Perceptron neural network 3.3.1. Network topology

The Perceptron, the simplest neural network, is only able to classify data into two classes [START_REF] Rosenblatt | The perceptron: a perceiving and recognizing automaton[END_REF]. Basically it consists of a single neuron with a number of adjustable weights [START_REF] Mccullogh | A logical calculus of the ideas immanent in nervous activity[END_REF]. It uses an adaptative learning rule. Given a problem which calls for more than two classes, several Perceptrons can be combined: the simplest form of a layered network just has an input layer of source nodes that connect to an output layer of neurons. The single-layer Perceptron can only classify linearly separable problems. For non-separable problems it is necessary to use more layers. A multi-layer network has one or more hidden layers whose neurons are called hidden neurons. The network is fully connected, every node in a layer is connected to all nodes in the next layer. According to the previous remarks, the network used for the present work is a multilayer Perceptron. It consists of one layer of linear output neurons and one hidden layer of nonlinear neurons [START_REF] Hornik | Multi-layer feedforward networks are universal approximation[END_REF] (Figure 7). According to previous tests, more than one hidden layer proved to cause slower convergence during the learning phase because intermediate neurons not directly connected to output neurons learn very slowly. Based on the principle of generalization versus convergence, both number of hidden neurons and iterations completed during the training phase were optimized [START_REF] Grieu | Prediction of parameters characterizing the state of a pollution removal biologic process[END_REF]. The multi-layer Perceptron neural network learns using an algorithm called backpropagation. During this iterative process, input data are repeatedly presented to the network. With each presentation, the network output is compared to the desired output and an error is computed. This error is then fed back to the network and used to adjust the weights such that it decreases with each iteration and the model gets closer and closer to produce the desired output [START_REF] Charalambous | Conjugate gradient algorithm for efficient training of artificial neural networks[END_REF].

The Levenberg-Marquardt algorithm

Several training methods were used, but the Levenberg-Marquardt algorithm [START_REF] Hagan | Training feedforward networks with the Marquardt algorithm[END_REF] proved to be the fastest and the most robust. It is particularly adapted for networks of moderate size and has memory reduction feature for use when the training set is large. Like the quasi-Newton methods, the Levenberg-Marquardt algorithm was designed to approach second-order training speed without having to compute the Hessian matrix. When the performance function has the form of a sum of squares, then the Hessian matrix can be approximated as:

The gradient can be computed as:

where is the Jacobian matrix that contains first derivatives of the network errors with respect to the weights and biases, and is a vector of network errors. The Jacobian matrix can be computed through a standard backpropagation technique that is much less complex than computing the Hessian matrix. The Levenberg-Marquardt algorithm uses this approximation to the Hessian matrix in the following Newton-like update:

When the scalar is zero, this is just Newton's method, using the approximate Hessian matrix. When μ is large, this becomes gradient descent with a small step size. Newton's method is faster and more accurate near an error minimum, so the aim is to shift towards Newton's method as quickly as possible. Thus, is decreased after each successful step and is increased only when a tentative step would increase the performance function. In this way, the performance function will always be reduced at each iteration of the algorithm. The main drawback of the Levenberg-Marquardt algorithm is that it requires the storage of some matrices that can be quite large for certain problems. The size of the Jacobian matrix is , where is the number of training sets and is the number of weights and biases in the network. It turns out that this matrix does not have to be computed and stored as a whole. For example, if we were to divide the Jacobian into two equal submatrices we could compute the approximate Hessian matrix as follows:

Therefore, the full Jacobian does not have to exist at one time. The approximate Hessian can be computed by summing a series of subterms. Once one subterm has been computed, the corresponding submatrix of the Jacobian can be cleared.

(7) (6)
(5) (4)

Forecast methodology

This section deals with the description of the proposed forecast methodology, based on both the concept of time series and the use of a wavelet-based multi-resolution analysis and multi-layer artificial neural networks. One speaks of MRA-ANN methodology.

Overall approach

A time series is a sequence of data points, measured typically at successive times and spaced at, often uniform, time intervals [START_REF] Brockwell | Time series: theory and methods[END_REF]. Time series forecasting is the use of a model to forecast future events based on known past events [START_REF] Brockwell | Introduction to time series and forecasting[END_REF]. The aim of the proposed MRA-ANN methodology is forecasting, thanks to artificial neural networks, a sequence of data points, using past sequences, also of length . According to the respective values of , and of the sampling time (set initially to 30 minutes; forecasted sequences are then upsampled to be in agreement with the sampling time of all the parameters measured at the district boiler), the results can be more or less accurate. Indeed, the further the forecasting horizon (set to 4h30, as requested by the boiler operators and taking into account the main objective of the OptiEnR research project, this leads to

), the more inaccurate it can be. Inaccuracy can also be the result of a lack of usable information provided by past sequences or the result of a too low number of considered past sequences. Past sequences are examples to be learned and used to forecast future values of the considered parameters (in this case, the outdoor temperature and the thermal power consumption of the hot water distribution network). To help the model to place itself in time, past sequences are completed, for each of their components, by the minute of the day and the day of the year. The proposed neuronal model has been developed (artificial neural networks have been trained) using half of the available sequences and validated thanks to the remaining sequences (Figure 8).

Wavelet-based multi-resolution analysis contribution

Since the data usually consist of low and high frequencies components including noise, obtaining the essential frequency-domain features from these data becomes an important aspect of data analysis, visualization and forecasting. As a consequence, understanding and being able to learn from data, using artificial neural networks, is not as easy as one may think. Usually, tools such as artificial neural networks learn very well, optimizing their topologies and all the training parameters, from limited frequency domain data. When the data consist of low and high frequencies components, a useful approach can be to develop one forecast model for each of the frequency domains and then to combine the information provided by all of the models to build the desired signal. That is why one needs a tool allowing, on one hand, decomposing a signal according to different scales of frequency, while preserving its temporal characteristics (which is fundamental when proposing a forecast methodology based on time series), and, on the other hand, rebuilding this signal. This tool is the wavelet-based multi-resolution analysis. It allows decomposing, according to both a decomposition level ( ) and a wavelet order ( ), a signal into approximations and details, using a filter bank composed of low-pass and high-pass filters. Keeping in mind the previously-mentioned characteristics of both the outdoor temperature and the thermal power consumption of the hot water distribution network, one can easily understand that the waveletbased multi-resolution analysis allows isolating the overall trend and the 24-hour pseudo period characterizing the considered time series (approximation coefficients) from the variability caused by climatic phenomena (detail coefficients). As previously-mentioned, the wavelet coefficients are of lower variability than the original and to-be-forecasted time series, which favors the estimation of future values. All the series of past sequences (of length ) considered to develop the model are so decomposed into details and an approximation of level . For each series, the target sequence is also decomposed in the same way (Figure 9).

Multi-layer Perceptrons for forecasting details and approximations of future sequences

Substituting the prediction task of an original time series of high variability by the prediction of its detail and approximation coefficients on different levels of lower variability is the key-point of the proposed MRA-ANN methodology. In this sense, a multi-layer Perceptron neural network is needed, and it has of course to be trained using the coefficients of the considered past sequences, for estimating each of the details and the approximation of level of the target sequence (the sequence to be forecasted). Let us remember that a sequence is composed of data points as well as temporal indicators (for each observation, the minute of the day and the day of the year). So, for a decomposition of level , the overall forecast model is composed of multi-layer Perceptrons. As previously mentioned, summing up the estimated coefficients allows rebuilding the to-be-forecasted sequence (Figure 10). Similar approaches have already been applied to the forecasting of parameters such as electric consumption [START_REF] Tran | Wavelets decomposition and Box-Jenkins methodology or neural networks for short-term electric consumptions forecast[END_REF][START_REF] Nengling | Techniques of applying wavelet transform into combined model for short-term load forecasting[END_REF], solar irradiation [START_REF] Grieu | Transformée en ondelettes discrète et réseaux de neurones artificiels pour la prédiction d'irradiation solaire à court terme, Energies renouvelables et éco-conception en génie électrique[END_REF][START_REF] Tran | Forecasting of solar irradiation using wavelets decomposition and cascade-correlation neural networks[END_REF] or wind speed [START_REF] Tran | Forecasting of wind speed using wavelets analysis and cascade-correlation networks[END_REF]. The particularities of the proposed methodology lies in the use of sequences (to train the artificial networks or to be forecasted), in the addition of temporal information allowing a better understanding of how the considered parameters ( and ) evolve in time, and in the use of a specific artificial neural network for estimating one of the wavelet coefficients of the future values we want to forecast. Of course, and this is the main weakness of such an approach, we have to train and find the right topology for each of the neural networks used.

Figure 10. Multi-layer Perceptrons for forecasting details and approximations of future sequences.

Rebuilding of these sequences (MRA-ANN methodology).

Selection of the most accurate models

As previously mentioned, the respective values of (the number of considered past sequences), (the length of the sequences), (the sampling time), (the wavelet decomposition level) and (the wavelet order) impact the forecast accuracy. The networks topology, basically the number of hidden neurons, as well as training parameters, such as the number of iterations, the learning rate and the error goal, also impact the networks generalization capability. That is why the present paper deals with the study and the quantification of all of these influences. It focuses also on the impact on the model performance of the random initialization of all the networks synaptic weights and on the relevance of developing "specialized" models (each model is developed to forecast outdoor temperature during only a specific period of the year). As it will be highlighted in the results section and because the Mean Relative Error (MRE) distribution (considering the networks validation phase) can be considered as a Gaussian function, with more or less skewness, one can easily select the most accurate models. According to statistical considerations, a model is selected, for a given configuration related with all the above-mentioned parameters and with a random initialization of all the network synaptic weights, when the validation error is inferior to the mean error (considering all the tested configurations) minus twice the deviation standard. The selected models define the overall forecast model.

Outdoor temperature forecasting

This section of the paper presents experimental results in forecasting outdoor temperature. All the parameters mentioned in section 3.4.4 have been optimized to obtain the most accurate models. With the aim of reducing the complexity of this optimization problem, all the considered parameters have not been studied at the same time. As a consequence, the models configurations obtained are probably not optimal configurations but, most certainly, sub-optimal configurations. According to the common topology (the number of hidden neurons) of the multi-layer Perceptrons used, the influence of both the wavelet order and the wavelet decomposition level is, first, studied (section 4.1). Next, the impact of the number of both the considered past sequences and specialized models (section 4.2) as well as of the sampling time (section 4.3), considering optimal values of , and , is studied. As mentioned in section 3.4.1, the forecasting horizon has been set to 4h30. It is a good compromise between accuracy and computation time, keeping in mind that forecasted sequences will be used by a model predictive controller to manage a thermal storage unit.

Influence of both the wavelet order and the wavelet decomposition level

As just above-mentioned, the two first parameters studied are both the wavelet order and the wavelet decomposition level. Only the number of hidden neurons varies, the sampling time is set to 30 minutes while only one past sequence of 9 points is considered ( and ). Figure 11 and 12 present the validation MRE according to , and (the networks common number of hidden neurons), ranging respectively between 2 and 10, between 1 and 9 and between 2 and 50. Let us note that the wavelet decomposition level has been averaged over its variation interval to obtain Figure 10. It is interesting to note that the best results (the most accurate forecasts) are achieved for a wavelet order equal to 4, using 5 hidden neurons (Figure 11). This can be explained by both the specific shape and frequency spectrum of a wavelet. This spectrum can be found, or not, in the analyzed signal. Considering, for example, a two-order wavelet, its shape is too different from the signal form and this leads to bad results. Having a quick look at Figure 12 (the wavelet order has been averaged over its variation interval), it can be noticed that when the wavelet decomposition level is high, only few hidden neurons are needed for accurately forecasting outdoor temperature. More hidden neurons are needed when the decomposition level is low. One can consider, and this is a very interesting result, that the number of hidden neurons to be used for accurately forecasting outdoor temperature is inversely proportional to the decomposition level. Indeed, for a too low decomposition level, each coefficient contains an important part of the initial information and, as a consequence, many hidden neurons are required to correctly model these coefficients using artificial neural networks. For a high decomposition level, the detail and approximation coefficients contain a small amount of information; that is why both their modeling and forecasting are easier and require, for each artificial neural network used, fewer hidden neurons. On the other side, let us remember that the higher the decomposition level, the more the networks needed (one for each coefficient to be estimated). The most accurate forecasts are achieved for a 5-level decomposition and 5 hidden neurons.

Figure 12. Influence of the wavelet decomposition level according to the networks topology (the number of hidden neurons).

Figure 13 depicts the linked influence of the wavelet order and the wavelet decomposition level (the common number of hidden neurons for all of the neural networks used has been averaged over its variation interval). It highlights specific domains of low MRE, one can speak of "valleys", and specific domains of higher MRE. Taking a look at Figure 13, it becomes apparent that, beyond an adequate wavelet decomposition level, the wavelet order is more influential than the decomposition level itself. Two valleys are clearly visible when the wavelet order is set to 4 or when it is set to 7 (the best choice seems to be 4). These valleys are very slightly dependent on the decomposition level.

Figure 13. Linked influence of the wavelet order and the wavelet decomposition level.

As a result of this first part of the study, one can conclude that the most accurate forecasts are achieved when decomposing the considered sequences into 5 levels, using Daubechies wavelets of order 4 and multi-layer Perceptrons which respective hidden layers are composed of 5 neurons. These results were taken into account to study the influence of the number of considered past sequences. Let us note that similar results were obtained with a sampling time set to 10 minutes.

Influence of both the number of considered past sequences and the number of specialized models

In this section, the influence of both the number of considered past sequences and the number of specialized models is analyzed. Several simulations have been carried out, considering a number of past sequences, a number of specialized models and a number of hidden neurons ranging respectively between 1 and 8, 1 and 3 and 3 and 21. Figure 14 depicts the influence of the number of considered pas sequences, considering only one overall model, according to the common topology of the neural networks used to model all of the wavelet decomposition coefficients. The results highlight the necessity of considering at least 4 past sequences (the optimal number) to obtain good forecasts. This clearly reduces the validation MRE observed. However, for a too high number of considered past sequences, the forecasts accuracy deteriorates gradually. Again, the common optimal number of hidden neurons for all of the multi-layer Perceptrons used is equal to 5. With this topology, considering 4 past sequences instead of 1 for developing and validating the overall model, the MRE is reduced of about 12.5% (from 4.8% to 4.2%).

Figure 14. Influence of the number of considered past sequences according to the networks topology (one overall model).

Figure 15. Influence of the number of considered past sequences according to the networks topology (three specialized models).

Figure 15 depicts the influence of the number of considered past sequences, but now considering three specialized models (each model is developed to forecast outdoor temperature during a specific period of the year only), according to the common topology of the neural networks used to model all of the wavelet decomposition coefficients. Each of the three specialized models being trained using less (three times less) information than when developing and validating an overall model, less hidden neurons are needed for accurately estimating, using multi-layer neural networks, approximation and detail coefficients. As a consequence, the influence of a bad networks sizing (too many hidden neurons favors overfitting) is more visible.

Considering the right number of hidden neurons for all of the neural networks used to model the wavelet decomposition coefficients, Figure 16 depicts the influence of the number of considered past sequences according to the number of specialized models. Analyzing this figure, one can easily conclude that using several specialized models, instead of a unique overall model, does not improve the forecasts accuracy, quite the opposite… Whereas the main objective of developing specialized models was to dissociate and better understand (learn) the various behaviors outdoor temperature might have, one can notice that both the too restrictive historical data considered (from mid January to early April) and the networks training using less examples do not allow reducing the MRE. The networks generalization capability seems to be negatively impacted. So, one can suppose that considering, for example, one year or historical data, developing specialized models according to the various seasons of the year (temperature evolution is quite different from a season to another) would be very useful. Finally, the influence of the sampling time on the models performance has also been studied, considering a wavelet order and a wavelet decomposition level equal to 4 and 5 respectively. One can conclude that reducing the sampling time from 30 minutes to 10 minutes does not impact significantly the forecasts accuracy. Figures 17 and18 depict, with a sampling time set to 10 minutes, the influence of the number of considered past sequences according to the topology of the neural networks used and the number of specialized models (in this case, the right number of hidden neurons for all of the neural networks used to model the wavelet decomposition coefficients is used) respectively. Both figures highlight that four past sequences and a unique overall model lead to the best results (in this case, the MRE observed is about 4.3%). A reduced sampling time impacts significantly on the common topology of the neural networks used to model all of the wavelet decomposition coefficients only. The common number of hidden neurons is about 15 (it was only of 5 with a sampling time set to 30 minutes). One can suppose that, because a reduced sampling time increases the length of the sequences (in this case, up to 27 points), more information has to be learned and, as a consequence, requires more hidden neurons for all of the neural networks used.

Optimal configuration and experimental results

As a conclusion of the study about outdoor temperature forecasting, Table 3 depicts the optimal configuration to meet performance needs when applying the MRA-ANN methodology, according to the wavelet order ( ), the wavelet decomposition level ( ), the topology (the common number of hidden neurons ) of the artificial neural networks used for estimating the detail and approximation coefficients, the number of considered past sequences ( ), the number of specialized models ( ) and, finally, the sampling time ( ). ), forecasted at the meteorological station ( ) and forecasted thanks to the MRA-ANN methodology based on a multi-resolution analysis and the use of artificial neural networks (

). Clearly, the proposed methodology provides better forecasts than the meteorological station. The temperature peaks are far better modeled (as previously mentioned, it was one of the major objectives of the work, the forecasted temperature being used by a MPC for optimizing the boiler functioning). Indeed, one can note that frequent differences of about 5°C between the measured and forecasted (at the meteorological station) outdoor temperatures are reduced, using the developed model, to only 1°C. Table 4 specifies, for the considered period (from February 15, 2009 to March 02, 2009) and the two forecasted outdoor temperatures (by the meteorological station and using the MRA-ANN methodology), the Mean Relative Error (MRE), the Mean Absolute Error (MAE) as well as the curve fitting (FIT) obtained. First, one can remark, considering the temperature forecasted at the meteorological station as a reference, that the developed model allows reducing both the MRE and the MAE by 33% (from 6.28% to 4.14%) and 29.4% (from 1.75°C to 1.15°C) respectively. The curve fitting is improved of about 40% (from 42.6% to 60.6%). As a conclusion, one can highlight the model performance and the validity of the proposed MRA-ANN methodology. ), forecasted at the meteorological station ( ) and forecasted using the MRA-ANN methodology ( ).

Table 4. Forecasting performance criteria (outdoor temperature).

Forecasted outdoor temperature MRE [%] MAE [°C] FIT [%]

6.28 1. 75 42.6 4.14 1.15 60.6

Thermal power consumption forecasting

This section of the paper presents experimental results in forecasting thermal power consumption. Two different ways of doing were considered: first, using the previously forecasted outdoor temperature (both parameters are strongly correlated) and, secondly, as we did for outdoor temperature, applying the MRA-ANN methodology. Finally, a hybrid approach is proposed.

Use of the forecasted outdoor temperature

As previously mentioned (section 3.1), outdoor temperature and thermal power consumption are inversely proportional. That is why the first of the two proposed approaches deals with a simple linear curve fitting, optimizing both the parameters and of Equation 8 thanks to the minimization of the difference between (the calculated thermal power consumption) and (the thermal power consumed by the hot water distribution network and estimated using the outdoor temperature forecasted at the meteorological station (

)) or (the thermal power consumed by the hot water distribution network and estimated using the outdoor temperature forecasted thanks to the MRA-ANN methodology (

)). Let us note that and if is used while and with .

(

) 8
Table 5 (lines 1 to 2) and Figure 20 present the results of the thermal power consumption ( ) forecasting. These results validate the first proposed approach: one can easily predict the thermal power consumption of the hot water distribution network using outdoor temperature. Whatever the forecasted temperature used, the results are quite similar. For the previously-mentioned period (from February 15, 2009 to March 02, 2009), the MRE, the MAE and the FIT are about 6.97%, 761.9 kW and 37.6% when using and about 6.94%, 758.3 kW and 38.4% with . However, having a quick look at Figure 20, one can remark, first, that the amplitude is not always well estimated and, secondly, that the curves are not always completely in phase. So, it could be interesting to use the proposed MRA-ANN methodology for directly forecasting the thermal power consumption of the hot water distribution network.

Table 5. Forecasting performance criteria (thermal power consumption).

Thermal power consumption ( ) forecasted using: MRE [%] MAE [kW] FIT [%]

6.97 761.9 37.6 ) and the hybrid approach ( ). Comparison with the calculated .

The MRA-ANN methodology proposed to forecast outdoor temperature, based on a multiresolution analysis and artificial neural networks, has also been used for directly estimating the thermal power consumption of the hot water distribution network (

). Due to the high correlation between the just-mentioned parameters, the optimal configuration found for outdoor temperature forecasting (Table 3) has been re-used to develop and validate the new model. Table 5 (line 3) and Figure 21 present the results of the thermal power consumption forecasting, using the MRA-ANN methodology. The MRE, the MAE and the FIT are about 6.07%, 663.6 kW and 44.7% respectively. These results are better than when dealing with a linear curve fitting, using the forecasted outdoor temperature ( ): one can highlight that the model allows reducing both the MRE and the MAE by 12.5%, from 6.94% to 6.07% and from 758.3 kW to 663.6 kW respectively. The curve fitting is improved of about 16.4% (from 38.4% to 44.7%). However, taking a look at Figure 21, one can remark that the peaks of consumption are again not always well estimated. That is why a last approach is proposed, combining time series forecasting (the MRA-ANN methodology) and linear curve fitting.

Hybrid approach

As just mentioned, a last approach is proposed, combining the MRA-ANN methodology and linear curve fitting, for predicting the thermal power consumption of the hot water distribution network. One can highlight that time series forecasting (section 5.2) provides better results than linear curve fitting (using the forecasted outdoor temperature, section 5.1) but sometimes underestimates peaks of consumption. Both approaches are valid but have different characteristics. That is why it appears that they can be used in tandem to improve the forecasts accuracy. The idea is so to combine the forecasted sequences provided by each of the two considered approaches ( and ) using a weighted mean whose coefficients and are optimized in the following way (equation 9):

The optimal values of and are about 0.28 and 0.71 respectively. As it was expected, the weight of the most accurate forecasted sequence, provided by the MRA-ANN methodology, is the highest weight. However, one can note that the influence (weight) of the sequence estimated by the "linear curve fitting" method, using the forecasted outdoor temperature, is regardless significant. Table 5 (line 4) and Figure 21 present the results of the thermal power consumption forecasting, using the just-described hybrid approach. The MRE, the MAE and the FIT are about 5.85%, 639.1 kW and 46.9% respectively. Taking as a reference the performance of the MRA-ANN methodology, one can observe that the proposed hybrid approach allows reducing both the MRE and the MAE by 3.7%, from 6.07% to 5.85% and from 663.6 kW to 639.1 kW respectively. The curve fitting is also improved of about 5% (from 44.7% to 46.9%). Although these improvements do not appear to be significant, they validate the last proposed hybrid approach for estimating the thermal power consumption of the hot water distribution network ( ).

Conclusion and perspectives

As part of the OptiEnR project, the present paper, the first of a series of three dealing with the modeling and the optimization of a district boiler, situated at La Rochelle (west coast of France), focuses on forecasting both the outdoor temperature and the thermal power consumption of the hot water distribution network. This district boiler is composed of a wood boiler, a gas-fuel oil boiler, a breaking pressure bottle and a cogeneration plant. With the aim of optimizing the use of a hot water thermal storage unit and reducing the gas-fuel oil boiler coverage rate (i.e. the fossil energy consumption), a model predictive controller will be developed and implemented. It requires forecasting, over the next 4h30, the two above-mentioned parameters.

The proposed MRA-ANN methodology deals with the concept of time series and uses a waveletbased multi-resolution analysis and artificial neural networks. Multi-layer Perceptrons were trained and validated, using the wavelet decomposition coefficients of past sequences, for estimating both the detail and approximation coefficients of the considered target sequences. Substituting the prediction task of an original time series of high variability by the estimation, using artificial [START_REF] Khotanzad | An artificial neural network hourly temperature forecaster with applications in load forecasting[END_REF] neural networks, of its wavelet coefficients on different levels of lower variability (according to different frequency domains), while preserving the temporal characteristics of the series, is the main idea of the present work. Then, the reconstruction of the target sequences is performed by simply summing up the estimated coefficients. One can note that the wavelet-based multi-resolution analysis allows isolating both the overall trend and the 24-hour pseudo period characterizing the considered time series (approximation coefficients) from the variability caused by climatic phenomena (detail coefficients). Because these coefficients are of lower variability than the original and to-beforecasted time series, the MRA-ANN methodology makes the estimation of future values easier.

The number of considered past sequences ( ), the sampling time ( ), the wavelet decomposition level ( ) and the wavelet order ( ) impact the forecast accuracy. The networks topology, basically the number of hidden neurons ( ), as well as training parameters, such as the number of iterations, the learning rate and the error goal, also impacts the networks generalization capability. That is why the present paper deals with the study and the quantification of all of these influences. It focuses also on the impact of the random initialization of all the networks synaptic weights on the performance of the models and on the suitability of developing specialized models ( ). As a result of the study about outdoor temperature forecasting, the optimal configuration, leading to the most accurate forecasts (the MRE, the MAE and the FIT are about 4.2%, 1.2°C and 60% respectively), is defined as follows: , , , , and . These results validate the proposed MRA-ANN methodology and highlight both the impact of working with wavelet decomposed sequences and the generalization capability of multi-layer neural networks.

Next, two different ways of doing have been considered for forecasting the thermal power consumption of the hot water distribution network: first, using the previously forecasted outdoor temperature (both parameters are strongly correlated) and, secondly, applying the MRA-ANN methodology used to forecast outdoor temperature again (in this case, the optimal configuration highlighted during the first part of the work has been re-used to develop and validate the new model). Finally, a hybrid approach is proposed, combining time series forecasting and linear curve fitting. Whatever the considered approach, the results are satisfactory. The MRE, the MAE and the FIT are ranging between 5.85% and 6.97%, 639.1 kW and 761.9 kW and 37.6% and 46.9% respectively. The hybrid approach provides the best results. However, one can note that, comparing the two forecasted parameters, forecasts accuracy is not as good as that achieved when predicting outdoor temperature. It can be explained by hardly-quantifiable influences leading to amplitude variations. Some peaks of consumption may be due to an unusual functioning of big buildings heating systems (for example, an exhibition hall connected to the hot water distribution network). Some other peaks may be the consequence of the way the thermal power consumption is specified, using both the temperatures of the water entering the distribution network and coming back from this network. So, if the gas-fuel oil boiler starts working, the temperature of the water entering the network will increase quickly but no significant impact will be noticed on the temperature of the cold water coming back to the plant until a few minutes. However, this leads to variations of the thermal power consumption, induced by both the functioning of the district boiler and the thermal inertia of the hot water distribution network. One can suppose that adding to the plant a hot water thermal storage unit (one of the next steps listed in the OptiEnR research project) will limit the presence of hardly-predictable peaks of consumption.

As just mentioned, future work will now focus on optimizing the performance of the district boiler of La Rochelle, adding to the plant a thermal storage unit which will be managed thanks to a model-based predictive controller. Using forecasted outdoor temperature and thermal power consumption sequences, such a controller calculates an optimal command sequence to be applied to the thermal storage feed pump. This allows adjusting the flow of the water passing through this unit. Let us remember that the main objective of the work deals with the minimization of the fossil energy consumption, stocking renewable energy during low-demand periods and using it during high-demand periods. Future papers will describe both the modeling and the way the district boiler functioning will be optimized.

Figure 1 .

 1 Figure 1. Diagram of the OptiEnR project. In red, the forecast module presented in the paper.

Figure 2 .

 2 Figure 2. Synopsis of the district boiler of La Rochelle. Its current situation (in black) and future perspectives (in red, blue and green).

Figure 3 .

 3 Figure 3. The outdoor temperatures forecasted at the meteorological station ( ) and measured at the district boiler ( ) (from early January to late March).

Figure 4

 4 Figure 4 depicts the calculated thermal power consumption of the hot water distribution network () during the first months of 2009 (from early January to late March). This parameter is not measured in situ and has an exogenous influence on the district boiler functioning because it impacts on the temperature of the water back from the distribution network.is calculated using the difference between (the temperature of the hot water leaving the breaking pressure bottle and entering the distribution network) and(the temperature of the cold water coming back from the distribution network) as well as the flow of the water passing through this network (). is estimated using the opening percentage of the distribution network feed pump valve ( ), related to the water differential pressure ( ), its set-point ( ) and outdoor temperature (Figure2). As depicted by Figure4, the thermal power consumption and the outdoor temperature are inversely proportional.

Figure 5 .

 5 Figure 5. Multi resolution analysis leading to the 3-level decomposition of a signal x.Because translating the mother wavelet leads to redundant information, one can use, instead of the continuous wavelet transform, the discrete one. In this case, both scale ( ) and translation ( ) factors are restricted to only discrete values such that and , where , , and. Basically, the discrete wavelet multi-resolution analysis, commonly based on Daubechies orthogonal wavelet basis[START_REF] Daubechies | Ten lectures on wavelets[END_REF], allows decomposing a signal into approximations (i.e. low frequency coefficients) and details (i.e. high frequency coefficients) using a filter bank composed of both low-pass ( ) and high-pass ( ) filters[START_REF] Mallat | A wavelet tour of signal processing[END_REF]. This process can be repeated times, producing levels of decomposition, but decomposing (downsampling) the approximations

Figure 6 .

 6 Figure 6. The Daubechies wavelets family, according to the wavelets order. In black: the wavelets, in grey: the associated scaling functions.

Figure 7 .

 7 Figure 7. Topology of the multi-layer Perceptrons used for estimating the wavelet coefficients of forecasted sequences: one hidden layer of nonlinear neurons and an output layer of linear neurons.

Figure 8 .

 8 Figure 8. Development phase of the forecast neuronal model.

Figure 9 .

 9 Figure 9. Decomposition of past and target sequences into details and an approximation of level .

Figure 11 .

 11 Figure 11. Influence of the wavelet order according to the networks topology.

Figure 16 .Figure 17 .

 1617 Figure[START_REF] Hippert | Combining neural networks and ARIMA models for hourly temperature forecasts[END_REF]. Influence of the number of considered past sequences according to the number of specialized models.4.3. Influence of the sampling time

Figure 18 .

 18 Figure 18. Influence of the number of considered past sequences according to the number of specialized models (sampling time = 10 minutes).

Figure 19

 19 Figure 19 allows comparing, from February 15, 2009 to March 02, 2009, the outdoor temperatures measured at the district boiler (), forecasted at the meteorological station ( ) and forecasted thanks to the MRA-ANN methodology based on a multi-resolution analysis and the use of artificial neural networks (). Clearly, the proposed methodology provides better forecasts than the meteorological station. The temperature peaks are far better modeled (as previously mentioned, it was one of the major objectives of the work, the forecasted temperature being used by a MPC for optimizing the boiler functioning). Indeed, one can note that frequent differences of about 5°C between the measured and forecasted (at the meteorological station) outdoor temperatures are reduced, using the developed model, to only 1°C.

Figure 19 .

 19 Figure 19. Outdoor temperature measured at the district boiler (), forecasted at the meteorological station ( ) and forecasted using the MRA-ANN methodology ( ).

Figure 20 . 2 .Figure 21 .

 20221 Figure 20. Thermal power consumption forecasted using ( ) or ( ). Comparison with the calculated .5.2. Use of the MRA-ANN methodology

Table 1 .

 1 Coverage rates of the three heat generators.

	Heat generator	Coverage rate
	Wood boiler	50%
	Gas-fuel boiler	15 to 20%
	Cogeneration plant	30 to 35%

Table 3 .

 3 Optimal configuration for outdoor temperature forecasting (MRA-ANN methodology).

	Parameter	Symbol	Optimal value
	Wavelet order		4
	Wavelet decomposition level		5
	Common number of hidden neurons		5
	Number of considered past sequences		4
	Number of specialized models		1
	Sampling time (minutes)		30


Ismael Mejía

Mario Südholt

Luis Daniel

Benavides Navarro
email: lbenavidesnavarro@acm.org

Invasive composition for the evolution of a health information system

Keywords: D.3.3 [Programming Languages]: Language Constructs and Features-Frameworks, Patterns Design, Languages Aspect-oriented programming, Distributed systems, Health information systems, Invasive software composition

In this paper we show that some of the evolution tasks in OpenMRS, a health information system, may require the invasive modification of interfaces and implementations in order to offer an appropriate modularization. We introduce a new composition framework in Java that supports the definition of expressive pattern-based invasive compositions. Furthermore, we show that the composition framework allows us to concisely define an evolution scenario of OpenMRS that supports the consolidation of patient data from different remote instances.

INTRODUCTION

Evolution of large-scale distributed systems is a fundamental challenge of current information systems, notably web-based ones. Two problems are particularly difficult to handle:

• The modification of heterogeneous distribution and communication requirements, for example, if previously colocated parts of an application are to be executed in a remote location.

• The need for invasive modifications [START_REF] Aßmann | Invasive Software Composition[END_REF], i.e., modifications to implementations, e.g., to enrich composition interfaces by exposing previously hidden behaviors.

We have studied these problems in the context of an information system for the health domain, the OpenMRS system [START_REF]OpenMRS home[END_REF], an open-source web-based medical record management system 1 . In this paper we consider an evolution of this system that enables the consolidation of information of a patient's medical record history from different sources. Such consolidation is needed, for instance, if different instances of OpenMRS are used in different locations. This system and evolution scenario is subject to the two problems mentioned above: the consolidation relies on data gathered from different sources (i.e., requires modification of communication requirements) and requires modifications to the handling of a patient's medical history that is not expressible in terms of the OpenMRS's corresponding interfaces but need (limited) access to implementations (requires invasive modifications to enrich composition interfaces).

Traditionally, such evolutions are implemented "by hand" without tool support for the required structural and behavioral modifications. Existing tools are limited to the generation of parts of the implementation from higher-level specifications, e.g., using model-driven engineering, or low-level manipulation of implementations using refactoring techniques. Frequently there is no higher-level specification to derive implementations from and evolutions are only expressible in terms of many non-trivial and low-level program refactorings. The evolution scenarios of OpenMRS that we consider are of this kind. Furthermore, design patterns and implementation patterns that provide, in principle, a general solution to such kinds of evolution scenarios cannot be applied easily to distributed applications that have heterogeneous communication and computation requirements (in contrast to more regular or 'simpler' systems, such as massively parallel applications [START_REF] Bromling | Pattern-based parallel programming[END_REF], workflow definitions [START_REF] Van Der Aalst | Workflow patterns[END_REF], or system integration in terms of message manipulations [START_REF] Hohpe | Enterprise Integration Patterns: Designing, building, and deploying messaging solutions[END_REF]).

Our approach to support heterogeneous evolutions scenarios builds on pattern-like solutions like algorithm skeletons [START_REF] Cole | Algorithmic skeletons: structured management of parallel computation[END_REF] for parallel algorithms. The corresponding patterns include farming out the same computation on many code or a gather pattern that blocks the execution of a task until it has received a set of responses from different nodes. Invasive distributed patterns [START_REF] Daniel | Invasive patterns for distributed applications[END_REF] extend algorithmic skeletons by the In this paper, we sketch an extension to the approach of invasive patterns [START_REF] Daniel | Invasive patterns for distributed applications[END_REF][START_REF] Mejía | Structured and flexible gray-box composition: application to task rescheduling for grid benchmarking[END_REF] that allows evolutions to be defined in terms of a language for, potentially invasive, composition. Concretely, we present the following contributions and issues for discussion:

• We briefly introduce a first design and implementation in Java of a composition framework for invasive patterns that instantiates the abstract language defined in [START_REF] Mejía | Structured and flexible gray-box composition: application to task rescheduling for grid benchmarking[END_REF].

• We motivate an evolution scenario of OpenMRS in which invasive composition is appropriate. We show how to implement it using our the concept of invasive patterns with our composition framework.

CONSOLIDATION OF ENCOUNTER IN-FORMATION IN OpenMRS

We have investigated evolution problems in the context of health information systems (HIS). Non-anticipated evolution tasks occur frequently in HIS when new technologies or applications need to be integrated, as well as when new business needs appear (e.g., new treatments, new legislation rules or new administrative processes).

We are investigating such evolution tasks in the context of the OpenMRS HIS [START_REF]OpenMRS home[END_REF], a well known HIS created to track medical records (patients, visits, diseases, etc) that has been deployed and used in several communities, mostly in developing countries. We present results related to one evolution task in this paper: enriching medical encounter information of patients with information concerning that patient from different sources, i.e., remote instances of OpenMRS where the patient has been treated. Such consolidation is useful e.g., to implement strategies to identify and handle epidemic situations. Technically, this task requires the addition of functionality for the management of distribution between different OpenMRS instances and consolidating patient information in the resulting distributed system (in the remainder, we denote this task by ConsolPatData). classes that represent the main functionalities of the system: patient and location management, consultations (aka medical encounters) with their respective diagnostic findings (aka observations) and orders, as well as administrative tasks like the creation of users, concepts (e.g., new drugs) and treatments. The part of the implementation relevant to our study can be roughly structured in two parts (we do not consider the data tier, i.e., how data is stored): the implementation of the web-tier and the health business functionality. The latter is implemented based on several core services, the business services, and an application-level service.

The implementation of the ConsolPatData evolution task basically requires two sets of modifications: (i) new functionality has to be inserted for the exchange of information between distributed instances of OpenMRS and (ii) the existing data handling processes have to be modified in order to correctly integrate the new information. Note that the latter requirement does not mean to only store the remote data and then handle it as local one: the fact that remote and local data co-exist means that existing procedures have to be modified, e.g., to make decisions only after the evaluation of the patient's health records at the remote and local site have been analyzed and compared.

In Fig. 1 gray stars and circles mark the parts of the Open-MRS implementation that are affected by these two sets of modifications. The stars mark invasive modifications, i.e., modifications to the implementations of the corresponding functionalities; the circle marks an interface-level modification.

The ConsolPatData evolution can be represented by the pattern composition shown in Fig. 2. Read from left to right, it defines how the information of a patient's medical record history is consolidated. First a medical center farms out (by means of pattern f arm1) an execution request to different nodes, e.g., other hospitals and clinics, in order to get the encounter information corresponding to a patient. Then some of the institutions relay requests to some of its internal departments (f arm2). As some of the departments do not have information for the patient we allow for partial compositions using a partial gather pattern (gather3). Finally, the answer of all the hosts is consolidated using the gather pattern gather4 to the requesting facility in order to construct the distributed medical record .

Implementing this pattern composition requires invasive modifications to be performed: in the initial OpenMRS in- stance patient data processes have to be modified in order to integrate and handle the remote data, in the other instances the information requested from the initial instance has to be extracted in order to be transformed. Furthermore, in all instances distribution code has to be introduced (since OpenMRS does not handle distribution natively). This composition task therefore cannot be performed solely in terms of standard patterns for distribution, such as described in [START_REF] Van Der Aalst | Workflow patterns[END_REF][START_REF] Hohpe | Enterprise Integration Patterns: Designing, building, and deploying messaging solutions[END_REF].

The necessary modifications are non-trivial and have to be applied at 16 different source code locations (and partially multiple times in different execution contexts at those locations). In previous work we have introduced invasive distributed patterns [START_REF] Daniel | Invasive patterns for distributed applications[END_REF][START_REF] Mejía | Structured and flexible gray-box composition: application to task rescheduling for grid benchmarking[END_REF] to handle such evolution tasks. That work introduced the conceptual framework and applied it using ad hoc implementations of pattern compositions. Here, we extend our previous work by presenting a first version of a Java framework for the systematic composition of invasive patterns. In the remainder of this paper, we briefly present this framework and illustrate its benefits for the evolution of OpenMRS.

A LANGUAGE FOR INVASIVE COMPO-SITION

In previous work [START_REF] Mejía | Structured and flexible gray-box composition: application to task rescheduling for grid benchmarking[END_REF] we have defined a kernel language for the composition of invasive patterns, a simplified version of which, sufficient for illustration purposes here, is shown in Fig. 3. The kernel language is based on the following key concepts: composition programs (P rog) are sequences of invasive patterns P , each of which is defined as a pair of operators O. The first operator (possibly invasively) extracts information from a set of source computations or hosts; the second applies (possibly invasive) modifications to a set of targets. Finally, operators are defined as pairs of event sequences, the first of which defining the context where the operator is to be applied, the second defining an adaptation to be applied when the operator context matches the current execution context.

The pattern-based composition program represented in Fig. 2 can be expressed with our language via the composition of different patterns characterized by event sequences as follows:

P := f arm1, f arm2, gather3, gather4
where . . . f arm2 := ( (e1@l2, ), ( , {e1,i@li} i∈{1,2,3} ) gather3 := ( ( , {e2,i@li} i∈{1,2,3} ), ( , e3@l2) ) . . . Here the composition program P is defined using a sequence of four composition operators. The operator gather3 is defined as an (outer) pair whose first component represents the source information to be gathered from the three nodes and the second pair the information to be integrated (invasively) in the target node. The sources and target are defined as (inner) pairs whose first components, both empty ( ) here, represent context information used to be generated via the events that are transmitted and integrated by means of the second components of respectively the first and second inner pair. Events define interesting points in the application e.g., e1 define a context to obtain the information of the encounters of a patient, but we are interested in a summarized view represented by the new service invocation e2 and finally consolidated in e5.

We have extended our previous work, by designing and implementing a framework for Java for invasive composition. This framework directly represents the basic concepts (composition programs, invasive patterns and operators), while providing a more powerful composition model than represented by the simplified kernel language above. Fig. 4 presents a high-level view of the main abstractions of our invasive operator library. Concretely we have modeled the main abstractions as interfaces that can be specialized. Localizable marks the elements which may belong to different locations. Class Context enables contexts to be defined by the matching of regular expressions of events (using a regular expression matcher implementing the IMatcher interface). Events are implemented using a parametric container class. Programmers may tailor contexts to their needs; we provide, however, a base implementation that is sufficient to match regular expressions over call-like pointcuts. An Adaptation defines a sort of task execution with an environment to receive and return arguments. Programmers should spe- cialize this class in order to add their own behavior if a given context is matched. Finally, an invasive Pattern is a composition of a context and an adaptation. The pattern class subsumes the operators presented in the kernel language, which provide simpler contexts and adaptations that consist solely of event sequences. The concrete implementation of the operators such as Execution, Farm and Gather, refine the semantics of the distributed characteristics of the given operations. And connectors such as: Program and Condi-tionalCompositor define the way patterns are composed: the execution of a successor pattern on a node is started as soon as the execution of the adaptation of a predecessor pattern has terminated.

OpenMRS EVOLUTION REVISITED

We have implemented the OpenMRS ConsolPatData evolution using the composition framework. Fig. 5 shows the inner part of the pattern composition. Here we show only the code of the inner composition (f arm2 → gather3) corresponding to the patterns and language that were presented previously in Fig. 2 and3. The patterns involved in this composition are nested farms, as well as partial gathers which collect information by allowing for timeouts. The composition shows in line 24, in particular, the use of composition operators over primitive invasive patterns that have been constructed by defining appropriate contexts and adaptations, as exemplified for the primitive farm pattern in lines 8-13.

Table 1 summarizes the different modifications done in order to evolve the system to support the distributed medical examination (the table only includes non deprecated methods).

The implementation of the distributed medical record requirement has been achieved with the use of our Invasive Patterns library. It allows us to program directly the involved distribution patterns and their protocol (e. case of the partial gather with timeout). We have implemented the complete solution in a concise and relatively simple way that involves one class which defines the patterns and its interaction points, and two additional utility classes that deal with concrete details of the adaptations.

RELATED WORK

Evolution and integration of distributed systems are complex tasks that require careful planning and execution. Currently, engineers execute such tasks by means of careful application of design patterns and manual modification of interfaces and implementations. Patterns for distributed applications are typically only used as informal sets of best practices for the modification of communication and interconnection structures. Furthermore, the problem of adapting interfaces to make different systems and components compatible is ordinarily also addressed by means of manual manipulation of code. Not many of the existing works proposes a more structured approach to composition. Finally, the problem of evolution and integration of distributed systems is intimately linked with the evolution of individual components. Here again, manual manipulation of code is the preferred mechanism of evolution. However, recent work on the refactoring using transformations and tool support based on Model Driven Engineering are striving for a more structured approach. In the remainder of this section we discuss related work to each of the mentioned problems of evolution and compare them with our proposal.

A fair number of approaches to the refactoring of applications using transformations have been considered. A first set of approaches addresses the problem of defining general frameworks or metamodels [START_REF] Lämmel | Towards generic refactoring[END_REF][START_REF] Reimann | Role-based generic model refactoring[END_REF][START_REF] Moha | Generic model refactorings[END_REF] and allow common refactoring actions over different programming languages and models to be expressed. These approaches extract the commonalities of refactoring actions and allow developers to apply refactoring actions on their programs and models. Other approaches have proposed modeling and transformation of source code to address refactoring actions for specific languages [START_REF] Taentzer | Specifying domain-specific refactorings for andromda based on graph transformation[END_REF][START_REF] Biermann | Graphical definition of in-place transformations in the eclipse modeling framework[END_REF], or the manipulation of source code by means of the specification of refactoring examples [START_REF] Brosch | An example is worth a thousand words: Composite operation modeling by-example[END_REF]. All these approaches address the problem of code restructuring but do not, in contrast to our approach, address the problem of refactoring complex communication patterns found in distributed systems. Furthermore, no support for complex pattern compositions is provided, as we do with our composition framework.

Concerning interface modifications, invasive Composition as proposed by Aßmann [START_REF] Aßmann | Invasive Software Composition[END_REF] adapts and extends components at hooks by means of transformations. In this approach modifications can be performed at explicit hooks that are provided by the component developer and implicit hooks as common abstractions provided by the programming platform, e.g. method entry points. Note that this approach is similar to ours in the sense that it admits invasive modification of composition interfaces, however, our approach goes a step further by allowing explicit manipulation of complex communication patterns. Furthermore, pattern compositions have to be implemented in an ad hoc manner.

Finally, other approaches use distribution patterns and parallel patterns for integration of distributed systems, and the configuration of parallel algorithms respectively [START_REF] Cole | Algorithmic skeletons: structured management of parallel computation[END_REF]. Hoppe and Woolf proposed Enterprise Integration Patterns [START_REF] Hohpe | Enterprise Integration Patterns: Designing, building, and deploying messaging solutions[END_REF], a catalog of best practices to address common problems found during integration of distributed systems using messaging middleware. Their approach explicitly addresses the manipulation of asynchronous communication, and provides several patterns for the composition of program communication behavior. However, these techniques rely on ad hoc manipulation of applications involved in the integrated system. In the domain of massively parallel applications, patterns are used to compose parallel algorithms over homogeneous deployment environments [START_REF] Bromling | Pattern-based parallel programming[END_REF] and do not support invasive composition.

CONCLUSION

In this paper we have motivated that evolution tasks in OpenMRS, a health information systems, require the invasive composition of interfaces and implementations. Concretely, we have shown that the consolidation of patient data from remote sites needs several interfaces and classes at different levels of the OpenMRS architecture to be modified. We have presented a new composition framework in Java that supports the expressive definition of pattern-based invasive compositions. Finally, we have provided evidence that our framework enables the evolution of patient data consolidation in a systematic manner and that pattern compositions significantly facilitate this evolution task.

Figure 1 :

 1 Figure 1: OpenMRS Architecture Diagram

Fig. 1

 1 shows a representation of the multi-tier architecture of OpenMRS. The ovals in the figure denote individual

Figure 2 :

 2 Figure 2: Pattern compositions for electronic medical record retrieval

Figure 3 :

 3 Figure 3: Kernel language for invasion composition

Figure 4 :

 4 Figure 4: Framework for invasive composition

public List<Encounter> getEncounters( 2 Figure 5 :

 25 Figure 5: Distributed encounter consolidation code

Table 1 :

 1 Modified elements to implement the distributed medical record

	g. in the

In this paper we have analyzed and applied an evolution task to OpenMRS, version 1.7.