L - Farm Irrigation

Description

Benny has a spacious farm land to irrigate. The farm land is a rectangle, and is divided into a lot of samll squares. Water pipes are placed in these squares. Different square has a different type of pipe. There are 11 types of pipes, which is marked from A to K, as Figure 1 shows.


Figure 1


Benny has a map of his farm, which is an array of marks denoting the distribution of water pipes over the whole farm. For example, if he has a map 

ADC
FJK
IHE

then the water pipes are distributed like 


Figure 2


Several wellsprings are found in the center of some squares, so water can flow along the pipes from one square to another. If water flow crosses one square, the whole farm land in this square is irrigated and will have a good harvest in autumn. 

Now Benny wants to know at least how many wellsprings should be found to have the whole farm land irrigated. Can you help him? 

Note: In the above example, at least 3 wellsprings are needed, as those red points in Figure 2 show.
 
 
 

Input

There are several test cases! In each test case, the first line contains 2 integers M and N, then M lines follow. In each of these lines, there are N characters, in the range of 'A' to 'K', denoting the type of water pipe over the corresponding square. A negative M or N denotes the end of input, else you can assume 1 <= M, N <= 50.

Output

For each test case, output in one line the least number of wellsprings needed.
 

Sample Input

2 2 DK HF 3 3 ADC FJK IHE -1 -1
 

Sample Output

2
3

#include <cstdio>
#include <cstring>
char map[60][60];
int d[3600];
int path[11][4]=           //上下左右为顺序,1表示有水管,0表示没有 
   { {1,0,1,0},{1,0,0,1},{0,1,1,0},{0,1,0,1},{1,1,0,0},{0,0,1,1},{1,0,1,1},{1,1,1,0},{0,1,1,1},{1,1,0,1},{1,1,1,1} };
int n,m,sum;
int find(int x)
{
    while(x!=d[x]) 
       x=d[x];
    return x;
}
int merge(int x,int y)                 //如果并入一个集合返回1,sum-- 
{
    int fx=find(x);
    int fy=find(y);
    if(fx==fy) return 0;
    d[fx]=fy;
    return 1;
}
void search(int x,int y)
{
    if(x>0 && path[map[x][y]-'A'][0] && path[map[x-1][y]-'A'][1])   //该位置的"上"和同列上一行位置的"下"都是1时相连 
        sum-=merge(x*m+y,(x-1)*m+y);
    if(y>0 && path[map[x][y]-'A'][2] && path[map[x][y-1]-'A'][3])   //该位置的"左"和同行前一列位置的"右"都是1时相连 
        sum-=merge(x*m+y,x*m+y-1);
}
int main()
{
    while(scanf("%d%d",&n,&m),n!=-1&&m!=-1)
    {
    
    	for(int i=0;i<3600;i++)
    	   d[i]=i;
 	   
        for(int i=0;i<n;i++)  
            scanf("%s",map[i]);
        
        sum=n*m;                    //sum初始化为n*m表示每个正方形为一个集合 
       for(int i=0;i<n;i++)
            for(int j=0;j<m;j++)
            {
                search(i,j);       //重点与技巧 
            }

        printf("%d\n",sum);
    } 
}


--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[4], line 1011 1008 logger.error("优化失败:", result.message) 1010 if __name__ == "__main__": -> 1011 main() Cell In[4], line 985 982 logger.info(f"候选储水罐位置: {tank_positions}") 984 # 6. 网络化优化 --> 985 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( 986 sprinkler_df, irrigation_demand, tank_positions) 988 if result.success: 989 logger.info("优化成功!") Cell In[4], line 563 560 tank_mst = list(nx.minimum_spanning_edges(tank_G, weight='length', data=True)) 562 # 找到离河流最近的储水罐作为储水罐网络入口 --> 563 distances_to_river_tank = [np.sqrt((x - RIVER_POoint[0])**2 + (y - RIVER_POINT[1])**2) 564 for x, y in tank_positions] 565 root_tank_idx = np.argmin(distances_to_river_tank) 567 # 分配喷头到储水罐(基于最近距离) NameError: name 'RIVER_POoint' 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 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] # 创建标签(所有喷头都分配到最近的储水罐) 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_POoint[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.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) # 绘制喷头覆盖范围 ax.add_patch(Circle((sprinkler['x'], sprinkler['y']), SPRINKLER_RADIUS, color='blue', alpha=0.1)) # 绘制喷头之间的连接(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) # 绘制储水罐覆盖范围(问题2为15m) ax.add_patch(Circle((tank[0], tank[1]), TANK_COVERAGE_RADIUS, color='purple', alpha=0.15)) 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()修改返回完整代码
最新发布
08-25
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[1], line 370 367 print("结果已保存到 '灌溉系统优化结果.csv'") 369 if __name__ == "__main__": --> 370 main() Cell In[1], line 319 317 print("优化成功!") 318 # 5. 计算成本 --> 319 total_cost, tank_cost, pipe_cost = calculate_total_cost( 320 result, sprinkler_df, tank_positions, assignments, direct_pipe_mask) 322 print(f"总成本: {total_cost:.2f} 元") 323 print(f"储水罐成本: {tank_cost:.2f} 元") Cell In[1], line 234 231 center_y = np.mean(direct_sprinklers['y']) 232 L_direct = np.sqrt((RIVER_POINT[0] - center_x)**2 + 233 (RIVER_POINT[1] - center_y)**2) --> 234 Q_direct = np.sum(sprinkler_df[direct_pipe_mask]['area']) * max(irrigation_demand) 235 pipe_cost += 50 * (L_direct ** 1.2) + 0.1 * (Q_direct ** 1.5) 237 total_cost = tank_cost + pipe_cost NameError: name 'irrigation_demand' is not definedimport numpy as np import pandas as pd from scipy.optimize import minimize, Bounds, LinearConstraint from sklearn.cluster import KMeans import matplotlib.pyplot as plt from scipy.spatial import Voronoi, voronoi_plot_2d # ==================== 常量定义 ==================== DEPTH = 0.05 # 5 cm DRY_DENS = 1500 # kg/m³ dry_mass_per_m2 = DRY_DENS * DEPTH # 75 kg/m² FARM_SIZE = 100 # 农场大小 100m × 100m SPRINKLER_RADIUS = 15 # 喷头半径 RIVER_POINT = (50, 100) # 河流取水点 # ==================== 数据加载与处理 ==================== def load_soil_moisture_data(): """ 加载土壤湿度数据 """ try: soil_data = pd.read_excel('附件\土壤湿度数据.xlsx') july_data = soil_data[(soil_data['日期'] >= '2021-07-01') & (soil_data['日期'] <= '2021-07-31')] daily_min_moisture = july_data.groupby('日期')['5cm_SM'].min() return daily_min_moisture except: # 如果无法加载数据,使用示例数据 print("使用示例土壤湿度数据") dates = pd.date_range('2021-07-01', '2021-07-31') return pd.Series(np.random.uniform(0.18, 0.25, len(dates)), index=dates) def calculate_irrigation_demand(daily_min_moisture): """ 计算每日灌溉需求 """ irrigation_demand = [] for moisture in daily_min_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=FARM_SIZE, radius=SPRINKLER_RADIUS): """ 生成六角形喷头布局 """ spacing = radius * 1.5 # 喷头间距 sprinklers = [] # 六角形网格生成 rows = int(farm_size / (spacing * np.sqrt(3)/2)) + 2 cols = int(farm_size / 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 and 0 <= y <= farm_size: sprinklers.append({ 'id': len(sprinklers), 'x': x, 'y': y, 'radius': radius }) return pd.DataFrame(sprinklers) def calculate_sprinkler_coverage(sprinkler_df, farm_size=FARM_SIZE): """ 计算每个喷头的覆盖面积(使用Voronoi图划分) """ points = sprinkler_df[['x', 'y']].values # 创建Voronoi图 vor = Voronoi(points) # 计算每个Voronoi单元的面积 from scipy.spatial import ConvexHull 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(0) else: areas.append(0) # 确保面积不为负且不超过农场总面积 areas = np.clip(areas, 0, farm_size**2) 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) 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): """ 计算总成本 """ 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=(10, 10)) # 绘制农场边界 ax.plot([0, FARM_SIZE, FARM_SIZE, 0, 0], [0, 0, FARM_SIZE, FARM_SIZE, 0], 'k-', linewidth=2) # 绘制河流 ax.plot([0, FARM_SIZE], [FARM_SIZE, FARM_SIZE], '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, 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_min_moisture = load_soil_moisture_data() irrigation_demand = calculate_irrigation_demand(daily_min_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) 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) else: print("优化失败:", 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] # 创建结果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) 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()
08-24
--------------------------------------------------------------------------- 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()修改返回完整代码
08-24
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值