From b83b895e2b8a2b05a9021c1b17ffebc2f311d3d7 Mon Sep 17 00:00:00 2001 From: Jiang Date: Fri, 6 Mar 2026 15:27:59 +0800 Subject: [PATCH] =?UTF-8?q?=E6=96=B0=E5=A2=9E=E7=88=86=E7=AE=A1=E4=BD=8D?= =?UTF-8?q?=E7=BD=AE=E6=A3=80=E6=B5=8B=E6=A8=A1=E5=9D=97=E5=8F=8A=E7=9B=B8?= =?UTF-8?q?=E5=85=B3API=E6=8E=A5=E5=8F=A3?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- app/algorithms/burst_location/__init__.py | 3 + .../burst_location/burst_location.py | 255 ++++++ .../burst_location/burst_locator.py | 609 +++++++++++++ .../burst_location/leak_simulator.py | 489 +++++++++++ .../burst_location/network_model.py | 169 ++++ .../burst_location/network_partitioner.py | 375 ++++++++ .../burst_location/noise_generator.py | 234 +++++ .../burst_location/similarity_calculator.py | 830 ++++++++++++++++++ app/api/v1/endpoints/burst_location.py | 33 + app/api/v1/router.py | 4 + app/services/burst_location.py | 83 ++ 11 files changed, 3084 insertions(+) create mode 100644 app/algorithms/burst_location/__init__.py create mode 100644 app/algorithms/burst_location/burst_location.py create mode 100644 app/algorithms/burst_location/burst_locator.py create mode 100644 app/algorithms/burst_location/leak_simulator.py create mode 100644 app/algorithms/burst_location/network_model.py create mode 100644 app/algorithms/burst_location/network_partitioner.py create mode 100644 app/algorithms/burst_location/noise_generator.py create mode 100644 app/algorithms/burst_location/similarity_calculator.py create mode 100644 app/api/v1/endpoints/burst_location.py create mode 100644 app/services/burst_location.py diff --git a/app/algorithms/burst_location/__init__.py b/app/algorithms/burst_location/__init__.py new file mode 100644 index 0000000..c7ec7a1 --- /dev/null +++ b/app/algorithms/burst_location/__init__.py @@ -0,0 +1,3 @@ +from .burst_location import run_burst_location + +__all__ = ["run_burst_location"] diff --git a/app/algorithms/burst_location/burst_location.py b/app/algorithms/burst_location/burst_location.py new file mode 100644 index 0000000..0f24c02 --- /dev/null +++ b/app/algorithms/burst_location/burst_location.py @@ -0,0 +1,255 @@ +"""爆管定位部署入口（基于外部 SCADA 实测数据）。""" + +import argparse +import json +from pathlib import Path +from typing import Iterable + +import pandas as pd + +from .burst_locator import ( + DN_search_multi_simple_add_flow_count_new, +) +from .network_model import ( + _build_node_pipe_maps, + cal_node_coordinate, + construct_graph, + load_inp, + read_inf_inp, + read_inf_inp_other, +) + + +def _read_id_list_json(path): + if path is None: + return None + data = json.loads(Path(path).read_text(encoding="utf-8")) + if isinstance(data, list): + return [str(item) for item in data] + if isinstance(data, dict): + if "ids" in data and isinstance(data["ids"], list): + return [str(item) for item in data["ids"]] + raise ValueError(f"ID JSON must be list or dict with key 'ids': {path}") + raise ValueError(f"Unsupported ID JSON format: {path}") + + +def _read_series_csv(path): + if path is None: + return None + df = pd.read_csv(path) + if df.shape[1] < 2: + raise ValueError(f"CSV must contain at least two columns (id,value): {path}") + if {"id", "value"}.issubset(df.columns): + id_col, value_col = "id", "value" + else: + id_col, value_col = df.columns[0], df.columns[1] + series = pd.Series( + df[value_col].values, index=df[id_col].astype(str).values, dtype=float + ) + return series + + +def _align_scada_series( + series: pd.Series, ids: Iterable[str], series_name: str +) -> pd.Series: + ids = [str(item) for item in ids] + aligned = series.copy() + aligned.index = aligned.index.map(str) + missing_ids = [item for item in ids if item not in aligned.index] + if missing_ids: + preview = ", ".join(missing_ids[:10]) + raise ValueError(f"{series_name} missing IDs: {preview}") + aligned = pd.to_numeric(aligned.loc[ids], errors="coerce") + invalid_ids = aligned[aligned.isna()].index.tolist() + if invalid_ids: + preview = ", ".join(invalid_ids[:10]) + raise ValueError( + f"{series_name} contains non-numeric values for IDs: {preview}" + ) + return aligned + + +def run_burst_location( + wn_inp_path, + pressure_scada_ids, + burst_pressure, + normal_pressure, + burst_leakage, + flow_scada_ids=None, + burst_flow=None, + normal_flow=None, + min_dpressure=2.0, + basic_pressure=10.0, +): + if pressure_scada_ids is None or len(pressure_scada_ids) == 0: + raise ValueError("pressure_scada_ids cannot be empty.") + if burst_pressure is None or normal_pressure is None: + raise ValueError("burst_pressure and normal_pressure are required.") + + has_any_flow = any( + value is not None for value in [flow_scada_ids, burst_flow, normal_flow] + ) + has_all_flow = all( + value is not None for value in [flow_scada_ids, burst_flow, normal_flow] + ) + if has_any_flow and not has_all_flow: + raise ValueError( + "flow_scada_ids, burst_flow, and normal_flow must be provided together." + ) + + inp_path = Path(wn_inp_path) + wn = load_inp( + inp_name=inp_path.name, + inp_location=str(inp_path.parent) + "/", + inp_time=0, + driven_mode="PDD", + require_p=float(basic_pressure), + minimum_p=0.0, + ) + + all_node, _, node_coordinates, candidate_pipe, _, _, pipe_length, _ = read_inf_inp( + wn + ) + _, pipe_start_node_all, pipe_end_node_all = read_inf_inp_other(wn) + node_x, node_y = cal_node_coordinate(all_node, node_coordinates) + G0 = construct_graph(wn) + node_pipe_dic, couple_node_length = _build_node_pipe_maps( + all_node, candidate_pipe, pipe_start_node_all, pipe_end_node_all, pipe_length + ) + all_node_series = pd.Series(range(len(all_node)), index=all_node) + + pressure_ids = [str(item) for item in pressure_scada_ids] + normal_pressure_aligned = _align_scada_series( + normal_pressure, pressure_ids, "normal_pressure" + ) + burst_pressure_aligned = _align_scada_series( + burst_pressure, pressure_ids, "burst_pressure" + ) + pressure_normal = normal_pressure_aligned.to_frame().T + pressure_monitor = burst_pressure_aligned.to_frame().T + pressure_predict = pressure_normal.copy() + timestep_list = list(pressure_normal.index) + + if has_all_flow: + flow_ids = [str(item) for item in flow_scada_ids] + if len(flow_ids) == 0: + raise ValueError( + "flow_scada_ids cannot be empty when flow data is provided." + ) + normal_flow_aligned = _align_scada_series(normal_flow, flow_ids, "normal_flow") + burst_flow_aligned = _align_scada_series(burst_flow, flow_ids, "burst_flow") + flow_normal = normal_flow_aligned.to_frame().T + flow_monitor = burst_flow_aligned.to_frame().T + flow_predict = flow_normal.copy() + similarity_mode = "CDF" + max_flow = flow_normal.iloc[0, :].abs() + else: + flow_normal = pd.DataFrame(index=timestep_list) + flow_monitor = pd.DataFrame(index=timestep_list) + flow_predict = pd.DataFrame(index=timestep_list) + similarity_mode = "CAD_new_gy" + max_flow = pd.Series(dtype=float) + + located_pipe, elapsed_seconds, simulation_times, _, similarity_series = ( + DN_search_multi_simple_add_flow_count_new( + wn=wn, + G0=G0, + all_node=all_node, + node_x=node_x, + node_y=node_y, + pipe_start_node_all=pipe_start_node_all, + pipe_end_node_all=pipe_end_node_all, + couple_node_length=couple_node_length, + node_pipe_dic=node_pipe_dic, + all_node_series=all_node_series, + top_group_ratio=0.3, + top_pipe_num_max=80, + top_pipe_num_min=10, + candidate_pipe_input_initial=candidate_pipe, + similarity_mode=similarity_mode, + pressure_monitor=pressure_monitor, + pressure_predict=pressure_predict, + pressure_normal=pressure_normal, + pressure_leak_all=None, + flow_monitor=flow_monitor, + flow_predict=flow_predict, + flow_normal=flow_normal, + flow_leak_all=None, + timestep_list=timestep_list, + max_flow=max_flow, + group_basic_num=30, + Top_sensor_num=min(5, len(pressure_ids)), + if_gy=0, + pressure_threshold=float(min_dpressure), + leak_mag=float(burst_leakage), + ) + ) + + return { + "located_pipe": located_pipe, + "burst_leakage": float(burst_leakage), + "elapsed_seconds": elapsed_seconds, + "simulation_times": int(simulation_times), + "top_candidates": list(similarity_series.index[:10]), + "similarity_mode": similarity_mode, + } + + +def _parse_args(): + parser = argparse.ArgumentParser(description="爆管定位主函数入口") + parser.add_argument("--wn-inp", required=True, help="EPANET inp 文件路径") + parser.add_argument( + "--pressure-ids-json", required=True, help="压力SCADA ID列表 JSON 文件" + ) + parser.add_argument( + "--flow-ids-json", default=None, help="（可选）流量SCADA ID列表 JSON 文件" + ) + parser.add_argument( + "--burst-pressure-csv", required=True, help="爆管时压力 CSV（id,value）" + ) + parser.add_argument( + "--normal-pressure-csv", required=True, help="正常时压力 CSV（id,value）" + ) + parser.add_argument( + "--burst-flow-csv", default=None, help="（可选）爆管时流量 CSV（id,value）" + ) + parser.add_argument( + "--normal-flow-csv", default=None, help="（可选）正常时流量 CSV（id,value）" + ) + parser.add_argument( + "--burst-leakage", type=float, required=True, help="爆管漏损流量" + ) + parser.add_argument( + "--min-dpressure", + type=float, + default=2.0, + help="（可选）最小压降阈值，默认 2.0", + ) + parser.add_argument( + "--basic-pressure", + type=float, + default=10.0, + help="（可选）基础服务压力，默认 10.0", + ) + return parser.parse_args() + + +def main(): + args = _parse_args() + result = run_burst_location( + wn_inp_path=args.wn_inp, + pressure_scada_ids=_read_id_list_json(args.pressure_ids_json), + burst_pressure=_read_series_csv(args.burst_pressure_csv), + normal_pressure=_read_series_csv(args.normal_pressure_csv), + burst_leakage=args.burst_leakage, + flow_scada_ids=_read_id_list_json(args.flow_ids_json), + burst_flow=_read_series_csv(args.burst_flow_csv), + normal_flow=_read_series_csv(args.normal_flow_csv), + min_dpressure=args.min_dpressure, + basic_pressure=args.basic_pressure, + ) + print(json.dumps(result, ensure_ascii=False)) + + +if __name__ == "__main__": + main() diff --git a/app/algorithms/burst_location/burst_locator.py b/app/algorithms/burst_location/burst_locator.py new file mode 100644 index 0000000..77c4af9 --- /dev/null +++ b/app/algorithms/burst_location/burst_locator.py @@ -0,0 +1,609 @@ +"""爆管定位主模块。""" + +import copy +import math +import sys +from datetime import datetime + +import networkx as nx +import numpy as np +import pandas as pd + +from .leak_simulator import cal_signature_pipe_multi_pf +from .network_partitioner import cal_group_num, metis_grouping_pipe_weight +from .similarity_calculator import ( + adjust_ratio, + cal_similarity_all_multi_new_sq_improve_double_lzr, + decode_mode, + extra_judge, + update_similarity, +) + + +def _ensure_signatures_for_centers( + wn, + center_list, # 本轮要用到的中心（list[str]） + pressure_leak_all, + flow_leak_all, # 全量缓存（可为空 DF） + timestep_list, # 你现有的时序列表 + pressure_monitor, + flow_monitor, # 用来推断传感器列名 + leak_mag, # 泄漏量，比如 400/3600 +): + """ + 只为缺失的中心补算 SLF（调用你现有的 cal_signature_pipe_multi_pf）， + 并把补算结果并回缓存。返回： + pressure_leak_subset, flow_leak_subset, pressure_leak_all_new, flow_leak_all_new + 其中 subset 只包含 center_list 的行（顺序与 center_list 保持一致）。 + """ + + # 1) 推断传感器列名（与现有数据保持一致） + sensor_name_all = list(pressure_monitor.columns) + sensor_f_name_all = ( + list(flow_monitor.columns) + if (flow_monitor is not None and hasattr(flow_monitor, "columns")) + else [] + ) + + # 2) 取出缓存里已经有的中心（考虑 MultiIndex 的第 0 层为 pipe） + def _existing_pipes(df): + if df is None or len(df) == 0: + return set() + idx = df.index + if isinstance(idx, pd.MultiIndex): + return set(idx.get_level_values(0)) + else: + return set(idx) + + exist_p = _existing_pipes(pressure_leak_all) + need = [p for p in center_list if p not in exist_p] + + # 3) 若有缺失中心，仅为这些中心补算一次 + if len(need) > 0: + p_new, _ = cal_signature_pipe_multi_pf( + wn, leak_mag, need, timestep_list, sensor_name_all + ) + # 初始化空缓存时，做一次“同构化” + if pressure_leak_all is None or len(pressure_leak_all) == 0: + pressure_leak_all = p_new + else: + pressure_leak_all = pd.concat([pressure_leak_all, p_new], axis=0) + # if (flow_leak_all is None or len(flow_leak_all) == 0) and f_new is not None: + # flow_leak_all = f_new + # elif f_new is not None: + # flow_leak_all = pd.concat([flow_leak_all, f_new], axis=0) + + # 去重（如果既有缓存里不小心有重复中心） + if isinstance(pressure_leak_all.index, pd.MultiIndex): + pressure_leak_all = pressure_leak_all[ + ~pressure_leak_all.index.duplicated(keep="last") + ] + if flow_leak_all is not None and len(flow_leak_all) > 0: + flow_leak_all = flow_leak_all[ + ~flow_leak_all.index.duplicated(keep="last") + ] + else: + pressure_leak_all = pressure_leak_all[ + ~pressure_leak_all.index.duplicated(keep="last") + ] + if flow_leak_all is not None and len(flow_leak_all) > 0: + flow_leak_all = flow_leak_all[ + ~flow_leak_all.index.duplicated(keep="last") + ] + + # 4) 从更新后的缓存里，取出这轮需要的中心子集（顺序与 center_list 一致） + if isinstance(pressure_leak_all.index, pd.MultiIndex): + pressure_subset = pressure_leak_all.loc[center_list] + flow_subset = ( + flow_leak_all.loc[center_list] + if (flow_leak_all is not None and len(flow_leak_all) > 0) + else None + ) + else: + pressure_subset = pressure_leak_all.loc[center_list, :] + flow_subset = ( + flow_leak_all.loc[center_list, :] + if (flow_leak_all is not None and len(flow_leak_all) > 0) + else None + ) + + return pressure_subset, flow_subset, pressure_leak_all, flow_leak_all + + +def area_output_num_ki_improve( + candidate_center, + candidate_group, + similarity, + new_all_node, + top_group_ratio, + top_pipe_num_max, + top_pipe_num_min, + cut_ratio, +): + final_area = [] + final_center = [] + all_node_iter = [] + if similarity.index.is_unique == False: + total_center_num = len(list(set(similarity.index))) + else: + total_center_num = len(similarity.index) + + next_group_num = min( + total_center_num, math.ceil(total_center_num / cut_ratio * top_group_ratio) + ) + + for i in range(next_group_num): + top_center = similarity.index[i] + top_center_index = find_list_repeat(candidate_center, top_center) + for j in range(len(top_center_index)): + final_area = final_area + candidate_group[top_center_index[j]] + all_node_iter = all_node_iter + list(new_all_node[top_center_index[j]]) + final_center.append(top_center) + final_area = list(set(final_area)) + + if len(final_area) > top_pipe_num_max: + if_end = 0 + elif len(final_area) > top_pipe_num_min: + if_end = 1 + elif total_center_num == next_group_num: + if_end = 1 + else: + if_end = 1 + for i in np.arange(next_group_num, total_center_num, 1): + + before_list = copy.deepcopy(final_area) + top_center = similarity.index[i] + top_center_index = candidate_center.index(top_center) + temp_group = final_area + candidate_group[top_center_index] + temp_area = list(set(temp_group)) + + if len(temp_area) < top_pipe_num_min: + + final_center.append(top_center) + all_node_iter = all_node_iter + list(new_all_node[top_center_index]) + final_area = temp_area + + elif len(temp_area) < top_pipe_num_max: + final_center.append(top_center) + all_node_iter = all_node_iter + list(new_all_node[top_center_index]) + final_area = temp_area + break + else: + a = len(temp_area) - top_pipe_num_max + b = top_pipe_num_min - len(before_list) + + if a >= b: + final_area = before_list + else: + final_center.append(top_center) + all_node_iter = all_node_iter + list(new_all_node[top_center_index]) + final_area = temp_area + break + + final_center = list(set(final_center)) + all_node_iter = list(set(all_node_iter)) + return final_area, final_center, all_node_iter, if_end + + +def find_list_repeat(candidate_center, target): + repeated_list = [] + for index, nums in enumerate(candidate_center): + if nums == target: + repeated_list.append(index) + return repeated_list + + +def cal_DtoTop1( + G0, pipe_leak, located_pipe, pipe_start_node_all, pipe_end_node_all, pipe_length +): + if pipe_leak == located_pipe: + result_DtoTop1 = 0 + result_DtoTop1_num = 0 + else: + pipe_leak_start_node = pipe_start_node_all[pipe_leak] + pipe_leak_end_node = pipe_end_node_all[pipe_leak] + located_pipe_start_node = pipe_start_node_all[located_pipe] + located_pipe_end_node = pipe_end_node_all[located_pipe] + DtoTop1_series = pd.Series(dtype=object) + DtoTop1_num_series = pd.Series(dtype=object) + DtoTop1_series["ss"] = nx.shortest_path_length( + G0, pipe_leak_start_node, located_pipe_start_node, weight="weight" + ) + DtoTop1_series["se"] = nx.shortest_path_length( + G0, pipe_leak_start_node, located_pipe_end_node, weight="weight" + ) + DtoTop1_series["es"] = nx.shortest_path_length( + G0, pipe_leak_end_node, located_pipe_start_node, weight="weight" + ) + DtoTop1_series["ee"] = nx.shortest_path_length( + G0, pipe_leak_end_node, located_pipe_end_node, weight="weight" + ) + DtoTop1_num_series["ss"] = nx.shortest_path_length( + G0, pipe_leak_start_node, located_pipe_start_node + ) + DtoTop1_num_series["se"] = nx.shortest_path_length( + G0, pipe_leak_start_node, located_pipe_end_node + ) + DtoTop1_num_series["es"] = nx.shortest_path_length( + G0, pipe_leak_end_node, located_pipe_start_node + ) + DtoTop1_num_series["ee"] = nx.shortest_path_length( + G0, pipe_leak_end_node, located_pipe_end_node + ) + + if DtoTop1_num_series.min() == 0: + result_DtoTop1_num = 1 + result_DtoTop1 = DtoTop1_series.max() / 2 + else: + result_DtoTop1_num = DtoTop1_num_series.min() + 1 + DtoTop1_type = DtoTop1_series.argmin() + result_DtoTop1 = ( + DtoTop1_series[DtoTop1_type] + + (pipe_length[pipe_leak] + pipe_length[located_pipe]) / 2 + ) + return result_DtoTop1, result_DtoTop1_num + + +def cal_RR(located_pipe, similarity_sp): + if located_pipe in similarity_sp.index: + rank = similarity_sp.index.get_loc(located_pipe) + RR = rank / len(similarity_sp.index) + else: + RR = 1.1 + return RR + + +def cal_cover(similarity, leak_pipe): + if leak_pipe in list(similarity.index): + cover = 1 + else: + cover = 0 + return cover + + +def cal_SD(located_pipe, real_pipe, pipe_x, pipe_y): + dx = pipe_x[located_pipe] - pipe_x[real_pipe] + dy = pipe_y[located_pipe] - pipe_y[real_pipe] + SD = math.sqrt(dx * dx + dy * dy) + return SD + + +def DN_search_multi_simple_add_flow_count_new( + wn, + G0, + all_node, + node_x, + node_y, + pipe_start_node_all, + pipe_end_node_all, + couple_node_length, + node_pipe_dic, + all_node_series, + top_group_ratio, + top_pipe_num_max, + top_pipe_num_min, + candidate_pipe_input_initial, + similarity_mode, + pressure_monitor, + pressure_predict, + pressure_normal, + pressure_leak_all, + flow_monitor, + flow_predict, + flow_normal, + flow_leak_all, + timestep_list, + max_flow, + group_basic_num, + Top_sensor_num, + if_gy, + pressure_threshold, + leak_mag=400 / 3600, +): + iter_count = 0 + all_node_iter = copy.deepcopy(all_node) + candidate_pipe_input = copy.deepcopy(candidate_pipe_input_initial) # 可能漏损管段 + t1 = datetime.now() + if_flow, if_only_cos, if_only_flow = decode_mode(similarity_mode) # 定位方法 + # threshold + if if_only_flow == 1: + dpressure = (flow_predict - flow_monitor).mean() + dpressure = dpressure.abs() + effective_sensor = list(dpressure.index) + + else: + dpressure = (pressure_predict - pressure_monitor).mean() + dpressure = dpressure.abs() + dpressure = dpressure[dpressure > pressure_threshold] + effective_sensor = list(dpressure.index) + simulation_times = 0 # 模拟次数 + if len(dpressure) > 0: + + cos_h = 0 + dis_h = 0 + dis_f_h = 0 + if_compalsive = 0 + record_center_dataset = [] + # iter + while 1: + final_area = [] + final_center = [] + group_num = cal_group_num(candidate_pipe_input, group_basic_num) + + # group 分组，得出候选漏损中心 + candidate_center_list, candidate_group_list, new_all_node = ( + metis_grouping_pipe_weight( + G0, + wn, + all_node_iter, + candidate_pipe_input, + group_num, + node_x, + node_y, + pipe_start_node_all, + pipe_end_node_all, + node_pipe_dic, + all_node_series, + couple_node_length, + ) + ) + simulation_times = simulation_times + len(candidate_center_list) + # pick_pressure_leak + # pressure_leak = pressure_leak_all.loc[candidate_center_list].loc[:, :] + # flow_leak = flow_leak_all.loc[candidate_center_list].loc[:, :] + # —— 新增泄漏量（保持你现在的一致，或从外部传入）—— + # —— 只为缺失中心补算，然后取本轮需要的中心子集 —— + pressure_leak, flow_leak, pressure_leak_all, flow_leak_all = ( + _ensure_signatures_for_centers( + wn=wn, + center_list=candidate_center_list, + pressure_leak_all=pressure_leak_all, + flow_leak_all=flow_leak_all, + timestep_list=timestep_list, + pressure_monitor=pressure_monitor, + flow_monitor=flow_monitor, + leak_mag=leak_mag, + ) + ) + + # pressure_leak_f= pressure_leak.swaplevel() + + # -------------------------------------------------------- + + add_center = [] + leak_center_dict = dict() + for i in range(len(candidate_center_list)): + houxuan_center = [] + for each_center in record_center_dataset: + if ( + each_center in candidate_group_list[i] + and each_center != candidate_center_list[i] + ): + houxuan_center.append(each_center) + add_center = add_center + houxuan_center + houxuan_center.append(candidate_center_list[i]) + leak_center_dict[candidate_center_list[i]] = houxuan_center + for each_center in candidate_center_list: + if each_center not in record_center_dataset: + record_center_dataset.append(each_center) + + # -------------------------------------------------------- + # -------------------------------------------------------- + # if len(add_center) > 0: + # s3 = pressure_leak_all.loc[add_center] + # pressure_leak = pd.concat([pressure_leak, s3]) + # s4 = flow_leak_all.loc[add_center] + # flow_leak = pd.concat([flow_leak, s4]) + # -------------------------------------------------------- + # 只为 add_center 里还没算过的中心补算，并与本轮中心合并 + if len(add_center) > 0: + pressure_add, flow_add, pressure_leak_all, flow_leak_all = ( + _ensure_signatures_for_centers( + wn=wn, + center_list=add_center, + pressure_leak_all=pressure_leak_all, + flow_leak_all=flow_leak_all, + timestep_list=timestep_list, + pressure_monitor=pressure_monitor, + flow_monitor=flow_monitor, + leak_mag=leak_mag, # 与上面一致 + ) + ) + pressure_leak = pd.concat([pressure_leak, pressure_add], axis=0) + if (flow_leak is not None) and (flow_add is not None): + flow_leak = pd.concat([flow_leak, flow_add], axis=0) + # -------------------------------------------------------- + # + if len(candidate_pipe_input) < 1.2 * top_pipe_num_max / top_group_ratio: + if_compalsive = 1 + cos_h, dis_h, dis_f_h = adjust_ratio(similarity_mode, cos_h, dis_h, dis_f_h) + candidate_center_list_sup = candidate_center_list + add_center + similarity, cos_h, dis_h, dis_f_h, break_flag = ( + cal_similarity_all_multi_new_sq_improve_double_lzr( + candidate_center_list_sup, + similarity_mode, + pressure_leak, + pressure_monitor, + pressure_predict, + pressure_normal, + if_flow, + if_only_cos, + if_only_flow, + flow_leak, + flow_monitor, + flow_predict, + flow_normal, + timestep_list, + Top_sensor_num, + if_gy, + effective_sensor, + cos_h, + dis_h, + dis_f_h, + if_compalsive, + max_flow, + ) + ) + if break_flag == 1: + break + + new_similarity = update_similarity( + candidate_center_list, similarity, leak_center_dict + ) + if len(candidate_pipe_input) > top_pipe_num_max / top_group_ratio: + cut_ratio, new_similarity = extra_judge(new_similarity) + else: + cut_ratio = 1 + final_area_t, final_center_t, all_node_new_1, if_end = ( + area_output_num_ki_improve( + candidate_center_list, + candidate_group_list, + new_similarity, + new_all_node, + top_group_ratio, + top_pipe_num_max, + top_pipe_num_min, + cut_ratio, + ) + ) + final_area = final_area + final_area_t + final_center = final_center + final_center_t + + final_area = list(set(final_area)) + final_center = list(set(final_center)) + if if_end == 1: + break + elif len(candidate_pipe_input) == len(final_area): + break + else: + candidate_pipe_input = final_area + all_node_iter = all_node_new_1 + iter_count += 1 + sys.stdout.write( + "\r" + + "已经完成" + + str(iter_count) + + "次迭代计算" + + "候选节点" + + str(len(final_area)) + + "个" + ) + # if break_flag == 0: + # final_area_pipe = copy.deepcopy(final_area) + # simulation_times = simulation_times + len(final_area) + # pressure_leak_sp = pressure_leak_all.loc[final_area_pipe].loc[:, :] + # flow_leak_sp = flow_leak_all.loc[final_area_pipe].loc[:, :] + # similarity_sp, cos_h, dis_h, dis_f_h, break_flag = cal_similarity_all_multi_new_sq_improve_double_lzr( + # final_area_pipe, similarity_mode, pressure_leak_sp, + # pressure_monitor, pressure_predict, pressure_normal, if_flow, + # if_only_cos, if_only_flow, + # flow_leak_sp, flow_monitor, flow_predict, flow_normal, + # timestep_list, Top_sensor_num, if_gy, effective_sensor, cos_h, dis_h, dis_f_h, if_compalsive, max_flow) + if break_flag == 0: + final_area_pipe = list(final_area) # 确保是 list + + # 只为还没算过的管段补齐 SLF（按需计算） + pressure_leak_sp, flow_leak_sp, pressure_leak_all, flow_leak_all = ( + _ensure_signatures_for_centers( + wn=wn, + center_list=final_area_pipe, # 这次要用的“最终区域里的所有管段” + pressure_leak_all=pressure_leak_all, # 累积缓存（会被更新） + flow_leak_all=flow_leak_all, + timestep_list=timestep_list, + pressure_monitor=pressure_monitor, + flow_monitor=flow_monitor, + leak_mag=leak_mag, + ) + ) + + # 如果你要精确统计模拟次数，这里可以加上“本次新补的数量”， + # 做法：让 _ensure_signatures_for_centers 额外返回 need_cnt，再 simulation_times += need_cnt + + similarity_sp, cos_h, dis_h, dis_f_h, break_flag = ( + cal_similarity_all_multi_new_sq_improve_double_lzr( + final_area_pipe, + similarity_mode, + pressure_leak_sp, + pressure_monitor, + pressure_predict, + pressure_normal, + if_flow, + if_only_cos, + if_only_flow, + flow_leak_sp, + flow_monitor, + flow_predict, + flow_normal, + timestep_list, + Top_sensor_num, + if_gy, + effective_sensor, + cos_h, + dis_h, + dis_f_h, + if_compalsive, + max_flow, + ) + ) + + else: + dpressure = (pressure_predict - pressure_monitor).mean() + dpressure = dpressure.abs() + simulation_times = simulation_times + len(dpressure.index) + similarity_sp = pd.Series(dtype=object) + for each_node in dpressure.index: + pipe = node_pipe_dic[each_node][0] + similarity_sp.loc[pipe] = dpressure.loc[each_node] + similarity_sp = similarity_sp.sort_values(ascending=False) + t2 = datetime.now() + final_area_pipe = [] + + sys.stdout.write( + "\r" + + "已经完成" + + str(iter_count + 1) + + "次迭代计算" + + "候选节点" + + str(len(final_area_pipe)) + + "个" + ) + t2 = datetime.now() + dt = (t2 - t1).seconds + else: + dpressure = (pressure_predict - pressure_monitor).mean() + dpressure = dpressure.abs() + + similarity_sp = pd.Series(dtype=object) + for each_node in dpressure.index: + pipe = node_pipe_dic[each_node][0] + similarity_sp.loc[pipe] = dpressure.loc[each_node] + similarity_sp = similarity_sp.sort_values(ascending=False) + t2 = datetime.now() + dt = (t2 - t1).seconds + return similarity_sp.index[0], dt, simulation_times, wn, similarity_sp + + +class BurstLocator: + @staticmethod + def DN_search_multi_simple_add_flow_count_new(*args, **kwargs): + return DN_search_multi_simple_add_flow_count_new(*args, **kwargs) + + @staticmethod + def area_output_num_ki_improve(*args, **kwargs): + return area_output_num_ki_improve(*args, **kwargs) + + @staticmethod + def cal_DtoTop1(*args, **kwargs): + return cal_DtoTop1(*args, **kwargs) + + @staticmethod + def cal_RR(*args, **kwargs): + return cal_RR(*args, **kwargs) + + @staticmethod + def cal_cover(*args, **kwargs): + return cal_cover(*args, **kwargs) + + @staticmethod + def cal_SD(*args, **kwargs): + return cal_SD(*args, **kwargs) diff --git a/app/algorithms/burst_location/leak_simulator.py b/app/algorithms/burst_location/leak_simulator.py new file mode 100644 index 0000000..b7616c6 --- /dev/null +++ b/app/algorithms/burst_location/leak_simulator.py @@ -0,0 +1,489 @@ +"""漏损模拟模块。""" + +import math +import sys + +import pandas as pd +import wntr + +_PIPE2LEAKNODE = None + + +def simple_add_leak(wn, leak_mag, leak_pipe): + whole_inf = dict() + leak_pipe_self = wn.get_link(leak_pipe) + pipe_diameter = leak_pipe_self.diameter + pipe_length = leak_pipe_self.length + pipe_roughness = leak_pipe_self.roughness + pipe_minor_loss = leak_pipe_self.minor_loss + # pipe_status = leak_pipe_self.status + # pipe_check_valve = leak_pipe_self.check_valve + pipe_start_node = leak_pipe_self.start_node_name + pipe_end_node = leak_pipe_self.end_node_name + + # close the pipe + # leak_pipe_self.status = 'Closed' + wn.remove_link(leak_pipe) + # add the pipe + add_pipe1 = leak_pipe + "A" + add_pipe2 = leak_pipe + "B" + add_node = leak_pipe + "_" + + start_n = wn.get_node(pipe_start_node) + end_n = wn.get_node(pipe_end_node) + if start_n.node_type == "Reservoir": + end_n_elevation = end_n.elevation + start_n_elevation = end_n_elevation + elif end_n.node_type == "Reservoir": + start_n_elevation = start_n.elevation + end_n_elevation = start_n_elevation + else: + end_n_elevation = end_n.elevation + start_n_elevation = start_n.elevation + elevation_self = (start_n_elevation + end_n_elevation) / 2 + coordinates_self = ( + (start_n.coordinates[0] + end_n.coordinates[0]) / 2, + (start_n.coordinates[1] + end_n.coordinates[1]), + ) + + wn.add_junction( + add_node, base_demand=0, elevation=elevation_self, coordinates=coordinates_self + ) + leak_node = wn.get_node(add_node) + wn.add_pipe( + add_pipe1, + start_node_name=pipe_start_node, + end_node_name=add_node, + length=pipe_length / 2, + diameter=pipe_diameter, + roughness=pipe_roughness, + minor_loss=pipe_minor_loss, + ) + wn.add_pipe( + add_pipe2, + start_node_name=pipe_end_node, + end_node_name=add_node, + length=pipe_length / 2, + diameter=pipe_diameter, + roughness=pipe_roughness, + minor_loss=pipe_minor_loss, + ) + + # simulation + leak_node.add_demand(base=leak_mag, pattern_name="add_leak") + + whole_inf["leak_node_name"] = add_node + whole_inf["add_pipe1"] = add_pipe1 + whole_inf["add_pipe2"] = add_pipe2 + whole_inf["leak_pipe"] = leak_pipe + whole_inf["pipe_start_node"] = pipe_start_node + whole_inf["pipe_end_node"] = pipe_end_node + whole_inf["pipe_length"] = pipe_length + whole_inf["pipe_diameter"] = pipe_diameter + whole_inf["pipe_roughness"] = pipe_roughness + whole_inf["pipe_minor_loss"] = pipe_diameter + return wn, whole_inf, add_pipe1 + + +def simple_recover_wn(wn, whole_inf): + leak_node = wn.get_node(whole_inf["leak_node_name"]) + + del leak_node.demand_timeseries_list[-1] + # update + wn.remove_link(whole_inf["add_pipe1"]) + wn.remove_link(whole_inf["add_pipe2"]) + wn.remove_node(whole_inf["leak_node_name"]) + # open the pipe + # leak_pipe_self.status = 'Open' + wn.add_pipe( + whole_inf["leak_pipe"], + start_node_name=whole_inf["pipe_start_node"], + end_node_name=whole_inf["pipe_end_node"], + length=whole_inf["pipe_length"], + diameter=whole_inf["pipe_diameter"], + roughness=whole_inf["pipe_roughness"], + minor_loss=whole_inf["pipe_minor_loss"], + ) + + return wn + + +def disable_all_controls_temporarily(wn): + """返回(控制名, 控制对象)的列表，之后可用 restore_controls 还原。""" + removed = [] + # WNTR 的控制都在 wn.control_name_list / wn.get_control / wn.remove_control + for cname in list(wn.control_name_list): + ctrl = wn.get_control(cname) + removed.append((cname, ctrl)) + wn.remove_control(cname) + return removed + + +def restore_controls(wn, removed): + """把先前禁用的控制全部加回去。""" + for cname, ctrl in removed: + wn.add_control(cname, ctrl) + + +def set_pipe2leaknode_mapping(mapping): + global _PIPE2LEAKNODE + _PIPE2LEAKNODE = mapping + + +def _get_or_create_leak_demand_ts(leak_node): + """ + 返回：泄漏专用 demand 的下标 idx。 + 若不存在，以 category='leak' 新建一条 base=0.0 的 demand。 + """ + # 先尝试找到已有的 'leak' 分类 + for i, ts in enumerate(leak_node.demand_timeseries_list): + # WNTR 的 Demand object 存在 category 属性 + if getattr(ts, "category", None) == "leak": + return i + # 没有则新建（base=0.0，后续临时改 base_value） + leak_node.add_demand(base=0.0, pattern_name=None, category="leak") + return len(leak_node.demand_timeseries_list) - 1 + + +def ensure_mid_node(wn, leak_pipe): + add_pipe1 = f"{leak_pipe}A" + add_pipe2 = f"{leak_pipe}B" + add_node = f"{leak_pipe}__mid" + + if add_node in wn.node_name_list: + return add_node + + if leak_pipe in wn.link_name_list: + leak_pipe_self = wn.get_link(leak_pipe) + pipe_diameter = leak_pipe_self.diameter + pipe_length = leak_pipe_self.length + pipe_roughness = leak_pipe_self.roughness + pipe_minor_loss = leak_pipe_self.minor_loss + pipe_start_node = leak_pipe_self.start_node_name + pipe_end_node = leak_pipe_self.end_node_name + + start_n = wn.get_node(pipe_start_node) + end_n = wn.get_node(pipe_end_node) + if start_n.node_type == "Reservoir": + end_elev = end_n.elevation + start_elev = end_elev + elif end_n.node_type == "Reservoir": + start_elev = start_n.elevation + end_elev = start_elev + else: + end_elev = end_n.elevation + start_elev = start_n.elevation + elev_mid = (start_elev + end_elev) / 2.0 + x_mid = (start_n.coordinates[0] + end_n.coordinates[0]) / 2.0 + y_mid = (start_n.coordinates[1] + end_n.coordinates[1]) / 2.0 + + wn.remove_link(leak_pipe) + wn.add_junction( + add_node, base_demand=0.0, elevation=elev_mid, coordinates=(x_mid, y_mid) + ) + wn.add_pipe( + add_pipe1, + start_node_name=pipe_start_node, + end_node_name=add_node, + length=pipe_length / 2.0, + diameter=pipe_diameter, + roughness=pipe_roughness, + minor_loss=pipe_minor_loss, + ) + wn.add_pipe( + add_pipe2, + start_node_name=add_node, + end_node_name=pipe_end_node, + length=pipe_length / 2.0, + diameter=pipe_diameter, + roughness=pipe_roughness, + minor_loss=pipe_minor_loss, + ) + return add_node + + # 若 A/B 已存在但中点不在，建议确认网络一致性 + raise KeyError(f"Cannot ensure mid node for pipe '{leak_pipe}'.") + + +def leak_simulation_pipe_dd_multi_pf(wn, leak_mag, leak_pipe, sensor_name): + """ + 优化版： + - 不再 remove/add link/node + - 直接在预插入的中点泄漏节点上设置 base_demand = leak_mag；仿真后设回 0 + """ + wn.options.hydraulic.demand_model = "DD" + + # 确保中点节点存在 + leak_node_name = ensure_mid_node(wn, leak_pipe) + leak_node = wn.get_node(leak_node_name) + + # 拿到泄漏专用的 demand time-series 下标 + leak_idx = _get_or_create_leak_demand_ts(leak_node) + ts_obj = leak_node.demand_timeseries_list[leak_idx] + + # 记录原值（通常是 0.0） + orig_base = ts_obj.base_value + try: + # 打开泄漏：只改 base_value，不碰 base_demand（只读） + ts_obj.base_value = float(leak_mag) + + # 仿真 + sim = wntr.sim.EpanetSimulator(wn) + results = sim.run_sim() + + # 输出（保持列顺序） + pressure_output = results.node["pressure"].loc[:, sensor_name] + # flow_output = results.link['flowrate'].loc[:, sensor_f_name] + return wn, pressure_output + finally: + # 关闭泄漏：还原 base_value + ts_obj.base_value = orig_base + + +def prepare_leak_infrastructure(wn, candidate_pipes): + """ + 把 candidate_pipes 每条管段切成两段，并在中点插入一个泄漏节点（base_demand=0）。 + 返回一个映射：pipe_id -> leak_node_name + 注意：只做一次；后续仿真通过在该节点设置 base_demand 实现“打开泄漏”，结束后恢复为 0。 + """ + pipe2leaknode = {} + for leak_pipe in candidate_pipes: + if leak_pipe in pipe2leaknode: + continue + + leak_pipe_self = wn.get_link(leak_pipe) + pipe_diameter = leak_pipe_self.diameter + pipe_length = leak_pipe_self.length + pipe_roughness = leak_pipe_self.roughness + pipe_minor_loss = leak_pipe_self.minor_loss + pipe_start_node = leak_pipe_self.start_node_name + pipe_end_node = leak_pipe_self.end_node_name + + # 计算中点高程/坐标（与原逻辑一致） + start_n = wn.get_node(pipe_start_node) + end_n = wn.get_node(pipe_end_node) + if start_n.node_type == "Reservoir": + end_elev = end_n.elevation + start_elev = end_elev + elif end_n.node_type == "Reservoir": + start_elev = start_n.elevation + end_elev = start_elev + else: + end_elev = end_n.elevation + start_elev = start_n.elevation + elev_mid = (start_elev + end_elev) / 2.0 + x_mid = (start_n.coordinates[0] + end_n.coordinates[0]) / 2.0 + y_mid = (start_n.coordinates[1] + end_n.coordinates[1]) / 2.0 + + # 先删原管，再加中点与两段半长管（只做一次） + wn.remove_link(leak_pipe) + + add_pipe1 = f"{leak_pipe}A" + add_pipe2 = f"{leak_pipe}B" + add_node = f"{leak_pipe}__mid" # 唯一命名，后面直接用它当泄漏节点 + + wn.add_junction( + add_node, base_demand=0.0, elevation=elev_mid, coordinates=(x_mid, y_mid) + ) + wn.add_pipe( + add_pipe1, + start_node_name=pipe_start_node, + end_node_name=add_node, + length=pipe_length / 2.0, + diameter=pipe_diameter, + roughness=pipe_roughness, + minor_loss=pipe_minor_loss, + ) + wn.add_pipe( + add_pipe2, + start_node_name=add_node, + end_node_name=pipe_end_node, + length=pipe_length / 2.0, + diameter=pipe_diameter, + roughness=pipe_roughness, + minor_loss=pipe_minor_loss, + ) + + pipe2leaknode[leak_pipe] = add_node + + return pipe2leaknode + + +def normal_simulation_pf( + wn, drive_mode, sensor_name, sensor_f_name, inp_time, require_p, minimum_p +): + # inp_time = 0 + if drive_mode == "PDD": # 需水量根据节点压力动态调整 + wn.options.hydraulic.demand_model = "PDD" + wn.options.hydraulic.required_pressure = require_p + wn.options.hydraulic.minimum_pressure = minimum_p + elif drive_mode == "DD": # 需水量固定，与压力无关 + wn.options.hydraulic.demand_model = "DD" + sim = wntr.sim.EpanetSimulator(wn) + results = sim.run_sim() + pressure_all = results.node["pressure"][sensor_name] + pressure = pressure_all.iloc[inp_time] + demand_all = results.node["demand"] + demand = demand_all.iloc[inp_time] + sum_demand = cal_sum_demand(demand) + flow_all = results.link["flowrate"][sensor_f_name] + flow = flow_all.iloc[inp_time] + + top_sensor = pressure.idxmin() + basic_p = results.node["pressure"] + basic_p = basic_p.iloc[inp_time] + return pressure, flow, basic_p, top_sensor, sum_demand + + +def normal_simulation_multi_pf( + wn, drive_mode, sensor_name, sensor_f_name, require_p, minimum_p +): + # inp_time = 0 + if drive_mode == "PDD": + wn.options.hydraulic.demand_model = "PDD" + wn.options.hydraulic.required_pressure = require_p + wn.options.hydraulic.minimum_pressure = minimum_p + elif drive_mode == "DD": + wn.options.hydraulic.demand_model = "DD" + sim = wntr.sim.EpanetSimulator(wn) + results = sim.run_sim() + pressure_all = results.node["pressure"][sensor_name] + pressure = pressure_all + demand_all = results.node["demand"] + demand = demand_all + flow = results.link["flowrate"][sensor_f_name] + sum_demand = pd.Series(dtype=object) + for i in range(len(demand.index)): + sum_demand[str(demand.index[i])] = cal_sum_demand(demand.iloc[i]) + if type(pressure) == pd.core.series.Series: + top_sensor = pressure.idxmin() + else: + mean_pressure = pressure.mean() + top_sensor = mean_pressure.idxmin() + basic_p = results.node["pressure"] + return pressure, flow, basic_p, top_sensor, sum_demand + + +def simple_simulation_pf(wn, sensor_name, sensor_f_name, leak_pipe, add_pipe1): + sim = wntr.sim.EpanetSimulator(wn) + results = sim.run_sim() + pressure_all = results.node["pressure"][sensor_name] + if len(leak_pipe) > 0 and leak_pipe in sensor_f_name: + + f_sensor_name = [add_pipe1 if i == leak_pipe else i for i in sensor_f_name] + flow_all = results.link["flowrate"][f_sensor_name] + flow_all.columns = sensor_f_name + else: + flow_all = results.link["flowrate"][sensor_f_name] + return pressure_all, flow_all + + +def cal_sum_demand(demand): + sum_demand = 0 + for i in range(len(demand)): + if demand.iloc[i] > 0: + sum_demand += demand.iloc[i] + return sum_demand + + +def cal_signature_pipe_multi_pf( + wn, leak_mag, candidate_center, timestep_list, sensor_name +): + candidate_center_num = len(candidate_center) + pressure_leak = pd.DataFrame( + index=pd.MultiIndex.from_product([candidate_center, timestep_list]), + columns=sensor_name, + ) + # flow_leak = pd.DataFrame(index=pd.MultiIndex.from_product([candidate_center, timestep_list]), + # columns=sensor_f_name) + pressure_leak = pressure_leak.sort_index() + # flow_leak = flow_leak.sort_index() + for i in range(candidate_center_num): + wn, pressure_output = leak_simulation_pipe_dd_multi_pf( + wn, leak_mag, candidate_center[i], sensor_name + ) + # leak_or_not_list.append(leak_or_not) + pressure_leak.loc[candidate_center[i]].loc[:, :] = pressure_output + # flow_leak.loc[candidate_center[i]].loc[:, :] = flow_output + sys.stdout.write("\r" + "已经完成计算" + str(i + 1) + "个特征中心") + return pressure_leak, candidate_center + + +def pick_pipe(all_pipes, pipe_diameter, limited_diameter): + candidate_pipe = [] + for each_pipe in all_pipes: + if pipe_diameter[each_pipe] >= limited_diameter: + candidate_pipe.append(each_pipe) + return candidate_pipe + + +def cal_possible_pipe(leak_flow, all_pipe, pipe_diameter): + basic_pressure = 10 # 基础压力 + discharge_coeff = 0.6 # 经验系数 + break_area_ratio = 1 # 爆管面积比 0.5 1.25 + break_area = leak_flow / ( + discharge_coeff * math.sqrt(2 * basic_pressure * 9.81) + ) # 爆管面积 m3/h + """break_area_diameter = math.sqrt(4 * break_area / math.pi) + min_diameter = (math.ceil(1000 * break_area_diameter / break_area_ratio)) / 1000""" + break_area_diameter = math.sqrt( + 4 * break_area / math.pi / break_area_ratio + ) # 爆管直径 + min_diameter = (math.ceil(1000 * break_area_diameter)) / 1000 # 向上取整 + new_all_pipe = pick_pipe(all_pipe, pipe_diameter, min_diameter) + return new_all_pipe, min_diameter + + +def extract_links(data, link_types, direction): + return [ + link + for res_data in data.values() + for link_type in link_types + for link in res_data[link_type][direction] + ] + + +class LeakSimulator: + @staticmethod + def simple_add_leak(*args, **kwargs): + return simple_add_leak(*args, **kwargs) + + @staticmethod + def simple_recover_wn(*args, **kwargs): + return simple_recover_wn(*args, **kwargs) + + @staticmethod + def leak_simulation_pipe_dd_multi_pf(*args, **kwargs): + return leak_simulation_pipe_dd_multi_pf(*args, **kwargs) + + @staticmethod + def normal_simulation_pf(*args, **kwargs): + return normal_simulation_pf(*args, **kwargs) + + @staticmethod + def normal_simulation_multi_pf(*args, **kwargs): + return normal_simulation_multi_pf(*args, **kwargs) + + @staticmethod + def simple_simulation_pf(*args, **kwargs): + return simple_simulation_pf(*args, **kwargs) + + @staticmethod + def cal_sum_demand(*args, **kwargs): + return cal_sum_demand(*args, **kwargs) + + @staticmethod + def cal_signature_pipe_multi_pf(*args, **kwargs): + return cal_signature_pipe_multi_pf(*args, **kwargs) + + @staticmethod + def cal_possible_pipe(*args, **kwargs): + return cal_possible_pipe(*args, **kwargs) + + @staticmethod + def pick_pipe(*args, **kwargs): + return pick_pipe(*args, **kwargs) + + @staticmethod + def extract_links(*args, **kwargs): + return extract_links(*args, **kwargs) diff --git a/app/algorithms/burst_location/network_model.py b/app/algorithms/burst_location/network_model.py new file mode 100644 index 0000000..4c22513 --- /dev/null +++ b/app/algorithms/burst_location/network_model.py @@ -0,0 +1,169 @@ +"""管网模型读取与图构建模块。""" + +import copy + +import numpy as np +import networkx as nx +import pandas as pd +import wntr + + +def load_inp(inp_name, inp_location, inp_time, driven_mode, require_p, minimum_p): + inp_file = inp_location + inp_name + wn = wntr.network.WaterNetworkModel(inp_file) + if driven_mode == "PDD": + wn.options.hydraulic.demand_model = "PDD" + wn.options.hydraulic.required_pressure = require_p + wn.options.hydraulic.minimum_pressure = minimum_p + else: + wn.options.hydraulic.demand_model = "DD" + return wn + + +def read_inf_inp(wn): + all_node = wn.node_name_list + node_elevation = wn.query_node_attribute("elevation") + node_coordinates = wn.query_node_attribute("coordinates") + + all_pipe = wn.pipe_name_list + # 改_wz__________________________________ + n_pipe = [] + for p in all_pipe: + pipe = wn.get_link(p) + if pipe.initial_status == 0: # 状态为'Closed' + n_pipe.append(p) + candidate_pipe_init = list(set(all_pipe) - set(n_pipe)) + pipe_start_node = wn.query_link_attribute( + "start_node_name", link_type=wntr.network.model.Pipe + ) + pipe_end_node = wn.query_link_attribute( + "end_node_name", link_type=wntr.network.model.Pipe + ) + pipe_length = wn.query_link_attribute("length") + pipe_diameter = wn.query_link_attribute("diameter") + + return ( + all_node, + node_elevation, + node_coordinates, + candidate_pipe_init, + pipe_start_node, + pipe_end_node, + pipe_length, + pipe_diameter, + ) + + +def read_inf_inp_other(wn): + all_link = wn.link_name_list + pipe_start_node_all = wn.query_link_attribute("start_node_name") + pipe_end_node_all = wn.query_link_attribute("end_node_name") + + return all_link, pipe_start_node_all, pipe_end_node_all + + +def construct_graph(wn): + length = wn.query_link_attribute("length") + G = wn.get_graph(wn, link_weight=length) + # 转为无向图 + G0 = G.to_undirected() + # A0 = np.array(nx.adjacency_graph(G0).todense()) + return G0 # , A0 + + +def cal_pipe_coordinate(all_pipe, pipe_start_node, pipe_end_node, node_coordinates): + pipe_num = len(all_pipe) + pipe_coordinates = np.zeros([pipe_num, 2]) + pipe_x = copy.deepcopy(pipe_start_node) + pipe_y = copy.deepcopy(pipe_start_node) + for i in range(pipe_num): + temp_pipe = all_pipe[i] + pipe_x[temp_pipe] = ( + node_coordinates[pipe_start_node[temp_pipe]][0] + + node_coordinates[pipe_end_node[temp_pipe]][0] + ) / 2 + pipe_y[temp_pipe] = ( + node_coordinates[pipe_start_node[temp_pipe]][1] + + node_coordinates[pipe_end_node[temp_pipe]][1] + ) / 2 + return pipe_x, pipe_y + + +def cal_node_coordinate(all_node, node_coordinates): + node_x = copy.deepcopy(node_coordinates) + node_y = copy.deepcopy(node_coordinates) + for i in range(len(node_x)): + temp_node = all_node[i] + node_x[temp_node] = node_coordinates[temp_node][0] + node_y[temp_node] = node_coordinates[temp_node][1] + return node_x, node_y + + +def produce_pattern_value(wn, all_node): + wn_o = copy.deepcopy(wn) + # 改_wz_____________________________ + # sample_node = wn_o.get_node(all_node[0]) + # num_categories = len(sample_node.demand_timeseries_list) + num_categories = 1 + columns = [f"D{i}" for i in range(num_categories)] + basic_demand_pd = pd.DataFrame(index=all_node, columns=columns) + for each in all_node: + node = wn_o.get_node(each) + for i in range(num_categories): + basic_demand_pd.loc[each, columns[i]] = node.demand_timeseries_list[ + i + ].base_value + + return basic_demand_pd + + +def _build_node_pipe_maps( + all_nodes, candidate_pipes, pipe_start_node, pipe_end_node, pipe_length +): + node_pipe_dic = {node: [] for node in all_nodes} + couple_node_length = {} + for pipe in candidate_pipes: + start_node = pipe_start_node[pipe] + end_node = pipe_end_node[pipe] + if start_node in node_pipe_dic: + node_pipe_dic[start_node].append(pipe) + if end_node in node_pipe_dic: + node_pipe_dic[end_node].append(pipe) + length = float(pipe_length[pipe]) + couple_node_length[f"{start_node},{end_node}"] = length + couple_node_length[f"{end_node},{start_node}"] = length + return node_pipe_dic, couple_node_length + + +class NetworkModelReader: + @staticmethod + def load_inp(*args, **kwargs): + return load_inp(*args, **kwargs) + + @staticmethod + def read_inf_inp(*args, **kwargs): + return read_inf_inp(*args, **kwargs) + + @staticmethod + def read_inf_inp_other(*args, **kwargs): + return read_inf_inp_other(*args, **kwargs) + + @staticmethod + def construct_graph(*args, **kwargs): + return construct_graph(*args, **kwargs) + + @staticmethod + def cal_pipe_coordinate(*args, **kwargs): + return cal_pipe_coordinate(*args, **kwargs) + + @staticmethod + def cal_node_coordinate(*args, **kwargs): + return cal_node_coordinate(*args, **kwargs) + + @staticmethod + def produce_pattern_value(*args, **kwargs): + return produce_pattern_value(*args, **kwargs) + + @staticmethod + def build_node_pipe_maps(*args, **kwargs): + return _build_node_pipe_maps(*args, **kwargs) diff --git a/app/algorithms/burst_location/network_partitioner.py b/app/algorithms/burst_location/network_partitioner.py new file mode 100644 index 0000000..5cd63b9 --- /dev/null +++ b/app/algorithms/burst_location/network_partitioner.py @@ -0,0 +1,375 @@ +"""管网分区模块。""" + +import math + +import matplotlib.pyplot as plt +import networkx as nx +import networkx as networkx +import numpy as np +import pandas as pd +import pymetis +from scipy.sparse import coo_matrix, csr_matrix +from scipy.sparse.csgraph import connected_components + + +def pick_center_pipe(node_x, node_y, candidate_pipe, pipe_start_node, pipe_end_node): + data_set_t = pd.DataFrame(dtype=object) + data_set_t["x"] = ( + node_x[pipe_start_node[candidate_pipe]].values + + node_x[pipe_start_node[candidate_pipe]].values + ) / 2 + data_set_t["y"] = ( + node_y[pipe_end_node[candidate_pipe]].values + + node_y[pipe_end_node[candidate_pipe]].values + ) / 2 + data_set_t.index = list(candidate_pipe) + mean_x = data_set_t["x"].mean() + mean_y = data_set_t["y"].mean() + data_set_t["d"] = abs(data_set_t["x"] - mean_x) + abs(data_set_t["y"] - mean_y) + distance_t = data_set_t["d"].sort_values(ascending=True, inplace=False) + """if distance_t.index==[]: + print(candidate_pipe) + else:""" + center_t = distance_t.index[0] + return center_t + + +def find_new_center_pipe( + node_x, node_y, candidate_pipe, pipe_start_node, pipe_end_node, record_center +): + new_candidate_pipe = list(set(candidate_pipe) - set(record_center)) + if new_candidate_pipe == []: + new_candidate_pipe = candidate_pipe + center_t = pick_center_pipe( + node_x, node_y, new_candidate_pipe, pipe_start_node, pipe_end_node + ) + return center_t + + +def cal_area_node_linked_pipe(nodeset, node_pipe_dic): + pipeset = [] + nodeset = list(nodeset) + for i in range(len(nodeset)): + temp_node = nodeset[i] + pipe = node_pipe_dic[temp_node] + pipeset = pipeset + pipe + return pipeset + + +def metis_grouping_pipe_weight( + G0, + wn, + all_node_iter, + candidate_pipe_input, + group_num, + node_x, + node_y, + pipe_start_node_all, + pipe_end_node_all, + node_pipe_dic, + all_node_series, + couple_node_length, +): + all_node_iter_series_new = all_node_series[all_node_iter] + all_node_iter_series_new = all_node_iter_series_new.sort_values(ascending=True) + all_node_iter_new = list(all_node_iter_series_new.index) + + G1 = G0.subgraph(all_node_iter_new) + delimiter = " " + adjacency_list = [] + node_dict = {} + c_new = 0 + for each_node in all_node_iter_new: + node_dict[each_node] = c_new + c_new = c_new + 1 + correspond_dic = {} + count_node = 0 + w = [] + for line in generate_adjlist_with_all_edges(G1, delimiter): + temp_node_name = line.split(sep=delimiter) + w_temp = [] + for i in range(len(temp_node_name) - 1): + temp_name_1 = temp_node_name[0] + "," + temp_node_name[i + 1] + w_temp.append(couple_node_length[temp_name_1]) + w.append(w_temp) + n_t = [] + for each_node in temp_node_name: + n_t.append(node_dict[each_node]) + correspond_dic[n_t[0]] = count_node + count_node = count_node + 1 + # del n_t[0] + adjacency_list.append(n_t) + adjacency_list_new = [[] * 1 for i in range(len(adjacency_list))] + w_new = [[] * 1 for i in range(len(adjacency_list))] + for i in range(len(adjacency_list)): + adjacency_list_new[int(adjacency_list[i][0])] = adjacency_list[i] + w_new[int(adjacency_list[i][0])] = w[i] + + for i in range(len(adjacency_list)): + del adjacency_list_new[i][0] + xadj = [0] + w_f = [] + final_adjacency_list = [] + for i in range(len(adjacency_list_new)): + final_adjacency_list = final_adjacency_list + adjacency_list_new[i] + xadj.append(len(final_adjacency_list)) + w_f = w_f + w_new[i] + + # (edgecuts, parts) = pymetis.part_graph(nparts=group_num, adjacency=adjacency_list_new) + (edgecuts, parts) = pymetis.part_graph( + nparts=group_num, adjncy=final_adjacency_list, xadj=xadj, eweights=w_f + ) + # (edgecuts, parts) = pymetis.part_graph(nparts=group_num, adjacency=adjacency_list_new) + candidate_group_list = [[] * 1 for i in range(group_num)] + for i in range(len(all_node_iter_new)): + candidate_group_list[parts[i]].append(all_node_iter_new[i]) + + """parts_new = np.zeros(len(candidate_node_input), dtype=int) + for i in range(len(candidate_group_list)): + temp_group = candidate_group_list[i] + for each_node in temp_group: + parts_new[node_dict[each_node]] = i + parts_new = list(parts_new)""" + + new_center = [] + new_group = [] + new_all_node = [] + candidate_pipe_set = set(candidate_pipe_input) + all_grouped_pipe = [] + for i in range(group_num): + # 构建子图 + G_sub = G0.subgraph(candidate_group_list[i]) + + # 计算联通子图 + sub_graphs = networkx.connected_components(G_sub) + if networkx.number_connected_components(G_sub) == 1: + # 求交集 + nodeset = G_sub.nodes() + pipeset_set = set(cal_area_node_linked_pipe(nodeset, node_pipe_dic)) + candidate_pipe = list(pipeset_set.intersection(candidate_pipe_set)) + + # 判断集合是否保留 + if len(candidate_pipe) > 0: + # 保留 计算中心 + center_t = pick_center_pipe( + node_x, + node_y, + candidate_pipe, + pipe_start_node_all, + pipe_end_node_all, + ) + + # 更新 + new_center.append(center_t) + new_group.append(candidate_pipe) + new_all_node.append(nodeset) + all_grouped_pipe = all_grouped_pipe + candidate_pipe + else: + for c in sub_graphs: + G_temp = G0.subgraph(c) + nodeset = G_temp.nodes() + pipeset = cal_area_node_linked_pipe(nodeset, node_pipe_dic) + pipeset_set = set(pipeset) + + # 求交集 + candidate_pipe = list(pipeset_set.intersection(candidate_pipe_set)) + # print(len(candidate_node)) + # 判断集合是否保留 + if len(candidate_pipe) > 0: + # 保留 计算中心 + center_t = pick_center_pipe( + node_x, + node_y, + candidate_pipe, + pipe_start_node_all, + pipe_end_node_all, + ) + # 更新 + new_center.append(center_t) + new_group.append(candidate_pipe) + new_all_node.append(nodeset) + all_grouped_pipe = all_grouped_pipe + candidate_pipe + record_center = [] + c_g = 0 + for each_group in new_group: + if len(each_group) < 3: + record_center.append(new_center[c_g]) + c_g += 1 + c_g = 0 + for each_group in new_group: + if len(each_group) >= 3: + if new_center[c_g] in record_center: + new_center[c_g] = find_new_center_pipe( + node_x, + node_y, + each_group, + pipe_start_node_all, + pipe_end_node_all, + record_center, + ) + record_center.append(new_center[c_g]) + c_g += 1 + + # visualize_metis_partition( + # G0, new_center, new_group, + # node_x, node_y, + # pipe_start_node_all, pipe_end_node_all + # ) + return new_center, new_group, new_all_node + + +def visualize_metis_partition( + G, center_pipes, pipe_groups, node_x, node_y, pipe_start_node_all, pipe_end_node_all +): + """ + 可视化METIS分区结果（单图模式） + 参数： + G: 原始管网图(nx.Graph) + center_pipes: 中心管道列表(list) + pipe_groups: 分组管道列表(list of lists) + node_x: 节点X坐标字典(dict) + node_y: 节点Y坐标字典(dict) + pipe_start_node_all: 管道起点字典(dict) + pipe_end_node_all: 管道终点字典(dict) + """ + plt.figure(figsize=(9, 10)) + + # 生成颜色映射（自动扩展颜色数量） + colors = plt.cm.tab20(np.linspace(0, 1, len(pipe_groups))) + + # --- 绘制背景管网（灰色半透明） --- + for edge in G.edges(): + start_node, end_node = edge + plt.plot( + [node_x[start_node], node_x[end_node]], + [node_y[start_node], node_y[end_node]], + color="lightgray", + linewidth=0.5, + alpha=0.3, + zorder=1, # 确保背景在底层 + ) + + # --- 绘制各分区管道（彩色）--- + legend_handles = [] # 用于图例的句柄 + for i, (group, center) in enumerate(zip(pipe_groups, center_pipes)): + color = colors[i % len(colors)] # 循环使用颜色 + + # 绘制分组管道 + for pipe in group: + start = pipe_start_node_all[pipe] + end = pipe_end_node_all[pipe] + line = plt.plot( + [node_x[start], node_x[end]], + [node_y[start], node_y[end]], + color=color, + linewidth=2.5, + alpha=0.8, + zorder=2, + ) + # 只为每个分组的第一个管道添加图例句柄 + if pipe == group[0]: + legend_handles.append(line[0]) + + # 高亮中心管道（红色虚线） + if center in pipe_start_node_all and center in pipe_end_node_all: + start = pipe_start_node_all[center] + end = pipe_end_node_all[center] + plt.plot( + [node_x[start], node_x[end]], + [node_y[start], node_y[end]], + color="red", + linewidth=4, + linestyle="--", + dash_capstyle="round", + zorder=3, # 确保中心管道在最顶层 + ) + + # --- 添加图例和标注 --- + # 分组图例 + group_labels = [f"Group {i + 1}" for i in range(len(pipe_groups))] + plt.legend( + legend_handles, + group_labels, + loc="upper right", + title="Partitions", + fontsize=8, + title_fontsize=10, + ) + + # 中心管道标注（可选） + for i, center in enumerate(center_pipes): + if center in pipe_start_node_all: + x = ( + node_x[pipe_start_node_all[center]] + node_x[pipe_end_node_all[center]] + ) / 2 + y = ( + node_y[pipe_start_node_all[center]] + node_y[pipe_end_node_all[center]] + ) / 2 + plt.text( + x, + y, + f"C{i + 1}", + color="red", + fontsize=10, + ha="center", + va="center", + bbox=dict(facecolor="white", alpha=0.8, edgecolor="none"), + ) + + # --- 图形美化 --- + plt.title("Water Network Partitioning Overview", fontsize=14, pad=20) + plt.xlabel("X Coordinate", fontsize=10) + plt.ylabel("Y Coordinate", fontsize=10) + plt.grid(True, alpha=0.2, linestyle=":") + plt.tight_layout() + + # 显示图形 + plt.show() + + +def generate_adjlist_with_all_edges(G, delimiter): + for s, nbrs in G.adjacency(): + line = str(s) + delimiter + for t, data in nbrs.items(): + line += str(t) + delimiter + yield line[: -len(delimiter)] + + +def cal_group_num(candidate_node_input, cal_group_num): + candidate_node_num = len(candidate_node_input) + if candidate_node_num > 100: + group_num_input = cal_group_num # 30 + else: + group_num_input = 10 + + return group_num_input + + +class NetworkPartitioner: + @staticmethod + def metis_grouping_pipe_weight(*args, **kwargs): + return metis_grouping_pipe_weight(*args, **kwargs) + + @staticmethod + def visualize_metis_partition(*args, **kwargs): + return visualize_metis_partition(*args, **kwargs) + + @staticmethod + def generate_adjlist_with_all_edges(*args, **kwargs): + return generate_adjlist_with_all_edges(*args, **kwargs) + + @staticmethod + def pick_center_pipe(*args, **kwargs): + return pick_center_pipe(*args, **kwargs) + + @staticmethod + def find_new_center_pipe(*args, **kwargs): + return find_new_center_pipe(*args, **kwargs) + + @staticmethod + def cal_area_node_linked_pipe(*args, **kwargs): + return cal_area_node_linked_pipe(*args, **kwargs) + + @staticmethod + def cal_group_num(*args, **kwargs): + return cal_group_num(*args, **kwargs) diff --git a/app/algorithms/burst_location/noise_generator.py b/app/algorithms/burst_location/noise_generator.py new file mode 100644 index 0000000..b3fa9eb --- /dev/null +++ b/app/algorithms/burst_location/noise_generator.py @@ -0,0 +1,234 @@ +"""噪声生成模块。""" + +import copy +import random + +import numpy as np +import pandas as pd + +from .leak_simulator import simple_add_leak, simple_recover_wn, simple_simulation_pf + + +def add_noise_pd(data, noise_type, noise_para): + output_data = copy.deepcopy(data) + if type(output_data) == pd.core.frame.Series: + if noise_type == "uni": + for x in output_data.index: + noise = (np.random.random() - 0.5) * 2 + output_data[x] = output_data[x] + noise * noise_para + + elif noise_type == "gauss": + noise = np.random.normal(loc=0, scale=noise_para, size=output_data.shape) + output_data = output_data + noise + elif type(output_data) == pd.core.frame.DataFrame: + if noise_type == "uni": + noise = (np.random.random(size=output_data.shape) - 0.5) * 2 + output_data = output_data + noise * noise_para + + elif noise_type == "gauss": + noise = np.random.normal(loc=0, scale=noise_para, size=output_data.shape) + output_data = output_data + noise + return output_data + + +def add_noise_number(data, noise_type, noise_para): + output_data = copy.deepcopy(data) + if noise_type == "uni": + noise = (np.random.random() - 0.5) * 2 + output_data = output_data + noise * noise_para + + elif noise_type == "gauss": + noise = random.gauss(0, noise_para) + output_data = output_data + noise + return output_data + + +def add_noise_number_flow(data, noise_para_mean, noise_para_std1, noise_para_std2): + output_data = copy.deepcopy(data) + noise_flag1 = np.random.random() - 0.5 + if noise_flag1 < 0: + noise = noise_para_mean - abs(np.random.normal(loc=0, scale=noise_para_std1)) + else: + noise = noise_para_mean + abs(np.random.normal(loc=0, scale=noise_para_std2)) + noise_flag2 = np.random.random() - 0.5 + if noise_flag2 < 0: + noise_f = noise * (-1) + else: + noise_f = noise + output_data = output_data + noise_f + return output_data + + +def produce_noise_number(noise_type, noise_para): + if noise_type == "uni": + noise = (np.random.random() - 0.5) * 2 + noise = noise * noise_para + elif noise_type == "gauss": + noise = random.gauss(0, noise_para) + else: + noise = 0 + return noise + + +def add_noise_percentage_pd(data, noise_type, noise_para): + output_data = copy.deepcopy(data) + + if type(output_data) == pd.core.frame.Series: + if noise_type == "uni": + for x in output_data.index: + noise = (np.random.random() - 0.5) * 2 + output_data[x] = output_data[x] * (1 + noise * noise_para / 100) + + elif noise_type == "gauss": + for x in output_data.index: + noise = np.random.gauss(0, noise_para) + output_data[x] = output_data[x] * (1 + noise / 100) + # std_noise = noise.std() + elif type(output_data) == pd.core.frame.DataFrame: + if noise_type == "uni": + noise = (np.random.random(size=output_data.shape) - 0.5) * 2 + output_data = output_data * (1 + noise * noise_para / 100) + + elif noise_type == "gauss": + noise = np.random.normal(loc=0, scale=noise_para, size=output_data.shape) + output_data = output_data * (1 + noise / 100) + # std_noise = noise.std().mean() + return output_data + + +def add_noise_in_wn_pf( + wn, + pipe_c_noise, + timestep_list, + pipe_coefficient, + sensor_name, + sensor_f_name, + all_node, + basic_demand_pd, + noise_type, + noise_para, + leak_pipe, + leak_flow, +): + + wn.options.time.duration = 0 + pipe_roughness_change = add_noise_pd(pipe_coefficient, noise_type, pipe_c_noise) + wn = change_para_of_wn(wn, pipe_roughness_change) + record_pressure = pd.DataFrame(index=timestep_list, columns=sensor_name) + record_flow = pd.DataFrame(index=timestep_list, columns=sensor_f_name) + record_noise_all = pd.DataFrame( + index=pd.MultiIndex.from_product([timestep_list, all_node]), + columns=basic_demand_pd.columns, + ) + record_noise_all = record_noise_all.sort_index() + + # normal 获取添加噪声后的监测点数据 + for i in range(len(timestep_list)): + wn, record_noise = change_node_demand( + wn, basic_demand_pd, all_node, noise_type, noise_para + ) + record_noise_all.loc[timestep_list[i]].loc[:, :] = record_noise + pressure_temp, flow_temp = simple_simulation_pf( + wn, sensor_name, sensor_f_name, [], [] + ) + record_pressure.iloc[i, :] = pressure_temp + record_flow.iloc[i, :] = flow_temp + + # leak_simulation 获取添加漏损后的监测点数据 + record_pressure_leak = pd.DataFrame(index=timestep_list, columns=sensor_name) + record_flow_leak = pd.DataFrame(index=timestep_list, columns=sensor_f_name) + # 改_wz_________________________________________ + # add leak + wn, whole_inf, add_pipe1 = simple_add_leak(wn, leak_flow, leak_pipe) + # simulation + for i in range(len(timestep_list)): + record_noise = record_noise_all.loc[timestep_list[i]] + wn = change_node_demand_leak(wn, record_noise, all_node) + pressure_temp, flow_temp = simple_simulation_pf( + wn, sensor_name, sensor_f_name, leak_pipe, add_pipe1 + ) + record_pressure_leak.iloc[i, :] = pressure_temp + record_flow_leak.iloc[i, :] = flow_temp + # delete leak + wn = simple_recover_wn(wn, whole_inf) + return wn, record_pressure, record_flow, record_pressure_leak, record_flow_leak + + +def change_node_demand(wn, basic_demand_pd, all_node, noise_type, noise_para): + # 改_wz_____________________________________ + record_noise = pd.DataFrame(index=all_node, columns=basic_demand_pd.columns) + for each_node in all_node: + node = wn.get_node(each_node) + num_columns = len(basic_demand_pd.columns) + # 处理前N-1列（如果有） + for i in range(num_columns - 1): + # 获取原始值并添加噪声 + record_noise.loc[each_node].iloc[i] = ( + 1 + produce_noise_number(noise_type, noise_para) + ) * basic_demand_pd.loc[each_node].iloc[i] + node.demand_timeseries_list[i].base_value = record_noise.loc[ + each_node + ].iloc[i] + # 处理最后一列（当列数>=1时） + if num_columns >= 1: + last_col = basic_demand_pd.columns[-1] + original_last = basic_demand_pd.loc[each_node, last_col] + record_noise.loc[each_node, last_col] = original_last + node.demand_timeseries_list[-1].base_value = original_last + + return wn, record_noise + + +def change_node_demand_leak(wn, record_noise, all_node): + sample_node = wn.get_node(all_node[0]) + # num_categories = len(sample_node.demand_timeseries_list) + num_categories = 1 + for each in all_node: + node = wn.get_node(each) + for i in range(num_categories): + node.demand_timeseries_list[i].base_value = record_noise.loc[each].iloc[i] + return wn + + +def change_para_of_wn(wn, pipe_roughness_change): + for pipe_name, pipe in wn.pipes(): + pipe.roughness = pipe_roughness_change[pipe_name] + return wn + + +class NoiseGenerator: + @staticmethod + def add_noise_pd(*args, **kwargs): + return add_noise_pd(*args, **kwargs) + + @staticmethod + def add_noise_number(*args, **kwargs): + return add_noise_number(*args, **kwargs) + + @staticmethod + def add_noise_number_flow(*args, **kwargs): + return add_noise_number_flow(*args, **kwargs) + + @staticmethod + def add_noise_percentage_pd(*args, **kwargs): + return add_noise_percentage_pd(*args, **kwargs) + + @staticmethod + def produce_noise_number(*args, **kwargs): + return produce_noise_number(*args, **kwargs) + + @staticmethod + def add_noise_in_wn_pf(*args, **kwargs): + return add_noise_in_wn_pf(*args, **kwargs) + + @staticmethod + def change_node_demand(*args, **kwargs): + return change_node_demand(*args, **kwargs) + + @staticmethod + def change_node_demand_leak(*args, **kwargs): + return change_node_demand_leak(*args, **kwargs) + + @staticmethod + def change_para_of_wn(*args, **kwargs): + return change_para_of_wn(*args, **kwargs) diff --git a/app/algorithms/burst_location/similarity_calculator.py b/app/algorithms/burst_location/similarity_calculator.py new file mode 100644 index 0000000..ff42828 --- /dev/null +++ b/app/algorithms/burst_location/similarity_calculator.py @@ -0,0 +1,830 @@ +"""相似性计算模块。""" + +import math + +import numpy as np +import pandas as pd + + +def cal_similarity_simple_return_dd( + similarity_mode, + monitor_p, + predict_p, + normal_p, + leak_p, + monitor_p_all, + predict_p_all, + normal_p_all, + leak_p_all, + important_sensor, + mean_dpressure, + dpressure_std, + dpressure_std_all, + if_gy=0, + cos_or_flow=1, +): + # cos_or_flow 用于 CAF + dpressure_s = normal_p - leak_p + dpressure = predict_p - monitor_p + act_dpressure = pd.Series(dtype=object) + for i in range(len(leak_p.index)): + if dpressure_std.iloc[i] > -200: # 0.001: + if if_gy == 1: + act_dpressure[leak_p.index[i]] = ( + leak_p.iloc[i] - monitor_p.iloc[i] + ) / dpressure_std.iloc[i] + else: + act_dpressure[leak_p.index[i]] = leak_p.iloc[i] - monitor_p.iloc[i] + + if similarity_mode == "COS" or (similarity_mode == "CAF" and cos_or_flow == 1): + + """if leak_p.min()<0: + none_flag = 1 + similarity_cos = 0 + similarity_dis = 0 + else:""" + none_flag = 0 + sensor_for_cos = list( + set(dpressure_s.index).intersection(set(act_dpressure.index)) + ) + """if len(dpressure_s) ==0 or len(dpressure) ==0: + jj=9 + else:""" + try: + s1 = np.dot( + np.transpose(dpressure_s.loc[sensor_for_cos]), + dpressure.loc[sensor_for_cos], + ) + s2 = np.linalg.norm(dpressure_s.loc[sensor_for_cos]) * np.linalg.norm( + dpressure.loc[sensor_for_cos] + ) + if s2 == 0: + s2 = s2 + 0.0001 + similarity_cos = s1 / s2 + similarity_dis = 0 + except Exception as e: + print(dpressure_s) + print(sensor_for_cos) + print(act_dpressure) + print(dpressure_std) + print(dpressure) + + elif similarity_mode == "DIS" or (similarity_mode == "CAF" and cos_or_flow == 2): + """if leak_p.min()<0: + none_flag = 1 + else:""" + none_flag = 0 + important_sensor = list( + set(important_sensor).intersection(set(act_dpressure.index)) + ) + part_dpressure = dpressure_s[important_sensor] - dpressure[important_sensor] + similarity_pre_DIS = np.linalg.norm(part_dpressure) + # similarity_pre_DIS_later = 1 / (1 + similarity_pre_DIS) + similarity_dis = similarity_pre_DIS + similarity_cos = 0 + + elif similarity_mode == "CAD_new": + act_dpressure = leak_p - monitor_p + """if leak_p.min() < 0: + none_flag = 1 + similarity_cos = 0 + similarity_dis =0 + else:""" + none_flag = 0 + # cos + s1 = np.dot(np.transpose(dpressure_s), dpressure) + s2 = np.linalg.norm(dpressure_s) * np.linalg.norm(dpressure) + + if s2 == 0: + s2 = s2 + 0.0001 + similarity_cos = s1 / s2 + + # DIS + part_dpressure = act_dpressure.loc[important_sensor] + similarity_pre_DIS = np.linalg.norm(part_dpressure) + similarity_pre_DIS_later = 1 / (1 + similarity_pre_DIS) + similarity_dis = similarity_pre_DIS + + elif similarity_mode == "CAD_new_gy" or similarity_mode == "CDF": + # cos + sensor_for_cos = list( + set(dpressure_s.index).intersection(set(act_dpressure.index)) + ) + if len(sensor_for_cos) == 0 and len(dpressure_s) == 0: + similarity_cos = 0 + elif len(sensor_for_cos) == 0 and len(dpressure_s) > 0: + sensor_for_cos = list(dpressure_s.index) + none_flag = 0 + s1 = np.dot( + np.transpose(dpressure_s.loc[sensor_for_cos]), + dpressure.loc[sensor_for_cos], + ) + s2 = np.linalg.norm(dpressure_s.loc[sensor_for_cos]) * np.linalg.norm( + dpressure.loc[sensor_for_cos] + ) + + if s2 == 0: + s2 = s2 + 0.0001 + similarity_cos = s1 / s2 + else: + none_flag = 0 + s1 = np.dot( + np.transpose(dpressure_s.loc[sensor_for_cos]), + dpressure.loc[sensor_for_cos], + ) + s2 = np.linalg.norm(dpressure_s.loc[sensor_for_cos]) * np.linalg.norm( + dpressure.loc[sensor_for_cos] + ) + + if s2 == 0: + s2 = s2 + 0.0001 + similarity_cos = s1 / s2 + + # DIS + important_sensor_new = list( + set(important_sensor).intersection(set(act_dpressure.index)) + ) + if len(important_sensor_new) == 0: + important_sensor_new = important_sensor + act_dpressure = pd.Series(dtype=object) + for i in range(len(leak_p_all.index)): + # if dpressure_std.iloc [i] > -200: # 0.001: + if if_gy == 1: + act_dpressure[leak_p_all.index[i]] = ( + leak_p_all.iloc[i] - monitor_p_all.iloc[i] + ) / dpressure_std_all.iloc[i] + else: + act_dpressure[leak_p_all.index[i]] = ( + leak_p_all.iloc[i] - monitor_p_all.iloc[i] + ) + # part_dpressure = act_dpressure.loc[important_sensor_new] + + part_dpressure = ( + dpressure.loc[important_sensor_new] - dpressure_s.loc[important_sensor_new] + ) + similarity_pre_DIS = np.linalg.norm(part_dpressure) ## chang test + # part_dpressure = dpressure_s.loc[important_sensor]-dpressure.loc[important_sensor] + # similarity_pre_DIS = np.linalg.norm(part_dpressure) + # similarity_pre_DIS_later = 1 / (1 + similarity_pre_DIS) + similarity_dis = similarity_pre_DIS + elif similarity_mode == "OF": + # cos + similarity_cos = 0 + none_flag = 0 + # DIS + important_sensor_new = list( + set(important_sensor).intersection(set(act_dpressure.index)) + ) + if len(important_sensor_new) == 0: + important_sensor_new = important_sensor + act_dpressure = pd.Series(dtype=object) + for i in range(len(leak_p_all.index)): + # if dpressure_std.iloc [i] > -200: # 0.001: + if if_gy == 1: + act_dpressure[leak_p_all.index[i]] = ( + leak_p_all.iloc[i] - monitor_p_all.iloc[i] + ) / dpressure_std_all.iloc[i] + else: + act_dpressure[leak_p_all.index[i]] = ( + leak_p_all.iloc[i] - monitor_p_all.iloc[i] + ) + # part_dpressure = act_dpressure.loc[important_sensor_new] + + part_dpressure = ( + dpressure.loc[important_sensor_new] - dpressure_s.loc[important_sensor_new] + ) + similarity_pre_DIS = np.linalg.norm(part_dpressure) ## chang test + # part_dpressure = dpressure_s.loc[important_sensor]-dpressure.loc[important_sensor] + # similarity_pre_DIS = np.linalg.norm(part_dpressure) + # similarity_pre_DIS_later = 1 / (1 + similarity_pre_DIS) + similarity_dis = similarity_pre_DIS + + return similarity_cos, similarity_dis, none_flag + + +def adjust( + similarity_cos, + similarity_dis, + record_success_candidate, + record_success_no_candidate, +): + if len(record_success_no_candidate) > 0: + for each in record_success_no_candidate: + similarity_cos[each] = similarity_cos[record_success_candidate].min() * 0.9 + similarity_dis[each] = similarity_dis[record_success_candidate].max() * 1.1 + return similarity_cos, similarity_dis + + +def cal_sq_all_multi( + similarity_cos, + similarity_dis, + similarity_f, + candidate_pipe, + timestep_list_spc, + if_flow, + if_only_cos, + if_only_flow, + cos_h_input, + dis_h_input, + dis_f_h_input, + if_compalsive, + cos_sensor_num, + flow_sensor_num, +): + if if_only_flow == 1: + similarity_f, h_f = cal_sq_single_array( + similarity_f.values.reshape((-1, 1)), if_direct=2 + ) + sq_cos = 0 + sq_dis = 0 + sq_f = 1 + similarity_all = similarity_f * sq_f + output_similarity = similarity_all.reshape((-1, len(candidate_pipe))) + output_similarity_pd = pd.DataFrame( + output_similarity, index=timestep_list_spc, columns=candidate_pipe + ) + else: + if if_only_cos == 0: + if if_flow == 1: + # standerdize + similarity_cos, h_cos = cal_sq_single_array( + similarity_cos.values.reshape((-1, 1)), if_direct=1 + ) + similarity_dis, h_dis = cal_sq_single_array( + similarity_dis.values.reshape((-1, 1)), if_direct=2 + ) + similarity_f, h_f = cal_sq_single_array( + similarity_f.values.reshape((-1, 1)), if_direct=2 + ) + if if_compalsive == 1: + sq_cos = cos_h_input + sq_dis = dis_h_input + sq_f = dis_f_h_input + else: + """sq_cos = h_cos/(h_cos +h_dis +h_f ) + sq_dis = h_dis/(h_cos +h_dis +h_f ) + sq_f = h_f/(h_cos +h_dis +h_f )""" + sq_cos, sq_dis, sq_f = add_weight_for_SQ( + h_cos, h_dis, h_f, cos_sensor_num, flow_sensor_num + ) + + """if cos_sensor_num == 2 and sq_cos>0.2: + sq_cos = 0.2 + sq_dis = 0.8*h_dis / (h_dis + h_f) + sq_f = 0.8*h_f / (h_dis + h_f) + if cos_sensor_num == 1 and sq_dis > 0.3: + sq_cos = 0.1 + sq_dis = 0.3 + sq_f = 0.6""" + sq_cos, sq_dis, sq_f = adjust_ratio("CDF", sq_cos, sq_dis, sq_f) + if cos_sensor_num <= 1: + sq_cos = 0 + # similarity + + similarity_all = ( + similarity_cos * sq_cos + + similarity_dis * sq_dis + + similarity_f * sq_f + ) + output_similarity = similarity_all.reshape((-1, len(candidate_pipe))) + output_similarity_pd = pd.DataFrame( + output_similarity, index=timestep_list_spc, columns=candidate_pipe + ) + else: + # standerdize + similarity_cos, h_cos = cal_sq_single_array( + similarity_cos.values.reshape((-1, 1)), if_direct=1 + ) + similarity_dis, h_dis = cal_sq_single_array( + similarity_dis.values.reshape((-1, 1)), if_direct=2 + ) + if if_compalsive == 1: + sq_cos = cos_h_input + sq_dis = dis_h_input + else: + sq_cos = h_cos / (h_cos + h_dis) + sq_dis = h_dis / (h_cos + h_dis) + if cos_sensor_num == 2 and sq_cos > 0.5: + sq_cos = 0.5 + sq_dis = 0.5 + sq_cos, sq_dis, sq_f = adjust_ratio("CAD_new_gy", sq_cos, sq_dis, 0) + sq_f = 0 + # similarity + similarity_all = similarity_cos * sq_cos + similarity_dis * sq_dis + output_similarity = similarity_all.reshape((-1, len(candidate_pipe))) + output_similarity_pd = pd.DataFrame( + output_similarity, index=timestep_list_spc, columns=candidate_pipe + ) + elif if_only_cos == 1: + if if_flow == 1: + # standerdize + similarity_cos, h_cos = cal_sq_single_array( + similarity_cos.values.reshape((-1, 1)), if_direct=1 + ) + similarity_f, h_f = cal_sq_single_array( + similarity_f.values.reshape((-1, 1)), if_direct=2 + ) + if if_compalsive == 1: + sq_cos = cos_h_input + sq_f = dis_f_h_input + else: + sq_cos = h_cos / (h_cos + h_f) + sq_f = h_f / (h_cos + h_f) + sq_cos, sq_dis, sq_f = adjust_ratio("CAF", sq_cos, 0, sq_f) + sq_dis = 0 + # similarity + similarity_all = similarity_cos * sq_cos + similarity_f * sq_f + output_similarity = similarity_all.reshape((-1, len(candidate_pipe))) + output_similarity_pd = pd.DataFrame( + output_similarity, index=timestep_list_spc, columns=candidate_pipe + ) + + else: + sq_cos = cos_h_input + sq_dis = dis_h_input + sq_f = dis_f_h_input + output_similarity_pd = similarity_cos + else: + sq_cos = cos_h_input + sq_dis = dis_h_input + sq_f = dis_f_h_input + output_similarity_pd = 1 / (similarity_dis + 1) + return output_similarity_pd, sq_cos, sq_dis, sq_f + + +def add_weight_for_SQ(h_cos, h_dis, h_f, sensor_cos_num, sensor_f_num): + h_f_new = h_f * sensor_f_num + if sensor_cos_num <= 1: + h_cos_new = 0 + h_dis_new = h_dis * sensor_cos_num + else: + h_cos_new = h_cos * sensor_cos_num # / 2 + h_dis_new = h_dis * sensor_cos_num # / 2 + cos_sq = h_cos_new / (h_cos_new + h_dis_new + h_f_new) + dis_sq = h_dis_new / (h_cos_new + h_dis_new + h_f_new) + f_sq = h_f_new / (h_cos_new + h_dis_new + h_f_new) + + if sensor_cos_num == 2 and cos_sq > 0.2: + cos_sq = 0.2 + dis_sq = 0.8 * h_dis_new / (h_dis_new + h_f_new) + f_sq = 0.8 * h_f_new / (h_dis_new + h_f_new) + """if sensor_cos_num == 1: + if dis_sq / f_sq > sensor_cos_num/sensor_f_num: + dis_sq = sensor_cos_num/sensor_f_num + f_sq=1-dis_sq""" + # if h_dis_new/h_f_new > sensor_cos_num/sensor_f_num + return cos_sq, dis_sq, f_sq + + +def cal_sq_single_array(similarity_pre, if_direct): + if similarity_pre.max() - similarity_pre.min() == 0: + similarity_pre = np.ones(similarity_pre.shape) + else: + if if_direct == 1: + similarity_pre = ( + 0.998 + * (similarity_pre - similarity_pre.min()) + / (similarity_pre.max() - similarity_pre.min()) + + 0.002 + ) + else: + similarity_pre = ( + 0.998 + * (similarity_pre.max() - similarity_pre) + / (similarity_pre.max() - similarity_pre.min()) + + 0.002 + ) + # calculate pij + similarity_p = similarity_pre / similarity_pre.sum() + # cal xinxishang + similarity_lnp = np.zeros((len(similarity_pre), 1)) + for j in range(len(similarity_p)): + similarity_lnp[j] = -similarity_p[j] * math.log(similarity_p[j], math.e) + h = 1 - 1 / math.log(len(similarity_pre), math.e) * similarity_lnp.sum() + return similarity_pre, h + + +def cal_similarity_all_multi_new_sq_improve_double_lzr( + candidate_pipe, + similarity_mode, + pressure_leak, + monitor_p, + predict_p, + normal_p, + if_flow, + if_only_cos, + if_only_flow, + flow_leak, + monitor_f, + predict_f, + normal_f, + timestep_list, + Top_sensor_num, + if_gy, + effective_sensor, + cos_h, + dis_h, + dis_f_h, + if_compalsive, + max_flow, +): + similarity = pd.Series(dtype=float, index=candidate_pipe) + important_p_sensor = cal_top_sensors(monitor_p, predict_p, Top_sensor_num) + # important_f_sensor, basic_f = cal_top_f_sensor(normal_f) + important_f_sensor = monitor_f.columns + if ( + len(important_p_sensor) > 0 or len(important_f_sensor) > 0 + ): # if len(important_p_sensor) > 0 + break_flag = 0 + pressure_leak_new = pressure_leak.swaplevel() + # flow_leak_new = flow_leak.swaplevel() + if isinstance(flow_leak, pd.DataFrame) and len(flow_leak) > 0: + flow_leak_new = flow_leak.swaplevel() + else: + flow_leak_new = None + total_similarity_cos = pd.DataFrame(index=timestep_list, columns=candidate_pipe) + total_similarity_dis = pd.DataFrame(index=timestep_list, columns=candidate_pipe) + total_similarity_dis_f = pd.DataFrame( + index=timestep_list, columns=candidate_pipe + ) + + for timestep in timestep_list: + # cal p_cos, p_dis, f_dis + if if_only_flow != 1: + pressure_leak_temp = pressure_leak_new.loc[timestep].loc[ + :, effective_sensor + ] + monitor_p_temp = monitor_p.loc[timestep, effective_sensor] + predict_p_temp = predict_p.loc[timestep, effective_sensor] + normal_p_temp = normal_p.loc[timestep, effective_sensor] + + ( + total_similarity_cos.loc[timestep, :], + total_similarity_dis.loc[timestep, :], + ) = cal_similarity_all_cos_dis( + candidate_pipe, + pressure_leak_temp, + similarity_mode, + monitor_p_temp, + predict_p_temp, + normal_p_temp, + pressure_leak_new.loc[timestep].loc[:, monitor_p.columns], + monitor_p.loc[timestep, :], + predict_p.loc[timestep, :], + normal_p.loc[timestep, :], + important_p_sensor, + if_gy, + cos_or_flow=1, + ) + + if if_flow == 1: + if len(timestep_list) == 1: + leak_f_temp = flow_leak_new.loc[timestep].loc[:, important_f_sensor] + monitor_f_temp = monitor_f.loc[timestep, important_f_sensor] + predict_f_temp = predict_f.loc[timestep, important_f_sensor] + normal_f_temp = normal_f.loc[timestep, important_f_sensor] + basic_normal_f_temp = abs(max_flow.loc[important_f_sensor]) + + leak_f_temp = leak_f_temp / basic_normal_f_temp + monitor_f_temp = monitor_f_temp / basic_normal_f_temp + predict_f_temp = predict_f_temp / basic_normal_f_temp + normal_f_temp = normal_f_temp / basic_normal_f_temp + + else: + basic_f = abs(max_flow.loc[important_f_sensor]) + leak_f_temp = ( + flow_leak_new.loc[timestep].loc[:, important_f_sensor] / basic_f + ) + monitor_f_temp = ( + monitor_f.loc[timestep, important_f_sensor] / basic_f + ) + predict_f_temp = ( + predict_f.loc[timestep, important_f_sensor] / basic_f + ) + normal_f_temp = normal_f.loc[timestep, important_f_sensor] / basic_f + _, total_similarity_dis_f.loc[timestep, :] = cal_similarity_all_cos_dis( + candidate_pipe, + leak_f_temp, + similarity_mode, + monitor_f_temp, + predict_f_temp, + normal_f_temp, + flow_leak_new.loc[timestep].loc[:, monitor_f.columns], + monitor_f.loc[timestep, :], + predict_f.loc[timestep, :], + normal_f.loc[timestep, :], + important_f_sensor, + if_gy, + cos_or_flow=2, + ) + else: + total_similarity_dis_f = [] + similarity_all, cos_h, dis_h, dis_f_h = cal_sq_all_multi( + total_similarity_cos, + total_similarity_dis, + total_similarity_dis_f, + candidate_pipe, + timestep_list, + if_flow, + if_only_cos, + if_only_flow, + cos_h, + dis_h, + dis_f_h, + if_compalsive, + len(important_p_sensor), + len(important_f_sensor), + ) + if len(timestep_list) == 1: + similarity = similarity_all.iloc[0] + elif len(timestep_list) > 3: + for each_candidate in candidate_pipe: + similarity[each_candidate] = remove_3_sigma( + similarity_all.loc[:, each_candidate] + ) + else: + for each_candidate in candidate_pipe: + similarity[each_candidate] = similarity_all.loc[ + :, each_candidate + ].mean() + similarity = similarity.sort_values(ascending=False) + else: + break_flag = 1 + similarity = 0 + cos_h = 0 + dis_h = 0 + dis_f_h = 0 + return similarity, cos_h, dis_h, dis_f_h, break_flag + + +def cal_similarity_all_cos_dis( + candidate_pipe, + pressure_leak, + similarity_mode, + monitor_p, + predict_p, + normal_p, + pressure_leak_all, + monitor_p_all, + predict_p_all, + normal_p_all, + important_sensor, + if_gy, + cos_or_flow, +): + similarity_cos = pd.Series(dtype=float, index=candidate_pipe) + similarity_dis = pd.Series(dtype=float, index=candidate_pipe) + dpressure = normal_p - pressure_leak + + # 无用 ---------------------------------------------- + mean_dpressure = dpressure.mean() + + monitor_new = pd.DataFrame(index=["monitor"], columns=monitor_p.index) + monitor_new.iloc[0] = monitor_p + add_m_leak_pressure = [pressure_leak, monitor_p] + add_m_leak_pressure = pd.concat(add_m_leak_pressure) + pressure_leak_std = np.std(add_m_leak_pressure, ddof=1) + pressure_leak_std = pd.Series(pressure_leak_std, index=pressure_leak.columns) + + add_m_leak_pressure_all = [pressure_leak_all, monitor_p_all] + add_m_leak_pressure_all = pd.concat(add_m_leak_pressure_all) + pressure_leak_std_all = np.std(add_m_leak_pressure_all, ddof=1) + pressure_leak_std_all = pd.Series( + pressure_leak_std_all, index=pressure_leak.columns + ) + # 无用 ---------------------------------------------- + + monitor_p_temp = monitor_p + predict_p_temp = predict_p + normal_p_temp = normal_p + monitor_p_temp_all = monitor_p_all + predict_p_temp_all = predict_p_all + normal_p_temp_all = normal_p_all + record_success_candidate = [] + record_success_no_candidate = [] + for i in range(len(candidate_pipe)): + leak_p = pressure_leak.iloc[i, :] + leak_p_all = pressure_leak_all.iloc[i, :] + similarity_cos.iloc[i], similarity_dis.iloc[i], none_flag = ( + cal_similarity_simple_return_dd( + similarity_mode, + monitor_p_temp, + predict_p_temp, + normal_p_temp, + leak_p, + monitor_p_temp_all, + predict_p_temp_all, + normal_p_temp_all, + leak_p_all, + important_sensor, + mean_dpressure, + pressure_leak_std, + pressure_leak_std_all, + if_gy, + cos_or_flow, + ) + ) + if none_flag == 0: + record_success_candidate.append(candidate_pipe[i]) + else: + record_success_no_candidate.append(candidate_pipe[i]) + similarity_cos, similarity_dis = adjust( + similarity_cos, + similarity_dis, + record_success_candidate, + record_success_no_candidate, + ) + return similarity_cos, similarity_dis + + +def cal_top_f_sensor(normal_f): + if type(normal_f) == pd.core.frame.DataFrame: + mean_f = normal_f.mean() + else: + mean_f = normal_f + output_sensor = [] + output_normal_f = pd.Series(dtype=object) + for i in range(len(mean_f.index)): + if abs(mean_f.iloc[i]) > 0.01 / 3600: + output_sensor.append(mean_f.index[i]) + output_normal_f[mean_f.index[i]] = mean_f.iloc[i] + return output_sensor, output_normal_f + + +def cal_top_sensors(monitor_p, predict_p, Top_sensor_num): + dpressure = abs(predict_p - monitor_p) + if type(dpressure) == pd.core.frame.DataFrame: + dpressure = dpressure.mean() + dpressure_rank = dpressure.sort_values(ascending=False) + return list(dpressure_rank.index[:Top_sensor_num]) + + +def remove_3_sigma(similarity_t): + all_sample = len(similarity_t.index) + apart_sample = math.ceil(all_sample * 0.6) + similarity = similarity_t.astype("float") + mean_t = similarity.mean() + std_t = similarity.std() + new_similarity = similarity[ + (similarity <= mean_t + 3 * std_t) & (similarity >= mean_t - 3 * std_t) + ] + mean_t_new = new_similarity.mean() + return mean_t_new + + +def update_similarity(leak_candidate_center, similarity, leak_center_dict): + similarity_new = pd.Series(dtype=object) + for each_center in leak_candidate_center: + houxuan_center = leak_center_dict[each_center] + if len(houxuan_center) > 1: + temp_similarity = similarity[houxuan_center] + similarity_new[each_center] = temp_similarity.max() + else: + if type(similarity[each_center]) == pd.core.series.Series: + similarity_new[each_center] = similarity[each_center].mean() + else: + similarity_new[each_center] = similarity[each_center] + similarity_new = similarity_new.sort_values(ascending=False) + return similarity_new + + +def extra_judge(similarity): + mean_similarity = float(similarity.mean()) + record_sensor = [] + record_value = [] + for i in range(len(similarity.index)): + if similarity.iloc[i] >= mean_similarity: + record_value.append(similarity.iloc[i]) + record_sensor.append(similarity.index[i]) + out_put_similarity = pd.Series(record_value, index=record_sensor, dtype=object) + cut_ratio = len(out_put_similarity.index) / len(similarity.index) + return cut_ratio, out_put_similarity + + +def adjust_ratio(similarity_mode, cos_h, dis_h, dis_f_h, low_limit=0.1): + if similarity_mode == "CAF": + if cos_h < low_limit: + cos_h = low_limit + dis_f_h = 1 - cos_h + elif dis_f_h < low_limit: + dis_f_h = low_limit + cos_h = 1 - dis_f_h + elif similarity_mode == "CAD_new_gy": + if dis_h < low_limit: + dis_h = low_limit + cos_h = 1 - dis_h + elif cos_h < low_limit: + cos_h = low_limit + dis_h = 1 - cos_h + elif similarity_mode == "CDF": + normal_index = [0, 1, 2] + h_list = [cos_h, dis_h, dis_f_h] + if cos_h < low_limit: + h_list[0] = low_limit + normal_index.remove(0) + if dis_h < low_limit: + h_list[1] = low_limit + normal_index.remove(1) + if dis_f_h < low_limit: + h_list[2] = low_limit + normal_index.remove(2) + + if len(normal_index) == 1: + h_list[normal_index[0]] = h_list[normal_index[0]] - (sum(h_list) - 1) + elif len(normal_index) == 2: + sum_list = sum(h_list) + multiper = 1 - (sum_list - 1) / ( + h_list[normal_index[0]] + h_list[normal_index[1]] + ) + h_list[normal_index[0]] = h_list[normal_index[0]] * multiper + h_list[normal_index[1]] = h_list[normal_index[1]] * multiper + + cos_h, dis_h, dis_f_h = h_list[0], h_list[1], h_list[2] + + return cos_h, dis_h, dis_f_h + + +def decode_mode(similarity_mode): + if similarity_mode == "COS": + if_flow = 0 + if_only_cos = 1 + if_only_flow = 0 + elif similarity_mode == "CAD_new_gy": + if_flow = 0 + if_only_cos = 0 + if_only_flow = 0 + elif similarity_mode == "CDF": + if_flow = 1 + if_only_cos = 0 + if_only_flow = 0 + elif similarity_mode == "CAF": + if_flow = 1 + if_only_cos = 1 + if_only_flow = 0 + elif similarity_mode == "DIS": + if_flow = 1 + if_only_cos = 2 + if_only_flow = 0 + elif similarity_mode == "OF": + if_flow = 1 + if_only_cos = 0 + if_only_flow = 1 + return if_flow, if_only_cos, if_only_flow + + +class SimilarityCalculator: + @staticmethod + def cal_similarity_simple_return_dd(*args, **kwargs): + return cal_similarity_simple_return_dd(*args, **kwargs) + + @staticmethod + def cal_similarity_all_cos_dis(*args, **kwargs): + return cal_similarity_all_cos_dis(*args, **kwargs) + + @staticmethod + def cal_similarity_all_multi_new_sq_improve_double_lzr(*args, **kwargs): + return cal_similarity_all_multi_new_sq_improve_double_lzr(*args, **kwargs) + + @staticmethod + def cal_sq_all_multi(*args, **kwargs): + return cal_sq_all_multi(*args, **kwargs) + + @staticmethod + def cal_sq_single_array(*args, **kwargs): + return cal_sq_single_array(*args, **kwargs) + + @staticmethod + def add_weight_for_SQ(*args, **kwargs): + return add_weight_for_SQ(*args, **kwargs) + + @staticmethod + def adjust_ratio(*args, **kwargs): + return adjust_ratio(*args, **kwargs) + + @staticmethod + def adjust(*args, **kwargs): + return adjust(*args, **kwargs) + + @staticmethod + def cal_top_sensors(*args, **kwargs): + return cal_top_sensors(*args, **kwargs) + + @staticmethod + def cal_top_f_sensor(*args, **kwargs): + return cal_top_f_sensor(*args, **kwargs) + + @staticmethod + def remove_3_sigma(*args, **kwargs): + return remove_3_sigma(*args, **kwargs) + + @staticmethod + def update_similarity(*args, **kwargs): + return update_similarity(*args, **kwargs) + + @staticmethod + def extra_judge(*args, **kwargs): + return extra_judge(*args, **kwargs) + + @staticmethod + def decode_mode(*args, **kwargs): + return decode_mode(*args, **kwargs) diff --git a/app/api/v1/endpoints/burst_location.py b/app/api/v1/endpoints/burst_location.py new file mode 100644 index 0000000..c3898cf --- /dev/null +++ b/app/api/v1/endpoints/burst_location.py @@ -0,0 +1,33 @@ +from typing import Any + +from fastapi import APIRouter, Depends, HTTPException +from pydantic import BaseModel + +from app.auth.keycloak_dependencies import get_current_keycloak_username +from app.services.burst_location import run_burst_location_by_network + +router = APIRouter() + + +class BurstLocationRequest(BaseModel): + network: str + pressure_scada_ids: list[str] + burst_pressure: dict[str, float] | list[dict[str, Any]] + normal_pressure: dict[str, float] | list[dict[str, Any]] + burst_leakage: float + flow_scada_ids: list[str] | None = None + burst_flow: dict[str, float] | list[dict[str, Any]] | None = None + normal_flow: dict[str, float] | list[dict[str, Any]] | None = None + min_dpressure: float = 2.0 + basic_pressure: float = 10.0 + + +@router.post("/locate/") +async def locate_burst( + data: BurstLocationRequest, + _username: str = Depends(get_current_keycloak_username), +) -> dict[str, Any]: + try: + return run_burst_location_by_network(**data.model_dump()) + except (TypeError, ValueError) as exc: + raise HTTPException(status_code=400, detail=str(exc)) diff --git a/app/api/v1/router.py b/app/api/v1/router.py index bc13d8d..42a0d0d 100644 --- a/app/api/v1/router.py +++ b/app/api/v1/router.py @@ -13,6 +13,7 @@ from app.api.v1.endpoints import ( risk, cache, leakage, + burst_location, user_management, # 新增：用户管理 audit, # 新增：审计日志 meta, @@ -85,6 +86,9 @@ api_router.include_router(misc.router, tags=["Misc"]) api_router.include_router(risk.router, tags=["Risk"]) api_router.include_router(cache.router, tags=["Cache"]) api_router.include_router(leakage.router, prefix="/leakage", tags=["Leakage"]) +api_router.include_router( + burst_location.router, prefix="/burst-location", tags=["Burst Location"] +) # Database Routers api_router.include_router(timescaledb_router, tags=["TimescaleDB"]) diff --git a/app/services/burst_location.py b/app/services/burst_location.py new file mode 100644 index 0000000..702f27a --- /dev/null +++ b/app/services/burst_location.py @@ -0,0 +1,83 @@ +from __future__ import annotations + +from pathlib import Path +from typing import Any +from uuid import uuid4 + +import pandas as pd + +from app.algorithms.burst_location import run_burst_location +from app.services.tjnetwork import dump_inp + +SeriesInput = pd.Series | dict[str, Any] | list[dict[str, Any]] + + +def _normalize_series(data: SeriesInput, field_name: str) -> pd.Series: + if isinstance(data, pd.Series): + series = data.copy() + elif isinstance(data, dict): + series = pd.Series(data, dtype=float) + elif isinstance(data, list): + if len(data) == 0: + return pd.Series(dtype=float) + frame = pd.DataFrame(data) + if not {"id", "value"}.issubset(frame.columns): + raise ValueError(f"{field_name} list item must include 'id' and 'value'.") + series = pd.Series( + frame["value"].values, index=frame["id"].astype(str).values, dtype=float + ) + else: + raise ValueError(f"Unsupported data format for {field_name}.") + + series.index = series.index.map(str) + return pd.to_numeric(series, errors="raise") + + +def run_burst_location_by_network( + *, + network: str, + pressure_scada_ids: list[str], + burst_pressure: SeriesInput, + normal_pressure: SeriesInput, + burst_leakage: float, + flow_scada_ids: list[str] | None = None, + burst_flow: SeriesInput | None = None, + normal_flow: SeriesInput | None = None, + min_dpressure: float = 2.0, + basic_pressure: float = 10.0, +) -> dict[str, Any]: + if not network: + raise ValueError("network is required.") + + tmp_filename = f"burst_location_{network}_{uuid4().hex}.inp" + inp_path = Path.cwd() / tmp_filename + + try: + dump_inp(network, tmp_filename) + + burst_pressure_series = _normalize_series(burst_pressure, "burst_pressure") + normal_pressure_series = _normalize_series(normal_pressure, "normal_pressure") + burst_flow_series = ( + _normalize_series(burst_flow, "burst_flow") if burst_flow is not None else None + ) + normal_flow_series = ( + _normalize_series(normal_flow, "normal_flow") + if normal_flow is not None + else None + ) + + return run_burst_location( + wn_inp_path=str(inp_path), + pressure_scada_ids=pressure_scada_ids, + burst_pressure=burst_pressure_series, + normal_pressure=normal_pressure_series, + burst_leakage=burst_leakage, + flow_scada_ids=flow_scada_ids, + burst_flow=burst_flow_series, + normal_flow=normal_flow_series, + min_dpressure=min_dpressure, + basic_pressure=basic_pressure, + ) + finally: + if inp_path.exists(): + inp_path.unlink()