E Zio

G Sansavini
email: giovanni.sansavini@polimi.it.

Keywords: Complex Networks, System of Systems, Critical Infrastructures, Smart Grids, Interdependencies, Vulnerability Analysis, Uncertainty I

This paper looks into the vulnerabilities of the electric power grid and associated communication network, in the face of intermittent power generation and uncertain demand within a complex network framework of analysis of smart grids. The perspective is typical for the system of systems analysis of interdependencies in a critical infrastructure (CI), i.e. the smart grid for electricity distribution. We assess how the integration of the two systems copes with requests to increase power generation due to enhanced power consumption at a load bus. We define adequate measures of vulnerability to identify the most limiting communication time delays. We quantify the probability that a reduction in the functionality of the communication system yields a faulty condition in the electric power grid, and find that a factual indicator to quantify the coupling strength between the two networks is the frequency of load-shedding actions due to excessive communication time delay. We evaluate safety margins with respect to communication specifications, i.e. the data rate of the network, to comply with the safety requirements in the electric power grid. Finally, we find a catastrophic phase transition with respect to this parameter, which affects the safe operation of the CI.

This brings the need of assessing and controlling the influences and limitations which interacting CIs impose on their operating conditions [START_REF] Zimmerman | Social implications of infrastructure network interactions[END_REF].

In particular, a focal subject of interest today is the evolution of the networks for electrical energy supply and their conception/renovation as "smart" grids [START_REF] Heckel | Smart Substation and Feeder Automation for a SMART Distribution Grid[END_REF]- [START_REF] Ilić | Modeling of Future Cyber-Physical Energy Systems for Distributed Sensing and Control[END_REF] with distributed generation, as opposed to the centralized power generation structure of the existing electric power grids. One key aspect characterizing smart grids is the combination of the electric power grid with an information and communication network to "improve the observability and controllability of the distribution grid and thereby convert it from a static infrastructure to be operated as designed to a flexible, dynamic infrastructure operated proactively" [START_REF] Heckel | Smart Substation and Feeder Automation for a SMART Distribution Grid[END_REF].

The concepts and configurations of smart grids vary sensibly with respect to the implementation at the different levels of the electrical infrastructure, i.e. the transmission and distribution systems (the level of the individual costumer, and the related pricing issues, lies beyond the scope of the study presented in this paper) [START_REF] Granger Morgan | The many meanings of "Smart Grid[END_REF].

At the level of the transmission system, the efficiency and stability of power system operation could be improved substantially by phase angle measurements at many key locations and the use of Flexible AC Transmission System control devices (FACTS), combined with distributed and autonomous control strategies and high-speed communication.

The level of the distribution system is the one usually referred to when introducing concepts of smart grids. Unlike the transmission grids that are arranged as meshed networks, the distribution networks form tree-like, radial structures, characterized by a central distribution hub where power enters and is routed to the loads along the parallel branches of the distribution feeders. If a malfunctioning occurs, circuit breakers automatically disconnect the entire feeder. The use of multiple power injection points, and sensors and remote control switches that can isolate and cut off the problem would guarantee the continuous power supply to the other buses of the feeder [START_REF] Ababnah | Optimal Control-Based Strategy for Sensor Deployment[END_REF]. Furthermore, at present stage the disconnection of the entire feeder is the only strategy to balance a load excess when generation is suddenly lost due to an emergency. Conversely, a combination of smart meters and advanced distribution automation would allow controlling individual loads along a distribution feeder so that critical services (such Vulnerability of Smart Grids with Variable Generation and Consumption: a System of Systems Perspective Enrico Zio, Senior Member, IEEE, and Giovanni Sansavini T as police stations, hospitals, emergency services) can remain connected, while loads that provide less critical services can be dropped. Furthermore, the deployment of small-size diffuse generation would relieve stresses on transmission and distribution systems, and allow the disconnection of the distribution feeder from the main power system to run it as an isolated island, serving only a few critical loads in case of an emergency (this is not currently allowed for technical, legal, safety, and regulatory reasons). Irrespective of the system level under consideration, the smart grid is an auto-balancing, self-monitoring power grid that has the ability to sense when a part of its system is overloaded and reroute power to reduce overload and prevent a potential outage situation. It enables real-time communication between users and utility, allowing optimal energy usage, and the introduction of renewable energy sources for minimizing the collectivity environmental footprint [START_REF] Heckel | Smart Substation and Feeder Automation for a SMART Distribution Grid[END_REF].

The use of renewable energy sources poses challenges to the design and the operations of the electric power grid. The availability of these sources is strictly dependent on the climatic conditions, i.e. cloud cover or wind speed and steadiness. Then, the power output of solar and wind electric generators is affected by uncertainty and intermittency [START_REF] Conti | Probabilistic load flow using Monte Carlo techniques for distribution networks with photovoltaic generators[END_REF]- [START_REF] Wang | Effects of Uncertainty in Both Component Reliability and Load Demand on Multistate System Reliability[END_REF]: paradoxically, due to its intrinsic stochastic nature, electric power from renewable sources could be abundantly available when it is not needed, while scarce in other cases of necessity. Moreover, electrical energy from renewable sources is typically produced in remote areas, e.g., where strong and steady winds are available or solar power plants can be conveniently located. This power has to be routed through the transmission network and delivered to distant locations where it is consumed. Yet, the existing power transmission infrastructure is not currently designed to carry intermittent electrical flows with possibly large peak flows, and in some cases it is inadequate to serve new remote power plants. As an example, a recent study [START_REF] Zanellini | Rete e vento -Lo sviluppo della rete elettrica italiana per la connessione e l'integrazione della fonte eolica[END_REF] commissioned by the Italian Association of Energy Producers from Renewable Sources (APER) has estimated that between 2008 and 2009 the wind farms located in Southern Italy have lost 700 GWh (25% of the maximum producible energy) due to the limited currentcarrying capacity of the high-voltage transmission lines and the scarce power demand in the area near the plants, corresponding to an economic loss of 144 millions of Euros.

The deployment of small-size renewable power generation at the electric distribution level is also rapidly changing the nature of this infrastructure. Traditionally a passive recipient of a unidirectional flow from the transmission network, the distribution network should now support injections of active power from its users. In this perspective, both the transmission and the distribution networks require a real-time control approach with an increasing level of complexity due to the broad space distribution of the generators and to the high degree of variability that characterizes electrical power from renewable sources [START_REF]Smart Grid Laboratory -Energylab, Smart Grid. Le reti elettriche di domani[END_REF].

Electric transmission networks already have some instrumentation that allows control centers to automatically monitor power flows and open or close circuit breakers at substations as a result of human-and computer-generated control decisions [START_REF] Granger Morgan | The many meanings of "Smart Grid[END_REF]. This communication level is based on heterogeneous technologies and platforms which are poorly integrated [START_REF]SmartGrids' -European Technology Platform for the Electricity Network of the Future[END_REF]. On the contrary, distribution networks lack instruments for the information exchange and the system control.

Considering the present lack of standards, the communication network will have to evolve with the developing Smart Grid. Communication links will need to use all kinds of resources, varying from hard-wired links to fiber optics, power line communication [START_REF] Galli | For the Grid and Through the Grid: The Role of Power Line Communications in the Smart Grid[END_REF], wireless [START_REF] Lai | Wireless network performance for residential demand-side participation[END_REF], satellites and micro-wave links, in order to serve the large geographical territory covered by the smart grid [START_REF] Sood | Developing a communication infrastructure for the Smart Grid[END_REF]. All of these communication links constitute opportunities but also introduce vulnerabilities, due to functional failures and exposition to malevolent attacks aiming at causing systemwide instabilities and blackouts, especially if they are not protected in dedicated supports but can be accessed over the Internet.

The coupling between the electric power and communication networks in Smart Grids raises safety and security issues [START_REF] Ten | Cybersecurity for Critical Infrastructures: Attack and Defense Modeling[END_REF]. In the smart grids, the continuous delivery of electric power is supported by the communication infrastructure and, vice versa, the exchange of information relies on the energy input of the electric power grid. Faults originating in one of the two networks impact on the other interdependent one. For example, the interruption or even the delay of the information flow communicating a control action upon the electric power grid may result in power unbalance and system instability that could trigger a blackout [START_REF] Sood | Developing a communication infrastructure for the Smart Grid[END_REF], [START_REF] Peng | Output Feedback Control of Discrete-Time Systems in Networked Environments[END_REF]. Vulnerability analysis must then be concerned with the interconnected, mutually dependent systems to properly account for the influences and limitations which the interacting electrical power grids and communication networks 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.

This paper looks into the vulnerabilities of the electric power grid and associated communication network, in the face of intermittent power generation and uncertain demand within a complex network framework of analysis of smart grids. Complex network approaches for the analysis of distributed systems [START_REF] Dobson | A loading-dependent model of probabilistic cascading failure[END_REF]- [START_REF] Yeh | An Improved Method for Multistate Flow Network Reliability With Unreliable Nodes and a Budget Constraint Based on Path Set[END_REF] take a holistic perspective in which the overall system behavior emerges as the result of the complex interactions of its constituents. The two main objectives of the analysis are (i) to identify preliminary vulnerabilities by topology-driven and dynamic analysis, and (ii) to guide and focus further detailed analyses of critical areas based on physical codes (e.g. incorporating line impedances and Kirchhoff's laws to assess the power flows from the generators to the loads of an electric network). Purely-topological, complex network models of the electric power grid [START_REF] Dobson | A loading-dependent model of probabilistic cascading failure[END_REF]- [START_REF] Yeh | An Improved Method for Multistate Flow Network Reliability With Unreliable Nodes and a Budget Constraint Based on Path Set[END_REF], which abstract the physical details of the provided service and characterize the topology of network systems, are here extended to capture essential operating features, i.e. the power is routed through the least resistant paths, the components are physically specialized in "generators" and "distributors", the lines are subjected to limits on their carrying-flow capacities and the power generation and consumption levels vary within a range of uncertain values.

In order to account for the interdependence effects due to the coupling [START_REF] Kallitsis | A Framework for Optimizing Measurement-Based Power Distribution under Communication Network Constraints[END_REF], we adopt a framework in which the communication delays are integrated in the electric CI. The communication network is assumed to be made of dedicated control signal channels. The coupling between the two interdependent networks occurs at the components level. Each node of the electric network is coupled with a router that ensures communication with the other elements of the network. The time required to exchange data from the measurement location to a control center or data concentrator, and the time required ultimately to communicate these data to the control devices, is collectively denoted as communication delay or latency. As stated before, excessive communication delay may result in power unbalance and system instability that could trigger a blackout. The calculation method of the communication delay encompasses several sources, including routing delays which are modeled as a series of M/M/1 queues, which symbolizes the cumulative time delay along the routing path from the measurement site to the control unit [START_REF] Stahlhut | Latency Viewed as a Stochastic Process and its Impact on Wide Area Power System Control Signals[END_REF]. According to the shorthand notation developed for queuing systems, the three-part descriptor M/M/1 denotes a one-server queuing system with exponentially-distributed inter-arrival and service times [START_REF] Kleinrock | Queuing g systems[END_REF].

We assess the vulnerability of the two interdependent networks subject to the specific hazard defined by the request of increase of power generation due to enhanced power consumption of one load. To this aim, we introduce vulnerability measures that identify the most critical communication time delays, and assess the probability that a reduction in the functionality of the communication system will generate a faulty condition in the electric power grid.

The paper is organized as follows: the embraced interpretation of vulnerability is detailed in Section II; the modeling of the smart grid CI with interdependencies is presented in Section III; Section III.A describes the model of power flow in the electrical transmission network and Section III.B describes the model of information flow in the communication network; Section III.C details the simulation procedure devised for analyzing the vulnerabilities of the two coupled CIs. In Section IV, the proposed model is applied to two interdependent networks whose structures are based on the 380 kV Italian power transmission network [START_REF] Terna | Rete Elettrica Nazionale, Dati statistici sull'energia elettrica Italia[END_REF], [START_REF] Rosato | Topological properties of highvoltage electrical transmission networks[END_REF], and Section IV.A discusses the numerical results. Conclusions are drawn in Section V.

II. VULNERABILITY ANALYSIS OF CRITICAL INFRASTRUCTURES

While the concept of risk is fairly mature and consensually agreed, the concept of vulnerability is still evolving and not yet established. In general terms, risk refers to a combination of the probability of occurrence of a specific event leading to loss, damage or injury and its extent. These quantities and their associated uncertainties are regarded as being numerically quantifiable.

The term vulnerability has been introduced as the concept of risk is somehow limited in portraying the hazard-centric perception of disasters. A hazard of low intensity could have severe consequences, while a hazard of high intensity could have negligible consequences: the level of vulnerability is making the difference [START_REF] White | Natural hazards: local, national and global[END_REF]. The concept of vulnerability seen as a global system property focuses on three elements [START_REF] Bouchon | The vulnerability of interdependent critical infrastructures systems: epistemological and conceptual state-of-the-art[END_REF]: (i) degree of loss and damages due to the impact of a hazard; (ii) 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 loss and damages; (iii) degree of resilience, i.e., the ability of a system to anticipate, cope with/absorb, resist and recover from the impact of a hazard or disaster.

Vulnerability is here interpreted as a flaw in the design, implementation or operation of an infrastructure system, or its elements, that renders it susceptible to destruction or incapacitation when exposed to a hazard or threat, or reduces its capacity to resume new stable conditions [START_REF] Kröger | Vulnerable Systems[END_REF], [START_REF] Cramer | Metric Optimization-Based Design of Systems Subject to Hostile Disruptions[END_REF]. The latter can be provided with a likelihood (frequency) while a measure for destruction or incapacitation needs specific elaborations. In this study, the vulnerability of the two interdependent networks is expressed in terms of frequency of load-shedding actions due to excessive communication time delay, with respect to several ranges of power increase request.

The goals of vulnerability analysis, and the associated modeling and simulation efforts, could be: (i) given a system and the end state of interest, identify the set of events and event sequences that can cause damages and loss effects; (ii) identify the relevant set of ''initiating events'' and evaluate their cascading impact on a subset of elements, or the system as a whole; (iii) given a system and the end state of interest, identify the set of events or respective event sequences that would cause this effect; (iv) given the set of initiating events and observed outcomes, determine and elaborate on (inter)dependencies (within the system and among systems) and on coupling effects of different orders [START_REF] Kröger | Vulnerable Systems[END_REF].

The achievement of these goals relies on the analysis of the system, its parts and their interactions within the system; the analysis must account for the environment which the system lives in and operates, and finally for the objectives the system is expected to achieve. The ultimate goal is to identify the vulnerabilities for managing and reducing them.

III. DYNAMIC MODEL OF INTERDEPENDENT NETWORKS AND SIMULATION PROCEDURE

The electric power grid and the communication network are represented by two graphs G 1 (N 1 ,K 1 ) and G 2 (N 2 ,K 2 ) with N 1 , N 2 nodes connected by K 1 , K 2 links, respectively. In our case, it has been assumed (without loss of generality) that N 1 =N 2 =N and K 1 =K 2 =K. The networks are specified by NN adjacency (connection) matrices whose entries are 1 if there is an edge joining node i to node j or 0, otherwise, and by connection length matrices whose entries are the physical length of the line connecting node i to node j if there is an edge joining node i to node j or 0, otherwise.

In the electric power grid, the nodes correspond to buses. They are further specified into N G generators, i.e. the providers of the service, and N D distributors, i.e. the recipient of the service. The edges correspond to electrical connections.

In the communication network, the nodes correspond to routers and the edges to dedicated control signal channels that connect them.

The interdependencies within the CI are due to the presence of routers on each electrical bus. Routers send and receive control commands that try to match the generation profile with the power demand in the electrical network according to the increase of electrical load.

The system, then, does not rely on a centralized control center. Control decisions are locally taken at the router/bus level and communicated to the other interacting buses/routers [START_REF] Cramer | Metric Optimization-Based Design of Systems Subject to Hostile Disruptions[END_REF]. During this process, the total time required to the system to balance generation and demand, i.e. the communication time delay, T, is monitored.

A. The functional model of power flow in the electrical transmission network

The functional model of the electric power grid abstracts the physical details while at the same time capturing its essential operating features. The elements of the network are characterized by loading values normalized in the interval [0, 1] that identify how close to the limit capacity they are operating. For generator and distributor nodes, L i  [0, 1] identifies the fraction of maximum power that generator i is injecting in the network, or alternatively, the fraction of maximum power that distributor i is absorbing from the network. The value L i = 1 identifies the maximum power that a generator i can inject in the system, or alternatively, the maximum power that distributor i has agreed to receive by contractual agreement. For lines, L ij  [0, 1] identifies the fraction of the maximum flow-carrying capacity that the line connecting node i to node j is serving. The value L ij = 1 indicates that the line connecting node i to node j, hereafter also called line ij, is operating at its maximum flow carryingcapacity.

The loading scenario is identified when the N values of L i and the K values of L ij are known. To account for the variability that characterizes generation and consumption, the L i values are sampled from uniform probability distributions in the intervals [LG min , LG max ] for generators and [LD min , LD max ] for distributors.

We assume that the flow of electrical power from generator i to distributor j follows the least resistant path between i and j. We further assume line resistance to be proportional to line length; hence the least resistant paths and the shortest paths coincide. We find the shortest paths that connect any generator i to any distributor j, and rank in descending order the lines with respect to their participation to the shortest paths. The L ij values are then also sampled from uniform probability distributions. In particular, two probability distributions are assumed that describe two possible states [ ] if line ij belongs to the second half of the rank. This sampling procedure assigns on the average a higher power flow to those lines that are traversed by many generator-distributor shortest paths, i.e. that have a high betweenness centrality [START_REF] Tizghadam | Betweenness Centrality and Resistance Distance in Communication Networks[END_REF].

In the following, we assess how the CI copes with an increase in the electrical power demand. To this aim, starting from an initial loading scenario, we simulate a request of additional power, W, by distributor i. We assume that W is described by a uniform probability distribution [W min , W max ]. The overall load L i + W may also possibly exceed the maximum contractual power that i has agreed to receive from the system, i.e. 1.

B. The model of information flow in the communication network

In [START_REF] Stahlhut | Latency Viewed as a Stochastic Process and its Impact on Wide Area Power System Control Signals[END_REF], Stahlhut et al. assess how the communication delay for measurements and control signals in a power system impacts the control system response. In the following, we embrace their calculation of the communication delay, and assume that the data transmitted are in the form of packets. The packets are a formatted block of information and are typically arranged in three sections: the header, the payload, and the trailer. The header has the following information: packet length, origin and destination address, packet type, and packet number (if a sequence of packets is being sent). The payload carries the data taken from the measurement. The trailer is at the end of the packet and carries information which permits the receiving device to identify the end of the packet. The time delay calculation encompasses several different delays that occur in communication systems. These delays typically are:

• serial delays: the delay of having bits being sent one after another;

• "between packet" serial delays: the time after a packet is sent to when the next packet is sent;

• propagation delays: the time required to transmit data over a particular communication medium;

• routing delays: the time required for data to be sent through a router, and resent to another location.

The total signal time delay may be represented as

s b p r T T T T T     (1) s s r T P D  (2) P T l v  (3)
where T s is the serial delay, T b is the between packet delay, T p is the propagation delay, T r is the routing delay, P s is the size of the packet (bits/packet), D r is the data rate of the network, l is the length of the communication medium, and v is the velocity at which the data are sent through the communication medium (e.g., 0.6 to c, where c is the speed of light).

Routers are a major part of the communication infrastructure in a communication network. There are a number of methods by which routing delay can be approximated [START_REF] Irvine | Data Communications and Networks[END_REF]- [START_REF] Lian | Performance evaluation of control networks: Ethernet, ControlNet, and DeviceNet[END_REF]. A simple and fast calculation method of routing delay can be used to approximate the delay at a router: this calculation method is based on a series of M/M/1 queues [START_REF] Kleinrock | Queuing g systems[END_REF], in which a path from the measurement to the control center is traced, and all of the routing delays are added up to represent the total routing delay for the measurement. For this calculation,

  1 N i rr i TT    (4)
where T r is the total routing delay,   i r T is the routing delay for a single router at location i, and N is the total number of routers. The M/M/1 queue is used to approximate the value of

  i r

T

for each router i, and it is modeled as a one-server queuing system with exponentially-distributed inter-arrival and service times [START_REF] Kleinrock | Queuing g systems[END_REF]. The M/M/1 queue has several associated performance measures such as: the average number of customers in line, the average number of customers in the system, the time waiting in line, and the time required to pass through the system [START_REF] Daigle | Queuing Theory with Applications to Packet Telecommunication[END_REF], [START_REF] Gross | Fundamentals of Queuing Theory[END_REF]. In this context, the customers are the sensory messages in the communication network. A noteworthy performance measure is the total waiting time (system time), and this encompasses both the amount of time waiting in line (waiting time), as well as the amount of time being served (service time). The total waiting time is quantified as

       ( 5
)
where  is the rate at which objects come in the system (e.g., packets/s) and µ is the rate at which the objects are being served (e.g., packets/s). The routing time   i r T for each router i can be estimated from [START_REF] Granger Morgan | The many meanings of "Smart Grid[END_REF] and used in (4) to calculate the total routing delay

          1 i N r i i i i T         (6)
We assume that information from router i to router j flows along the shortest path connecting i and j. Therefore, the total routing time, T r , from router i to router j is evaluated as the sum of the routing delays of the encountered routers along the shortest path connecting i and j.

The maximum time delay that the smart grid can tolerate depends on the system task to be executed and ranges from tenths of milliseconds up to tenths of seconds [START_REF]Smart Grid Laboratory -Energylab, Smart Grid. Le reti elettriche di domani[END_REF], [START_REF] Sood | Developing a communication infrastructure for the Smart Grid[END_REF]. For example, the standard specifies that for an unintentional island in which the distributed generation energizes a portion of the distribution network, the control system shall detect the island and cease to energize the portion within 2s of the formation of an island. Moreover, latencies in the tens of ms allow rapid detection of faults and are the accepted fault detection times [START_REF] Sood | Developing a communication infrastructure for the Smart Grid[END_REF]. On the contrary, voltage and active power regulation actions should be accomplished within 10s and some minutes, respectively [START_REF]Smart Grid Laboratory -Energylab, Smart Grid. Le reti elettriche di domani[END_REF].

C. Simulation procedure

To test the vulnerability of the CI with respect to the stability upon a power increase request, W, we proceed in successive stages as follows: 1. All N components of the electric power grid are working under independent uniformly random initial loads. In particular,

L i  [LG min , LG max ] for generator i = 1,…N G , L i  [LD min , LD max ] for distributor i = 1+ N G ,…N D + N G , and L ij  [L high min , L high max ], or L ij  [L low min , L low max ]
, if line ij belongs to the first or second half of the betweenness centrality rank (Section III.A), respectively, i , j = 1,…, N. These values identify the loading scenario. 2. A distributor node i is chosen randomly and its power demand is increased from L i to L i + W, where W is sampled from the uniform probability distribution [W min , W max ].

The N

G generators that are connected to i are sorted by ascending distance from i, i.e. ascending resistance. The closest generator j = 1 is selected. Communication delay T is set to zero. 4. Distributor i queries generator j for additional power W.

Generator j answers providing the value of the power available to i, i.e.

   

min ,1 j SG j W L 

, where 1 -L j is the maximum additional power that j can provide. Communication time delay, T, is updated. 5. Check whether the shortest path between i and j can accommodate the additional power SG(j). To this aim, every line yz on the shortest path between i and j is tested for infringement of the current-carrying capacity limit in the new configuration, i.e. L yx + SG(j)  1. If at least one limit is exceeded the additional power is updated to the maximum power value that satisfies all the line constraints, i.e.

, for every line yz on the shortest path between distributor i and generator j, SG(j) = 0 in case that a line yz is already operating at its maximum flow-carrying capacity. The residual power that distributor i still requires, U(i), and the power that generator j could not supply, UG(j), are evaluated as

    1 () j s U i UG j W SG s     
. They are recorded for the vulnerability evaluation. Communication time delay, T, is updated. 6. The generator index j is incremented by 1, and the procedure is returned to step 4 until U(i) = 0, in case of a successful system state, or j = N G and U(i) > 0, in case the power request is not satisfied. The final communication time delay, T, is recorded. The value of U(i) at the end of the simulation is the fraction of the additional load W that has to be shed in order to balance the power generation and the power consumption. Therefore, we assume that load can be continuously shed.

Furthermore, it can be noticed that the power that generator j could not supply, UG(j) can be greater than 0 for two different causes: either j is operating at its maximum capacity, or there is no path for the power produced by j for reaching the distributor i where it is needed. This feature of the algorithm is consistent with the behavior of smart grids that include intermittent and remotely-generated electricity from renewable energy sources, as discussed in Section I.

To incorporate the effects of load and generation uncertainties, the above algorithm is embedded in a Monte Carlo simulation framework, in which a large number of additional load requests, e.g., 100000 in this study, is simulated for the same range of [LG min During the simulations, relevant statistics are recorded to evaluate different vulnerability indexes, i.e. the distribution of the communication time delay, T, in case of a successful (no load-shedding is required to balance consumption and generation) or unsuccessful (load-shedding is required to balance consumption and generation) final system state.

The effects of the interdependencies between the electric power grid and the communication network are evaluated with respect to the communication time delay, T , that elapses from the instant in which the request for power increase is done, until the request is fulfilled. As seen in Section III.B, the maximum time delay requirements, T MAX , for the smart grid range from tenths of milliseconds up to tenths of seconds depending on the system task to be executed. In this respect, various communication network configurations differing for the value of the data rate, D r , can be tested to quantify the probability that a malfunction in the communication system, i.e. a communication delay, will generate a faulty condition in the electric power grid. The compliance with T MAX , i.e. T < T MAX , is evaluated at the end of the simulations.

IV. CASE STUDY

The approach to evaluate the vulnerability of the smart grid CI with interdependencies between power and communication networks introduced in Section III, is exemplified with reference to the topological network of the 380 kV Italian power transmission network (Fig. 1), focusing only on its structure with no further reference on the electrical properties, coupled with a dedicated communication network. For simplicity, but with no loss of generality, we assume that the two networks have the same topologies. The developed methodology can be applied to interdependent networks with different topologies. 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 generator and N D =97 distributor nodes) connected by K=171 links [START_REF] Terna | Rete Elettrica Nazionale, Dati statistici sull'energia elettrica Italia[END_REF], [START_REF] Rosato | Topological properties of highvoltage electrical transmission networks[END_REF], defined by its N×N adjacency (connection) matrix a ij and connection length matrix l ij . We assume that the communication network consists of dedicated control signal channels, which follow the same layout of the overhead electrical lines. Therefore, both the electric power grid and the communication network are described by the same a ij and l ij matrices (without loss of generality): the topology and the geographical location of the 380 kV Italian power transmission network components serve as reference.

Fig. 1. The 380 kV Italian power transmission network [START_REF] Terna | Rete Elettrica Nazionale, Dati statistici sull'energia elettrica Italia[END_REF], [START_REF] Rosato | Topological properties of highvoltage electrical transmission networks[END_REF].

We test the vulnerability of the two interdependent networks with respect to the stability upon a power increase request, W, for eight increasing ranges of [W min , W max ] reported in Table II, while the other system parameters are constant and specified in Table I. We assume that the electric power grid operates with an average 20% safety margin for all the components except for the highly loaded lines which operate with an average 15% safety margin. The data for the communication network have been taken from [START_REF] Stahlhut | Latency Viewed as a Stochastic Process and its Impact on Wide Area Power System Control Signals[END_REF] and assume a dedicated communication network. Then, we test the vulnerability of the two interdependent networks with respect to the stability upon a power increase request, W  [W min , W max ], for thirteen values of data rate of the network, D r , in Table III which identify different dedicated network configurations for the communication CI. We aim at quantifying the probability that a malfunction in the communication system, i.e. a communication delay, will generate a faulty condition in the electric power grid.

A. Numerical results

Table II reports the frequency of unsatisfied power request, i.e. load-shedding intervention, with respect to the eight sampling intervals [W min , W max ]. As W max increases, a greater amount of additional load W is requested on the average and the generators closest to the increasing load cannot provide the total required amount. Thus, distant generators are queried, but power cannot always find a path to reach the increasing load due to growing likelihood of line limit infringement along the path, and the system has little chances to fully satisfy W. [Fig. 2-7 in the previous manuscript were removed in order to comply with the 12 page limit for IEEE Transactions on Systems, Man and Cybernetics, Part A -Systems and Humans regular papers. Table 3-6 in the previous manuscript were removed in order to comply with the 12 page limit for IEEE Transactions on Systems, Man and Cybernetics, Part A -Systems and Humans regular papers. The comments to the removed results were also removed.] The frequency of unsatisfied power request in Table II characterizes the self-limitations that the electric power grid imposes on its own operations due to its structural configuration, the limits on flow-carrying capacities of lines, the intermittent nature of power generation from renewable energy sources and the uncertainties in the load demand profile. We now devote our attention to the temporal aspects that characterize the power balance procedure in the CI. In this respect, Fig. 2 shows the probability distribution function of the communication time delay, PDF(T), when the request to increase load is not satisfied by generators, i.e. U > 0, for the eight sampling intervals [W min , W max ] in Table II. Fig. 2 is relative to situations in which the additional load request W could not be met by generators, and the communication signal spans the entire extension of the network to query for additional generation capability. This appears from the larger T values in the PDF(T) plot, which are the communication time delays needed to route the power to the increasing load bus from remote regions of the network. II.

The differences in the communication time delay, T, distribution for normal operations and load shedding appear in Fig. 3, where the CDF(T)'s are plotted for these two situations and encompass the eight sampling intervals [W min , W max ] in Table II. We expect that CDF(T) for normal operations (solid curve in Fig. 3) is shifted to smaller T values than the one relative to load-shedding interventions (dashed curve in Fig. 3). When the extra power request is satisfied, U = 0, only a fraction of the total number of generators is queried for additional power. On the contrary, when the extra power request is not satisfied, U > 0, the communication signal spans the entire generation set looking for additional generation capability that is not found. Given the maximum admissible communication time delay, T MAX , Fig. 3 allows identifying the frequency of load-shedding due to the excess of communication latency, i.e. T > T MAX , even for those loading scenarios that would eventually result in U = 0 otherwise (solid line in Fig. 3). As detailed in Section III.B, T MAX depends on the system task to be executed and ranges from tenths of milliseconds up to tenths of seconds. In this context, T MAX can be thought of as the maximum time delay before instabilities arise in the electric power grid due to the imbalance between generation and consumption. If the value T MAX equal to 100 ms is assumed, we evaluate CDF(T = T MAX = 100 ms) = 96.77% ± 0.06%, and we discover that due to latency constraints the frequency of faulty conditions is 3.23% ± 0.06% even for those scenarios that we classified as normal operations because the additional power request is fully satisfied, i.e. U = 0 (solid line in Fig. 3). If a narrower latency constraint is assumed, e.g., T MAX = 50 ms, we evaluate CDF(T = T MAX = 50 ms) = 87.85% ± 0.10%, and conclude that the frequency of faulty conditions is 12.2% ± 0.10% even for those scenarios that we classified as normal operations because the additional power request is fully satisfied, i.e. U = 0.

The frequency of faulty conditions due to the excess of communication time delay, T, quantifies the limitations that the communication system imposes on the operations of the interdependent electric system. It can be used to measure the strength of the interdependency or the degree of the coupling between these two systems within the CI.

After having found an index that quantifies the extent to which the delay in the communication network impacts on the operations of the electric power grid, we investigate how the frequency of faulty conditions due to the excess of communication time delay, T, varies with respect to different configurations of the communication network. We test several networks characterized by decreasing values of the data rate, D r (Table III first column), that is defined as the number of bits that are conveyed or processed per unit of time. Since higher investments correspond to higher D r , this kind of analysis allows identifying the optimal communication network that respects the constraint relative to the frequency of faulty conditions due to the excess of communication time delay, T. Alternatively, we can evaluate the extent to which a malfunction that reduces the speed of communication, D r , will impact on the power grid due to the interdependencies between the two systems.

Fig. 4 shows the cumulative distribution functions of the communication time delay, T, for decreasing values of the data rate of the communication network, D r , reported in Table III first column, and W  [0.10, 0.25]. These results correspond to scenarios of the electric power grid for which a final balance between generation and consumption could be found, i.e. U = 0. Fig. 4 allows evaluating the extent to which D r values influence the degree of coupling between the two CIs. Smaller D r values produce CDF(T) that are shifted to larger T values. If we assume a maximum admissible communication time delay, T MAX , from Fig. 4 we identify the frequency of load-shedding due to the excess of communication latency, i.e. T > T MAX , even for those loading scenarios that would eventually result in U equal to 0 otherwise. As an example, if T MAX is equal to 100 ms, the scenarios that are found to the right of the vertical dash-dotted line in Fig. 4 will require load-shedding activity. We quantify the effects of the interdependency between the two networks by two complementary indexes that provide safety margins with respect to the coupled operations of the two CIs. The 95 th percentile of the distribution of T quantifies the minimum admissible time delay, T MAX , that limits to 5% the frequency of load-shedding due to excessive latency. From Table III second column and Fig. 5, we see how the 95 th percentile varies with D r . For small D r , e.g., D r = 2‧ 10 4 bps, the system must tolerate a minimum latency T MAX greater than or equal to 2.069 s ± 0.140 s, in order to limit to 5% the frequency of load-shedding due to excessive latency. In other words, power increase operations that must be safely carried out within a time smaller than 2.069 s ± 0.140 s should not be attempted, if one wants to comply with the safety margin of 5% load-shedding frequency. Conversely, if D r is equal to 5‧ 10 6 bps, the system can admit a minimum latency as small as T MAX equal to or greater than 0.6678 s ± 0.0149 s, and the frequency of load-shedding due to excessive latency would still be limited to 5%. Because the admissible latency T MAX is a constraint in the system, a somewhat more intuitive safety margin is provided with respect to P(T  T MAX ) = CDF(T = T MAX ), i.e. the probability of final safe conditions (no load-shedding is required), given the admissible latency T MAX . In Table III III third column, we see that when D r = 1‧ 10 5 bps, the probability of final safe conditions is 31.18% ± 0.15% and, consequently, the frequency of load-shedding due to excessive latency is 68.82% ± 0.15%. The information in Fig. 6 can be used during the design phase of the CI with the aim of identifying the communication specifications that comply with the required safety standards, given T MAX , i.e. the smallest admissible time delay, and [W min , W max ], i.e. the electric load profile. In our example, a communication network with D r  1.26‧ 10 6 bps has to be selected to comply with the 95% safety margin, i.e. 5% probability of load-shedding, when T MAX = 100ms and W  [0.10, 0.25]. Furthermore, the information in Fig. 6 can be used during operations to quantify the extent to which a possible loss of quality in the communication, i.e. a D r reduction, impacts on the electric power grid. In our example, a D r reduction from D r = 5‧ 10 6 bps (normal operations) to D r = 1.94‧ 10 5 bps (anomalous operations), entails a reduction of the safety margin from 96.77% ± 0.06% to approximately 70%, i.e. the probability of load-shedding increases from 3.23% ± 0.06% to approximately 30%. Yet, most important is the fact that an additional small D r reduction to D r = 7‧ 10 4 bps entails the complete disappearance of the safety margin. Therefore, the knowledge of the relation in Fig. 6 is vital for operating the CI with interdependencies in a region of the system parameters that is far from the catastrophic phase transition point.

V. CONCLUSIONS

In this paper, we have looked into the vulnerabilities that arise when complex networks are coupled, with reference to smart grids modeled as a combination of the electric power grid and the communication network. We have adopted a modeling framework in which the communication delays are integrated in the electric CI. The complex network approach proposed extends the purely-topological models of the electric power grid to capture essential operating features, i.e. the power is routed through the least resistant paths, the components are physically specialized in "generators" and "distributors", and the lines are subjected to limits on their carrying-flow capacities; moreover, it includes the further

,

  LG max ], [LD min , LD max ], [W min , W max ], in order to obtain statistically significant results for various realizations of the same average loading condition.

Fig. 2 .

 2 Fig. 2. Probability distribution function of the communication time delay, PDF(T), when the request to increase load is not satisfied by generators, i.e. U > 0, for the eight sampling intervals [W min , W max ] in TableII.

Fig. 3 .

 3 Fig.3. Cumulative distribution function of the communication time delay, T, when the request to increase load is not satisfied, i.e. U > 0 (dashed line), and when the request to increase load is fully satisfied by generators, i.e. U = 0 (solid line), for the eight sampling intervals [W min , W max ] in TableII.

  Fig.3. Cumulative distribution function of the communication time delay, T, when the request to increase load is not satisfied, i.e. U > 0 (dashed line), and when the request to increase load is fully satisfied by generators, i.e. U = 0 (solid line), for the eight sampling intervals [W min , W max ] in TableII.

  The communication time delay, T, is distributed within the intervals T[62.[START_REF]Smart Grid Laboratory -Energylab, Smart Grid. Le reti elettriche di domani[END_REF], 951] ms and T[643, 2000] ms, in case of normal operations and of load-shedding actions, respectively.

Fig. 4 .

 4 Fig. 4. Cumulative distribution functions of the communication time delay, T, for decreasing values of the data rate of the communication network, Dr, reported in Table III and W  [0.10, 0.25]. Smaller Dr values yield distributions shifted towards larger T values. The vertical dash-dotted line indicates the maximum allowable communication time delay, T MAX = 100 ms.

Fig. 5 .

 5 Fig. 5. 95 th percentile of the communication time delay, T, vs. the values of the data rate of the communication network, Dr, and W  [0.10, 0.25].

  third column and Fig. 6, P(T  T MAX = 100ms) = CDF(T = T MAX = 100ms) values are presented as functions of D r , given W  [0.10, 0.25] and T MAX = 100ms. From Table

Fig. 6

 6 displays how the probability of final safe conditions decreases as D r progressively decreases from D r = 1‧ 10 6 bps to D r = 5‧ 10 3 bps. P(T  T MAX = 100ms) exhibits a catastrophic transition in the range D r  [3‧ 10 5 , 7‧ 10 4 ] bps, as D r values decrease below 3‧ 10 5 bps. No final safe condition is realized below D r = 7‧ 10 4 bps, i.e. every additional power request, W, involves some load-shedding actions due to excessive latency in the communication, T > T MAX .

Fig. 6 .

 6 Fig. 6. Probability of final safe conditions, i.e. P(T  T MAX ) = CDF(T = T MAX = 100ms), given the admissible latency T MAX , vs. data rate of the communication network, Dr, and W  [0.10, 0.25].

TABLE II FREQUENCY

 II OF UNSATISFIED POWER REQUEST, I.E. LOAD-SHEDDING INTERVENTION, AND THE ASSOCIATED STANDARD ERROR VS. THE WIDTHS OF THE EIGHT UNIFORM SAMPLING INTERVALS FOR THE ADDITIONAL POWER, W. THE NUMBER OF SIGNIFICANT DIGITS THAT WE REPORT THROUGHOUT THE PAPER IS EVALUATED BY YONEDA'S RULE[START_REF] Song | Omitting Meaningless Digits in Point Estimates: The Probability Guarantee of Leading-Digit Rules[END_REF]. Sampling interval for the additional power request, W

		Frequency of unsatisfied request	Standard error
	[0.10, 0.15]	0.6670 %	0.1252 %
	[0.10, 0.20]	1.1950 %	0.1263 %
	[0.10, 0.25]	2.6000 %	0.1717 %
	[0.10, 0.30]	6.565 %	0.384 %
	[0.10, 0.35]	11.788 %	0.316 %
	[0.10, 0.40]	17.587 %	0.554 %
	[0.10, 0.45]	23.083 %	0.425 %
	[0.10, 0.50]	28.835 %	0.553 %

TABLE III 95

 III TH PERCENTILE OF THE COMMUNICATION TIME DELAY, T, (SECOND COLUMN) AND PROBABILITY THAT T  T MAX (THIRD COLUMN) WITH STANDARD ERRORS, FOR DECREASING VALUES OF THE DATA RATE OF THE COMMUNICATION NETWORK, DR, (FIRST COLUMN) AND W  [0.10, 0.25].

	Data rate of the network, Dr [bps]	95 th Percentile of the communication time delay, T [ms]	P(T  T MAX ) = CDF(T = T MAX = 100ms)
	5‧ 10 6	66.789 ± 1.489	96.77% ± 0.06%
	2.5‧ 10 6	80.94 ± 2.95	96.47% ± 0.06%
	1‧ 10 6	108.986 ± 1.029	94.69% ± 0.07%
	5‧ 10 5	148.97 ± 4.31	88.89% ± 0.10%
	3‧ 10 5	202.22 ± 3.25	81.14% ± 0.13%
	2‧ 10 5	268.91 ± 6.27	72.48% ± 0.14%
	1‧ 10 5	468.9 ± 12.5	31.18% ± 0.15%
	9‧ 10 4	513.33 ± 3.31	30.99% ± 0.15%
	7‧ 10 4	640 ± 53	0%
	5‧ 10 4	869 ± 41	0%
	2‧ 10 4	2069 ± 140	0%
	9‧ 10 3	4513 ± 300	0%
	5‧ 10 3	8069 ± 200	0%

ACKNOWLEDGMENT

The authors are very thankful to the four anonymous reviewers for helping in improving the paper through their observations and suggestions.

complexity that arises from the variability of power generation from renewable energy sources and intermittent load consumption. However, given the somewhat abstract level of the modeling supporting complex network analysis, the insights gained with respect to the vulnerable areas in the system (first findings) may not be clear-cut, and major hidden vulnerabilities may still be expected. Then, to accurately model the physical behavior of the interconnected systems, computational frameworks that propagate the flows of power and information based on physical laws have to be embraced, and more detailed information about the system and its operating environment may be needed. For large, real-scale systems, this requires the development of adequate computational methodologies through clustering computing, co-operative simulation, or similar architectures [START_REF]Advances in Computers -Parallel, Distributed, and Pervasive Computing[END_REF]- [START_REF] Ferreira | A Cooperative Simulation Framework for Traffic and Transportation Engineering[END_REF].

We have quantified the limitations that both the electric power grid and the communication network impose on the operations of the CI with interdependencies, when a load bus demands additional power from the generators. In the application to the 380 kV Italian power transmission network, with the assumptions made and the numerical data used, we have observed that a factual indicator for quantifying the coupling strength between the two integrated systems is the frequency of load-shedding actions due to excessive communication time delay. By means of this indicator, we have evaluated the extent to which a loss of quality in the communication impacts on the electric power grid, and we have selected appropriate communication specifications, i.e. the data rate of the network, that comply with the required safety requirements in the electric power grid.

Finally, we have detected a catastrophic phase transition point in the frequency of faulty conditions with respect to the data rate of the communication network. The CI with interdependencies safely operate only in a region of the system parameters that are far from this catastrophic transition point. To this aim, the introduction of adequate safety margins with respect to the data rate of the communication network is suggested.

In summary, the main contributions to the knowledge generated for the network general case are: (i) the extension of the purely-topological complex network modeling framework by inclusion of relevant physical aspects of the system; (ii) the use of the frequency of load-shedding actions due to excessive communication time delay as a quantitative indicator for quantifying the vulnerability of the smart grid; (iii) the detection of a catastrophic phase transition with respect to the data rate of the communication network in the smart grid.


F Navarro
email: fabien.navarro@math.unicaen.fr

M J Fadili

C Chesneau

Adaptive parameter selection for block wavelet-thresholding deconvolution

Keywords: Inverse problem, Estimation parameters, Deconvolution, Signal reconstruction, Coloured noise

In this paper, we propose a data-driven block thresholding procedure for waveletbased non-blind deconvolution. The approach consists in appropriately writing the problem in the wavelet domain and then selecting both the block size and threshold parameter at each resolution level by minimizing Stein's unbiased risk estimate. The resulting algorithm is simple to implement and fast. Numerical illustrations are provided to assess the performance of the estimator.

INTRODUCTION

Problem statement

Let A : H → K be a known linear operator between the Hilbert spaces H and K endowed with the corresponding inner products and associated norms. Assume that we observe the data Y according to Y = Af + ǫW, (1) where f ∈ H is the unknown function we wish to recover, ǫ > 0 is the noise level and W is a zero-mean Gaussian white noise. We aim to estimate a function f from the observations Y .

Motivating example

In this paper, we consider the problem of deconvolution, but other examples can also be considered as well, e.g., density with error in variable, tomography, etc., if the operator possesses K the so-called Wavelet-Vaguelette Decomposition (WVD) of Donoho [1995]. The degradation process can be modelled as a blurring operator A which is characterised by the kernel g (i.e. also often referred to as a point spread function (PSF) or impulse response)

Af : Ω → R t → Af (t) = g ⊛ f (t) = Ω g(s)f (t -s)ds, t ∈ Ω
where Ω = [0, 1] and both f and g belong to H = L 2 per (Ω) = {h; h is 1-periodic on Ω and Ω |h(t)| 2 dt < ∞}. Deconvolution is a challenging ill-posed inverse problem.

Overview of previous work and contributions

There is already an extensive literature on deconvolution for the case where A is known (see, e.g., Donoho [1995],

Cavalier [2008[ ], Cavalier et al. [2004]], Cavalier [2002], Chesneau [2007], Johnstone [1999]).

In this work, we consider adaptive block estimators based on the Wavelet-Vaguelette Decomposition (WVD) of Donoho [1995]. Minimaxity of James-Stein block thresholding over Besov balls with fixed block size has been established for deconvolution in Chesneau et al. [2010]. The performance of Block thresholding estimators strongly depends on the block size L and threshold level λ. In the nonparametric regression setting, in order to choose these two key parameters in an optimal way, Cai and Zhou [2009] proposed an adaptive James-Stein block thresholding estimator whose parameters minimize the Stein's unbiased risk estimate (SURE) and established its minimax rates of convergence under the mean squared error over Besov balls. Johnstone [1999] described a SURE-based leveldependent term-by-term thresholding estimator for a class of linear inverse problems possessing a WVD. Recently, the SURE has proved to be a powerful tool for signal/image restoration Luisier et al. [2007], Pesquet et al. [2009], Vaiter et al. [2012]. For image denoising problems, a Stein risk estimator have been proposed in Peyré et al. [2011] to both adapt the block-sparsity structure and the threshold.

In this work, we extend several of the previous works to solve linear inverse problems, mainly deconvolution here, with wavelet block thresholding estimators. Our estimator is constructed using a periodized Meyer wavelet basis and two popular block thresholding estimators; namely Block-Soft and Block James-Stein (JS) (defined in the sequel). The corresponding optimal block size and threshold parameter are determined by minimizing the SURE after rewriting the model appropriately in the wavelet domain.

The paper is organized as follows. Section 2 briefly reviews wavelets and Besov balls and describes the block thresholding-based deconvolution estimator. Section 3 dis-cusses numerical results, before conclusions are drawn in Section 4.

NONLINEAR ESTIMATION VIA BLOCK THRESHOLDING

Wavelets and Besov balls

For the purpose of this paper, we use the periodized Meyer wavelet bases on L 2 per (Ω). Any function f ∈ L 2 per (Ω) can be expanded into a wavelet series

f (t) = 2 j 0 -1 k=0 α j0,k φ j0,k (t) + ∞ j=j0 2 j -1 k=0 β j,k ψ j,k (t), t ∈ Ω, where α j0,k = f, φ j0,k , β j,k = f, ψ j,k (2)
with f, g = Ω f g, and φ j0,k and ψ j,k are scaled and translated versions of the father and mother wavelets φ and ψ. The family {φ j0,k , k = 0, . . . , 2 j0 -1; ψ j,k , j ≥ j 0 , k = 0, . . . , 2 j -1} forms an orthonormal basis of L 2 per (Ω). Meyer wavelets are compactly supported in the Fourier domain with

D j = supp (φ j,k ) ℓ ⊂ [-2 j+1 c, 2 j+1 c], C j = supp (ψ j,k ) ℓ ⊂ [-2 j+2 c, -2 j c] ∪ [2 j c, 2 j+2 c],
where c = 2π/3, g ℓ = g, e ℓ with e ℓ = e 2πiℓt , and supp denotes the support. See Meyer [1992] for a detailed account.

We say that a function f belongs to the Besov ball B s p,r (M ) if and only if Ω |f (t)| p dt ≤ M , and there exists a constant M * > 0 (depending on M ) such that the associated wavelet coefficients (2) satisfy

2 j0(1/2-1/p)   2 j 0 -1 k=0 |α j0,k | p   1/p +    ∞ j=j0   2 j(s+1/2-1/p)   2 j -1 k=0 |β j,k | p   1/p    r    1/r ≤ M * ,
with a smoothness parameter 0 < s < q, and the norm parameters p and r. Besov spaces capture a variety of smoothness features in a function including spatially inhomogeneous behavior, see Meyer [1992].

Smoothness of the kernel g

In this study, we focus on the standard context of ordinary smoothness of g. That is, that there exist three constants, c g > 0, C g > 0 and δ > 1 such that, for any ℓ ∈ Z,

c g 1 (1 + ℓ 2 ) δ/2 ≤ |g ℓ | ≤ C g 1 (1 + ℓ 2 ) δ/2 .
(3)

This assumption controls the decay of the Fourier coefficients of g. It is a standard hypothesis usually adopted in the field of nonparametric estimation for deconvolution problems, see e.g. [START_REF] Cai | On adaptive wavelet estimation of a derivative and other related linear inverse problems[END_REF], Donoho [1995], Pensky and Vidakovic [1999]. k 2 j . . .

2 j 1 . . . 2 j 2 A j 0 2 j -1 K = 1 b j,1 L . . . K = 2 j /L
Fig. 1. Illustration of the blocks partition at a given scale.

Wavelet deconvolution in the Fourier domain

We set, f ℓ = f, e ℓ , W ℓ = W, e ℓ , where W ℓ are i.i.d. standard (complex-valued) Gaussian random variables.

Then, re-writing the model (1) in the Fourier domain

y ℓ = f ℓ g ℓ + ǫW ℓ , ℓ ∈ Z.
(4) Moreover, we can compute the wavelet transform via the Parseval formula

β j,k = ℓ∈Z f ℓ (ψ j,k ) ℓ ,
where (• j,k ) ℓ = • j,k , e ℓ , Then, combining the latter with (4), we get

ℓ∈Cj y ℓ g ℓ (ψ j,k ) ℓ = ℓ∈Cj f ℓ (ψ j,k ) ℓ + ǫ ℓ∈Cj e ℓ g ℓ (ψ j,k ) ℓ .
Thus, we can perform wavelet transform and deconvolution simultaneously and (4) reduces to a general sequence model with correleted noise

z j,k = β j,k + ǫW j,k , (5)
where

W j,k = ℓ∈Cj e ℓ g ℓ (ψ j,k ) ℓ ∼ N 0, ǫ 2 ℓ∈Cj |(ψ j,k ) ℓ | 2 |g ℓ | 2 .
Let us denote by σ j the standard deviations of the noise process W j,k and let λ j = λσ j be a sequence of thresholds to be applied to the wavelet coefficients at level j. In view of (3), we have at level j that σ j ∼ ǫ2 jδ . We define the estimator of the wavelet coefficients of f as

β * j,k = η * (z j,k , λ j ), (6)
where η * is either Block-Soft or Block-JS thresholding defined in (8).

Block deconvolution estimator

Let the observed signal Y be the vector of n observations Y 1 , . . . , Y n . Let j 1 = ⌊ln 2 (ln n)⌋ be the coarsest resolution level, and j 2 = ⌊(1/(2δ + 1)) ln 2 (n/ ln n)⌋, where, for any a ∈ R, ⌊a⌋ denotes the integer part of a. For any j ∈ {j 1 , . . . , j 2 }, let L ≥ 1 be the block size.

Let A j = {1, . . . , ⌊2 j L -1 ⌋} be the set indexing the blocks at resolution j. For each j, let {b j,K } K∈Aj be a uniform and disjoint open covering of {0, . . . , 2 j -1}, i.e. K∈Aj b j,K = {0, . . . , 2 j -1} and for any (K, K ′ ) ∈ A 2 j with K = K ′ , b j,K ∩ b j,K ′ = ∅ and Card(b j,K ) = L, where b j,K = {k ∈ {0, . . . , 2 j -1}; (K -1)L ≤ k ≤ KL -1} is the Kth block (see Fig. 1).

We define our block estimator f λ of f by

f λ (t) = 2 j 1 -1 k=0 α j1,k φ j1,k (t) + j2 j=j1 K∈Aj k∈bj,K β * j,k ψ j,k (t), (7) t ∈ [0, 1],
where for any resolution level j and position k ∈ b j,K within the Kth block, the wavelet coefficients of f are estimated via (6) with one of the rules

η Soft (z j,k , λ) = 1 - λ z bj,K 2 + z j,k , η JS (z j,k , λ) = 1 - λ 2 z bj,K 2 2 + z j,k , (8) with, z bj,K 2 2 = k∈bj,K |z j,k | 2
, and (a) + = max(a, 0).

α j1,k = ℓ∈Dj 1 y ℓ g ℓ (φ j,k ) ℓ is the empirical approximation coefficient.
Note that ( 7) is simple to implement and fast (thanks to the Fast Fourier Transform algorithm) which allows us to perform the selection procedure in a reasonable time.

Block thresholding estimators depend on both the block size L and threshold parameter λ which mainly determine the performance of the resulting estimator. So it is important to select these parameters in an optimal way.

Unbiased risk estimation for automatic parameter selection

Since the discrete wavelet transform is orthogonal, the risk of an estimator f of f will be the same as that of its wavelets coefficients. The risk R( f , f ) = E( ff 2 ) will measure the estimation error in the Mean Square Error sense. Since the blocks are also non-overlapping, an estimator of the risk is derived by summing the risk estimates over blocks and resolution levels j.

More precisely, since each of the thresholding rules in ( 8) is weakly differentiable with an essentially bounded gradient, Stein lemma allows to get an estimator of the risk in each block, which solely depends on the observation. Denoting z bj,K = (z j,k ) k∈bj,K , the SURE on each block b j,K is given by

J Soft (z bj,K , λ j , σ j ) = Lσ 2 j + (L z bj,K 2 -2Lσ 2 j )1 { z b j,K <λj } + (Lλ 2 j -2σ 2 j (L -1) λ j z bj,K )1 { z b j,K λj } J JS (z bj,K , λ j , σ j ) = Lσ 2 j + (L z bj,K 2 -2Lσ 2 j )1 { z b j,K <λj } + Lλ 2 j -2σ 2 j (L -2) z bj,K 2 /λ 2 j 1 { bj,K λj } ,
where 1 {•} denotes the indicator function of the event in its argument. Thus, the overall SURE of β * j,k is computed by simple summation

J * ({z j,k } j,k , λ) = j2 j=j1 K∈Aj k∈bj,K J * (z bj,K , λ j , σ j ) . (9)
This is indeed an unbiased estimator of the risk. Consequently, our approach can be used as a principled way to objectively choose the optimal parameters λ and L that minimize (9).

Deconvolution algorithm

The deconvolution algorithm can be summarized as follows:

Algorithm 1. deconvolution algorithm Parameters: The observed blurred and noisy signal Y , the PSF g. Initialization:

• Block size: L = ⌊ln n⌋.

• Coarsest decomposition scale j 1 = ⌊ln 2 (ln n)⌋.

Step one: Fourier and wavelet domain.

• Apply an inverse filtering on the data using FFT (i.e. 1/g ℓ ).

• Perform one Forward wavelet transform:

z j,k = ℓ∈Cj y ℓ g ℓ (ψ j,k ) ℓ .
• Compute and store the variance σ 2 j .

Step two: Compute the risk.

• For each possible dyadic block b j,K compute J * (z bj,K , λ j , σ j ) and sum over the scale to get J * ({z bj,K } , λ j , σ j ).

Step three: Optimal parameters.

• Get the optimal parameters λ and L by a grid search. Output: Get deconvolved signal f λ at the optimal λ and L.

SIMULATION EXPERIMENTS

Three test functions, representing different degrees of irregularity, were used. The signals have been blurred with a Laplacian PSF and were corrupted by a zero-mean white Gaussian noise. Fig. 2 depicts the SURE and true risk curves as a function of the block size and the threshold λ, for each thresholding rule, plotted at each scale with the optimal block size (1st row) and summed over scales for each block size used (2nd row). It can be observed that the SURE gives a very reliable estimate for the risk, and in turn, also a high-quality estimate of the optimal λ. The restoration quality can be assessed both visually and quantitatively in Fig. 3, which shows the restored signals at the optimally chosen λ and L, with an input BSNR = 30 dB.

CONCLUSION

In this paper, a data-driven block wavelet-based deconvolution estimator was presented. Its usefulness has been illustrated on automatic selection of the both the block size L and threshold level λ for signal deconvolution. Although we focused on one-dimensional convolution, the approach can handle other operators A and easily extended to a multi-dimensional setting (e.g. images). We reported some simulation experiments to support our findings. Our efforts are now directed towards establishing the theoretical minimax performance of this estimator.

  Fig. 2. SURE and risk as a function of the thresholdfor Soft (left) and JS (right) block thresholding both at each scale with the optimal block size (top) and summed over scales for each block size used (bottom).

Fig. 3 .

 3 Fig. 3. Top panel: Original and observed signal (top);Bottom panels: Block-Soft (left) and Block-JS (right) estimates at the optimal λ and L.

T. Cai and H. Zhou. A data-driven block thresholding approach to wavelet estimation. Ann. Stat., volume 37, pages 569-595, 2009. L. Cavalier. Nonparametric statistical inverse problems. Inverse Problems, volume 24 (2008). L. Cavalier and A.B. Tsybakov. Sharp adaptation for inverse problems with random noise. Probab. Theory and Related Fields, volume 123, pages 323-354, 2002. L. Cavalier, Y. Golubev, O. Lepski and A.B. Tsybakov. Block thresholding and sharp adaptive estimation in severely ill-posed inverse problems. Probab. Theory and Related Fields, volume 48, pages 426-446, 2004. C. Chesneau. Wavelet estimation via block thresholding: A minimax study under the L p risk. Stat. Sinica, volume 18, pages 1007-1024, 2008. C. Chesneau, M.J. Fadili, and J.-L. Starck. Stein block thresholding for wavelet-based image deconvolution. Elec. J. Statistics, volume 4, pages 415-435, 2010. D.L. Donoho. Nonlinear Solution of Inverse Problems by Wavelet-Vaguelette Decomposition. Appl. Comput. Harmon. Anal., volume 2, pages 101-126, 1995. J. Fan and J.Y. Koo. Wavelet deconvolution. IEEE Trans. Inf. Theory, volume 48, pages 734-747, 2002. J. Johannes. Deconvolution with unknown error distribution. Ann. Stat., volume 38, pages 2301-2323, 2009. I. Johnstone. Wavelet Shrinkage For Correlated Data And Inverse Problems: Adaptivity Results. Statistica Sinica, volume 9, pages 51-83, 1999. F. Luisier, T. Blu, and M. Unser, A new SURE approach to image denoising: Inter-scale orthonormal wavelet thresholding. IEEE Trans. Image Process., volume 16, pages 593-606, 2007. Y. Meyer. Wavelets and Operators. Cambridge University Press, Cambridge, 1992. M. Pensky and B. Vidakovic. Adaptive wavelet estimator for nonparametric density deconvolution. Ann. Stat., volume 27, pages 2033-2053, 1999. J.C. Pesquet, A. Benazza-Benyahia, and C. Chaux. A SURE approach for digital signal/image deconvolution problems. IEEE Trans. Signal Process., volume 57, pages 4616-4632, 2009. G. Peyré, J. Fadili and C. Chesneau. Adaptive Structured Block Sparsity Via Dyadic Partitioning. Proc. EUSIPCO 2011, pages 303-307, 2011. S. Vaiter, C. Deledalle, G. Peyré, J. Fadili, and C. Dossal. Local behavior of sparse analysis regularization: Applications to risk estimation. to appear in Appl. Comput. Harmon. Anal., 2012.


Piero Baraldi
email: piero.baraldi@polimi.it

Francesco Di Maio

Marco Rigamonti

Enrico Zio

Redouane Seraoui

Transients Analysis of a Nuclear Power Plant Component for Fault Diagnosis

We analyze signal data collected during 148 shut-down transients of a nuclear power plant (NPP) turbine for fault diagnosis. The aim is to identify groups of transients with similar characteristics, under the conjecture that faults and malfunctioning of the same type lead to a similar behavior of the measured signal transients. The fault diagnosis method is based on the combined use of an original fuzzy-based slope and similarity analysis with spectral clustering.

Introduction

Fault diagnosis can be seen as a classification problem in which specific classes of fault are associated to specific values of measured signals [START_REF] Zio | Neuro-fuzzy Pattern Classification for Fault Diagnosis in Nuclear Components[END_REF]. In general terms, fault diagnosis methods can be divided into two categories: model-based and data-driven (Venkatasubramian et al., 2003[START_REF] Sheng | Decision Tree-Based Methodology for High Impedance Fault Detection[END_REF]. Model-based methods need expert specific domain knowledge for model building while not guaranteeing the timely detection of faults [START_REF] Tian | Fault Detection of Oil Pump Based on Classify Support Vector Machine[END_REF], whereas data driven methods can satisfy a number of practical requirements, such as short calculation time and high accuracy (Zio, 2012;Baraldi, 2012a). We develop a new data driven method for identifying specific signal transient behaviors associated to different types of fault. A first task to be done is the extraction of signal features relevant for the transients classification. Several techniques of feature extraction can be used, such as pointwise difference (Baraldi et al., 2012b), correlation [START_REF] Guyon | An Introduction to Variable and Feature Selection[END_REF], Principal Component Analysis (PCA) [START_REF] Baraldi | Local Fusion of an Ensemble of Models for the Reconstruction of Faulty Signals[END_REF]. Problems may arise when dealing with misaligned signals [START_REF] Secchi | K-mean alignment for curve clustering[END_REF], i.e., amplified and/or delayed transients data. To overcome this, we develop an effective feature extraction technique for misaligned data, which we call Fuzzy-based slope analysis. For the clustering, we then propose to compute the similarity between signal trajectories by measuring the fuzzy similarity of the extracted features [START_REF] Zio | A Fuzzy Similarity-Based Method for Failure Detection and Recovery Time Estimation[END_REF]. A similarity graph (von [START_REF] Luxburg | A Tutorial on Spectral Clustering[END_REF] is built, in which each vertex represents a trajectory and the weight associated to the edge connecting two vertices is the value of (fuzzy) similarity between the two corresponding trajectories. Spectral clustering is finally applied in order to find an optimal partition of the graph (von [START_REF] Luxburg | A Tutorial on Spectral Clustering[END_REF]. The developed framework of analysis has been tested on a real industrial case regarding the shut-down transients of a turbine of a NPP. The rest of the paper is organized as follows: Section 2 states the problem and describes the available NPP turbine data; in Section 3 the proposed framework of analysis is described in details along with the application to the case study; in Section 4 the obtained results are shown and discussed in detail; finally, in Section 5, some conclusions and remarks are drawn.

NPP turbine shut-down transients

We consider the measured values of 70 signals taken at T=4500 time steps with a fixed sampling frequency during N=148 NPP turbine shut-down transients. The generic i-th transient, i=1,...,148, can be seen as a multidimensional trajectory into the 70-dimensional signal space, and represented by the matrix of values i X whose component i lk x represents the value of signal k taken at time t l . For fault diagnosis, we want to partition the 148 multidimensional trajectories into an a priori unknown number of clusters, C, each one containing trajectories of similar functional behavior, with the assumption that they have been originated by the same physical setting (nominal or faulty).

Data pre-treatment: reduction of the dimensionality of the signals space

Each i-th transient entails a very large amount of data (4500 simultaneous measurements of 70 signals) that challenge the identification of similarities between the transients in a signal space characterized by such a high dimensionality (Baraldi et al., 2012d). To overcome this problem, we have calculated the correlation among the signals and retained those with large correlation values, i.e., containing relevant information for clustering the transients. A spectral clustering algorithm applied to a correlation matrix of size [70,70] has allowed identifying six groups of signals characterized by a high degree of correlation among them and a low degree with the signals of other groups. The largest group is composed by 27 signals and includes the active power signal, two vibration signals and twenty four temperature signals.

For limitation of paper length, the analyses described in Sections 3 and 4 will take into account only this group of signals.

The framework for unsupervised clustering

In this Section, we introduce the basic method for fault diagnosis developed and applied to the shut-down transients of NPP turbines. In Section 3.1, the fuzzy-based slope analysis feature extraction technique is applied to the raw signal data; Section 3.2 describes the calculation of the fuzzy-similarity measure between transients represented by the extracted features; finally, Section 3.3 describes the application of the spectral clustering method to the obtained similarity matrix.

Feature extraction: Fuzzy-based slope analysis

The developed feature extraction technique is properly devised to overcome problems of misaligned transients. It consists of five steps, applied separately to each k-th signal, k=1,...,27:

1. Signal slope computation: divide the i-th transient of the k-th signal in R=9 intervals of length L= T/R=4500/9=500. The number of intervals has been found by trial-and-error. For each r-th interval, r=1,..,9, the slope i k r  of the signal is computed using the least squares regression method on the signal values in the considered interval [START_REF] Frank | Intermediate least squares regression method[END_REF].

2. Signal slope distribution: collect the slope values of the k-th signal in all the 148 transients and build separately the signal slope empirical distributions for the positive and the negative values. 3. Percentiles computation: calculate the 1 st and the 50 th percentiles of the negative slopes distribution (N1 and N50), and the 50 th and the 99 th percentiles of the positive slopes distribution (P50 and P99).

4. Fuzzy sets and membership functions creation: consider the signal slope,  , as a linguistic variable that can be approximated within a fuzzy framework defined by the fuzzy sets (terms) high negative, low negative, low positive and high positive. In this case, the membership functions which define the fuzzy sets are S=5 asymmetric, unevenly spaced triangular functions centered on the percentile values computed at step (3) and on the zero value (V0): we chose asymmetric, unevenly spaced triangular fuzzy sets, based on the computed percentiles values, because they only depend on the distribution of the available signal slope values, and allow to describe the signals slope behavior better than symmetric triangular fuzzy sets, which are centered on a priori established values and may lead to the creation of empty fuzzy sets [START_REF] Baldwin | Asymmetric Triangular Fuzzy Sets for Classification Models[END_REF].

Once the membership functions have been defined, each computed signal slope can be characterized by its degrees of membership to the fuzzy sets which represent the level to which the slope verifies the extent properties of the set, as shown in Figure 2.

   R r i k r s i k s R 1 1   (1)
In this way, since five different membership values are obtained for each signal, the matrix

Pattern matching by fuzzy-similarity

The similarity between transients i and j can be computed by considering the corresponding vectors of extracted features i Y and j Y . To this aim, a fuzzy similarity measure is considered to determine the degree of closeness of the two trajectories with reference to the pointwise difference between the corresponding feature values [START_REF] Zio | A Fuzzy Similarity-Based Method for Failure Detection and Recovery Time Estimation[END_REF]Baraldi et al., 2012b). Without loss of generality, let us define the pointwise difference ij  between the trajectories i Y and j Y :

      135 1 2 f j f i f ij y y  (2)
The pointwise difference of the two trajectories is then evaluated with reference to an "approximately zero" fuzzy set (FS) specified by a function which maps ij  into a value ij  of membership to the condition of "approximately zero": values of ij  close to 0 indicate that the signal evolutions in the two transients i and j are very different, whereas values close to 1 indicate high similarity. In this work, the following bell-shaped function is used:

2 2    ij ij e   (3)
The arbitrary parameter σ be set by the analyst to shape the desired interpretation of similarity into the fuzzy set: the larger the value of σ, the narrower the fuzzy set and the stronger the definition of similarity [START_REF] Zio | A Fuzzy Similarity-Based Method for Failure Detection and Recovery Time Estimation[END_REF]. Applying eq. ( 2) and eq. ( 3) with an optimized value of the bell-shaped function parameter = 0.55, we have built the similarity matrix W of size [148,148], whose values represent the functional similarity between transients.

Transient clustering by spectral clustering

From the matrix W of size [148,148] a similarity graph G = (V,E) is constructed, where each vertex vi represents the i-th trajectory and the weight associated to the edge pij connecting the two vertices i and j is the similarity value ij [START_REF] Alpert | Spectral partitioning: the more eigenvectors, the better[END_REF]. The original problem of identifying groups of similar trajectories can then be reformulated as that of finding a partition of the similarity graph such that the edges connecting elements of different groups have small weights and the edges connecting elements within a group have large weights. Figure 3 shows the 148 eigenvalues of the Random walk Laplacian matrix rw L obtained by applying the spectral analysis method to matrix W (Baraldi et al., 2012b). Since the first four eigenvalues are closer to zero and the gap between the fourth and the fifth is larger, the number of clusters C is set equal to 4, according to the eigengap theory (von [START_REF] Luxburg | A Tutorial on Spectral Clustering[END_REF].

The relevant information to be used for partitioning the graph (i.e., clustering the transients) consists into the eigenvectors 4 , 2 , 1 ...,u u u associated to the 4 smallest eigenvalues [START_REF] Baraldi | Unsupervised Clustering for Fault Diagnosis in Nuclear Power Plant Components and Systems[END_REF]. The problem of clustering the 148 trajectories i X is now reduced to the problem of finding four clusters among the 148 4-dimensional vectors, where for each i-th transient

  i i i i u u u u
- BASED SLOPE ANALYSIS Signal Measurement                i Z i i i i i i i i i i i x x x x x x x x x x x X ,
27 5 2 3 2 2 2 1 1 5 1 4 1 3 1 2 1 1 i i i i i i i i i i Y           i=1,…, 148 i=1,…, 148
th transient to the c-th cluster, c=1, 2,3,4 (Baraldi et al., 2012c). The transient with the largest value of membership to a cluster is named prototypical trajectory, and its functional behavior can be taken as most characteristic of the cluster.

Clustering Results

With respect to the matrix of values ic  , two categories of transients have been distinguished:

-Transients assigned to a cluster with high confidence (the largest membership value is larger than 0.7) are hereafter labelled "assigned". -Transients not assigned to any cluster with enough confidence are hereafter labelled "not assigned". As an example, Figure 4 shows the evolution of the signal 3 in the "assigned" transients of the four identified clusters. Clusters 1 and 4 mainly differ in the rate of decrease of the signals: those belonging to cluster 1 are characterized by a sharper, less smooth decrease than those of cluster 4. Although the distinction between these two groups depends on the rate of decrease of the signals, their functional behaviour is almost the same. Furthermore, it is possible to note that, even if the distinction between the more characteristic transients belonging to cluster 1 and cluster 4 is clearly marked, there are several transients at the border between the two clusters for which the signal behavior seems to be very similar. On the contrary, clusters 2 and 3 show peculiar functional behaviors which allow distinguishing them from other transients: those belonging to cluster 2 exhibit a smooth rise followed by a decrease of the signal, whereas those belonging to cluster 3 show a short and more marked decrease followed by a rise. The low degrees of membership of the "not assigned" transients can be due to two different causes:

the "border-line" transients could belong to more than one cluster, i.e., being at the border of two clusters they share some characteristics with both. -The "outliers" transients are very different to those of the four identified clusters, i.e. their signal behaviours are not similar to those of the transients of any of the identified clusters. In order to identify if the "not assigned" transients are "border-line" or "outliers", we have used an autoassociative method, previously developed in the context of fault detection (Baraldi et al., 2012d). The idea is to develop an auto-associative model using a training set formed by the patterns of the "assigned" transient (Figure 5). When the "not assigned" transients are fed to the developed auto-associative model AAKR (Baraldi et al., 2012d) two cases may arise: the signal reconstructions (output of the AAKR model) are similar to the signal behaviour of the "not assigned" transient (input to the AAKR model) or they are different. In the former case, we can conclude that the "not assigned" transient is similar to those used for the model training, i.e. the transients of the four clusters, and, thus, the transient is at the border of two clusters. Contrarily, in case of remarkable difference between the reconstruction and the "not assigned" transient, we can conclude that the transient is atypical with respect to the training transients and, thus, it is an outlier characterized by anomalous signal behaviors. According to this analysis, 6 out of the 15 transients are "outliers", whereas the remaining 9 transients are "border-line". -

Conclusions

In this work, we have developed a framework of unsupervised classification of transients for identifying different behaviours associated to different faults. We have preceded the clustering task with a feature extraction capable of treating misaligned data. The proposed fuzzy-based slope analysis for feature extraction is able to capture functional shapes in case of misaligned transients and seems to be a general and reliable technique. We have then computed the similarity between the trajectories by considering the fuzzy similarity of the extracted features, obtaining a squared similarity matrix. The transients clustering is then obtained by the application of the spectral analysis and the following FCM algorithm to the similarity matrix. The application of the proposed framework to a real industrial case study concerning 148 shutdown transients of a NPP turbine has demonstrated the feasibility of the approach. Four different groups of transients have been identified: two groups mainly differ for the rate of decrease of signals correlated to the turbine speed, whereas the other two groups are characterized by transients with anomalous behaviors of some signals, which should be analyzed by NPPs experts to understand whether they are associated to malfunctioning, faults or peculiar working conditions. The applied methodology has also identified the presence of 6 atypical transients (outliers) which have not been associated to any group with enough confidence, which can be considered also as a symptom of malfunctionings or faults.

5.

  Feature extraction: for each i-th transient compute the mean membership i k s  of the k-th signal to each of the s fuzzy sets:

Figure 1 :

 1 Figure 1: Sketch of the proposed technique for the fuzzy-based slope analysis for the feature extraction

  representation of i X . The Fuzzy C-Means (FCM) clustering algorithm fed with i u provides the memberships ic  of the i-

FUZZY

Figure 2 :Figure 3 :

 23 Figure 2: Membership functions of the linguistic variable slope

Figure 4 :Figure 6 ;

 46 Figure 4: Evolutions of the signal 3 in the transients assigned to the four clusters