import numpy as np
import pandas as pd
from scipy.optimize import minimize, Bounds
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi, voronoi_plot_2d, ConvexHull
from geopy.distance import geodesic
import warnings
# 忽略特定警告
warnings.filterwarnings("ignore", category=UserWarning, module="sklearn")
# ==================== 常量定义 ====================
DEPTH = 0.05 # 5 cm
DRY_DENS = 1500 # kg/m³
dry_mass_per_m2 = DRY_DENS * DEPTH # 75 kg/m²
# 根据坐标计算农场实际尺寸
def calculate_farm_dimensions():
"""
根据采样点坐标计算农场实际尺寸
"""
# 坐标点(经度, 纬度)
points = [
(125.621498, 44.791779), # 坐标1
(125.622031, 44.791988), # 坐标2
(125.622885, 44.790893), # 坐标3
(125.622354, 44.790613) # 坐标4
]
# 计算东西方向和南北方向的距离
east_west = geodesic((points[0][1], points[0][0]), (points[2][1], points[2][0])).meters
north_south = geodesic((points[0][1], points[0][0]), (points[1][1], points[1][0])).meters
return east_west, north_south
# 计算农场尺寸
FARM_SIZE_X, FARM_SIZE_Y = calculate_farm_dimensions()
print(f"农场实际尺寸: {FARM_SIZE_X:.1f}m × {FARM_SIZE_Y:.1f}m")
SPRINKLER_RADIUS = 15 # 喷头半径
RIVER_POSITION = 'south'
# 根据农场尺寸调整河流取水点
RIVER_POINT = (FARM_SIZE_X / 2, FARM_SIZE_Y) # 河流取水点
# ==================== 数据加载与处理 ====================
def load_soil_moisture_data():
"""
加载多个采样点的土壤湿度数据并计算平均值
"""
try:
# 尝试加载所有采样点的数据
soil_data_list = []
for i in range(1, 5):
sheet_name = f'Sample area coordinate {i}'
data = pd.read_excel('附件\该地土壤湿度数据.xlsx', sheet_name='JingYueTan')
soil_data_list.append(data)
# 合并数据并计算平均值
combined_data = pd.concat(soil_data_list)
july_data = combined_data[(combined_data['DATE'] >= '2021-07-01') &
(combined_data['DATE'] <= '2021-07-31')]
# 按日期分组并计算平均土壤湿度
daily_avg_moisture = july_data.groupby('DATE')['5cm_SM'].mean()
return daily_avg_moisture
def calculate_irrigation_demand(daily_avg_moisture):
"""
计算每日灌溉需求
"""
irrigation_demand = []
for moisture in daily_avg_moisture:
demand_per_m2 = max(0.0, (0.22 - moisture) * dry_mass_per_m2)
irrigation_demand.append(demand_per_m2)
return irrigation_demand
# ==================== 喷头布局生成 ====================
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:
sprinklers.append({
'id': len(sprinklers),
'x': x,
'y': y,
'radius': radius
})
return pd.DataFrame(sprinklers)
def calculate_sprinkler_coverage(sprinkler_df, farm_size_x=FARM_SIZE_X, farm_size_y=FARM_SIZE_Y):
"""
计算每个喷头的覆盖面积(使用Voronoi图划分)
"""
points = sprinkler_df[['x', 'y']].values
# 创建Voronoi图
vor = Voronoi(points)
# 计算每个Voronoi单元的面积
areas = []
for region in vor.regions:
if not -1 in region and len(region) > 0:
polygon = vor.vertices[region]
try:
hull = ConvexHull(polygon)
area = hull.volume # 对于2D,volume是面积
areas.append(area)
except:
# 对于无效区域使用默认面积
areas.append((farm_size_x * farm_size_y) / len(sprinkler_df))
else:
areas.append((farm_size_x * farm_size_y) / len(sprinkler_df))
# 确保面积不为负且不超过农场总面积
areas = np.clip(areas, 0, farm_size_x * farm_size_y)
sprinkler_df['area'] = areas[:len(sprinkler_df)]
return sprinkler_df
# ==================== 多储水罐优化 ====================
def generate_candidate_tanks(sprinkler_df, num_tanks=3):
"""
生成候选储水罐位置
"""
coordinates = sprinkler_df[['x', 'y']].values
kmeans = KMeans(n_clusters=num_tanks, random_state=42, n_init=10)
kmeans.fit(coordinates)
return kmeans.cluster_centers_
def calculate_distance_matrix(sprinkler_df, tank_positions):
"""
计算喷头到储水罐的距离矩阵
"""
num_sprinklers = len(sprinkler_df)
num_tanks = len(tank_positions)
distances = np.zeros((num_sprinklers, num_tanks))
for i, sprinkler in sprinkler_df.iterrows():
for j, tank in enumerate(tank_positions):
dist = np.sqrt((sprinkler['x'] - tank[0])**2 +
(sprinkler['y'] - tank[1])**2)
distances[i, j] = dist
return distances
def two_stage_optimization(sprinkler_df, irrigation_demand, tank_positions):
"""
两阶段优化:先分配喷头,再优化容量和管道
"""
# 阶段1: 喷头分配到储水罐或管道
distances = calculate_distance_matrix(sprinkler_df, tank_positions)
num_sprinklers = len(sprinkler_df)
num_tanks = len(tank_positions)
# 计算每个喷头的最大日需求
max_demand = max(irrigation_demand)
sprinkler_demands = sprinkler_df['area'] * max_demand
# 分配喷头到最近的储水罐(如果在15m内),否则标记为直接连接管道
assignments = np.argmin(distances, axis=1)
direct_pipe_mask = np.min(distances, axis=1) > SPRINKLER_RADIUS
# 阶段2: 优化储水罐容量和管道流量
def objective(vars):
V = vars[:num_tanks] # 储水罐容量
Q = vars[num_tanks:num_tanks+num_tanks] # 到每个储水罐的管道流量
# 储水罐成本
tank_cost = 5 * np.sum(V)
# 管道成本
pipe_cost = 0
for j in range(num_tanks):
L = np.sqrt((RIVER_POINT[0] - tank_positions[j][0])**2 +
(RIVER_POINT[1] - tank_positions[j][1])**2)
pipe_cost += 50 * (L ** 1.2) + 0.1 * (Q[j] ** 1.5)
# 惩罚项:储水罐容量不足
penalty = 0
for j in range(num_tanks):
tank_demand = np.sum(sprinkler_demands[(assignments == j) & ~direct_pipe_mask])
if V[j] < tank_demand:
penalty += 1000 * (tank_demand - V[j])
return tank_cost + pipe_cost + penalty
# 约束条件
def constraints(vars):
V = vars[:num_tanks]
Q = vars[num_tanks:num_tanks+num_tanks]
cons = []
# 流量约束:管道流量应大于等于储水罐需求
for j in range(num_tanks):
tank_demand = np.sum(sprinkler_demands[(assignments == j) & ~direct_pipe_mask])
cons.append(Q[j] - tank_demand)
# 容量约束
for j in range(num_tanks):
cons.append(V[j]) # V[j] ≥ 0
return np.array(cons)
# 初始值
initial_V = [np.sum(sprinkler_demands[(assignments == j) & ~direct_pipe_mask])
for j in range(num_tanks)]
initial_Q = initial_V.copy()
x0 = initial_V + initial_Q
# 边界
bounds = Bounds([0] * (2 * num_tanks), [np.inf] * (2 * num_tanks))
# 优化
result = minimize(objective, x0, bounds=bounds,
constraints={'type': 'ineq', 'fun': constraints},
method='SLSQP', options={'disp': True})
return result, assignments, direct_pipe_mask
# ==================== 成本计算与可视化 ====================
def calculate_total_cost(result, sprinkler_df, tank_positions, assignments, direct_pipe_mask, irrigation_demand):
"""
计算总成本
"""
num_tanks = len(tank_positions)
V_opt = result.x[:num_tanks]
Q_opt = result.x[num_tanks:2*num_tanks]
# 储水罐成本
tank_cost = 5 * np.sum(V_opt)
# 管道成本
pipe_cost = 0
for j in range(num_tanks):
L = np.sqrt((RIVER_POINT[0] - tank_positions[j][0])**2 +
(RIVER_POINT[1] - tank_positions[j][1])**2)
pipe_cost += 50 * (L ** 1.2) + 0.1 * (Q_opt[j] ** 1.5)
# 直接连接管道的喷头成本
direct_sprinklers = sprinkler_df[direct_pipe_mask]
if len(direct_sprinklers) > 0:
# 计算主管道长度(到直接连接喷头的中心)
center_x = np.mean(direct_sprinklers['x'])
center_y = np.mean(direct_sprinklers['y'])
L_direct = np.sqrt((RIVER_POINT[0] - center_x)**2 +
(RIVER_POINT[1] - center_y)**2)
Q_direct = np.sum(sprinkler_df[direct_pipe_mask]['area']) * max(irrigation_demand)
pipe_cost += 50 * (L_direct ** 1.2) + 0.1 * (Q_direct ** 1.5)
total_cost = tank_cost + pipe_cost
return total_cost, tank_cost, pipe_cost
def visualize_system(sprinkler_df, tank_positions, assignments, direct_pipe_mask):
"""
可视化灌溉系统
"""
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei'] # 微软雅黑支持更多符号
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
fig, ax = plt.subplots(figsize=(12, 10))
# 绘制农场边界
ax.plot([0, FARM_SIZE_X, FARM_SIZE_X, 0, 0],
[0, 0, FARM_SIZE_Y, FARM_SIZE_Y, 0], 'k-', linewidth=2)
# 绘制河流
if RIVER_POSITION == 'north':
ax.plot([0, FARM_SIZE_X], [FARM_SIZE_Y, FARM_SIZE_Y], 'b-', linewidth=3, label='河流')
# 绘制喷头
colors = ['red', 'green', 'blue', 'orange', 'purple']
for i, sprinkler in sprinkler_df.iterrows():
if direct_pipe_mask[i]:
# 直接连接管道的喷头
ax.plot(sprinkler['x'], sprinkler['y'], 'ko', markersize=5)
ax.add_patch(plt.Circle((sprinkler['x'], sprinkler['y']),
SPRINKLER_RADIUS, color='gray', alpha=0.1))
else:
# 由储水罐覆盖的喷头
tank_id = assignments[i]
color = colors[tank_id % len(colors)]
ax.plot(sprinkler['x'], sprinkler['y'], 'o', color=color, markersize=5)
ax.add_patch(plt.Circle((sprinkler['x'], sprinkler['y']),
SPRINKLER_RADIUS, color=color, alpha=0.1))
# 绘制储水罐
for j, tank in enumerate(tank_positions):
color = colors[j % len(colors)]
ax.plot(tank[0], tank[1], 's', color=color, markersize=10, label=f'储水罐 {j+1}')
ax.add_patch(plt.Circle((tank[0], tank[1]), SPRINKLER_RADIUS*2,
color=color, alpha=0.2))
# 绘制管道
for j, tank in enumerate(tank_positions):
ax.plot([RIVER_POINT[0], tank[0]], [RIVER_POINT[1], tank[1]],
'r-', linewidth=2, label='管道' if j == 0 else "")
ax.set_xlabel('X坐标 (m)')
ax.set_ylabel('Y坐标 (m)')
ax.set_title('农业灌溉系统布局优化')
ax.legend()
ax.grid(True)
plt.savefig('灌溉系统布局图.png', dpi=300)
plt.show()
# ==================== 主函数 ====================
def main():
"""
主优化流程
"""
print("开始农业灌溉系统优化...")
# 1. 加载和处理数据
daily_avg_moisture = load_soil_moisture_data()
irrigation_demand = calculate_irrigation_demand(daily_avg_moisture)
print(f"最大日灌溉需求: {max(irrigation_demand):.2f} L/m²")
# 2. 生成喷头布局
sprinkler_df = generate_sprinkler_layout()
sprinkler_df = calculate_sprinkler_coverage(sprinkler_df)
print(f"生成喷头数量: {len(sprinkler_df)}")
# 3. 生成候选储水罐位置
tank_positions = generate_candidate_tanks(sprinkler_df, num_tanks=3)
print(f"候选储水罐位置: {tank_positions}")
# 4. 两阶段优化
result, assignments, direct_pipe_mask = two_stage_optimization(
sprinkler_df, irrigation_demand, tank_positions)
if result.success:
print("优化成功!")
# 5. 计算成本
total_cost, tank_cost, pipe_cost = calculate_total_cost(
result, sprinkler_df, tank_positions, assignments, direct_pipe_mask, irrigation_demand)
print(f"总成本: {total_cost:.2f} 元")
print(f"储水罐成本: {tank_cost:.2f} 元")
print(f"管道成本: {pipe_cost:.2f} 元")
# 6. 可视化
visualize_system(sprinkler_df, tank_positions, assignments, direct_pipe_mask)
# 7. 输出结果表格
output_results(sprinkler_df, tank_positions, result, assignments, direct_pipe_mask, irrigation_demand)
else:
print("优化失败:", result.message)
def output_results(sprinkler_df, tank_positions, result, assignments, direct_pipe_mask, irrigation_demand):
"""
输出优化结果表格
"""
num_tanks = len(tank_positions)
V_opt = result.x[:num_tanks]
Q_opt = result.x[num_tanks:2*num_tanks]
# 创建结果DataFrame
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],
'管道流量(L/天)': [f"{q:.0f}" for q in Q_opt],
'管道长度(m)': [f"{np.sqrt((RIVER_POINT[0]-tank[0])**2 + (RIVER_POINT[1]-tank[1])**2):.1f}"
for tank in tank_positions]
})
# 计算覆盖统计
coverage_stats = []
for j in range(num_tanks):
covered = np.sum((assignments == j) & ~direct_pipe_mask)
coverage_stats.append(covered)
results_df['覆盖喷头数'] = coverage_stats
results_df['直接连接喷头数'] = np.sum(direct_pipe_mask)
# 添加直接连接喷头信息
direct_sprinklers = sprinkler_df[direct_pipe_mask]
if len(direct_sprinklers) > 0:
center_x = np.mean(direct_sprinklers['x'])
center_y = np.mean(direct_sprinklers['y'])
L_direct = np.sqrt((RIVER_POINT[0] - center_x)**2 +
(RIVER_POINT[1] - center_y)**2)
Q_direct = np.sum(sprinkler_df[direct_pipe_mask]['area']) * max(irrigation_demand)
results_df['直接连接中心X'] = f"{center_x:.1f}"
results_df['直接连接中心Y'] = f"{center_y:.1f}"
results_df['直接连接管道长度'] = f"{L_direct:.1f}"
results_df['直接连接流量'] = f"{Q_direct:.0f}"
print("\n优化结果详情:")
print(results_df.to_string(index=False))
# 保存结果
results_df.to_csv('灌溉系统优化结果.csv', index=False, encoding='utf-8-sig')
print("结果已保存到 '灌溉系统优化结果.csv'")
if __name__ == "__main__":
main() Cell In[9], line 70
def calculate_irrigation_demand(daily_avg_moisture):
^
SyntaxError: expected 'except' or 'finally' block,修改错误,并返回
最新发布