Francesco Di Maio

Piercesare Secchi

Simone Vantini

Enrico Zio
email: enrico.zio@polimi.it

FUZZY C-MEANS CLUSTERING OF SIGNAL FUNCTIONAL PRINCIPAL COMPONENTS FOR POST-PROCESSING DYNAMIC SCENARIOS OF A NUCLEAR POWER PLANT DIGITAL INSTRUMENTATION AND CONTROL SYSTEM

Keywords: dynamic reliability, fault scenarios classification, functional principal components analysis, fuzzy clustering, Fuzzy C-Means, digital instrumentation and control faults, nuclear power plants

This paper addresses the issue of the classification of accident scenarios generated in a dynamic safety and reliability analyses of a Nuclear Power Plant (NPP) equipped with a Digital Instrumentation and Control system (I&C). More specifically, the classification of the final state reached by the system at the end of an accident scenario is performed by Fuzzy C-Means clustering the Functional Principal Components (FPCs) of selected relevant process variables. The approach allows capturing the characteristics of the process evolution determined by the occurrence, timing, and magnitudes of the fault events. An illustrative case study is considered, regarding the fault scenarios of the digital I&C system of the Lead Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS). The results obtained are compared with those of the Kth Nearest Neighbor (KNN) and Classification and Regression Tree (CART) classifiers.

Introduction

Nuclear Power Plants (NPPs) are replacing and upgrading their aging and obsolete Instrumentation and Control (I&C) components, with a transition from analog to digital technology. Furthermore, new NPP designs involve an increased use of computer-based systems, expected to improve the plant's safety, reliability and failure detection capability. To assess the impact of digital I&C on NPPs safety, quantifiable reliable models are needed, along with data for digital systems that are compatible with the Probabilistic Risk Assessment Methodology (PRA). Due to the many unique attributes of these systems (e.g., software, fault tolerant features, different human-system interfaces, dynamic interaction between the plant system and processes), several challenges exist in modeling and data collection. The feature addressed in this paper is limited to the modeling of the dynamic interactions between the plant systems and processes, which in many instances regards both analog and digital I&C systems. Indeed, it has been shown that the order and timing of the fault events occurring along an accident sequence and the magnitude of the process variables at the time of event occurrence can be critical in determining the evolution of the accident and thus the risk associated with the system operation [Siu, 1994;[START_REF] Aldemir | A Benchmark Implementation of Two Dynamic Methodologies for the Reliability Modeling of Digital Instrumentation and Control Systems[END_REF]. In the case of digital I&C systems the modeling is rendered more difficult by the complex interactions of the software with the hardware and human components that are not easily captured by the existing PRA modeling tools, e.g., the Event-Tree (ET)/Fault-Tree (FT) approach which do not take explicitly into account the sequencing of the events nor their timing in the system stochastic evolution [START_REF] Rutt | Distributed Dynamic Event Tree Generation for Reliability and Risk Assessment[END_REF], Zio et al., 2009]. Methodological approaches have been studied and practical tools are being developed, in an effort to complement the current static PRA methods and tools. These so-called dynamic reliability [Devooght, 1997;[START_REF] Labeau | Dynamic Reliability: Towards an Integrated Platform for Probabilistic Risk Assessment[END_REF][START_REF] Dufour | Dynamic Reliability: A New Model, 13 th ESREL2002 European Conference[END_REF] or probabilistic dynamics [Devooght et al., 1992a;Devooght et al., 1992b] methods aim at giving explicit account to the interactions among the physical parameters of the process (temperature, pressure, speed, etc.), the human operators actions and the failures of the hardware and software components. Some extensions of classical methods have been proposed but in most instances these remain unable to take into account the system state dynamic changes and reconfigurations along the accident scenarios and their application is limited to problems with a small number of scenarios [START_REF] Medjoudj | A Method for Deriving Feared Scenarios in Hybrid Systems[END_REF]. Some methods allow a visual presentation of the sequence of events ordered in time, e.g. the Events Sequence Diagrams (ESD) [START_REF] Swaminathan | The Event Sequence Diagram Framework for Dynamic PRA[END_REF], Petri Nets [Peterson, 1977] and Dynamic Flowgraph Methodology (DFM) [START_REF] Guarro | Development of Tools for Safety Analysis of Control Software in Advanced Reactors[END_REF][START_REF] Garret | Automated Hazard Analysis of Digital Control Systems[END_REF]. Most of these methods can be regarded as semi-dynamic, because they represent system dynamics qualitatively (e.g., Petri Nets) or in a coarse partitioning of the system state space (i.e., in terms of large, small, medium changes in the controlled process variables such as in DFM). Discrete methods have also been developed to generate dynamic failure scenarios, e.g. the Discrete Dynamic Event Trees (DDET) [START_REF] Marchand | DDET and Monte Carlo Simulation to Solve Some Dynamic Reliability Problems[END_REF]], the DYnamic Logical Analytical Methodology (DYLAM) [START_REF] Amendola | DYLAM-1, A Software Package for Event Sequence and Consequence Spectrum Methodology[END_REF]Cojazzi, 1996] and the Dynamic Event Tree Analysis Method (DETAM) [START_REF] Acosta | Dynamic Event Trees in Accident Sequence Analysis: Application to Steam Generator Tube Rupture[END_REF]. The advantage of the dynamic event-tree generation techniques (such as DDET, DYLAM and DETAM) is that they are compatible with the existing PRA structure and are able to generate possible scenarios of the system evolution exhaustively. A major challenge of the dynamic reliability techniques is their computational complexity, both in model construction and implementation. Indeed, the number of dynamic scenario branches increases according to the power law with the number of occurring events and thus is much larger than in the classical FT/ET approach; the a posteriori information retrieval then becomes quite burdensome and difficult [START_REF] Labeau | Dynamic Reliability: Towards an Integrated Platform for Probabilistic Risk Assessment[END_REF].

The work in this paper reconsiders the problem of processing dynamic reliability scenarios as formulated in [START_REF] Zio | Processing Dynamic Scenarios from a Reliability Analysis of a Nuclear Power Plant Digital Instrumentation and Control System[END_REF], i.e., in terms of the classification of dynamic patterns on the basis of the similarity of both their stochastic features (i.e., times and magnitudes of failure event occurrences) and deterministic features (i.e., process variable values, such as temperatures and pressures). In general terms, the nature and complexity of the scenarios prevent the possibility of identifying crisp clusters of similarity in the features space (Figure 1); uncertainties and overlaps exist in the association of the accident sequences to the three end state classes. Such situation can be effectively tackled within a fuzzy clustering paradigm for classification [START_REF] Klir | Fuzzy Sets and Fuzzy Logic: Theory and Application[END_REF]. In Figures 2 and3, the logic framework and operative steps for the classification of the scenarios are sketched. Figure 2 (left) shows that the larger the number of failure modes, the larger the number of branches in a traditional ET, which can become infinite if considering the continuous stochastic times and magnitudes of the failure modes within a continuous ET framework [START_REF] Amendola | DYLAM-1, A Software Package for Event Sequence and Consequence Spectrum Methodology[END_REF]Smidts, 1994]; the dynamic scenario post-processing here proposed is aimed at classifying these branches into a small number of representative clusters of scenarios with common evolution features and leading to one of the pre-defined end states (right). The classification algorithm is built on a limited number of accident scenarios and can then be used to evaluate new accident scenarios for identifying their end states; this leads to savings in the computational time which makes feasible the dynamic analysis, irrespective of the number of system configurations to be considered.

In more details, the first step of the approach is the selection of the features of the scenarios relevant for the classification; the second step is the evolutionary training of the clustering algorithm on patterns of known classes. The latter one essentially involves the search for geometric clusters in the feature space based on a metric iteratively optimized so that the clusters finally obtained are close to the actual classes of scenarios of same end state [START_REF] Jain | Algorithms for Clustering Data[END_REF].

Process variables values

Stochastic variables values

End state 1

End state 2

End state 3 As in [START_REF] Zio | Processing Dynamic Scenarios from a Reliability Analysis of a Nuclear Power Plant Digital Instrumentation and Control System[END_REF], the supervised Fuzzy C-Means (FCM) clustering algorithm developed in [Bezdek, 1981;[START_REF] Zio | Identification of Nuclear Transients via Optimized Fuzzy Clustering[END_REF][START_REF] Mercurio | Identification and classification of dynamic event tree scenarios via possibilistic clustering: application to a steam generator tube rupture event[END_REF] is utilized as classifier of the accident scenarios. The original objective of the present work is the investigation of the feasibility of applying Functional Principal Component analysis for representing in a finite dimensional space with uncorrelated components, the characteristics of the accident scenarios that are typically nonfinite nor immediately describable by means of uncorrelated components [START_REF] Ramsay | Functional Data Analysis[END_REF].

As additional original contribution with respect to the previous work in [START_REF] Zio | Processing Dynamic Scenarios from a Reliability Analysis of a Nuclear Power Plant Digital Instrumentation and Control System[END_REF], the present study offers a comparison with a Kth Nearest Neighbor (KNN) [START_REF] Hastie | The elements of statistical learning. Data mining, inference, and prediction[END_REF] and a Classification and Regression Tree (CART) [START_REF] Breiman | Classification and regression trees[END_REF] algorithms for classification of the end state classes of the dynamic accident scenarios simulated.

The case study considered for illustration is also the same as in the previous cited work, i.e., the dynamic scenarios are taken from the analysis of the Lead Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS) equipped with digital I&C [START_REF] Cammi | Modelling and control strategy of the Italian LBE-XADS[END_REF]. This work focuses on scenario classification. For this reason, the analysis does not cover the study of the software and its possible failure modes, nor the benefit of fault tolerant features, nor the interactions of the software with the hardware and human components. The simplified dynamic scenario modeling is then intended to serve the purpose of showing the feasibility of effectively post-processing accident scenarios resorting to a FCM clustering algorithm based on the information provided by the FPCs. The classification results could then be used for the identification of the system prime implicants, i.e., the dynamic equivalent of the minimal cut sets of the system under analysis [START_REF] Aldemir | A Benchmark Implementation of Two Dynamic Methodologies for the Reliability Modeling of Digital Instrumentation and Control Systems[END_REF]. Moreover, it is worth pointing out that the actual implementation of the approach as a qualified tool in support of a complete quantitative dynamic reliability analysis and its inclusion in PRA would need to be supported by full dynamic accident calculations, inclusive of comprehensive models of hardware, software and human failure modes and their interactions.

The paper is organized as follows. In Section 2, the presentation made in [START_REF] Zio | Processing Dynamic Scenarios from a Reliability Analysis of a Nuclear Power Plant Digital Instrumentation and Control System[END_REF] of the mechanistic model used to generate the scenarios for the dynamic reliability analysis of the LBE-XADS is repeated, for completeness of the paper. In Section 3, some limitations of traditional static PRA modeling tools (ET/FT) are shown; these limits can be overcome by means of a classification method of dynamic scenarios based on FCM clustering of FPCs. In Section 4, the results of the application of the approach to the scenarios of the LBE-XADS are presented and supported with a comparison to the KNN and CART classifiers. Conclusions and remarks are given in Section 5.

Analysis of the LBE-XADS

The model

The Lead-Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS) is a subcritical, fast reactor in which the fission process for providing thermal power ( ) P t is sustained by an external neutron source through spallation reaction by a proton beam ( ) Q t accelerated by a synchrotron on a lead-bismuth eutectic target [START_REF] Bowman | Nuclear energy generation and waste transmutation using an accelerator-driven intense thermal neutron source[END_REF]Van Tuyle et al., 1993;[START_REF] Venneri | Accelerator-driven Transmutation of Waste (ATW) -A New Method for Reducing the Long-term Radioactivity of Commercial Nuclear Waste[END_REF][START_REF] Carminati | An Energy Amplifier for Cleaner and Inexhaustible Nuclear Energy Production Driven by a Particle Beam Accelerator[END_REF][START_REF] Rubbia | Conceptual Design of a Fast Neutron Operated High Power Energy Amplifier[END_REF]. In the current design [Ansaldo, 2001], the core contains Uranium and Plutonium dioxide fuel rods; future developments are aimed at housing also long-lived transuranic elements. A simplified scheme of the plant is sketched in Figure 4. The primary cooling system is of pooltype with Lead-Bismuth Eutectic (LBE) liquid metal coolant leaving the top of the core, at full power nominal conditions, at temperature A dedicated, dynamic simulation model has been implemented in SIMULINK for providing a simplified, lumped and zero-dimensional description of the coupled neutronic and thermo-hydraulic evolution of the system [START_REF] Cammi | Modelling and control strategy of the Italian LBE-XADS[END_REF]. The model allows the simulation of the system controlled dynamics as well as of the free dynamics when the control module is deactivated and the air cooler flow is kept constant. Both feedforward and feedback digital control schemes have been adopted for the operation of the system. The feedback controller is a PID (Proportional, Integral and Derivative)-based configuration with low values of both the proportional and the integral gains: it contributes for 70% of the load to fulfill the control strategy. On the other hand, the remaining 30% of the load is provided by the feedforward action which consists of a monotonically increasing function relating the required air mass flow rate to the reactor power. The logic of operation of the discrete-state control block can be represented as in Figure 5.

Figure 5 Block diagram representing the logic structure of the LBE-XADS control

The control is set to keep a steady state value of approximately 300 °C of the average temperature of the diathermic oil , av S o T : this value represents the optimal working point of the diathermic oil at the steady state, full nominal power of 80 MW th . On the contrary, an oil temperature beyond the upper threshold , th u o T =340 °C would lead to degradation of its physical and chemical properties, whereas a temperature below the lower threshold , th l o T =280 °C could result in thermal shocks for the primary fluid and, eventually, for the structural components [START_REF] Cammi | Modelling and control strategy of the Italian LBE-XADS[END_REF].

The Monte Carlo-driven fault injection engine

Multiple component failures can occur during the system life. To simulate this, the model has been embedded within a Monte Carlo (MC) sampling procedure for injecting faults at random times and of random magnitudes. This allows generating transients representative of the system dynamic accident scenarios.

The set of faults considered are (Figure 8): The combination of the selected faults gives rise to 64 possible system configurations depending on the number of failed components and the order of occurrence of the failures (Table 1). It is important to underline that the procedure implemented in this work for sampling the fault events is not intended to reproduce the actual stochastic failure behavior of the system components; rather, the choices and hypotheses for modeling the faults (i.e., the mission time, the number and typology of faults, the distributions of failure times and magnitudes) have been arbitrarily made with the aim of favoring multiple failures in the sequences and capturing the dynamic influence of their order, timing, and magnitude.

In particular, times and magnitudes of faults are obtained by a stratified sampling with respect to the 64 possible ordered accident sequences of Table 1, assumed equally probable. Within each stratum (i.e., given one of the 64 possible ordered sequences), the corresponding failure times are generated by means of a "stick-breaking" strategy [Halmos, 1944]: the first failure time is sampled from the uniform distribution [0,3000] [s] and the successive failure times from the conditional distributions, uniform from the last sampled time to 3000 [s]. Failure magnitudes are instead independently sampled from uniform distributions in [0,797], [0,1000] and [0,797] [kg/s] for m 1 , m 2 , and m 3, respectively. This sampling strategy models a wearing system, since the average failure rate is increasing along time. Moreover, even rare multiple fault events are included in the set of accident scenarios. Thus, it can be easily stated that the experimental set-up is a very conservative assumption and in this sense it further tests the robustness of the classification procedure. On one hand, in real applications these hypotheses can be relaxed, due to the high reliability of the equipment usually adopted in real NPPs; on the other hand, the complete analysis of a realistic system may involve dozens of components and failure modes, whose permutation of occurrence may lead to a combinatorial explosion of the system failure sequences; this, linked with the need of simulating a representative number of scenarios for each sequence, can make the analysis impracticable. To overcome this problem in practice, i) the system failure sequences to be considered will need to be identified by experienced-developed ETs which eliminate some permutations, but in a wise and controlled manner avoiding to ignore some relevant dynamic behaviors or collapsing together accident scenarios which are actually different and lead to different end states [START_REF] Aldemir | A Benchmark Implementation of Two Dynamic Methodologies for the Reliability Modeling of Digital Instrumentation and Control Systems[END_REF]; ii) the simulation of the required number of dynamic accident scenarios for the considered system failure sequences might need to be performed by fast-running empirical models providing sufficiently accurate response surfaces of the system dynamic behavior. The evolution of the accident scenarios may lead to three different end states, within the mission time of 3000 Table 1 Failure sequences: -= safe component, 1 st = first event of the sequence, 2 nd = second event of the sequence, 3 rd = third event of the sequence, 4 th = fourth event of the sequence

Actuators-PID communication interruption

1 1 st - - - 2 - 1 st - - 3 - - 1 st - 4 - - 1 st 5 1 st 2 nd - - 6 2 nd 1 st - - 7 1 st - 2 nd - 8 2 nd - 1 st - 9 1 st - - 2 nd 2 nd - - 1 st - 1 st 2 nd - - 2 nd 1 st - - 1 st - 2 nd - 2 nd - 1 st - - 1 st 2 nd - - 2 nd 1 st 1 st 2 nd 3 rd - 1 st 3 rd 2 nd - 2 nd 1 st 3 rd - 2 nd 3 rd 1 st - 3 rd 2 nd 1 st - 3 rd 1 st 2 nd - 1 st - 2 nd 3 rd 1 st - 3 rd 2 nd 2 nd - 1 st 3 rd 2 nd - 3 rd 1 st 3 rd - 2 nd 1 st 3 rd - 1 st 2 nd - 1 st 2 nd 3 rd - 1 st 3 rd 2 nd - 2 nd 1 st 3 rd - 2 nd 3 rd 1 st - 3 rd 2 nd 1 st - 3 rd 1 st 2 nd 1 st 2 nd - 3 rd 1 st 3 rd - 2 nd 2 nd 1 st - 3 rd 2 nd 3 rd - 1 st 3 rd 2 nd - 1 st 3 rd 1 st -

Limitation of traditional approaches

In order to show i) the limitation of traditional PRA tools in analyzing dynamic accident scenarios and ii) the feasibility of the approach here proposed for classifying dynamic accident scenarios, the system has been classically analyzed to identify its minimal cut sets with respect to the low-and high-temperature failure modes (the interested reader may refer to [START_REF] Zio | Processing Dynamic Scenarios from a Reliability Analysis of a Nuclear Power Plant Digital Instrumentation and Control System[END_REF]). To investigate this, the MC-driven fault injection engine introduced in Section 2.2 is used to sample for the 64 system configurations of Table 1, the time and magnitude of the components failures from their probability distributions. With the aim of generating a sufficiently wide variety of system dynamic behaviors, to guarantee the necessary coverage of system failure variability coverage, the number of simulated accident scenarios has been taken equal to 120 for each one of the 64 failure sequences; for each sequence, the random times and magnitudes of the components faults have been sampled from the respective assumed distributions and the SIMULINK model of the LBE-XADS in the sampled system configuration has been run. Figure 9, reports the frequency of the three system end states (high-temperature, safe, lowtemperature) for the 64 system configurations of Table 1. It is clear that the order, timing and magnitude of the events occurring along an accident scenario determine the evolution of the accident towards safe or fault end states [START_REF] Aldemir | A Benchmark Implementation of Two Dynamic Methodologies for the Reliability Modeling of Digital Instrumentation and Control Systems[END_REF][START_REF] Zio | Processing Dynamic Scenarios from a Reliability Analysis of a Nuclear Power Plant Digital Instrumentation and Control System[END_REF]. All but one configurations lead unequivocally to one and only one end state: on one side, this means that none of those system configurations is a minimal cut set of the system; on the other side, it also means that even when the order of the events in the sequence is accounted for this is not sufficient to unequivocally determine the consequent system end state, which depends also on the timing and magnitude of the occurring failures. The difficulties of the traditional approach to safety and reliability analysis appear clearly.

Processing dynamic scenarios for classification

The processing of the dynamic scenarios here performed aims at identifying classes of behavioural similarity in the scenario evolutions in terms of the characteristic features of the accident sequences and at relating them with the three possible end states of the system, which define the classes. The characteristics of these classes constitute information that could allow identifying the prime implicants of the system. The classification approach adopted is founded on a FCM clustering algorithm fed by FPCs of a selected safety variable (Appendix A at the end of the paper) [START_REF] Ramsay | Functional Data Analysis[END_REF] and optimized by mean of a supervised evolutionary search scheme (Appendix B at the end of the paper) [Bezdek, 1981;[START_REF] Zio | Identification of Nuclear Transients via Optimized Fuzzy Clustering[END_REF][START_REF] Mercurio | Identification and classification of dynamic event tree scenarios via possibilistic clustering: application to a steam generator tube rupture event[END_REF].

The supervised evolutionary clustering classifier based on functional principal components

In mathematical terms, the target of the supervised optimization is to find optimal Mahalanobis metrics which define c geometric clusters of the available data set which minimize the distance ( ) ( ) The overall iterative training scheme for the supervised optimization leading to the identification of the optimal Mahalanobis metrics and the corresponding clusters is given in Appendix B at the end of the paper. Note that the use of individual specific Mahalanobis metrics for defining the different clusters allows obtaining different ellipsoidal shapes and orientations of the clusters that can more adequately fit the actual data partition than the traditional Euclidean metric which leads to spherical clusters [START_REF] Yuan | Data driven identification of key variables[END_REF].

* 1 1 1 | | ( , ) ( , ) k k t t c c N i i t i i i i k x x D D c N c µ µ = = = - Γ Γ Γ Γ = = ⋅ ∑ ∑∑ (
When fed with a new pattern of feature values x of a given dynamic scenario, the trained classification algorithm provides the values of the membership functions ( )

* i x µ , 1 ,

LBE-XADS fault scenario classification by FCM

The three classes (c=3) which supervise the construction of the clusters in the feature space correspond to the three system end states which may arise in the dynamic scenarios generated by the accident sequences: Once constructed, the classifier can be used to classify any pattern of dynamic scenarios. In the present work, a total of M=1280 newly simulated accident scenarios of the three classes have been fed for classification by the fuzzy clustering classifier. These new scenarios are generated from the 64 possible accident sequences of Table 1 by sampling new realizations of the times of occurrence and magnitudes of the failure events (20 for each ordered sequence). The error rate of classification of the test patterns provided by the FCM classifier is summarized in Table 2: patterns are assigned to the class for which the membership value of the pattern is the highest. For comparison, the error rate when ) is considered as the most important feature for the post-processing clustering analysis is also reported as in [START_REF] Zio | Processing Dynamic Scenarios from a Reliability Analysis of a Nuclear Power Plant Digital Instrumentation and Control System[END_REF]. The analysis of the misclassification performance highlights the improved performance of the FCM based on FPC with respect to the methodology originally developed and presented in [START_REF] Zio | Processing Dynamic Scenarios from a Reliability Analysis of a Nuclear Power Plant Digital Instrumentation and Control System[END_REF]. The reason stands in the capability of FPCs in capturing relevant features of the stochastic sequences of failures, which are otherwise partially hidden when considering the safety relevant parameter , av C LB T and almost completely missing when analyzing it by resorting to static approaches (Figure 6).

Comparisons with KNN and CART classifiers

In this Section, the misclassification error rates of the test patterns associated to the FCM classifier are compared with those of the KNN and CART classifiers [START_REF] Hastie | The elements of statistical learning. Data mining, inference, and prediction[END_REF][START_REF] Breiman | Classification and regression trees[END_REF]. The KNN classifier is a non-parametric classification method; the central idea of this method is to determine the unknown class of a new pattern looking at the known classes of its neighbors. More specifically, the classification of a new pattern is done by assignment to the class whose frequency is the highest among the kth nearest neighbors of the new pattern, where the kth nearest neighbors are the k patterns of the training sample that are closest to the new pattern according to the Euclidean distance [START_REF] Hastie | The elements of statistical learning. Data mining, inference, and prediction[END_REF]. A CART classifier is essentially a recursive dicotomic partition of the input space in which each element of the partition (namely, the leaves of the tree) is associated to a class; classification of a new pattern is done by assignment to the class whose frequency is the highest among the patterns of the training sample within the leaf of the tree which the new pattern belongs to; maximization of the decrement of the Gini index is used as growing criterion [START_REF] Breiman | Classification and regression trees[END_REF], and minimization of the 10-fold crossvalidation misclassification error is used as pruning criterium [START_REF] Breiman | Classification and regression trees[END_REF]. In Table 3, the misclassification error rate on M=1280 newly simulated patterns are reported (columns refer to the different classifiers, rows refer to the different input vectors). As expected, because of the massive size of the training sample, the KNN classifier (that is fully non-parametric) performs similarly or better than both FCM and CART classifiers for all possible input vectors; this might not be confirmed for smaller training samples. Functional principal components appear to be effective not only in terms of descriptiveness (more of the 99% of the total variability is explained by the first three principal components, as shown in Figure 11), as expected, but also in terms of discriminating power. For instance, focusing on the FCM classifier, using the first three principal components of the average temperature of the LBE coolant (

Conclusions

Dynamic approaches to system safety and reliability analysis embed the physical models of process dynamics into the stochastic models governing the hardware, software, human components failure behavior. To capture the different system behaviors which arise, a large number of time-dependent scenarios of possible system evolutions are simulated and then interpreted for drawing conclusions on the system safety and reliability characteristics. This paper has addressed the problem of processing the numerous scenarios arising from a dynamic system safety and reliability analysis. With reference to a literature case study regarding the LBE-XADS system, a Fuzzy C-Means clustering approach has been presented for classifying scenarios with similar characteristics described by the Functional Principal Components (FPCs) of a selected safety parameter whose behavior is affected by the occurred events. The classification approach is based on the use of FPCs within an optimized fuzzy clustering scheme. The system end states of interest and the relative safety parameter are identified a priori by expert judgment on the basis of the analysis of the system design, logic and dynamics of the LBE-XADS system. In the case study presented, the safety parameter is the diathermic oil secondary coolant temperature which cannot exceed lower and upper thresholds otherwise the system enters low-and high-temperature failure modes, respectively. The application of the approach to the case study considered has demonstrated:

• the feasibility of the proposed approach of dynamic accident scenarios post-processing for retrieving safety-relevant information, and • the increased classification performance when relying on FPCs of the selected safety relevant parameters, with respect to other characteristic features. A comparison with the KNN and CART classification approaches has confirmed these findings. Future research could investigate to what extent the results are robust to changes in the Monte Carlo (MC) sampling procedure, and further develop the methodology aiming at the identification of the prime implicants of the system under analysis (dynamic analogue of minimal cut sets) to be included into PRA.

APPENDIX A

The basic concepts pertaining to Functional Principal Component Analysis (FPCA) as a tool for dimensional reduction are here presented. Let F be a random function such that

2 ( ) b a E F t dt   < +∞     ∫ (
Φ + = ∑ = µ
be the projection of F on the q-dimensional affine space centered on the mean function generated by the first q principal components. It can be proven that the approximation obtained by means of this projection is "statistically optimal". Indeed, identifying with Ρ the orthogonal projector on this space, it can be proven that:

( ) ( )       Ρ - =       - ∫ ∫ Ρ b a S b a q dt t F t F E dt t F t F E q q S 2 2 ) ( ) ( min arg ) ( ) ( (4)
Hence, the approximation obtained by representing F through its first q principal components minimizes the expected quadratic error over any other q-dimensional representation.

The effectiveness of this approximation can be easily quantified by means of the so called "fraction of explained total variance", i.e.:

( )

( ) ( ) ( ) ∑ ∑ ∫ ∫ ∫ ∫ ∞ = = =       -       - - =       -       - 1 1 2 2 2 2 ) ( ) ( ) ( ) ( 1 ) ( ) ( ) ( ) ( k k q k k b a b a q b a b a q dt t t F E dt t F t F E dt t t F E dt t t F E λ λ µ µ µ (5)
In the end, FPCA provides a useful tool to interpret and represent in a finite dimensional space with uncorrelated components, phenomena that are intrinsically non-finite and not immediately describable by means of uncorrelated components. matrices. Notice that the elements of the i -th Mahalanobis matrix are updated proportionally to the distance i D + between the i -th a priori known class and the i -th cluster found. In this way, only the matrices of those clusters which are not satisfactory for the classification purpose are modified. 6. Return to step 2.

Figure 1

 1 Figure 1 Schematics of the clustering of the accident sequences for classification of their end states, based on the combination of process and stochastic variables information

Figure 2 Figure 3

 23 Figure 2 Dynamic scenario classification into a low number of branches corresponding to the clusters of scenarios sharing common stochastic and deterministic evolution features which lead to one of the pre-defined end states

  τ equal to 400 °C pushed by natural circulation enhanced by argon gas injection into the heat exchangers of the secondary cooling circuit and then re-entering the core from the bottom through the down-comer at temperature , system is a flow of an organic diathermic oil at 290-320 °C, at full power conditions. Cooling of the diathermic oil is obtained through an air flow

Figure 4

 4 Figure 4 LBE-XADS simplified schematics. A = Accelerator; C = core; P = primary heat exchanger; S = secondary heat exchanger

  representing the SIMULINK model of the LBE-XADS is shown in Figure6: the controlled variable is the average temperature of diathermic oil ( , av S o T ), whereas the control variable is the mass flow rate of air ( a Γ ) in the air coolers battery. In Figure7, the profile of the average temperature of diathermic oil ( , av S o T ) at full power nominal conditions is shown: even if the system is stable at nominal conditions (303.85 °C), the discrete-state regulation of the air coolers causes visible ripples of the diathermic oil temperature.

Figure 6

 6 Figure 6 Block diagram representing the SIMULINK model of the LBE-XADS [Cammi et al., 2006]

Figure 7

 7 Figure 7 Diathermic oil temperature profile for the LBE-XADS in the stable, nominal conditions

Figure 8

 8 Figure 8 Sketch of the faults that can be injected into the system: the PID controller fails stuck at a random output value, the air coolers fails stuck at a random position, the feedforward control fails stuck at a random output value, the communication between air coolers actuators and the PID controller is interrupted

Figure 9

 9 Figure 9 Stacked bar chart of the frequency of end states for each of the 64 accident scenarios listed in Table 1

  this case c corresponds to the possible end states):

  the a priori known membership of the k -th pattern k x to the i -th physical class (possibly not known with absolute precision, in which case it has a membership less than one) the membership to the corresponding geometric cluster in the feature space.

class 1 :T

 1 low-temperature failure mode, with the safety parameter (i.e., the diathermic oil temperature , From the training sample, the main uncorrelated modes of variability in time (i.e., the FPCs) of the average temperature of the LBE coolant , of the underlying stochastic sequence of failures; on the other hand, they may lack a complete interpretation in reliability terms. In this case, the first principal time (Figure11); on the contrary, the interpretation of the second and third principal components are less clear.

Figure

  Figure 10 A sub-sample of sixty-four , av C LB T

  value over the transient duration of 3000 [s] of the average temperature of the LBE coolant , av C LB T

  is able to correctly classify nearly 19 patterns out of 20 (94.75%), while using the mean value of the average temperature of the LBE coolant , av C LB T it is able to correctly classify nearly 17 patterns out of 20 (83.35%).

•

  The PID controller fails stuck at time t 1 in [0,3000] [s], with a flow rate output of magnitude m 1 in [0,797] [kg/s]. • The air coolers fail stuck at time t 2 in [0,3000] [s] in a random position that provides a corresponding air flow mass m 2 in [0,1000] [kg/s]. • The feedforward controller fails stuck at time t 3 in [0,3000] [s] with a corresponding flow rate value m 3 in [0,797] [kg/s]. • The communication between air coolers actuators and PID controller fails at time t 4 in [0,3000] [s] so that the PID is provided with the same input value of the previous time step.

Table 2 Results of the classifier for different sets of FCP of the process variable

 2

	Misclassification error rate	FCM
	(	T	, av C LB	)					16.65%
	(	1 LB PC , , av C T	)				11.64%
	(	, , av C LB PC T	2	,	, , av C LB PC T	2	)	9.06%
	(	1 LB PC , , av C T	,	, , av C LB PC T	2	,	3 LB PC , , av C T	)	5.25%
										, av C LB T	and its mean value	T	, av C LB

Table 3 Results of the classifiers for different sets of FCP of the process variable

 3

	, av C LB T	and its mean value	T	, av C LB

  ; thus, if the sample size is sufficiently large, they can be used to estimate functional principal components.It can be proven that under assumption (2), the set { }, i.e. the space which the realizations of the random function F belong to. Thus, F can be decomposed as follows:The function k Φ is known as the kth functional principal component and k F as the score relative to the kth functional principal component. Moreover, it can be proven that

																																	2)
	with thus the kernel of its covariance operator. Often [ ] ) ( ) ( t F E t = µ being its mean function and estimated. Sample mean and sample covariance function provide consistent estimates of [ ] ) ( ) ( ) , ( s F t F E s t = its covariance function, and σ ) (t µ and ) , ( s t are not known and thus need to be σ ) (t and µ ) , ( s t σ
	covariance operator forms an orthonormal basis for	, 2 1 Φ [ ] Φ ( ) b a L ,	,...	of the eigenfunctions of the
						F	( t	)	=	µ	( t	)	+	∑	1 = ∞ k	F k Φ	k	( t	)	with	k F	=	∫ a b	(	F	t (	)	-	µ	t (	)	) Φ	k	t (	)	dt	(3)
																								[ ] k k F λ = var	with k λ being the eigenvalue associated to the
	eigenvector	k Φ (i.e.	k λ		Φ	k	(	s	)	=	∫ a b	σ	t (	,	s	)	Φ	k	t (	)	dt	) and	cor	[ F k F , k	'	]	=	δ	kk	'	(i.e. scores related to
	different functional principal components are uncorrelated scalar random variable). Note that ,... , , 3 2 1 Φ Φ Φ are ordered such that 0 ... 2 1 ≥ ≥ ≥ λ λ .
	Let	F	q	( t	)	( t	)			k q	1	F k			k	( t	)

2

APPENDIX B

The iterative training scheme for the supervised optimization leading to the identification of the optimal Mahalanobis metrics and the corresponding cluster can be summarized as follows:

1. At the first iteration ( 1 τ = ), initialize the metrics of all the c clusters to the Euclidean

, where I is the identity matrix.

2. At the generic iteration step τ , run the FCM clustering algorithm [Bezdek, 1981] Cholesky factors [Labeau, 1996],

+ is a lower triangular matrix with positive entries on the main diagonal. More precisely, at iteration τ , the entries ( ) ( )

, , max 10 , 0,


Enrico Zio
email: enrico.zio@ecp.fr

Giovanni Sansavini
email: giovanni.sansavini@polimi.it

Modeling Interdependent Network Systems for Identifying Cascade-Safe Operating Margins

Keywords: Critical infrastructures, interdependent systems, complex networks, failure cascade. ACRONYMS CIs critical infrastructures CI critical infrastructure

Infrastructure interdependency stems from the functional and logical relations among individual components in different distributed systems. To characterize the extent to which a contingency affecting an infrastructure is going to weaken, and possibly disrupt, the safe operation of an interconnected system, it is necessary to model the relations established through the connections linking the multiple components of the involved infrastructures. In this work, the modeling of interdependencies among network systems and of their effects on failure propagation is carried out within the simulation framework of a failure cascade process. The sensitivity of the critical loading value (the lower bound of the cascading failure region) and of the average cascade size with respect to the coupling parameters defining the interdependency strength is investigated as a means to arrive at the definition and prescription of cascade-safe operating margins.

N i number of nodes in system i P load transferred to neighboring components after failure in the two identically operating systems P i load transferred to neighboring components after failure in system i S average cascade size in the two identically operating systems S i average cascade size in system i S cr maximum allowable cascade size in the two identically operating systems S i cr maximum allowable cascade size in system i

I. INTRODUCTION

Understanding, modeling, and assessing the normal functioning and the possible faulty conditions of critical infrastructures (CIs) is essential to safely design these complex systems, and to effectively operate the services they provide. Distributed infrastructures have been modeled as graphs to unravel their structure and dynamics, and to evaluate to what extent the structure may impact the dynamics [START_REF] Strogatz | Exploring complex networks[END_REF]. However, it is evident that infrastructures do not exist in isolation of one another, and the relations among them must be identified to perform realistic and applicable analyses [START_REF] Rinaldi | Modeling and simulating critical infrastructures and their interdependencies[END_REF]. The focus of the research on critical infrastructures must then be shifted from single, isolated systems to multiple, interconnected, mutually dependent systems. The aim becomes that of assessing the influences and limitations which interacting infrastructures impose on the individual system operating conditions, for avoiding fault propagation by designing redundancies and alternative modes of operations, and by detecting and recognizing threats [START_REF] Zimmerman | Social implications of infrastructure network interactions[END_REF].

In modern society, the linking among service infrastructures is required for optimal and economical operation.

Yet, these interconnections introduce weaknesses in the systems due to the fact that failures may cascade from one system to other interdependent systems until, possibly, affecting the overall functioning, if proper protection of interdependencies is not considered [START_REF] Little | Controlling cascading failure: understanding the vulnerabilities of interconnected infrastructures[END_REF]. The role of dependencies among infrastructures (so called interdependencies), and the intrinsic difficulties arising in their modeling, have been highlighted in empirical studies. For example, the database in [START_REF] Luiijf | Empirical findings on critical infrastructure dependencies in Europe[END_REF] has been built from public reports of disruptions of CIs from open sources like newspapers and internet news. Events have been classified as "cascade initiating" (i.e., an event that causes an event in another CI), "cascade resulting" (i.e., an event that results from an event in another CI), and "independent" (i.e., an event that is neither a cascade initiating nor a cascade resulting event). The information in the database shows in particular that: i) "cascade resulting" events are more frequent than generally believed, and that "cascade initiators" are about half as frequent; ii) the dependencies are more focused and directional than often thought; and

iii) energy and telecommunication are the main "cascading initiating" sectors.

Interdependencies among telecommunication systems, transportation systems, and power distribution grids played a negative role in the small telecommunication blackout which took place in the suburbs of Rome in January 2004 [START_REF]Tools and techniques for interdependency analysis[END_REF]. The chain of events originated in a major telecommunication service node when a metallic pipe carrying cooling water for the air conditioning broke. The resulting flood led to a communication blackout in the area, which in turn hit the country's biggest printed news agency transmissions; stopped the check-in, ticketing and luggage acceptance at the International Fiumicino Airport; disturbed post offices and banks operations; and caused problems to the communication network connecting the main Italian research institutions. Furthermore, the telecommunication blackout had effects on the operation of the power distribution grid causing the disconnection of two control centers, and resulting in the loss of control on a number of remote substations in the area of Rome.

Similar occurrences [START_REF] Zimmerman | Social implications of infrastructure network interactions[END_REF] point to the need for identifying infrastructure interdependencies [START_REF] Rinaldi | Identifying, understanding and analyzing critical infrastructures interdependencies[END_REF], and determining the design and operational conditions which can prevent the onset and propagation of cascading failures.

Due to the highlighted intrinsic complexities related to interdependencies among CIs, their modeling is far from trivial, and abstractions can be of help for preliminary investigation and broad understanding.

Many models and computer simulation techniques exist for analyzing the behavior of individual infrastructures (e.g., load flow and stability codes for electric power networks, connectivity and hydraulic analyses for pipeline systems, traffic management models for transportation networks), but modeling and simulation frameworks that allow the coupling of multiple interdependent infrastructures are needed to address infrastructure protection, mitigation, response, and recovery issues. In developing these, it is important to know that simply linking together existing infrastructure models into a reduced model may fail to capture the emergent behavior arising in interdependent infrastructures [START_REF] Gautam | Multi-agent based models of social contagion and emergent collective behavior[END_REF].

In the last ten years, a number of researchers have focused on the study of interdependent infrastructures, and developed a variety of methodologies to be applied for ensuring the safe, reliable, continuous operation of interdependent infrastructures. According to [START_REF] Rinaldi | Modeling and simulating critical infrastructures and their interdependencies[END_REF], these modeling and simulation techniques can be grouped into six broad categories: (i) aggregate supply and demand tools, which evaluate the total demand for infrastructure services in a region and the ability to supply those services [START_REF] Min | Toward modeling and simulation of critical national infrastructure interdependencies[END_REF]- [START_REF] Helseth | Structural vulnerability of energy distribution systems: incorporating infrastructural dependencies[END_REF]; (ii) dynamic simulations, which examine infrastructure operations, the effects of disruptions, and the associated downstream consequences [START_REF] Newman | Risk assessment in complex interacting infrastructure systems[END_REF]- [START_REF] Ouyang | A methodological approach to analyze vulnerability of interdependent infrastructures[END_REF]; (iii) agent-based models, which allow the analysis of the operational characteristics and physical states of infrastructures [START_REF] Barton | An agent-based microsimulation of critical infrastructure systems[END_REF]- [START_REF] Schläpfer | Reliability analysis of electric power systems using an object-oriented hybrid modeling approach[END_REF];

(iv) physics-based models, which analyze physical aspects of infrastructures with standard engineering techniques (e.g., power flow and stability analyses for electric power grids or hydraulic analyses on pipeline systems) [START_REF] Chen | Identifying high risk N-k contingencies for online security assessment[END_REF], [START_REF] Kodsi | Application of a stability-constrained optimal power flow to tuning of oscillation controls in competitive electricity markets[END_REF];

(v) population mobility models, which examine the movement of entities through geographical regions [START_REF] Hong | A group mobility model for ad hoc wireless networks[END_REF]- [START_REF] Casalicchio | MobileOnRealEnvironment-GIS: a federated mobile network simulator of mobile nodes on real geographic data[END_REF];

and (vi) Leontief input-output models, which in the basic case provide a linear, aggregated, time-independent analysis of the generation, flow, and consumption of various commodities among infrastructure sectors [START_REF] Haimes | Leontief-based model of risk in complex interconnected infrastructures[END_REF]- [START_REF] Cagno | Risk analysis of underground infrastructures in urban area: time-dependent interoperability analysis[END_REF].

The present work falls in the dynamic simulations category, and aims at developing a simulation framework which abstracts the physical details of the services provided by the infrastructures while at the same time capturing their essential operating features. Such type of simulation models is appealing because of its simplicity and feasibility of use, which allows identifying the factors which most contribute to the safe, efficient design and operation of interdependent infrastructures. This knowledge is achieved through the assessment of the extent to which the presence of interdependencies affects the performance of the individual infrastructures. In particular, this scheme of analysis can be used for a first evaluation of the operational margins of interdependent systems which ensure the services provided by the involved infrastructures with cascade-safe conditions and prevent the propagation of failures over large areas.

A model of cascading failures in interdependent network systems is developed as an extension of the work in [38] and [START_REF] Newman | Risk assessment in complex interacting infrastructure systems[END_REF], considering the local propagation of the overload originated from a failure to first-neighbors, and to the interdependent set of the failed component. The proposed model differs from similar studies such as [START_REF] Motter | Cascade-based attacks on complex networks[END_REF] because it does not rely on the assumption that the load on a component is proportional to the number of shortest paths running through it, which do not always reflect the actual routes of the flow in a network; on the contrary, it shares similarity with the fiber bundle model in complex networks which has been applied to the blackout scenarios for cascading breakdown of power plants from overloading [START_REF] Kim | Universality class of fiber bundle model on complex networks[END_REF]. Similar to the fiber bundle model, if the load in a node exceeds its capacity, the node gets disconnected from the network, and the demand is transferred locally to neighboring nodes through the edges of the network. Moreover, it is assumed that the network connections cannot fail while the vertices are damaged and removed from the systems.

The current model is applied on a modified literature case study [START_REF] Grigg | The IEEE reliability test system -1996[END_REF] with the aim of identifying the interdependency features most critical for the cascade-safe operations of interdependent infrastructures, and defining related cascade-safe operational margins. The findings are in agreement with those of [START_REF] Carreras | Interdependent risk in interacting infrastructure systems[END_REF], in which a parallel study on cascading models in interdependent networks is presented, based on the forest fire model [START_REF] Bak | A forest-fire model and some thoughts on turbulence[END_REF], and mean-field theory analysis [START_REF] Grassberger | On a forest fire model with supposed self-organized criticality[END_REF]. A percolation threshold [START_REF] Stauffer | Introduction to percolation theory[END_REF] is identified above which cascading failures of all sizes can propagate in the interdependent systems, and a parameter is found to characterize the cascading threshold in the coupled systems. Following [START_REF] Carreras | Interdependent risk in interacting infrastructure systems[END_REF], the present work adopts the loading level in the interdependent systems as a metric for practical application to characterize the percolation threshold in more realistic systems.

The paper is structured as follows. The research framework and the motivations for the proposed analysis are stated in Section I. The modeling assumptions of interdependencies among network systems are detailed in Section II. The resulting algorithm for modeling cascading failures in interdependent networks is illustrated in Section III. In Section IV, a crude sensitivity analysis is performed with respect to the main model parameters related to interdependency, and conclusions are drawn in Section V.

II. INTERDEPENDENT NETWORK SYSTEMS

The modeling development carried out in this work considers two interdependent network systems, each one made by a number of components (nodes) connected by links (arcs) representing the physical and/or logical relations among them. The interdependencies are modeled as links connecting nodes of the two systems; these links are conceptually similar to those of the individual systems, and are bidirectional with respect to the "flow" between the two interdependent networks.

The number of interdependency links among interdependent components in the two systems, M, and the fixed amount of load transferred over the interdependency links upon a component failure, I, are essential features characterizing the "coupling energy" between the two interdependent systems.

In the following analysis, the interest is on systems with fixed "interdependency energy", i.e., such that each node in system 1 could be interdependent with any other node in system 2, but M is fixed. Communication systems, in which each node in system 1 can interact with any other node in system 2 but there is a maximum amount of connecting energy between the two systems, could be an example of such energy-limited systems.

To account for the dynamics of changing connections between the two network systems under developing failure cascade processes, Monte Carlo simulations are performed in which M, and I are kept constant, but the interdependency connections among components are randomly rewired at each Monte Carlo trial. This approach allows evaluating the average systems responses to cascading failures because it accounts for multiple connectivity patterns among the M interdependent nodes in the two systems.

Fig. 1 shows the systems considered. They are two identical networks which are an abstract topological model of the IEEE Reliability Test System -96 [START_REF] Grigg | The IEEE reliability test system -1996[END_REF]. M interdependency links are drawn between them as explained above (dashed lines in Fig. 1). In the proposed analysis, interest is on cascade onset and propagation over the bare topological structure of the test systems; no reference is made to the specific electrical properties which characterizes these electrical infrastructures. The topology corresponds to the IEEE Reliability Test System -96; however, no other aspect of this paper is related to the real network. For simplicity, but with no loss of generality, all components in the i-th system are assumed to have the same limit of operation L i fail , beyond which they are failed (L 1 fail = L 2 fail = L fail =1 in this study, upon normalization of all loads relative to the failure load value). In the reference cascading failure model of [START_REF] Newman | Risk assessment in complex interacting infrastructure systems[END_REF], when a component fails, a fixed and positive amount of load is transferred to each of the system's components; on the contrary, in the proposed extended model, the overload is propagated locally, to the first-neighbors of the failed node within the network structure it belongs to (a fixed and positive amount of load, P i , i = 1, 2; P 1 = P 2 = P = 0.07 in this study), and to the interdependent components which the failed component is connected to in the other network system (a fixed and positive amount of load, I), if any. If there is no working node in the neighborhood of a failed component or among the interdependent nodes connected to it, the cascade spreading along that path is stopped. The case of two fully connected systems, where all nodes are first-neighbors and every component in system 1 is interdependent to every component in system 2, coincides with the original model proposed in [START_REF] Newman | Risk assessment in complex interacting infrastructure systems[END_REF].

The interdependency links between the two network systems are treated in the same way as the individual system links. They are bidirectional connections and upon the failure of a node in system 1 or 2, the fixed amount of load, I, is propagated locally to the nodes in the interdependent network system 2 or 1, if any interdependency is present for the failed node. This transfer accomplishes the coupling between the two systems.

To start the cascade, an initial disturbance imposes an additional load The algorithm for simulating the cascading failures proceeds in successive stages as follows:

0 At stage m = 0, all N 1 +N 2 components j i = 1, 2,…, N i , i = 1, 2
in systems 1 and 2 are working under independent uniformly random initial loads

i j L  [L i min , L i max ], with L i max  L i fail , j i = 1, 2,…, N i , i = 1, 2.
1 M interdependency links between system 1 and system 2 are generated, connecting two randomly selected components; multiple interdependency links are allowed for each component.

2 An initial disturbance i j D , j i = 1, 2,…, N i , i = 1, 2
, is added to the load of each component in the two systems. 3 Each component that has not failed is tested for failure: for j i = 1, ..., N i , if component j i has not failed, and its load i j L > L i fail , then component j i in system i fails, i = 1, 2. 4 The component loads are incremented, taking into account the network topology, i.e. the failed component neighborhood: for each failed node in system i, the load of its first-neighbors is incremented by an amount P i , i = 1, 2. If the working neighborhood set of the failed node is empty, the associated failure propagation into the system comes to an end. [START_REF] Luiijf | Empirical findings on critical infrastructure dependencies in Europe[END_REF] The component loads are incremented, taking into account the interdependency pattern, i.e., the nodes interdependent to the failed component: for each failed node in system 1 or 2, the load of its interdependent nodes in system 2 or 1 is incremented by an amount I. If the interdependency set of the working components of the failed node is empty, the associated failure propagation to the interdependent system comes to an end. [START_REF]Tools and techniques for interdependency analysis[END_REF] The stage counter m is incremented by 1, and the algorithm is returned to step 3.

The algorithm stops when failures are no further propagated within or between the two systems.

Various initial system loading levels are evaluated, varying the uniform sampling ranges [L i min , L i max ], i = 1, 2, whose midpoints, L i , are indicators of the average initial systems loading levels. Large L i values relate to operating conditions in which the systems are more stressed. In this study, when the systems are operating at varying average initial load, L i , its range of variation is [0.5, 1] at steps of 0.005.

IV. SENSITIVITY ANALYSIS

The effects of the interdependencies between the two systems are shown in Fig. 2 in terms of the average cascade size, S i , i.e., the number of failed components in the i-th system at the end of the cascade spread, versus L i , which represents the system operating level, i = 1, 2. For each value of L i , varying in the range [0.5, 1] at steps of 0.005, 100,000 Monte Carlo simulations have been repeated. In these simulations, L i fail =1, i j D =2%, P i =I=7%, K i =M=34 for

j i = 1, 2,…, N i , i = 1, 2.
Each simulation corresponds to a different sampled pattern of the M interdependency links. The triangles in Fig. 2 represent S i in system i = 1, 2, as a function of L i =L, for the isolated single system i = 1, 2, i.e., when no interdependency is present. The overlapping squares and circles represent the same quantity S 1 = S 2 for the identical and identically operating systems 1, and 2, respectively.

As expected, the interdependencies cause a shift to lower values of the loading threshold for which the cascading phenomenon starts appearing with significance (from approximately 0.9 for the individual isolated system to approximately 0.8 for the interdependent systems). As L i =L increases, the systems are increasingly vulnerable to cascading failures. The transition, to the region where cascade of significant size appear, is smoother for the two interdependent networks than for the individual isolated system. This result is due to the fact that cascades appear at lower L i =L in the interdependent networks, which thus are less stressed, and prone to their propagation. As found in [START_REF] Carreras | Interdependent risk in interacting infrastructure systems[END_REF], the interdependencies between the two networks make access to the critical point easier, which means that the systems when coupled are more susceptible to large-scale failures, and a failure in one system can cause a similar failure in the coupled system.

To quantitatively assess the effects of the interdependency in Fig. 2, a threshold representing the maximum allowable cascade size, S i cr , in system i = 1, 2, can be set. The intersection with the horizontal line, S i = S i cr , identifies the critical load, L i cr , beyond which the threshold is exceeded in system i = 1, 2. S i cr is interpreted as the maximum number of components which can be lost in system i without affecting the global service provided by the infrastructure. This threshold can vary from system to system, and is a distinguishing feature of the provided service.

In the following, for simplicity but with no loss of generality, S 1 cr = S 2 cr = S cr = 15% is assumed, which identifies L 1 cr = L 2 cr = L cr = 0.8662 for the individual systems in isolated conditions and L 1 cr = L 2 cr = L cr = 0.7266, for the two interdependent systems (Fig. 2).

L i

cr is a relevant feature of a network system because it identifies, together with the continuous change in gradient, a type-two transition [START_REF] Huang | Statistical mechanics[END_REF] between the cascade-safe region and the onset of disrupting cascades in terms of the loading conditions, L i . Along with S i , it gives essential information on the system vulnerability towards cascading failures, and it can help identify safety margins of system operation. On the other hand, see from Fig. 2 that S i has small values until it reaches the proximity of L i cr , beyond which the transition to cascade region occurs. Thus, its use as a vulnerability indicator needs to be gauged against the sensitivity of detecting small changes in its value. Then, appropriate countermeasures of reducing the working load for limiting the failure propagation can be taken as the operation conditions are approaching L i cr . From Fig. 2, observe that the transition from cascade-safe to cascade regime is smoother for interdependent networks, allowing more time for reacting against the cascading failure onset. Fig. 3 shows the probability distribution function for S 1 in the interdependent system 1, when L 1 = L 1 cr . The point is characterized by the L 1 cr value at which the percolation transition occurs [START_REF] Stauffer | Introduction to percolation theory[END_REF], and cascades of all sizes can propagate in the interdependent systems. The power-law tail for large S 1 values developed at the critical point L 1 cr with an exponent of approximately -0.4 is consistent with the forest tree model analyzed in [START_REF] Carreras | Interdependent risk in interacting infrastructure systems[END_REF]. It is important to stress that even at the critical point L 1 cr , which marks the upper bound of the cascade-safe region, cascades of significant size (S 1 =22, i.e. 92% of the system size N 1 =24) appear, even if their frequency is extremely low. Once more, this finding points to the fact that the use of L i cr as a vulnerability indicator needs to be gauged against the sensitivity of detecting small changes in its value. Then, appropriate countermeasures of reducing the working load for limiting the failure propagation can be taken as the operation conditions are approaching L i cr .

Fig. 3. The probability distribution of the cascade size S i at the critical load L i cr . In these simulations L i fail =1, i j D =2%,

P i =I=7%, K i =M=34 for i=1, 2.
To understand the effects on the cascade process of the parameters characterizing the interdependency between the two systems, we conduct a further sensitivity analysis in several operating conditions which reflect real system operations.

A. Sensitivity with respect to the load of interdependent system 2

The following analysis aims at assessing the extent to which an interdependent system working at different, fixed load levels influences the coupled network system with respect to its vulnerability towards cascading failures. To this aim, the variation of L 1 cr is assessed while system 2 is working at fixed constant loads; the analysis is performed crudely for fixed values of average initial loads of system 2 ranging between L 2 =0.5 and L 2 =1, in steps of 0.05. Fig. 4 shows the results of this analysis. As expected, the coupling between the two systems is such to weaken the resistance of system 1 to failure cascade, forcing it to be operated at increasingly lower levels as L 2 increases. The emerging functional dependence, however, could not be easily anticipated. L 1 cr decreases linearly for system 2 loading levels below L 2 =0.75, and suddenly drops when system 2 loading levels rise above L 2 =0.75, until it saturates for system 2 average initial loads beyond L 2 =0.9. This behavior indicates that, when L 2 0.9, system 1 experiences unbearable cascades irrespective of its loading level. There is no cascade-safe region for system 1 when system 2 is operating beyond 90% of L 2 fail . The trend found indicates that the coupling between the two interdependent systems is such that under given loading conditions and beyond certain thresholds, the effects of the system nonlinearities become relevant, and an emergent behavior arises in the interdependent systems [START_REF] Rinaldi | Identifying, understanding and analyzing critical infrastructures interdependencies[END_REF]. This result is in agreement with [START_REF] Carreras | Interdependent risk in interacting infrastructure systems[END_REF], in which it was found that symmetric coupling of network systems actually decreases nonlinearly the percolation threshold above which cascading failures of all sizes can propagate in the interdependent systems. In the present model, the percolation threshold for the propagation of cascading failures is given by L 1 cr , which has a nonlinear relation with the load of the interdependent system, L 2 . Fig. 4. Critical load, L 1 cr , in system 1 for constant average initial load levels, L 2 , of system 2.

This result points to the fact that the loading level of the interdependent system plays a role of paramount importance in the failure behavior of its companion coupled system because it affects the cascading failure behavior of the latter in an unpredictable way. When designing and operating interdependent infrastructures, it is then necessary to control the operating levels of the systems, and assess the values beyond which nonlinearities start governing the system cascading failure behavior.

B. Sensitivity with respect to the number of interdependency links, M

In the following analysis, the effects of M on the vulnerability to cascading failures are assessed in two different system operating conditions.

In the first case, both systems are operating at the same varying L 1 = L 2 = L[0.5, 1] at steps of 0.005, and the variation of L 1 cr is analyzed with respect to M (Fig. 5). L 1 cr is assessed starting the cascade simulation from L=0.5, and increasing this value until the onset of the cascade regime is encountered. Because both systems are identical and operate at the same loading conditions, they will show identical trends of L 1 cr = L 2 cr = L cr , similarly to what is observed in Fig. 2. See from Fig. 5 that there is an approximately linear functional relationship between L cr and M, up to the value M=70 for which the system cascade-safe region disappears. If M  70, the systems are going to experience unbearable cascades irrespective of the loading level. There is no cascade-safe region when more than M=70 interdependency links are present between the two systems. Thus, if one were to try to protect the interdependent systems from cascade failure by controlling M, it appears that nonlinearities do not play a significant role. A linear decrease of the cascade-safe region is to be characterized with respect to the addition of interdependency links between the two systems instead.

The characterization of this relationship is relevant in the definition of cascade-safe operating regimes for the interdependent systems: for a fixed M, a critical loading level can be identified below which the systems can be safely operated. In the present example, it turns out that there is no safety margin when more than M=70 interdependency links are present between the two systems, which is more than twice the number of links in each system, K i =34, i =1, 2. This means that beyond M=70, cascades of significant size spread in the network irrespective of L. In the second case, both systems are operating at the same constant L 1 = L 2 = L = 75% of L 1 fail = L 2 fail = L fail , and S 1 = S 2 = S is assessed with respect to the variation of M (Fig. 6).

As expected, S 1 = S 2 = S increases as M increases until a saturation value is reached, which is a function of the load transferred over the interdependency links, I, and the constant L.

Cascade-safe operating regimes can be identified with respect to S; once the operating level is known (L = 75% in this case), the systems can be operated or designed to limit the maximum S. As an example, from Fig. 6 see that to have cascades involving less than 15% of the system components, no more than M=29 interdependency links can be operated between the two systems.

C. Sensitivity with respect to load transferred over the interdependency links, I

In the following analysis, the effects of the load transferred upon failure over the interdependency links of the failed component, are assessed in the same two system operating conditions introduced in Section B.

In the first case of both systems operating at the same varying average initial load, the variation of L cr is investigated as a function of I (Fig. 7). Similar to the findings of the previous analysis, an approximately linear functional relationship emerges. The system cascade-safe region is dropping slowly, and only when I=25% of L fail is transferred over each interdependency link, cascades of significant size spread in the network irrespective of L. This slow drop points to the fact that I is a less critical parameter in designing and operating interdependent network systems against cascading failures.

As before, cascade-safe operating regimes for the systems can be identified with respect to L cr : for a fixed value of I, a critical loading level can be determined below which the systems can be safely operated. In the present example (Fig. 7), if I=9%, operating the systems beyond L cr =70% of L fail , could result in cascading sizes larger than the maximum allowable size (15% of the system size, N=24). In the second case of both systems operating at the same constant L 1 = L 2 = L = 75% of L 1 fail = L 2 fail = L fail = 1, S 1 = S 2 = S is assessed with respect to the variation of I (Fig. 8).

As expected, S increases as I increases until a saturation value is reached, which is a function of M and the constant L. Cascade-safe operating regimes can be identified with respect to S: given the operating level (L = 75% in this case), the systems can be operated or designed to limit the maximum S. As an example, from Fig. 8 it can be understood that to have cascades involving less than 15% of the system components, no more than I=5.4% of L fail should be allowed to flow over each interdependency link connecting the two systems.

V. CONCLUSION

An abstract model of cascading failures in interdependent network systems has been developed. Its application on a modified literature case study network structure has shown the usefulness of this kind of models for identifying the relevant factors affecting the failure cascade process, and for capturing their relationships with the effects of the cascade.

The simulation results show that the critical load, L cr , and the average cascade size, S, give essential information on the system vulnerability towards cascading failures, and help identify cascade-safe regions for the system operations.

To investigate the effects of the parameters characterizing the interdependency between the two systems, a sensitivity analysis has been performed in several operating conditions which reflect real system operations.

The analysis has revealed that the interdependent system loading level is a vital feature to control because it affects the interdependent systems cascading failure behavior in an unpredictable way. In particular, when designing and operating interdependent infrastructures, it is important to determine the threshold level of loading beyond which nonlinearities emerge in the system cascading failure behavior.

The results of analyses of this kind allow defining cascade-safe operating regimes with respect to L cr for the network systems: given the number of interdependency links in the system, M, and the load transferred over the interdependency links, I, loading levels can be identified below which the systems can be safely operated.

Furthermore, cascade-safe operating regimes can be identified with respect to S. Once the operating level is known, the systems can be designed or operated, tweaking M and/or I, to limit the maximum S.

An interesting direction currently under investigation is the extension to coupled infrastructures of an analytical model for computing thresholds of criticality. The thresholds are functions of the average degree of the network, and of the amount of load transferred in the uniform network base case.

As a concluding remark, a recent NSF workshop report points at the fact that uncertainty is pervasive in complex systems [START_REF] Guckenheimer | Foundations for complex systems research in the physical sciences and engineering[END_REF]. Thus, developing reliable predictions about system behavior in the face of the large numbers of uncertain parameters in models of actual complex systems is a major challenge. Quantifying this uncertainty, and determining how it propagates throughout the system is a key aspect of reliable prediction and control of cascading failure in critical infrastructures. Indeed, the choice of the model parameters can be critical for the system response, and for the full representation of the variability associated with the system response, including that due to the uncertainty in the model parameters (e.g. I). Therefore, the evaluation of the cascade-safe operating margins will be complemented with a full quantification of uncertainty in a following study. This future study will add a further level of confidence to the findings of the sensitivity analysis.

  the components in the two identically operating systems i j D initial disturbance on component j i in system i I amount of load transferred over the interdependency links K number of links in the two identical systems K i number of links in system i L average initial load in the two identically operating systems L i average initial load in system i L cr critical load in the two identically operating systems L i cr critical load in system i L fail limit of operation for components in the two identically operating systems L i fail limit of operation for components in system i i j L load of component j i in system i L i max maximum value of the initial load sampling interval in system i L i min minimum value of the initial load sampling interval in system i M number of interdependency links N number of nodes in the two identical systems

Fig. 1 .

 1 Fig. 1. The topological model of the IEEE Reliability Test System -96 [41]. Each system i has N i =24 nodes (circles), and K i =34 links (solid lines), i=1, 2. For visual clarity, only six interdependency links connecting the two systems are shown (dashed lines).

2 jDD

 2 = D = 0.02 in this study). If the sum of the initial load of component j i in system i = 1, 2 is larger than the component load threshold L i fail , component j i fails. This failure occurrence leads to the redistribution of additional loads, P i on the neighboring nodes, and I on the interdependent nodes, which may, in turn, get overloaded, and thus fail in a cascade which follows the connection and interdependency patterns of the network systems. As the components become progressively more loaded, the cascade proceeds.

Fig. 2 .

 2 Fig. 2. The average cascade size vs. the average initial load. Triangles: isolated single system. Squares and circles: identical interdependent systems.

Fig. 5 .

 5 Fig. 5. Critical load, L 1 cr = L 2 cr = L cr , in systems 1 and 2 vs. M. Both systems are working at the same varying average initial loads, L 1 = L 2 =L. In these simulations I=0.07.

Fig. 6 .

 6 Fig. 6. Average cascade size, S 1 = S 2 = S, in systems 1 and 2 vs. M, while both systems operate at the same constant working load, L 1 = L 2 =L = 0.75. In these simulations, I=0.07.

Fig. 7 .

 7 Fig. 7. Critical load, L 1 cr = L 2 cr = L cr , in systems 1 and 2 vs. I, while both systems work at the same varying average initial load, L 1 = L 2 =L. In these simulations, M=34.

Fig. 8 .

 8 Fig. 8. Average cascade size, S 1 = S 2 = S, in systems 1 and 2 vs. I, while both systems operate at the same constant working load, L 1 = L 2 = L = 0.75. In these simulations M=34.

ACKNOWLEDGMENT

The authors wish to thank the anonymous referees for the numerous constructive comments and suggestions which have led to a significantly improved work.