655 lines
30 KiB
Python
655 lines
30 KiB
Python
# 改进灵敏度法
|
||
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
|
||
import project_info
|
||
|
||
# 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=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)
|