Local Fusion of an Ensemble of Models for the Reconstruction of Faulty Signals

Piero Baraldi, Antonio Cammi, Francesca Mangili, and Enrico E. Zio

Abstract-Sensors are placed at various locations in a production plant to monitor the state of its components and accordingly operate its control and protection. For the plant state monitoring to be effective, the sensors themselves must be monitored for detecting

I. INTRODUCTION

S ENSORS contribute to the safe and productive operation of a nuclear power plant by conveying information on its state to the automated controls and the operators. For this to be effective, it is important to timely detect sensor malfunctions and reconstruct the incorrect signals before using them in the operation, control and protection of the plant [START_REF] Hoffmann | On-Line Monitoring for Calibration Reduction HWR-784[END_REF], [START_REF] Hoffmann | Signal grouping algorithm for an improved on-line calibration monitoring system[END_REF].

An approach to sensor signal validation relies on auto-associative models [START_REF] Holbert | An integrated signal validation for nuclear power plants[END_REF]- [START_REF] Roverso | Solutions for plant-wide on-line calibration monitoring[END_REF] e.g. based on Principal Component Analysis (PCA) [START_REF] Tharrault | Fault detection and isolation with robust principal component analysis[END_REF], AutoAssociative Kernel Regression (AAKR) [START_REF] Singer | Model-based nuclear power plant monitoring and fault detection: Theoretical foundations[END_REF] and AutoAssociative Neural Network (AANN) [START_REF] Hoffmann | Signal grouping algorithm for an improved on-line calibration monitoring system[END_REF] methods. The practical problem, however, is that a single reconstruction model cannot handle the multiplicity of signals measured on a real plant [START_REF] Zio | Genetic algoritms for grouping of signals for system monitoring and diagnostics[END_REF]- [START_REF] Fantoni | Multiple-failure signal validation in nuclear power plants using artificial neural networks[END_REF]. A possible way to overtake this limitation is to subdivide the signals into small overlapping groups, develop an ensemble of reconstruction models, one for each group, and finally combine their outcomes. Key to building of the ensemble is the diversity of the individual models. In the approach investigated in this work, diversity is promoted by randomly generating the signals groups according to the Random Feature Selection Ensemble (RFSE) technique [START_REF] Bryll | Attribute bagging: Improving accuracy of classifiers ensembles by using random feature subsets[END_REF]; this is a completely random technique in which no optimization of the composition of the individual groups is sought, i.e., no relevance is given, for example, to the correlation between the signals in the groups or to their capability of reconstruction. The groups thereby created are then used to develop a corresponding number of signal reconstruction PCA models [START_REF] Jolliffe | Principal Component Analysis[END_REF]- [START_REF] Moore | Principal component analysis in linear systems: Controllability, observability, and model reduction[END_REF].

The investigation of different methods for aggregating the outcomes of the individual models of the ensemble is the main objective of the present work. In particular, LF methods are explored, in which the aggregation is guided by the local performance of each model, i.e., its reconstruction accuracy on signal patterns of training similar (and for this reason also called neighbors) to those to be reconstructed (also called test patterns) [START_REF] Baraldi | A procedure for the reconstruction of faulty signals by means of an ensemble of regression models based on principal components analysis[END_REF]. These methods rely on the idea that each model can perform well in some regions of the signals space and poorly in others.

The general LF process is based on the following three steps: 1) retrieve neighbors of the test pattern from the set of training patterns; 2) associate a weight and a bias to each individual model of the ensemble depending on its reconstruction accuracy on the retrieved neighbors; 3) aggregate the outputs, accounting for the models weights and bias. With respect to step 1), two strategies have been compared for the identification of the neighborhood of a test pattern: the first one considers all the patterns of the training set located in a hyper-rectangle centered on the test pattern of interest; the second one includes only its nearest neighbors.

With respect to step 2), for a given test pattern the bias correction of each model in the ensemble is taken equal to the local mean error which the model makes on the neighbor patterns of training, whereas different weighting strategies have been considered: a) weight proportional to the inverse of the mean absolute error made by the model on the neighboring patterns of training; b) weight proportional to the logarithm of the inverse of the mae; c) the borda-count method [START_REF] Polikar | Ensemble based system in decision making[END_REF]. With respect to step 3), the output aggregation is performed by a weighted average with the weights computed at step 2 modifying the model outcomes by subtracting the bias.

In order to enhance the robustness of the method, a variation of these strategies has been investigated in which the is computed on training patterns to which a noise has been added. Also, to improve the accuracy of the reconstruction, past signal measurements are used as further input to the reconstruction models and the reconstruction of the faulty signals is iterated until satisfactory convergence. The remaining part of the paper is organized as follows. Section II states the problems of fault detection and signal reconstruction; Section III recalls briefly the RFSE technique; Section IV describes the methodologies for fusing multiple models outcomes, while Section V presents the results from a set of experiments based on the reconstruction of the signals measured by the sensors of a simulated PWR pressurizer; finally, Section VI presents the conclusions and describes potential future work.

II. FAULT DETECTION AND SIGNAL RECONSTRUCTION

A sensor measurement provides an estimate of the true value of a physical quantity; ideally, it should be , but in practice the measurement is affected by random noise and sometimes it also deviates from the true value because of a sensor fault.

Under normal conditions, sensor measurements are strongly correlated by the physical relations governing the process. Then, validation of sensor measurements is possible by taking advantage of the redundant information coming from the correlated measurements; this allows detecting incorrect measurements, identifying faulty sensors and reconstructing the correct values of measurements. More specifically, fault detection and identification aim at detecting if a measured signal is no longer a reliable estimation of the true signal value , due to an identified sensor malfunctioning. The detection of a sensor fault is often achieved by comparison between the actual sensor measurement and the measurement estimate provided by a physical or empirical model reproducing the response of the sensor in non-faulty conditions; many techniques, like the sequential probability ratio test (SPRT) technique [START_REF] Wald | Sequential Analysis[END_REF], [START_REF] Yang | Monitoring and Uncertainty Analysis of Feedwater Flow Rate Using Data-Based Modeling Method[END_REF], are based on the monitoring of the residual signal . Faulty signal reconstruction uses the information contained in other correlated signals to infer a more satisfactory estimate of the true value of the physical quantity: we call this estimate the reconstructed value of the faulty signal . The main objective of the present paper is the investigation of techniques for reducing the deviation between the estimate of the quantity measured by a faulty sensor and its real value . The improvement of the estimate brings two main benefits:

i. a more reliable estimate of the residual which leads to a more prompt and accurate detection of the sensor faults; ii. a more accurate reconstruction of the faulty signal in case of sensor fault. The first beneficial effect above mentioned is not further developed here, since the definition of an explicit procedure for sensor faults detection and identification is outside the scope of the present work. , is available for building a model for the reconstruction of a detected faulty signal. Given the typically large value of , a single model cannot perform the reconstruction task with the desired accuracy and reliability. Then the signal set must be partitioned into subsets, for each of which a reconstruction model is built. Actually in the ensemble approach here proposed the subsets of signals are overlapping (i.e., two subsets may contain signals which are the same), a reconstruction model is built for each subset and the reconstruction of a signal contained in different models is obtained by aggregating their outcomes, within an ensemble approach.

III. THE RANDOM FEATURE SELECTION ENSEMBLE APPROACH

In this work, the individual models employed in the ensemble for signal reconstruction are based on the PCA technique; the motivation is the need for a trade-off between computational time and accuracy of the reconstruction. Other techniques such as AANN and AAKR have given satisfactory accuracies [START_REF] Chevalier | Assessment of statistical and classification models for monitoring EDF's assets[END_REF] but are unfit to be used in an ensemble approach, due to their very high computational cost (long training time for AANN, long on-line execution time for AAKR). PCA is based on the transformation of the data to a smaller set of variables which are linear combinations of the original variables; the transformation is made so as to retain as much information as possible. The principal components correspond to the directions in which the projected observations have the largest variance: the largest variance comes to lie on the first principal component, the second largest variance on the second principal component, and so on. Principal components are included in the transformation until the cumulative variance associated to the retained components reaches a fixed fraction of the total variance. Appendix B gives some details on the application of PCA to signal reconstruction.

In this work, signal grouping is performed by the RFSE technique which consists in randomly sampling from the available signals, with replacement, subsets , , each constituted by signals [START_REF] Bryll | Attribute bagging: Improving accuracy of classifiers ensembles by using random feature subsets[END_REF]. This guarantees high signal diversity in the overlapping groups upon which the PCA models are built and allows for rapid construction of the signal groups. Furthermore, randomly selecting the signals in the groups with a reasonable choice of the group size parameters and , can basically guarantee coverage of all the signals in the ensemble with adequate redundancy [START_REF] Fantoni | On-line monitoring of instrument channel performance in nuclear power plant using PEANO[END_REF], [START_REF] Atiya | Novel methods for the feature subset ensembles approach[END_REF].

As mentioned above, the diverse signal groups generated are used as bases for developing a corresponding number of PCA reconstruction models. To do this, the data set of signal patterns available is partitioned into a training set (made of patterns) and a test set (made of patterns). The former is used to train the individual models, whereas the latter is used to verify the ensemble performance in the signal reconstruction task.

Each signal is present in a number of groups and thus a corresponding number of individual PCA models provide its reconstruction. Different methods can be used to aggregate the outcomes of these individual models in the ensemble, to get the final reconstructed value. Some of these will be discussed in the next Section IV.

After the ensemble of PCA models has been trained, its performance is verified on the test data: the entire set of n signals is fed into the ensemble of models which returns their reconstructed values. In particular, each model receives in input and reconstructs only the signals belonging to the corresponding subset .

IV. LOCAL FUSION STRATEGIES FOR MODELS OUTCOME AGGREGATION

A. Motivations

Let be the reconstruction of the -th signal of a test pattern made by a generic model built on the training set , and its real value. Since the model depends on the actual data used to build it, the randomness of the training set implies that the response of the model at any point is a random variable. In order to evaluate in a probabilistic sense the quality of as estimator of the mean squared error over all possible training datasets is commonly used:

(1)

The can be decomposed in a variance term plus a bias component [START_REF] Geman | Neural networks and the bias/variance dilemma[END_REF]:

(2) where is the expected value of over all possible training sets . Notice that the bias component represents the mean squared error that is obtained in the ideal situation where the outcomes of infinite models built over all possible training sets are combined by means of the arithmetic average, since this combination technique compensates the variance component, i.e., the fluctuations around the mean value of the reconstructions obtained by models built on different datasets.

In [START_REF] Bonissone | Fast meta-models for local fusion of multiple predictive models[END_REF], two issues are raised in this regard: 1) although an ensemble of models built on different training sets actually boosts the accuracy of the reconstruction by reducing the variance component, the simple average fusion strategy of the single model outcomes cannot effectively reduce the bias; 2) the mean squared error also depends on the test pattern under reconstruction since all terms in (2) depend on it. Therefore, different patterns in the signals space will usually have different error profiles of the model reconstructions. In this context, a procedure is defined in [START_REF] Bonissone | Fast meta-models for local fusion of multiple predictive models[END_REF] for the aggregation of an ensemble of models trained over different datasets generated through the bootstrap technique. The proposed fusion strategy aims at favoring the best performing models in each region of the input space and, at the same time, effectively correcting not only the variance component of the error, but also the bias.

In this work, this local strategy for the aggregation of multiple outcomes is applied in the context of multiple models in which diversity is injected through the RFSE technique, since the considerations about the dependence of the model error on the input pattern remains true for each model built on the signal subset .

To further delve into the aspect of local variation of the error profile, Fig. 2 and Fig. 3 show some examples of models reconstructions for the water level signal in the PWR pressurizer application. Fig. 2 represents the time evolution of the water level in the pressurizer during a plant transient (top) and the residual error made by model 12 in its reconstruction (bottom). Notice that the reconstruction error varies in the different time regions and seems correlated to the signal value. For example, the dotted boxes highlight regions corresponding to flat peaks of the level measurement signal, to which correspond similar trends of the residual error. Thus, in order to correct the error in the reconstruction of a test pattern falling in the region of a flat peak, it would seem more effective to compute the mean error of model 12 considering only the similar (neighbor) training patterns of the highlighted regions, instead of using the mean error on all the available training patterns.

Fig. 3 shows the performances of models 8 and 11 in the reconstruction of the same signal of Fig. 2; although

B. Models Outcome Aggregation

In general, the aggregation of the models outcomes requires to associate a weight and a bias correction to the reconstruction of each model . The idea is to correct the values of by subtracting the estimated bias and to combine with the other models estimates by means of a weighted average:

(3)
According to the considerations of Section IV-A, the ensemble performance can be increased if both the bias corrections and the weights are computed locally. Fig. 4 shows the overall framework for the reconstruction of the generic signal .

The following two subsections provide some details on how to select the neighboring training patterns to be used for the computation of the bias and weight parameters and how to associate them to the models of the ensemble. is formed by the training patterns nearest to the test pattern (Fig. 6). Notice that in this way the cardinality of is fixed. A drawback of this approach is that when the neighborhood of the test pattern has a low density of patterns, its dimensions increase and the precision of the estimated parameters is reduced.

2) Local Weights and Bias Estimation: In order to evaluate the performance of the generic model in the reconstruction of the -th signal of a given test pattern, the model reconstructions of the training patterns in the neighborhood set are compared to the measured signal , which, for a fault-free training set, is supposed to be equal to the true signal value , except for the presence of noise. The true signal value cannot be used in real applications, since it is not available. To this purpose, the following two indicators are considered:

• the local mean error (4) In this work three different relationships between the weights and the local mean error are investigated:

1) The weight is the inverse of the mean absolute error:

(7)

2) The weight is proportional to the logarithm of the inverse of the normalized mean absolute error:

(8)
The normalization factor is the maximum value of the residual error over all patterns of the training set and all models ; being always greater than the , it forces the argument of the logarithm to be greater than 1, thus maintaining the value of the weight positive. This strategy, which is inspired by the assignment of the weights in the adaboost algorithm [START_REF] Polikar | Bootstrap-inspired techniques in computational intelligence[END_REF], is here employed to attenuate the difference in the weights of groups with very different local mean absolute errors; this can be useful, in particular, when an outstanding performance of a group on the training set is not due to the real quality of the model, but derives from overfitting.

3) The Borda Count method [START_REF] Polikar | Ensemble based system in decision making[END_REF]. The estimated local error is used to make a ranking of the different models and to assign them a score , , according to their position in the ranking, being 1 the score associated to the worst performing model and that to the best performing one: [START_REF] Baraldi | Genetic algorithms for signal grouping in sensor validation: A comparison of the filter and wrapper approaches[END_REF] Another possibility to overcome the overfitting problem is to compute the model performance not on the training data but on a different dataset. In this respect, an artificial noise has been added to the training data before their reconstruction so that the training error used to estimate the weights will also account for the robustness of the model, reducing in this way the risk of giving high weights to models that are overfitting.

All the three proposed strategies defined to assign a weight to the different models (( 7), ( 8), ( 9)) can be used with the local mean absolute error computed without bias and with the disturbed training error instead of the .

V. APPLICATION

The case study used to test the LF strategies presented in the previous Sections concerns the reconstruction of signals measured in a controlled PWR pressurizer. The pressurizer maintains the pressure in the primary system of the plant at an opportune value; in particular, for safety and operational reasons, it is fundamental to have a pressure sufficiently high to avoid boiling of the coolant while the plant is operating.

The data for the analysis (training and test patterns) were generated by means of a SIMULINK model purposedly developed. The SIMULINK model is based on a system of non-linear differential equations derived from the mass and energy conservation laws applied to the two regions of vapour and liquid of the pressurizer; exchanges between the two regions, due to evaporation of liquid and condensation of steam, are taken into account. The water and steam properties used in the SIMULINK model are based on the Industrial Formulation IAPWS-IF97 [26].

The variables controlled in the pressurizer are the level and the pressure. Their control is achieved by acting on the heaters, the charging/letdown flows and the relief valve (Fig. 7). The relief valve is a safety device which opens only when the pressure exceeds a limit value (here set equal to 165 bar). Since this limit

TABLE I INITIAL CONDITIONS OF THE PRESSURIZER

pressure value has never been reached during the simulated transients, the relief valve flow signal has not been considered in this work.

Fig. 8 reports the control scheme adopted. The measured level is compared with the level set point: in case of low level, the charging/letdown flows correct the mismatches and the heaters are turned off independently from the value of the pressure. Besides, when the measured pressure is lower than its set point, the heaters are turned on. The control signal sent to charging/letdown flows, relief valve and heaters is obtained by using a proportional-integral-derivative (PID) controller.

Starting from the initial conditions reported in Table I, 50 transients of 100 time instants each have been simulated by randomly changing the value of the surge line flow entering the pressurizer. Fig. 9 reports the evolution of the signals in one of the simulated transients.

The surge line flow depends on the reactor power; in particular, when the Primary Heat Transport (PHT) system pressure rises, the higher pressure pushes some coolant from the PHT into the pressurizer (in-surge flow); on the contrary, in case of low PHT system pressure, some liquid moves from the pressurizer into the PHT system circuit (out-surge flow). In this application, the transients have been generated by assuming surge flows in the range of In order to test the method on pseudo-realistic data, white noise has been added to each signal according to engineering considerations on the sensors accuracy [START_REF] Hashemian | Sensor Performance and Reliability[END_REF]- [START_REF] Johnson | Nuclear Reactor Controls and Instrumentation[END_REF]. Table II reports the standard deviations of the considered noises.

The reconstruction task has been carried out both for clean and noisy data.

A. Sensor Measurement Reconstruction

In order to catch the dynamic evolution of the system, the set of input measurements which are given as input to the reconstruction model is constituted not only by the current values of the 7 signal measurements but also by a sliding window that contains measurements at previous times: in particular, for each one of the 7 signals, the current and the previous 9 measurements have been considered, leading to a total number of 70 input measurements (Fig. 10).

In this application, 2730 patterns, corresponding to 30 transients, have been used to train the PCA models and to compute the single model local mean errors and local mean absolute errors which are the quantities used for the settings of the ensemble local weights and bias, whereas the remaining 1820 patterns, taken from 20 different transients, are used to test its performance. Notice that, in order to avoid any overfitting of the data, the ensemble performance is evaluated on test data that are different from the training data used to evaluate the single model performances which are necessary for setting the ensemble parameters.

TABLE II STANDARD DEVIATION OF SENSORS NOISE

The optimal number of models to be used in the ensemble and the optimal number of input signals in each PCA model have been established through optimization of the performances on fault-free noisy validation data (different from both the training and test data). The best results have been obtained by using 20 PCA models, each one with 25 input measurements. In this way, all signals are reconstructed by at least 7 models. Notice that in order to reduce the computational time, redundancies higher than 7 have not been considered although they may increase the ensemble performance.

The fraction of retained variance for each PCA model is .

Finally notice that in the evaluation of the performance of the ensemble model, only the reconstruction of the 7 signal measurements at the present time and not their past values have been considered for the estimation of the mean absolute error. The values refer to the signals normalized in the range [0.2; 1].

The following Section V-B considers the problem of reconstructing the true signal values in the case in which there are no failures in the sensors. Since the comparison between the reconstructed signal values and the sensor measurements is the base for the identification of anomalies in the sensors, the reconstruction is fundamental for a correct sensor fault detection. In Section V-C, the problem of signal reconstruction in the case of sensor failure will be tackled.

In both cases, the use of a set of signals containing also previous measurements (9 in this work) has demonstrated to be effective. In case of noisy data the of a single PCA model built on the 7 present signals is

, while considering the past measurements the obtained is , with a reduction of the error of 5.5%. Moreover, as shown in Fig. 11, referring to the case of a constant bias fault of the pressure sensor, a significant improvement is obtained in the signals reconstruction during the first instants after the initiation of a fault. When past information is considered, the signals used for the reconstruction, in the 9 instants just after the fault, have in part been measured before the beginning of the faulty condition; this allows an improvement in the reconstruction of the faulty signal, which facilitates the fault detection and identification. Finally, if the pressure signal is reconstructed by the iteration procedure explained in Section V-C, a of is obtained using the past measurements in input, which is 60% lower than the of obtained when only the current measurements are used.

B. Signal Validation for Sensor Fault Detection

In this Section, the results obtained by using the proposed LF strategies are compared with those obtained by using two global approaches: 1) the Simple Mean (SM) and 2) the Globally Weighted Average (GWA). The latter is realized by assigning to the -th model reconstruction of signal a weight proportional to the inverse of its , computed on the whole training dataset [START_REF] Bonissone | Fast meta-models for local fusion of multiple predictive models[END_REF].

Table III shows the performances obtained in the reconstruction of clean (top) and noisy (bottom) data with the different fusion strategies.

Several runs have been carried out to set the values of the hyper-rectangle length and the number of nearest neighbors , in order to get the best results for fair comparison. The value of has to be high enough to avoid empty neighborhoods in correspondence of any test pattern. Notice that in the case of clean data, the value which leads to the best results is the smallest one respecting the condition. This is in accordance with the fact that the optimal number of nearest neighbors is very low in case of clean data. When noise is added to the test data measurements, the optimal values of and increase. This is due to the fact that, when noise is present, the error committed by a model in the reconstruction of signal of a generic pattern is the result of two contributions: the modeling error and the The reported percentages refer to the improvement of the performance with respect to the SM. Hyper-rectangle-based neighborhood; k-nn-based neighborhood.

noise error , i.e. the error resulting from having inaccurate inputs to the model. The error component cannot be properly estimated using the training data because of the random nature of the noise. The results from the average of the local errors made on the neighbors ((10)).

(10) When we consider a large number of neighbors, the last term of (10) can be approximated to 0, so that the contribution of noise to the error estimate is eliminated. However, as a consequence of the increased value of , the local character of the procedure is reduced because a wider region of the input space is considered; thus the quality of the parameters estimation is decreased and the performance of local and global techniques tend to become more similar.

In order to distinguish the contribution to the improvement of the performance of the bias correction from that of the local weights, Table III reports in the third row the performance obtained by assigning a constant weight equal to to all the models and correcting the bias locally, according to (3). In the fourth row, Table III reports the results obtained by applying the complete LF scheme of (3), which includes bias correction and local weighting, using a weight proportional to the inverse of the (( 7)). Notice that the main contribution to the improved performances of the local fusion is due to the local weights, whereas the bias correction conveys only a small amelioration, particularly in case of noisy signals.

Finally, notice that the best strategy for the selection of the neighbors seems to be the -nn-based neighborhood, given that it leads to lower reconstruction errors than the hyper-rectanglebased neighborhood, in both cases of clean and noisy data.

C. Reconstruction in Case of Faulty Sensor

When one of the plant sensors is faulty, a faulty signal is sent in input to the PCA models which include that signal. The ensemble should however still provide a good estimate of the true value of the signal, by exploiting the correlated information coming from the non-faulty signals. To verify the performance of the different local fusion techniques in dealing with sensor faults, a number of tests have been performed as described in this Section.

Two faults types are considered: constant bias faults (CF) and noise faults (NF), affecting both the pressure and level sensors. These two sensors have been chosen since their measurements influence the controller decisions and thus are particularly critical for plant operation. The faults are supposed to start in all the transients at , with various intensities. Fig. 12 shows some examples of sensor faults analyzed.

The LF signal reconstruction approach based on the -nnbased neighborhood (which has given the best results in the previous Section), has been considered with in case of clean training data and in case of noise. The different techniques introduced in the previous Section IV-B-2 for the estimation of the local weights have been investigated and compared to the SM global approach. Table IV reports the obtained in the reconstruction of the true value of the signal measured by the faulty sensor.

The results show that in the case of faults affecting the pressure sensor, the SM performs better than the LF techniques. Also, among the local techniques, the best performing weighting mechanism is the one based on logarithmic weights. The reason of this is that some models achieve a very high performance in the reconstruction of the training data just because they tend to exactly reproduce the input measurements and thus they are associated with very high weights; however, such models are not robust in the case of sensor faults; nevertheless, according to the good performance on the training set, the LF Fig. 13. Weights associated to the 7 models of the ensemble using the training data and corresponding reconstruction of the pressure during a plant transient. mechanism will assign to them weights excessively high compared to the other models. A logarithmic weight mechanism reduces the weights associated to these non robust models and thus improve the local fusion performance with respect to the technique, but its performance is still worse than the SM weighting.

In order to investigate this effect, Fig. 13 (left) shows the reconstruction of the pressure in a transient obtained by the seven models of the ensemble which use as input measurements the actual pressure, while Fig. 13 (right) shows the weights assigned to the models according to their performance on the training dataset. Notice that from approximately , model 1 has an exceptionally high weight even if, when the fault is present, it is the worst performing one. In this example, the use of logarithmic weights permits to reduce the predominance of model 1 in the reconstruction, thus improving the ensemble performance. However, even with a logarithmic weight mechanism, the weight associated to model 1 would be greater than the weight associated to the actual best performing model.

To overtake this limitation of the local fusion techniques, it is possible to estimate the actual robustness of the models by computing the reconstruction error not on the training dataset but on a new data set obtained by adding noise on the training data. The results shown in the 4th and 8th columns of Table IV confirm our expectation: the new weights assigned to the pressure reconstructions of Fig. 13 (left) are shown in Fig. 14.

Finally, notice that once the fault has been detected, in order to correctly reconstruct the true signal value, it is possible to repeat the reconstruction of the faulty signal times, using as input to the ensemble model at the -th iteration the value of the reconstructed measurement at the previous iteration [START_REF] Chevalier | Assessment of statistical and classification models for monitoring EDF's assets[END_REF]. This technique has proven very effective in the reconstruction of the measurements of faulty sensors as it can be seen in Table V, which reports the results obtained by repeating the reconstruction . First of all, notice the reduction of the obtained by repeating the reconstruction. Furthermore, it is interesting to observe that computing the weight by adding a noise on the training data (NT technique) leads to the best signal reconstruction in the case of noisy data, whereas, as expected, it is not effective in case of clean data where the best results are obtained by the logarithmic weighting technique.

D. Results Discussion

The results reported in Sections V-B and V-C show that the local fusion approach performs better in the signal validation task, than the two global techniques considered, SM and GWA.

Concerning the choice of the local neighborhood for the fusion, the best results are obtained using a nearest neighbors approach. When data are affected by noise, the number of neighbors to be considered increases significantly.

In case of sensor failure, the best results are obtained by assigning the model weights proportional to the inverse of the absolute error computed on training data to which a noise has been added. This strategy permits to avoid assigning high weights to overfitting models, i.e., models that perform very well on the training data but are not suitable for the reconstruction of faulty signals.

A significant improvement of the results can also be obtained by repeating the reconstruction task iteratively several times, each time feeding the reconstructed value of the previous iteration.

VI. CONCLUSIONS

The ensemble approach has been shown to be an effective way for increasing the accuracy and robustness of signal reconstruction.

The goodness of the results depends on the methods used for aggregating the multiple outcomes of the ensemble models. This work has shown that local approaches can increase the performance of the ensemble in the signal reconstruction task. The best results have, in fact, been obtained when the parameters of the aggregation procedure dynamically change depending on the values of the signals to be reconstructed, i.e., the location of the test pattern in the signal space relative to the set of training patterns.

Realistic noise and sensor faults have been considered in order to show that the proposed local fusion approach is effective for boosting model prediction accuracy, in practical settings.

Furthermore, the use of a set of signals containing also previous measurements has been demonstrated to be an effective way to reduce the detrimental contribution of noise or faults.

In the case of sensor faults, the proposed procedure of faulty signal reconstruction, based on iterations of the reconstruction step and memorization of the reconstructed signal for use at the next iteration, improves the accuracy of the reconstruction of the value of the signal.

Although the methodology has been demonstrated in the context of multiple models based on PCA and constructed on groups of signals selected by the random feature selection technique, it seems to be broadly applicable to other empirical modeling paradigms. Investigations in this direction could be of interest, considering in particular other fast empirical regression techniques.

Finally, the local performance identification and parameters computation have to be performed online; this could be a costly procedure that might not be suitable when response time is an issue. Future work will be devoted to the set up of a procedure for the offline assignment of the local fusion parameters values.

  anomalies in their functioning and for reconstructing the correct values of the signals measured. In this work, the task of sensor monitoring and signal reconstruction is tackled with an ensemble of Principal Component Analysis (PCA) models. The novelty of the work consists in the investigation of local fusion (LF) strategies for the aggregation of the outcomes of the different models of the ensemble. In the reconstruction of a signal, each model of the ensemble is assigned a weight and a bias related to the error committed in the reconstruction of training patterns similar to the one under reconstruction. Iteration of the reconstruction procedure and use of past measurements of the signals are introduced for improved performance. The proposed methodology is applied to a case study concerning the reconstruction of seven signals in the pressurizer of a Pressurized Water Reactor (PWR) nuclear power plant. Index Terms-Local fusion, pressurizer, random feature selection ensemble, signal monitoring, signal reconstruction.

Fig. 1 .

 1 Fig. 1. Approach to the reconstruction of signal f contained in groups F ; F ; . . . ; F .

Fig. 1

 1 Fig. 1 reports a sketch of the flow of modeling for the reconstruction of a faulty signal. A set of sensors' signals ,, is available for building a model for the reconstruction of a detected faulty signal. Given the typically large value of , a single model cannot perform the reconstruction task with the desired accuracy and reliability. Then the signal set must be partitioned into subsets, for each of which a reconstruction model is built. Actually in the ensemble approach here proposed the subsets of signals are overlapping (i.e., two subsets may contain signals which are the same), a reconstruction model is built for each subset and the reconstruction of a signal contained in different models is obtained by aggregating their outcomes, within an ensemble approach.In this work, the individual models employed in the ensemble for signal reconstruction are based on the PCA technique; the motivation is the need for a trade-off between computational time and accuracy of the reconstruction. Other techniques such as AANN and AAKR have given satisfactory accuracies[START_REF] Chevalier | Assessment of statistical and classification models for monitoring EDF's assets[END_REF] but are unfit to be used in an ensemble approach, due to their very high computational cost (long training time for AANN, long on-line execution time for AAKR). PCA is based on the transformation of the data to a smaller set of variables which are linear combinations of the original variables; the transformation is made so as to retain as much information as possible. The principal components correspond to the directions in which the

Fig. 2 .

 2 Fig. 2. Top, true signal value; bottom, reconstruction error of model 12.

Fig. 3 .

 3 Fig. 3. Comparison between model 8 and model 11 reconstruction errors. Dotted boxes highlight some regions where the model built on the 11th group of signals performs better than the one built on the 8th group.

Fig. 4 .

 4 Fig. 4. Framework of the locally weighted fusion of model outcomes.

1 )Fig. 5 .

 15 Fig. 5. Hyper-rectangle-based neighborhood. The crossed circle represents the test pattern; the other circles are training patterns.

Fig. 6 . 5 )

 65 Fig. 6. k-nn-based neighborhood. The crossed circle represents the test pattern; the other circles are training patterns.

Fig. 7 .

 7 Fig. 7. Input, output and control variables of the SIMULINK model of the pressurizer.

Fig. 8 .

 8 Fig. 8. Pressurizer control scheme.

  [ 40; 40] kg/s, positive values indicating in-surge flow and negative values out-surge flow. During the simulations, the following seven signals are recorded: • the level in the pressurizer • the liquid temperature • the steam temperature • the pressure • the charging/letdown flows • the surge line flow • the heaters power.

Fig. 9 .

 9 Fig. 9. Evolution of the seven signals during a plant transient.

Fig. 10 .

 10 Fig. 10. Input pattern used for the reconstruction of the current signal measurement.

Fig. 11 .

 11 Fig. 11. Comparison of the reconstruction of the pressure signal undergoing a constant bias fault, obtained from a single PCA model built on the current values of the 7 signals and the one built on the set containing also the 9 most recent measurements.

Fig. 12 .

 12 Fig. 12. Examples of signals in fault conditions.

mae( 10 )

 10 OF THE FAULTY SIGNAL RECONSTRUCTION SM = simple mean; LF = local fusion, k-nn-based neighborhood; 1=" = weights computed with eq:(7); Log = weights computed with eq:(8); Borda = weights computed with eq:(9); NT = weights computed with eq:(7), adding noise on training data; CF = constant fault, NF = noise fault.

Fig. 14 .

 14 Fig. 14. Weights associated to the reconstructions of Fig. 13, obtained adding a noise to the training data.

mae( 10 )

 10 OF THE FAULTY SIGNAL RECONSTRUCTION AFTER 5 ITERATIONS SM = simple mean; LF = local fusion; 1=" = weights computed with eq:(7); Log = weights computed with eq:(8); NT = weights computed with eq:(7), adding noise on training data; CF = constant fault, NF = noise fault.

TABLE III MEAN

 III ABSOLUTE ERROR (mae) OBTAINED BY USING DIFFERENT COMBINATION STRATEGIES

TABLE IV

 IV

TABLE V

 V

Manuscript received November 05, 2009; revised January 21, 2010 and January 29, 2010. Current version published April 14, 2010.The authors are with the Department of Energy, Polytechnic of Milan, 20133 Milano, Italy (e-mail: piero.baraldi@polimi.it; antonio.cammi@polimi.it; francesca.mangili@mail.polimi.it; enrico.zio@polimi.it).Digital Object Identifier 10.1109/TNS.2010.2042968

ACKNOWLEDGMENT

The authors wish to thank the reviewers for their constructive comments which have allowed improving the paper.

In this Appendix, the procedural steps of Principal Component Analysis (PCA) are briefly sketched [START_REF] Baraldi | Genetic algorithms for signal grouping in sensor validation: A comparison of the filter and wrapper approaches[END_REF]. Let be the data set matrix with row vector , , which gathers the signals patterns measured at different time instants . The purpose of the PCA is to identify a -dimensional subspace in which the most of the data set variation is retained and the least information is lost.

Let be a matrix constituted by orthonormal principal vectors which are the eigenvectors associated to the eigenvalues of the covariance matrix of :

where it is assumed, without loss of generality, that is a diagonal matrix with diagonal elements in decreasing order of magnitude.

represents an orthonormal basis for the data set and the transformation (B2) projects the data in the new basis. The data set has now two representations: when intended in the original basis, the -th pattern is the vector with components ; when intended in the orthonormal basis, the same -th pattern is represented by the -th rows of the principal component matrix with component along . To get the inverse relation, we right-multiply by to obtain the data set in the original basis:

(B3)

Eqs. (B2) and (B3) represent an optimal linear transformation of the data matrix in the sense that the new set of variables, the principal components (PCs), are uncorrelated and ordered so that the first few retain the most of the variation present in all of the original variables [START_REF] Jolliffe | Principal Component Analysis[END_REF]. Up to this point the equations are exact and the data values are transformed in both senses without any loss of information.

The PCA approximation consists in mapping the observations in a subspace identified by vectors. For that, the matrix is partitioned as follows:

where is the submatrix constituted by the first columns of and is the submatrix constituted by the last columns of . The column vectors in and constitute the bases of the two mutually orthogonal subspaces and in which has been divided. In terms of the above submatrices, matrix can be rewritten as:

where is the principal part of the data explained by the first eigenvectors:

(B6)

The -th row of , namely , is the orthonormal projection of onto . If all the information about the data set essentially lies in a -dimensional space (apart from small components in given by as stated in (B5), then each observation vector can be approximated by its orthonormal projection plus a residual vector in which is postulated to be independent of :

The best residual vector is the one which, on average, minimizes the absolute value of the square error between the real and approximated data patterns, i.e.:

(B8)

The calculations for the minimization can be found in [START_REF] Baraldi | Genetic algorithms for signal grouping in sensor validation: A comparison of the filter and wrapper approaches[END_REF]; the resulting optimal residual vector takes the form:

(B9) where the vector is the arithmetic average of the observation vectors, i.e. the average value of the signals. In particular, the -th component of is the arithmetic average of the -th column of .

Then, from (B7), the expression for the PCA approximation of the data pattern is given by: (B10) or in matrix form, (B11)

As demonstrated and explained in details in [START_REF] Jolliffe | Principal Component Analysis[END_REF], [START_REF] Diamantaras | Principal Component Neural Networks: Theory and Applications[END_REF] and [START_REF] Marseguerra | Lecture Notes on Principal Compponents Analysis[END_REF], the minimum error corresponding to the optimal residual vector can be written as:

and the fact of choosing as the eigenvectors of the covariance matrix corresponds to minimizing the minimal error of (B12).

Finally, for the specific problem of signal validation and reconstruction the orthonormal basis has been obtained by (B1) from the covariance matrix of the pairwise signal correlations between the signals. The number of eigenvectors to be retained is determined as follows: components corresponding to the largest variance of the data are retained until the cumulative variance associated to the retained components reaches a fraction of the total variance. Furthermore, in order to simplify the calculations, the time trends of the signals have been previously normalized so that their mean is zero and their standard deviation equals 1. This allows skipping the computation of the residuals, since and, according to (B10) (B13)


Marco Lenci

Serge Troubetzkoy

Infinite-horizon Lorentz tubes and gases: recurrence and ergodic properties

Keywords: MSC 2010 37D50, 37A40, 60K37, 37B20, 36A25

We construct classes of two-dimensional aperiodic Lorentz systems that have infinite horizon and are 'chaotic', in the sense that they are (Poincaré) recurrent, uniformly hyperbolic and ergodic, and the first-return map to any scatterer is K-mixing.

In the case of the Lorentz tubes (i.e., Lorentz gases in a strip), we define general measured families of systems (ensembles) for which the above properties occur with probability 1. In the case of the Lorentz gases in the plane, we define families, endowed with a natural metric, within which the set of all chaotic dynamical systems is uncountable and dense.

Introduction

A Lorentz system is a dynamical system of a point particle moving inertially in an unbounded domain (in this paper we consider only planar domains) endowed with an infinite number of dispersing (i.e., locally convex) scatterers. When the particle hits a scatterer, which is regarded as infinitely massive, it undergoes an elastic collision: the angle of reflection equals the angle of incidence.

The most popular such system is undoubtedly the Lorentz gas, devised by Lorentz in 1905 [Lo] to study the dynamics of an electron in a crystal; the term 'gas' was introduced later in the century, when versions of the Lorentz model were used to give a statistical description of the motion of a molecule in a gas. 1 Date: July 12, 2011.

1

To our knowledge, the first appearance of this model within the scope of the kinetic theory of gases is in a 1932 textbook of theoretical physics [J]. The model is used briefly for the estimation of the number of collisions per unit time of a molecule in a gas; however, no mention of Lorentz is made. The first occurrence of the phrase 'Lorentz gas' in the scientific literature seems to date back to 1941 [Gr]. Curiously, the first occurrence of the phrase 'Lorentzian gas' that we are aware of is found in the article preceding [Gr] in the same issue of the same journal [HIV]. Finally, it is interesting to notice that, although Lorentz's original papers [Lo] are

In the past century, Lorentz systems have been preferred models in the fields of statistical physics, optics, acoustics, and generally anywhere the diffusive properties of a chaotic motion were to be investigated. Throughout this history, for reasons of mathematical convenience, the models that were studied most often and most deeply were periodic (i.e., the configuration of scatterers was invariant for the action of a discrete group of translations) and with finite horizon (i.e., the free flight was bounded above). Only recently have aperiodic Lorentz systems come to the fore [START_REF] Lenci | Aperiodic Lorentz gas: recurrence and ergodicity, Ergodic Theory Dynam[END_REF][START_REF] Lenci | Typicality of recurrence for Lorentz gases, Ergodic Theory Dynam[END_REF][START_REF] Dolgopyat | Limit theorems for locally perturbed planar Lorentz processes[END_REF][START_REF] Cristadoro | Recurrence for quenched random Lorentz tubes[END_REF][START_REF] Seri | Recurrence and higher ergodic properties for quenched random Lorentz tubes in dimension bigger than two[END_REF][START_REF] Troubetzkoy | Typical recurrence for the Ehrenfest wind-tree model[END_REF] (within the scope of dynamical systems, that is, aperiodic models had already been studied in other contexts; e.g., [Ga, BBS, BBP]). However, while there is a considerable literature on infinite-horizon periodic Lorentz gases, almost nothing is known, at least to these authors, on aperiodic Lorentz systems with infinite horizon.

In this note we consider 2D Lorentz gases and also Lorentz tubes, that is, Lorentz systems confined to a strip of R 2 [CLS, SLDC]. We construct billiards that have infinite horizon and possess the ergodic properties that one would expect of these chaotic systems, such as ergodicity and strong mixing properties. For these infinite-measure preserving dynamical systems, it turns out that (Poincaré) recurrence is not only a necessary but also a sufficient condition for ergodicity; this is a consequence of the hyperbolic structure that our systems can be shown to have [START_REF] Lenci | Aperiodic Lorentz gas: recurrence and ergodicity, Ergodic Theory Dynam[END_REF][START_REF] Cristadoro | Recurrence for quenched random Lorentz tubes[END_REF]. As for mixing, since a universally accepted definition of mixing is not available in infinite ergodic theory (see, e.g., the discussion in [START_REF] Lenci | On infinite-volume mixing[END_REF]), we characterize this aspect of the dynamics by proving that certain first-return maps are K-automorphisms (which implies strong mixing). This is again a consequence of recurrence and hyperbolicity.

Another important question that we aim to discuss is that of the typicality of such ergodic properties. One would expect that, when the effective dimension is one (Lorentz tubes) or two (Lorentz gases), recurrence, and all the stochastic properties that it entails, hold for "most" systems (and one would expect the former case to be more easily worked out than the latter).

We address this question by introducing the following two classes of billiards. A strip and the plane are partitioned into infinitely many congruent cells. In each cell we place a configuration of dispersing scatterers, i.e., a finite union of piecewise smooth closed sets, whose smooth boundary components are seen as convex from the exterior.

about electrons in a metal and not molecules in a gas, he treats the problem by deriving and solving a Boltzmann-like transport equation; see also [K].

The set of all the global configurations of scatterers gives rise to a family of dynamical systems of the same type. A natural distance between two configurations can be defined, which makes the above set a metric space. Furthermore, if the configuration is chosen according to a probability law, the family becomes a measured family, or an ensemble, of dynamical systems. This structure is often referred to as a quenched random dynamical system.

In the case of the tubes we give fairly explicit sufficient conditions for the above-mentioned ergodic properties, thus proving that, for many reasonable random laws on the "disorder" (including all non-degenerate Bernoulli measures), such properties hold almost surely in the ensemble; we state these results in Section 2. For the harder case of the gases, we prove that the set of the ergodic systems is uncountable and dense within the whole space; this is described in Section 3. Outlines of the proofs are given in Section 4.

It is worthwhile to mention that the billiards we construct have infinite but locally finite horizon, that is, though the free flight has no upper bound, no straight line exists that intersects no scatterers.

Lorentz tubes

Let C denote the unit square, which will be henceforth referred to as the cell. Let G 1 denote an open segment along one of the sides of C, say the left one, and G 2 the corresponding segment on the opposite side (via the natural orthogonal projection). G 1 and G 2 are called the gates of C.

A local configuration of scatterers is a "fat" closed subset Γ ⊂ C (this means that Γ is the closure of its interior) such that:

(A1) ∂Γ is made up of a finite number of C 3 -smooth curves γ i , which may only intersect at their endpoints (γ i is always considered a closed set). (A2) ∂C \ (G 1 ∪ G 2 ) ⊂ ∂Γ; and G 1 , G 2 do not intersect ∂Γ. (A3) Either γ i is part of ∂C or its curvature is bounded below by a positive constant (with the convention that positive curvature means that γ i bends towards the inside of Γ).

(A4) The angle formed by γ i and each intersecting γ j (or G j ) is nonzero.

Remark 1. Observe that (A1)-(A3) imply that the flat parts of ∂Γ are 'external' boundaries, in the sense that they can never be reached by a particle (the contrary would violate the hypothesis on the "fatness" of Γ). So we are dealing with dispersing, not semi-dispersing, billiards.

We consider a finite number Γ 1 , Γ 2 , . . . , Γ m of such local configurations (see Figs. 1 and2).

A Lorentz tube (LT) is a chain of cells C n (n ∈ Z) such that G 2 n , the right gate of C n , coincides with G 1 n+1 , the left gate of C n+1 . More precisely, call C n := [n, n + 1] × [0, 1]
the particular copy of C immersed in R 2 as indicated, and denote Ω := {1, 2, . . . , m}. Then, for ℓ := (ℓ n ) n∈Z ∈ Ω Z , we define the billiard table

(1) Q = Q ℓ := n∈Z C n \ Γ ℓn n , Figure 3. A Lorentz tube
where Γ ℓn n is the configuration Γ ℓn translated to the cell C n (see Fig. 3). The collection (Γ ℓn n ) n∈Z -equivalently ℓ-is called the global configuration of scatterers.

We henceforth say 'the LT ℓ' to mean both the table Q ℓ and the billiard dynamics defined on it. By this we mean, precisely, the dynamical system (M ℓ , T ℓ , µ ℓ )-more concisely, (M, T, µ)-where:

• M is the collection of all the line elements of the dynamics, i.e., all the position-velocity pairs (q, v), where q ∈ ∂Q and v is a unit vector based in q and pointing toward the interior of Q. (q, v) is meant to represent the dynamical variables of the particle right after a collision (v can be chosen unitary because in this Hamiltonian system the conservation of energy equals the conservation of speed). • T is the map that takes (q, v) into the next post-collisional line element (q ′ , v ′ ), along the billiard trajectory of (q, v); this map fails to be well defined only at a negligible set of line elements, which are called singular, cf. below. T is usually called the (standard) billiard map. • µ is the invariant measure induced on the Poincaré section M by the Liouville measure; it is well known that dµ(q, v) = n q , v dqdv, where n q is inner unit normal to ∂Q in q. It is easy to verify that the set of all singular points in phase space is null w.r.t. µ Finally, as is evident, µ(M) < ∞ if and only if the total length of ∂Q is finite. Notice that, by the definition of Q, a trajectory that intersects a gate G i n crosses it. (The trajectory that runs along the segment G i n will always be considered singular.)

We assume that there are two types of local configurations: the blocking configurations, corresponding to the index set Ω B := {1, 2, . . . , m ′ } (m ′ < m), and the non-blocking confugurations, corresponding to the set Ω N B := {m ′ + 1, m ′ + 2, . . . , m}. The former type verifies the following condition:

(A5) If a ∈ Ω B , any billiard trajectory that enters a cell with configuration Γ a must experience a collision before leaving it (Fig. 1).

An example of a configuration Γ a , with a ∈ Ω N B , is shown in Fig. 2.

Clearly, an LT ℓ ∈ Ω Z B has finite horizon, i.e., the free flight between two successive collisions has an upper bound. An arbitrary LT in Ω Z might not have this property.

A word a 1 a 2 • • • a l , with a i ∈ Ω, is called a factor of ℓ ∈ Ω Z if there exists n such that ℓ n+i = a i , for i = 1, 2, . . . , l. The factor is called blocking (respectively, non-blocking) if all the a i belong to Ω B (respectively, Ω N B ); it is called constant if they are all equal. The positive integer l is called the length of the factor.

The technique used by Troubetzkoy in [START_REF] Troubetzkoy | Recurrence in generic staircases[END_REF], when applied to the present models, easily implies that any LT in Ω Z which has arbitrarily long blocking constant factors (both forwards and backwards) is recurrent. Furthermore, Cristadoro, Lenci and Seri have shown that, for LTs in Ω Z B , ergodicity and recurrence are equivalent [CLS], and they imply K-mixing for suitable return maps (this last result is actually stated in [SLDC] but its proof applies as well to the models of [CLS]).

Our results on the Lorentz tubes are a combination and an extension of these ideas. To describe them we introduce some notation that will appear obscure at first, but will be explained shortly. For fixed ℓ, define g + 0 := g - 0 := 0 and, recursively for j > 0,

(2)

g ± j+1 := min k ∈ Z + ℓ ± j i=0 g ± i ± k ∈ Ω B .
In other words, g + j -1 (respectively, g - j -1) is the length of the j th nonblocking factor to the right (respectively, to the left) of the cell C 0 (with the convention that between two blocking cells there is a non-blocking factor of length 0). Notice that, if ℓ 0 ∈ Ω B , this coding reflects exactly the sequence of non-blocking factors of ℓ. Otherwise, there is a little difference which, see Remark 3 below, does not affect the upcoming statement.

Theorem 2. Assume (A1)-(A5). For any ℓ ∈ Ω Z which has arbitrarily long blocking constant factors, both forward and backwards, and such that both sequences (g ± j ) j∈N grow at most like a power-law, the corresponding dynamical system (M, T, µ) is:

(a) uniformly hyperbolic, in the sense that local stable and unstable manifolds exist at a.e. point of M, and the corresponding (Tinvariant) laminations are absolutely continuous w.r.t. µ (see, e.g., [START_REF] Lenci | Aperiodic Lorentz gas: recurrence and ergodicity, Ergodic Theory Dynam[END_REF] for details); the contraction coefficient of T n , along the stable direction, is bounded above by Cλ n and its expansion coefficient, along the unstable direction, is bounded below by C -1 λ -n , where C > 0 and λ ∈ (0, 1) are uniform constants;

(b) recurrent in the sense of Poincaré, i.e., given a measurable A ⊂ M, the orbit of a.e. point in A comes back to A infinitely many times; (c) ergodic, that is, if T (A) = A mod µ, then either A or its complement has measure zero.

Furthermore, the first-return map to any smooth component of ∂Q, of the type γ i as in (A1), is K-mixing.

Remark 3. In the above theorem, both hypotheses are shift-invariant in Ω Z , as a relocation of the origin on ℓ will produce at most a shift in (g ± j ) and a change of the first few terms. Therefore, when convenient-to do away with the problem mentioned in the previous paragraph-we make the convention that any ℓ is shifted to the left the minimum amount of times for C 0 to have a blocking configuration.

An LT as described in Theorem 2 need not have infinite horizon. To start with, it is necessary that the non-blocking configurations let free trajectories through, something that was not postulated. But more is needed as well. For instance, if one designs the local configurations Γ a , with a ∈ Ω N B , so that, for any l ∈ Z + , all non-blocking factors of length l admit a free flight of length ≥ l (and this is easy, cf. Fig. 2), then a necessary and sufficient condition for an LT to have infinite horizon is that at least one of the two sequences (g ± j ) is unbounded. Or one might ask for something less: for example, that just the constant nonblocking factors admit long free flights. This is enough to guarantee that very many LTs have an infinite horizon, cf. Corollary 4.

An important question is whether the LTs to which Theorem 2 applies are typical in Ω Z . This of course depends on the definition of 'typical'. One strong notion of typicality is the measure-theoretic notion, provided a probability measure Π is put on Ω Z (endowed with the natural σ-algebra generated by the cylinders). In this case, (Ω Z , Π) becomes a measured family (in jargon, an ensemble) of dynamical systems, which we call quenched random Lorentz tube.

Many reasonable measures Π ensure that the assertions of Theorem 2 hold Π-almost surely. Here is an example: Corollary 4. Let (p 1 , p 2 , . . . , p m ) be a stochastic vector, with p a > 0, and let Π be the Bernoulli measure on Ω Z relative to that vector (i.e., the unique measure for which Π(ℓ n = a) = p a , for all n). Then Π-a.e. LT in Ω Z has the properties stated in Theorem 2. Furthermore, if any non-blocking constant factor of length l (∀l ∈ Z + ) admits a free flight of length ≥ l, then a.e. LT has infinite horizon as well. a=1 p a and p N B := m a=m ′ +1 p a . By hypothesis, p B , p N B ∈ (0, 1). It is evident that, w.r.t. Π, the random variables (g ± j ) are i.i.d., with distribution Π(g

± j = k) = p B (p N B ) k-1 (k ∈ Z + ).
We claim that, for Π-a.e. ℓ, there exists K = K(ℓ) such that g ± j ≤ Kj. In fact, observe that (3) Π(g ± j ≥ j) = (p N B ) j-1 . This implies that the probabilities of the 'events' {g ± j ≥ j} form a summable sequence. Thus, by Borel-Cantelli, for a.e. ℓ, the inequality g ± j /j ≥ 1 is verified only for a finite number of j's. Setting K := max j (g ± j /j) proves the claim. To finish the proof of Corollary 4, observe that, for a non-degenerate Bernoulli measure, a.e. ℓ contains arbitrarily long blocking and nonblocking constant factors.

Lorentz gases

We now turn to the truly two-dimensional case. We tile the plane with Z 2 copies of the unit square C. In this case C is endowed with four gates: G 1 and G 2 , congruent and opposite open segments (say, the left and the right gate, respectively); G 3 and G 4 , again congruent and opposite open segments (the lower and the upper gates).

Apart from this, we have the same structure as in Section 2: a finite number of local configurations indexed by the set Ω = Ω B ∪Ω N B , where Ω B denotes the blocking and Ω N B the non-blocking configurations, as in Figs. 4 and5. Again, these sets verify five assumptions, (A1), (A3), (A4), (A5), and the counterpart of (A2): (A2 ′ ) ∂C \ (G 1 ∪ . . . ∪ G 4 ) ⊂ ∂Γ; and G 1 , . . . , G 4 do not intersect ∂Γ.

:= (ℓ n ) n∈Z 2 ∈ Ω Z 2 , with n := (n 1 , n 2 ), and the billiard table is (4) Q = Q ℓ := n∈Z 2 C n \ Γ ℓn n ,
where

C n := [n 1 , n 1 + 1] × [n 2 , n 2 + 1],
and Γ ℓn n is the configuration Γ ℓn translated to C n . We call Q (and the billiard map thereon) a Lorentz gas (LG)-see a realization in Fig. 6.

A word a 1 a 2 • • • a l is called a horizontal factor of ℓ ∈ Ω Z 2 if there exists n = (n 1 , n 2 ) such that ℓ n 1 +i,n 2 = a i , for i = 1, 2, . . . , l. The analogous definition is given for a vertical factor. As in Section 2, a factor is called non-blocking if a i ∈ Ω N B and constant if a i = a, for all i.

Since we are interested in infinite-horizon billiards, we discuss a sufficient condition for obtaining this property. It will be seen in Section 4 that the LGs that we construct possess arbitrarily long, both horizontal and vertical, non-blocking constant factors. They have infinite horizon if all non-blocking, say, horizontal factors of length l (∀l ∈ Z + ) admit a free flight of length ≥ l.

In the two-dimensional setting no one has succeeded in proving that recurrence is a typical property in a measure-theoretic sense (this is actually an important open problem, cf. [START_REF] Lenci | Typicality of recurrence for Lorentz gases, Ergodic Theory Dynam[END_REF][START_REF] Chernov | Hyperbolic billiards and statistical physics[END_REF]). The known results are with respect to a topological notion of typicality. For this, we make Ω Z 2 a metric space by endowing it with the distance

(5) dist(ℓ, ℓ ′ ) := n∈Z 2 2 -|n 1 |-|n 2 | |ℓ n -ℓ ′ n | .
It was shown by Lenci that, relative to the above metric, the Bairetypical LG in Ω Z 2 B is recurrent [START_REF] Lenci | Typicality of recurrence for Lorentz gases, Ergodic Theory Dynam[END_REF] and ergodic [START_REF] Lenci | Aperiodic Lorentz gas: recurrence and ergodicity, Ergodic Theory Dynam[END_REF]. On the other Figure 6. A Lorentz gas hand, Troubetzkoy has shown that, under the assumption that long factors of non-blocking configurations admit long free flights, the Bairetypical LG in Ω Z 2 is recurrent and has infinite horizon [START_REF] Troubetzkoy | Typical recurrence for the Ehrenfest wind-tree model[END_REF]. Once again, we combine and extend these ideas to show that a great number of infinite-horizon LGs are recurrent and chaotic in the sense of Theorem 2.

Theorem 5. Assuming (A1), (A2 ′ ), ( A3)-(A5), the metric space Ω Z 2 contains a dense uncountable set of LGs that are hyperbolic, recurrent and ergodic in the sense of Theorem 2, and such that the return map to any smooth component of the type γ i is K-mixing. Furthermore, if, for all l ∈ Z + , any non-blocking horizontal constant factor of length l admits a free flight of length ≥ l, then those LGs have infinite horizon as well.

Proofs

4.1. Sketch of the proof of Theorem 2. Let us fix an LT as in the statement of Theorem 2. Its recurrence is proved essentially in the same way as for the staircase billiards of [START_REF] Troubetzkoy | Recurrence in generic staircases[END_REF]. For the sake of completeness, we give an outline of the argument.

Here and in the remainder we are going to need a notation for the portion of the phase space pertaining to the cell C n :

(6)

M n := {(q, v) ∈ M | q ∈ C n } .
It is clear by the hypotheses that the LT contains, both forward and backwards, arbitrarily long constant factors of the same blocking configuration, say Γ 1 . In this paragraph and the next, when we say blocking factor, we will always mean a factor of configurations Γ 1 . Consider a blocking factor beginning with the cell C n 1 and ending with the cell C n 2 (n 2 > n 1 ). Let us call A n 2 n 1 the set of all the line elements of M n 1 whose trajectories visit a cell to the right of C n 2 , before coming back to C n 1 . Likewise, let A n 1 n 2 be the set of all the line elements of M n 2 whose trajectories visit a cell to the left of C n 1 , before coming back to C n 2 . Now fix k ∈ Z + . Since an LT made up entirely of cells Γ 1 is recurrent, there exists a positive integer l k such that, if the length of the factor (namely, n 2 -n 1 + 1) is bigger than or equal to l k , both µ(A n 2 n 1 ) and µ(A n 1 n 2 ) are smaller than 1/k. Consider a wandering set W and, for any n ∈ Z, set W n := W ∩ M n . W n is also a wandering set. For all k ∈ Z + , there exist a blocking factor of length ≥ l k (say, from the cell C n 1 to the cell C n 2 ) to the right of C n and a blocking factor of length ≥ l k (say, from C n 3 to C n 4 ) to the left of C n . By definition, the orbits of the points of W n are all disjoint and unbounded, which means they must intersect either A n 2 n 1 or A n 3 n 4 in distinct points. Therefore, using the invariance of the measure,

µ(W n ) ≤ µ(A n 2 n 1 ∪ A n 3 n 4 ) = 2/k.
Since k and n are arbitrary, W is a null set.

Turning to the hyperbolicity and ergodicity, once the recurrence is known, one proceeds as in [START_REF] Lenci | Aperiodic Lorentz gas: recurrence and ergodicity, Ergodic Theory Dynam[END_REF] or [CLS]. (A fairly accurate summary of the whole proof is given in [SLDC], although that article refers to LTs in dimension higher than two.) Here we limit ourselves to mentioning which parts of the proof can be worked out using the standard techniques for classical hyperbolic billiards [KS, LW] and which ones need to be adapted to our particular infinite-measure system.

Focusing on the hyperbolicity first, it can be seen that all the arguments used in the proof of Theorem 2(a), with one exception, are local, that is, depend on the value of T (its invariant cones, its distortion coefficients, etc.) on a neighborhood of a given point, or a given orbit. In other words, they cannot distinguish between a finite-measure dispersing billiard-for which everything works well [KS]-or an infinitemeasure one.

The only argument that is not local is the one whereby µ-almost all orbits stay sufficiently far away from the singular points of T (so that the singularities of the map do not interfere with the construction of the local stable and unstable manifolds). The singular points are organized in smooth curves, called singularity lines. Each such line corresponds to all the line elements whose first collision occurs at a given corner of ∂Q, or tangentially to a certain smooth component of it. Hence, there are a countable number of singularity lines. They can be counted (or at least overestimated) in a way that, in each region of the phase space M γ := {(q, v) ∈ M | q ∈ γ}, where γ is a smooth portion of the boundary as in (A1), there are at most two singularity lines for every source of singularity (a vertex or a tangency) "seen" from γ. (The factor 2 comes from the fact that there are two ways to be tangent to a smooth boundary, one for each orientation.) Also, the length of each singularity line is bounded above by a universal constant having to do with the size of the cell.

Let us indicate with S the singular set of T , that is, the union of all the singularity lines described above. It turns out that, if Lemma 6 below holds, almost all orbits approach S no faster than a negative power-law in time (where by time we mean 'number of collisions'). This is enough to make the sought argument work [START_REF] Lenci | Aperiodic Lorentz gas: recurrence and ergodicity, Ergodic Theory Dynam[END_REF][START_REF] Cristadoro | Recurrence for quenched random Lorentz tubes[END_REF].

Moving on to statement (c), as explained in [START_REF] Lenci | Aperiodic Lorentz gas: recurrence and ergodicity, Ergodic Theory Dynam[END_REF], one can exploit a suitable Local Ergodicity Theorem for billiards (say, the version of [LW]), which uses only local arguments except in its most delicate part, the so-called Tail Bound. It turns out, however, that a version of the Tail Bound can be proved for our LTs too, if the following result holds (cf. Sec. 3 of [SLDC]). Lemma 6. There exist constants C, α > 0 such that, for any n ∈ Z, any U ⊆ M n , and any t ∈ Z + , t j=0 T j (U) intersects at most Ct α singularity lines from S. (In more descriptive terms, each trajectory with initial conditions in U may approach at most Ct α singularity lines of T , in time t, and this bound is uniform in U, if U is not too large.)

Proof. Let us define f ± j := ± j i=0 g ± i , so that f + j is the location of the blocking cell to the right of the j th non-blocking factor, on the right side of the tube w.r.t. C 0 ; while f - j is the location of the blocking cell to the left of the j th non-blocking factor, on the left side of the tube; cf. definition (2). It follows from the hypotheses of Theorem 2 that both sequences (|f ± j |) j∈N are bounded by a power-law in j. Coming to the statement of the lemma, we may assume that U = M 0 (in fact, proving the result for U ′ := M n automatically proves it for U; then, choosing n = 0 is no loss of generality because a translation of the tube does not change our argument-cf. Remark 3).

By (A5), the configuration space trajectories of length t, with initial conditions in M 0 , cannot go further left than the cell C f - t or further right than C f + t . The corresponding phase space orbits are thus contained in (7)

A t := f + t n=f - t M n .
A line element in A t can "see" at most those cells that range from C f - t+1 to C f + t+1 . Therefore, the number of singularity lines of T in each set of the type M γ ⊂ A t (M γ was defined earlier) does not exceed C ′ (f + t+1f - t+1 + 1), for some C ′ > 0. But there are at most C ′′ (f + t -f - t + 1) sets of that type. The product of these two estimates, which, as shown earlier, grows no faster than a power-law in t, is an upper bound for the number of singularity lines from S in A t .

Once we have local ergodicity (all but countably many points in phase space have a neighborhood contained in one ergodic component), global ergodicity is easily shown. Also, the assertion about the Kmixing first-return map is proved as in [SLDC], Sec. 4.

L i := ℓ = (ℓ n ) ∈ Ω Z 2 ℓ n = 1, ∀n ∈ Z i .
In other words, L i is the set of all the global configurations which have "blocking circles" (namely, circles filled with cells of type Γ 1 ) at all radii j 2 , with j ≥ i. Clearly, L i ∼ = Ω Z i in a natural sense. Note that L i ⊂ L i+1 and that i L i is dense in Ω Z 2 . In each L i we apply the method of [START_REF] Troubetzkoy | Typical recurrence for the Ehrenfest wind-tree model[END_REF] to construct a G δ dense set R i of recurrent Lorentz gases. Let us sketch this method.

In what follows, whenever we mention circles, balls, annuli, we will always mean circles, balls, annuli in Z 2 , relative to the norm • and centered in the origin. Fix i and denote by ξ a configuration of cells in a ball of Z i , equivalently, a vector of Ω B , where B is the restriction to Z i of a ball in Z 2 ; B will be referred to as the support of ξ. Clearly, there are countably many such (finite) configurations, so we can index them as (ξ k ) k∈Z + (the dependence on i is suppressed). For each such k, let us construct a finite configuration η k such that:

• the support of η k is a ball of radius ρ > ρ 1 , where ρ 1 is radius of the support of ξ k , and the restriction of η k to the support of ξ k is ξ k ; • there exists a positive integer ρ 2 ∈ (ρ 1 , ρ) such that, if ρ 1 < n ≤ ρ 2 , ℓ n = 1, i.e., there is a "blocking annulus" (of type Γ 1 ) of radii ρ 1 , ρ 2 ; ρ 2 must be so large that the phase space measure of all the line elements based in the inner circle of the annulus, whose trajectories reach the outer circle before coming back to inner circle, does not exceed 1/k; • ρ -ρ 2 ≥ k and, for n ∈ Z i with ρ 2 < n ≤ ρ, ℓ n = m, i.e., the outer part of the configuration η k (except for the blocking circles D j , which do not belong to Z i ) is a non-blocking annulus (of type Γ m ) with thickness ≥ k. Then, let us denote by A k i the cylinder in L i defined by all the configurations in Z i that coincide with η k on its support. It is not hard to show that A k i is open in L i w.r.t. the metric (5). Hence (10)

R i := n∈Z + k≥n
A k i is a G δ set that is clearly dense in L i . The recurrence of any LG ℓ ∈ R i is proved essentially as in Section 4.1, by showing that, if W is a wandering set, then µ(W ∩ M n ) ≤ 1/k, for all n ∈ Z 2 and k large enough. Also, the construction of the non-blocking annulus in each η k and the hypothesis on the non-blocking horizontal constant factors (cf. Theorem 5) imply that ℓ has infinite horizon. The presence of the blocking circles D j is necessary to ensure that the set of all the trajectories with initial positions, say, in C n , stay confined, within time t ∈ N, to a portion of the LG that comprises at most Ct α cells, for some C, α > 0. This makes the equivalent of Lemma 6 hold, which in turn yields hyperbolicity, ergodicity and the other statements of Theorem 5.

Finally, R := i R i is a dense uncountable set of LGs in Ω Z 2 that possess all the sought properties.

Figure 1 .

 1 Figure 1. A blocking cell for the Lorentz tube

Figure 2 .

 2 Figure 2. A non-blocking cell for the Lorentz tube

Figure 4 .

 4 Figure 4. A blocking cell for the Lorentz gas

Figure 5 .

 5 Figure 5. A non-blocking cell for the Lorentz gas

4. 2 .

 2 Sketch of the proof of Theorem 5. Let us endow Z 2 with the norm n = (n 1 , n 2 ) := |n 1 | + |n 2 |. For j ∈ Z + , set D j := {n ∈ Z 2 | n = j 2 } (this set resembles the border of a rhombus in Z 2 ). Given i ∈ Z + ,

Centre de physique théorique, Federation de Recherches des Unites de Mathematique de Marseille, Institut de mathématiques de Luminy and, Université de la Méditerranée, Luminy, Case 907, F-13288 Marseille Cedex 9, France E-mail address: troubetz@iml.univ-mrs.fr URL: http://iml.univ-mrs.fr/ ~troubetz/

Acknowledgments. We thank A. J. Kox and J. Lebowitz for helping us with the history of the Lorentz gas, and two anonymous referees for their careful reading of the first version of this paper. M. L. is partially supported by the FIRB-"Futuro in Ricerca" Project RBFR08UH60 (MIUR, Italy). S. T. is partially supported by Projet ANR "Perturbations" (France).

Dipartimento di Matematica, Università di Bologna, Piazza di Porta San Donato 5, 40126 Bologna, Italy E-mail address: marco.lenci@unibo.it URL: http://www.dm.unibo.it/ ~lenci/