917 lines
31 KiB
Python
917 lines
31 KiB
Python
"""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 wntr
|
|
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
|
|
|
|
|
|
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.
|
|
"""
|
|
|
|
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
|
|
|
|
|
|
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)
|
|
|
|
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]
|
|
|
|
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,
|
|
)
|
|
|
|
|
|
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
|
|
|
|
|
|
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],
|
|
)
|
|
)
|
|
quotas[active_components[order[:sensor_num]]] = 1
|
|
return candidate_components, quotas
|
|
|
|
quotas[active_components] = 1
|
|
remaining = sensor_num - len(active_components)
|
|
while remaining > 0:
|
|
available = active_components[
|
|
quotas[active_components] < capacities[active_components]
|
|
]
|
|
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],
|
|
)
|
|
whole_count = int(whole.sum())
|
|
if whole_count:
|
|
quotas[available] += whole
|
|
remaining -= whole_count
|
|
continue
|
|
|
|
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
|
|
|
|
return candidate_components, quotas
|
|
|
|
|
|
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.
|
|
"""
|
|
|
|
if sensor_num <= 0:
|
|
raise ValueError("监测点数量必须大于 0")
|
|
if min_diameter < 0:
|
|
raise ValueError("最小管径不能小于 0")
|
|
|
|
total_started = perf_counter()
|
|
simulation_started = total_started
|
|
results = _run_hydraulic_simulation(wn)
|
|
simulation_seconds = perf_counter() - simulation_started
|
|
|
|
preparation_started = perf_counter()
|
|
prepared = _prepare_network(wn, results, min_diameter=min_diameter)
|
|
preparation_seconds = perf_counter() - preparation_started
|
|
|
|
sensitivity_started = perf_counter()
|
|
log_sensitivity = _estimate_log_pressure_sensitivity(prepared)
|
|
sensitivity_seconds = perf_counter() - sensitivity_started
|
|
|
|
distance_started = perf_counter()
|
|
hydraulic_distance_sums = _estimate_hydraulic_distance_sums(prepared)
|
|
distance_seconds = perf_counter() - distance_started
|
|
|
|
selection_started = perf_counter()
|
|
selected = _select_sensor_nodes(
|
|
prepared,
|
|
log_sensitivity,
|
|
hydraulic_distance_sums,
|
|
sensor_num=sensor_num,
|
|
)
|
|
selection_seconds = perf_counter() - selection_started
|
|
|
|
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
|
|
|
|
|
|
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,
|
|
)
|
|
|
|
|
|
def get_ID(name: str, sensor_num: int, min_diameter: int) -> list[str]:
|
|
"""Compatibility entry point used by the sensor placement service."""
|
|
|
|
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,
|
|
)
|