---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 727
724 logger.info("结果已保存到 '灌溉系统优化结果.csv'")
726 if __name__ == "__main__":
--> 727 main()
Cell In[15], line 653
650 logger.info(f"直接连接喷头比例: {np.sum(direct_pipe_mask)/len(sprinkler_df)*100:.1f}%")
652 # 8. 可视化
--> 653 visualize_system(sprinkler_df, tank_positions, assignments, direct_pipe_mask, vor)
654 plot_cost_breakdown(tank_cost, pipe_cost)
656 # 9. 储水罐容量分析
Cell In[15], line 556
554 mid_y = (RIVER_POINT[1] + tank[1]) / 2
555 length = np.sqrt((RIVER_POINT[0]-tank[0])**2 + (RIVER_POINT[1]-tank[1])**2)
--> 556 ax.text(mid_x, mid_y, f'{length:.1f}m\n{Q_opt[j]:.0f}L/d', fontsize=8,
557 bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.7))
559 ax.set_xlabel('X坐标 (m)')
560 ax.set_ylabel('Y坐标 (m)')
NameError: name 'Q_opt' is not definedimport 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
from scipy.spatial import Voronoi, voronoi_plot_2d, ConvexHull
import warnings
import logging
from matplotlib.patches import Polygon
import matplotlib.cm as cm
from sklearn.preprocessing import StandardScaler
# 设置日志记录
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 = {
'sorghum': {
'area': 0.5 * 10000, # 公顷转换为平方米
'sowing_water': 5, # 播种期需水量 (L/m²)
'flowering_water': 10, # 开花期需水量 (L/m²)
'mature_water': 8, # 成熟期需水量 (L/m²)
'sowing_days': 20, # 播种期天数
'flowering_days': 50 # 开花期天数
},
'corn': {
'area': 0.3 * 10000,
'sowing_water': 6,
'flowering_water': 12,
'mature_water': 10,
'sowing_days': 32,
'flowering_days': 50
},
'soybean': {
'area': 0.2 * 10000,
'sowing_water': 4,
'flowering_water': 8,
'mature_water': 6,
'sowing_days': 40,
'flowering_days': 40
}
}
# 作物分布（假设按区域分布）
CROP_REGIONS = {
'sorghum': {'x_min': 0, 'x_max': FARM_SIZE_X, 'y_min': 0, 'y_max': FARM_SIZE_Y/2},
'corn': {'x_min': 0, 'x_max': FARM_SIZE_X/2, 'y_min': FARM_SIZE_Y/2, 'y_max': FARM_SIZE_Y},
'soybean': {'x_min': FARM_SIZE_X/2, '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, FARM_SIZE_Y)
# 调整成本公式参数
PIPE_LENGTH_COEFF = 50 # 管道长度系数
PIPE_FLOW_COEFF = 0.1 # 管道流量系数
PIPE_LENGTH_EXP = 1.1 # 管道长度指数 (降低)
PIPE_FLOW_EXP = 1.3 # 管道流量指数 (降低)
TANK_COST_PER_LITER = 5 # 储水罐单位容积成本
# ==================== 数据加载与处理 ====================
def load_soil_moisture_data():
"""
加载土壤湿度数据
"""
try:
# 尝试加载数据
data = pd.read_excel('附件/该地土壤湿度数据.xlsx', sheet_name='JingYueTan')
logger.info(f"土壤湿度数据结构: {data.info()}")
# 转换日期格式
data['DATE'] = pd.to_datetime(data['DATE'])
# 提取7月数据
july_data = data[(data['DATE'] >= '2021-07-01') &
(data['DATE'] <= '2021-07-31')]
# 按日期分组并计算平均土壤湿度
daily_avg_moisture = july_data.groupby('DATE')['5cm_SM'].mean()
# 绘制土壤湿度变化图
plt.figure(figsize=(12, 6))
plt.plot(daily_avg_moisture.index, daily_avg_moisture.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 daily_avg_moisture, july_data
except Exception as e:
logger.error(f"加载土壤湿度数据时出错: {e}")
# 创建模拟数据
dates = pd.date_range('2021-07-01', periods=31)
default_moisture = pd.Series([0.15 + 0.1 * np.sin(i/10) for i in range(31)], index=dates)
# 绘制模拟土壤湿度变化图
plt.figure(figsize=(12, 6))
plt.plot(dates, default_moisture.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 default_moisture, pd.DataFrame({'DATE': dates, '5cm_SM': default_moisture})
def calculate_daily_irrigation_demand(daily_moisture, crop_type=None):
"""
计算每日每平方米灌溉需求，考虑作物类型
"""
base_demand = max(0.0, (MIN_SOIL_MOISTURE - daily_moisture) * dry_mass_per_m2)
# 根据作物类型调整需求
if crop_type and crop_type in CROPS:
# 7月假设为开花期
crop_factor = CROPS[crop_type]['flowering_water'] / 10 # 以高粱的开花需水量10为基准
base_demand *= crop_factor
return base_demand
def calculate_irrigation_demand(daily_avg_moisture, sprinkler_df):
"""
计算每日灌溉需求，考虑作物类型
"""
irrigation_demand = {}
for crop in CROPS.keys():
crop_demand = []
for moisture in daily_avg_moisture:
demand_per_m2 = calculate_daily_irrigation_demand(moisture, crop)
crop_demand.append(demand_per_m2)
irrigation_demand[crop] = crop_demand
# 计算每个喷头的需求
sprinkler_demands = []
for _, sprinkler in sprinkler_df.iterrows():
crop = sprinkler['crop_type']
crop_idx = list(CROPS.keys()).index(crop)
max_demand = max(irrigation_demand[crop])
sprinkler_demands.append(sprinkler['area'] * max_demand)
sprinkler_df['max_demand'] = sprinkler_demands
# 绘制灌溉需求变化图
plt.figure(figsize=(12, 6))
for crop, demand in irrigation_demand.items():
plt.plot(daily_avg_moisture.index, demand, label=f'{crop}需求')
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 irrigation_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_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, vor
# ==================== 多储水罐优化 ====================
def determine_optimal_clusters(sprinkler_df, max_clusters=10):
"""
使用肘部法和轮廓系数确定最佳聚类数量，考虑作物类型
"""
# 创建特征矩阵：位置 + 作物类型（one-hot编码）
coordinates = sprinkler_df[['x', 'y']].values
crop_types = pd.get_dummies(sprinkler_df['crop_type']).values
# 标准化特征
features = np.hstack([coordinates, crop_types])
scaler = StandardScaler()
scaled_features = scaler.fit_transform(features)
# 计算不同聚类数量的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_)
if k > 1: # 轮廓系数需要至少2个聚类
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的下降率
sse_diff = np.diff(sse)
sse_ratio = sse_diff[:-1] / sse_diff[1:]
# 找到肘点 (SSE下降率最大的点)
elbow_point = np.argmax(sse_ratio) + 2 # 从k=2开始
# 找到轮廓系数最大的点
best_silhouette = np.argmax(silhouette_scores[1:]) + 2
# 综合选择最佳聚类数量
optimal_clusters = min(max(3, elbow_point), max(3, best_silhouette)) # 至少3个储水罐
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):
"""
生成候选储水罐位置，考虑作物类型
"""
# 创建特征矩阵：位置 + 作物类型（one-hot编码）
coordinates = sprinkler_df[['x', 'y']].values
crop_types = pd.get_dummies(sprinkler_df['crop_type']).values
# 标准化特征
features = np.hstack([coordinates, crop_types])
scaler = StandardScaler()
scaled_features = scaler.fit_transform(features)
kmeans = KMeans(n_clusters=num_tanks, random_state=42, n_init=10)
kmeans.fit(scaled_features)
# 获取聚类中心的位置部分（忽略作物类型部分）
cluster_centers = kmeans.cluster_centers_[:, :2]
# 反标准化位置
cluster_centers = scaler.inverse_transform(
np.hstack([cluster_centers, np.zeros((num_tanks, crop_types.shape[1]))])
)[:, :2]
return cluster_centers, kmeans.labels_
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)
# 计算每个喷头的最大日需求
sprinkler_max_demands = sprinkler_df['max_demand'].values
# 分配喷头到最近的储水罐（如果在15m内），否则标记为直接连接管道
assignments = np.argmin(distances, axis=1)
direct_pipe_mask = np.min(distances, axis=1) > SPRINKLER_RADIUS
# 限制直接连接管道的比例不超过10%
direct_ratio = np.sum(direct_pipe_mask) / num_sprinklers
if direct_ratio > 0.1:
# 如果直接连接比例超过10%，重新分配最远的喷头
sorted_indices = np.argsort(np.min(distances, axis=1))[::-1] # 从最远到最近
for i in sorted_indices:
if direct_pipe_mask[i] and direct_ratio > 0.1:
direct_pipe_mask[i] = False
direct_ratio = np.sum(direct_pipe_mask) / num_sprinklers
if direct_ratio <= 0.1:
break
# 阶段2: 优化储水罐容量和管道流量
def objective(vars):
V = vars[:num_tanks] # 储水罐容量
Q = vars[num_tanks:num_tanks+num_tanks] # 到每个储水罐的管道流量
# 储水罐成本
tank_cost = TANK_COST_PER_LITER * 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 += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q[j] ** PIPE_FLOW_EXP)
# 惩罚项：储水罐容量不足
penalty = 0
for j in range(num_tanks):
# 计算该储水罐覆盖喷头的总最大需求
tank_demand = np.sum(sprinkler_max_demands[(assignments == j) & ~direct_pipe_mask])
if V[j] < tank_demand:
penalty += 1000 * (tank_demand - V[j])
# 惩罚项：直接连接管道过多
direct_penalty = 10000 * max(0, np.sum(direct_pipe_mask) / num_sprinklers - 0.1)
return tank_cost + pipe_cost + penalty + direct_penalty
# 约束条件
constraints = []
# 流量约束：管道流量应大于等于储水罐需求
for j in range(num_tanks):
tank_demand = np.sum(sprinkler_max_demands[(assignments == j) & ~direct_pipe_mask])
constraints.append({'type': 'ineq', 'fun': lambda x, j=j, td=tank_demand: x[num_tanks+j] - td})
# 容量约束
for j in range(num_tanks):
constraints.append({'type': 'ineq', 'fun': lambda x, j=j: x[j]}) # V[j] ≥ 0
# 初始值
initial_V = [np.sum(sprinkler_max_demands[(assignments == j) & ~direct_pipe_mask])
for j in range(num_tanks)]
initial_Q = [1.2 * v for v in initial_V] # 流量比容量多20%以提供缓冲
x0 = initial_V + initial_Q
# 边界
bounds = Bounds([0] * (2 * num_tanks), [np.inf] * (2 * num_tanks))
# 使用SLSQP进行优化
logger.info("开始优化...")
result = minimize(
objective,
x0,
bounds=bounds,
constraints=constraints,
method='SLSQP',
options={'disp': True, 'ftol': 1e-6, 'maxiter': 100}
)
return result, assignments, direct_pipe_mask
# ==================== 成本计算与可视化 ====================
def calculate_total_cost(result, sprinkler_df, tank_positions, assignments, direct_pipe_mask):
"""
计算总成本
"""
num_tanks = len(tank_positions)
V_opt = result.x[:num_tanks]
Q_opt = result.x[num_tanks:2*num_tanks]
# 储水罐成本
tank_cost = TANK_COST_PER_LITER * 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 += PIPE_LENGTH_COEFF * (L ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_opt[j] ** PIPE_FLOW_EXP)
# 直接连接管道的喷头成本
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]['max_demand'])
pipe_cost += PIPE_LENGTH_COEFF * (L_direct ** PIPE_LENGTH_EXP) + PIPE_FLOW_COEFF * (Q_direct ** PIPE_FLOW_EXP)
total_cost = tank_cost + pipe_cost
return total_cost, tank_cost, pipe_cost
def visualize_system(sprinkler_df, tank_positions, assignments, direct_pipe_mask, vor):
"""
可视化灌溉系统
"""
fig, ax = plt.subplots(figsize=(14, 12))
# 绘制农场边界
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 == 'south':
ax.plot([0, FARM_SIZE_X], [FARM_SIZE_Y, FARM_SIZE_Y], 'b-', linewidth=4, label='河流')
# 绘制作物区域
colors = {'sorghum': 'lightgreen', 'corn': 'lightyellow', 'soybean': 'lightblue'}
for crop, region in CROP_REGIONS.items():
rect = plt.Rectangle((region['x_min'], region['y_min']),
region['x_max'] - region['x_min'],
region['y_max'] - region['y_min'],
alpha=0.3, color=colors[crop], label=crop)
ax.add_patch(rect)
# 绘制Voronoi图
voronoi_plot_2d(vor, ax=ax, show_points=False, show_vertices=False, line_colors='gray', line_width=0.5, alpha=0.5)
# 绘制喷头
tank_colors = cm.rainbow(np.linspace(0, 1, len(tank_positions)))
for i, sprinkler in sprinkler_df.iterrows():
if direct_pipe_mask[i]:
# 直接连接管道的喷头
ax.plot(sprinkler['x'], sprinkler['y'], 'ko', markersize=4)
ax.add_patch(plt.Circle((sprinkler['x'], sprinkler['y']),
SPRINKLER_RADIUS, color='gray', alpha=0.1))
else:
# 由储水罐覆盖的喷头
tank_id = assignments[i]
color = tank_colors[tank_id]
ax.plot(sprinkler['x'], sprinkler['y'], 'o', color=color, markersize=4)
ax.add_patch(plt.Circle((sprinkler['x'], sprinkler['y']),
SPRINKLER_RADIUS, color=color, alpha=0.1))
# 绘制储水罐
for j, tank in enumerate(tank_positions):
color = tank_colors[j]
ax.plot(tank[0], tank[1], 's', color=color, markersize=10, markeredgecolor='black', markeredgewidth=2)
ax.add_patch(plt.Circle((tank[0], tank[1]), SPRINKLER_RADIUS*2,
color=color, alpha=0.2))
ax.text(tank[0], tank[1], f'T{j+1}', fontsize=12, ha='center', va='center', fontweight='bold')
# 绘制管道
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 "")
# 标注管道长度和流量
mid_x = (RIVER_POINT[0] + tank[0]) / 2
mid_y = (RIVER_POINT[1] + tank[1]) / 2
length = np.sqrt((RIVER_POINT[0]-tank[0])**2 + (RIVER_POINT[1]-tank[1])**2)
ax.text(mid_x, mid_y, f'{length:.1f}m\n{Q_opt[j]:.0f}L/d', fontsize=8,
bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.7))
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 plot_cost_breakdown(tank_cost, pipe_cost):
"""
绘制成本分解图
"""
labels = ['储水罐成本', '管道成本']
sizes = [tank_cost, pipe_cost]
colors = ['lightblue', 'lightcoral']
plt.figure(figsize=(8, 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 plot_tank_capacity_analysis(V_opt, tank_demands):
"""
绘制储水罐容量分析图
"""
num_tanks = len(V_opt)
x = np.arange(num_tanks)
width = 0.35
fig, ax = plt.subplots(figsize=(10, 6))
rects1 = ax.bar(x - width/2, V_opt, width, label='设计容量')
rects2 = ax.bar(x + width/2, tank_demands, width, label='需求容量')
ax.set_xlabel('储水罐编号')
ax.set_ylabel('容量 (L)')
ax.set_title('储水罐容量与需求对比')
ax.set_xticks(x)
ax.set_xticklabels([f'T{i+1}' for i in range(num_tanks)])
ax.legend()
# 添加数值标签
for i, v in enumerate(V_opt):
ax.text(i - width/2, v + max(V_opt)*0.01, f'{v:.0f}', ha='center')
for i, v in enumerate(tank_demands):
ax.text(i + width/2, v + max(V_opt)*0.01, f'{v:.0f}', ha='center')
plt.tight_layout()
plt.savefig('储水罐容量分析.png', dpi=300)
plt.show()
# ==================== 主函数 ====================
def main():
"""
主优化流程
"""
logger.info("开始农业灌溉系统优化...")
# 1. 加载和处理数据
daily_avg_moisture, july_data = load_soil_moisture_data()
# 2. 生成喷头布局
sprinkler_df = generate_sprinkler_layout()
sprinkler_df, vor = calculate_sprinkler_coverage(sprinkler_df)
logger.info(f"生成喷头数量: {len(sprinkler_df)}")
# 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, direct_pipe_mask = two_stage_optimization(
sprinkler_df, irrigation_demand, tank_positions)
if result.success:
logger.info("优化成功!")
# 7. 计算成本
total_cost, tank_cost, pipe_cost = calculate_total_cost(
result, sprinkler_df, tank_positions, assignments, direct_pipe_mask)
logger.info(f"总成本: {total_cost:.2f} 元")
logger.info(f"储水罐成本: {tank_cost:.2f} 元 ({tank_cost/total_cost*100:.1f}%)")
logger.info(f"管道成本: {pipe_cost:.2f} 元 ({pipe_cost/total_cost*100:.1f}%)")
logger.info(f"直接连接喷头比例: {np.sum(direct_pipe_mask)/len(sprinkler_df)*100:.1f}%")
# 8. 可视化
visualize_system(sprinkler_df, tank_positions, assignments, direct_pipe_mask, vor)
plot_cost_breakdown(tank_cost, pipe_cost)
# 9. 储水罐容量分析
num_tanks = len(tank_positions)
V_opt = result.x[:num_tanks]
tank_demands = []
for j in range(num_tanks):
tank_demand = np.sum(sprinkler_df['max_demand'][(assignments == j) & ~direct_pipe_mask])
tank_demands.append(tank_demand)
plot_tank_capacity_analysis(V_opt, tank_demands)
# 10. 输出结果表格
output_results(sprinkler_df, tank_positions, result, assignments, direct_pipe_mask)
else:
logger.error("优化失败:", result.message)
def output_results(sprinkler_df, tank_positions, result, assignments, direct_pipe_mask):
"""
输出优化结果表格
"""
num_tanks = len(tank_positions)
V_opt = result.x[:num_tanks]
Q_opt = result.x[num_tanks:2*num_tanks]
# 计算每个储水罐的需求
tank_demands = []
for j in range(num_tanks):
tank_demand = np.sum(sprinkler_df['max_demand'][(assignments == j) & ~direct_pipe_mask])
tank_demands.append(tank_demand)
# 创建结果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"{d:.0f}" for d in tank_demands],
'安全余量(%)': [f"{(v-d)/d*100:.1f}" if d > 0 else "N/A" for v, d in zip(V_opt, tank_demands)],
'管道流量(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]['max_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}"
logger.info("\n优化结果详情:")
print(results_df.to_string(index=False))
# 保存结果
results_df.to_csv('灌溉系统优化结果.csv', index=False, encoding='utf-8-sig')
logger.info("结果已保存到 '灌溉系统优化结果.csv'")
if __name__ == "__main__":
main()修改返回完整代码
最新发布
扫一扫
2471
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
