Piero Baraldi
email: piero.baraldi@polimi.it

Enrico Zio

Francesca Mangili

Giulio Gola

Bent H Nystad

An hybrid ensemble based approach for process parameter estimation in offshore oil platforms

The valve flow coefficient is commonly used as a parameter to assess the erosion state of choke valves in offshore oil platforms. In particular, the difference between the theoretical value of the valve flow coefficient and its actual value calculated during operation is retained as the valve health indicator. The actual valve flow coefficient is analytically calculated from the oil, water and gas mass flow rates. These quantities, which are allocated on a daily basis based on the measured total production from a number of wells, on physical parameters (pressures and temperatures) related to the specific well, and on a physical model of the process, can be affected by large uncertainties. Based on such values, the evaluation of the health indicator becomes unreliable and undermines the possibility of using it for prognostic purposes. Similar situations arise every time health monitoring rely on unreliable measurement taken by sensors subject to hard working condition, as often happen in the nuclear industry. This paper proposes a method to obtain more accurate daily estimates of the actual values of the oil, water and gas flow rates, from which improved estimates of the flow coefficient will follow. In this respect, an hybrid ensemble aggregating the physical model with data-driven models built using the Kernel Regression (KR) method has been used. Ensemble diversity is ensured by using different training sets;a local procedure based on the historical performance of the models is adopted to aggregate their predictions. The method is verified on real measurements performed on a number of similar offshore choke valves.

Introduction

In this paper, we consider the degradation of choke valves located topside at wells on the Norwegian Continental Shelf [START_REF] Nystad | Technical Condition Assessment and Remaining Useful Life Estimation of Choke Valves subject to Erosion[END_REF] [START_REF] Haugen | Sand Erosion of Wear Resistant Materials: Erosion in Choke Valves[END_REF]. The difference between the actual valve flow coefficient and its theoretical value is retained as the indicator of the choke valve health state and is used to assess the degree of erosion affecting the choke. While the theoretical value of the valve flow coefficient depends only on the choke opening, the actual valve flow coefficient is analytically calculated on a daily basis as a function of the pressure drop through the choke which is directly measured and oil, gas and water flow rates which are allocated based on the measured total production from a number of wells and on physical parameters (pressures and temperatures) related to the single well. Such flow rates are actually measured only during a number of well tests carried out throughout the valve life.

In practice, the resulting indicator of the choke valve state is very noisy and lacks the physical monotonicity of the erosion process. In [START_REF] Baraldi | Ensemble of kernel regression models for assessing the health state of choke valves in offshore oil platforms[END_REF] it has been shown that the allocated values of oil, gas and water flow rates can cause large inaccuracies and uncertainties in the calculation of the actual valve flow coefficient. In this work, an ensemble of Kernel Regression (KR) models has been devised to correct these values based on the relations among all parameters and aggregated to the outcome of the physical model to increase its accuracy. KR is a distance-based regression algorithm [START_REF] Nadaraya | On Estimating Regression[END_REF] [START_REF] Atkeson | Locally Weighted Learning[END_REF]; an ensemble of multiple KR models is used to avoid the need of selecting the optimal model and to increase the robustness and reduce the uncertainty of the estimate [START_REF] Perrone | When networks disagree: ensemble methods for hybrid neural networks[END_REF]. Diversity is injected in the ensemble by differentiating the training procedure for each KR model. The aggregation of the KR model outcomes is obtained based on the models performance on historical data closed to the test data under reconstruction. This approach can be found in literature under the name of local fusion (Barald, 2010).

The paper is framed as follows. The traditional procedure for the construction of a health indicator assessing the choke valve erosion state is presented in Section 2; in Section 3, an ensemble of KR models is proposed to improve the accuracy of the allocated flow rates; Section 4 shows the results of the application of the method to the choke dataset; finally, conclusion and potential perspectives for future work are drawn in Section 5.

Choke valve erosion assessment

In oil and gas industries, choke valves are normally located on top of each well and are used to balance the pressure on several wells into a common manifold to control flow rates and protect the equipment from unusual pressure fluctuations.

In Figure 1, a choke valve is sketched. The throttle mechanism consists of two circular disks, each with a pair of circular openings to create variable flow areas. One of the disks is fixed in the valve body, whereas the other is rotated either by manual operation or by actuator, to vary or close the opening. For large pressure drops, the well streams which contain gas, liquid and sand particles can reach 400-500 m/s and produce heavy metal loss mainly due to solids, liquid droplets, cavitation and combined mechanisms of erosion-corrosion, resulting in choke lifetimes of less than a year. Erosion management is vital to avoid failures that may result in loss of containment, production being held back, and increased maintenance costs. Moreover, several chokes are located subsea, where the replacement cost is high. Then, the need has increased for reliable models to estimate erosion and lifetime of choke valves, in order to allow implementing effective maintenance strategies [START_REF] Wold | New Generation Software Integrates Intrusive and Non-intrusive Systems for Corrosion and Sand/erosion Monitoring[END_REF] [9] [START_REF] Birchenough | Assessment of Choke Valve Erosion in a High-Pressure, High-Temperature Gas Condensate Well Using TLA[END_REF].

Figure 1. Typical choke valve of rotating disk type (http://www.vonkchokes.nl/).

A common indicator of the valve flow capacity is the flow coefficient C V , which is related to the effective flow cross-section of the valve.

For a specific valve opening, erosion produces a gradual increase of the valve area available for the flow transit, thus determining an increase of C V (eq. 1). For this reason, knowing the value of the flow coefficient is fundamental for assessing the health state of the choke. During operation, C V is not directly measured but computed for a two-phase flow as [START_REF]Flow Control Manual[END_REF]:

              o o w w g g p V f f J f P F m C    2 27.3  (1)
where

g w o m m m m       
is the total mass flow rate of the oil-water-gas mixture,

m m f g w o g w o   / , , , ,



is the fraction of the oil, water and gas mass flow rates, respectively, ρ o,w,g are the corresponding densities, J is the gas expansion factor, F p (θ) is the piping geometry factor accounting for the geometry of the valve/pipe reducer assembly and ΔP is the pressure drop through the choke. Eq. ( 1) and the values of ρ o,w,g , J, F p (θ) and N 6 are derived from fluid dynamics; parameters ΔP, θ, o m  , w m  and g m  are measured or allocated during operation.

Choke valve dataset

For a correct assessment of the choke erosion state, it is fundamental to obtain frequent and reliable measurements or estimates of the parameters ΔP, θ, o m  , w m  and g m  used to compute the flow coefficient C V . Nevertheless, only the pressure drop ΔP and the valve opening θ are measured during standard daily inspections (SI), whereas measures of water, oil and gas flows rates are taken downstream of the choke only during well tests (WT) with a multiphase flow separator. On a daily basis, the values of o m  , w m  and g m  are allocated for a single well by a software based on the measured total production from a number of wells and on physical parameters (pressures and temperatures) related to the specific well. In this work, we consider data regarding valve degradation collected on a daily basis from five different wells. Table 1 and

Improving the quality of the allocated parameters

Since the allocated values of o m  , w m  , and g m  are noisy and unreliable [START_REF] Baraldi | Ensemble of kernel regression models for assessing the health state of choke valves in offshore oil platforms[END_REF], an on-line procedure for improving the accuracy of the estimates of those parameters values is here proposed. The procedure is based on empirical models which learn from a training set the relationships between the parameters, and provides as output an estimate x ˆ of the input parameters x. Different regression techniques such as those based on the use of principal component analysis [START_REF] Baraldi | Aggregation of randomized model ensemble outcomes for reconstructing nuclear signals from faulty sensors[END_REF], artificial neural networks [START_REF] Fantoni | Multiple-Failure Signal Validation in Nuclear Power Plants using Artificial Neural Networks[END_REF] [14], support vector machines [START_REF] Sun | Lidar Signal Denoising Using Least-Squares Support Vector Machine[END_REF], evolving clustering methods [START_REF] Chevalier | Assessment of Statistical and Classification Models For Monitoring EDF's Assets[END_REF] have been applied to this purpose. In this work, Kernel Regression models [START_REF] Nadaraya | On Estimating Regression[END_REF] [5] have been chosen.

Kernel regression models

Nonparametric Kernel Regression (KR) is used to build a model for improving the quality of the allocated values of oil, water and gas mass flow rates. Compared with parametric methods, which are defined by sets of parameters and predefined functional relationships, nonparametric methods have the advantage that they do not require any assumption about the mathematical structure of the regression model [START_REF] Nadaraya | On Estimating Regression[END_REF].

KR
      trn trn N k k N k RG k k RG tst w w 1 1 ˆx x (2)
The kernel function K must be such that training patterns with small distances from the test pattern are assigned large weights and vice versa. Among the several functions which satisfy this criterion, the Gaussian kernel is commonly used [START_REF] Wand | Gaussian-based kernels for curve estimation and window width selection[END_REF]:

          2 2 2 2 exp 2 1 ) ( h d h d K PG PG  (3)
where the parameter h defines the kernel bandwidth and is used to control how close training patterns must be to the test pattern to be assigned a large weight. In order to compute d PG , the PG parameters are normalized to mean equal to 0 and standard deviation equal to 1.

In the present case study, the choice of training dataset and predictor parameters is critical. In this respect, four different models can be devised by differentiating the training set as listed in Table 2.

Table 3. Model training procedures.

Model Training set

Predictor parameters

1 Well test data X WT Measured ] , [  P PG k   x 2 Standard inspections data X SI Measured ] , [  P PG k   x 3 Well test data X WT Measured & allocated ] , , , , [ g w o PG k m m m P       x 4 Standard inspections data X SI Measured & allocated ] , , , , [ g w o PG k m m m P       x
In any case, for estimating the test pattern x tst collected for a specific choke, only patterns concerning the same choke are considered. This is done because, although valves are similar, the well behavior, and thus the relation among the observed parameters, tends to vary with time and from one well to another. Thus, data collected from other chokes do not provide useful information about the behavior of the choke under study.

In practice, when estimating the test pattern x k for a specific choke c, only the patterns x j , j=1,…,k-1 previously collected during the life of the c-th choke are used as training patterns.

In order to verify the performance of the process parameter estimation models, we have considered a test set formed by the N WT patterns of X SI collected the same day of a well test. Since an accurate measure of the process parameters under estimation is available for these patterns, they can be used to assess the performance of the estimation models.

The response group is formed by the unreliable parameters that need to be estimated ] , , [

g w o RG k m m m     x .

Ensemble approach

Since the performance of the models depends on the characteristics of the parameter to be estimated and the intensity of the noise, it is difficult to identify a single best model [START_REF] Baraldi | Ensemble of kernel regression models for assessing the health state of choke valves in offshore oil platforms[END_REF].

Using an ensemble of models allows overcoming this dilemma. Indeed, the general idea underlying ensembles is to create many models and combine their outputs in order to achieve a performance better than that provided by each individual model in the ensemble [START_REF] Perrone | When networks disagree: ensemble methods for hybrid neural networks[END_REF]. Models' prediction diversity plays a fundamental role when ensemble approaches are devised. In fact, individual models committing diverse errors can be opportunely combined in such a way that the error of the aggregated prediction is smaller than the error of any of the individual models.

Different techniques for the aggregation of the outcomes of individual models have been proposed in literature, the most common being statistics methods like the simple mean, the median and the trimmed mean [12] [18]. Other aggregation techniques, which allow improving the ensemble performance, consider weighted averages of the model outcomes with weights proportional to the performance of the individual models. In this respect, both global approaches (in which the performance is computed on all the available patterns) and local approaches (which measure the performance only on the patterns closed to the test pattern) have been proposed [START_REF] Baraldi | Local Fusion of an Ensemble of Models for the Reconstruction of Faulty Signals[END_REF]. Since in the choke valve case study a complete input-output set of patterns is not available, model weighting cannot be based on a measure of the performance of the individual model. For this reason, a new strategy is here proposed based on the use of the Analytic Hierarchy Process (AHP) [START_REF] Saaty | The analytic Hierarchy Process, Planning, Priority Setting, Resource Allocation[END_REF].

Outcome aggregation with Analytic Hierarchy Process

AHP is used to assign performance weights to the models of the ensemble. The procedure allows ranking different models outcomes using relative performance measurements, without resorting to an absolute measurement of the model performance. AHP is a multi-criteria decision method that uses hierarchic structures to represent a decision problem and provides ranking of different choices [START_REF] Saaty | The analytic Hierarchy Process, Planning, Priority Setting, Resource Allocation[END_REF]. It consists of two main steps: 1) structuring a hierarchy; 2) assigning priorities to the elements of each hierarchy level by comparative judgments of the elements based on a pre-defined scale.

In this application, the hierarchy structure sketched in Figure 2 is used. The four models on level 3 are compared with respect to the two criteria Z 1 and Z 2 of the level 2 towards the goal (level 1) of obtaining high model accuracy. The basic tool for assigning priorities to the elements of a level of the hierarchy are matrices of pairwise comparisons based on the criteria defined at the previous level. For the hierarchy of Figure 2, two matrices of comparisons A is defined, the vector of priorities l Z π of the models in level 3 of the hierarchy with respect to criterion Z l is given by the eigenvector associated to the maximum eigenvalue of matrix l Z A . The priority vectors obtained for each criterion are weighted with the priority assigned to the corresponding criterion and averaged to obtain the overall priority vector π=[π 1 , π 2 , π 3, π 4 ] assigning the priority π m to model m.

In the proposed aggregation method, the priorities assigned to each model are used as weights to aggregate the models' outcomes through a weighted average: In this application, the first criterion Z

∏ ≠i,j m tst m ij ij a a ) (x  (6)
According to the AHP method, the quality of a matrix of comparison can be evaluated considering its consistency. Matrix

1 Z
A is consistent if the following equation is satisfied for any i, j and k [START_REF] Saaty | The analytic Hierarchy Process, Planning, Priority Setting, Resource Allocation[END_REF]:

ik k i k j j i jk ij a π π π π π π a a    (7)
In our case, substituting eqs. ( 5) and ( 6) in eq. ( 7) gives:

∏ ∏ ∏ ∏ ∏ ∏ ≠ ≠ ≠ ≠ ≠ ≠ k i h ik im km k j i m im km ji jk i k j m jm km ji ki k j i m im jm ik jk k j m jm km j i m im jm jk ij a d d d d d d d d d d d d d d d d d d a a , , , , , , , , ,  
                   where ) , ( , , j RG tst i RG tst ij d d x x =
and, by definition,

ji ij d d =
. This shows that, in the proposed approach, matrix

1 Z A is consistent.
A second criterion Z 2 for evaluating the performance of a model takes into account the RMSE in reconstructing the reliable parameters ΔP and θ, i.e. the root mean square difference between the reconstructed and measured values. This second criterion takes into account the fact that robust and reliable models should be able to correctly reconstruct the reliable parameters despite the noise on the mass flow rates. Since all model performances are evaluated with respect to the same reference, i.e., the reliable measurements of ΔP and θ, the pair-wise comparison is not needed, and the vector of priorities Finally, the two criterions Z 1 and Z 2 of level 2 of the hierarchy are given the same importance and thus the priority vector π is given by:

       2 1 ] 5 . 0 5 . 0 [ Z Z π π π (8)

Results

The KR models in Table 2 are applied to the choke valve case study to obtain estimates of the mass flow rates o m  , w m  and g m  . The performance of the models is evaluated by considering the mean square error (MSE) between the estimates of the mass flow rates o m  , w m  and g m  and the corresponding well test measurements. Table 4 compares the MSE of the four models and that which is obtained by considering the allocations as mass flow rate estimates. Notice that only model 2 allows achieving estimates more accurate than the allocations. Furthermore, using all five parameters as predictor parameters produces, in general, best results. 4, a second ensemble considering only the three best performing models (i.e., models 1, 2 and 4) and a third one considering both the allocations and the three best performing models. Results show that all the three AHP-based ensembles outperform the allocations; moreover, ensembles 2 and 3 outperform all single models.

In Figure 3, the MSE of the allocations and of Ensemble 3 are compared for every choke valve.

The obtained results show that the ensemble estimates do not always outperform the allocation. This can be due to two main reasons: 1.

for some choke valve the allocations are highly accurate, and thus they cannot be improved through the KR estimates; 2.

in some cases characterized by abrupt changes in the operating conditions of the choke (e.g., large variations of the choke opening), the pattern used for the model training may not cover the range of parameter values of interest for the test pattern.

To this purpose, future research work should be devoted to the a priori identification of those cases in which the KR ensemble will not outperform the allocations. In particular, with respect to the situation at point 1) one can have an a-priori indiction of the allocations and the KR ensemble accuracy by considering the corresponding MSE obtained in the estimates of the previous well test performed on the same choke, whereas with respect to the situation at point 2) an estimate of the confidence of the ensemble outcome should be provided by analyzing the position of the test pattern with respect to the distribution of the training patterns.

Conclusions

In this paper, we have tackled the problem of improving the quality of the estimates of some process parameters used in offshore oil platforms for assessing the health state of choke valves in which undergoes erosion. To this purpose, we have proposed a method which estimates the allocated process parameters based on available measurements of other parameters, which are conjectured to be reliable, and on few reliable measurements of the process parameters collected during a number of well tests performed throughout the valve life.

The method is based on an ensemble of Kernel Regression models trained using different procedures in order to inject diversity into the models ensemble. To aggregate the outcomes of the individual models, an original technique based on the Analytic Hierarchy Process (AHP) method has been used. The results obtained on a number of similar eroding choke valves have confirmed the improved performances of the ensemble with respect to any of the single KR models, allowing significant improvement of the oil, water and gas mass flow rates estimates.

It has been shown, however, that in some cases the allocated value of the uncertain parameters are more accurate than those estimated by the proposed KR models ensemble; two main situations which causes this have been identified and a possible strategy for the improvement of the flow rate estimates, based on the identification of these situations, has been briefly sketched.

The development of a method for the implementation of this strategy will be the object of future research work.

Althoght the methodology has been developed in the oil and gas context, a general application of the proposed approach is envisioned in situations in which unreliable parameters' measurements can be improved by resorting to a set of reliable parameters. The goal of the research is the application of such methodology to obtain reliable component health state indicators and prognostic models for remaining useful life estimation also in the nuclear energy field where diagnostics and prognostics are critical issues.

6.

Figure 2 .

 2 Figure 2. Model weighting hierarchy structure.

  provided by model m of the response group parameters in RG tst x .

2 π

 2 Z is computed by taking for each model h=1,…,4, the inverse of its RMSE, i.e.,

Figure 3 .

 3 Figure 3. Comparison of the allocation and estimation performance for each choke.

Table 1 .

 1 Available information

	2

m  real values measured during well tests. Since degraded valves are replaced, data collected for a single well refer to different chokes.

Table 2 .

 2 Number of SI and WT patterns for each choke.

	Well	N SI	N WT	Nval
	1	1854	87	68
	2	2143	96	59
	3	657	39	20
	4	1859	96	54
	5	1678	71	36

  models provide estimates by developing local models in the neighborhoods of the test patterns they are fed with. Estimates are obtained as weighted averages of the training patterns, with weights decreasing as the distance between the test and the training pattern increases. In this view, training patterns closer to the test pattern are conjectured to be more similar to it, thus giving the most relevant contribution to its estimate.

Let X trn ={x k }, k=1,…, N trn be the training set used for the estimate of the test pattern x tst . To develop the KR model, parameters are divided into a predictor group (PG) and a response group (RG) (with the two groups possibly overlapping). For the estimate of x tst , the KR algorithm assigns to each training pattern x k a weight w k =K[d PG (x tst ,x k )], where K is the kernel function which produces the weight for a given distance d PG (x tst ,x k ), between the training and the test patterns, computed considering only the parameters of the predictor group. The estimate RG tst x ˆ of the RG parameters of the test patterns is obtained as a weighted average of the RG parameters of the training patterns:

  Assuming that the model outcomes of the models left out of the pair-wise comparison are distributed around the correct value, this criterion assigns larger weights to the model (i or j) whose outcome is more similar to that of the models left out.

	and the entry a ij of the comparison matrix values ) ( tst m ij a x m=1,…,4, m≠i, j:	1 A is given by the product of the relative importance Z
										Level 1 Goal			High Accuracy
								Level 2 Performance indexes	Z 1	Z 2
										Level 3 Models	Model 1	Model 2	Model 3	Model 4
													1 chosen to evaluate the relative importance	( tst a x ij	)	of model
	x i with respect to model j in the reconstruction of a test pattern tst two models outcomes i RG tst , x and j RG tst , x to the remaining models outcome is the relative similarity of the m RG tst , x , m≠i, j. The similarity of two patterns i RG tst , x and m RG , x has been estimated by the inverse of their tst Euclidean distance ) , ( , , m RG tst i RG tst d x x ; the relative importance ) ( tst m ij a x of a model i with respect to
	model j when model m is taken as reference is defined by:
	a	m ij	(	x	tst	)		d d	( ( x x	, , i RG tst j RG tst	, , x x	, , m RG tst m RG tst	) )		(5)

Table 4 .

 4 Comparison of the performance of the SI allocations with the KR models estimates.Table5compares the performance of the allocations with respect to those obtained by three different AHP ensembles: a first ensemble aggregating all four KR models in Table

			SI allocations	Model 1	Model 2	Model 3	Model 4
	m 	o	0.1205	0.0925	0.0794	0.1086	0.1267
	m 	w	0.1547	0.1967	0.1445	0.2658	0.2060
	m 	g	0.1993	0.2875	0.2068	0.3087	0.2096
	Average	0.1582	0.1923	0.1436	0.2277	0.1808

Table 5 .

 5 Comparison of the performance of the SI allocations with those of three different ensemble estimates

			Allocations	Ensemble 1	Ensemble 2	Ensemble 3
	m 	o	0.1205	0.0672	0.0646	0.0647
	m 	w	0.1547	0.1666	0.1302	0.1271
	m 	g	0.1993	0.2103	0.1915	0.1862
	Average	0.1582	0.1480	0.1288	0.1260
				.

Acknowledgments

The authors wish to thank Erling Lunde and Morten Løes at Statoil ASA for proving us with the operational choke valve data and the IO Center for Integrated Operations in the Petroleum Industry (www.ntnu.no/iocenter) for partly funding this research project.

7.


INTRODUCTION

We consider network protection (NP) by line switching by the system operator, following an attack or a failure to given components. Line switching is common practice for handling line overloads [START_REF] Granelli | Optimal network reconfiguration for congestion management by deterministic and genetic algorithms[END_REF], avoiding loss and/or cost reductions [START_REF] Schnyder | Security enhancement using an optimal switching power flow[END_REF], improving system security [START_REF] Schnyder | Integrated security control using an optimal power flow and switching concepts[END_REF], and the combination of these [START_REF] Shao | Corrective switching algorithm for relieving overloads and voltage violations[END_REF]. NP aims at hindering failure propagation by cutting off the possible "directions" along which failures can spread within the network. The strategy for such protection stands on the possibility of finding the optimal set of lines to be switched-off in order to limit the consequences of cascade failures; in large networks, this poses a combinatorial optimization problem.

In this paper, we address the NP optimization problem by a modified version of the binary differential evolution (MBDE) algorithm [START_REF] Wang | A Modified Binary Differential Evolution Algorithm, Life System Modeling and Intelligent Computing[END_REF]. Differential Evolution (DE) is an emerging population-based stochastic optimization technique [START_REF] Storn | Differential evolution -A simple and efficient adaptive scheme for global optimization over continuous spaces[END_REF]. The standard DE algorithm is simple and efficient and has been successfully applied in various scientific and engineering fields (Greem & Kim 2001[START_REF] Mahdavi | An improved harmony search algorithm for solving optimization problems[END_REF][START_REF] Coelho | An improved harmony search algorithm for synchronization of discrete-time chaotic systems[END_REF][START_REF] Saka | Optimum design of steel sway frames to BS5950 using harmony search algorithm[END_REF], often with superior performance than alternative optimization algorithms, e.g. particle swarm optimization and GAs [START_REF] Ponsich | Differential Evolution performances for the solution of mixed integer constrained Process Engineering problems[END_REF]. MBDE is a novel binary version of DE developed to tackle binary-coded optimization problems [START_REF] Wang | A Modified Binary Differential Evolution Algorithm, Life System Modeling and Intelligent Computing[END_REF]. The goodness of the line switching solutions proposed during the MBDE search is evaluated by a model of cascading failure process which relies on a topological representation of the network connection pattern (Zio & Sansavini 2009). This modeling en-tails considering the interplay between structural characteristics and dynamical aspects and provides indications on the elements critical for the propagation process [START_REF] Zio | Modeling failure cascades in critical infrastructures with physically-characterized components and interdependencies[END_REF], and on the actions that prevent or mitigate undesired effects [START_REF] Motter | Cascade control and defense in complex Networks[END_REF]). Moreover, these models have shed light on the way complex networks react to faults and attacks, evaluating their consequences when propagation dynamics are taken into account. The system response often results in a dramatic cascade phenomenon due to avalanches of node outages. The rest of the paper is organized as follows. Section 2 describes the basis of the model of cascading failures. The description of the MBDE algorithm is provided in Section 3. Section 4 presents the suggested NP strategy. The results of the NP optimization framework applied to a realistic case study, i.e. the Italian 380kV power transmission grid, are summarized in Section 5. Section 6 provides conclusions and points out directions for future work.

THE MODEL OF CASCADING FAILURES IN NETWORK SYSTEMS

We adopt a topological representation of the power grid as a network of N nodes (substations or buses) and K edges (transmission lines). Two types of buses are distinguished here: N G generators, i.e. sources for power, and N D distribution nodes (loads or primary substations), i.e. at the outer edge of the transmission grid. In order to abstractly account for the flow of the physical quantities in the network, we assume that the power is transmitted along the generatordistributor shortest paths [START_REF] Batagelj | Semirings for social networks analysis[END_REF][START_REF] Motter | Cascade-based attacks on complex networks[END_REF][START_REF] Albert | Structural vulnerability of the North American power grid[END_REF] ABSTRACT: A modified binary differential evolution scheme is devised to identify an optimal protection strategy that minimizes the consequences of cascading failures on an electrical grid. The 380 kV Italian power transmission network is employed as a realistic test case for the suggested protection strategy.

We model the load on a component, e.g. a bus or transmission line, as dependent on the number of shortest paths transiting through it, when flow is sent from each available generation node to each distribution node. In this view, the load, L j , or amount of flow passing through a node, is measured by the node betweenness [START_REF] Sabidussi | The centrality index of graphs[END_REF][START_REF] Nieminen | On the centrality in a graph[END_REF], Freeman 1978[START_REF] Freeman | Centrality in valued graphs: A measure of betweenness based on network flow[END_REF][START_REF] Little | Controlling Cascading Failure: Understanding the Vulnerabilities of Interconnected Infrastructures[END_REF], calculated as the number of generator-distributor shortest paths that pass through that node. Assuming that power is routed through the most direct path, the betweenness of a bus gives an approximation of the amount of power flowing through it [START_REF] Albert | Structural vulnerability of the North American power grid[END_REF]. Likewise, the load, L ij , or amount of flow passing through an edge, is measured by the edge betweenness [START_REF] Girvan | Community structure in social and biological networks[END_REF]. Node/edge betweenness measures the influence of a node/edge with respect to the flow of information in the network, especially in cases where the flow of information follows the available generator-distributor shortest path [START_REF] Albert | Structural vulnerability of the North American power grid[END_REF]).

In order to assess the status of a component, its load is compared with the component capacity, i.e. the maximum load that it can process. The capacity is limited by technological constraints and economic factors. The capacity of node j is dimensioned proportionally to its nominal load L j at which it is designed to operate, C j = (1+ α)L j . Likewise, the capacity of link ij is dimensioned proportionally to its nominal load L ij at which it is designed to operate,

C ij = (1+ α)L ij ,
where α > 0 is the tolerance parameter of the power transmission network. This parameter can be regarded as an operating margin allowing safe operations of the component under possible load increments. When α = 0, the system is working at its limit capacity, its operating margin being null: any further load added to a component would result in its failure and, possibly, in the propagation of a cascading failure involving a large part of the system. When a component fails, it ceases to provide its service and is removed from the system. This leads to a redistribution of the shortest paths in the network and, consequently, to a change in the loads of the components remained active. If upon redistribution of the shortest paths, the load on a component increases beyond its capacity, the component fails and a new redistribution follows. The cascade process of successive failures is a function of the nodes states and also of the edge states, extending a common reference model [START_REF] Batagelj | Semirings for social networks analysis[END_REF][START_REF] Motter | Cascade-based attacks on complex networks[END_REF][START_REF] Albert | Structural vulnerability of the North American power grid[END_REF][START_REF] Newman | Finding and evaluating community structure in networks[END_REF], Latora & Marchiori 2005, Duenas-Osorio & Vemuru 2009[START_REF] Zio | Modeling failure cascades in critical infrastructures with physically-characterized components and interdependencies[END_REF].

In this work, the effects of cascading failures are summarized by two indicators, the connectivity loss, C L , and the cascade size, S. The former quantifies the decrease of the ability of distribution substations to receive power from generators; its calculation re-lies on the topological structure of the network and the available "least-resistance" pathways:

1 1 1 D i N G L i DG N C NN    (1)
where N G and N D represent the numbers of generator and distribution nodes in the unperturbed state of the network, and i G N represents the number of generation units able to supply flow to distribution node i after failures take place. While the failure of lines only contributes to the further propagation of the cascade through load redirection, the failure of buses implies also the inability to provide electrical load to areas of the network. For this reason, we select to measure the disruption caused by the failure cascade by the number of nodes made inactive by the cascade, through S, and by the inability of generator nodes to power distributor nodes, through C L .

MODIFIED BINARY DIFFERENTIAL EVOLUTION ALGORITHM

An improved binary DE [START_REF] Wang | A Modified Binary Differential Evolution Algorithm, Life System Modeling and Intelligent Computing[END_REF]) is used for optimizing the NP strategy. DE has been originally developed as a populationbased global optimization algorithm for real-valued optimization problems [START_REF] Storn | Differential evolution -A simple and efficient adaptive scheme for global optimization over continuous spaces[END_REF]. The original idea of DE lies in the way of perturbing the current population: this is obtained by adding to a chromosome the weighted difference between two others randomly selected from the population. The initial population is randomly generated by sampling values uniformly distributed between the lower and upper bounds of each of the decision variables. Three evolutionary operations, i.e., mutation, crossover and selection are used to update the population.

During the first operation, the mutated individual v i , called noisy vector, is generated according to the following equation

1 2 3 1 () t t t t i r r r v x F x x      (2)
where F, a positive real constant, is the scaling factor, t is the index of generation and 1 2 3

, and

t t t r r r
x x x are three randomly chosen individuals, i.e. target vectors, with indexes r 1 ≠r 2 ≠r 3 ≠i. After mutation, the noisy vector is further modified by the crossover process, in which the noisy and target vectors are mixed according to some rule to create the trial vector u i . The trial vector inherits from the noisy and the target vectors different pieces of chromosome, as regulated by the crossover rate CR. The commonly used binomial crossover is defined as:

  1 1 , if (0,1] or irand( )
, if (0,1] and irand( ) 1, 2,...,

t ji t ji t ji v U CR j P u x U CR j P jK            (3)
where U(0,1] is a uniform continuous random value ϵ[0,1], irand(P) is a uniform discrete random number in the set ϵ{1, 2, …, P}, P is the population size, j is the index of the dimensionality and K is the length of the chromosome. The selection operator is defined as:

  11 1 , if f( ) f( ) , otherwise 1, 2,..., t t t i i i t i t i u u x x x iP           (4)
Referring to a single objective (SO) minimization, if the fitness value of the trial individual is better than the target one, it replaces the target individual; otherwise, the target individual is preserved in the next generation. This criterion ensures that the next generation is better than or at least equal to the previous one.

In the NP optimization problem here of interest, the decision variable vector x is composed of K elements corresponding to the network links; each element contains 1 if the corresponding line is to be switched-off, 0 otherwise. The objective function of the optimization is C L (x).

We define the attack/failure scenario whose consequences have to be minimized via an optimal NP strategy. We consider the attack to/failure of a single component whose following cascading process maximizes the damage to the network.

NETWORK PROTECTION STRATEGY

In practice, we identify the component whose removal yields the largest disruption in the network. NP aims at minimizing the effects of the cascade failure propagation by switching lines after the cascade is triggered. Within the modeling framework here adopted, the cascade unfolds rapidly, through the shortest paths.

The NP strategy consists in taking action only once after the cascade is triggered, by opening several lines to prevent failure propagation. The issue is to decide which set of lines to disconnect in order to protect most from the cascade propagation.

The protection optimization seeks interventions that minimize the objective function C L (x), i.e. the connectivity loss of the network configuration x = {x j | x j ϵ {0,1}, j = 1,…,K}, where x j = 1 if link j is operating, or 0 otherwise. Initially, all components, N nodes and K links, are characterized by initial loads L i and L ij , respectively, and maximum capacities, C i and C ij , respectively. The optimization is carried out by the MBDE algorithm, that uses binary variables to characterize the line switching decisions, as detailed in Section 3. The algorithm for simulating the combined effects of the cascading failures and protection strategy proceeds in successive propagation steps as follows.

0. The initiating event occurs, i.e. the failure of the most critical component, which is removed from the network. Each component that is operating is tested for failure: for i = 1,…, N, if node i is working and L i > C i then node i fails; if link ij is working and L ij > C ij then link ij fails. 1. The operator intervenes by switching off a set of lines, in order to limit the cascade failure effects, expressed by the objective function. The MBDE optimization algorithm returns the most appropriate network configuration that minimizes the objective function at the end of the current step t = 1.

2. The components loads are re-distributed taking into account the new, protected network topology. Each component that is operating is tested for failure: for i = 1,..., N, if node i is working and

L i > C i then node i fails; if link ij is working and L ij > C ij then link ij fails.
3. The stage counter t is incremented by 1 and the algorithm is returned to step 2.

4. The algorithm stops when no more working nodes fail.

RESULTS AND DISCUSSION

The NP optimization framework of Sections 3 and 4 has been applied to the topological network of the 380 kV Italian power transmission network (Figure 1). The 380 kV Italian power transmission network is a branch of a high voltage level transmission, which can be modeled as a network of N = 127 nodes (N G = 30 generators and N D = 97 distributors) and K = 171 edges. Figures 2 and3 show the load distribution on each component when the network operates at nominal conditions. Node 88 and link 104 (the line connecting nodes 75 and 88) handle the largest load in the system, i.e. the largest number of generatordistribution shortest paths under our modeling scheme. The initiating event is identified as described in Section 4. Contrary to what can be supposed, it is found that the largest damages to the network are due to removals of components that are not among the most loaded (Figures 4567). The removals of nodes 86 and 78 have the largest impacts on the network with respect to both C L and S. The failure of the most loaded node, i.e. node 88, leads to cascades involving 28 components, compared to 48 components that fail when node 86 triggers the cascade. Also, C L = 0.77 when the system stabilizes if node 88 is removed, but C L = 0.94 when node 78 fails. The criticality of nodes 86 and 78 with respect to cascading failures is due to their position in the network, in proximity of highly loaded nodes. The effects of cascading failures are investigated also with respect to the loss of lines. The removal of link 107 connecting nodes 78 and 81, results in the largest damage in terms of S (57 nodes fail at the end of the cascade failure propagation) and C L = 0.96, while the failure of the most loaded link, i.e. 104, leads to S = 27 components and C L = 0.71. The average initial load has little influence on the criticality, whereas the neighbors load and the failure of the connection between node 78 and 81, have a large impact (node 78 is found to be a critical node in triggering the cascading failure propagation). We observe that, on average, the loss of links results in a more severe damage than that caused by the loss of nodes, both in terms of C L and S. This is due to the fact that failure of a link causes the redistribution of the components loads to other working components (nodes and links). The disconnection of the transmission line (edge) 107 is identified as the worst-case scenario.

Protection

The cascade propagation is triggered at t = 0 by an attack/failure consisting in the disruption of link 107, and produces the loss of six nodes by the end of the first propagation step, t = 1. We aim at reducing the effects of the continuation of the cascade of failures in the network through the application of the protection strategy. The optimization parameters are summarized in Table 1. The optimal NP strategy requires disabling 27 lines, i.e. 21% of the network connections (Table 2, column 2). This represents a feasible number of lines to be switched-off and it is consistent with the size of other protection strategies reported in the literature.

As an example, [START_REF] Salmeron | Analysis of Electric Grid Security Under Terrorist Threat[END_REF]) suggests a protection strategy that involves the hardening of 50% of lines in a network characterized by 24 nodes and 40 lines.

The optimal set of lines identified in the protection strategy is represented in Figure 8. The lines are spread over the entire network with a major concentration in the northern part. By adopting this protection strategy, the operator intervention maintains a good level of working components with acceptable C L .

Figure 8. The dotted lines in this graphic representation of the network are the optimal set of 27 edges to be switched-off to optimally contain the failure cascade Figure 9 presents the effects of the cascade after the initial failure, without any operator intervention (dotted line) and with the operator intervention of switching off the optimal set of 27 lines, both in terms of S (left graph) and C L (right graph). The effects of the NP intervention can be evinced by comparing the evolutions in Figure 9. The optimal line switching intervention at t = 1 reduces S and C L by 87.7% and by 34.6%, respectively. The failures keep propagating in the network and the cascade ends at t = 3. The results obtained when the system is stabilized are satisfactory: S is reduced from 57 nodes to 10 nodes, while C L is improved by approximately 21%. There are four nodes that fail after the operator intervenes with the protection strategy, i.e. Pieve A substation, situated in the Lombardy region and Rondissone, Trino and Turbigo substations, situated in the Piemonte region (Figure 1). The results obtained at the end of the optimization are summarized in Table 2. -8, 5-6, 14-73, 15-17, 17-18, 24-28, 25-27, 28-30, 32-35, 34-38, 41-64, 41-66, 43-46, 48-53, 51-54, 55-56, 79-81, 89-91, 91-92, 92-93, 96-98, 97-103, 103-104, 110-111, 113-120, 115-118, 121-123 0.608 0.758

CONCLUSIONS

We have considered the optimal protection of power networks by extending a previously developed model for cascading failures propagation and embedding it in a MBDE scheme of optimization. The worstcase initiating failure scenario is determined by looking at the effects of the cascading failures triggered by the removal of single nodes or links. The optimal protection strategy of line switching is identified by MBDE so as to limit the cascade failure consequences.

Future research steps will extend the abstract model of cascading failures beyond network abstraction to incorporate physical properties and time dynamics of power transmission. Within this framework, other interdiction scenarios and multiple protection strategies can be developed. In particular, an iterative NP strategy can be implemented, that incorporates realtime constraints arising from the network and load configuration updates given previous protective actions.

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. The 380 kV Italian power transmission network (TERNA 2002, Rosato et al. 2007)

Figure 4 .Figure 5 .

 45 Figure 4. The cascade size, S, caused by the removal of each node (abscissa) in the power transmission system of Figure 1

Figure 6 .Figure 7 .

 67 Figure6. The cascade size, S, caused by the removal of each link (abscissa) in the power transmission system of Figure1

Figure 9 .

 9 Figure9. The cascade evolution in terms of cascade size, S, and connectivity loss, C L , after the failure of line 107 and with no intervention at t = 1 (solid lines) or optimal operator intervention at t = 1 (dotted lines)

Table 1

 1

	. The parameters of the MBDE algorithm
	MBDE parameters
	Population size P	40
	Dimensionality of solutions K	171
	Crossover rate CR	0.8
	Scaling factor F	0.2
	Minimum fitness error	10 -4
	Maximum number of generation	1500

Table 2 .

 2 Best optimal switching scheme

	Number	Optimal line switching	Connectivity loss, C L
	of lines	scheme	computational	end of
			time t = 1	cascade
		2
	27

ACKNOWLEDGEMENTS

This work was funded by the Foundation pour une Culture de Securite Industrielle of Toulouse, France, under the research contract AO2006-01.


INTRODUCTION

Critical infrastructures (CIs) like the electricity, oil & gas supply, rail, road, air, sea transport, internet networks, are highly interconnected and mutually dependent in complex ways, both physically and through a multitude of information and communication technologies (so called cyber-based systems) used for data acquisition and control. The coupling of CIs leads to the concept of "systems-of-systems", which implies that single CIs cannot be studied in isolation from other CIs; rather, it is necessary to assess the limitations that interacting CIs pose on the operating conditions of the individual infrastructures so as to implement adequate protections for preventing failures in one CI from cascading to other dependent CIs.

The 2001 prolonged power crisis in California demonstrates the importance of coupling in interdependent CIs (Rinaldi et al., 2001). The crisis took place when electric power disruptions at various times curtailed natural gas production (first order effects); the latter generated a shutdown of steam injection in heavy oil production (second order effects).

A common denominator is needed to assess the vulnerability of a system that is exposed to natural and accidental hazards, and threats of malevolent acts. The need is to capture the CI vulnerability sources and issues, given its technical and physical features, and the dependencies and interdependencies on other systems. This requires an evaluation of the exposure to different hazards, including threats of malevolent acts.

An analysis aimed at identifying the causes of damage or disruption of services in CIs needs to embrace an all-hazard approach (Waugh, 2004), [START_REF] Pollet | All-Hazard approach for Assessing Readiness of Critical Infrastructure[END_REF], encompassing a general view on the hazards targeting a given system.

We propose a framework for an All-HAZard ANalysis (A-HAZAN) which relies on tailored tabular procedures to organize the qualitative and quantitative features of the system relevant for revealing and highlighting its vulnerabilities. The A-HAZAN framework is intended as a tool for managers, analysts and stakeholders of CIs to carry out the identification of all the sources of vulnerability in an allhazard perspective. The paper is organized as follows. In Section 2, some methodologies to assess the vulnerabilities of CIs are reviewed. In Section 3, the concepts of vulnerability and susceptibility are outlined from the perspective of CIs. In Section 4, the A-HAZAN tabular procedure and a methodology are presented. In Section 5, the A-HAZAN methodology is applied to a literature case study. Conclusions are drawn in Section 6.

METHODOLOGIES OF VULNERABILITY ASSESSMENT

Screening methodologies for prioritizing scenarios of terrorist threats and identifying vulnerabilities in single and interdependent CIs have been proposed. Apostolakis and Lemon (Apostolakis and An All-Hazard Approach for the Vulnerability Analysis of Critical Infrastructures

E.Zio

Ecole Centrale Paris-Supélec, Paris, France, Politecnico di Milano

R.Piccinelli and G.Sansavini

Politecnico di Milano

ABSTRACT: In this paper, a framework is proposed for the All-HAZard ANalysis (A-HAZAN) of Critical Infrastructures (CIs). Starting from the identification of the task of each component in the infrastructure, we use tabular procedures to organize the information on the susceptibility to attacks, to single and cascading failures. All variables and states are identified that may impact on the component's role as a possible source of vulnerability within the CI and towards interdependent CIs. This is a starting point for a quantitative evaluation of the degree of exposure to intentional acts. A case study of literature is taken as an exemplary demonstration of the procedures of the analysis. [START_REF] Apostolakis | A screening methodology for the identification and ranking of infrastructure vulnerabilities due to terrorism[END_REF] and [START_REF] Patterson | Identification of critical locations across multiple infrastructures for terrorist actions[END_REF] focus on the identification of critical locations in infrastructures; these are seen as geographical points that are exposed to intentional attacks. Critical locations are not limited to individual infrastructures but may affect multiple infrastructures: for example, water and electrical distribution systems may occupy the same service tunnels. In the proposed scheme, the conditional probabilities that the terrorists will successfully exploit a vulnerability need to be evaluated. The procedure for this relies on extensive use of expert judgment and may be challenging in practice. The impacts of attacks are treated without including the probabilities of various levels of damage, and without consideration of any intervention by first responders: for this reason, it can be defined conservative. The vulnerabilities and their ranking according to potential impact are eventually obtained by Multi-Attribute Utility Theory (MAUT) [START_REF]Categorizing Risks for Risk Ranking[END_REF]. [START_REF] Konce | Bulk power risk analysis: Ranking infrastructure elements according to their risk significance[END_REF] have proposed a methodology for ranking components of a bulk power system with respect to its risk significance to the involved stakeholders; the likelihood and the extent of power outages when components fail to perform their designed functions are analyzed; the consequences associated with the failures are determined by considering the type and number of customers affected.

Johansson and Hassel (Johansson and H. Hassel, 2010) have proposed a framework for considering structural and functional properties of interdependent systems and developing a predictive model in a vulnerability analysis context. Piwowar et al. [START_REF] Piwowar | An efficient process to reduce infrastructure vulnerabilities facing malevolence[END_REF] have proposed a systemic analysis which accounts for malevolence, i.e., the willingness to cause damage.

The aim of the present work is to develop a systematic framework of system analysis for identifying the vulnerable elements of CIs, considering natural hazards, random failures and intentional attacks. While the first two types of vulnerability are characterized by stochastic uncertainties and can be analyzed by traditional safety analysis tools, intentional attacks require a new way of analysis.

HAZARDS, THREATS AND VULNERA-BILITY

In the United Nations' view, hazard is defined as "a potentially damaging physical event, phenomenon and/or human activity, which may cause loss of life or injury, property damage, social and economic disruption or environmental degradation. Hazards can be single, sequential or combined in their origin and effects" [START_REF] Turner | A framework for vulnerability analysis in sustainability science[END_REF]. The European Cooperation for Space Standardization (ECSS) defines hazard as "an existing or potential condition of an item that can result in an accident" (White, 1974); the condition is associated with the design, fabrication, operation or environment of the item, and has the potential for accidents. These two definitions encompass the idea of hazard described as a "condition prerequisite to a mishap" (UNISDR, 2004) and as "a source of potential harm" (ECSS, 2004). The concept of hazard is strictly tied to the presence of a potential source of difficulty both natural or manmade.

On the other hand, the concept of threat is defined as "a potential intent to cause harm or damage to the system by adversely changing its state" (ISO guide, 2009). This definition is strictly linked to intentional and malevolent acts [START_REF] Apostolakis | A screening methodology for the identification and ranking of infrastructure vulnerabilities due to terrorism[END_REF], and it seems in contrast with the one by the US Homeland Security, which describes threat as a "natural or manmade occurrence, individual, entity, or action that has or indicates the potential to harm life, information, operations, the environment, and/or property" (OHS, 2009). In the latter, little reference is made to the idea of intention embedded in the former definition of threat. From these definitions, the broader concept of hazard emerges as a general condition of potential source of harm. Therefore, it encompasses the intentionality of threats, e.g., terrorist acts that are distinguished by a malevolent intelligence directed toward maximum social disruption.

In the all-hazard approach, malevolent acts, accidental and natural occurrences are all considered. Yet, they require a different analytical treatment. Random accidents, natural failures and unintentional man-made hazards are typically known and categorized by emergency planners. Their occurrence can be typically addressed within a probabilistic framework (Figure 1). Conversely, terrorism is a hazard that eludes a quantification by probability theory due to the intentional and malevolent planning it implies (Figure 1). The concept of vulnerability follows the degree of impact that an hazard has on the CI. In [START_REF] Konce | Bulk power risk analysis: Ranking infrastructure elements according to their risk significance[END_REF], vulnerability is defined as "the degree to which a system, a subsystem or a system component is likely to experience harm due to exposure to a hazard, either a perturbation or stress" or, equivalently, in [START_REF] Apostolakis | A screening methodology for the identification and ranking of infrastructure vulnerabilities due to terrorism[END_REF] as the "manifestation of inherent states of the system (e.g., physical, technical, organizational, cultural) that can be exploited by an adversary to harm or damage the system". Along the same line of thought, vulnerability is also defined as a physical feature or operational attribute that renders an entity open to exploitation or susceptible to a given hazard (OHS, 2009). The United Nations define vulnerability as the conditions determined by physical, social, economic, and environmental factors or processes, which increase the susceptibility of a community to the impact of hazards (UNISDR, 2004).

Quantitatively, vulnerability focuses on three aspects [START_REF] Konce | Bulk power risk analysis: Ranking infrastructure elements according to their risk significance[END_REF]:

-degree of losses and damages due to the impact of a hazard; -degree of exposure to the hazards, i.e., likelihood of being exposed to hazards of a certain degree and susceptibility of an element at the risk of suffering losses and damages; -degree of resilience, i.e., ability of a system to anticipate, cope with/absorb, resist and recover from the impact of a hazard or disaster.

Practically, vulnerability comes from flaws or weaknesses in the design, implementation, operation and management of an infrastructure that makes it susceptible to destruction or incapacitation when exposed to a hazard or threat, or reduces its capacity to resume stable conditions. For example, the vulnerability of the electric power system might be quantified in terms of changes in network characteristics following attacks on nodes, and the scale (e.g., number of nodes/lines lost) or duration of the associated losses. In a somewhat more intuitive interpretation, vulnerability characterizes a system component or an aspect of a system [START_REF] Jönsson | Identifying critical components in technical infra-structure networks[END_REF]. A component is said to be a vulnerability of a system if its failure causes large negative consequences. In this sense, the component is said to be critical and the term vulnerability describes a property which can be employed for ranking the system's components with respect to their criticality.

Given the above, the goals of vulnerability analysis becomes then [START_REF] Kroger | Vulnerable Systems[END_REF]:

1. identifying the set and sequences of events that cause damages and losses;

2. identifying the relevant set of "initiating events" and evaluate their cascading impact on a subset of elements, or the system as a whole; 3. determining and elaborating on (inter)dependencies (within the system and among systems) and on coupling of different orders, given the set of initiating events and observed outcomes. The ultimate goal is to identify obvious and, most important, hidden vulnerabilities in infrastructure systems to act for managing and reducing them. The achievement of these goals relies on the analysis of the system, of its parts and of their interactions. The analysis must account for the environment where the system operates, and for the objectives the system is designed to achieve. During the development of such basic system understanding, first vulnerabilities may already emerge.

In this paper, an All-HAZard ANalysis (A-HAZAN) framework is proposed to grasp the complementary aspects of random failures, or unintentional or natural occurrences, and malevolent attacks (Figure 2).

THREAT HAZARD All-Hazard Vulnerability Analysis

THE A-HAZAN ANALYSIS

We propose a tabular methodology for the allhazard vulnerability analysis of CIs. It aims at identifying the features, operating conditions and failure modes relevant to CI vulnerability. A tabular procedure is developed to collect and organize the relevant information to the vulnerability characterization of the system's components (Table 1). Starting from their functional role, the components of the system are broadly divided into three main categories: namely, user, transmitter and provider. A user is the target or the recipient of a particular task or service, e.g., a load bus in the power transmission network (or the consumer in a water supply system). A transmitter functions as a spreader of the task or service, e.g., the transmission lines in the power grid (or the pipelines in the water supply system). A provider is the component which originates that particular task or service, e.g., the generating units in the power transmission network (or the waterworks in the water supply system).

For each of the components, the relevant features that impact on its role as source of vulnerability are listed. The size is the physical dimension (if the component is large, it may need large protections). The level of supply, of transmission or of demand are the quantities of power supplied, transmitted or consumed by the component, and the fraction of the service produced, transmitted or consumed by the component with respect to the overall power produced, transmitted or demanded by the whole system. The level of protection refers to the logical and physical barriers deployed to prevent or discourage malevolent acts. The social criticality anticipates the impact on public opinion of the effects of the intentional attack, given that it is successfully accomplished. The most relevant consequences here considered are measured in terms of human losses and geographic extension. Other features are i.e., perma-nent outages and transient outages, the effect of environmental conditions and temperature, or maintenance operations.

Other critical aspects of a system's component are its logical position and geographic location (Patterson and [START_REF] Patterson | Identification of critical locations across multiple infrastructures for terrorist actions[END_REF]. It is important to know the position of the component with respect to the system and to the environment, and the connections and interconnections between the component and other systems.

Other specific features related to the physics of the provided service are accounted for in this qualitative analysis step. Examples are given in Section 5 with reference to a power transmission grid.

Then, vulnerability characterization in the allhazard approach considers the following (Table 1):

 Threats. These are potential events characterized by the act of a malevolent intelligence directed towards maximum social disruption (Section 3).  Random failures. These are typically permanent or transient outages due to components' failures and may be identified by standard risk analysis techniques (e.g., FMECA, Hazop, and others).  INTRA-system failures. These are failures within the system, typically dependent failures.

 External causes. This term considers natural hazards, such as meteorological or seismic phenomena, and unintentional humaninduced hazards, such as the processing or the storage of potentially hazardous materials, or nearby military installations. External causes may be further grouped into local external causes, e.g., a lightening striking a building or an aircraft crash, or diffuse external causes, such as earthquakes, hurricanes, flooding or the leakage of explosive or toxic materials.

A-HAZAN OF THE IEEE RTS-96

The IEEE 1996 reliability test system (RTS-96) (IEEE RTS-96, 1999) is considered (Figure 3). It contains 24 buses and 38 transmission lines. The buses consist of nine load-only buses, eight load/generation buses, three generation-only buses and four transmission buses (no load or generation on the bus). The test system does not refer to any particular geographic location; however, in the following, we suppose and characterize the locations of its components in order to contextualize the allhazard framework. Five different types of generating units are considered: oil/steam, coal/steam, oil/combustion turbine (CT), hydro and nuclear generating units. In Table 2, the A-HAZAN Table of a generating unit, i.e., the nuclear power plant sited at bus (ID)18, is exemplified. It is supposed to be located in a plain, near a river or a lake, possibly in a non-seismic area, far from urbanized areas and highways. In general, it can be supposed that the power plant does not lie under commercial air routes. This is a typical location for European nuclear power plants.

In the second column, the features of the power plant are summarized. Along with the generated power (size), i.e., 400 MW of active power and 137 MVAR of reactive power, and the percentage of the overall generated power, i.e., 12% of the entire power generated by the grid (level of supply), the type of generating unit, i.e., nuclear plant, and the total number of units of that type in the system, i.e., 1, are listed. A nuclear power plant has a high social criticality: the impact of an accident in the plant is high, both on public opinion and on public health. Due to its intrinsic dangerousness, a nuclear power plant is provided with maximum security measures and is classified as completely secure. The identification of environmental conditions is more suitable for a component than for an entire plant. For an entire plant, environmental conditions are replaced by the age of the infrastructure, or the characterizing ambiance: brackish air or degree of humidity. Maintenance accounts for all the joint operations regarding fuel supply and waste disposal, as well as components maintenance, e.g., turbines, thermal exchangers, or steam generators, as well as maintenance of the buildings and of the non-operational part of the plant, e.g., offices, air filtering systems.

In the third column, the specific area is detailed. In a general view, it could be assumed as a flat area.

Following the definitions given in Section 4, in the fifth column, threats are grouped under the label "sabotage" and "terroristic attacks", where sabotage is meant as a deliberate action intended to ''damage, disrupt, or subvert the organization's operations for the personal purposes of the saboteur by creating unfavorable publicity, embarrassment, delays in production, damage to property, the destruction of working relationships, or the harming of employees or customers'' (Crino, 1994). For example, in a workplace setting, sabotage is the conscious withdrawal of efficiency generally directed at causing some change in workplace conditions. On the other hand, terroristic attacks are actions intended to cause a strong psychological effect by means of disruption and death.

Random failures are identified via FMEA/FMECA analysis on the power plant. Considering the function of the plant as a provider element, random failure rates are given in (IEEE RTS-96, 1999).

Intra-system features account for the connection of the component to other elements of the system, for example, lines T-29, T-31-1, T-31-2 outages which would prevent the power to flow from the provider to the users along the involved path.

Inter-system features encompass all interdependencies between the system and other infrastructures. We consider the possible interdependencies between the power plant and other infrastructure systems: water and fuel supply, transport network or the interdependencies between the communication system and the provider. Natural hazards and human activities can be found under the definition of external causes. This is strictly tied with the location of the component and of the system in general. Natural hazards can be roughly divided as follows: seismic phenomena (i.e., earthquakes, landslides, etc.), volcano eruptions, meteorological phenomena (frequent meteorological events, i.e., wind, precipitations, snow pack, temperature, etc. and rare meteorological phenomena such as: tornadoes, storms, etc.) and flooding (i.e., hazards resulting from flooding of river, sea, lakes and semi-enclosed water bodies, groundwater, etc.). In Table 6, considering a frat terrain as the site, specific natural hazards are: flooding, storm and earthquakes.

The identification of potential sources of hazard due to unintentional human activities in the proximate areas are considered: processing or storing of potentially hazardous materials (such as explosive, flammable, corrosive, toxic or radioactive materials), aircraft crashes, railway rolling stocks and road traffic, military installations and future human activities in the planning stage, e.g., lands with potential for commercial development.

The external causes can be specified further in local and diffuse hazards, depending on the affected portion of the system. An example of local natural hazard is a stroke lightening on a transmission line. Diffuse hazards are mudslides due to heavy rains or flooding of rivers or lakes, for power plants and transmission towers.

In order to proceed in the analysis, a small portion of the RTS-96 (IEEE RTS-96, 1999) has been selected (Figure 4). It consists of a user (load bus (ID)101), a provider (generating unit located at bus (ID)1), a transmission line (line T-2) and a transmission line with a transformer (line T-7). The A-HAZAN Tables describing the characteristics of these components are here not reported, due to limitation of space. However, the salient aspects are described in the following.

In the A-HAZAN Table for the load bus, several features are highlighted, i.e., the absolute quantity of power requested (size), i.e., 108 MW of active power and 22 MVAR of reactive power needed by the bus (ID)101, and the percentage of power flow required by the particular bus with respect to the entire transmission network power requirements (level of demand), i.e., 3.8% of the total 2850 MW requested by the network, (IEEE RTS-96, 1999). Depending on the type of load bus, the level of protection is taken into account: a transmission grid load bus is typically located far from densely-populated areas and it is usually locked, surrounded by fences but no other complex barriers. The impact of the component role is assumed moderate, since it is expected that if a load bus cannot receive power, the generation of power is easily modulated, the power excess is eliminated, and the overall infrastructure still provides its service. The presence of special costumers on the load bus should be verified, e.g., hospitals, airports, energy-intensive factories. Special attention is devoted to the Environmental conditions, other than natural hazards, characterizing the component in a network. For example, the aging or the degradation of the components, due to specific conditions of the location site, i.e., the humidity of the air or the brackish ambience may cause corrosion or damage to the buildings. The temperature range is also an important issue to account for. An exceptionally hot summer or an extremely cold winter impose additional strain on the component.

In the Table, a description of the hazards is also given. Some of the items have been discussed in Section 4, and in reference to Table 2. For this particular component, the intra-system hazards are referred to failures that may occur when the connections (T-1, T-2 and T-3) between the component and the network are damaged.

A specific A-HAZAN contains the description of the generating unit sited at bus (ID)1. Additional features are, along with the generated power (size) and the percentage of generated power in the grid (level of supply), the type of generating unit and the number of units of the particular type of provider. The level of protection is also considered: the plant is isolated from the urbanized areas and it is usually guarded with security patrols, video surveillance of the entire power plant and alarms. A provider is assumed to have a high level of criticality, because in general its entire power supply cannot be readily replaced by alternative generation. Environmental conditions should also be taken into account; age and degradation affect the power plant (or its constituents) full functionality.

For this particular component, the intra-system hazards are referred to failures that may occur when the connections (T-1, T-2 and T-3) between the plant and the network are damaged.

Two tables are used to report the features pertaining to components referred to as transmitters because they are not the recipients of the electrical service but perform the propagation of the service. Transmitters include transmission lines and transformers. The highlighted physical parameters are: the direction of the lines, from bus 101 to bus 103, their capacities, their lengths, 55 miles, and the electrical characteristics of the transmitter, i.e., resistance (R = 0.055 Ω), reactance (X = 0.211 Ω) and susceptance (B = 0.057 S). Transmission towers are usually located in isolated sites, e.g., open country and they are not provided with any particular fence or barriers. Nor are they watched by patrol. Transmission lines environmental issues are salt pollution depositing on insulators on overhead lines and on substations, or floods and fires adjacent to electrical equipment, e.g., beneath overhead lines. The pruning of trees sited along transmission lines is also crucial: for example, the Italian 2003 blackout was triggered and accrued by two consecutive flashovers towards a tree of two overhead lines (Sforna and Delfanti, 2006).

Particular features of the transformer that connects bus 103 to bus 124: a zero miles length line, a very small resistance, R=0.002 Ω, a reactance X= 0.084 Ω and no susceptance, B = 0.000 S.

CONCLUSIONS

A practical all-hazard analysis framework has been presented for merging two different perspectives on vulnerability of critical infrastructures. On one hand, it captures the vulnerabilities due to random failures and natural hazard; on the other hand, it includes vulnerabilities due to malevolent acts. In this sense, it extends common approaches of system hazard identification.

A general organization of the relevant information on the system components is offered on the basis of their tasks and of the features that characterize them as potential sources of vulnerability.

For the characterization, inter-and intradependencies are considered. The A-HAZAN framework is intended as a tool for managers, analysts and stakeholders of CIs to carry out the identification of the sources of vulnerability in an allhazard perspective. It can serve as an entry point into the quantitative evaluation of the degree of exposure of CIs to hazards of different nature.

The future step of the analysis will be the development of a decision logic framework for evaluating the susceptibility to the hazards that loom over a CI, i.e., random failures, unintentional acts and natural hazards, but also malevolent acts.

Figure 1 .

 1 Figure 1. All-Hazard Approach overview.

Figure 2 .

 2 Figure2. All-hazard vulnerability analysis including the assessment of random, natural and unintentional occurrences, as well as threats of intentional, malevolent acts.

Figure 3 .

 3 Figure 3. Single area IEEE RTS-96 grid (IEEE RTS-96, 1999).

Figure 4 .

 4 Figure 4. Portion of the RTS-96 (IEEE RTS-96, 1999) power grid described in Section 5 including bus (ID)101, generating unit (ID)1 and transmission lines T-2 and T-7.

Table 1 .

 1 The A-HAZAN Table.

	Task (Funtion) Component	Features	Geographic Location	Hazards
	USER	Load Bus	Size	Area	THREATS
			Level of demand Bus Name	RANDOM FAILURES
			Environmental conditions		INTRA-System	Degree Centralities Communities
			Maintenance			Physical
			Level of protec-tion		INTER-System	Logical
			Social criticality			Cyber
						Geographic
					EXTERNAL	Diffuse
					CAUSES
					Natural hazards	Local
					Human activity

Table 2 .

 2 Qualitative part of the A-HAZAN Table for a nuclear generating unit in theRTS-96 (IEEE RTS-96, 1999).

	Task	Compnent	Features			Geographic		Hazards
	(Funtion)					Location
			ID	MW MVAR
							INTENTIONAL ATTACKS
						Area: 1	Sabotage
	USER	Load Bus	101 Size:	108	22	(zone11)	Terroristic attacks
			Level of de-			Bus Name:	RANDOM FAILURES from
			mand: 3,8%			Abel	PRA/FMEA/FMECA
			Level of pro-					Degree
			tection:				INTRA-System:	Communities
			locked, non-				lines outage : T-
			complex bar-				1, T-2, T-3	Centralities
			riers, fences
			Social critical-
			ity: moderate
			Environmental
			conditions: aging and deg-radation of the				INTER-System	Physical: transport network (rail, road and air)
			plant, humidi-
			ty, brackish
			ambience,
			temperature
			Maintenance					Logical: economics, political
								Cyber: Scada System/ Power
								for switches
								Geographic: co-location
							EXTERNAL	Diffuse: earthquakes, flooding,
							CAUSES	tornadoes, storms
								Local: stroke of lightning, land-
							Natural Hazards Human Activity	slides, snow, military maneu-vers, aircraft crashes