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
# 设置日志记录
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 # 喷头半径
TANK_COVERAGE_RADIUS = SPRINKLER_RADIUS # 储水罐覆盖半径与喷头相同
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 # 管道流量指数
BASE_TANK_CAPACITY = 40000 # 基础免费容量(L)
EXTRA_TANK_COST = 5 # 超出基础容量的单位成本(元/L)
# 单位转换系数
L_TO_M3 = 0.001 # 1L = 0.001m³
# 系统参数
DAILY_WATER_SOURCE_RATIO = 0.8 # 日常水源中河水的比例
EMERGENCY_WATER_SOURCE_RATIO = 0.0 # 不需要应急储备水源
MIN_DISTANCE_FROM_BORDER = 10 # 喷头离农场边线的最小距离
MIN_DISTANCE_BETWEEN_SPRINKLER_TANK = 0 # 喷头和储水罐位置相同
# ==================== 数据加载与处理 ====================
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.rcParams['font.sans-serif'] = ['Microsoft YaHei'] # 微软雅黑支持更多符号
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
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并保持喷头数量为23个"""
spacing = radius * 1.3 # 减小间距以保持喷头数量
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 = 5 # 固定行数以控制喷头数量
cols = 5 # 固定列数以控制喷头数量
# 生成喷头位置
for i in range(rows):
for j in range(cols):
# 计算位置,确保在缩小的区域内分布均匀
x = min_x + j * (max_x - min_x) / (cols - 1) if cols > 1 else min_x
y = min_y + i * (max_y - min_y) / (rows - 1) if rows > 1 else min_y
# 六角形布局偏移
if i % 2 == 1:
x += (max_x - min_x) / (cols - 1) / 2 if cols > 1 else 0
# 确保位置在农场内
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,
'is_tank': False # 初始都不是储水罐
})
# 如果喷头数量不足23,在中间区域补充
while len(sprinklers) < 23:
# 在中心区域添加额外喷头
center_x = (min_x + max_x) / 2
center_y = (min_y + max_y) / 2
# 围绕中心生成额外喷头
angle = (len(sprinklers) - rows*cols) * (2 * np.pi / 12)
dist = (max_x - min_x) / 6
x = center_x + dist * np.cos(angle)
y = center_y + dist * np.sin(angle)
# 确保在边界内
x = max(min_x, min(x, max_x))
y = max(min_y, min(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,
'is_tank': False
})
# 如果喷头数量超过23,移除边缘的喷头
if len(sprinklers) > 23:
# 计算每个喷头到中心的距离
center_x = (min_x + max_x) / 2
center_y = (min_y + max_y) / 2
distances = []
for idx, sprinkler in enumerate(sprinklers):
dist = np.sqrt((sprinkler['x'] - center_x)**2 + (sprinkler['y'] - center_y)** 2)
distances.append((dist, idx))
# 排序并保留距离中心最近的23个喷头
distances.sort(reverse=True) # 从远到近排序
to_remove = [idx for (dist, idx) in distances[:len(sprinklers)-23]]
sprinklers = [s for idx, s in enumerate(sprinklers) if idx not in to_remove]
# 重新编号
for i, s in enumerate(sprinklers):
s['id'] = i
return pd.DataFrame(sprinklers)
def adjust_specific_rows_spacing(sprinkler_df, add_distance=10):
"""
调整第二行和倒数第二行喷头间距,总增加距离为add_distance
以中间的喷头为参照物,保持其位置不变,左右两侧喷头分别向外移动
对第二行和倒数第二行的右侧第一喷头再向右移动6米
对第二行和倒数第二行的右侧第二个喷头再向右移动2米(原1米基础上再加1米)
对第二行和倒数第二行的左侧第一喷头再向左移动2米
"""
# 按y坐标分组,识别各行喷头
rows_dict = {}
for _, row in sprinkler_df.iterrows():
y_coord = row['y']
if y_coord not in rows_dict:
rows_dict[y_coord] = []
rows_dict[y_coord].append(row.to_dict()) # 转字典方便修改
# 按y坐标排序行,方便取第二行、倒数第二行
sorted_row_ys = sorted(rows_dict.keys())
num_rows = len(sorted_row_ys)
# 调整第二行(索引1)和倒数第二行(索引-2)间距
if num_rows >= 3: # 至少需要3行才能有第二行和倒数第二行
# 处理第二行
second_row_y = sorted_row_ys[1]
second_row_sprinklers = rows_dict[second_row_y]
second_row_sprinklers.sort(key=lambda s: s['x']) # 按x排序
if len(second_row_sprinklers) > 1:
# 计算每个间隔需要增加的距离
spacing_increase = add_distance / (len(second_row_sprinklers) - 1)
# 确定中间喷头的索引(参照物)
mid_idx = len(second_row_sprinklers) // 2
# 左侧喷头向左移动(以中间喷头为参照)
for j in range(mid_idx):
# 左侧偏移量:从中间向左侧递增
offset = (mid_idx - j) * spacing_increase
second_row_sprinklers[j]['x'] -= offset
# 右侧喷头向右移动(以中间喷头为参照)
for j in range(mid_idx + 1, len(second_row_sprinklers)):
# 右侧偏移量:从中间向右侧递增
offset = (j - mid_idx) * spacing_increase
second_row_sprinklers[j]['x'] += offset
# 对右侧第一喷头额外再向右移动6米
if len(second_row_sprinklers) > 0:
rightmost_idx = len(second_row_sprinklers) - 1
second_row_sprinklers[rightmost_idx]['x'] += 6
# 对右侧第二个喷头额外再向右移动2米(原1米基础上增加1米)
if len(second_row_sprinklers) > 1:
second_right_idx = len(second_row_sprinklers) - 2
second_row_sprinklers[second_right_idx]['x'] += 2
# 对左侧第一喷头额外再向左移动2米
if len(second_row_sprinklers) > 0:
leftmost_idx = 0
second_row_sprinklers[leftmost_idx]['x'] -= 2
# 处理倒数第二行
second_last_row_y = sorted_row_ys[-2] # 倒数第二行
second_last_row_sprinklers = rows_dict[second_last_row_y]
second_last_row_sprinklers.sort(key=lambda s: s['x']) # 按x排序
if len(second_last_row_sprinklers) > 1:
spacing_increase = add_distance / (len(second_last_row_sprinklers) - 1)
mid_idx = len(second_last_row_sprinklers) // 2 # 中间喷头索引
# 左侧喷头向左移动
for j in range(mid_idx):
offset = (mid_idx - j) * spacing_increase
second_last_row_sprinklers[j]['x'] -= offset
# 右侧喷头向右移动
for j in range(mid_idx + 1, len(second_last_row_sprinklers)):
offset = (j - mid_idx) * spacing_increase
second_last_row_sprinklers[j]['x'] += offset
# 对右侧第一喷头额外再向右移动6米
if len(second_last_row_sprinklers) > 0:
rightmost_idx = len(second_last_row_sprinklers) - 1
second_last_row_sprinklers[rightmost_idx]['x'] += 6
# 对右侧第二个喷头额外再向右移动2米(原1米基础上增加1米)
if len(second_last_row_sprinklers) > 1:
second_right_idx = len(second_last_row_sprinklers) - 2
second_last_row_sprinklers[second_right_idx]['x'] += 2
# 对左侧第一喷头额外再向左移动2米
if len(second_last_row_sprinklers) > 0:
leftmost_idx = 0
second_last_row_sprinklers[leftmost_idx]['x'] -= 2
# 把字典转回DataFrame
adjusted_sprinklers = []
for y in rows_dict:
adjusted_sprinklers.extend(rows_dict[y])
adjusted_df = pd.DataFrame(adjusted_sprinklers)
# 重新按id排序
adjusted_df.sort_values(by='id', inplace=True)
adjusted_df.reset_index(drop=True, inplace=True)
return adjusted_df
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 select_tank_positions(sprinkler_df, num_tanks=3):
"""从喷头位置中选择储水罐位置"""
# 使用K-means聚类确定最优位置
coords = sprinkler_df[['x','y']].values
kmeans = KMeans(n_clusters=num_tanks, n_init=10, random_state=42).fit(coords)
cluster_centers = kmeans.cluster_centers_
# 找到距离聚类中心最近的喷头作为储水罐位置
tank_indices = []
tank_positions = []
for center in cluster_centers:
min_dist = float('inf')
best_idx = -1
for idx, row in sprinkler_df.iterrows():
dist = np.sqrt((center[0]-row['x'])**2 + (center[1]-row['y'])** 2)
if dist < min_dist:
min_dist = dist
best_idx = idx
if best_idx not in tank_indices:
tank_indices.append(best_idx)
tank_positions.append([sprinkler_df.loc[best_idx, 'x'], sprinkler_df.loc[best_idx, 'y']])
sprinkler_df.loc[best_idx, 'is_tank'] = True
return tank_positions, tank_indices, sprinkler_df
# ==================== 网络连接优化 ====================
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 = [] # 误差平方和
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 network_optimization(sprinkler_df, irrigation_demand, tank_positions, tank_indices):
"""
网络优化:储水罐与喷头形成统一网络,总管直接连接到最近的喷头
"""
# 确定总管连接的喷头(离河流最近的喷头)
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)
logger.info(f"总管连接到最近的喷头: 喷头 #{root_sprinkler_idx}")
# 构建统一网络(包括所有喷头和储水罐)
all_points = sprinkler_df[['x', 'y']].values
# 构建完全连接网络图
G = nx.Graph()
n_points = len(sprinkler_df)
for i in range(n_points):
G.add_node(i, demand=sprinkler_df.iloc[i]['max_demand'])
# 创建所有点之间的连接
for i in range(n_points):
for j in range(i+1, n_points):
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 = list(nx.minimum_spanning_edges(G, weight='length', data=True))
# 计算网络中各边的流量
edge_flows, _ = calculate_network_flows(sprinkler_df, root_sprinkler_idx)
# 分配喷头到储水罐(基于最近距离)
num_tanks = len(tank_positions)
distances = np.zeros((len(sprinkler_df), num_tanks))
for i, row in sprinkler_df.iterrows():
for j, tank in enumerate(tank_positions):
distances[i, j] = np.sqrt((row['x']-tank[0])**2 + (row['y']-tank[1])** 2)
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)
# 总管成本(从河流到连接的喷头)
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)
# 网络管道成本(所有喷头和储水罐之间的连接)
network_pipe_cost = 0
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
network_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP)
# 储水罐成本(仅计算超出基础容量的部分)
tank_cost = 0
for v in V:
if v > BASE_TANK_CAPACITY:
tank_cost += (v - BASE_TANK_CAPACITY) * EXTRA_TANK_COST
# 惩罚项
penalty = 0
# 总管流量约束
total_demand = np.sum(sprinkler_df['max_demand'])
if Q_total < total_demand * DAILY_WATER_SOURCE_RATIO:
penalty += 1000 * (total_demand * DAILY_WATER_SOURCE_RATIO - Q_total)
# 储水罐约束:总容量需满足非日常水源部分
total_tank_capacity = sum(V)
required_capacity = total_demand * (1 - DAILY_WATER_SOURCE_RATIO) * 1.2 # 增加20%的安全余量
if total_tank_capacity < required_capacity:
penalty += 1000 * (required_capacity - total_tank_capacity)
return tank_cost + main_pipe_cost + network_pipe_cost + penalty
# 初始值
total_demand = np.sum(sprinkler_df['max_demand'])
# 按比例分配储水罐容量
tank_ratios = [d / sum(tank_demands) if sum(tank_demands) > 0 else 1/num_tanks for d in tank_demands]
required_capacity = total_demand * (1 - DAILY_WATER_SOURCE_RATIO) * 1.2
initial_V = [ratio * required_capacity for ratio in tank_ratios]
initial_Q_total = total_demand * DAILY_WATER_SOURCE_RATIO
x0 = initial_V + [initial_Q_total]
# 边界
bounds = Bounds([0] * (num_tanks + 1), [np.inf] * (num_tanks + 1))
# 优化
logger.info("开始优化...")
result = minimize(
objective,
x0,
bounds=bounds,
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_tank_network_opt = [tank_demands[j] * (1 - DAILY_WATER_SOURCE_RATIO) for j in range(num_tanks)]
return result, assignments, tank_demands, None, mst, root_sprinkler_idx, V_opt, Q_total_opt, Q_tank_network_opt
# ==================== 成本计算与可视化 ====================
def calculate_total_cost(result, sprinkler_df, tank_positions, assignments, tank_mst, root_sprinkler_idx):
num_tanks = len(tank_positions)
V_opt = result.x[:num_tanks]
Q_total_opt = result.x[num_tanks]
# 计算储水罐需求
tank_demands = []
sprinkler_max_demands = sprinkler_df['max_demand'].values
distances = np.zeros((len(sprinkler_df), num_tanks))
for i, row in sprinkler_df.iterrows():
for j, tank in enumerate(tank_positions):
distances[i, j] = np.sqrt((row['x']-tank[0])**2 + (row['y']-tank[1])** 2)
assignments = np.argmin(distances, axis=1)
for j in range(num_tanks):
demand = np.sum(sprinkler_max_demands[assignments == j])
tank_demands.append(demand)
# 计算每个储水罐的网络流量(与network_optimization保持一致)
Q_tank_network_opt = [tank_demands[j] * (1 - DAILY_WATER_SOURCE_RATIO) for j in range(num_tanks)]
# 总管成本(从河流到连接的喷头)
main_pipe_cost = 0
distances_to_river = [np.sqrt((x - RIVER_POINT[0])**2 + (y - RIVER_POINT[1])** 2)
for _, (x, y) in sprinkler_df[['x', 'y']].iterrows()]
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)
# 网络管道成本
network_pipe_cost = 0
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
network_pipe_cost += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_m3 ** PIPE_FLOW_EXP)
# 储水罐成本(仅计算超出基础容量的部分)
tank_cost = 0
for v in V_opt:
if v > BASE_TANK_CAPACITY:
tank_cost += (v - BASE_TANK_CAPACITY) * EXTRA_TANK_COST
total_cost = tank_cost + main_pipe_cost + network_pipe_cost
cost_breakdown = {
'tank_cost': tank_cost,
'main_pipe_cost': main_pipe_cost,
'network_pipe_cost': network_pipe_cost
}
return total_cost, cost_breakdown, Q_total_opt, Q_tank_network_opt
def visualize_network_system(sprinkler_df, tank_positions, tank_indices, assignments, tank_mst, network_mst,
root_sprinkler_idx, Q_total_opt, Q_tank_network_opt):
"""可视化网络连接的灌溉系统"""
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei'] # 微软雅黑支持更多符号
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
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], 'co', markersize=15, label='取水点')
# 绘制10米边界线(用于可视化)
ax.plot([MIN_DISTANCE_FROM_BORDER, FARM_SIZE_X-MIN_DISTANCE_FROM_BORDER,
FARM_SIZE_X-MIN_DISTANCE_FROM_BORDER, MIN_DISTANCE_FROM_BORDER,
MIN_DISTANCE_FROM_BORDER],
[MIN_DISTANCE_FROM_BORDER, MIN_DISTANCE_FROM_BORDER,
FARM_SIZE_Y-MIN_DISTANCE_FROM_BORDER, FARM_SIZE_Y-MIN_DISTANCE_FROM_BORDER,
MIN_DISTANCE_FROM_BORDER], 'k--', linewidth=1, 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, row in sprinkler_df.iterrows():
if i not in tank_indices:
ax.plot(row['x'], row['y'], 'bo', markersize=6)
circle = Circle((row['x'], row['y']), SPRINKLER_RADIUS,
color='blue', alpha=0.1, fill=True)
ax.add_patch(circle)
# 绘制储水罐(无主次之分)
for i, idx in enumerate(tank_indices):
row = sprinkler_df.iloc[idx]
ax.plot(row['x'], row['y'], 'rs', markersize=10,
markeredgecolor='black', markeredgewidth=1.5, label='储水罐' if i == 0 else "")
# 绘制储水罐覆盖范围
circle = Circle((row['x'], row['y']), TANK_COVERAGE_RADIUS,
color='red', alpha=0.15, fill=True)
ax.add_patch(circle)
ax.text(row['x'], row['y']+3, f'T{i+1}', fontsize=10,
ha='center', va='bottom', fontweight='bold')
# 标记总管连接的喷头
root_row = sprinkler_df.iloc[root_sprinkler_idx]
ax.plot(root_row['x'], root_row['y'], 'go', markersize=10,
markeredgecolor='black', markeredgewidth=1.5, label='总管连接喷头')
ax.text(root_row['x'], root_row['y']-3, f"S{root_sprinkler_idx}", fontsize=9,
ha='center', va='top', fontweight='bold')
# 绘制网络连接(所有喷头和储水罐之间的连接)
for i, j, data in network_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 "")
# 绘制总管(河流到连接的喷头)
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))
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['network_pipe_cost']
]
colors = ['lightblue', 'lightcoral', 'lightgreen']
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, tank_indices, V_opt, assignments,
Q_total_opt, Q_tank_network_opt, tank_demands, cost_breakdown):
"""输出优化结果表格"""
num_tanks = len(tank_positions)
results_df = pd.DataFrame({
'储水罐编号': range(1, num_tanks+1),
'原喷头编号': [f"S{tank_indices[i]}" for i in range(num_tanks)],
'X坐标': [f"{tank_positions[i][0]:.1f}" for i in range(num_tanks)],
'Y坐标': [f"{tank_positions[i][1]:.1f}" for i in range(num_tanks)],
'容量(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
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)
# 计算储水罐成本明细
base_capacity_used = min(BASE_TANK_CAPACITY * num_tanks, tank_supply)
extra_capacity = max(0, tank_supply - BASE_TANK_CAPACITY * num_tanks)
extra_cost = extra_capacity * EXTRA_TANK_COST
logger.info(f"\n系统总体信息:")
logger.info(f"总管连接喷头编号: {root_sprinkler_idx}")
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" - 免费基础容量: {min(BASE_TANK_CAPACITY * num_tanks, tank_supply):.0f}L")
logger.info(f" - 额外付费容量: {extra_capacity:.0f}L (成本: {extra_cost:.2f}元)")
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['network_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()
# 3. 调整第二行和倒数第二行喷头间距,以中间喷头为参照,增加10米
sprinkler_df = adjust_specific_rows_spacing(sprinkler_df, add_distance=10)
# 4. 计算喷头覆盖面积
sprinkler_df = calculate_sprinkler_coverage(sprinkler_df)
logger.info(f"生成喷头数量: {len(sprinkler_df)}")
# 5. 验证喷头间距
spacing_ok = validate_sprinkler_spacing(sprinkler_df)
if not spacing_ok:
logger.warning("喷头间距不符合要求,可能需要调整布局!")
# 6. 计算灌溉需求
irrigation_demand, sprinkler_df = calculate_irrigation_demand(daily_avg_moisture, sprinkler_df)
# 7. 确定最佳聚类数量
optimal_clusters = determine_optimal_clusters(sprinkler_df, max_clusters=8)
# 8. 从喷头中选择储水罐位置
tank_positions, tank_indices, sprinkler_df = select_tank_positions(sprinkler_df, optimal_clusters)
logger.info(f"储水罐位置 (来自喷头布局): {tank_positions}")
# 9. 网络化优化(修改后的版本)
global root_sprinkler_idx # 声明为全局变量,供output_results使用
result, assignments, tank_demands, tank_mst, network_mst, root_sprinkler_idx, V_opt, Q_total_opt, Q_tank_network_opt = network_optimization(
sprinkler_df, irrigation_demand, tank_positions, tank_indices)
if result.success:
logger.info("优化成功!")
# 10. 计算成本
total_cost, cost_breakdown, Q_total_opt, Q_tank_network_opt = calculate_total_cost(
result, sprinkler_df, tank_positions, assignments, tank_mst, root_sprinkler_idx)
# 11. 可视化
visualize_network_system(sprinkler_df, tank_positions, tank_indices, assignments, tank_mst, network_mst,
root_sprinkler_idx, Q_total_opt, Q_tank_network_opt)
plot_cost_breakdown(cost_breakdown)
# 12. 输出结果
output_results(sprinkler_df, tank_positions, tank_indices, V_opt, assignments,
Q_total_opt, Q_tank_network_opt, tank_demands, cost_breakdown)
# 13. 最终验证报告
logger.info("\n最终系统验证报告:")
logger.info(f"1. 喷头间距验证: {'通过' if spacing_ok else '失败'}")
logger.info(f"2. 储水罐位置: 全部来自喷头布局")
logger.info(f"3. 网络结构: 储水罐与喷头形成统一网络")
logger.info(f"4. 总管连接: 直接连接到最近的喷头 #{root_sprinkler_idx}")
logger.info(f"5. 系统总成本: {total_cost:.2f} 元")
else:
logger.error("优化失败:", result.message)
if __name__ == "__main__":
main()
第二题已经求解完毕,接下来基于第二题来求解第三题
最新发布