diff --git a/sensitivity.py b/sensitivity.py new file mode 100644 index 0000000..3678548 --- /dev/null +++ b/sensitivity.py @@ -0,0 +1,654 @@ +# 改进灵敏度法 +import networkx +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 tjnetwork import * +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 + + +# 2025/03/12 +# Step1: 获取节点坐标 +def getCoor(wn: wntr.network.WaterNetworkModel) -> pandas.DataFrame: + """ + 获取管网模型的节点坐标 + :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 + + +# 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 skater_partition(G, n_clusters): + """ + 使用 SKATER 算法对输入的无向图 G 进行区域划分， + 保证每个划分区域在图论意义上是连通的， + 同时依据节点坐标的空间信息进行划分。 + + 参数: + 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) + 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): + + # 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)) + # 泵、起终点为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)) + + # === 改动新增部分：筛选管径小于 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()] + 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) + 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 + + def Jaco(self, hL: pandas.DataFrame): + """ + 计算灵敏度矩阵（节点压力对粗糙度变化的响应） + :param hL: 水力距离矩阵 + :return: + """ + # global result + # A：numpy.matrix, 节点-管道关系矩阵 + A = self.A2 + wn = self.wn + + 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 + + +# 2025/03/12 +# Step4: 传感器布置优化 +# Sensorplacement +# weight：分配权重 +# sensor：传感器布置的位置 +class Sensorplacement(wn_func): + """ + Sensorplacement 类继承了 wn_func 类，并且用于计算和优化传感器布置的位置。 + """ + def __init__(self, wn: wntr.network.WaterNetworkModel, sensornum: int, min_diameter: int): + """ + + :param wn: 由wntr生成的模型 + :param sensornum: 传感器的数量 + :param min_diameter: 安装的最小管径 + """ + wn_func.__init__(self, wn, min_diameter=min_diameter) + self.sensornum = sensornum + + # 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') # 存储节点总灵敏度 + + # 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 = {} + + # 改动 + 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 + + # 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]) + + for key, value in group.items(): + if Smaxnode in value: + sumSS = sumSS.drop(index=group[key]) + continue + + sumSS = sumSS.sort_values(by=[0], ascending=[False]) + + return sensorindex, sensorindex_2 + + +# 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() + + # 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='bb', 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)