我自己对第二问进行了求解,现在我需要你把我修改一下需求: 确定储水罐位置和容量
初始假设储水罐容量为40,000L,初始成本为0,只有当容量超过40000L时,超出部分才按每升5元的成本计算。
根据农场布局和作物分布,确定储水罐的最佳位置
储水罐应位于能够覆盖尽可能多作物的位置
2. 基于储水罐位置确定喷头布局
确保喷头与储水罐的覆盖范围不重叠
喷头布局应确保全覆盖农场且满足最小间距要求
3. 管道网络优化
从河流引水到储水罐和喷头网络
优化管道布局以最小化成本: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")
# 作物面积比例 (公顷)
CROP_AREAS = {
'高粱': 0.5,
'玉米': 0.3,
'大豆': 0.2
}
# 计算垂直分布的高度边界
total_height = FARM_SIZE_Y
gaoliang_height = total_height * CROP_AREAS['高粱'] # 高粱区域高度
yumi_height = total_height * CROP_AREAS['玉米'] # 玉米区域高度
dadou_height = total_height * CROP_AREAS['大豆'] # 大豆区域高度
# 作物区域垂直分布 (高粱在下,玉米在中,大豆在上)
CROP_REGIONS = {
'高粱': {'x_min': 0, 'x_max': FARM_SIZE_X, 'y_min': 0, 'y_max': gaoliang_height},
'玉米': {'x_min': 0, 'x_max': FARM_SIZE_X, 'y_min': gaoliang_height, 'y_max': gaoliang_height + yumi_height},
'大豆': {'x_min': 0, 'x_max': FARM_SIZE_X, 'y_min': gaoliang_height + yumi_height, '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.0 # 问题2不需要应急储备水源
TANK_COVERAGE_RADIUS = 15 # 储水罐覆盖半径(问题2为15m)
MIN_DISTANCE_FROM_BORDER = 10 # 喷头离农场边线的最小距离
MIN_DISTANCE_BETWEEN_SPRINKLER_TANK = SPRINKLER_RADIUS + TANK_COVERAGE_RADIUS + 5 # 喷头和储水罐之间的最小距离
# ==================== 数据加载与处理 ====================
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):
"""计算每日灌溉需求"""
daily_demand_per_m2 = [calculate_daily_irrigation_demand(m) for m in daily_avg_moisture]
max_demand_per_m2 = max(daily_demand_per_m2)
avg_demand_per_m2 = np.mean(daily_demand_per_m2)
# 使用最大需求作为设计值(保守设计)
sprinkler_df['max_demand'] = sprinkler_df['area'] * max_demand_per_m2
# 记录日志
logger.info(f"最大日灌溉需求: {max_demand_per_m2:.2f} L/m²")
logger.info(f"平均日灌溉需求: {avg_demand_per_m2:.2f} L/m²")
plt.figure(figsize=(12, 6))
plt.plot(daily_avg_moisture.index, daily_demand_per_m2, 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_per_m2, sprinkler_df
# ==================== 喷头布局生成 ====================
def generate_sprinkler_layout(farm_size_x=FARM_SIZE_X, farm_size_y=FARM_SIZE_Y, radius=SPRINKLER_RADIUS):
"""生成六角形喷头布局,确保离边界至少10m"""
spacing = radius * 1.5
sprinklers = []
# 确保喷头离边界至少10m
min_x = MIN_DISTANCE_FROM_BORDER
max_x = farm_size_x - MIN_DISTANCE_FROM_BORDER
min_y = MIN_DISTANCE_FROM_BORDER
max_y = farm_size_y - MIN_DISTANCE_FROM_BORDER
rows = int((max_y - min_y) / (spacing * np.sqrt(3)/2)) + 2
cols = int((max_x - min_x) / spacing) + 2
for i in range(rows):
for j in range(cols):
x = min_x + j * spacing
y = min_y + i * spacing * np.sqrt(3)/2
if i % 2 == 1:
x += spacing / 2
if min_x <= x <= max_x and min_y <= y <= max_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'}
# 高粱区域(底部)
rect_gaoliang = Rectangle(
(CROP_REGIONS['高粱']['x_min'], CROP_REGIONS['高粱']['y_min']),
CROP_REGIONS['高粱']['x_max'] - CROP_REGIONS['高粱']['x_min'],
CROP_REGIONS['高粱']['y_max'] - CROP_REGIONS['高粱']['y_min'],
alpha=0.3, color=colors['高粱'], label=f'高粱 ({CROP_AREAS["高粱"]}公顷)'
)
plt.gca().add_patch(rect_gaoliang)
# 玉米区域(中部)
rect_yumi = Rectangle(
(CROP_REGIONS['玉米']['x_min'], CROP_REGIONS['玉米']['y_min']),
CROP_REGIONS['玉米']['x_max'] - CROP_REGIONS['玉米']['x_min'],
CROP_REGIONS['玉米']['y_max'] - CROP_REGIONS['玉米']['y_min'],
alpha=0.3, color=colors['玉米'], label=f'玉米 ({CROP_AREAS["玉米"]}公顷)'
)
plt.gca().add_patch(rect_yumi)
# 大豆区域(顶部)
rect_dadou = Rectangle(
(CROP_REGIONS['大豆']['x_min'], CROP_REGIONS['大豆']['y_min']),
CROP_REGIONS['大豆']['x_max'] - CROP_REGIONS['大豆']['x_min'],
CROP_REGIONS['大豆']['y_max'] - CROP_REGIONS['大豆']['y_min'],
alpha=0.3, color=colors['大豆'], label=f'大豆 ({CROP_AREAS["大豆"]}公顷)'
)
plt.gca().add_patch(rect_dadou)
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_network_flows(sprinkler_df, root_idx):
"""计算喷头网络中每条边的流量(使用BFS遍历)"""
# 构建完全连接的网络图
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 in range(n):
for j in range(i+1, n):
x1, y1 = sprinkler_df.iloc[i][['x', 'y']]
x2, y2 = sprinkler_df.iloc[j][['x', 'y']]
L = np.sqrt((x1 - x2)**2 + (y1 - y2)**2)
G.add_edge(i, j, length=L)
# 计算最小生成树(形成实际连接)
mst_edges = list(nx.minimum_spanning_edges(G, weight='length', data=True))
# 计算每条边的流量
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:
# 检查这条边是否在MST中
edge_in_mst = any((min(node, neighbor) == min(i, j) and
max(node, neighbor) == max(i, j))
for i, j, _ in mst_edges)
if edge_in_mst:
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, mst_edges
def determine_optimal_clusters(sprinkler_df, max_clusters=10):
"""确定最佳聚类数量"""
coordinates = sprinkler_df[['x', 'y']].values
scaler = StandardScaler()
scaled_features = scaler.fit_transform(coordinates)
sse = [] # 这里定义了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_) # 这里应该是sse而不是se
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):
"""生成候选储水罐位置,确保不与喷头位置重合,并且协同覆盖整个农场"""
# 计算农场的四个角落和中心点作为初始候选位置
candidate_positions = [
[MIN_DISTANCE_FROM_BORDER, MIN_DISTANCE_FROM_BORDER], # 左下角
[FARM_SIZE_X - MIN_DISTANCE_FROM_BORDER, MIN_DISTANCE_FROM_BORDER], # 右下角
[FARM_SIZE_X - MIN_DISTANCE_FROM_BORDER, FARM_SIZE_Y - MIN_DISTANCE_FROM_BORDER], # 右上角
[MIN_DISTANCE_FROM_BORDER, FARM_SIZE_Y - MIN_DISTANCE_FROM_BORDER], # 左上角
[FARM_SIZE_X / 2, FARM_SIZE_Y / 2] # 中心点
]
# 如果需要的储水罐数量多于预设位置,使用KMeans生成额外位置
if num_tanks > len(candidate_positions):
coordinates = sprinkler_df[['x', 'y']].values
scaler = StandardScaler()
scaled_features = scaler.fit_transform(coordinates)
kmeans = KMeans(n_clusters=num_tanks - len(candidate_positions), random_state=42, n_init=10)
kmeans.fit(scaled_features)
additional_centers = scaler.inverse_transform(kmeans.cluster_centers_)
# 将额外位置添加到候选位置
for center in additional_centers:
candidate_positions.append([center[0], center[1]])
# 只保留前num_tanks个位置
candidate_positions = candidate_positions[:num_tanks]
# 调整储水罐位置,确保不与喷头位置重合
for i in range(len(candidate_positions)):
tank_pos = candidate_positions[i]
min_dist = float('inf')
closest_sprinkler_idx = -1
# 找到最近的喷头
for j, sprinkler in sprinkler_df.iterrows():
dist = np.sqrt((tank_pos[0] - sprinkler['x'])**2 +
(tank_pos[1] - sprinkler['y'])**2)
if dist < min_dist:
min_dist = dist
closest_sprinkler_idx = j
# 如果距离太近,调整储水罐位置
if min_dist < MIN_DISTANCE_BETWEEN_SPRINKLER_TANK:
closest_sprinkler = sprinkler_df.iloc[closest_sprinkler_idx]
dx = tank_pos[0] - closest_sprinkler['x']
dy = tank_pos[1] - closest_sprinkler['y']
angle = np.arctan2(dy, dx)
# 将储水罐移动到最小距离之外
new_x = closest_sprinkler['x'] + np.cos(angle) * MIN_DISTANCE_BETWEEN_SPRINKLER_TANK
new_y = closest_sprinkler['y'] + np.sin(angle) * MIN_DISTANCE_BETWEEN_SPRINKLER_TANK
# 确保新位置在农场范围内
new_x = max(MIN_DISTANCE_FROM_BORDER, min(FARM_SIZE_X - MIN_DISTANCE_FROM_BORDER, new_x))
new_y = max(MIN_DISTANCE_FROM_BORDER, min(FARM_SIZE_Y - MIN_DISTANCE_FROM_BORDER, new_y))
candidate_positions[i] = [new_x, new_y]
# 确保储水罐之间也有足够距离
for i in range(len(candidate_positions)):
for j in range(i+1, len(candidate_positions)):
dist = np.sqrt((candidate_positions[i][0] - candidate_positions[j][0])**2 +
(candidate_positions[i][1] - candidate_positions[j][1])**2)
if dist < MIN_DISTANCE_BETWEEN_SPRINKLER_TANK:
# 调整位置使它们分开
dx = candidate_positions[j][0] - candidate_positions[i][0]
dy = candidate_positions[j][1] - candidate_positions[i][1]
angle = np.arctan2(dy, dx)
# 将第二个储水罐移动到最小距离之外
new_x = candidate_positions[i][0] + np.cos(angle) * MIN_DISTANCE_BETWEEN_SPRINKLER_TANK
new_y = candidate_positions[i][1] + np.sin(angle) * MIN_DISTANCE_BETWEEN_SPRINKLER_TANK
# 确保新位置在农场范围内
new_x = max(MIN_DISTANCE_FROM_BORDER, min(FARM_SIZE_X - MIN_DISTANCE_FROM_BORDER, new_x))
new_y = max(MIN_DISTANCE_FROM_BORDER, min(FARM_SIZE_Y - MIN_DISTANCE_FROM_BORDER, new_y))
candidate_positions[j] = [new_x, new_y]
# 创建标签(所有喷头都分配到最近的储水罐)
labels = []
for i, sprinkler in sprinkler_df.iterrows():
min_dist = float('inf')
closest_tank_idx = -1
for j, tank_pos in enumerate(candidate_positions):
dist = np.sqrt((sprinkler['x'] - tank_pos[0])**2 +
(sprinkler['y'] - tank_pos[1])**2)
if dist < min_dist:
min_dist = dist
closest_tank_idx = j
labels.append(closest_tank_idx)
return candidate_positions, np.array(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_G = nx.Graph()
n_sprinklers = len(sprinkler_df)
for i in range(n_sprinklers):
sprinkler_G.add_node(i, demand=sprinkler_df.iloc[i]['max_demand'])
# 创建所有喷头之间的连接
for i in range(n_sprinklers):
for j in range(i+1, n_sprinklers):
x1, y1 = sprinkler_df.iloc[i][['x', 'y']]
x2, y2 = sprinkler_df.iloc[j][['x', 'y']]
L = np.sqrt((x1 - x2)**2 + (y1 - y2)**2)
sprinkler_G.add_edge(i, j, length=L)
# 计算喷头最小生成树(形成实际连接)
sprinkler_mst = list(nx.minimum_spanning_edges(sprinkler_G, weight='length', data=True))
# 找到离河流最近的喷头作为喷头网络入口
distances_to_river = []
for i in range(n_sprinklers):
x, y = sprinkler_df.iloc[i][['x', 'y']]
dist = np.sqrt((x - RIVER_POINT[0])**2 + (y - RIVER_POINT[1])**2)
distances_to_river.append(dist)
root_sprinkler_idx = np.argmin(distances_to_river)
# 计算喷头网络中各边的流量
edge_flows, _ = calculate_network_flows(sprinkler_df, root_sprinkler_idx)
# 阶段2: 构建储水罐网络(完全连接)
num_tanks = len(tank_positions)
# 构建储水罐完全连接网络图
tank_G = nx.Graph()
for j in range(num_tanks):
tank_G.add_node(j)
for i in range(num_tanks):
for j in range(i+1, num_tanks):
x1, y1 = tank_positions[i]
x2, y2 = tank_positions[j]
L = np.sqrt((x1 - x2)**2 + (y1 - y2)**2)
tank_G.add_edge(i, j, length=L)
# 计算储水罐最小生成树
tank_mst = list(nx.minimum_spanning_edges(tank_G, weight='length', data=True))
# 找到离河流最近的储水罐作为储水罐网络入口
distances_to_river_tank = [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_tank)
# 分配喷头到储水罐(基于最近距离)
distances = calculate_distance_matrix(
sprinkler_df[['x', 'y']].values,
np.array(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)
# 优化目标函数
def objective(vars):
# 解析变量
V = vars[:num_tanks] # 储水罐容量(L)
Q_total = vars[num_tanks] # 总管流量(L/day)
Q_sprinkler = vars[num_tanks+1] # 分配到喷头网络的流量(L/day)
Q_tank_network = vars[num_tanks+2:2*num_tanks+2] # 分配到储水罐网络的流量(L/day)
# 总管成本(从河流到分流点)
main_pipe_cost = 0
L_main = distances_to_river[root_sprinkler_idx] # 使用到最近喷头的距离作为总管长度
Q_m3 = Q_total * L_TO_M3
main_pipe_cost += PIPE_LENGTH_COEFF * (L_main ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP)
# 喷头网络管道成本
sprinkler_pipe_cost = 0
# 分流点到喷头网络入口的管道
L_sprinkler = 0 # 分流点就是喷头网络入口,所以长度为0
Q_m3 = Q_sprinkler * L_TO_M3
sprinkler_pipe_cost += PIPE_LENGTH_COEFF * (L_sprinkler ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP)
# 喷头之间的管道成本(只计算MST中的边)
for (i, j), flow in edge_flows.items():
x1, y1 = sprinkler_df.iloc[i][['x', 'y']]
x2, y2 = sprinkler_df.iloc[j][['x', 'y']]
L = np.sqrt((x1 - x2)**2 + (y1 - y2)**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)
# 储水罐成本
tank_cost = TANK_COST_PER_LITER * np.sum(V)
# 分流点到储水罐网络的管道成本
tank_pipe_cost = 0
# 分流点到储水罐网络入口的管道
L_tank = np.sqrt((sprinkler_df.iloc[root_sprinkler_idx]['x'] - tank_positions[root_tank_idx][0])**2 +
(sprinkler_df.iloc[root_sprinkler_idx]['y'] - tank_positions[root_tank_idx][1])**2)
Q_m3 = np.sum(Q_tank_network) * L_TO_M3
tank_pipe_cost += PIPE_LENGTH_COEFF * (L_tank ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP)
# 储水罐之间的管道成本(只计算MST中的边)
for i, j, data in tank_mst:
L = data['length']
# 计算两个储水罐之间的流量(取最大值)
Q_avg = max(Q_tank_network[i], Q_tank_network[j])
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)
# 惩罚项
penalty = 0
# 总管流量约束
total_demand = np.sum(sprinkler_df['max_demand'])
if Q_total < total_demand:
penalty += 1000 * (total_demand - Q_total)
# 喷头网络流量约束
if Q_sprinkler < total_demand * DAILY_WATER_SOURCE_RATIO:
penalty += 1000 * (total_demand * DAILY_WATER_SOURCE_RATIO - Q_sprinkler)
# 储水罐约束
for j in range(num_tanks):
# 储水罐容量只需要满足部分需求,因为河流提供80%
required_capacity = tank_demands[j] * (1 - DAILY_WATER_SOURCE_RATIO) * 1.2 # 增加20%的安全余量
if V[j] < required_capacity:
penalty += 1000 * (required_capacity - V[j])
if Q_tank_network[j] < tank_demands[j] * (1 - DAILY_WATER_SOURCE_RATIO):
penalty += 1000 * (tank_demands[j] * (1 - DAILY_WATER_SOURCE_RATIO) - Q_tank_network[j])
return tank_cost + main_pipe_cost + sprinkler_pipe_cost + tank_pipe_cost + penalty
# 约束条件
constraints = []
# 初始值
total_demand = np.sum(sprinkler_df['max_demand'])
initial_V = [tank_demands[j] * (1 - DAILY_WATER_SOURCE_RATIO) * 1.2 for j in range(num_tanks)] # 增加20%安全余量
initial_Q_total = total_demand
initial_Q_sprinkler = total_demand * DAILY_WATER_SOURCE_RATIO
initial_Q_tank_network = [tank_demands[j] * (1 - DAILY_WATER_SOURCE_RATIO) for j in range(num_tanks)]
x0 = initial_V + [initial_Q_total, initial_Q_sprinkler] + initial_Q_tank_network
# 边界
bounds = Bounds([0] * (2 * num_tanks + 2), [np.inf] * (2 * num_tanks + 2))
# 优化
logger.info("开始优化...")
result = minimize(
objective,
x0,
bounds=bounds,
constraints=constraints,
method='SLSQP',
options={'disp': True, 'ftol': 1e-6, 'maxiter': 100}
)
# 提取优化结果
V_opt = result.x[:num_tanks] # 储水罐容量(L)
Q_total_opt = result.x[num_tanks] # 总管流量(L/day)
Q_sprinkler_opt = result.x[num_tanks+1] # 分配到喷头网络的流量(L/day)
Q_tank_network_opt = result.x[num_tanks+2:2*num_tanks+2] # 分配到储水罐网络的流量(L/day)
return result, assignments, tank_demands, tank_mst, sprinkler_mst, root_sprinkler_idx, root_tank_idx, V_opt, Q_total_opt, Q_sprinkler_opt, Q_tank_network_opt
# ==================== 成本计算与可视化 ====================
def calculate_total_cost(result, sprinkler_df, tank_positions, assignments, tank_mst, root_sprinkler_idx, root_tank_idx):
num_tanks = len(tank_positions)
V_opt = result.x[:num_tanks]
Q_total_opt = result.x[num_tanks]
Q_sprinkler_opt = result.x[num_tanks+1]
Q_tank_network_opt = result.x[num_tanks+2:2*num_tanks+2]
# 总管成本(从河流到分流点)
main_pipe_cost = 0
distances_to_river = []
for i in range(len(sprinkler_df)):
x, y = sprinkler_df.iloc[i][['x', 'y']]
dist = np.sqrt((x - RIVER_POINT[0])**2 + (y - RIVER_POINT[1])**2)
distances_to_river.append(dist)
root_sprinkler_idx = np.argmin(distances_to_river)
L_main = distances_to_river[root_sprinkler_idx]
Q_m3 = Q_total_opt * L_TO_M3
main_pipe_cost += PIPE_LENGTH_COEFF * (L_main ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP)
# 喷头网络管道成本
sprinkler_pipe_cost = 0
# 分流点到喷头网络入口的管道(长度为0)
L_sprinkler = 0
Q_m3 = Q_sprinkler_opt * L_TO_M3
sprinkler_pipe_cost += PIPE_LENGTH_COEFF * (L_sprinkler ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP)
# 喷头之间的管道成本
edge_flows, _ = calculate_network_flows(sprinkler_df, root_sprinkler_idx)
for (i, j), flow in edge_flows.items():
x1, y1 = sprinkler_df.iloc[i][['x', 'y']]
x2, y2 = sprinkler_df.iloc[j][['x', 'y']]
L = np.sqrt((x1 - x2)**2 + (y1 - y2)**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)
# 储水罐成本
tank_cost = TANK_COST_PER_LITER * np.sum(V_opt)
# 分流点到储水罐网络的管道成本
tank_pipe_cost = 0
# 分流点到储水罐网络入口的管道
L_tank = np.sqrt((sprinkler_df.iloc[root_sprinkler_idx]['x'] - tank_positions[root_tank_idx][0])**2 +
(sprinkler_df.iloc[root_sprinkler_idx]['y'] - tank_positions[root_tank_idx][1])**2)
Q_m3 = np.sum(Q_tank_network_opt) * L_TO_M3
tank_pipe_cost += PIPE_LENGTH_COEFF * (L_tank ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP)
# 储水罐之间的管道成本
for i, j, data in tank_mst:
L = data['length']
Q_avg = max(Q_tank_network_opt[i], Q_tank_network_opt[j])
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)
total_cost = tank_cost + main_pipe_cost + sprinkler_pipe_cost + tank_pipe_cost
cost_breakdown = {
'tank_cost': tank_cost,
'main_pipe_cost': main_pipe_cost,
'sprinkler_pipe_cost': sprinkler_pipe_cost,
'tank_pipe_cost': tank_pipe_cost
}
return total_cost, cost_breakdown, Q_total_opt, Q_sprinkler_opt, Q_tank_network_opt
def visualize_network_system(sprinkler_df, tank_positions, assignments, tank_mst, sprinkler_mst,
root_sprinkler_idx, root_tank_idx, Q_total_opt, Q_sprinkler_opt, Q_tank_network_opt):
"""可视化网络连接的灌溉系统"""
plt.figure(figsize=(16, 14))
ax = plt.gca()
# 设置坐标轴等比例,确保圆形显示为圆形
ax.set_aspect('equal')
# 绘制农场边界和河流
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'}
# 高粱区域(底部0-50m)
rect_gaoliang = Rectangle(
(0, 0), FARM_SIZE_X, 50,
alpha=0.2, color=colors['高粱'], label='高粱 (0.5公顷)'
)
ax.add_patch(rect_gaoliang)
# 玉米区域(中部50-80m)
rect_yumi = Rectangle(
(0, 50), FARM_SIZE_X, 30,
alpha=0.2, color=colors['玉米'], label='玉米 (0.3公顷)'
)
ax.add_patch(rect_yumi)
# 大豆区域(顶部80-100m)
rect_dadou = Rectangle(
(0, 80), FARM_SIZE_X, 20,
alpha=0.2, color=colors['大豆'], label='大豆 (0.2公顷)'
)
ax.add_patch(rect_dadou)
# 绘制喷头
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)
# 绘制喷头覆盖范围(确保是圆形)
circle = Circle((sprinkler['x'], sprinkler['y']), SPRINKLER_RADIUS,
color='blue', alpha=0.1, fill=True)
ax.add_patch(circle)
# 绘制喷头之间的连接(MST边)
for i, j, data in sprinkler_mst:
x1, y1 = sprinkler_df.iloc[i][['x', 'y']]
x2, y2 = sprinkler_df.iloc[j][['x', 'y']]
ax.plot([x1, x2], [y1, y2], '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)
# 绘制储水罐覆盖范围(确保是圆形)
circle = Circle((tank[0], tank[1]), TANK_COVERAGE_RADIUS,
color='purple', alpha=0.15, fill=True)
ax.add_patch(circle)
ax.text(tank[0], tank[1], f'T{j+1}', fontsize=10,
ha='center', va='center', fontweight='bold')
# 绘制储水罐之间的连接(MST边)
for i, j, data in tank_mst:
x1, y1 = tank_positions[i]
x2, y2 = tank_positions[j]
ax.plot([x1, x2], [y1, y2], 'g-', linewidth=2, label='储水罐间管道' if i == 0 and j == 1 else "")
# 标注管道长度
mid_x = (x1 + x2) / 2
mid_y = (y1 + y2) / 2
ax.text(mid_x, mid_y, f'{data["length"]:.1f}m', fontsize=8,
bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.7))
# 绘制总管(河流到分流点)
connection_point = sprinkler_df.iloc[root_sprinkler_idx][['x', 'y']].values
ax.plot([RIVER_POINT[0], connection_point[0]],
[RIVER_POINT[1], connection_point[1]],
'r-', linewidth=3, label='总管')
# 标注总管信息
mid_x = (RIVER_POINT[0] + connection_point[0]) / 2
mid_y = (RIVER_POINT[1] + connection_point[1]) / 2
length = np.sqrt((RIVER_POINT[0]-connection_point[0])**2 +
(RIVER_POINT[1]-connection_point[1])**2)
Q_m3 = Q_total_opt * L_TO_M3
ax.text(mid_x, mid_y,
f'{length:.1f}m\n{Q_total_opt:.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))
# 绘制分流点到储水罐网络的管道
tank_pos = tank_positions[root_tank_idx]
ax.plot([connection_point[0], tank_pos[0]],
[connection_point[1], tank_pos[1]],
'm--', linewidth=2, label='分流点到储水罐网络')
# 标注分流点到储水罐管道信息
mid_x = (connection_point[0] + tank_pos[0]) / 2
mid_y = (connection_point[1] + tank_pos[1]) / 2
L = np.sqrt((connection_point[0]-tank_pos[0])**2 +
(connection_point[1]-tank_pos[1])**2)
Q_total_tank = np.sum(Q_tank_network_opt)
Q_m3 = Q_total_tank * L_TO_M3
ax.text(mid_x, mid_y,
f'{L:.1f}m\n{Q_total_tank:.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))
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['main_pipe_cost'],
cost_breakdown['sprinkler_pipe_cost'],
cost_breakdown['tank_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, V_opt, assignments, tank_mst,
Q_total_opt, Q_sprinkler_opt, Q_tank_network_opt, tank_demands, cost_breakdown):
"""输出优化结果表格"""
num_tanks = len(tank_positions)
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_tank_network_opt],
'分配到储水罐流量(m³/天)': [f"{q * L_TO_M3:.2f}" for q in Q_tank_network_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, data in tank_mst:
tank_pipe_info.append(f'T{i+1}-T{j+1}: {data["length"]:.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 = Q_total_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"分配到喷头网络流量: {Q_sprinkler_opt:.0f} L/天 ({Q_sprinkler_opt * L_TO_M3:.2f} m³/天)")
logger.info(f"分配到储水罐网络流量: {np.sum(Q_tank_network_opt):.0f} L/天 ({np.sum(Q_tank_network_opt) * 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['main_pipe_cost']:.2f} 元")
logger.info(f"喷头网络管道成本: {cost_breakdown['sprinkler_pipe_cost']:.2f} 元")
logger.info(f"储水罐管道成本: {cost_breakdown['tank_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)
# 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, V_opt, Q_total_opt, Q_sprinkler_opt, Q_tank_network_opt = two_stage_optimization(
sprinkler_df, irrigation_demand, tank_positions)
if result.success:
logger.info("优化成功!")
# 7. 计算成本
total_cost, cost_breakdown, Q_total_opt, Q_sprinkler_opt, Q_tank_network_opt = calculate_total_cost(
result, sprinkler_df, tank_positions, assignments, tank_mst, root_sprinkler_idx, root_tank_idx)
# 8. 可视化
visualize_network_system(sprinkler_df, tank_positions, assignments, tank_mst, sprinkler_mst,
root_sprinkler_idx, root_tank_idx, Q_total_opt, Q_sprinkler_opt, Q_tank_network_opt)
plot_cost_breakdown(cost_breakdown)
# 9. 输出结果
output_results(sprinkler_df, tank_positions, V_opt, assignments, tank_mst,
Q_total_opt, Q_sprinkler_opt, Q_tank_network_opt, tank_demands, cost_breakdown)
# 10. 最终验证报告
logger.info("\n最终系统验证报告:")
logger.info(f"1. 喷头间距验证: {'通过' if spacing_ok else '失败'}")
logger.info(f"2. 系统可靠性: {total_cost:.2f} 元")
else:
logger.error("优化失败:", result.message)
if __name__ == "__main__":
main()
最新发布