-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathSensitivityMatrix_02.R
158 lines (113 loc) · 5.71 KB
/
SensitivityMatrix_02.R
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
#...............................................................................
# 1. - Initialize Session ####
#...............................................................................
cat("\014")
rm(list=ls())
# Installs libraries
library(tidyverse)
library(dplyr)
library(zoo)
library(lubridate)
library(epanetReader)
library(epanet2toolkit)
library(ggfortify)
library(ggthemes)
library(scales)
library(purrr)
library(visNetwork)
# Initialize params
params <- list(base_network = "base_dma_02",
new_network = "base_dma_w_leaks",
functs_name = "epanet_api_functions",
inlet_valves = "PRV_",
jt_to_analyze = "^JT_0[A-K]", # RegExp
pipe_to_analyze = "PS_", # RegExp
leak_rate = 0.01, # Percentage of the network with leaks
demad_factor = list( names =c( "wd_spring_summer",
"hw_spring_summer",
"wd_summer_break",
"hw_summer_break",
"wd_fall_winter",
"hw_fall_winter"),
factors = c( 0.92, 1.00, 1.09,
0.81, 0.66, 0.95)),
work_folders = list( dir_work = getwd(),
dir_report = file.path(getwd(),"reports"),
dir_data = file.path(getwd(),"data"),
dir_bin = file.path(getwd(),"reports"),
dir_func = file.path(getwd(),"func")))
params$f_names <- list( base_file_inp = file.path(params$work_folders$dir_data,
paste0(params$base_network,".inp")),
base_file_report = file.path(params$work_folders$dir_report,
paste0(params$base_network,".rpt")),
new_file_inp = file.path(params$work_folders$dir_data,
paste0(params$new_network,".inp")),
new_file_report = file.path(params$work_folders$dir_report,
paste0(params$new_network,".rpt")),
file_func = file.path(params$work_folders$dir_func,
paste0(params$functs_name,".R")))
# Load Functions Standard
source(params$f_names$file_func)
# load model
net_input_01 <- read.inp(params$f_names$base_file_inp)
# initialise th Juction for the leakage analyse
junctions <- net_input_01$Junctions %>%
filter(grepl(params$jt_to_analyze,ID))
# initialise data frames
leak_node <- tibble( "ID" = character(),
"FlowCoef" = double(),
"leakflow" = double())
# Generate Leack
emitters_ID <- junctions$ID[sample(1:nrow(junctions),1)]
# coef ~= y = 0.1436*(l/s) - 0.0026
# Leack(l/s) = 2.0 then coef = 0.285
coefficient <- 0.2846
net_input_01$Emitters <- data.frame(ID = emitters_ID, FlowCoef = coefficient)
write.inp(net_input_01, params$f_names$new_file_inp)
rm(net_input_01)
#...............................................................................
# 3. Running a Full Simulation ####
# The function ENepanet() runs a full simulation and
# writes the results to a file.
#...............................................................................
ENepanet(params$f_names$base_file_inp,
params$f_names$base_file_report)
ENepanet(params$f_names$new_file_inp,
params$f_names$new_file_report)
net_input_01 <- read.inp(params$f_names$new_file_inp)
report_base <- read.rpt(params$f_names$base_file_report)
report_leack <- read.rpt(params$f_names$new_file_report)
# 4.1.- the average inflow f(t) was calculated at each hour t
base_inletflow <- inlet_flows(report_base, params$inlet_valves, group = FALSE)
leak_inletflow <- inlet_flows(report_leack, params$inlet_valves, group = FALSE)
inletflow <- full_join(base_inletflow, leak_inletflow,
by = "timeInSeconds") %>%
select(timeInSeconds, inflow.x, inflow.y) %>%
mutate(leakflow = inflow.y - inflow.x)
inletflow <- median(inletflow$leakflow)
#-----
leak_node <- add_row(leak_node, "ID" = emitters_ID,
"FlowCoef" = coefficient,
"leakflow" = inletflow)
#-------------------------
rm(base_inletflow,leak_inletflow)
# residual vector
rep01 <- eval_nodes (report_base, node_type = "",
id_nodes = params$jt_to_analyze,
group = FALSE)
rep02 <- eval_nodes (report_leack, node_type = "",
id_nodes = params$jt_to_analyze,
group = FALSE)
residual_vector <- full_join(rep01, rep02, by = c("timeInSeconds","ID")) %>%
mutate(D_Pressure = Pressure.y - Pressure.x) %>%
select(timeInSeconds,ID, D_Pressure) %>%
group_by(ID) %>%
summarise( rv_min = min(D_Pressure))
residual_matrix <- 1
rm(rep01,rep02)
# sensitivity vector
sensitivity_vector <- residual_vector %>%
mutate(sensitivity = rv_min/inletflow)
# detection capability matrix Mdc Where :
# - the rows represent the pontential sensor locations and
# - the columns represent the potential leak location