A Boolean Networks Approach to Modeling and Resilience Analysis of Interdependent Critical Infrastructures

In this article, we propose a dynamic Boolean network‐type mathematical representation of networked critical infrastructures under stress. Our formulation is obtained as a threshold‐based approximation of a discrete‐time switched system meant to describe failure and recovery phases that may span considerably different time periods. The formalism can incorporate mitigating factors, resource constraints, and allocation strategies affecting the response of the system and the unfolding of restoration processes. As such, the tool may inform both contingency planning and consequence assessment. We also discuss an illustrative case study concerning the recovery process of selected infrastructures of the Thessaloniki Port in an earthquake scenario. Simulations are performed for assigned parameter distributions associated with the different operational modes. Finally, a performance loss indicator is computed to assess the overall service perturbation affecting the system. For the needs of this analysis, fragility curves have been selected from the literature, whereas restoration curves, recovery priorities, classification of components, and other modeling variables have been defined in collaboration with the local Port Authority.


INTRODUCTION
The vulnerability of critical infrastructures (CIs), their protection, and management throughout adverse cir-This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. cumstances have been addressed by a number of policy documents, including President's Commission on Critical Infrastructure Protection (1997), European Commission (2004), and Obama (2013). The research community has recognized complexity and heterogeneity of mutual interactions among the key aspects related to the analysis of modern CIs (Rinaldi et al., 2001;Kröger, 2008;Vespignani, 2010;Ouyang, 2014). Relevant factors investigated by recent studies also include human behavior (Subotic et al., 2007;Oxstrand and Sylvander, 2010;Boring and Gertman, 2010;Blythe, 2012) and systems' adaptivity (Brown et al., 2004;Negenborn et al., 2010). Assessing emerging behaviors is significant in this process: Bak and Paczuski (1995) introduced the concept of self-organized criticality, referring to large-scale systems reaching highly interactive critical states where even minor perturbations may lead to sizable damages and avalanches; Pescaroli and Alexander (2016) acknowledged the importance of vulnerability paths toward a more thorough understanding of cascading disasters and consequences.
A concept overarching many of these topics is that of resilience, defined by the National Academy of Sciences as the ability "to plan and prepare for, absorb, respond to, and recover from disasters and adapt to new conditions" (National Research Council, 2012). Further definitions can be found in today's steadily growing literature; among them, resilience has also been referred to as "the ability to reduce effectively both the magnitude and duration of the deviation from targeted system performance levels" (Vugrin et al., 2010). The emphasis on the importance of a dynamic and operational assessment of resilience is predominant in recent research. It is mirrored, for instance, in the robustness, redundancy, resourcefulness, and rapidity factors found in Bruneau et al. (2003). It interrelates with the absorptive, adaptive, and restorative capacity indicators studied in Vugrin et al. (2013). Moreover, a number of novel contributions are addressing aspects such as the development of resilience indexes and optimization techniques, see, for instance, Bozza et al. (2017), Ouyang (2017), Ouyang and Fang (2017), and Wang et al. (2017). Resilience analysis and enhancement also display interesting relationships with sustainability and sustainable design (Redman, 2014;Wang and Adeli, 2014;Xu et al., 2015;Rafiei and Adeli, 2016;Wang et al., 2017).
To derive resilience metrics useful to assess a system's performance under stress, the literature often refers to functionality and affine concepts, see also Nayak and Turnquist (2016). Building on the resilience triangle representation considered in Bruneau et al. (2003) and successive literature, we can illustrate the idea through a functionality measure F(t) having a nominal value F 0 and deviating from it under the influence of a perturbation. This may, for instance, lead to some minimum level of performance F min and be followed by recovery actions aiming to correct the performance deterioration and report the system to the pre-event status, see Figure 1. Based on this type of representation, a number of resilience indicators have been proposed. These include, for instance, the well-established joint assessment of systemic impact ("the difference between a targeted system performance level and an actual system performance following a disruptive event") and total recovery effort ("the amount of resources expended during recovery processes following the disruption") (Vugrin et al., 2010). Ganin et al. (2016) map planning, absorption, recovery, and adaptation to different time phases associated to this type of representation. Moreover, it claims that resilience assessment should "identify the critical functionality of a system and evaluate the tem-poral profile of system recovery in response to adverse events," and introduces an associated resilience metric. Notably, as stated in Nayak and Turnquist (2016), the deterioration phase and recovery phase durations ((t 1 − t 0 ) and (t 2 − t 1 ) in Figure 1, respectively) may differ consistently, as in many practical situations, the damage processes targeted by this representation are much more rapid than the subsequent repair actions, see also Ouyang et al. (2012) and related references.
As the scale of the system under study grows and the represented entities become heterogeneous, the use of network-based approaches is gaining increasing attention in today's literature. Indeed, a vast array of such techniques has been developed to address different aspects of CI modeling, analysis, and optimization. In this respect, Ouyang (2014) identifies the two categories of topology-based methods and flow-based methods and provides an account of recent achievements in the two domains. Interestingly, to the purpose of the present article, often resilience assessment of the resulting network can be managed through a discretization of the functionality status of nodes/edges. Correspondingly, resilience metrics can focus, for instance, on the number of nodes having a given status over the analysis horizon (Ganin et al., 2016) and take into consideration graph-mediated dependencies. Furthermore, discretization can allow focalizing the attention on mutual dependencies and enriching their description, see, for instance, the case of operational/nonoperational nodes with conjunctive and disjunctive dependencies introduced in Sen et al. (2014) and expressed through Boolean formulas.
Another trait often found in network-based CI analysis literature is the conceptual distinction between the internal functionality of nodes and their ability to support the working of the rest of the network, both of which, in turn, may also depend on the status of upstream nodes. For instance, in Trucco et al. (2012), functionality and operability variables interplay in determining the ability of each node to provide service, also allowing dynamical models to smoothly embed further resilience factors such as buffering components and other mitigating elements. Another example can be found in the case of dynamic input-output inoperability models, see, for instance, the operability and production operability notions proposed in Barker and Santos (2010). In related work (Barker and Santos, 2010;Galbusera et al., 2016), the presence of inventories loosens the coupling among the inoperability states of different industries and consequently influences shock propagation. More broadly, in network analysis, the articulation of internal (continuous/discrete) dynamics and external (discrete) behavior has been exploited in various ways to study propagation phenomena, for instance, through threshold-related inner-outer states. In the kinetic logic approach developed in Thomas (1973) for the modeling of genetic control circuits, internal node dynamics is absorbed into a discrete, symbolic description through thresholding, whereas node interactions are mediated by a graph with an associated set of logical equations. Another notable example is provided by the Abelian sandpile model (or Bak-Tang-Wisenfeld model) proposed in Bak et al. (1987). This represents a cellular automaton used to stylize cascading phenomena through the symbolism of grains of sand dropped over the nodes in a graph. A thresholding rule is then fixed so that as the number of grains on a given node exceeds a fixed amount, they topple to adjacent nodes. Applications of this model range from banking (Haldane and May, 2011) to natural catastrophes (Malamud et al., 1998;Hergarten, 2003) and critical events on the electrical grid (Dobson et al., 2007;Brummitt et al., 2012). Furthermore, thresholding is taken into account in the Demon model discussed in Carreras et al. (2014), where operating, failing, and failed modes are specified for network components and probabilistic rules are associated to establish propagation if at least one of a node's nearest components is failing. The latter reference also highlights the ability of the model to incorporate mono-and bidirectional couplings, fully symmetric or asymmetric coupling strengths, homogeneous and heterogeneous strengths, negative or positive reinforcement. In Di Muro et al. (2017), the concepts of internal and external functionality rules are specified for intra-and internetwork dependencies. Moreover, a threshold-based rule is exploited to relate a node's functionality to the number of functional neighbors.
Interestingly, thresholding has also been used to formulate dynamic Boolean network (BN) representations starting from dynamical systems described in terms of ordinary differential equations (ODEs). In Ghil et al. (2008), thresholding is introduced as a way to bring a continuous real-valued state variables representation to a Boolean domain. Furthermore, in Stötzel et al. (2015), an Euler-based method is proposed to pass from an ODE-based description of a dynamical system to a dynamical BN approximation. As emphasized in latter reference, some key advantages of the operation reside in the fact that, in spite of some coarseness introduced by this approximation, BNs can facilitate global analysis. This may be especially useful when the transformation is done in such a way to maintain key structural features of the original model (e.g., the presence of cyclic attractors).
As we will discuss in this article, dynamical BN can also support the study of the different phases a CI system goes through during and after a critical event. This aspect has also been observed, for instance, in the domain of economic damage propagation assessment, where Boolean delay equations have been employed (Coluzzi et al., 2011). In our study, we will start from a switched continuous representation of the system, wherein nodes assume different operating modes according to the current functional status (e.g., decay and recovery phases illustrated in Figure 1). The formalism incorporates the ability to specify binary dependencies between nodes, mediated by the concept of thresholding. As a further step, we will propose a method to pass to an approximated delayed BN representation, able to take into account the dynamical features of the different switching modes on a node basis. This allows coping with the differences between the phases mentioned above. At the same time, the discrete-domain treatment introduced by such a BN representation enables to describe insightfully various types of logical relationships, to introduce mitigating factors, and to represent the recovery actions in the awareness of priorities and resource allocation constraints. It has to be observed that the introduced BN formulation opens opportunities for the application of interesting techniques proposed in the recent literature. Notably, these include an extensive set of BN analysis and control methods formulated on top of algebraic state-space representations based on semitensor matrix products (Cheng et al., 2010).
We also illustrate the application of our approach to a case study wherein we describe the postearthquake recovery process of a complex, multidomain port infrastructure network endowed with buffering elements and subject to resource constraints. Thus, the proposed analysis is also related to the systemic risk assessment of interconnected infrastructures exposed to different natural hazards. In this context, different methodologies have been proposed, which can be distinguished according to the scale (city, regional, and national) and the objectives of infrastructure stakeholders (emergency, recovery, or mitigation planning), see, e.g., Adachi andEllingwood (2009), Cavalieri et al. (2014), Esposito et al. (2015), Argyroudis et al. (2015), and Franchin and Cavalieri (2015). The main computational modules of these approaches involve models for hazard, physical damage, and system's performance analysis. The outcomes of these studies include the distribution of damage and performance loss for single scenario events or for stochastically generated sets of events in case of probabilistic analysis. The results provide input for the resilience analysis of the infrastructure on the basis of efficient allocation of resources and disaster management.
The rest of the article is organized as follows: in Section 2, we propose the formulation of the switched, thresholded dynamical model upon which we base our analysis; Section 3 introduces a method for its approximated representation in terms of BNs; then, in Section 4, we discuss the considered case study; Section 5 provides associated numerical results and is followed by some concluding remarks.
Notation. Symbols ¬, ∧, and ∨ denote Boolean NOT, AND, and OR functions, respectively; and represent n-ary AND and OR,i.e.,i∈{1,...,n} is the unit step function, conventionally equal to 1 for a nonnegative argument and 0 otherwise. For discrete-time index k, we will use the notation k + = k + 1.

MODEL FORMULATION
In this section, we introduce our model formulation, which takes into consideration the following aspects: node functionality, articulated in inner and outer terms to reflect the internal status of each node as well as the node-to-node relationships; serviceability, which expresses the service capability of the nodes with respect to the end consumers.

Inner and outer functionality
Network theory has been used extensively for CI representation, to describe both the internal structure of distributed infrastructures and mutual interdependencies. On the latter aspect, a recent review of various treatments of CI interdependencies involving network theory concepts has been proposed in Banerjee et al. (2017).
To formalize our modeling approach, we first introduce the notion of directed functionality graph G = (V, E), where V is a set of n nodes and E ⊆ V × V is the edge set. Term e i j ∈ E describes an edge from node i to j, and self-loops are also admitted. The set of neighbors In our study, nodes in G will be used to represent infrastructures as dynamical systems and their dependence on other components of the CI network will be specified through the corresponding edge topology. Dependencies, in particular, will be mediated by Boolean relationships determining the switching of nodes to different operational phases. Therefore, they will be used to encode both perturbation propagation patterns and decisions such as recovery strategies. In this sense, the introduction of self-loops allows, for instance, the implementation of logical rules depending on the status of the node itself.
To each node, we associate the following key properties: r Inner functionality: This may generally be intended as a composite indicator of structural integrity and, for instance, the integrity of internal processes of a specified asset; it will be expressed by means of a (normalized) functionality state variable x i (t) ∈ [0, 1], where 1 denotes full functionality and 0 its opposite; this variable will reflect, at the level of each node, the notion of functionality presented in the introduction.
r Outer functionality: This quantity discriminates the ability or inability of a node to ensure an adequate level of internal functionality to be active with respect to the rest of the network; in our representation, it will be described by means of a variable y i (t) ∈ {0, 1}, per node, obtained from x i (t) through thresholding.
Observe that the internal and external functionality rules found in Di Muro et al. (2017) map each node's functionality to the status of neighbors belonging to either the same network or another interdependent one. Instead, here, inner functionality reflects the internal status of a node, whereas outer functionality relates to the status that the node exposes. At the same time, similarly to the approach found in the latter reference, a node's exposed functionality is discrete in our formulation. Moreover, both inner and outer functionality can depend on the (outer) functionality of other nodes, as detailed below.
We are now in a position to introduce the following discrete-time switched model for each node i ∈ V: Herein, . . , N i } is a switching signal that selects at each time instant the active mode and will mediate node-to-node dependencies based on the following assignment: where binary signal w i (·) : N 0 → {0, 1} discriminates the absence or presence (0 or 1, respectively) of an external perturbation on node i at time k. Correspondingly, φ i (·, ·) : {0, 1} V(i)+1 → K i will determine the active mode at each time instant as a function of both node i's upstream graph dependencies and external perturbations.
Observe that, for simplicity, we assumed that threshold values are specific to each node, but not dependent on downstream nodes. More complex situations can be handled by techniques such as node splitting and arraytype definitions of outer functionality.
From now on, we shall assume |K i | = 2, ∀i ∈ V, where mode 1 will be used to describe the performance degradation phase of a node subject to a perturbation and mode 2 for the normal operation and recovery phases. Accordingly, we will have x * i,1 ∈ [0, 1), to describe different node-dependent perturbed equilibria, and x * i,2 = 1, ∀i ∈ V, to characterize a fully functional inner status for the network. (1) and (2) target a CI network subject to sudden changes in the dynamics under abrupt external perturbations on one side and to changes in the operational logic applied to control the network on the other. Function φ i (·, ·), in particular, can be used to encode both aspects. In this sense, observe that, through (2), both inner and outer functionalities of each node i can be affected by the outer functionality of upstream nodes. This allows, for instance, describing control actions taking into account factors such as resource limitations, prioritization constraints, and service dependencies among different CIs.

Serviceability
Based on the observations made in Remark 1, (2) may embed corrective measures, including, for instance, some reconstruction/reconfiguration component. At the same time, to complete our representation of the CI network in terms of its service provision to the end consumers, we introduce the concept of serviceability z i (t), ∀i ∈ V. This quantity is kept separated from outer functionality to describe, in particular, the presence of compensation measures meant to somehow mitigate the functionality drop and to damp the effects of perturbation on the delivered service. Therefore, serviceability allows some degree of separation between the internal working of the CI network and the output performance.
We now introduce a service graphḠ = (V,Ē), wherē V = V,Ē defines the service chains active in the considered CI system, and the set of neighbors of i ∈V inḠ is denoted asV(i). Then, for all i ∈V, we express serviceability as We will denote the model for node i expressed by Equations (1) and (4) as M i .
Mitigating factors may include, for instance, buffering elements, allowing a given asset to (temporarily) provide service under perturbed operational statuses. In this perspective, formula (2) may be useful to express, for instance, the activation of specific corrective control actions such as a recovery process subject to the availability of specific recovery resources. On the other side, (4) might be associated to the representation of compensative actions complementing corrective actions.

BOOLEAN NETWORK REPRESENTATION
The model structure specified in the previous section allows to accommodate logical constraints through the definition of switching rules (2) for each node. At the same time, the inner dynamics of each node has a dependency on other nodes mediated by the respective switching signal. These features allow the conversion of the proposed model to an approximate BN representation. In more detail, the objective is to start from M i and formulate the following approximated delayed BN modelM i , for each node i ∈ V: where D i = {0, . . . , δ max i } and δ max i ∈ N 0 is a maximum delay to be specified. D i will affect the approximation of functionality in terms of our delayed BN representation; thus, it is kept separate from T i as this relates to serviceability, instead. Observe that the structure of functions φ i (·, ·) and ψ i (·) is preserved in this transformation.
A method to derive the mentioned representation M i from M i will be now introduced. The objective, in particular, is the synthesis of function γ i (·), ∀i. To this end, for each node i ∈ V, we consider the set of initial states x seeds i corresponding to the steady states defined by f i, j (·), ∀ j ∈ K i .
A search routine to accomplish this task is outlined in Algorithms 1 and 2. Function BN forestscan is meant to determine, through subroutine BN treescan , switching signal patterns leading to threshold crossing for each node in the network and each root in x seeds i . These will be then exploited as conditions for the firing of a mutation in the Boolean state of the node, in the corresponding BN approximation.
In the proposed algorithms, constant δ max i specifies the upper bound to the tree search depth and is subject to tuning based on the system under consideration. Furthermore, ] serves as the threshold crossing indicator for node i and for a spec-ified seed. In particular, it discriminates whether the computed state sequence at a given tree exploration stage crossed the thresholdx i , for a given i ∈ V. Furthermore, function BN treebound(x i , x (+) i , x seed i ) introduces termination criteria to bound exploration when based on already found solutions and prefix switching patterns. Optimizations are also possible to accommodate reduced search spaces based on criteria such as path constraints and node-specific dynamical features. In particular, the search procedure can be simplified taking into account monotonicity considerations on the modes associated to the switched models.
The outcome of the search process is represented by ∈ X i is an ordered set of constrained-length state sequences leading to threshold crossing at node i, with X i +1 specify the associated tree depths and switching sequences, respectively.
The identified set of switching patterns leading to threshold crossing can be exploited toward the definition of γ i (·) in (5) based on present and past values of Boolean variablesŷ i (·) and ofû i (·), ∀i ∈ V. In particular, once triplet (X i , i , S i ) has been computed, the approximate BN model (5a) with x i (0) ∈ x seeds i , ∀i, can be specified as follows: Remark 2. Representation (5), together with the specification of term (5a) provided in (7) In the next sections, we will present the application of the proposed modeling approach and BN approximation framework to the analysis of recovery of selected infrastructures of the Thessaloniki Port in an earthquake scenario.

CASE STUDY: POSTEARTHQUAKE RECOVERY OF THESSALONIKI PORT INFRASTRUCTURES
The Thessaloniki Port is one of the most important ports in Southeast Europe and the largest transit-trade port in Greece. It occupies a total surface of 1.5 km 2 and includes six piers (five of them provide services to vessels), spreading over A representative seismic scenario is examined based on the results of a detailed microzonation study corresponding to a return period of 475 years or 10% probability of being exceeded in 50 years (Pitilakis et al., 2010), considering both ground shaking and ground failure due to liquefaction. Seismic forces are provided in terms of peak ground acceleration (PGA), whereas ground failure is expressed in terms of permanent ground deformation (PGD) at the components site. PGA values in the port area vary from 0.20 g to 0.30 g, whereas PGD is estimated in the order of 6.0-19.0 cm. Given the spatial distribution of the earthquake intensity and typological features of CI elements, damage probability at component level is assessed through fragility functions (National Institute of Building Sciences (NIBS), 2004). The functionality level (i.e., 0/1) for each component is assigned based on the associated damage level (i.e., the component is assumed to be fully functional if it sustains no or minor damage).
In this context, we consider a failure scenario involving the disruption of electrical substations and cranes. The different typologies of elements and their dependencies have been identified in collaboration with the local Port Authority using the SYNER-G taxonomy (Pitilakis et al., 2014). Threshold values, importance classifications, recovery priorities, and buffering capabilities have also been defined in collaboration with the Port Authority according to its practice or/and based on empirical choices.

Description of the components and their interaction
In this subsection, we introduce a description of the system's components (substations and cranes) in terms of their grouping, functional interlinking, recovery priorities, and recovery resource allocation. The considered infrastructure network consists of and V c = {12, . . . , 51} index electrical substations and cranes, respectively, whereas E will be specified through the switching rules below in this section; r a service graphḠ = (V,Ē) with no loops, illustrated in Figure 2 together with the geographic displacement of assets.
Service interdependencies are defined by specifying the substation (node) that supplies electric power to each crane or end-user substation. These are intended as substations used to supply power to other uses, e.g., container storage area, building facilities. In most cases, given the limited availability of resources, a prioritization of postearthquake restoration efforts is necessary to enhance system resilience and minimize the time needed for full recovery. This procedure has to take into account specific characteristics and functioning of components as well as the whole system layout features. Next, in this section, we will introduce the criteria associated to the implementation of recovery actions in the present case study.

Electrical substations.
They are partly used to serve the cranes and partly for other purposes related to the general operation of the port. Some of them (enduser substations) do not serve cranes, whereas they have an upstream dependency from other substations. It is noted that the distribution lines between the substations are installed underground and in this case study are considered as nonvulnerable. In addition, interdependencies with the external network are not considered (i.e., a continuous power supply to the substations is assumed).
Part of the substations have backup power supply capabilities able to guarantee some buffering time τ b i > 0 before disservice. The importance of different substations is classified by an integer-valued function cl(i), ∀i ∈ V s , taking values from 1 to 3 (highest to lowest importance) and partitioning V s into subsets V ( j) s , j = 1, 2, 3. As restoration times are in general short (few  days) compared to the case of cranes, according to expert opinion, the recovery of nonfunctional components should be prioritized by a combination of the classification of importance and the existence or not of backup appliances, as specified in Table 1.
s the partitions of set V s by recovery priority, ∀ j ∈ {1, . . . , 4}. In Table 2, we list the nodes belonging to each of these subsets. Finally, the number of available crews is assumed to be sufficient to restore in parallel all components within the same priority level and the buffering elements have by assumption negligible restoration time.

Cranes.
Cargo handling equipment has nonanchored components without backup power supply. Four gantry cranes are used for container loading-unloading services and located in the western part of the sixth pier. Cranes have only upstream functional dependencies from substations, due to their need for electricity to operate. In this case, the importance of cranes is classified by function cl(·) with values ranging from 1 to 5, and the corresponding partitions of V c are denoted as V ( j) c , j = 1, . . . , 5. In each of the five classes, there is a minimum number of cranes that should be functional to ensure a minimum performance level of the port in the case of an unexpected event and to prevent a total breakdown of the port's functionality (in terms of served ships and transported cargos). This requirement is translated here into assuming that the recovery process will include two stages: in the first one, the objective is to ensure that the number of operational cranes in each classification group j reaches the minimum threshold θ 1 j = 2 3 |V ( j) c | ; the second stage is devoted to the completion of the recovery process, i.e., to reaching the full operability threshold θ 2 j = |V ( j) c |, ∀ j. Accordingly, we can specify V c as the node sets associated to the different importance classes in the first and second stages.  Table 3 reports the established recovery priorities for the different crane classes throughout the two considered recovery stages. It is noted that the prioritization considers the location of cranes, the operational depth of the pier, and the type of crane. In both stages, thus, higher priority is given to components with higher importance, until the target threshold for the current stage is achieved.
Recovery priority groups are then established in terms of sets V [ j] c , j ∈ {1, . . . , 10}, as described in Table 4.
In order to further specify both intra-and intergroup recovery prioritization, the following set of rules is established: r Within the same classification group, the prioritization will be based on the ordering of the cranes' node indexes (low to high). r There are two crews available to work simultaneously toward recovery. This means that up to two cranes, Table 4 Definition of crane recovery priority groups and associated number of elements  not necessarily belonging to the same classification group, can be repaired simultaneously.
r The Port Authority has two mobile cranes that can be used during the recovery operations to substitute any damaged cranes, with the exception of those with classification of importance equal to 1. Therefore, mobile cranes can help compensating the lack of functionality during the stages of the recovery process, still taking into account the established classification of importance toward allocation.

Restoration curves.
The ability of all components (substations, cranes) to promptly recover after the impacting event is defined using restoration curves, which quantify the percentage of the component's functionality conditional on the damage state and the time after the earthquake occurrence. Site-specific restoration curves proposed in previous national projects ( SRM-LIFE, 2003UPGRADE, 2013UPGRADE, -2016 are used (Figures 3 and 4).

Switching rules
In this and the next subsections, we will introduce, for the described case study, a mathematical representation of r switching rules, involving both the external perturbation and the operational logics governing the recovery activities; r serviceability, wherein we will explicitly consider the role of buffering in electrical substations and temporary replacement by mobile cranes.
Next, switching rules are introduced for both substations and cranes.

Substations.
We can express the switching rules associated to each node i ∈ V [ j] s with priority j ∈ {1, . . . , 4} as follows: where c s ∈ {0, 1} is used to specify the availability of a crew to perform recovery of nodes in V s . Expression (8) specifies three alternative reasons leading to the activation of mode σ i (k) = 1, namely: r the presence of an external perturbation w i (k) affecting node i at time k; r the unavailability of recovery crew (term σ α i, j (k)), whereas node i is in need at time k; r the presence of higher priority nodes needing recovery actions (term σ β i, j (k)), whereas node i is in need at time k.
A graphical representation of the recovery sequencing implications expressed by these conditions for the substations subsystem is provided in Figure 5.

Cranes.
To express the more complex recovery resource allocation logics associated to the case of the cranes, we introduce auxiliary variables r [ j] i (k), for i = 1, 2 and j = 1, . . . , 10. These specify the allocation of the ith recovery resource to the crane group V [ j] c at time k, see Table 4. In particular, r [ j] 1 (k) determines the allocation of one of the two crews considered in our case study, whereas r [ j] 2 (k) will describe the simultaneous allocation of both to the same group.
In particular, we assign the following logical rules: and where c c,1 , c c,2 ∈ {0, 1} specify the availability of the two crews devoted to the recovery of nodes in V c . Terms r 1 [ j,α] (k) and r 2 [ j,α] (k) constrain the availability of resources at a given recovery priority group j (see Table 4) to the nonoccupancy status at higher priority layers. Furthermore, terms r 1 [ j,β] (k) and r 2 [ j,β] (k) link the assignment of resources to the different layers j to actual need, based on comparison with the relevant thresholds.
We are now in a position to introduce a switching rule Therein, observations similar to the case of (8) hold, with the notable exception that recovery resources here are staged, which affects σ α i (k) and σ β i (k); moreover, the latter term also takes into account the presence of multiple (two) recovery resources, allowing to test the allocation of the second resource when the first one is allocated.

Serviceability variables
In our representation, serviceability refers to the ability of the considered infrastructure network to provide service to end users. In particular, in our case study, nodes i ∈ V with in-degree greater than zero and zero outdegree are collected in a set of end-user service provider nodesV ⊂V. SetV can be partitioned asV =V s ∪V c , whereV s ⊂ V s collects the nodes associated to end-user substations andV c = V c all those associated to cranes. Serviceability throughout the considered case study is enhanced by the presence of electrical buffering elements, present for some of the substations, and mobile cranes.
In the case of substations, in particular, serviceability will be computed as follows: Observe that the serviceability of end-user substations requires jointly their own (buffered) outer functionality and service by any of the neighboring supporting substations, which in the considered case study configuration reduces to a single element in all cases.
The assignment of the two considered mobile cranes obeys principles similar to those applied for the recovery resources as described in (9) and (10). Nevertheless, assignment of mobile cranes is excluded for cranes belonging to recovery priority group 1. Consequently, we introduce the following quantities: where μ [ j] 1 (k) determines the allocation of one of the two mobile cranes to crane recovery priority group j, whereas μ [ j] 2 (k) describes the simultaneous allocation of both mobile cranes to the same group.
Based on the group allocation rules specified above, node allocation of mobile cranes to a specific node i ∈ V (h) c can be done through the following formula: where ρ i (k) = 1 implies that a mobile crane is assigned to node i at time k.
The availability of mobile cranes serves to foster the serviceability of the CI nodes, according to priority principles. Also, in the case of cranes, neighbors are represented by substations, whose serviceability has to be taken into account. Consequently, we formulate the following expression for the serviceability of crane elements: assuming that mobile cranes are served by the same substation as the cranes they replace.

NUMERICAL RESULTS
In this section, we report some simulation results obtained by implementing the proposed case study and applying the BN model approximation. All considered restoration scenarios refer to a fixed failure description involving a complete failure of the considered CI network, i.e., ∃k i > 0 such that x i (k i ) <x i , ∀i ∈ V. In particular, an abrupt perturbation to the nodes' performance was described by selecting models f i,1 (·) characterized by fast dynamics with perturbed equilibrium at x * i,1 <x i for each of them, consistently with the representations provided in Figures 3 and 4. Full availability of the recovery crew was assumed for both substations and cranes, so that c s , c c,1 , c c,2 = 1. Similarly, we assumed full availability of the mobile cranes, i.e., μ c,1 , μ c,2 = 1.
A Monte Carlo analysis was performed by considering a set of 500 scenarios for the considered system. These were formulated by introducing randomization over the dynamical modes f i,2 (·) associated to each node i in the CI network. In particular, the randomization was considered over suitable parameter distributions to take into account some possible variability in the duration of the recovery processes associated to the different nodes. In all scenarios, the system was initialized as fully functional.
Focusing on the crane compartment, we take into account the classification of the corresponding nodes into five classification groups (see Table 3). As all considered scenarios are fully perturbed, we have to expect that each group of cranes undergoes a recovery process involving two stages, to meet the threshold requirements specified in Table 4. Accordingly, in Figure 6, we report the distribution of the return times T [ j] c , j = 1, . . . , 10 (i.e., the time horizon needed for each crane recovery priority group to reach the associated threshold of operational nodes after being subject to the effects of the perturbation) over the random scenarios and organized per crane classification group. The figure confirms that  Table 4; red: second phase, i.e., j = 6, . . . , 10). higher priority is given, in each phase, to cranes belonging to lower indexed classification groups. The results are consistent with Figure 4 and involve different time scales with respect to the case of substations, whose recovery times are generally much shorter. In Figures 7 and 8, we detail, for a single scenario, further simulation results organized per crane category. These concern the allocation of the two available crane crews and substituting mobile cranes to different crane classes, respectively. Observe that the model correctly represents the sequencing of the resource allocation, which follows the priority principle. Furthermore, the model correctly describes the constraint specifying that substituting mobile cranes cannot be assigned to class 1, as per our standing assumptions. Furthermore, in Figure 9, we report the sums of the Boolean states of crane nodes, in time and per importance class, for the same scenario. Observe here that recovery processes overlap both between different class groups and within the same group.
Finally, we calculate an overall performance loss indicator for the port infrastructure in terms of the performance loss emerging from the operation of the terminal service providing nodesV over time horizon The distribution of the performance loss indicator over the considered random scenarios is reported in Figure 10. Based on such indicator, the infrastructure stakeholder can have an overview of the expected losses for given scenarios and existing restoration capacities, which provides a preliminary assessment of some port's resilience attributes and can be exploited to quantify economic losses. This kind of result may support a more efficient decision-making and management of recovery resources and processes.

CONCLUDING REMARKS
This article proposed a method for the representation of a networked CI in terms of a switched dynamical system, wherein functionality and serviceability are qualified and expressed taking into account mutual dependencies among nodes. Furthermore, we introduced a method for the approximation of the mentioned model in terms of a delayed BN representation, which is able to take into account different dynamical features of the nodes. The proposed approach was also applied to a case study describing the recovery process of a port. Toward this purpose, several modeling variables have been defined in collaboration with the Thessaloniki Port Authority. As part of the case study, we also computed a performance loss indicator that takes into account the loss of service associated to the terminal service providing nodes.
In the framework of interdependence types laid down in Rinaldi et al. (2001), this article proposes the exploitation of BNs to describe the service logic of a set of mutually connected components, as well as failure and recovery processes subject to constrained resources. Despite the reduced resolution introduced by the use of Boolean variables and the limitations imposed by the approximation process, the resulting formulation is able to accommodate complex sequencing constraints and restrictions on available resources in the context of infrastructures with heterogeneous failure and recovery dynamics. As such, the proposed modeling approach may be part of broader suite of techniques enabling the joint assessment of different interdependence classes and supporting the response to critical events.