-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathSysrel.py
165 lines (148 loc) · 8.44 KB
/
Sysrel.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
# -*- coding: utf-8 -*-
"""
Python module for system reliability. Contains the functions directly called by the main file:
- Load Network Data: Creates a graph object based on the line exposure and node damage information (i.e. prob. of failure).
Network fragility defines which node taxonomy corresponds to source and consumer nodes
- Evaluate System Loads: estimates the loads at nodes and edges, based on shortest path algorithm between source and consumer nodes
- Assign Initial Capacities: assign capacities to nodes and edges based on precomputed loads and a given safety factor
- Monte Carlo Simulation: simulates the hazard action and cascading effects for multiple random nodal damage states
- Compute output: based on the samples, returns sample of global values (total affected population) and probability of affectation for
each consumer area
Created on Tue Aug 13 10:52:00 2019
@author: hfrv2
"""
import numpy as np
import networkx as nx
import copy
import Constants as cons
import Netsim as ns
##### ----------------------------- Functions called in the main file ---------------------------------########
# Create the graph from geojson files
def load_network_data(DamageNodes,ExposureLines,NetworkFragility):
node_types_keys=NetworkFragility[cons.META][cons.TAXONOMIES]
G=nx.Graph()
EdgesFea=ExposureLines[cons.FEATURES]
NodesFea=DamageNodes[cons.FEATURES]
NodeBunch={}
#Initialize attributes for nodes and edges
NodeBunch={NodesFea[i][cons.PROPERTIES][cons.NODE_NAME]:{
cons.TAXONOMY:'',
cons.COORDINATES:NodesFea[i][cons.GEOMETRY][cons.COORDINATES],
cons.NODE_DAMAGE:0.0,
cons.NODE_DELTADAMAGE:0.0,
cons.NODE_POF:0.0
} for i in range(0,len(NodesFea))}
EdgeBunch=[(EdgesFea[i][cons.PROPERTIES][cons.FROM],EdgesFea[i][cons.PROPERTIES][cons.TO],{
cons.LINE_NAME:'',
cons.LENGTH:1.0,
cons.COORDINATES:EdgesFea[i][cons.GEOMETRY][cons.COORDINATES],
cons.VOLTAGE:1.0,
cons.WEIGHT:1.0,
cons.RESISTANCE:1.0,
cons.REACTANCE:1.0,
cons.LINE_DAMAGE:0.0,
cons.LINE_DELTADAMAGE:0.0,
cons.LOAD:1.0,
cons.CAPACITY:1.0
}) for i in range(0,len(EdgesFea))]
G.add_edges_from(EdgeBunch)
del EdgeBunch
nx.set_node_attributes(G,NodeBunch)
#set node attributes (only from 'properties' key)
for attr in NodesFea[0][cons.PROPERTIES].keys():
if attr==cons.NODE_NAME:
continue
NodeAttr={NodesFea[i][cons.PROPERTIES][cons.NODE_NAME]:{attr:NodesFea[i][cons.PROPERTIES][attr]} for i in range(0,len(NodesFea))}
nx.set_node_attributes(G,NodeAttr)
#set edge attributes (only from 'properties' key)
for attr in EdgesFea[0][cons.PROPERTIES].keys():
EdgeAttr={(EdgesFea[i][cons.PROPERTIES][cons.FROM],EdgesFea[i][cons.PROPERTIES][cons.TO]):{attr:EdgesFea[i][cons.PROPERTIES][attr]} for i in range(0,len(EdgesFea))}
nx.set_edge_attributes(G,EdgeAttr)
#update the weights; depending on available data
EdgeWeights={}
#try desired case: weightL*X, (L*sqrt(R^2+X^2), R<<X), R: resistance, X: reactance
try:
#EdgeWeights={(EdgesFea[i][cons.PROPERTIES][cons.FROM],EdgesFea[i][cons.PROPERTIES][cons.TO]):
#{cons.WEIGHT:float(EdgesFea[i][cons.PROPERTIES][cons.LENGTH])*
# np.sqrt((float(EdgesFea[i][cons.PROPERTIES][cons.REACTANCE])**2+
# float(EdgesFea[i][cons.PROPERTIES][cons.RESISTANCE])**2))} for i in range(0,len(EdgesFea))}
EdgeWeights={(EdgesFea[i][cons.PROPERTIES][cons.FROM],EdgesFea[i][cons.PROPERTIES][cons.TO]):
{cons.WEIGHT:float(EdgesFea[i][cons.PROPERTIES][cons.LENGTH])*
float(EdgesFea[i][cons.PROPERTIES][cons.REACTANCE])} for i in range(0,len(EdgesFea))}
except:
#otherwise, try just with the resistance: weight=L*R
try:
EdgeWeights={(EdgesFea[i][cons.PROPERTIES][cons.FROM],EdgesFea[i][cons.PROPERTIES][cons.TO]):
{cons.WEIGHT:float(EdgesFea[i][cons.PROPERTIES][cons.LENGTH])*
float(EdgesFea[i][cons.PROPERTIES][cons.RESISTANCE])} for i in range(0,len(EdgesFea))}
#otherwise, just assume weight~length/voltage
except:
EdgeWeights={(EdgesFea[i][cons.PROPERTIES][cons.FROM],EdgesFea[i][cons.PROPERTIES][cons.TO]):
{cons.WEIGHT:float(EdgesFea[i][cons.PROPERTIES][cons.LENGTH])/(float(EdgesFea[i][cons.PROPERTIES][cons.VOLTAGE]))} for i in range(0,len(EdgesFea))}
nx.set_edge_attributes(G,EdgeWeights)
# create a dictionary of node lists by taxonomy
Types=nx.get_node_attributes(G,cons.NODE_TYPE)
#node_types={typ:[nod for nod in G.nodes() if Types[nod]==typ] for typ in node_types_keys}
source=NetworkFragility[cons.META][cons.SOURCE]
consumer=NetworkFragility[cons.META][cons.CONSUMER]
s_nodes=[nod for nod in G.nodes() if Types[nod] in source]
c_nodes=[nod for nod in G.nodes() if Types[nod] in consumer]
return G,s_nodes,c_nodes
# evaluation of system loads
# load is computed as the number of shortest paths that pass through the component (node or edge)
def evaluate_system_loads(G,s_nodes,c_nodes):
ns.evaluate_system_loads(G,s_nodes,c_nodes)
#Assign component capacities
def assign_initial_capacities(G,alpha):
node_caps={no:{cons.CAPACITY:alpha*G.nodes[no][cons.LOAD]} for no in G.nodes()}
edge_caps={ed:{cons.CAPACITY:alpha*G.edges[ed][cons.LOAD]} for ed in G.edges()}
nx.set_node_attributes(G,node_caps)
nx.set_edge_attributes(G,edge_caps)
# MONTE CARLO SIMULATION
def run_Monte_Carlo_simulation(Graph,s_nodes0,c_nodes0,ExposureConsumerAreas,mcs):
#component_state_dha=[]
#component_state_idp=[]
#store samples of affected areas in a list
affected_areas=[]
max_iteration=5#max number of iterations in simulation of cascading effects
for i in range(0,mcs):
i_component_state_dha={}
#i_component_state_casc={}
#modify the networks in a copy of them
NetG=copy.deepcopy(Graph)
s_nodes=copy.deepcopy(s_nodes0)
c_nodes=copy.deepcopy(c_nodes0)
## Direct Hazard Action
# Simulate effects of hazard action on components
ns.direct_hazard_action(NetG)
# store coponent damage states
i_component_state_dha={cons.NODES:nx.get_node_attributes(NetG,cons.NODE_DAMAGE),cons.EDGES:nx.get_edge_attributes(NetG,cons.LINE_DAMAGE)}
#component_state_dha.append(i_component_state_dha)
# Update network description (weights and capacities)
ns.update_network(NetG,s_nodes,c_nodes)
# at least one component has significant damage (>10% damage)
if max(i_component_state_dha[cons.NODES].values())>cons.MIN_DAMAGE or max(i_component_state_dha[cons.EDGES].values())>cons.MIN_DAMAGE:
#i_component_state_casc=i_component_state_dha
## Damage Propagation
# if there are surviving supply and consumtion nodes
if len(s_nodes)>0 and len(c_nodes)>0:
# cascading effects
ns.simulate_cascading_effects(NetG,s_nodes,c_nodes,max_iteration)
#else:
# loss of functionality, i.e. no more flow in network
ns.update_network(NetG,s_nodes,c_nodes)
# no disconnection or cascading failures in this network
# Save Results
#component_state_idp.append(i_component_state_casc)
i_affected_areas=ns.set_state_consumers(ExposureConsumerAreas,NetG)
affected_areas.append(i_affected_areas)
print("MCS iteration: "+str(i))
return affected_areas
#Post Processing
def compute_output(SampleDamageAreas,ExposureConsumerAreas,nmcs):
SampleDamageNetwork=[int(np.sum([SampleDamageAreas[i][i_area]*ExposureConsumerAreas[cons.FEATURES][i_area][cons.PROPERTIES][cons.AREA_POPULATION]
for i_area in range(0,len(SampleDamageAreas[0]))])) for i in range(0,nmcs)]
for i_area in range(0,len(ExposureConsumerAreas[cons.FEATURES])):
est_apof=np.mean([SampleDamageAreas[i][i_area]for i in range(0,nmcs)])
ExposureConsumerAreas[cons.FEATURES][i_area][cons.PROPERTIES][cons.AREA_POF]=est_apof
return ExposureConsumerAreas,SampleDamageNetwork