返回完整的代码流程，便于我运行

<think>我们按照要求修改代码，将喷头之间连接成网状（使用最小生成树），然后通过一根水管连接河流。储水罐也采用最小生成树连接，并且只从河流连接一根水管到最近的储水罐。 修改点： 1. 喷头网络：所有喷头通过最小生成树（MST）连接，并选择离河流最近的喷头作为整个喷头网络的入口（即连接到河流的取水点）。 2. 储水罐网络：储水罐之间也通过MST连接，并选择离河流最近的储水罐作为整个储水罐网络的入口（即河流水管连接到这个储水罐）。 3. 优化目标函数：成本计算将基于网络结构（管道长度和流量）进行。 注意：由于喷头网络和储水罐网络是分开的，我们需要将储水罐连接到喷头网络。这里，每个储水罐连接到其负责的喷头集群中最近的一个喷头（作为接入点），然后通过喷头网络将水分配到各个喷头。 具体步骤： a) 生成喷头布局，并计算每个喷头的覆盖面积和最大灌溉需求。 b) 构建喷头网络的最小生成树（MST），并计算每条边的流量（基于下游喷头的总需求）。 c) 确定喷头网络的入口（离河流最近的喷头）。 d) 对储水罐进行聚类（数量由肘部法确定），并构建储水罐的MST，确定储水罐网络的入口（离河流最近的储水罐）。 e) 将每个储水罐连接到其负责的喷头集群中最近的一个喷头（即接入喷头网络）。 f) 优化目标函数包括： - 储水罐成本（与容量成正比） - 河流到储水罐入口的管道成本（一根水管） - 储水罐之间管道的成本（MST中的边） - 储水罐到喷头接入点的管道成本（每个储水罐到其接入喷头的管道） - 喷头网络管道的成本（MST中的边，流量为下游喷头的总需求） 由于代码较长，我们将修改关键函数：`two_stage_optimization`，并添加新的函数（如计算喷头网络MST和流量分配）。同时，修改可视化函数以展示网络结构。 以下是修改后的完整代码（注意：由于问题复杂，我们只修改关键部分，其他部分如数据加载、喷头生成等保持不变）：</think>以下是根据您的要求修改后的完整代码流程，已实现喷头之间网状连接和储水罐最短连接优化： ```python import numpy as np import pandas as pd from scipy.optimize import minimize, Bounds from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score import matplotlib.pyplot as plt import warnings import logging import math from matplotlib.patches import Circle, Rectangle from sklearn.preprocessing import StandardScaler import os import networkx as nx from scipy.sparse.csgraph import minimum_spanning_tree # 设置日志记录 logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s') logger = logging.getLogger(__name__) # 忽略特定警告 warnings.filterwarnings("ignore", category=UserWarning, module="sklearn") warnings.filterwarnings("ignore", category=RuntimeWarning) # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei'] # 使用黑体 plt.rcParams['axes.unicode_minus'] = False # 正常显示负号 # ==================== 常量定义 ==================== DEPTH = 0.05 # 5 cm DRY_DENS = 1500 # kg/m³ dry_mass_per_m2 = DRY_DENS * DEPTH # 75 kg/m² MIN_SOIL_MOISTURE = 0.22 # 最低土壤湿度 # 农场尺寸 (1公顷正方形) FARM_SIZE = 100 # 米 FARM_SIZE_X, FARM_SIZE_Y = FARM_SIZE, FARM_SIZE logger.info(f"农场尺寸: {FARM_SIZE_X}m × {FARM_SIZE_Y}m") # 作物信息 CROPS = { '高粱': {'area': 0.5 * 10000}, # 公顷转换为平方米 '玉米': {'area': 0.3 * 10000}, '大豆': {'area': 0.2 * 10000} } # 作物分布 CROP_REGIONS = { '高粱': {'x_min': 0, 'x_max': FARM_SIZE_X, 'y_min': 0, 'y_max': FARM_SIZE_Y/2}, '玉米': {'x_min': 0, 'x_max': 60, 'y_min': FARM_SIZE_Y/2, 'y_max': FARM_SIZE_Y}, '大豆': {'x_min': 60, 'x_max': FARM_SIZE_X, 'y_min': FARM_SIZE_Y/2, 'y_max': FARM_SIZE_Y} } SPRINKLER_RADIUS = 15 # 喷头半径 RIVER_POSITION = 'south' RIVER_POINT = (FARM_SIZE_X / 2, 0) # 河流取水点（南侧中心） # 成本公式参数 PIPE_LENGTH_COEFF = 50 # 管道长度系数 PIPE_FLOW_COEFF = 0.1 # 管道流量系数 PIPE_LENGTH_EXP = 1.2 # 管道长度指数 PIPE_FLOW_EXP = 1.5 # 管道流量指数 TANK_COST_PER_LITER = 5 # 储水罐单位容积成本 # 单位转换系数 L_TO_M3 = 0.001 # 1L = 0.001m³ # 系统参数 DAILY_WATER_SOURCE_RATIO = 0.8 # 日常水源中河水的比例 EMERGENCY_WATER_SOURCE_RATIO = 0.2 # 应急水源比例 # ==================== 数据加载与处理 ==================== def load_soil_moisture_data(): """从Excel文件加载真实的土壤湿度数据""" try: file_path = '附件/该地土壤湿度数据.xlsx' if not os.path.exists(file_path): logger.error(f"土壤湿度数据文件不存在: {file_path}") dates = pd.date_range('2021-07-01', periods=31) moisture_values = 0.15 + 0.1 * np.sin(np.linspace(0, 2*np.pi, 31)) daily_avg_moisture = pd.Series(moisture_values, index=dates) logger.warning("使用模拟数据替代") return daily_avg_moisture logger.info(f"从Excel文件加载土壤湿度数据: {file_path}") data = pd.read_excel(file_path, sheet_name='JingYueTan') required_columns = ['DATE', '5cm_SM'] if not all(col in data.columns for col in required_columns): logger.error(f"Excel文件中缺少必要的列: {required_columns}") dates = pd.date_range('2021-07-01', periods=31) moisture_values = 0.15 + 0.1 * np.sin(np.linspace(0, 2*np.pi, 31)) daily_avg_moisture = pd.Series(moisture_values, index=dates) logger.warning("使用模拟数据替代") return daily_avg_moisture data['DATE'] = pd.to_datetime(data['DATE']) data.set_index('DATE', inplace=True) start_date = pd.Timestamp('2021-07-01') end_date = pd.Timestamp('2021-07-31') july_data = data.loc[(data.index >= start_date) & (data.index <= end_date)] if july_data.empty: logger.warning("2021年7月数据为空，使用全年数据") july_data = data.copy() july_data.sort_index(inplace=True) plt.figure(figsize=(12, 6)) plt.plot(july_data.index, july_data['5cm_SM'].values, 'b-', linewidth=2) plt.axhline(y=MIN_SOIL_MOISTURE, color='r', linestyle='--', label='最低土壤湿度阈值') plt.title('2021年7月土壤湿度变化') plt.xlabel('日期') plt.ylabel('土壤湿度 (5cm_SM)') plt.legend() plt.grid(True) plt.xticks(rotation=45) plt.tight_layout() plt.savefig('土壤湿度变化图.png', dpi=300) plt.show() return july_data['5cm_SM'] except Exception as e: logger.error(f"加载土壤湿度数据时出错: {e}") dates = pd.date_range('2021-07-01', periods=31) moisture_values = 0.15 + 0.1 * np.sin(np.linspace(0, 2*np.pi, 31)) daily_avg_moisture = pd.Series(moisture_values, index=dates) logger.warning("使用模拟数据替代") return daily_avg_moisture def calculate_daily_irrigation_demand(daily_moisture): """计算每日每平方米灌溉需求""" return max(0.0, (MIN_SOIL_MOISTURE - daily_moisture) * dry_mass_per_m2) def calculate_irrigation_demand(daily_avg_moisture, sprinkler_df): """计算每日灌溉需求""" max_demand_per_m2 = max(calculate_daily_irrigation_demand(m) for m in daily_avg_moisture) sprinkler_df['max_demand'] = sprinkler_df['area'] * max_demand_per_m2 plt.figure(figsize=(12, 6)) daily_demand = [calculate_daily_irrigation_demand(m) for m in daily_avg_moisture] plt.plot(daily_avg_moisture.index, daily_demand, label='灌溉需求') plt.title('2021年7月每日灌溉需求') plt.xlabel('日期') plt.ylabel('灌溉需求 (L/m²)') plt.legend() plt.grid(True) plt.xticks(rotation=45) plt.tight_layout() plt.savefig('灌溉需求变化图.png', dpi=300) plt.show() return daily_demand, sprinkler_df # ==================== 喷头布局生成 ==================== def generate_sprinkler_layout(farm_size_x=FARM_SIZE_X, farm_size_y=FARM_SIZE_Y, radius=SPRINKLER_RADIUS): """生成六角形喷头布局""" spacing = radius * 1.5 sprinklers = [] rows = int(farm_size_y / (spacing * np.sqrt(3)/2)) + 2 cols = int(farm_size_x / spacing) + 2 for i in range(rows): for j in range(cols): x = j * spacing y = i * spacing * np.sqrt(3)/2 if i % 2 == 1: x += spacing / 2 if 0 <= x <= farm_size_x and 0 <= y <= farm_size_y: crop_type = None for crop, region in CROP_REGIONS.items(): if region['x_min'] <= x <= region['x_max'] and region['y_min'] <= y <= region['y_max']: crop_type = crop break sprinklers.append({ 'id': len(sprinklers), 'x': x, 'y': y, 'radius': radius, 'crop_type': crop_type }) return pd.DataFrame(sprinklers) def calculate_circle_segment_area(radius, overlap): """计算圆形边界重叠部分的面积""" angle = 2 * math.acos(overlap / radius) segment_area = (radius**2) * (angle - math.sin(angle)) / 2 return segment_area def calculate_sprinkler_coverage(sprinkler_df, farm_size_x=FARM_SIZE_X, farm_size_y=FARM_SIZE_Y): """计算每个喷头的覆盖面积""" full_area = np.pi * SPRINKLER_RADIUS ** 2 areas = [] for _, sprinkler in sprinkler_df.iterrows(): x, y = sprinkler['x'], sprinkler['y'] effective_area = full_area if x < SPRINKLER_RADIUS: overlap = SPRINKLER_RADIUS - x segment_area = calculate_circle_segment_area(SPRINKLER_RADIUS, overlap) effective_area -= segment_area if x > farm_size_x - SPRINKLER_RADIUS: overlap = SPRINKLER_RADIUS - (farm_size_x - x) segment_area = calculate_circle_segment_area(SPRINKLER_RADIUS, overlap) effective_area -= segment_area if y < SPRINKLER_RADIUS: overlap = SPRINKLER_RADIUS - y segment_area = calculate_circle_segment_area(SPRINKLER_RADIUS, overlap) effective_area -= segment_area if y > farm_size_y - SPRINKLER_RADIUS: overlap = SPRINKLER_RADIUS - (farm_size_y - y) segment_area = calculate_circle_segment_area(SPRINKLER_RADIUS, overlap) effective_area -= segment_area areas.append(effective_area) sprinkler_df['area'] = areas return sprinkler_df def validate_sprinkler_spacing(sprinkler_df, min_spacing=15): """验证喷头间距是否≥15m""" points = sprinkler_df[['x', 'y']].values num_sprinklers = len(points) min_distance = float('inf') min_pair = (-1, -1) for i in range(num_sprinklers): for j in range(i+1, num_sprinklers): dist = np.sqrt((points[i][0] - points[j][0])**2 + (points[i][1] - points[j][1])**2) if dist < min_distance: min_distance = dist min_pair = (i, j) plt.figure(figsize=(12, 10)) plt.scatter(sprinkler_df['x'], sprinkler_df['y'], c='blue', s=50, label='喷头') if min_pair != (-1, -1): plt.plot([points[min_pair[0]][0], points[min_pair[1]][0]], [points[min_pair[0]][1], points[min_pair[1]][1]], 'r--', linewidth=2, label=f'最小间距: {min_distance:.2f}m') for i, row in sprinkler_df.iterrows(): plt.text(row['x'], row['y'], f"S{i}", fontsize=9, ha='center', va='bottom') plt.text(row['x'], row['y'], f"({row['x']:.1f},{row['y']:.1f})", fontsize=8, ha='center', va='top') colors = {'高粱': 'lightgreen', '玉米': 'lightyellow', '大豆': 'lightblue'} for crop, region in CROP_REGIONS.items(): width = region['x_max'] - region['x_min'] height = region['y_max'] - region['y_min'] rect = plt.Rectangle((region['x_min'], region['y_min']), width, height, alpha=0.3, color=colors[crop], label=crop) plt.gca().add_patch(rect) plt.plot([0, FARM_SIZE_X], [0, 0], 'b-', linewidth=4, label='河流') plt.title(f'喷头布局图 (最小间距: {min_distance:.2f}m)') plt.xlabel('X坐标 (m)') plt.ylabel('Y坐标 (m)') plt.grid(True) plt.legend() plt.tight_layout() plt.savefig('喷头布局验证图.png', dpi=300) plt.show() if min_distance >= min_spacing: logger.info(f"喷头间距验证通过! 最小间距: {min_distance:.2f}m ≥ {min_spacing}m") return True else: logger.warning(f"喷头间距验证失败! 最小间距: {min_distance:.2f}m < {min_spacing}m") return False # ==================== 网络连接优化 ==================== def calculate_minimum_spanning_tree(points): """计算点集的最小生成树""" n = len(points) dist_matrix = np.zeros((n, n)) for i in range(n): for j in range(n): if i != j: dist_matrix[i, j] = np.sqrt((points[i][0] - points[j][0])**2 + (points[i][1] - points[j][1])**2) mst = minimum_spanning_tree(dist_matrix) edges = [] rows, cols = mst.nonzero() for i, j in zip(rows, cols): if i < j: edges.append((i, j, dist_matrix[i, j])) return edges def calculate_network_flows(sprinkler_df, root_idx, mst_edges): """计算喷头网络中每条边的流量""" G = nx.Graph() n = len(sprinkler_df) for i in range(n): G.add_node(i, demand=sprinkler_df.iloc[i]['max_demand']) for i, j, L in mst_edges: G.add_edge(i, j, length=L) edge_flows = {} visited = set() def dfs(node): visited.add(node) total_flow = sprinkler_df.iloc[node]['max_demand'] for neighbor in G.neighbors(node): if neighbor not in visited: subtree_flow = dfs(neighbor) total_flow += subtree_flow edge_flows[(min(node, neighbor), max(node, neighbor))] = subtree_flow return total_flow total_flow = dfs(root_idx) return edge_flows def determine_optimal_clusters(sprinkler_df, max_clusters=10): """确定最佳聚类数量""" coordinates = sprinkler_df[['x', 'y']].values scaler = StandardScaler() scaled_features = scaler.fit_transform(coordinates) sse = [] silhouette_scores = [] k_range = range(2, max_clusters + 1) for k in k_range: kmeans = KMeans(n_clusters=k, random_state=42, n_init=10) kmeans.fit(scaled_features) sse.append(kmeans.inertia_) if k > 1: silhouette_scores.append(silhouette_score(scaled_features, kmeans.labels_)) else: silhouette_scores.append(0) plt.figure(figsize=(12, 5)) plt.subplot(1, 2, 1) plt.plot(k_range, sse, 'bo-') plt.xlabel('聚类数量') plt.ylabel('SSE') plt.title('肘部法') plt.grid(True) plt.subplot(1, 2, 2) plt.plot(k_range[1:], silhouette_scores[1:], 'ro-') plt.xlabel('聚类数量') plt.ylabel('轮廓系数') plt.title('轮廓系数法') plt.grid(True) plt.tight_layout() plt.savefig('聚类数量选择.png', dpi=300) plt.show() sse_diff = np.diff(sse) sse_ratio = sse_diff[:-1] / sse_diff[1:] elbow_point = np.argmax(sse_ratio) + 2 best_silhouette = np.argmax(silhouette_scores[1:]) + 2 optimal_clusters = min(max(3, elbow_point), max(3, best_silhouette)) logger.info(f"肘部法建议聚类数量: {elbow_point}") logger.info(f"轮廓系数法建议聚类数量: {best_silhouette}") logger.info(f"最终选择聚类数量: {optimal_clusters}") return optimal_clusters def generate_candidate_tanks(sprinkler_df, num_tanks): """生成候选储水罐位置""" coordinates = sprinkler_df[['x', 'y']].values scaler = StandardScaler() scaled_features = scaler.fit_transform(coordinates) kmeans = KMeans(n_clusters=num_tanks, random_state=42, n_init=10) kmeans.fit(scaled_features) cluster_centers = scaler.inverse_transform(kmeans.cluster_centers_) return cluster_centers, kmeans.labels_ def calculate_distance_matrix(points1, points2): """计算两点集之间的距离矩阵""" num_points1 = len(points1) num_points2 = len(points2) distances = np.zeros((num_points1, num_points2)) for i, point1 in enumerate(points1): for j, point2 in enumerate(points2): dist = np.sqrt((point1[0] - point2[0])**2 + (point1[1] - point2[1])**2) distances[i, j] = dist return distances def two_stage_optimization(sprinkler_df, irrigation_demand, tank_positions): """ 两阶段优化：喷头网络连接 + 储水罐网络连接 """ # 阶段1: 构建喷头网络 sprinkler_points = sprinkler_df[['x', 'y']].values sprinkler_mst = calculate_minimum_spanning_tree(sprinkler_points) # 找到离河流最近的喷头作为网络入口 distances_to_river = [np.sqrt((x - RIVER_POINT[0])**2 + (y - RIVER_POINT[1])**2) for x, y in sprinkler_points] root_sprinkler_idx = np.argmin(distances_to_river) # 计算喷头网络中各边的流量 edge_flows = calculate_network_flows(sprinkler_df, root_sprinkler_idx, sprinkler_mst) # 阶段2: 分配喷头到储水罐 distances = calculate_distance_matrix( sprinkler_df[['x', 'y']].values, np.array(tank_positions) ) num_sprinklers = len(sprinkler_df) num_tanks = len(tank_positions) # 分配喷头到最近的储水罐 assignments = np.argmin(distances, axis=1) # 计算每个储水罐的需求 sprinkler_max_demands = sprinkler_df['max_demand'].values tank_demands = [] for j in range(num_tanks): demand = np.sum(sprinkler_max_demands[assignments == j]) tank_demands.append(demand) # 构建储水罐网络 tank_mst = calculate_minimum_spanning_tree(tank_positions) # 找到离河流最近的储水罐作为入口 distances_to_river = [np.sqrt((x - RIVER_POINT[0])**2 + (y - RIVER_POINT[1])**2) for x, y in tank_positions] root_tank_idx = np.argmin(distances_to_river) # 优化目标函数 def objective(vars): # 解析变量 V = vars[:num_tanks] # 储水罐容量（L） Q_river = vars[num_tanks:2*num_tanks] # 从河流到储水罐的流量（L/day） Q_tank = vars[2*num_tanks:3*num_tanks] # 从储水罐到喷头网络的流量（L/day） # 储水罐成本 tank_cost = TANK_COST_PER_LITER * np.sum(V) # 河流到储水罐的管道成本（只连接到根储水罐） river_pipe_cost = 0 L = distances_to_river[root_tank_idx] Q_m3 = Q_river[root_tank_idx] * L_TO_M3 river_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP) # 储水罐之间的管道成本 tank_pipe_cost = 0 for i, j, L in tank_mst: Q_avg = (Q_river[i] + Q_river[j]) / 2 Q_m3 = Q_avg * L_TO_M3 tank_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP) # 喷头网络管道成本 sprinkler_pipe_cost = 0 for (i, j), flow in edge_flows.items(): L = np.sqrt((sprinkler_points[i][0] - sprinkler_points[j][0])**2 + (sprinkler_points[i][1] - sprinkler_points[j][1])**2) Q_m3 = flow * L_TO_M3 sprinkler_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP) # 储水罐到喷头网络的管道成本 for j in range(num_tanks): tank_pos = tank_positions[j] cluster_sprinklers = sprinkler_df[assignments == j] if not cluster_sprinklers.empty: dists = np.sqrt((cluster_sprinklers['x'] - tank_pos[0])**2 + (cluster_sprinklers['y'] - tank_pos[1])**2) connection_idx = dists.idxmin() connection_point = sprinkler_df.loc[connection_idx][['x', 'y']].values L = np.sqrt((tank_pos[0] - connection_point[0])**2 + (tank_pos[1] - connection_point[1])**2) Q_m3 = Q_tank[j] * L_TO_M3 sprinkler_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP) # 惩罚项 penalty = 0 for j in range(num_tanks): total_demand = tank_demands[j] * (1 + EMERGENCY_WATER_SOURCE_RATIO) if V[j] < total_demand: penalty += 1000 * (total_demand - V[j]) if Q_river[j] < tank_demands[j] * DAILY_WATER_SOURCE_RATIO: penalty += 1000 * (tank_demands[j] * DAILY_WATER_SOURCE_RATIO - Q_river[j]) if Q_tank[j] < tank_demands[j]: penalty += 1000 * (tank_demands[j] - Q_tank[j]) return tank_cost + river_pipe_cost + tank_pipe_cost + sprinkler_pipe_cost + penalty # 约束条件 constraints = [] # 初始值 initial_V = [tank_demands[j] * (1 + EMERGENCY_WATER_SOURCE_RATIO) for j in range(num_tanks)] initial_Q_river = [tank_demands[j] * DAILY_WATER_SOURCE_RATIO for j in range(num_tanks)] initial_Q_tank = [tank_demands[j] for j in range(num_tanks)] x0 = initial_V + initial_Q_river + initial_Q_tank # 边界 bounds = Bounds([0] * (3 * num_tanks), [np.inf] * (3 * num_tanks)) # 优化 logger.info("开始优化...") result = minimize( objective, x0, bounds=bounds, constraints=constraints, method='SLSQP', options={'disp': True, 'ftol': 1e-6, 'maxiter': 100} ) return result, assignments, tank_demands, tank_mst, sprinkler_mst, root_sprinkler_idx, root_tank_idx # ==================== 成本计算与可视化 ==================== def calculate_total_cost(result, sprinkler_df, tank_positions, assignments, tank_mst): num_tanks = len(tank_positions) V_opt = result.x[:num_tanks] Q_river_opt = result.x[num_tanks:2*num_tanks] Q_tank_opt = result.x[2*num_tanks:3*num_tanks] # 储水罐成本 tank_cost = TANK_COST_PER_LITER * np.sum(V_opt) # 河流到储水罐的管道成本 river_pipe_cost = 0 distances_to_river = [np.sqrt((x - RIVER_POINT[0])**2 + (y - RIVER_POINT[1])**2) for x, y in tank_positions] root_tank_idx = np.argmin(distances_to_river) L = distances_to_river[root_tank_idx] Q_m3 = Q_river_opt[root_tank_idx] * L_TO_M3 river_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP) # 储水罐之间的管道成本 tank_pipe_cost = 0 for i, j, L in tank_mst: Q_avg = (Q_river_opt[i] + Q_river_opt[j]) / 2 Q_m3 = Q_avg * L_TO_M3 tank_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP) # 喷头网络管道成本 sprinkler_pipe_cost = 0 sprinkler_points = sprinkler_df[['x', 'y']].values sprinkler_mst = calculate_minimum_spanning_tree(sprinkler_points) root_sprinkler_idx = np.argmin([np.sqrt((x - RIVER_POINT[0])**2 + (y - RIVER_POINT[1])**2) for x, y in sprinkler_points]) edge_flows = calculate_network_flows(sprinkler_df, root_sprinkler_idx, sprinkler_mst) for (i, j), flow in edge_flows.items(): L = np.sqrt((sprinkler_points[i][0] - sprinkler_points[j][0])**2 + (sprinkler_points[i][1] - sprinkler_points[j][1])**2) Q_m3 = flow * L_TO_M3 sprinkler_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP) # 储水罐到喷头网络的管道成本 for j in range(num_tanks): tank_pos = tank_positions[j] cluster_sprinklers = sprinkler_df[assignments == j] if not cluster_sprinklers.empty: dists = np.sqrt((cluster_sprinklers['x'] - tank_pos[0])**2 + (cluster_sprinklers['y'] - tank_pos[1])**2) connection_idx = dists.idxmin() connection_point = sprinkler_df.loc[connection_idx][['x', 'y']].values L = np.sqrt((tank_pos[0] - connection_point[0])**2 + (tank_pos[1] - connection_point[1])**2) Q_m3 = Q_tank_opt[j] * L_TO_M3 sprinkler_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP) total_cost = tank_cost + river_pipe_cost + tank_pipe_cost + sprinkler_pipe_cost cost_breakdown = { 'tank_cost': tank_cost, 'river_pipe_cost': river_pipe_cost, 'tank_pipe_cost': tank_pipe_cost, 'sprinkler_pipe_cost': sprinkler_pipe_cost } return total_cost, cost_breakdown, Q_river_opt, Q_tank_opt def visualize_network_system(sprinkler_df, tank_positions, assignments, tank_mst, sprinkler_mst, root_sprinkler_idx, root_tank_idx, Q_river_opt, Q_tank_opt): """可视化网络连接的灌溉系统""" plt.figure(figsize=(16, 14)) ax = plt.gca() # 绘制农场边界 ax.plot([0, FARM_SIZE_X, FARM_SIZE_X, 0, 0], [0, 0, FARM_SIZE_Y, FARM_SIZE_Y, 0], 'k-', linewidth=2) # 绘制河流（南侧） ax.plot([0, FARM_SIZE_X], [0, 0], 'b-', linewidth=4, label='河流') ax.plot(RIVER_POINT[0], RIVER_POINT[1], 'bo', markersize=10, label='取水点') # 绘制作物区域 colors = {'高粱': 'lightgreen', '玉米': 'lightyellow', '大豆': 'lightblue'} for crop, region in CROP_REGIONS.items(): width = region['x_max'] - region['x_min'] height = region['y_max'] - region['y_min'] rect = Rectangle((region['x_min'], region['y_min']), width, height, alpha=0.2, color=colors[crop], label=crop) ax.add_patch(rect) # 绘制喷头 sprinkler_points = sprinkler_df[['x', 'y']].values for i, sprinkler in sprinkler_df.iterrows(): if i == root_sprinkler_idx: ax.plot(sprinkler['x'], sprinkler['y'], 'o', color='red', markersize=8, markeredgecolor='black', markeredgewidth=1.5, label='喷头网络入口') else: ax.plot(sprinkler['x'], sprinkler['y'], 'o', color='blue', markersize=6) ax.add_patch(Circle((sprinkler['x'], sprinkler['y']), SPRINKLER_RADIUS, color='blue', alpha=0.1)) # 绘制喷头之间的连接 for i, j, L in sprinkler_mst: ax.plot([sprinkler_points[i][0], sprinkler_points[j][0]], [sprinkler_points[i][1], sprinkler_points[j][1]], 'b-', linewidth=1.5, alpha=0.7, label='喷头间管道' if i == 0 and j == 1 else "") # 绘制储水罐 for j, tank in enumerate(tank_positions): if j == root_tank_idx: ax.plot(tank[0], tank[1], 's', color='red', markersize=12, markeredgecolor='black', markeredgewidth=2, label='储水罐网络入口') else: ax.plot(tank[0], tank[1], 's', color='purple', markersize=10, markeredgecolor='black', markeredgewidth=1.5) ax.add_patch(Circle((tank[0], tank[1]), SPRINKLER_RADIUS*2, color='purple', alpha=0.15)) ax.text(tank[0], tank[1], f'T{j+1}', fontsize=10, ha='center', va='center', fontweight='bold') # 绘制储水罐之间的连接 for i, j, L in tank_mst: ax.plot([tank_positions[i][0], tank_positions[j][0]], [tank_positions[i][1], tank_positions[j][1]], 'g-', linewidth=2, label='储水罐间管道' if i == 0 and j == 1 else "") mid_x = (tank_positions[i][0] + tank_positions[j][0]) / 2 mid_y = (tank_positions[i][1] + tank_positions[j][1]) / 2 ax.text(mid_x, mid_y, f'{L:.1f}m', fontsize=8, bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.7)) # 绘制河流到储水罐的管道 ax.plot([RIVER_POINT[0], tank_positions[root_tank_idx][0]], [RIVER_POINT[1], tank_positions[root_tank_idx][1]], 'r-', linewidth=3, label='河流到储水罐管道') mid_x = (RIVER_POINT[0] + tank_positions[root_tank_idx][0]) / 2 mid_y = (RIVER_POINT[1] + tank_positions[root_tank_idx][1]) / 2 length = np.sqrt((RIVER_POINT[0]-tank_positions[root_tank_idx][0])**2 + (RIVER_POINT[1]-tank_positions[root_tank_idx][1])**2) Q_m3 = Q_river_opt[root_tank_idx] * L_TO_M3 ax.text(mid_x, mid_y, f'{length:.1f}m\n{Q_river_opt[root_tank_idx]:.0f}L/d\n({Q_m3:.2f}m³/d)', fontsize=9, fontweight='bold', bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.8)) # 绘制储水罐到喷头网络的连接 for j in range(len(tank_positions)): cluster_sprinklers = sprinkler_df[assignments == j] if not cluster_sprinklers.empty: tank_pos = tank_positions[j] dists = np.sqrt((cluster_sprinklers['x'] - tank_pos[0])**2 + (cluster_sprinklers['y'] - tank_pos[1])**2) connection_idx = dists.idxmin() connection_point = sprinkler_df.loc[connection_idx][['x', 'y']].values ax.plot([tank_pos[0], connection_point[0]], [tank_pos[1], connection_point[1]], 'm--', linewidth=2, label='储水罐到喷头网络' if j == 0 else "") ax.set_xlabel('X坐标 (m)') ax.set_ylabel('Y坐标 (m)') ax.set_title('网络化灌溉系统优化布局') handles, labels = ax.get_legend_handles_labels() by_label = dict(zip(labels, handles)) ax.legend(by_label.values(), by_label.keys(), loc='best') ax.grid(True) plt.tight_layout() plt.savefig('网络化灌溉系统布局.png', dpi=300) plt.show() def plot_cost_breakdown(cost_breakdown): """绘制成本分解图""" labels = ['储水罐成本', '河流到储水罐管道成本', '储水罐间管道成本', '喷头网络管道成本'] sizes = [ cost_breakdown['tank_cost'], cost_breakdown['river_pipe_cost'], cost_breakdown['tank_pipe_cost'], cost_breakdown['sprinkler_pipe_cost'] ] colors = ['lightblue', 'lightcoral', 'lightgreen', 'lightyellow'] plt.figure(figsize=(10, 8)) plt.pie(sizes, labels=labels, colors=colors, autopct='%1.1f%%', startangle=90) plt.axis('equal') plt.title('成本分解') plt.savefig('成本分解图.png', dpi=300) plt.show() def output_results(sprinkler_df, tank_positions, result, assignments, tank_mst, Q_river_opt, Q_tank_opt, tank_demands, cost_breakdown): """输出优化结果表格""" num_tanks = len(tank_positions) V_opt = result.x[:num_tanks] results_df = pd.DataFrame({ '储水罐编号': range(1, num_tanks+1), 'X坐标': [f"{tank[0]:.1f}" for tank in tank_positions], 'Y坐标': [f"{tank[1]:.1f}" for tank in tank_positions], '容量(L)': [f"{v:.0f}" for v in V_opt], '容量(m³)': [f"{v * L_TO_M3:.2f}" for v in V_opt], '需求(L/天)': [f"{d:.0f}" for d in tank_demands], '需求(m³/天)': [f"{d * L_TO_M3:.2f}" for d in tank_demands], '河流到储水罐流量(L/天)': [f"{q:.0f}" for q in Q_river_opt], '河流到储水罐流量(m³/天)': [f"{q * L_TO_M3:.2f}" for q in Q_river_opt], '储水罐到喷头流量(L/天)': [f"{q:.0f}" for q in Q_tank_opt], '储水罐到喷头流量(m³/天)': [f"{q * L_TO_M3:.2f}" for q in Q_tank_opt] }) coverage_stats = [] for j in range(num_tanks): covered = np.sum(assignments == j) coverage_stats.append(covered) results_df['覆盖喷头数'] = coverage_stats tank_pipe_info = [] for i, j, L in tank_mst: tank_pipe_info.append(f'T{i+1}-T{j+1}: {L:.1f}m') results_df['储水罐间管道'] = [', '.join(tank_pipe_info)] * num_tanks logger.info("\n优化结果详情:") print(results_df.to_string(index=False)) results_df.to_csv('灌溉系统优化结果.csv', index=False, encoding='utf-8-sig') logger.info("结果已保存到 '灌溉系统优化结果.csv'") total_demand = np.sum(sprinkler_df['max_demand']) river_supply = np.sum(Q_river_opt) tank_supply = np.sum(V_opt) logger.info(f"\n系统总体信息:") logger.info(f"总灌溉需求: {total_demand:.0f} L/天 ({total_demand * L_TO_M3:.2f} m³/天)") logger.info(f"河流供水能力: {river_supply:.0f} L/天 ({river_supply * L_TO_M3:.2f} m³/天)") logger.info(f"储水罐总容量: {tank_supply:.0f} L ({tank_supply * L_TO_M3:.2f} m³)") logger.info(f"系统可靠性: {(river_supply + tank_supply)/total_demand*100:.1f}%") logger.info(f"\n成本详情:") logger.info(f"储水罐成本: {cost_breakdown['tank_cost']:.2f} 元") logger.info(f"河流到储水罐管道成本: {cost_breakdown['river_pipe_cost']:.2f} 元") logger.info(f"储水罐间管道成本: {cost_breakdown['tank_pipe_cost']:.2f} 元") logger.info(f"喷头网络管道成本: {cost_breakdown['sprinkler_pipe_cost']:.2f} 元") total_cost = sum(cost_breakdown.values()) logger.info(f"总成本: {total_cost:.2f} 元") # ==================== 主函数 ==================== def main(): """主优化流程""" logger.info("开始农业灌溉系统优化...") # 1. 加载和处理数据 daily_avg_moisture = load_soil_moisture_data() # 2. 生成喷头布局 sprinkler_df = generate_sprinkler_layout() sprinkler_df = calculate_sprinkler_coverage(sprinkler_df) logger.info(f"生成喷头数量: {len(sprinkler_df)}") # 验证喷头间距 spacing_ok = validate_sprinkler_spacing(sprinkler_df) if not spacing_ok: logger.warning("喷头间距不符合要求，可能需要调整布局!") # 3. 计算灌溉需求 irrigation_demand, sprinkler_df = calculate_irrigation_demand(daily_avg_moisture, sprinkler_df) logger.info(f"最大日灌溉需求: {sprinkler_df['max_demand'].sum()/sprinkler_df['area'].sum():.2f} L/m²") # 4. 确定最佳聚类数量 optimal_clusters = determine_optimal_clusters(sprinkler_df, max_clusters=8) # 5. 生成候选储水罐位置 tank_positions, cluster_labels = generate_candidate_tanks(sprinkler_df, optimal_clusters) logger.info(f"候选储水罐位置: {tank_positions}") # 6. 网络化优化 result, assignments, tank_demands, tank_mst, sprinkler_mst, root_sprinkler_idx, root_tank_idx = two_stage_optimization( sprinkler_df, irrigation_demand, tank_positions) if result.success: logger.info("优化成功!") # 7. 计算成本 total_cost, cost_breakdown, Q_river_opt, Q_tank_opt = calculate_total_cost( result, sprinkler_df, tank_positions, assignments, tank_mst) logger.info(f"总成本: {total_cost:.2f} 元") # 8. 可视化 visualize_network_system(sprinkler_df, tank_positions, assignments, tank_mst, sprinkler_mst, root_sprinkler_idx, root_tank_idx, Q_river_opt, Q_tank_opt) plot_cost_breakdown(cost_breakdown) # 9. 输出结果 output_results(sprinkler_df, tank_positions, result, assignments, tank_mst, Q_river_opt, Q_tank_opt, tank_demands, cost_breakdown) # 10. 最终验证报告 logger.info("\n最终系统验证报告:") logger.info(f"1. 喷头间距验证: {'通过' if spacing_ok else '失败'}") logger.info(f"2. 系统可靠性: {(total_cost - result.fun)/total_cost*100:.1f}%") else: logger.error("优化失败:", result.message) if __name__ == "__main__": main() ``` ### 完整代码说明： 1. **网络化连接设计**： - 喷头之间通过最小生成树(MST)形成网状连接 - 储水罐之间也通过MST形成最短连接 - 河流只连接最近的储水罐，储水罐连接喷头网络 2. **关键优化点**： - 喷头网络入口选择离河流最近的喷头 - 储水罐网络入口选择离河流最近的储水罐 - 每个储水罐连接到其负责区域中最近的喷头 - 管道成本基于实际网络路径和流量计算 3. **可视化增强**： - 清晰显示喷头网络和储水罐网络 - 突出显示网络入口点 - 标注管道长度和流量信息 - 区分不同类型的管道 4. **成本分析**： - 详细分解储水罐、管道等各项成本 - 提供成本饼图直观展示 - 输出完整的优化结果表格 这个实现满足了您的需求：喷头之间形成网状连接，储水罐之间形成最短连接，整个系统通过一根水管连接河流，最大程度降低了管道成本。