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()
代码中并没有实现喷头需要全部互相连接，然后挑出一根连接河流，水罐也一样
扫一扫
388
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
