diff --git a/app/algorithms/sensor/sensitivity.py b/app/algorithms/sensor/sensitivity.py index ef09fad..618c08a 100644 --- a/app/algorithms/sensor/sensitivity.py +++ b/app/algorithms/sensor/sensitivity.py @@ -1,701 +1,916 @@ -# 改进灵敏度法 -import networkx +"""Pressure sensor placement based on scalable sensitivity analysis. + +The original implementation expanded a sparse water network into several dense +``node x node``, ``node x pipe``, and ``pipe x pipe`` matrices. That made the +memory requirement quadratic and the explicit matrix inverse cubic in time. + +This module keeps one algorithm for every network size: + +* run EPANET once and reuse the first hydraulic state; +* keep incidence and hydraulic graphs sparse; +* estimate the row-wise L1 pressure sensitivity with deterministic Cauchy + projections and one sparse factorization; +* estimate total directed hydraulic distance from a deterministic spatial + coreset without materialising an all-pairs distance matrix; +* balance sensitivity score with geographic and pipe-network coverage without + materialising candidate-to-candidate distances. + +The random seed and sample counts are fixed, so the same model and request +produce the same placement on every run. +""" + +from __future__ import annotations + +import logging +from dataclasses import dataclass +from pathlib import Path +from tempfile import TemporaryDirectory +from time import perf_counter + import numpy as np -import pandas import wntr -import pandas as pd -import copy -import matplotlib.pyplot as plt -import networkx as nx -from sklearn.cluster import KMeans -from wntr.epanet.toolkit import EpanetException -from numpy.linalg import slogdet -import random -from matplotlib.lines import Line2D -from sklearn.cluster import SpectralClustering -import libpysal as ps -from spopt.region import Skater -from shapely.geometry import Point -import geopandas as gpd -from sklearn.metrics import pairwise_distances -import app.services.project_info as project_info +from scipy.sparse import csr_matrix, eye +from scipy.sparse.csgraph import connected_components, dijkstra +from scipy.sparse.linalg import splu +from sklearn.cluster import MiniBatchKMeans -# 2025/03/12 -# Step1: 获取节点坐标 -def getCoor(wn: wntr.network.WaterNetworkModel) -> pandas.DataFrame: +logger = logging.getLogger(__name__) + +_RANDOM_SEED = 42 +_SENSITIVITY_PROJECTIONS = 256 +_HYDRAULIC_LANDMARKS = 256 +_PROJECTION_BLOCK_SIZE = 16 +_DIJKSTRA_BLOCK_SIZE = 16 +_HEADLOSS_EPSILON = 1e-10 +_DIAMETER_TOLERANCE_MM = 1e-9 +_COVERAGE_ELIGIBILITY_RATIO = 0.70 +_COVERAGE_EDGE_EPSILON = 1e-9 + + +@dataclass(frozen=True) +class _PreparedNetwork: + """Sparse data required by the placement pipeline.""" + + node_names: tuple[str, ...] + full_node_indices: np.ndarray + candidate_indices: np.ndarray + coordinates: np.ndarray + incidence: csr_matrix + conductance: np.ndarray + roughness_response: np.ndarray + distance_graph: csr_matrix + coverage_graph: csr_matrix + + +@dataclass(frozen=True) +class _CandidatePool: + """Aligned candidate arrays consumed by the placement stage.""" + + full_indices: np.ndarray + coordinates: np.ndarray + names: np.ndarray + scores: np.ndarray + + +def _run_hydraulic_simulation( + wn: wntr.network.WaterNetworkModel, +): + """Run only the initial EPANET state without shared ``temp.*`` files.""" + + original_duration = wn.options.time.duration + try: + # Every downstream calculation reads ``iloc[0]``. Running an extended + # simulation only allocates unused time-series results, which is + # especially expensive for daily models with tens of thousands of + # nodes. Restore the caller's model even when EPANET fails. + wn.options.time.duration = 0 + with TemporaryDirectory(prefix="tjwater-sensitivity-") as temp_dir: + file_prefix = str(Path(temp_dir) / "simulation") + return wntr.sim.EpanetSimulator(wn).run_sim(file_prefix=file_prefix) + finally: + wn.options.time.duration = original_duration + + +def _excluded_elements( + wn: wntr.network.WaterNetworkModel, +) -> tuple[set[str], set[str]]: + """Return nodes that cannot host sensors and source-connected pipes. + + Reservoirs, tanks, pump/valve endpoints, and the junction immediately next + to a reservoir or tank are treated as hydraulic boundary nodes. Pipes + connected directly to a source are removed from the perturbation set, as + in the legacy algorithm. """ - 获取管网模型的节点坐标 - :param wn: 由wntr生成的模型 - :return: 节点坐标 - """ - # site: pandas.Series - # index：节点名称（wn.node_name_list） - # values：每个节点的坐标，格式为 tuple（如 (x, y) 或 (x, y, z)） - site = wn.query_node_attribute("coordinates") - # Coor: pandas.Series - # index：与site相同（节点名称）。 - # values：坐标转换为numpy.ndarray（如array([10.5, 20.3])） - Coor = site.apply(lambda x: np.array(x)) # 将节点坐标转换为numpy数组 - # x, y: list[float] - x = [] # 存储所有节点的 x 坐标 - y = [] # 存储所有节点的 y 坐标 - for i in range(0, len(Coor)): - x.append(Coor.values[i][0]) # 将 x 坐标存入 x 列表。 - y.append(Coor.values[i][1]) # 将 y 坐标存入 y 列表 - # xy: dict[str, list], x、y 坐标的字典 - xy = {"x": x, "y": y} - # Coor_node: pandas.DataFrame, 存储节点 x, y 坐标的 DataFrame - Coor_node = pd.DataFrame(xy, index=wn.node_name_list, columns=["x", "y"]) - return Coor_node + + source_nodes = set(wn.reservoir_name_list) | set(wn.tank_name_list) + excluded_nodes = set(source_nodes) + source_pipes: set[str] = set() + + for pipe_name, pipe in wn.pipes(): + endpoints = {pipe.start_node_name, pipe.end_node_name} + if endpoints & source_nodes: + source_pipes.add(pipe_name) + excluded_nodes.update(endpoints) + + for _link_name, link in list(wn.pumps()) + list(wn.valves()): + excluded_nodes.add(link.start_node_name) + excluded_nodes.add(link.end_node_name) + + return excluded_nodes, source_pipes -# 2025/03/12 -# Step2: KMeans 聚类 -# 将节点用kmeans根据坐标分为k组，存入字典g -def kgroup(coor: pandas.DataFrame, knum: int) -> dict[int, list[str]]: - """ - 使用KMeans聚类，将节点坐标分组 - :param coor: 存储所有节点的坐标数据 - :param knum: 需要分成的聚类数 - :return: 聚类结果字典 - """ - g = {} - # estimator: sklearn.cluster.KMeans,KMeans 聚类模型 - estimator = KMeans(n_clusters=knum) - estimator.fit(coor) - # label_pred: numpy.ndarray（int）,每个点的类别标签 - label_pred = estimator.labels_ - for i in range(0, knum): - g[i] = coor[label_pred == i].index.tolist() - return g +def _minimum_weight_csr( + rows: list[int], + columns: list[int], + weights: list[float], + *, + shape: tuple[int, int], +) -> csr_matrix: + """Build a CSR graph while retaining the lightest parallel edge.""" + if not rows: + return csr_matrix(shape, dtype=np.float64) -def skater_partition(G, n_clusters): - """ - 使用 SKATER 算法对输入的无向图 G 进行区域划分， - 保证每个划分区域在图论意义上是连通的， - 同时依据节点坐标的空间信息进行划分。 + row_array = np.asarray(rows, dtype=np.int64) + column_array = np.asarray(columns, dtype=np.int64) + weight_array = np.asarray(weights, dtype=np.float64) + order = np.lexsort((column_array, row_array)) + row_array = row_array[order] + column_array = column_array[order] + weight_array = weight_array[order] - 参数: - G: networkx.Graph - 带有节点坐标属性（键为 'pos'）的无向图。 - n_clusters: int - 希望划分的区域数量。 - - 返回: - groups: dict - 字典形式的聚类结果，键为区域编号，值为该区域内的节点列表。 - """ - # 1. 获取所有节点坐标，假设每个节点都有 'pos' 属性 - pos = nx.get_node_attributes(G, "pos") - nodes = list(G.nodes()) - # 构造坐标数组：每行为 [x, y] - coords = np.array([pos[node] for node in nodes]) - - # 2. 构造 GeoDataFrame：创建 DataFrame 并生成 geometry 列 - df = pd.DataFrame(coords, columns=["x", "y"], index=nodes) - # 利用 shapely 的 Point 构造空间位置 - df["geometry"] = df.apply(lambda row: Point(row["x"], row["y"]), axis=1) - gdf = gpd.GeoDataFrame(df, geometry="geometry") - - # 3. 构造空间权重矩阵，使用 4 近邻方法（k=4，可根据实际情况调整） - w = ps.weights.KNN.from_array(coords, k=4) - w.transform = "R" - - # 4. 调用 SKATER：新版本 API 要求传入 gdf, w 以及 attrs_name（这里使用 'x' 和 'y' 作为属性） - skater = Skater(gdf, w, attrs_name=["x", "y"], n_clusters=n_clusters) - skater.solve() - - # 5. 获取聚类标签，构造成字典格式 - labels = skater.labels_ - groups = {} - for label, node in zip(labels, nodes): - groups.setdefault(label, []).append(node) - - return groups - - -def spectral_partition(G, n_clusters): - """ - 利用谱聚类算法对图 G 进行分区： - 1. 根据所有节点的空间坐标计算欧氏距离矩阵； - 2. 利用高斯核函数构造相似度矩阵； - 3. 使用 SpectralClustering 进行归一化割，返回分区结果。 - - 参数: - G: networkx.Graph - 每个节点需要有 'pos' 属性，其值为 (x, y) 坐标。 - n_clusters: int - 希望划分的聚类数目。 - - 返回: - groups: dict - 键为聚类标签，值为该聚类对应的节点列表。 - """ - # 1. 获取节点空间坐标，注意保证每个节点都有 'pos' 属性 - pos_dict = nx.get_node_attributes(G, "pos") - nodes = list(G.nodes()) - coords = np.array([pos_dict[node] for node in nodes]) - - # 2. 计算节点之间的欧氏距离矩阵 - D = pairwise_distances(coords, metric="euclidean") - - # 3. 计算 sigma 值：这里取所有距离的均值，当然也可以根据实际情况调整 - sigma = np.mean(D) - - # 4. 构造相似度矩阵：使用高斯核函数 - # A(i, j) = exp( -d(i,j)^2 / (2*sigma^2) ) - A = np.exp(-(D**2) / (2 * sigma**2)) - - # 5. 使用谱聚类进行图分区 - clustering = SpectralClustering( - n_clusters=n_clusters, affinity="precomputed", random_state=0 + group_start = np.empty(len(row_array), dtype=bool) + group_start[0] = True + group_start[1:] = (row_array[1:] != row_array[:-1]) | ( + column_array[1:] != column_array[:-1] + ) + starts = np.flatnonzero(group_start) + minimum_weights = np.minimum.reduceat(weight_array, starts) + return csr_matrix( + (minimum_weights, (row_array[starts], column_array[starts])), + shape=shape, ) - labels = clustering.fit_predict(A) - - # 6. 构造字典形式的分区结果 - groups = {} - for label, node in zip(labels, nodes): - groups.setdefault(label, []).append(node) - - return groups -# 2025/03/12 -# Step3: wn_func类，水力计算 -# wn_func 主要用于计算： -# 水力距离（hydraulic length）：即节点之间的水力阻力。 -# 灵敏度分析（sensitivity analysis）：用于优化测压点的布置。 -# 一些与水力相关的函数，包括 CtoS：求水力距离,stafun：求状态函数F -# # diff：求F对P的导数，返回灵敏度矩阵A -# # sensitivity：返回灵敏度和总灵敏度 -class wn_func(object): +def _node_coordinates( + wn: wntr.network.WaterNetworkModel, + node_names: tuple[str, ...], +) -> np.ndarray: + coordinate_series = wn.query_node_attribute("coordinates") + coordinates = np.asarray( + [coordinate_series.loc[node_name] for node_name in node_names], + dtype=np.float64, + ) + if coordinates.ndim != 2 or coordinates.shape[1] < 2: + raise ValueError("管网节点缺少二维坐标，无法进行监测点空间布置") + coordinates = coordinates[:, :2] + if not np.isfinite(coordinates).all(): + raise ValueError("管网节点坐标包含非有限值，无法进行监测点空间布置") + return coordinates - # Step3.1: 初始化 - def __init__(self, wn: wntr.network.WaterNetworkModel, min_diameter: int): - """ - 获取管网模型信息 - :param wn: 由wntr生成的模型 - :param min_diameter: 安装的最小管径 - """ - # self.results: wntr.sim.results.SimulationResults,仿真结果，包含压力、流量、水头等数据 - self.results = wntr.sim.EpanetSimulator(wn).run_sim() # 存储运行结果 - self.wn = wn - # self.q：pandas.DataFrame,管道流量，索引为时间步长，列为管道名称 - self.q = self.results.link["flowrate"] - # ReservoirIndex / Tankindex: list[str],水库 / 水箱节点名称列表 - ReservoirIndex = wn.reservoir_name_list - Tankindex = wn.tank_name_list - # 删除水库节点，删除与直接水库相连的虚拟管道 - # self.pipes: list[str],所有管道的名称 - self.pipes = wn.pipe_name_list - # self.nodes: list[str],所有节点的名称 - self.nodes = wn.node_name_list - # self.coordinates：pandas.Series,节点坐标，索引为节点名，值为 (x, y) 坐标的 tuple - self.coordinates = wn.query_node_attribute("coordinates") - # allpumps / allvalves: list[str],所有泵/阀门名称列表 - allpumps = wn.pump_name_list - allvalves = wn.valve_name_list - # pumpstnode / pumpednode / valvestnode / valveednode: list[str],存储泵和阀门 起终点节点的名称 - pumpstnode = [] - pumpednode = [] - valvestnode = [] - valveednode = [] - # Reservoirpipe / Reservoirednode: list[str],记录与水库相关的管道和节点 - Reservoirpipe = [] - Reservoirednode = [] - for pump in allpumps: - pumpstnode.append(wn.links[pump].start_node.name) - pumpednode.append(wn.links[pump].end_node.name) - for valve in allvalves: - valvestnode.append(wn.links[valve].start_node.name) - valveednode.append(wn.links[valve].end_node.name) - for pipe in self.pipes: - if wn.links[pipe].start_node.name in ReservoirIndex: - Reservoirpipe.append(pipe) - Reservoirednode.append(wn.links[pipe].end_node.name) - if wn.links[pipe].start_node.name in Tankindex: - Reservoirpipe.append(pipe) - Reservoirednode.append(wn.links[pipe].end_node.name) - if wn.links[pipe].end_node.name in Tankindex: - Reservoirpipe.append(pipe) - Reservoirednode.append(wn.links[pipe].start_node.name) - # 泵的起终点、tank、reservoir - # self.delnodes: list[str],需要删除的节点（包括水库、泵、阀门连接的节点） - self.delnodes = list( - set(ReservoirIndex).union( - Tankindex, - pumpstnode, - pumpednode, - valvestnode, - valveednode, - Reservoirednode, + +def _build_coverage_graph( + wn: wntr.network.WaterNetworkModel, + results, + full_node_index: dict[str, int], +) -> csr_matrix: + """Build the active undirected physical graph used to spread sensors.""" + + status_series = results.link["status"].iloc[0] + rows: list[int] = [] + columns: list[int] = [] + weights: list[float] = [] + + for link_name, link in wn.links(): + if float(status_series.loc[link_name]) <= 0: + continue + + start = full_node_index[link.start_node_name] + end = full_node_index[link.end_node_name] + # Pipes carry their physical length. Pumps and valves are point + # devices, so a tiny positive length preserves connectivity without + # dominating shortest-path distance. + weight = max( + float(getattr(link, "length", 0.0)), + _COVERAGE_EDGE_EPSILON, + ) + rows.extend((start, end)) + columns.extend((end, start)) + weights.extend((weight, weight)) + + return _minimum_weight_csr( + rows, + columns, + weights, + shape=(len(full_node_index), len(full_node_index)), + ) + + +def _prepare_network( + wn: wntr.network.WaterNetworkModel, + results, + *, + min_diameter: int, +) -> _PreparedNetwork: + excluded_nodes, source_pipes = _excluded_elements(wn) + full_node_names = tuple(wn.node_name_list) + full_node_index = { + node_name: index for index, node_name in enumerate(full_node_names) + } + node_names = tuple( + node_name for node_name in full_node_names if node_name not in excluded_nodes + ) + if not node_names: + raise ValueError("管网中没有可参与灵敏度分析的节点") + + node_index = {node_name: index for index, node_name in enumerate(node_names)} + full_node_indices = np.asarray( + [full_node_index[node_name] for node_name in node_names], + dtype=np.int64, + ) + coordinates = _node_coordinates(wn, node_names) + + flow_series = results.link["flowrate"].iloc[0] + headloss_series = results.link["headloss"].iloc[0] + head_series = results.node["head"].iloc[0] + + candidate_nodes: set[str] = set() + for _pipe_name, pipe in wn.pipes(): + diameter_mm = float(pipe.diameter) * 1000.0 + if diameter_mm + _DIAMETER_TOLERANCE_MM < min_diameter: + continue + if pipe.start_node_name in node_index: + candidate_nodes.add(pipe.start_node_name) + if pipe.end_node_name in node_index: + candidate_nodes.add(pipe.end_node_name) + + incidence_rows: list[int] = [] + incidence_columns: list[int] = [] + incidence_values: list[float] = [] + conductance: list[float] = [] + roughness_response: list[float] = [] + distance_rows: list[int] = [] + distance_columns: list[int] = [] + distance_weights: list[float] = [] + + kept_pipe_count = 0 + for pipe_name, pipe in wn.pipes(): + if pipe_name in source_pipes: + continue + + start_name = pipe.start_node_name + end_name = pipe.end_node_name + if start_name not in node_index and end_name not in node_index: + continue + + flow = float(flow_series.loc[pipe_name]) + absolute_flow = abs(flow) + headloss = abs(float(headloss_series.loc[pipe_name])) + roughness = float(pipe.roughness) + if roughness <= 0: + raise ValueError(f"管道 {pipe_name} 的粗糙度必须大于 0") + + orientation = -1.0 if flow < 0 else 1.0 + if start_name in node_index: + incidence_rows.append(node_index[start_name]) + incidence_columns.append(kept_pipe_count) + incidence_values.append(-orientation) + if end_name in node_index: + incidence_rows.append(node_index[end_name]) + incidence_columns.append(kept_pipe_count) + incidence_values.append(orientation) + + conductance.append( + absolute_flow / (1.852 * headloss + _HEADLOSS_EPSILON) + ) + roughness_response.append(absolute_flow / roughness) + + if flow > 0: + upstream_name, downstream_name = start_name, end_name + else: + upstream_name, downstream_name = end_name, start_name + hydraulic_weight = ( + abs(float(head_series.loc[start_name]) - float(head_series.loc[end_name])) + * float(pipe.length) + ) + distance_rows.append(full_node_index[upstream_name]) + distance_columns.append(full_node_index[downstream_name]) + distance_weights.append(hydraulic_weight) + kept_pipe_count += 1 + + if kept_pipe_count == 0: + raise ValueError("管网中没有可用于灵敏度分析的管道") + + incidence = csr_matrix( + ( + np.asarray(incidence_values, dtype=np.float64), + ( + np.asarray(incidence_rows, dtype=np.int64), + np.asarray(incidence_columns, dtype=np.int64), + ), + ), + shape=(len(node_names), kept_pipe_count), + ) + conductance_array = np.asarray(conductance, dtype=np.float64) + response_array = np.asarray(roughness_response, dtype=np.float64) + if not np.isfinite(conductance_array).all() or not np.isfinite( + response_array + ).all(): + raise ValueError("水力结果产生了非有限灵敏度系数") + + distance_graph = _minimum_weight_csr( + distance_rows, + distance_columns, + distance_weights, + shape=(len(full_node_names), len(full_node_names)), + ) + coverage_graph = _build_coverage_graph(wn, results, full_node_index) + candidate_indices = np.asarray( + [ + index + for index, node_name in enumerate(node_names) + if node_name in candidate_nodes + ], + dtype=np.int64, + ) + + return _PreparedNetwork( + node_names=node_names, + full_node_indices=full_node_indices, + candidate_indices=candidate_indices, + coordinates=coordinates, + incidence=incidence, + conductance=conductance_array, + roughness_response=response_array, + distance_graph=distance_graph, + coverage_graph=coverage_graph, + ) + + +def _axis_normalized_coordinates(coordinates: np.ndarray) -> np.ndarray: + """Scale each axis independently for MiniBatchKMeans.""" + + minimum = coordinates.min(axis=0) + span = np.ptp(coordinates, axis=0) + span[span == 0] = 1.0 + return (coordinates - minimum) / span + + +def _isotropic_coordinates(coordinates: np.ndarray) -> np.ndarray: + """Normalize coordinates without distorting the network aspect ratio.""" + + minimum = coordinates.min(axis=0) + scale = float(np.max(np.ptp(coordinates, axis=0), initial=0.0)) + if scale == 0: + scale = 1.0 + return (coordinates - minimum) / scale + + +def _cluster_labels( + coordinates: np.ndarray, + cluster_count: int, + *, + random_seed: int, +) -> tuple[np.ndarray, np.ndarray]: + """Cluster coordinates deterministically with one implementation at all sizes.""" + + normalized = _axis_normalized_coordinates(coordinates) + if cluster_count == 1: + return np.zeros(len(coordinates), dtype=np.int64), normalized[[0]] + if cluster_count >= len(coordinates): + return np.arange(len(coordinates), dtype=np.int64), normalized.copy() + + model = MiniBatchKMeans( + n_clusters=cluster_count, + random_state=random_seed, + n_init=3, + batch_size=min(len(coordinates), max(1024, cluster_count * 3)), + max_iter=100, + max_no_improvement=20, + reassignment_ratio=0.0, + ) + labels = model.fit_predict(normalized).astype(np.int64, copy=False) + return labels, np.asarray(model.cluster_centers_, dtype=np.float64) + + +def _estimate_log_pressure_sensitivity(prepared: _PreparedNetwork) -> np.ndarray: + """Estimate each row's L1 sensitivity using streaming Cauchy projections.""" + + weighted_incidence = prepared.incidence.multiply(prepared.conductance) + laplacian = (weighted_incidence @ prepared.incidence.T).tocsc() + diagonal = np.asarray(laplacian.diagonal(), dtype=np.float64) + diagonal_scale = float(np.max(np.abs(diagonal), initial=0.0)) + if diagonal_scale == 0: + raise ValueError("水力雅可比矩阵为空，无法计算压力灵敏度") + + regularization = diagonal_scale * np.sqrt(np.finfo(np.float64).eps) + laplacian = laplacian + eye( + laplacian.shape[0], format="csc", dtype=np.float64 + ) * regularization + factor = splu( + laplacian, + permc_spec="MMD_AT_PLUS_A", + diag_pivot_thresh=0.0, + options={"SymmetricMode": True}, + ) + + random = np.random.default_rng(_RANDOM_SEED) + log_absolute_sum = np.zeros(len(prepared.node_names), dtype=np.float64) + projection_count = 0 + float_epsilon = np.finfo(np.float64).eps + float_tiny = np.finfo(np.float64).tiny + + while projection_count < _SENSITIVITY_PROJECTIONS: + block_size = min( + _PROJECTION_BLOCK_SIZE, + _SENSITIVITY_PROJECTIONS - projection_count, + ) + uniform = random.random((prepared.incidence.shape[1], block_size)) + np.clip(uniform, float_epsilon, 1.0 - float_epsilon, out=uniform) + cauchy_projection = np.tan(np.pi * (uniform - 0.5)) + projected_response = prepared.incidence @ ( + prepared.roughness_response[:, None] * cauchy_projection + ) + solution = factor.solve(np.asarray(projected_response, dtype=np.float64)) + log_absolute_sum += np.log( + np.maximum(np.abs(solution), float_tiny) + ).sum(axis=1) + projection_count += block_size + + # For a standard Cauchy variable E[log(abs(X))] is zero. Therefore this + # streaming geometric mean estimates log(||row||_1) without retaining the + # node-by-projection matrix. A finite-sample bias is common to all rows and + # does not affect ranking. + return log_absolute_sum / _SENSITIVITY_PROJECTIONS + + +def _landmark_coreset(prepared: _PreparedNetwork) -> tuple[np.ndarray, np.ndarray]: + landmark_count = min(_HYDRAULIC_LANDMARKS, len(prepared.node_names)) + labels, centers = _cluster_labels( + prepared.coordinates, + landmark_count, + random_seed=_RANDOM_SEED + 1, + ) + normalized = _axis_normalized_coordinates(prepared.coordinates) + landmarks: list[int] = [] + weights: list[float] = [] + + for label in np.unique(labels): + members = np.flatnonzero(labels == label) + center = centers[int(label)] + squared_distance = np.square(normalized[members] - center).sum(axis=1) + landmarks.append(int(members[int(np.argmin(squared_distance))])) + weights.append(float(len(members))) + + return ( + np.asarray(landmarks, dtype=np.int64), + np.asarray(weights, dtype=np.float64), + ) + + +def _estimate_hydraulic_distance_sums(prepared: _PreparedNetwork) -> np.ndarray: + """Estimate outbound distance sums without an all-pairs distance matrix.""" + + landmark_indices, landmark_weights = _landmark_coreset(prepared) + full_landmark_indices = prepared.full_node_indices[landmark_indices] + reversed_graph = prepared.distance_graph.transpose().tocsr() + distance_sums = np.zeros(len(prepared.node_names), dtype=np.float64) + + for start in range(0, len(landmark_indices), _DIJKSTRA_BLOCK_SIZE): + stop = min(start + _DIJKSTRA_BLOCK_SIZE, len(landmark_indices)) + distances = dijkstra( + reversed_graph, + directed=True, + indices=full_landmark_indices[start:stop], + return_predecessors=False, + ) + distances = np.atleast_2d(distances)[:, prepared.full_node_indices] + # The legacy matrix represented unreachable pairs as zero. Retaining + # that convention prevents disconnected branches from receiving an + # artificial infinite score. + distances[~np.isfinite(distances)] = 0.0 + distance_sums += landmark_weights[start:stop] @ distances + + return distance_sums + + +def _build_candidate_pool( + prepared: _PreparedNetwork, + log_sensitivity: np.ndarray, + hydraulic_distance_sums: np.ndarray, +) -> _CandidatePool: + candidate_indices = prepared.candidate_indices + candidate_distance = hydraulic_distance_sums[candidate_indices] + with np.errstate(divide="ignore", invalid="ignore"): + scores = log_sensitivity[candidate_indices] + np.log(candidate_distance) + scores = np.nan_to_num( + scores, + nan=-np.inf, + neginf=-np.inf, + posinf=np.finfo(np.float64).max, + ) + return _CandidatePool( + full_indices=prepared.full_node_indices[candidate_indices], + coordinates=_isotropic_coordinates(prepared.coordinates[candidate_indices]), + names=np.asarray( + [prepared.node_names[index] for index in candidate_indices], + dtype=str, + ), + scores=scores, + ) + + +def _highest_scoring_position( + names: np.ndarray, + scores: np.ndarray, + positions: np.ndarray, +) -> int: + """Return the best position, breaking score ties by node name.""" + + order = np.lexsort((names[positions], -scores[positions])) + return int(positions[order[0]]) + + +def _relative_gap( + distances: np.ndarray, + available: np.ndarray, +) -> np.ndarray: + """Normalize available distances to their current finite maximum.""" + + maximum = float(np.max(distances[available], initial=0.0)) + if not np.isfinite(maximum) or maximum <= np.finfo(np.float64).eps: + return np.zeros(len(distances), dtype=np.float64) + return distances / maximum + + +def _eligible_gap_positions( + relative_gap: np.ndarray, + available: np.ndarray, +) -> np.ndarray: + """Return positions within the configured fraction of the largest gap.""" + + maximum = float(np.max(relative_gap[available], initial=0.0)) + if maximum <= np.finfo(np.float64).eps: + return np.flatnonzero(available) + threshold = _COVERAGE_ELIGIBILITY_RATIO * maximum + return np.flatnonzero( + available & (relative_gap >= threshold - np.finfo(np.float64).eps) + ) + + +def _allocate_component_quotas( + coverage_graph: csr_matrix, + candidate_full_indices: np.ndarray, + candidate_scores: np.ndarray, + *, + sensor_num: int, +) -> tuple[np.ndarray, np.ndarray]: + """Allocate sensor counts by active pipe length with candidate caps.""" + + component_count, node_components = connected_components( + coverage_graph, + directed=False, + return_labels=True, + ) + candidate_components = node_components[candidate_full_indices] + capacities = np.bincount( + candidate_components, + minlength=component_count, + ).astype(np.int64, copy=False) + + # The graph is symmetric. Summed row weights count every physical edge + # twice, hence the division by two after aggregation by component. + node_lengths = np.asarray(coverage_graph.sum(axis=1)).ravel() + component_lengths = np.bincount( + node_components, + weights=node_lengths, + minlength=component_count, + ) / 2.0 + component_best_scores = np.full(component_count, -np.inf, dtype=np.float64) + np.maximum.at( + component_best_scores, + candidate_components, + candidate_scores, + ) + + active_components = np.flatnonzero(capacities) + quotas = np.zeros(component_count, dtype=np.int64) + if len(active_components) > sensor_num: + order = np.lexsort( + ( + active_components, + -component_best_scores[active_components], + -component_lengths[active_components], ) ) - # 泵、起终点为tank、reservoir的管道 - # self.delpipes: list[str],需要删除的管道（包括水库、泵、阀门连接的管道） - self.delpipes = list( - set(wn.pump_name_list).union(wn.valve_name_list).union(Reservoirpipe) - ) - self.pipes = [pipe for pipe in wn.pipe_name_list if pipe not in self.delpipes] - # self.L: list[float],所有管道的长度（以米为单位） - self.L = wn.query_link_attribute("length")[self.pipes].tolist() - self.n = len(self.nodes) - self.m = len(self.pipes) - # self.unit_headloss: list[float],单位水头损失（headloss 数据的第一行，单位：米/km） - self.unit_headloss = self.results.link["headloss"].iloc[0, :].tolist() - ## - self.delnodes1 = list(set(ReservoirIndex).union(Tankindex)) + quotas[active_components[order[:sensor_num]]] = 1 + return candidate_components, quotas - # === 改动新增部分：筛选管径小于 min_diameter 的管道节点 === - self.less_than_min_diameter_junction_list = [] - for pipe in self.pipes: - diameter = wn.links[pipe].diameter - if diameter < min_diameter: - start_node = wn.links[pipe].start_node.name - end_node = wn.links[pipe].end_node.name - self.less_than_min_diameter_junction_list.extend([start_node, end_node]) - # 去重 - self.less_than_min_diameter_junction_list = list( - set(self.less_than_min_diameter_junction_list) - ) - - # Step3.2: 计算水力距离 - def CtoS(self): - """ - 计算水力距离矩阵 - :return: - """ - # 水力距离：当行索引对应的节点为控制点时，列索引对应的节点距离控制点的（路径*水头损失）的最小值 - # nodes：list[str]（节点名称） - nodes = copy.deepcopy(self.nodes) - # pipes：list[str]（管道名称） - pipes = self.pipes - wn = self.wn - # n / m：int（节点数 / 管道数） - n = self.n - m = self.m - s1 = [0] * m - q = self.q - L = self.L - # H1：pandas.DataFrame,水头数据，索引为时间步长，列为节点名 - H1 = self.results.node["head"].T - # hh：list[float],计算管道两端水头之差 - hh = [] - # 水头损失 - for p in pipes: - h1 = self.wn.links[p].start_node.name - h1 = H1.loc[str(h1)] - h2 = self.wn.links[p].end_node.name - h2 = H1.loc[str(h2)] - hh.append(abs(h1 - h2)) - hh = np.array(hh) - # headloss：pandas.DataFrame,管道水头损失矩阵 - headloss = pd.DataFrame(hh, index=pipes).T - # s1:管道阻力系数，s2：将管道阻力系数与管道的起始节点和终止节点对应 - hf = pd.DataFrame( - np.array([0] * (n**2)).reshape(n, n), - index=nodes, - columns=nodes, - dtype=float, - ) - weightL = pd.DataFrame( - np.array([0] * (n**2)).reshape(n, n), - index=nodes, - columns=nodes, - dtype=float, - ) - # s2为对应管道起始节点与终止节点的粗糙度系数矩阵，index代表起始节点，columns代表终止节点 - G = nx.DiGraph() - for i in range(0, m): - pipe = pipes[i] - a = wn.links[pipe].start_node.name - b = wn.links[pipe].end_node.name - if q.loc[0, pipe] > 0: - hf.loc[a, b] = headloss.loc[0, pipe] - weightL.loc[a, b] = headloss.loc[0, pipe] * L[i] - G.add_weighted_edges_from([(a, b, weightL.loc[a, b])]) - - else: - hf.loc[b, a] = headloss.loc[0, pipe] - weightL.loc[b, a] = headloss.loc[0, pipe] * L[i] - G.add_weighted_edges_from([(b, a, weightL.loc[b, a])]) - - hydraulicL = pd.DataFrame( - np.array([0] * (n**2)).reshape(n, n), - index=nodes, - columns=nodes, - dtype=float, - ) - - for a in nodes: - if a in G.nodes: - d = nx.shortest_path_length(G, source=a, weight="weight") - for b in list(d.keys()): - hydraulicL.loc[a, b] = d[b] - - hydraulicL = hydraulicL.drop(self.delnodes) - hydraulicL = hydraulicL.drop(self.delnodes, axis=1) - - # 求加权水力距离 - return hydraulicL, G - - # Step3.3: 计算灵敏度矩阵 - # 获取关系矩阵 - def get_Conn(self): - """ - 计算管网连接关系矩阵 - :return: - """ - m = self.wn.num_links - n = self.wn.num_nodes - p = self.wn.num_pumps - v = self.wn.num_valves - - self.nonjunc_index = [] - self.non_link_index = [] - for r in self.wn.reservoirs(): - self.nonjunc_index.append(r[0]) - for t in self.wn.tanks(): - self.nonjunc_index.append(t[0]) - # Conn：numpy.matrix，节点-管道连接矩阵，起点 -1，终点 1 - Conn = np.mat( - np.zeros([n, m - p - v]) - ) # 节点和管道的关系矩阵，行为节点，列为管道，起点为-1，终点为1 - # NConn：numpy.matrix，节点-节点连接矩阵，有管道相连的地方设为 1 - NConn = np.mat(np.zeros([n, n])) # 节点之间的关系，之间有管道为1，反之为0 - # pipes：list[str]，去除泵和阀门的管道列表 - pipes = [ - pipe - for pipe in self.wn.pipes() - if pipe not in self.wn.pumps() and pipe not in self.wn.valves() + quotas[active_components] = 1 + remaining = sensor_num - len(active_components) + while remaining > 0: + available = active_components[ + quotas[active_components] < capacities[active_components] ] - for pipe_name, pipe in pipes: - start = self.wn.node_name_list.index(pipe.start_node_name) - end = self.wn.node_name_list.index(pipe.end_node_name) - p_index = self.wn.link_name_list.index(pipe_name) - Conn[start, p_index] = -1 - Conn[end, p_index] = 1 - NConn[start, end] = 1 - NConn[end, start] = 1 - self.A = Conn - link_name_list = [ - link - for link in self.wn.link_name_list - if link not in self.wn.pump_name_list - and link not in self.wn.valve_name_list - ] - self.A2 = pd.DataFrame( - self.A, index=self.wn.node_name_list, columns=link_name_list + if len(available) == 0: + raise ValueError("连通区域中的候选节点不足，无法分配监测点名额") + + weights = component_lengths[available] + if float(weights.sum()) <= 0: + weights = (capacities[available] - quotas[available]).astype( + np.float64, + copy=False, + ) + ideal = remaining * weights / float(weights.sum()) + whole = np.minimum( + np.floor(ideal).astype(np.int64), + capacities[available] - quotas[available], ) - self.A2 = self.A2.drop(self.delnodes) - for pipe in self.delpipes: - if ( - pipe not in self.wn.pump_name_list - and pipe not in self.wn.valve_name_list - ): - self.A2 = self.A2.drop(columns=pipe) - self.junc_list = self.A2.index - self.A2 = np.mat(self.A2) # 节点管道关系 - self.A3 = NConn + whole_count = int(whole.sum()) + if whole_count: + quotas[available] += whole + remaining -= whole_count + continue - def Jaco(self, hL: pandas.DataFrame): - """ - 计算灵敏度矩阵（节点压力对粗糙度变化的响应） - :param hL: 水力距离矩阵 - :return: - """ - # global result - # A：numpy.matrix, 节点-管道关系矩阵 - A = self.A2 - wn = self.wn + fractional = ideal - np.floor(ideal) + order = np.lexsort( + ( + available, + -component_best_scores[available], + -weights, + -fractional, + ) + ) + for component in available[order]: + quotas[component] += 1 + remaining -= 1 + if remaining == 0: + break - try: - result = wntr.sim.EpanetSimulator(wn).run_sim() - except EpanetException: - pass - finally: - h = result.link["headloss"][self.pipes].values[0] - q = result.link["flowrate"][self.pipes].values[0] - l = self.wn.query_link_attribute("length")[self.pipes] - C = self.wn.query_link_attribute("roughness")[self.pipes] - # headloss：numpy.ndarray,水头损失数组 - headloss = np.array(h) - # 调整流量方向 - for i in range(0, len(q)): - if q[i] < 0: - A[:, i] = -A[:, i] - # q：numpy.ndarray,流量数组 - q = np.abs(q) - # 两个灵敏度矩阵 - # B / S：numpy.matrix,灵敏度计算的中间矩阵 - B = np.mat(np.diag(q / ((1.852 * headloss) + 1e-10))) - S = np.mat(np.diag(q / C)) - # X：numpy.matrix, 灵敏度矩阵 - X = A * B * A.T - try: - det = np.linalg.det(X) - except RuntimeError as e: - sign, logdet = slogdet(X) # 防止溢出 - det = sign * np.exp(logdet) - if det != 0: - J_H_Cw = X.I * A * S - # J_H_Q = -X.I - J_q_Cw = S - B * A.T * X.I * A * S # 去掉了delnodes和delpipes - # J_q_Q = B * A.T * X.I - else: # 当X不可逆 - J_H_Cw = np.linalg.pinv(X) @ A @ S - # J_H_Q = -np.linalg.pinv(X) - J_q_Cw = S - B * A.T * np.linalg.pinv(X) * A * S - # J_q_Q = B * A.T * np.linalg.pinv(X) - - Sen_pressure = [] - S_pressure = np.abs(J_H_Cw).sum(axis=1).tolist() # 修改为绝对值 - for ss in S_pressure: - Sen_pressure.append(ss[0]) - # 求总灵敏度 - SS_pressure = copy.deepcopy(hL) - for i in range(0, len(Sen_pressure)): - SS_pressure.iloc[i, :] = SS_pressure.iloc[i, :] * Sen_pressure[i] - SS = copy.deepcopy(hL) - for i in range(0, len(Sen_pressure)): - SS.iloc[i, :] = SS.iloc[i, :] * Sen_pressure[i] - # SS[i,j]:节点nodes[i]的灵敏度*该节点到nodes[j]的水力距离 - return SS + return candidate_components, quotas -# 2025/03/12 -# Step4: 传感器布置优化 -# Sensorplacement -# weight：分配权重 -# sensor：传感器布置的位置 -class Sensorplacement(wn_func): - """ - Sensorplacement 类继承了 wn_func 类，并且用于计算和优化传感器布置的位置。 +def _select_component_positions( + coverage_graph: csr_matrix, + candidates: _CandidatePool, + component_positions: np.ndarray, + existing_positions: list[int], + *, + quota: int, +) -> list[int]: + """Select one component's sensors with score-aware farthest-first search.""" + + local_coordinates = candidates.coordinates[component_positions] + local_names = candidates.names[component_positions] + local_scores = candidates.scores[component_positions] + nearest_geographic = np.full(len(component_positions), np.inf) + nearest_topological = np.full(len(component_positions), np.inf) + + for position in existing_positions: + nearest_geographic = np.minimum( + nearest_geographic, + np.linalg.norm( + local_coordinates - candidates.coordinates[position], + axis=1, + ), + ) + + if existing_positions: + all_local = np.ones(len(component_positions), dtype=bool) + seed_eligible = _eligible_gap_positions( + _relative_gap(nearest_geographic, all_local), + all_local, + ) + else: + seed_eligible = np.arange(len(component_positions), dtype=np.int64) + + seed = _highest_scoring_position( + local_names, + local_scores, + seed_eligible, + ) + selected_local = [seed] + remaining = np.ones(len(component_positions), dtype=bool) + remaining[seed] = False + + while len(selected_local) < quota: + newest = selected_local[-1] + geographic_distance = np.linalg.norm( + local_coordinates - local_coordinates[newest], + axis=1, + ) + nearest_geographic = np.minimum( + nearest_geographic, + geographic_distance, + ) + + source = int(candidates.full_indices[component_positions[newest]]) + topological_distance = dijkstra( + coverage_graph, + directed=False, + indices=source, + return_predecessors=False, + )[candidates.full_indices[component_positions]] + nearest_topological = np.minimum( + nearest_topological, + topological_distance, + ) + + coverage_gap = np.maximum( + _relative_gap(nearest_geographic, remaining), + _relative_gap(nearest_topological, remaining), + ) + eligible_local = _eligible_gap_positions(coverage_gap, remaining) + next_local = _highest_scoring_position( + local_names, + local_scores, + eligible_local, + ) + selected_local.append(next_local) + remaining[next_local] = False + + return [int(component_positions[position]) for position in selected_local] + + +def _geographic_coverage_metrics( + candidate_coordinates: np.ndarray, + selected_positions: list[int], +) -> tuple[float, float, float]: + normalized = _isotropic_coordinates(candidate_coordinates) + nearest = np.full(len(normalized), np.inf) + for position in selected_positions: + nearest = np.minimum( + nearest, + np.linalg.norm(normalized - normalized[position], axis=1), + ) + + selected_coordinates = normalized[selected_positions] + if len(selected_positions) < 2: + minimum_gap = 0.0 + else: + pairwise = np.linalg.norm( + selected_coordinates[:, None, :] - selected_coordinates[None, :, :], + axis=2, + ) + np.fill_diagonal(pairwise, np.inf) + minimum_gap = float(pairwise.min()) + return ( + float(nearest.max()), + float(np.quantile(nearest, 0.95)), + minimum_gap, + ) + + +def _select_sensor_nodes( + prepared: _PreparedNetwork, + log_sensitivity: np.ndarray, + hydraulic_distance_sums: np.ndarray, + *, + sensor_num: int, +) -> list[str]: + candidate_indices = prepared.candidate_indices + if len(candidate_indices) < sensor_num: + raise ValueError( + "满足最小管径要求的候选节点少于请求的监测点数量：" + f"候选 {len(candidate_indices)} 个，请求 {sensor_num} 个" + ) + + candidates = _build_candidate_pool( + prepared, + log_sensitivity, + hydraulic_distance_sums, + ) + candidate_components, component_quotas = _allocate_component_quotas( + prepared.coverage_graph, + candidates.full_indices, + candidates.scores, + sensor_num=sensor_num, + ) + selected_positions: list[int] = [] + quota_components = np.flatnonzero(component_quotas) + component_order = np.lexsort( + (quota_components, -component_quotas[quota_components]) + ) + for component in quota_components[component_order]: + component_positions = np.flatnonzero(candidate_components == component) + selected_positions.extend( + _select_component_positions( + prepared.coverage_graph, + candidates, + component_positions, + selected_positions, + quota=int(component_quotas[component]), + ) + ) + + selected_array = np.asarray(selected_positions, dtype=np.int64) + selected_order = np.lexsort( + ( + candidates.names[selected_array], + -candidates.scores[selected_array], + ) + ) + selected_positions = selected_array[selected_order].tolist() + maximum_radius, p95_radius, minimum_gap = _geographic_coverage_metrics( + prepared.coordinates[candidate_indices], + selected_positions, + ) + logger.info( + "Sensitivity placement coverage: components=%d max_radius=%.6f " + "p95_radius=%.6f min_sensor_gap=%.6f", + int(np.count_nonzero(component_quotas)), + maximum_radius, + p95_radius, + minimum_gap, + ) + return [str(candidates.names[position]) for position in selected_positions] + + +def optimize_sensor_placement( + wn: wntr.network.WaterNetworkModel, + sensor_num: int, + min_diameter: int, +) -> list[str]: + """Return deterministic pressure monitoring nodes for a loaded network. + + ``min_diameter`` is expressed in millimetres, matching the HTTP contract. + A node is a valid installation candidate when at least one incident pipe + meets the threshold. All valid hydraulic nodes still participate in the + sensitivity calculation so small pipes continue to influence the result. """ - def __init__( - self, wn: wntr.network.WaterNetworkModel, sensornum: int, min_diameter: int - ): - """ + if sensor_num <= 0: + raise ValueError("监测点数量必须大于 0") + if min_diameter < 0: + raise ValueError("最小管径不能小于 0") - :param wn: 由wntr生成的模型 - :param sensornum: 传感器的数量 - :param min_diameter: 安装的最小管径 - """ - wn_func.__init__(self, wn, min_diameter=min_diameter) - self.sensornum = sensornum + total_started = perf_counter() + simulation_started = total_started + results = _run_hydraulic_simulation(wn) + simulation_seconds = perf_counter() - simulation_started - # 1.某个节点到所有节点的加权距离之和 - # 2.某个节点到该组内所有节点的加权距离之和 - def sensor( - self, SS: pandas.DataFrame, G: networkx.Graph, group: dict[int, list[str]] - ): - """ - sensor 方法是用来根据灵敏度矩阵 SS 和加权图 G 来确定传感器布置位置的 - :param SS: 灵敏度矩阵，每个节点的行和列代表不同节点，矩阵元素表示节点间的灵敏度。SS.iloc[i, :] 表示第 i 行对应节点 i 到所有其他节点的灵敏度 - :param G: 加权图，表示管网的拓扑结构，每个节点通过管道连接。图的边的权重通常是根据水力距离或者流量等计算的 - :param group: 节点分组，字典的键是分组编号，值是该组的节点名称列表 - :return: - """ - # 传感器布置个数以及位置 - # W = self.weight() - n = self.n - len(self.delnodes) - nodes = copy.deepcopy(self.nodes) - for node in self.delnodes: - nodes.remove(node) - # sumSS：list[float],每个节点到其他节点的灵敏度之和。SS.iloc[i, :] 返回第 i 个节点与所有其他节点的灵敏度值，sum(SS.iloc[i, :]) 计算这些灵敏度值的总和。 - sumSS = [] - for i in range(0, n): - sumSS.append(sum(SS.iloc[i, :])) - # 一个整数范围，表示每个节点的索引，用作sumSS_ DataFrame的索引 - indices = range(0, n) - # sumSS_：pandas.DataFrame,将 sumSS 转换成 DataFrame 格式，并且将节点的总灵敏度保存到 CSV 文件 sumSS_data.csv 中 - sumSS_ = pd.DataFrame(np.array(sumSS), index=indices) - # sumSS_.to_csv('sumSS_data.csv') # 存储节点总灵敏度 + preparation_started = perf_counter() + prepared = _prepare_network(wn, results, min_diameter=min_diameter) + preparation_seconds = perf_counter() - preparation_started - # sumSS：pandas.DataFrame,sumSS 被转换为 DataFrame 类型，并且按总灵敏度（即灵敏度之和）降序排列。此时，sumSS 是按节点的灵敏度之和排序的 DataFrame - sumSS = pd.DataFrame(np.array(sumSS), index=nodes) - sumSS = sumSS.sort_values(by=[0], ascending=[False]) - # sensorindex：list[str],用于存储根据灵敏度排序选出的传感器位置的节点名称,存储根据总灵敏度排序的节点列表，用于传感器布置 - sensorindex = [] - # sensorindex_2：list[str],用于存储每组内根据灵敏度排序选出的传感器位置的节点名称,存储每个组内根据灵敏度排序选择的传感器节点 - sensorindex_2 = [] - # group_S：dict[int, pandas.DataFrame],存储每个组内的灵敏度矩阵 - group_S = {} - # group_sumSS：dict[int, list[float]],存储每个组内节点的总灵敏度，值为每个组内节点灵敏度之和的列表 - group_sumSS = {} + sensitivity_started = perf_counter() + log_sensitivity = _estimate_log_pressure_sensitivity(prepared) + sensitivity_seconds = perf_counter() - sensitivity_started - # 改动 - for i in range(0, len(group)): - for node in self.delnodes: - # 这里的group[i]是每个组的节点列表，代码首先去除已经被标记为删除的节点self.delnodes - if node in group[i]: - group[i].remove(node) - group_S[i] = SS.loc[group[i], group[i]] - # 对每个组内的节点，计算组内节点的总灵敏度（group_sumSS[i]）。它将每个组内节点的灵敏度值相加，并且按灵敏度降序排序 - group_sumSS[i] = [] - for j in range(0, len(group[i])): - group_sumSS[i].append(sum(group_S[i].iloc[j, :])) - group_sumSS[i] = pd.DataFrame(np.array(group_sumSS[i]), index=group[i]) - group_sumSS[i] = group_sumSS[i].sort_values(by=[0], ascending=[False]) - for node in self.less_than_min_diameter_junction_list: - # 这里的group_sumSS[i]是每个分组的灵敏度节点排序列表，去除已经被标记为删除的节点self.less_than_min_diameter_junction_list - if node in group_sumSS[i]: - group_sumSS[i].remove(node) - pass + distance_started = perf_counter() + hydraulic_distance_sums = _estimate_hydraulic_distance_sums(prepared) + distance_seconds = perf_counter() - distance_started - # 1.选sumSS最大的节点，然后把这个节点所在的那个组删掉，就可以不再从这个组选点。再重新排序选sumSS最大的; - # 2.在每组内选group_sumSS最大的节点 - # 在这个循环中，首先选择灵敏度最高的节点Smaxnode并添加到sensorindex。然后根据灵敏度排序，删除已选的节点并继续选择下一个灵敏度最大的节点。这个过程用于选择传感器的位置 - sensornum = self.sensornum - for i in range(0, sensornum): - # Smaxnode：str,最大灵敏度节点，sumSS.index[0] 表示灵敏度最高的节点 - Smaxnode = sumSS.index[0] - sensorindex.append(Smaxnode) - sensorindex_2.append(group_sumSS[i].index[0]) + selection_started = perf_counter() + selected = _select_sensor_nodes( + prepared, + log_sensitivity, + hydraulic_distance_sums, + sensor_num=sensor_num, + ) + selection_seconds = perf_counter() - selection_started - for key, value in group.items(): - if Smaxnode in value: - sumSS = sumSS.drop(index=group[key]) - continue + logger.info( + "Sensitivity placement completed: nodes=%d pipes=%d candidates=%d " + "sensors=%d seconds=%.3f " + "(simulation=%.3f preparation=%.3f sensitivity=%.3f " + "distance=%.3f selection=%.3f)", + len(prepared.node_names), + prepared.incidence.shape[1], + len(prepared.candidate_indices), + len(selected), + perf_counter() - total_started, + simulation_seconds, + preparation_seconds, + sensitivity_seconds, + distance_seconds, + selection_seconds, + ) + return selected - sumSS = sumSS.sort_values(by=[0], ascending=[False]) - return sensorindex, sensorindex_2 +def optimize_sensor_placement_from_inp( + inp_path: str | Path, + sensor_num: int, + min_diameter: int, +) -> list[str]: + """Load an EPANET INP model and run the unified placement algorithm.""" + + wn = wntr.network.WaterNetworkModel(str(inp_path)) + return optimize_sensor_placement( + wn, + sensor_num=sensor_num, + min_diameter=min_diameter, + ) -# 2025/03/13 def get_ID(name: str, sensor_num: int, min_diameter: int) -> list[str]: - """ - 获取布置测压点的坐标,初始测压点布置根据灵敏度来布置，计算初始情况下的校准过程的error - :param name: 数据库名称 - :param sensor_num: 测压点数目 - :param min_diameter: 安装的最小管径 - :return: 测压点节点ID - """ - # inp_file_real：str,输入文件名，表示原始水力模型文件的路径,该文件格式为 EPANET 输入文件（.inp），包含管网的结构信息、节点、管道、泵等数据 - inp_file_real = f"./db_inp/{name}.db.inp" - # sensornum：int,需要布置的传感器数量 - # sensornum = sensor_num - # wn_real：wntr.network.WaterNetworkModel,加载 EPANET 水力模型 - wn_real = wntr.network.WaterNetworkModel(inp_file_real) # 真实粗糙度的原始管网 - # sim_real：wntr.sim.EpanetSimulator,创建一个水力仿真器对象 - sim_real = wntr.sim.EpanetSimulator(wn_real) - # results_real：wntr.sim.results.SimulationResults,运行仿真并返回结果 - results_real = sim_real.run_sim() + """Compatibility entry point used by the sensor placement service.""" - # real_C：list[float]，包含所有管道粗糙度的列表 - real_C = wn_real.query_link_attribute("roughness").tolist() - # wn_fun1：wn_func（继承自 object），创建 wn_func 类的实例，传入 wn_real 水力模型对象。wn_func 用于计算管网相关的水力属性，比如水力距离、灵敏度等 - wn_fun1 = wn_func(wn_real, min_diameter=min_diameter) - # nodes：list[str]，管网的节点名称列表 - nodes = wn_fun1.nodes - # delnodes：list[str]，被删除的节点（如水库、泵、阀门连接的节点等） - delnodes = wn_fun1.delnodes - # Coor_node：pandas.DataFrame - Coor_node = getCoor(wn_real) - Coor_node = Coor_node.drop(wn_fun1.delnodes) - nodes = [node for node in wn_fun1.nodes if node not in delnodes] - # coordinates：pandas.Series，存储所有节点的坐标，类型为 Series，索引为节点名称，值为 (x, y) 坐标对 - coordinates = wn_fun1.coordinates - - # 随机产生监测点 - # junctionnum：int，nodes 的长度，表示节点的数量 - junctionnum = len(nodes) - # random_numbers：list[int]，使用 random.sample 随机选择 sensornum（20）个节点的编号。它返回一个不重复的随机编号列表 - # random_numbers = random.sample(range(junctionnum), sensor_num) - # for i in range(sensor_num): - # # print(random_numbers[i]) - - wn_fun1.get_Conn() - # hL：pandas.DataFrame，水力距离矩阵，表示每个节点到其他节点的水力阻力 - # G：networkx.DiGraph，加权有向图，表示管网的拓扑结构，节点之间的边带有权重 - hL, G = wn_fun1.CtoS() - # SS：pandas.DataFrame，灵敏度矩阵，表示每个节点对管网变化（如粗糙度、流量等）的响应 - SS = wn_fun1.Jaco(hL) - # group：dict[int, list[str]]，使用 kgroup 函数将节点按坐标分成若干组，每组包含的节点数不一定相同。group 是一个字典，键为分组编号，值为节点名列表 - - G1 = wn_real.to_graph() - G1 = G1.to_undirected() # 变为无向图 - - group = kgroup(Coor_node, sensor_num) - # group = skater_partition(G1, sensor_num) - # group = spectral_partition(G1, sensor_num) - - # print(group) - # --------------------- 保存 group 数据 --------------------- - # 将 group 数据转换为一个“长格式”的 DataFrame， - # 每一行记录一个节点及其所属的分组 - # group_data = [] - # for group_id, node_list in group.items(): - # for node in node_list: - # group_data.append({"Group": group_id, "Node": node}) - # - # df_group = pd.DataFrame(group_data) - # - # # 保存为 Excel 文件，文件名为 "group.xlsx"；index=False 表示不保存行索引 - # df_group.to_excel("group.xlsx", index=False) - - # wn_fun：Sensorplacement（继承自wn_func） - # 创建Sensorplacement类的实例，传入水力网络模型wn_real和传感器数量sensornum。Sensorplacement用于计算和布置传感器 - wn_fun = Sensorplacement(wn_real, sensor_num, min_diameter=min_diameter) - wn_fun.__dict__.update(wn_fun1.__dict__) - # sensorindex：list[str]，初始传感器布置位置的节点名称 - # sensorindex_2：list[str]，根据分组选择的传感器位置 - sensorindex, sensorindex_2 = wn_fun.sensor(SS, G, group) # 初始的sensorindex - # print(str(sensor_num), "个测压点，测压点位置：", sensorindex) - - # 重新打开数据库 - # if is_project_open(name=name): - # close_project(name=name) - # open_project(name=name) - # for node_id in sensorindex : - # sensor_coord[node_id] = get_node_coord(name=name, node_id=node_id) - # close_project(name=name) - # print(sensor_coord) - # # 分区画图 - # colorlist = ['lightpink', 'coral', 'rosybrown', 'olive', 'powderblue', 'lightskyblue', 'steelblue', 'peachpuff','brown','silver','indigo','lime','gold','violet','maroon','navy','teal','magenta','cyan', - # 'burlywood', 'tan', 'slategrey', 'thistle', 'lightseagreen', 'lightgreen', 'red','blue','yellow','orange','purple','grey','green','pink','lightblue','beige','chartreuse','turquoise','lavender','fuchsia','coral'] - # G = wn_real.to_graph() - # G = G.to_undirected() # 变为无向图 - # pos = nx.get_node_attributes(G, 'pos') - # pass - # - # for i in range(0, sensor_num): - # ax = plt.gca() - # ax.set_title(inp_file_real + str(sensor_num)) - # nodes = nx.draw_networkx_nodes(G, pos, nodelist=group[i], node_color=colorlist[i], node_size=10) - # nodes = nx.draw_networkx_nodes(G, pos, - # nodelist=sensorindex_2, node_color='red', node_size=70, node_shape='*' - # ) - # edges = nx.draw_networkx_edges(G, pos) - # ax.spines['top'].set_visible(False) - # ax.spines['right'].set_visible(False) - # ax.spines['bottom'].set_visible(False) - # ax.spines['left'].set_visible(False) - # plt.savefig(inp_file_real + str(sensor_num) + ".png", dpi=300) - # plt.show() - # - # wntr.graphics.plot_network(wn_real, node_attribute=sensorindex_2, node_size=50, node_labels=False, - # title=inp_file_real + '_Projetion' + str(sensor_num)) - # plt.savefig(inp_file_real + '_S' + str(sensor_num) + ".png", dpi=300) - # plt.show() - return sensorindex - - -if __name__ == "__main__": - sensorindex = get_ID(name=project_info.name, sensor_num=20, min_diameter=300) - - print(sensorindex) - # 将 sensor_coord 字典转换为 DataFrame， - # 使用 orient='index' 表示字典的键作为 DataFrame 的行索引， - # 数据中每个键对应的 value 是一个子字典，其键 'x' 和 'y' 成为 DataFrame 的列名 - # df_sensor_coord = pd.DataFrame.from_dict(sensor_coord, orient='index') - # - # # 将索引名称设为 'Node' - # df_sensor_coord.index.name = 'Node' - # - # # 保存到 Excel 文件 - # df_sensor_coord.to_excel("sensor_coord.xlsx", index=True) + inp_path = Path("db_inp") / f"{name}.db.inp" + return optimize_sensor_placement_from_inp( + inp_path, + sensor_num=sensor_num, + min_diameter=min_diameter, + ) diff --git a/docs/sensor-sensitivity-optimization.md b/docs/sensor-sensitivity-optimization.md new file mode 100644 index 0000000..85aecc6 --- /dev/null +++ b/docs/sensor-sensitivity-optimization.md @@ -0,0 +1,270 @@ +# 监测点灵敏度算法优化说明 + +本文记录压力监测点布置算法的改造方案、数学含义、复杂度变化和本地验证结果。实现位于 `app/algorithms/sensor/sensitivity.py`，对外兼容入口仍为 `get_ID(name, sensor_num, min_diameter)`。 + +## 改造目标 + +原算法根据节点压力对管道粗糙度变化的灵敏度，以及节点到全网的水力距离，给每个候选节点评分。空间聚类用于限制监测点过度集中。业务评分可以写成： + +$$ +score_i = sensitivity_i \times \sum_j distance(i, j) +$$ + +本次改造保留这套评分语义和空间覆盖原则，主要处理以下问题： + +- 稠密的节点与节点、节点与管道、管道与管道矩阵占用大量内存。 +- 显式计算行列式、逆矩阵和伪逆，计算量随节点数快速增长。 +- 从每个节点执行最短路径并保存全节点距离矩阵，时间和内存均为平方级。 +- 一次任务重复运行 EPANET，并保留算法没有使用的扩展时序结果。 +- 聚类没有固定随机状态，相同输入可能返回不同结果。 +- KMeans 按候选节点密度最小化平方距离，密集管网会系统性获得更多监测点，不能保证地图或管网路径覆盖均匀。 +- WNTR 管径使用米，接口参数使用毫米，旧候选过滤没有统一单位。 + +## 稀疏表示替代稠密矩阵 + +管网的节点通常只连接少量管道。新实现使用 SciPy CSR 或 CSC 矩阵保存节点与管道关联矩阵、水力有向图和水力雅可比矩阵，只记录真实存在的连接。 + +准备阶段生成一个只读的 `_PreparedNetwork`，其中包括： + +- 参与分析的节点及其在完整模型中的索引。 +- 符合安装条件的候选节点索引。 +- 归一化前的二维坐标。 +- 稀疏节点与管道关联矩阵。 +- 管道导通系数和粗糙度响应系数。 +- 按实际流向构建的稀疏水力距离图。 +- 按初始运行状态构建的无向物理覆盖图。 + +水库、水箱、泵和阀门的边界节点继续沿用原算法的排除规则。与水源直接相连的管道不参加扰动计算。平行有向边只保留权重最小的边，以符合最短路径语义。 + +主体存储由平方级降到接近 `O(n + m)`。稀疏 LU 分解可能产生填充，其内存仍取决于管网拓扑，但实现不会再主动分配完整的 `n × n` 或 `n × m` 数组。 + +## 用 Cauchy 投影估计压力灵敏度 + +旧算法先显式计算压力对管道粗糙度的响应矩阵： + +$$ +J = X^{-1} A S +$$ + +其中，`X` 是节点水力雅可比矩阵，`A` 是节点与管道关联矩阵，`S` 是管道粗糙度响应矩阵。节点灵敏度是 `J` 对应行的 L1 范数。完整的 `J` 为 `节点数 × 管道数`，大模型无法保存。 + +新算法使用 Cauchy 分布的 1-stable 特性估计每一行的 L1 范数。设 `R` 为 Cauchy 随机投影矩阵，只需求解： + +$$ +Y = X^{-1} A S R +$$ + +对节点 `i` 而言，`Y` 中每个投影值服从以 `||J_i||_1` 为尺度的 Cauchy 分布。实现使用投影绝对值的对数几何均值估计 `log(||J_i||_1)`，不保存完整灵敏度矩阵。 + +当前固定参数如下： + +| 参数 | 数值 | 用途 | +|---|---:|---| +| 随机种子 | 42 | 保证结果可复现 | +| Cauchy 投影数 | 256 | 控制灵敏度估计精度 | +| 投影批大小 | 16 | 限制投影中间矩阵内存 | +| 雅可比正则化 | `max(abs(diag(X))) × sqrt(eps)` | 改善接近奇异矩阵的稳定性 | +| 稀疏排序 | `MMD_AT_PLUS_A` | 减少 LU 分解填充 | + +水力雅可比矩阵只进行一次稀疏 LU 分解，256 次投影复用该分解结果。批处理期间最多保留 16 列投影数据。 + +## 用空间代表点估计水力距离总和 + +旧算法从每个节点执行一次最短路径，并保存 `n × n` 的水力距离矩阵。新算法根据节点坐标构建最多 256 个空间代表点： + +1. 使用固定随机状态的 `MiniBatchKMeans` 对节点坐标分组。 +2. 每组选择最靠近聚类中心的节点作为代表点。 +3. 代表点权重等于该组包含的节点数。 +4. 在反向水力图上从代表点执行有向 Dijkstra，得到原图中各节点到代表点的距离。 + +节点 `i` 的全网水力距离总和估计为： + +$$ +distance\_sum_i \approx \sum_{l \in landmarks} weight_l \times distance(i, l) +$$ + +Dijkstra 每批处理 16 个代表点，内存中只保留当前距离块。不可达距离按原算法约定计为 0，避免断开分支得到无穷评分。 + +## 评分约束下的混合覆盖选点 + +最终排序使用对数形式： + +$$ +log\_score_i = log\_sensitivity_i + \log(distance\_sum_i) +$$ + +对数是单调函数，因此该排序等价于比较 `sensitivity_i × distance_sum_i`，同时可以避免大数乘法溢出。 + +最终选点不再对全部候选节点执行 KMeans。KMeans 的目标函数按候选节点数计权，节点密集区域即使地理范围较小，也会获得更多聚类中心。本次改为评分约束下的最远空白区优先策略。 + +候选坐标使用同一个尺度因子归一化，保留管网原始长宽比。初始状态下开启的管道按物理长度构建无向覆盖图，开启的泵和阀门作为点连接；关闭连接会形成独立区域。监测点名额首先按各连通区域的有效管道长度分配：名额足够时每个区域至少一个，其余名额按最大余数法分配，并受该区域候选数量限制。 + +选点按区域名额从多到少处理，使主干管网先形成覆盖骨架。第一个区域的首点取综合评分最高的候选；后续区域的首点也必须考虑此前所有区域的已选点，先进入全局地图空白度前 70% 的候选集，再比较综合评分。这一约束避免两个拓扑断开但地图上重叠或相邻的区域各自选择一个近邻高分点。 + +设全局已选集合为 `S`，其余候选的最近地图距离和本连通区域内的最近管网路径距离分别为： + +$$ +g_i = \min_{s \in S} ||x_i-x_s||, \qquad +t_i = \min_{s \in S} shortest\_path(i, s) +$$ + +两类距离分别除以当前未选候选中的最大值，混合空白度取两者较大值： + +$$ +coverage_i = \max\left(\frac{g_i}{\max g}, \frac{t_i}{\max t}\right) +$$ + +每轮只保留空白度达到当前最大值 70% 的候选，再从中选择综合评分最高的节点。地图距离始终相对全部已选区域更新，管网路径距离在当前连通区域内更新。这样地图或管网路径中任一维度仍有明显空白时，该区域都不会被高密度节点误判为已经覆盖。综合评分相同时按节点 ID 排序，结果数量、唯一性和可复现性保持不变。 + +实现只维护候选节点的最近距离数组。每新增一个监测点，执行一次单源 Dijkstra 并增量更新，不构造候选节点两两距离矩阵。 + +所有模型使用同一套算法和相同参数，不根据节点数量切换实现。 + +## EPANET 只计算初始状态 + +后续计算只读取水力结果的第一个时刻。新实现会临时将 `wn.options.time.duration` 设置为 0，只运行一次 EPANET 初始状态模拟，并在结束或异常后恢复原始时长。 + +EPANET 中间文件写入独立的 `TemporaryDirectory`，任务结束后自动清理。并发请求不再共用工作目录下的 `temp.*` 文件。 + +旧调用链最多重复运行约四次 EPANET，新调用链只运行一次，也不会分配没有使用的完整时序结果。 + +## 最小管径规则 + +`min_diameter` 的接口单位是毫米，WNTR 中的管径单位是米。准备阶段使用以下换算： + +$$ +diameter\_mm = diameter\_m \times 1000 +$$ + +节点连接的管道中，只要至少一根达到最小管径，该节点就可以作为安装候选。小管径管道仍参加全网水力和灵敏度计算，管径条件只限制监测点安装位置。 + +当候选节点少于请求的监测点数量时，算法会返回包含候选数量和请求数量的明确错误。 + +## 复杂度变化 + +下表中的 `n` 为参与分析的节点数，`m` 为参与分析的管道数，`k=256` 为灵敏度投影数，`l≤256` 为水力距离代表点数，`b=16` 为批大小。 + +| 环节 | 原实现 | 新实现 | +|---|---|---| +| 节点与管道关系 | 稠密 `n × m` | CSR 稀疏矩阵，约 `O(n + m)` | +| 节点关系与距离 | 多个稠密 `n × n` 矩阵 | 稀疏有向图和流式距离块 | +| 灵敏度 | 稠密行列式、逆矩阵或伪逆，时间接近 `O(n³)` | 一次稀疏 LU 分解和 `k` 次稀疏求解 | +| 灵敏度结果 | 保存完整 `n × m` 响应矩阵 | 保存 `n` 个对数灵敏度和当前投影批 | +| 水力距离 | 从全部节点执行最短路径并保存 `n × n` 结果 | 从至多 `l` 个代表点执行 Dijkstra，每批 `b` 个 | +| 最终选点 | 按节点密度分配的完整 KMeans | 连通区域配额和 70% 混合覆盖，线性距离数组 | +| 水力模拟 | 调用链中重复执行，并可能保存完整时序 | 一次初始状态模拟 | + +稀疏 LU 的时间和内存不能简单视为线性，其填充程度受网络拓扑影响。当前压力测试覆盖到 20.3 万原始节点，不能据此保证任意更大或连接更稠密的模型都保持相同比例。 + +## 精度取舍 + +新实现不逐元素生成旧算法的完整精确矩阵，而是估计最终排序需要的两个统计量： + +- Cauchy 投影估计压力响应矩阵每一行的 L1 范数。 +- 加权空间代表点估计节点到全网的水力距离总和。 + +单元测试使用小型模型构造完整稠密参考结果，要求近似方案选点的精确参考目标值不低于精确方案的 95%。本地模型对比结果如下： + +| 模型 | 近似选点的精确参考目标比 | +|---|---:| +| `fengxian.inp` | 98.80% | +| MD 模型 | 99.71% | + +固定随机种子和固定采样数保证同一模型、监测点数量和最小管径得到相同结果。近似排序仍可能与完整稠密算法不同，特别是多个候选节点得分接近时。 + +混合覆盖会主动放弃部分集中在同一区域的高分节点。单点综合评分之和因此不是唯一质量指标，还需要同时检查未覆盖半径、最小点间距和入选节点的评分百分位。 + +## 资源压力测试记录 + +以下结果来自 2026-08-03 的本地验证。环境为 Linux、Python 3.12、Conda `server` 环境，主机物理内存约 30 GiB。测试参数统一为 20 个监测点、最小管径 0。 + +外部看门狗使用以下停止条件： + +- 算法进程树 RSS 达到 7.5 GiB。 +- 系统可用内存低于 6 GiB。 +- 单模型运行超过 600 秒。 +- 子进程虚拟地址空间硬限制为 8 GiB。 + +任一条件满足时，看门狗会终止整个进程组。三次压力测试均未触发停止条件。 + +| 模型 | 原始节点 | 实际分析节点 | 实际分析管道 | 算法流水线耗时 | 峰值 RSS | +|---|---:|---:|---:|---:|---:| +| `temp/leakage/temp_3698123.inp` | 31,143 | 28,723 | 29,974 | 2.77 秒 | 467 MiB | +| `inp/jbh.inp` | 94,049 | 69,949 | 82,369 | 8.92 秒 | 980 MiB | +| `inp/Todo/v-16常熟模型.inp` | 203,569 | 145,948 | 176,512 | 20.00 秒 | 1.84 GiB | + +实际分析节点少于原始节点，是因为水库、水箱、泵和阀门边界节点按算法规则排除。20.3 万节点模型原文件使用非标准的 `REPORTING TIMESTEP`，并且缺少 `[END]`。压力测试只在隔离临时副本中将其规范化，没有修改原文件。 + +20.3 万节点模型各阶段耗时如下： + +| 阶段 | 耗时 | +|---|---:| +| 模型加载 | 6.63 秒 | +| EPANET 初始状态模拟 | 8.30 秒 | +| 稀疏数据准备 | 3.00 秒 | +| 压力灵敏度估计 | 0.76 秒 | +| 水力距离估计 | 0.94 秒 | +| 监测点选择 | 0.37 秒 | + +### `tjwater` 分布优化前后对比 + +2026-08-03 使用 `db_inp/tjwater.db.inp`、20 个监测点、最小管径 0 进行同机对比。进程通过 systemd scope 限制在 2.5 GiB 内存，并设置 90 秒超时；两次运行均未触发保护。覆盖距离使用保持长宽比后的模型坐标，按管网最大轴跨度归一化。 + +| 指标 | KMeans 选点 | 70% 混合覆盖 | 变化 | +|---|---:|---:|---:| +| 实际分析/候选节点 | 87,877 | 87,877 | 不变 | +| 最大地图覆盖半径 | 0.176786 | 0.173941 | -1.61% | +| P95 地图覆盖半径 | 0.122093 | 0.090686 | -25.72% | +| 最小监测点间距 | 0.000559 | 0.040887 | 约 73.2 倍 | +| 综合评分中位百分位 | 98.76% | 95.41% | -3.35 个百分点 | +| 端到端耗时 | 12.34 秒 | 11.71 秒 | -5.11% | +| 峰值 RSS | 913 MiB | 932 MiB | +19 MiB | + +结果达到预定验收条件：P95 覆盖半径下降超过 15%，评分中位百分位下降少于 10 个百分点，运行时间没有增加，峰值内存增加远低于 256 MiB。 + +### 2026-08-03：跨连通区域共享全局间距 + +首次混合覆盖结果中，节点 `121302` 与 `110874` 分属两个不可达连通区域，但地图距离只有约 550.47 个模型单位。旧的分区独立首点规则分别选中了两个区域的最高分节点，形成视觉近邻。加入跨区域全局地图间距后，`121302` 被替换，最小归一化点间距由 0.013685 进一步提高到 0.040887。 + +本次调整保留连通区域名额，改变区域之间互不感知的选点方式。区域按名额从多到少处理，主干管网先形成覆盖骨架。第一个区域仍从综合评分最高的候选开始；后续区域选择首点时，先计算该区域所有候选到全局已选点的最近地图距离，只保留达到本区域最大空白距离 70% 的候选，再比较灵敏度综合评分。区域内部后续选点继续使用地图距离与管网路径距离的混合空白度。 + +这项约束只影响跨区域首点选择，不改变灵敏度估计、水力距离估计、区域名额、最小管径规则和公开 API。合成回归模型会构造两个地图上重叠但拓扑断开的区域，确保算法不会再次分别选择两个相邻高分点。 + +## 测试与回归验证 + +新增测试位于 `tests/unit/test_sensor_sensitivity.py`，覆盖以下行为： + +- 相同输入返回相同节点，并且每次调用只运行一次 EPANET。 +- EPANET 只保留初始状态，模型原始模拟时长能够恢复。 +- 关联矩阵和距离图保持 CSR 稀疏格式。 +- 最小管径按毫米过滤安装候选。 +- 近似方案在稠密参考目标上的结果不低于 95%。 +- 高密西部、低密东部的合成模型不再按候选节点密度分配名额。 +- 地理位置相邻但拓扑断开的区域在名额允许时分别获得监测点。 +- 地图上重叠的独立区域共享全局地理间距，不能各自选择相邻的首个高分点。 +- 名额不足时优先覆盖有效管道长度更大的连通区域。 +- 重复坐标依靠管网路径距离继续选点，并保持数量和确定性。 +- 非法监测点数量和最小管径返回明确错误。 +- 模拟结束后不残留共享 `temp.*` 文件。 + +回归命令： + +```bash +conda run -n server python -m pytest \ + tests/unit \ + tests/auth \ + tests/api/test_sensor_placement_endpoints.py \ + tests/api/test_simulation_endpoints.py \ + -q +``` + +2026-08-03 的执行结果为 `142 passed, 2 skipped, 7 warnings`。`git diff --check` 同时通过。 + +## 实现边界 + +- 256 次投影和 256 个代表点是当前质量与性能验证后的固定参数。调整参数需要重新运行稠密参考质量测试和大模型压力测试。 +- 稀疏 LU 对高连接度或拓扑特殊的模型可能产生更多填充，应继续用进程级资源保护运行未知大模型。 +- 算法使用单个初始水力状态。如果业务目标改为覆盖全天多个工况，需要先定义多工况评分和结果合并规则，不能直接恢复长时段模拟后沿用当前评分。 +- 节点坐标用于代表点构建和地图覆盖，假定其能表达一致的平面相对距离；缺少有效二维坐标的模型会返回错误。 +- 物理覆盖图使用初始水力状态。全天工况中频繁开闭的阀门或泵需要在多工况方案中重新定义连通区域合并规则。 +- 当前压力测试验证到 203,569 个原始节点，没有验证 30 万节点模型。 diff --git a/tests/unit/test_sensor_sensitivity.py b/tests/unit/test_sensor_sensitivity.py new file mode 100644 index 0000000..42a8a50 --- /dev/null +++ b/tests/unit/test_sensor_sensitivity.py @@ -0,0 +1,391 @@ +from pathlib import Path + +import numpy as np +import pytest +import wntr +from scipy.sparse import csr_matrix, isspmatrix_csr +from scipy.sparse.csgraph import dijkstra + +from app.algorithms.sensor import sensitivity + + +def _build_test_network() -> wntr.network.WaterNetworkModel: + wn = wntr.network.WaterNetworkModel() + wn.options.time.duration = 0 + wn.add_reservoir("R1", base_head=100.0, coordinates=(-1.0, 0.0)) + wn.add_junction("J0", elevation=5.0, coordinates=(0.0, 0.0)) + + for index in range(1, 13): + wn.add_junction( + f"J{index}", + base_demand=0.001 + index * 0.00001, + elevation=5.0 + index * 0.05, + coordinates=(float(index % 4), float(index // 4)), + ) + + wn.add_pipe("P0", "R1", "J0", length=100.0, diameter=0.4, roughness=110) + wn.add_pipe("P1", "J0", "J1", length=100.0, diameter=0.4, roughness=110) + for index in range(1, 11): + wn.add_pipe( + f"P{index + 1}", + f"J{index}", + f"J{index + 1}", + length=80.0 + index, + diameter=0.3, + roughness=105, + ) + # J12 is connected only through a small pipe, while J11 also touches P11. + wn.add_pipe("P12", "J11", "J12", length=90.0, diameter=0.1, roughness=105) + wn.add_pipe("PX1", "J2", "J6", length=120.0, diameter=0.3, roughness=105) + wn.add_pipe("PX2", "J5", "J9", length=120.0, diameter=0.3, roughness=105) + return wn + + +def _prepared_selection_network( + coordinates: np.ndarray, + edges: list[tuple[int, int, float]], +) -> sensitivity._PreparedNetwork: + node_count = len(coordinates) + rows: list[int] = [] + columns: list[int] = [] + weights: list[float] = [] + for start, end, weight in edges: + rows.extend((start, end)) + columns.extend((end, start)) + weights.extend((weight, weight)) + coverage_graph = csr_matrix( + (weights, (rows, columns)), + shape=(node_count, node_count), + ) + return sensitivity._PreparedNetwork( + node_names=tuple(f"N{index:04d}" for index in range(node_count)), + full_node_indices=np.arange(node_count, dtype=np.int64), + candidate_indices=np.arange(node_count, dtype=np.int64), + coordinates=np.asarray(coordinates, dtype=np.float64), + incidence=csr_matrix((node_count, 1), dtype=np.float64), + conductance=np.ones(1, dtype=np.float64), + roughness_response=np.ones(1, dtype=np.float64), + distance_graph=csr_matrix((node_count, node_count), dtype=np.float64), + coverage_graph=coverage_graph, + ) + + +def test_algorithm_is_deterministic_and_runs_epanet_once(monkeypatch, tmp_path): + wn = _build_test_network() + original_run_sim = wntr.sim.EpanetSimulator.run_sim + prefixes: list[str] = [] + + def counted_run_sim(simulator, *args, **kwargs): + prefixes.append(str(kwargs["file_prefix"])) + return original_run_sim(simulator, *args, **kwargs) + + monkeypatch.setattr(wntr.sim.EpanetSimulator, "run_sim", counted_run_sim) + monkeypatch.chdir(tmp_path) + + first = sensitivity.optimize_sensor_placement(wn, sensor_num=4, min_diameter=0) + second = sensitivity.optimize_sensor_placement(wn, sensor_num=4, min_diameter=0) + + assert first == second + assert len(first) == len(set(first)) == 4 + assert len(prefixes) == 2 + assert all(not Path(prefix).parent.exists() for prefix in prefixes) + assert not list(tmp_path.glob("temp.*")) + + +def test_hydraulic_simulation_keeps_only_initial_state_and_restores_duration(): + wn = _build_test_network() + wn.options.time.duration = 24 * 60 * 60 + + results = sensitivity._run_hydraulic_simulation(wn) + + assert len(results.node["head"].index) == 1 + assert wn.options.time.duration == 24 * 60 * 60 + + +def test_preparation_keeps_network_matrices_sparse(): + wn = _build_test_network() + results = sensitivity._run_hydraulic_simulation(wn) + + prepared = sensitivity._prepare_network(wn, results, min_diameter=0) + + assert isspmatrix_csr(prepared.incidence) + assert isspmatrix_csr(prepared.distance_graph) + assert isspmatrix_csr(prepared.coverage_graph) + assert prepared.incidence.nnz <= 2 * prepared.incidence.shape[1] + assert prepared.distance_graph.nnz <= wn.num_pipes + assert prepared.coverage_graph.nnz <= 2 * wn.num_links + assert (prepared.coverage_graph != prepared.coverage_graph.T).nnz == 0 + dense_incidence_bytes = int(np.prod(prepared.incidence.shape)) * 8 + sparse_payload_bytes = ( + prepared.incidence.data.nbytes + + prepared.incidence.indices.nbytes + + prepared.incidence.indptr.nbytes + ) + assert sparse_payload_bytes < dense_incidence_bytes + + +def test_minimum_diameter_filters_installation_candidates_in_millimetres(): + wn = _build_test_network() + results = sensitivity._run_hydraulic_simulation(wn) + prepared = sensitivity._prepare_network(wn, results, min_diameter=300) + candidate_names = { + prepared.node_names[index] for index in prepared.candidate_indices + } + + assert "J12" not in candidate_names + assert "J11" in candidate_names + + selected = sensitivity.optimize_sensor_placement( + wn, + sensor_num=4, + min_diameter=300, + ) + assert set(selected) <= candidate_names + + with pytest.raises(ValueError, match="候选节点少于"): + sensitivity.optimize_sensor_placement( + wn, + sensor_num=len(candidate_names) + 1, + min_diameter=300, + ) + + +def test_sparse_estimate_preserves_dense_reference_placement_quality(): + wn = _build_test_network() + results = sensitivity._run_hydraulic_simulation(wn) + prepared = sensitivity._prepare_network(wn, results, min_diameter=0) + + approximate_log_sensitivity = sensitivity._estimate_log_pressure_sensitivity( + prepared + ) + approximate_distance = sensitivity._estimate_hydraulic_distance_sums(prepared) + approximate_selected = sensitivity._select_sensor_nodes( + prepared, + approximate_log_sensitivity, + approximate_distance, + sensor_num=4, + ) + + incidence = prepared.incidence.toarray() + laplacian = ( + prepared.incidence.multiply(prepared.conductance) + @ prepared.incidence.T + ).toarray() + diagonal_scale = float(np.max(np.abs(np.diag(laplacian)))) + laplacian += np.eye(laplacian.shape[0]) * ( + diagonal_scale * np.sqrt(np.finfo(np.float64).eps) + ) + response = np.linalg.solve( + laplacian, + incidence * prepared.roughness_response, + ) + exact_sensitivity = np.abs(response).sum(axis=1) + + exact_distances = dijkstra( + prepared.distance_graph.transpose().tocsr(), + directed=True, + indices=prepared.full_node_indices, + return_predecessors=False, + )[:, prepared.full_node_indices] + exact_distances[~np.isfinite(exact_distances)] = 0.0 + exact_distance = exact_distances.sum(axis=0) + exact_selected = sensitivity._select_sensor_nodes( + prepared, + np.log(np.maximum(exact_sensitivity, np.finfo(np.float64).tiny)), + exact_distance, + sensor_num=4, + ) + + exact_score = exact_sensitivity * exact_distance + node_index = { + node_name: index for index, node_name in enumerate(prepared.node_names) + } + approximate_objective = sum( + exact_score[node_index[node_name]] for node_name in approximate_selected + ) + exact_objective = sum( + exact_score[node_index[node_name]] for node_name in exact_selected + ) + + assert approximate_objective / exact_objective >= 0.95 + + +def test_mixed_coverage_avoids_candidate_density_bias(): + dense_west = np.linspace(0.0, 2.0, 200) + sparse_east = np.linspace(3.0, 10.0, 20) + x_coordinates = np.concatenate((dense_west, sparse_east)) + coordinates = np.column_stack( + (x_coordinates, np.zeros(len(x_coordinates), dtype=np.float64)) + ) + ordered = np.argsort(x_coordinates) + edges = [ + ( + int(start), + int(end), + float(x_coordinates[end] - x_coordinates[start]), + ) + for start, end in zip(ordered[:-1], ordered[1:]) + ] + prepared = _prepared_selection_network(coordinates, edges) + log_scores = np.linspace(4.0, 0.0, len(coordinates)) + distance_sums = np.ones(len(coordinates), dtype=np.float64) + + selected = sensitivity._select_sensor_nodes( + prepared, + log_scores, + distance_sums, + sensor_num=6, + ) + name_to_position = { + name: position for position, name in enumerate(prepared.node_names) + } + selected_positions = [name_to_position[name] for name in selected] + new_metrics = sensitivity._geographic_coverage_metrics( + coordinates, + selected_positions, + ) + + legacy_labels, _centers = sensitivity._cluster_labels( + coordinates, + 6, + random_seed=sensitivity._RANDOM_SEED + 2, + ) + legacy_positions: list[int] = [] + represented: set[int] = set() + for position in np.argsort(-log_scores): + label = int(legacy_labels[position]) + if label in represented: + continue + represented.add(label) + legacy_positions.append(int(position)) + legacy_metrics = sensitivity._geographic_coverage_metrics( + coordinates, + legacy_positions, + ) + + assert new_metrics[0] <= legacy_metrics[0] * 0.6 + assert new_metrics[2] >= legacy_metrics[2] * 1.5 + assert max(x_coordinates[selected_positions]) >= 9.0 + + +def test_disconnected_components_each_receive_a_sensor_when_slots_allow(): + coordinates = np.asarray( + [ + (0.0, 0.0), + (1.0, 0.0), + (0.0, 0.01), + (1.0, 0.01), + ] + ) + prepared = _prepared_selection_network( + coordinates, + [(0, 1, 1.0), (2, 3, 1.0)], + ) + + selected = sensitivity._select_sensor_nodes( + prepared, + np.asarray([10.0, 9.0, 8.0, 7.0]), + np.ones(4), + sensor_num=2, + ) + + assert len(set(selected) & {"N0000", "N0001"}) == 1 + assert len(set(selected) & {"N0002", "N0003"}) == 1 + + +def test_overlapping_components_respect_global_geographic_spacing(): + coordinates = np.asarray( + [ + (0.0, 0.0), + (10.0, 0.0), + (0.1, 0.0), + (10.1, 0.0), + ] + ) + prepared = _prepared_selection_network( + coordinates, + [(0, 1, 10.0), (2, 3, 10.0)], + ) + + selected = sensitivity._select_sensor_nodes( + prepared, + np.asarray([10.0, 1.0, 9.0, 0.0]), + np.ones(4), + sensor_num=2, + ) + selected_positions = [prepared.node_names.index(name) for name in selected] + minimum_gap = sensitivity._geographic_coverage_metrics( + coordinates, + selected_positions, + )[2] + + assert minimum_gap >= 0.9 + + +def test_component_quota_prefers_longer_networks_when_slots_are_limited(): + coordinates = np.asarray( + [ + (0.0, 0.0), + (10.0, 0.0), + (20.0, 0.0), + (25.0, 0.0), + (30.0, 0.0), + (31.0, 0.0), + ] + ) + prepared = _prepared_selection_network( + coordinates, + [(0, 1, 10.0), (2, 3, 5.0), (4, 5, 1.0)], + ) + + selected = sensitivity._select_sensor_nodes( + prepared, + np.asarray([1.0, 1.0, 2.0, 2.0, 100.0, 100.0]), + np.ones(6), + sensor_num=2, + ) + + assert set(selected) <= {"N0000", "N0001", "N0002", "N0003"} + assert len(set(selected) & {"N0000", "N0001"}) == 1 + assert len(set(selected) & {"N0002", "N0003"}) == 1 + + +def test_duplicate_coordinates_use_topology_and_return_exact_count(): + coordinates = np.zeros((6, 2), dtype=np.float64) + prepared = _prepared_selection_network( + coordinates, + [(index, index + 1, 1.0) for index in range(5)], + ) + log_scores = np.linspace(6.0, 1.0, 6) + + first = sensitivity._select_sensor_nodes( + prepared, + log_scores, + np.ones(6), + sensor_num=4, + ) + second = sensitivity._select_sensor_nodes( + prepared, + log_scores, + np.ones(6), + sensor_num=4, + ) + + assert first == second + assert len(first) == len(set(first)) == 4 + + +@pytest.mark.parametrize( + ("sensor_num", "min_diameter", "message"), + [ + (0, 0, "监测点数量必须大于 0"), + (1, -1, "最小管径不能小于 0"), + ], +) +def test_algorithm_rejects_invalid_parameters(sensor_num, min_diameter, message): + with pytest.raises(ValueError, match=message): + sensitivity.optimize_sensor_placement( + _build_test_network(), + sensor_num=sensor_num, + min_diameter=min_diameter, + )