2Sum closest vs. 2Minus closest

最新推荐文章于 2019-08-23 14:45:44 发布
原创 最新推荐文章于 2019-08-23 14:45:44 发布 · 896 阅读 · 0 ·
CC 4.0 BY-SA版权
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

一 如果去掉closest,，即2sum 和2minus问题

1  这两个题都可以转化成查找问题，枚举其中一个数a，查找另一个，对于2sum 就是查找sum-a，对于2minus就是查找a-diff和a+diff，用hashmap实现O(n)。

2 也可以先排序再用双指针法，注意2sum的双指针是两头夹逼，2minus的双指针是同向移动（即滑动窗口法）


二 如果是closest 问题

hasmap法不可用，因为无法转化为查找问题，两个数都只能是枚举。

1 对于2Sum closest，排序后用双指针两头夹逼法枚举两个数，

2 对于2Minus closest，排序后用同向移动的双指针法枚举两个数，注意，如果diff > target, L++的时候，如果和R重合，要推着R走

如果是求差之最接近0的两个数，可以排序后枚举相邻元素的差，求差最小的两个相邻数。



2sum, 3sum, 4sum和3sum closest总结
IT成长之路
12-13 260
[size=medium] 2sum, 3sum, 4sum和3sum closest是四道类似的题目，后面三个都是2sum的延伸。[/size] [b]1，2sum[/b] 给定一个数组numbers[]和一个目标数值target，从数组中找到两个元素，使它们的和等于目标元素。返回这两个值的索引indices，并且index1 大于index2。（假设只有一个结果） 例如：输入: nu...
LeetCode - Two Sum I - II - III - IV、3Sum、3Sum Closest、4Sum、4Sum II 系列学习总结
Bob__yuan的博客
09-06 394
目录 1 - Two Sum 2 - Two Sum II - Input array is sorted 3 - Two Sum III - Data structure design 4 - Two Sum IV - Input is a BST 5 - 3Sum 6 - 3Sum Closest 7 - 3Sum Smaller 8 - 4Sum 9 - 4Sum II ...
Leetcode ——2sum, 3sum, 3sumcloset. 4sum
qy724728631的博客
08-23 205
1. Two Sum Given an array of integers, return indices of the two numbers such that they add up to a specific target. You may assume that each input would have exactly one solution, and you may not u...
K Sum(2 Sum,3 Sum,4 Sum,3-Sum Closest)
weixin_34352449的博客
10-10 122
算是经典算法问题了。这里主要针对只存在一个解或者只需要求一个解的情况描述一下解题思路。若需要找到所有可能解，方法需要略作调整。如有问题，欢迎指正。 2 sum： 如果已排序，可直接用夹逼法，即两指针从头尾向中间移动，使和靠近target。时间复杂度为O(n)。 若未排序，因为除非有特别限制，比较排序的时间复杂度为O(nlgn)。此时用hash更为合适，可以达到O(n)的时间复杂度。此方法实际...
2sum问题和3sum问题
kunpengku
02-21 1506
算法 这个书，190页介绍了 这两个问题 。2sum的意思是 在一组数中，找到 两个数的和为零。有多少个这样的组合。 3sum是 找 有多少三个数的组合 ，他们的和为零。常规的for循环，问题的规模分别是 平方级别 和立方级别。使用先排序，然后查找的方式，可以是 对数线性， 和 平方对数级别。 这样就快了一些。这个转变是 ，把求一个组合的问题 ，转化成 已知一个数，找另一个数的问题。 查找就先排
16. 3Sum Closest
Jingwang Li
09-24 184
16. 3Sum Closest 和之前3Sum那道题思路一样，刚开始想的是不用每次计算误差，只需要在小于目标和大于目标的转折点出计算两次然后临界的时候计算一次即可，可是实际上和每次都计算误差速度一样的。。。 class Solution { public: int threeSumClosest(vector<int>& nums, int target) { ...
LeetCode 16. 3Sum Closest
小白菜又菜
08-05 163
question find three numbers whose sum is closed to the target number. algorithm Enumerate the first number, then get two numbers whose sum is closed to the target number minus the first number. code c...
【LeetCode】16. 3Sum Closest
A_Huiii的博客
08-23 129
题目： Given an array nums of n integers and an integer target, find three integers in nums such that the sum is closest to target. Return the sum of the three integers. You may assume that each input wo...
import os import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib as mpl from scipy.stats import gaussian_kde, zscore from scipy.signal import savgol_filter from matplotlib.colors import Normalize from matplotlib.pyplot import MultipleLocator import matplotlib.font_manager as fm from matplotlib.colors import LinearSegmentedColormap # ========== 全局设置 ========== FONT_SIZE = 16 # 字体大小 FONT_NAME = 'SimSun' # 使用支持中文的字体 MARKER_SIZE = 50 # 标记大小 # ========== 列名配置（根据实际数据列名修改）========== COLUMN_NAMES = { 'rise_time': '上升时间', # 上升时间列名 'amplitude': '幅值', # 幅度列名 'ring_count': '计数', # 振铃计数列名 'duration': '持续时间' # 持续时间列名 } # 设置支持中文的字体 plt.rcParams['font.sans-serif'] = ['SimSun', 'Microsoft YaHei', 'SimHei', 'KaiTi', 'Arial Unicode MS'] plt.rcParams['axes.unicode_minus'] = False # 正确显示负号 # 设置字体和样式 config = { "font.family": 'sans-serif', "font.size": FONT_SIZE, "mathtext.fontset": 'stix' } plt.rcParams.update(config) def calculate_RA_AF(rise_time_us, amplitude_uV, ring_count, duration_us): """ 计算 RA 和 AF 的值，处理无效值 """ # 确保数据类型正确 amplitude_uV = amplitude_uV.astype(float) duration_us = duration_us.astype(float) # 创建掩码过滤无效值 valid_mask = (amplitude_uV > 0) & (duration_us > 0) # 计算 RA 和 AF RA_values = np.zeros_like(rise_time_us, dtype=float) * np.nan AF_values = np.zeros_like(rise_time_us, dtype=float) * np.nan RA_values[valid_mask] = rise_time_us[valid_mask] / amplitude_uV[valid_mask] AF_values[valid_mask] = (ring_count[valid_mask] / duration_us[valid_mask]) * 1000 return RA_values, AF_values, valid_mask def remove_outliers_zscore(df, threshold=30): """ 使用Z-score方法剔除异常值，避免精度损失问题 """ # 只对数值列计算Z-score numeric_cols = df.select_dtypes(include=[np.number]).columns # 处理常数列问题 z_scores = np.zeros_like(df[numeric_cols], dtype=float) for i, col in enumerate(numeric_cols): col_data = df[col].values if np.std(col_data) > 1e-10: # 避免几乎相同的列 z_scores[:, i] = np.abs(zscore(col_data)) else: z_scores[:, i] = 0 # 常数列Z-score设为0 # 标记异常点 outliers = (z_scores >= threshold).any(axis=1) print(f"剔除的异常点数量: {np.sum(outliers)}") if np.sum(outliers) > 0: print("部分异常点示例:") print(df[outliers].head()) return df[~outliers] def remove_outliers_IQR(df, factor=1.5): """ 使用IQR方法剔除异常值 """ # 只对数值列计算IQR numeric_cols = df.select_dtypes(include=[np.number]).columns # 初始化异常值掩码 outliers = np.zeros(len(df), dtype=bool) for col in numeric_cols: col_data = df[col].values Q1 = np.percentile(col_data, 25) Q3 = np.percentile(col_data, 75) IQR = Q3 - Q1 if IQR > 0: # 避免IQR为0的情况 lower_bound = Q1 - factor * IQR upper_bound = Q3 + factor * IQR col_outliers = (col_data < lower_bound) | (col_data > upper_bound) outliers |= col_outliers print(f"剔除的异常点数量: {np.sum(outliers)}") if np.sum(outliers) > 0: print("部分异常点示例:") print(df[outliers].head()) return df[~outliers] def find_knee_point(x_values, y_values): """ 使用kneedle算法计算曲线的拐点 基于曲率最大值确定拐点位置 """ # 平滑曲线以减少噪声影响 if len(y_values) > 5: window_length = min(11, len(y_values)) if window_length % 2 == 0: # 确保窗口长度为奇数 window_length += 1 smoothed_y = savgol_filter(y_values, window_length=window_length, polyorder=2) else: smoothed_y = y_values # 计算一阶导数（斜率） dx = x_values[1] - x_values[0] first_deriv = np.gradient(smoothed_y, dx) # 计算二阶导数（曲率） second_deriv = np.gradient(first_deriv, dx) # 计算曲率 (K = |f''| / (1 + (f')^2)^(3/2)) curvature = np.abs(second_deriv) / (1 + first_deriv ** 2) ** (3 / 2) # 找到最大曲率点（拐点） # 排除前5%和后5%的点，避免边界效应 n = len(x_values) start_idx = max(0, int(0.05 * n)) end_idx = min(n, int(0.95 * n)) valid_indices = np.arange(start_idx, end_idx) if len(valid_indices) == 0: # 如果没有有效点，使用整个范围 valid_indices = np.arange(n) max_curv_idx = valid_indices[np.argmax(curvature[valid_indices])] knee_x = x_values[max_curv_idx] knee_y = smoothed_y[max_curv_idx] return knee_x, knee_y, curvature, smoothed_y def find_b0_for_curve(b_values, percentages): """ 为单个曲线找到临界截距b0 确保 b0 > 0，并正确反映曲线特征： - b < b0 时：剪切裂纹占比随b增大而迅速增加 - b > b0 时：剪切裂纹占比随b增大增速放缓或水平变化 参数: b_values: 截距b值数组 percentages: 剪切裂纹占比数组 返回: b0: 临界截距值 (保证 b0 > 0) """ # 平滑曲线以减少噪声影响 if len(percentages) > 5: window_length = min(11, len(percentages)) if window_length % 2 == 0: # 确保窗口长度为奇数 window_length += 1 smoothed_p = savgol_filter(percentages, window_length=window_length, polyorder=2) else: smoothed_p = percentages # 计算一阶导数（变化率） dx = b_values[1] - b_values[0] first_deriv = np.gradient(smoothed_p, dx) # 计算二阶导数（变化率的变化） second_deriv = np.gradient(first_deriv, dx) # 只考虑 b > 0 的区域 positive_b_mask = b_values > 0 if not np.any(positive_b_mask): # 如果没有b>0的点，返回最小正值 return np.max(b_values) * 0.5 # 返回中间值 # 找到增速开始放缓的点（二阶导数由正变负的点） # 在 b>0 区域内寻找二阶导数小于0的点 negative_second_deriv = second_deriv < 0 # 确保只考虑b>0的区域 valid_mask = positive_b_mask & negative_second_deriv if not np.any(valid_mask): # 如果找不到增速放缓的点，使用b>0区域的最大值 return np.max(b_values[positive_b_mask]) # 找到第一个增速放缓的点（在b>0区域中第一个二阶导数小于0的点） slowdown_indices = np.where(valid_mask)[0] if len(slowdown_indices) == 0: return np.max(b_values[positive_b_mask]) first_slowdown_idx = slowdown_indices[0] b0 = b_values[first_slowdown_idx] # 确保 b0 > 0 if b0 <= 0: # 如果计算出的b0小于等于0，使用b>0区域的最小值 b0 = np.min(b_values[positive_b_mask]) return b0 def plot_all_k_b_curves_with_b0(k_values, b_values, all_percentages, knee_k, output_path, font2): """ 绘制所有斜率条件下剪切裂纹占比随截距b值的变化曲线，并标记临界截距b0 参数: k_values: 斜率k值列表 b_values: 截距b值列表 all_percentages: 字典，键为k值，值为对应b值的百分比列表 knee_k: 拐点k值 output_path: 输出文件路径 font2: 字体设置 返回: avg_b0: 所有临界截距b0的平均值 b0_values: 每个k值对应的临界截距b0列表 """ # 创建自定义颜色映射 colors = plt.cm.viridis(np.linspace(0, 1, len(k_values))) fig, ax = plt.subplots(figsize=(14, 10), dpi=100) # 存储每个k值对应的临界截距b0 b0_values = [] # 绘制每条曲线 for i, k in enumerate(k_values): # 特别标出拐点k值对应的曲线 if abs(k - knee_k) < 1e-5: # 浮点数比较容差 linewidth = 3.5 color = 'red' label = f'拐点 k={knee_k:.2f}' zorder = 10 # 确保在最上层 else: linewidth = 1.0 color = colors[i] label = f'k={k:.0f}' zorder = 1 # 获取当前k值对应的百分比曲线 percentages = all_percentages[k] # 计算当前曲线的临界截距b0 (确保b0>0) b0 = find_b0_for_curve(b_values, percentages) b0_values.append(b0) # 绘制曲线 line, = ax.plot(b_values, percentages, linewidth=linewidth, color=color, label=label, zorder=zorder) # 标记临界截距b0 (确保只标记b0>0的点) if b0 > 0: # 找到最接近b0的索引 idx = np.abs(b_values - b0).argmin() ax.scatter(b0, percentages[idx], s=100, marker='*', color=color, edgecolor='black', zorder=zorder + 1) # 计算所有临界截距b0的平均值 (只考虑b0>0的值) positive_b0_values = [b0 for b0 in b0_values if b0 > 0] if positive_b0_values: avg_b0 = np.mean(positive_b0_values) else: avg_b0 = np.max(b_values) * 0.5 # 默认值 # 在图中标记平均临界截距b0 ax.axvline(x=avg_b0, color='purple', linestyle='--', linewidth=2, alpha=0.8) ax.text(avg_b0 + 5, 10, f'平均临界截距 $b_0$ = {avg_b0:.2f}', fontsize=16, fontname=FONT_NAME, color='purple') # 添加区域标注 ax.text(avg_b0 * 0.3, 80, '区域 I: b < b₀\n剪切裂纹占比迅速增加', fontsize=14, fontname=FONT_NAME, bbox=dict(facecolor='white', alpha=0.8)) ax.text(avg_b0 * 1.5, 80, '区域 II: b > b₀\n剪切裂纹占比增速放缓', fontsize=14, fontname=FONT_NAME, bbox=dict(facecolor='white', alpha=0.8)) # 设置坐标轴 ax.set_xlabel('截距 b', fontdict=font2) ax.set_ylabel('剪切裂纹占比 (%)', fontdict=font2) ax.set_title('不同斜率条件下剪切裂纹占比随截距b的变化及临界截距$b_0$', fontdict=font2) # 添加网格 ax.grid(True, linestyle='--', alpha=0.6) # 添加图例（只显示部分k值） # 选择代表性的k值显示图例，避免过于拥挤 num_legend_items = min(15, len(k_values)) step = max(1, len(k_values) // num_legend_items) # 获取图例句柄和标签 handles, labels = ax.get_legend_handles_labels() # 创建新的图例项列表，包含拐点k值和部分其他k值 new_handles = [] new_labels = [] # 确保拐点k值在图例中 for handle, label in zip(handles, labels): if '拐点' in label: new_handles.append(handle) new_labels.append(label) break # 添加部分其他k值 for i in range(0, len(k_values), step): if i < len(handles) and '拐点' not in labels[i]: new_handles.append(handles[i]) new_labels.append(labels[i]) # 添加临界截距标记的图例 b0_marker = plt.Line2D([0], [0], marker='*', color='w', markerfacecolor='black', markersize=15, label='临界截距 $b_0$') new_handles.append(b0_marker) new_labels.append('临界截距 $b_0$') # 添加图例 ax.legend(new_handles, new_labels, loc='best', fontsize=12, title='图例说明', title_fontsize=14) # 保存图像 plt.savefig(output_path, dpi=800, bbox_inches='tight', pad_inches=0) plt.close() print(f"所有斜率b值曲线图（含临界截距$b_0$）已保存至: {output_path}") return avg_b0, b0_values def analyze_ra_af(input_file, output_dir, remove_method='zscore', k_range=(0, 50), k_step=0.1, b_analysis_k_range=(0, 300, 10), # 新增参数：用于b值分析的k范围 b_analysis_b_range=(-600, 600, 1), # 新增参数：b值范围 density_cmap='Spectral_r'): """ 执行完整的RA-AF分析流程： 1. 读取数据并计算RA和AF 2. 剔除异常值 3. 绘制RA-AF散点密度图（添加AF = knee_k + b0分界线） 4. 计算剪切裂纹占比随k值变化的曲线 5. 使用kneedle算法找到拐点k 6. 分析不同斜率条件下剪切裂纹占比随截距b值的变化 7. 对每个斜率计算临界截距b0 8. 输出包含原始数据、RA、AF、关键k值和分类结果的数据表格 9. 分别计算剪切和张拉裂纹比例 参数: input_file: 输入数据文件路径 output_dir: 输出目录路径 remove_method: 异常值剔除方法 ('zscore' 或 'IQR') k_range: k值的范围 (start, end) k_step: k值的步长 b_analysis_k_range: 用于b值分析的k范围 (start, stop, step) b_analysis_b_range: 用于b值分析的b范围 (start, stop, step) density_cmap: 密度图的颜色映射 """ # 确保输出目录存在 os.makedirs(output_dir, exist_ok=True) # 从文件名提取基本名称用于输出 base_name = os.path.splitext(os.path.basename(input_file))[0] # 1. 读取数据 df = pd.read_excel(input_file) # 使用列名配置提取数据 try: rise_time = df[COLUMN_NAMES['rise_time']].values amplitude = df[COLUMN_NAMES['amplitude']].values ring_count = df[COLUMN_NAMES['ring_count']].values duration = df[COLUMN_NAMES['duration']].values except KeyError as e: print(f"错误：数据文件中找不到配置的列名 - {e}") print("请检查COLUMN_NAMES配置与实际数据列名是否匹配") print("数据文件中的列名：", df.columns.tolist()) return # 2. 计算RA和AF RA, AF, valid_mask = calculate_RA_AF(rise_time, amplitude, ring_count, duration) # 创建包含原始数据和RA/AF的DataFrame ra_af_df = df.copy() ra_af_df['RA'] = RA ra_af_df['AF'] = AF # 只保留有效数据点 ra_af_df = ra_af_df[valid_mask] # 3. 剔除异常值 if remove_method == 'zscore': cleaned_df = remove_outliers_zscore(ra_af_df) elif remove_method == 'IQR': cleaned_df = remove_outliers_IQR(ra_af_df) else: print(f"警告: 未知的异常值剔除方法 '{remove_method}', 使用默认的zscore方法") cleaned_df = remove_outliers_zscore(ra_af_df) # 确保我们处理的是DataFrame的副本，避免SettingWithCopyWarning cleaned_df = cleaned_df.copy() # 提取处理后的数据用于绘图 x = cleaned_df['RA'].values y = cleaned_df['AF'].values # 4. 绘制散点密度图 xy = np.vstack([x, y]) z_density = gaussian_kde(xy)(xy) # 排序点以便更好地可视化 idx = z_density.argsort() x_sorted, y_sorted, z_density_sorted = x[idx], y[idx], z_density[idx] fig, ax = plt.subplots(figsize=(12, 9), dpi=100) scatter = ax.scatter(x_sorted, y_sorted, marker='o', c=z_density_sorted, edgecolors='None', s=15, label='LST', cmap=density_cmap) # 添加颜色条 cbar = plt.colorbar(scatter, shrink=1, orientation='vertical', extend='both', pad=0.015, aspect=30, label='点密度') cbar.ax.set_ylabel('点密度', fontsize=26, fontname=FONT_NAME) cbar.ax.tick_params(labelsize=26) # 设置坐标刻度 plt.tick_params(labelsize=26) labels = ax.get_xticklabels() + ax.get_yticklabels() [label.set_fontname(FONT_NAME) for label in labels] # 使用支持中文的字体 plt.rcParams['font.sans-serif'] = ['SimSun', 'Microsoft YaHei', 'SimHei', 'KaiTi', 'Arial Unicode MS'] font2 = {'family': FONT_NAME, 'weight': 'normal', 'size': 26} plt.ylabel("AF(kHz)", fontdict=font2) plt.xlabel("RA(ms/mv)", fontdict=font2) # 设置坐标轴范围和刻度 plt.axis([0, 40, 0, 600]) ax.xaxis.set_major_locator(MultipleLocator(20)) ax.yaxis.set_major_locator(MultipleLocator(50)) # 添加文本标签 plt.text(x=2, y=550, s=base_name, fontdict=font2) # 保存散点图 scatter_plot_path = os.path.join(output_dir, f'{base_name}_scatter.tiff') plt.savefig(scatter_plot_path, dpi=800, bbox_inches='tight', pad_inches=0) plt.close() print(f"散点图已保存至: {scatter_plot_path}") # 5. 计算剪切裂纹占比随k值的变化 k_values = np.arange(k_range[0], k_range[1] + k_step, k_step) percentages = [] for k in k_values: y_line = k * x mask = y < y_line percentage = np.sum(mask) / len(y) * 100 percentages.append(percentage) # 保存k值和百分比 output_df = pd.DataFrame({'k': k_values, 'Percentage': percentages}) csv_filename = os.path.join(output_dir, f'{base_name}_percentage_values.csv') output_df.to_csv(csv_filename, index=False) print(f"k值百分比数据已保存至: {csv_filename}") # 6. 使用kneedle算法找到拐点k knee_k, knee_percentage, curvature, smoothed_p = find_knee_point(k_values, percentages) # 绘制曲线图 fig2, (ax2, ax3) = plt.subplots(2, 1, figsize=(12, 12), dpi=100, sharex=True) # 使用支持中文的字体 plt.rcParams['font.sans-serif'] = ['SimSun', 'Microsoft YaHei', 'SimHei', 'KaiTi', 'Arial Unicode MS'] # 绘制原始百分比曲线 ax2.plot(k_values, percentages, 'b-', linewidth=1.5, alpha=0.6, label='原始数据') ax2.plot(k_values, smoothed_p, 'g--', linewidth=1.5, alpha=0.7, label='平滑后数据') # 标记拐点 ax2.axvline(x=knee_k, color='r', linestyle='--', alpha=0.7) ax2.plot(knee_k, knee_percentage, 'ro', markersize=8) ax2.annotate(f'拐点 k = {knee_k:.2f}', (knee_k + 0.5, knee_percentage + 5), fontsize=16, fontname=FONT_NAME) # 添加拐点说明 explanation_text = f"当斜率 k > {knee_k:.2f} 时，\n裂纹分类结果基本稳定" ax2.text(0.65, 0.25, explanation_text, transform=ax2.transAxes, fontsize=16, bbox=dict(facecolor='white', alpha=0.8), fontname=FONT_NAME) # 设置坐标轴 ax2.set_ylabel('剪切裂纹占比 (%)', fontdict=font2) ax2.set_title(f'剪切裂纹占比随分界线斜率 k 的变化 ({base_name})', fontdict=font2) # 设置坐标范围 ax2.set_xlim(k_range) ax2.set_ylim(0, 100) # 设置网格和图例 ax2.grid(True, linestyle='--', alpha=0.6) ax2.legend(loc='lower right', fontsize=16, prop={'family': FONT_NAME}) # 绘制曲率图 ax3.plot(k_values, curvature, 'm-', linewidth=2, label='曲率') ax3.axvline(x=knee_k, color='r', linestyle='--', alpha=0.7) ax3.plot(knee_k, curvature[np.where(k_values == knee_k)[0][0]], 'ro', markersize=8) # 设置曲率图 ax3.set_xlabel('分界线斜率 k', fontdict=font2) ax3.set_ylabel('曲率', fontdict=font2) ax3.set_title('曲率变化图', fontdict=font2) ax3.grid(True, linestyle='--', alpha=0.6) ax3.legend(loc='upper right', fontsize=16, prop={'family': FONT_NAME}) # 设置字体 for ax in [ax2, ax3]: labels = ax.get_xticklabels() + ax.get_yticklabels() [label.set_fontname(FONT_NAME) for label in labels] ax.tick_params(labelsize=20) # 保存百分比图 percentage_plot_path = os.path.join(output_dir, f'{base_name}_percentage_curve.tiff') plt.savefig(percentage_plot_path, dpi=800, bbox_inches='tight', pad_inches=0) plt.close(fig2) print(f"k值百分比曲线图已保存至: {percentage_plot_path}") # 7. 分析不同斜率条件下剪切裂纹占比随截距b值的变化 print("\n开始分析不同斜率条件下剪切裂纹占比随截距b值的变化...") # 准备用于b值分析的k值范围 k_start, k_stop, k_step_b = b_analysis_k_range k_values_b = np.arange(k_start, k_stop + k_step_b, k_step_b) # 准备b值范围 b_start, b_stop, b_step = b_analysis_b_range b_values = np.arange(b_start, b_stop + b_step, b_step) # 存储不同k值的临界截距b0 b0_values = [] k_b0_points = [] # 存储占比数据的字典 all_percentages = {} # 创建子目录存储b值分析结果 b_output_dir = os.path.join(output_dir, 'b_analysis') os.makedirs(b_output_dir, exist_ok=True) # 计算每个k值对应的b值占比 for i, k in enumerate(k_values_b): print(f"计算斜率 k={k:.2f} 的b值占比 ({i + 1}/{len(k_values_b)})") # 存储当前k值下不同b值的占比 percentages_b = [] for b in b_values: y_line = k * x + b mask = y < y_line percentage = np.sum(mask) / len(y) * 100 percentages_b.append(percentage) # 保存当前k值的占比数据 all_percentages[k] = percentages_b # 找到临界截距b0 b0 = find_b0_for_curve(b_values, percentages_b) # 存储临界截距b0 b0_values.append(b0) k_b0_points.append(k) # 保存数据 df_b = pd.DataFrame({ 'b': b_values, 'percentage': percentages_b, }) df_b.to_excel(os.path.join(b_output_dir, f'b_analysis_k_{k:.2f}.xlsx'), index=False) # 绘制所有斜率条件下剪切裂纹占比随截距b的变化曲线（含临界截距b0） all_k_b_path = os.path.join(output_dir, f'{base_name}_all_k_b_curves_with_b0.tiff') avg_b0, b0_values = plot_all_k_b_curves_with_b0( k_values_b, b_values, all_percentages, knee_k, all_k_b_path, font2 ) # 保存临界截距b0数据 b0_df = pd.DataFrame({ 'k': k_b0_points, 'b0': b0_values }) b0_csv = os.path.join(output_dir, f'{base_name}_critical_b0_values.csv') b0_df.to_csv(b0_csv, index=False) print(f"临界截距b0数据已保存至: {b0_csv}") # 单独绘制拐点k值时占比随b的变化 fig_b, ax_b = plt.subplots(figsize=(12, 8), dpi=100) # 获取拐点k值对应的数据 # 找到最接近拐点k值的实际k值 k_closest_idx = np.abs(k_values_b - knee_k).argmin() k_closest = k_values_b[k_closest_idx] percentages_b_closest = all_percentages[k_closest] # 找到临界截距b0 b0_closest = find_b0_for_curve(b_values, percentages_b_closest) # 绘制曲线 ax_b.plot(b_values, percentages_b_closest, 'b-', linewidth=2.5, label=f'k={k_closest:.2f} (拐点k值)') # 标记临界截距b0 ax_b.axvline(x=b0_closest, color='r', linestyle='--', alpha=0.7) ax_b.plot(b0_closest, percentages_b_closest[np.abs(b_values - b0_closest).argmin()], 'ro', markersize=8) ax_b.annotate(f'临界截距 $b_0$ = {b0_closest:.2f}', (b0_closest + 5, percentages_b_closest[np.abs(b_values - b0_closest).argmin()] + 5), fontsize=16, fontname=FONT_NAME) # 添加说明 ax_b.text(0.05, 0.95, f'斜率 k = {k_closest:.2f} (拐点k值)', transform=ax_b.transAxes, fontsize=14, bbox=dict(facecolor='white', alpha=0.8)) # 添加区域标注 ax_b.text(b0_closest * 0.3, 80, '区域 I: b < b₀\n剪切裂纹占比迅速增加', fontsize=14, fontname=FONT_NAME, bbox=dict(facecolor='white', alpha=0.8)) ax_b.text(b0_closest * 1.5, 80, '区域 II: b > b₀\n剪切裂纹占比增速放缓', fontsize=14, fontname=FONT_NAME, bbox=dict(facecolor='white', alpha=0.8)) # 设置坐标轴 ax_b.set_xlabel('截距 b', fontdict=font2) ax_b.set_ylabel('剪切裂纹占比 (%)', fontdict=font2) ax_b.set_title(f'拐点k值(k={k_closest:.2f})时剪切裂纹占比随截距b的变化', fontdict=font2) # 添加网格和图例 ax_b.grid(True, linestyle='--', alpha=0.6) ax_b.legend(loc='best', fontsize=14) # 保存图像 target_k_b_path = os.path.join(output_dir, f'{base_name}_knee_k_b_curve.tiff') plt.savefig(target_k_b_path, dpi=800, bbox_inches='tight', pad_inches=0) plt.close() print(f"拐点k值对应的b曲线图已保存至: {target_k_b_path}") # 8. 在散点图中添加AF = knee_k + b0分界线并重新保存 print("\n在散点图中添加AF = knee_k + b0分界线...") fig, ax = plt.subplots(figsize=(12, 9), dpi=100) # 重新绘制散点图 scatter = ax.scatter(x_sorted, y_sorted, marker='o', c=z_density_sorted, edgecolors='None', s=15, label='LST', cmap=density_cmap) # 添加颜色条 cbar = plt.colorbar(scatter, shrink=1, orientation='vertical', extend='both', pad=0.015, aspect=30, label='点密度') cbar.ax.set_ylabel('点密度', fontsize=26, fontname=FONT_NAME) cbar.ax.tick_params(labelsize=26) # 设置坐标轴 plt.tick_params(labelsize=26) labels = ax.get_xticklabels() + ax.get_yticklabels() [label.set_fontname(FONT_NAME) for label in labels] plt.ylabel("AF(kHz)", fontdict=font2) plt.xlabel("RA(ms/mv)", fontdict=font2) # 设置坐标轴范围和刻度 plt.axis([0, 40, 0, 600]) ax.xaxis.set_major_locator(MultipleLocator(20)) ax.yaxis.set_major_locator(MultipleLocator(50)) # 添加文本标签 plt.text(x=2, y=550, s=base_name, fontdict=font2) # 添加分界线 x_line = np.linspace(0, 40, 100) # 1. 拐点k值分界线: AF = knee_k * RA y_knee = knee_k * x_line ax.plot(x_line, y_knee, 'r-', linewidth=2.5, label=f'分界线: AF = {knee_k:.2f} * RA') # 2. 拐点k值 + 平均临界截距b0分界线: AF = knee_k * RA + avg_b0 y_b0 = knee_k * x_line + avg_b0 ax.plot(x_line, y_b0, 'b--', linewidth=2.5, label=f'分界线: AF = {knee_k:.2f} * RA + {avg_b0:.2f}') # 添加图例 ax.legend(loc='upper right', fontsize=14) # 保存带分界线的散点图 scatter_with_lines_path = os.path.join(output_dir, f'{base_name}_scatter_with_lines.tiff') plt.savefig(scatter_with_lines_path, dpi=800, bbox_inches='tight', pad_inches=0) plt.close() print(f"带分界线的散点图已保存至: {scatter_with_lines_path}") # 9. 分别计算剪切和张拉裂纹比例 print("\n计算剪切和张拉裂纹比例...") # 使用拐点k值分界线分类 mask_shear = y < knee_k * x mask_tensile = ~mask_shear # 计算比例 shear_percentage = np.sum(mask_shear) / len(y) * 100 tensile_percentage = 100 - shear_percentage # 使用拐点k值 + 平均临界截距b0分界线分类 mask_shear_b0 = y < (knee_k * x + avg_b0) mask_tensile_b0 = ~mask_shear_b0 # 计算比例 shear_percentage_b0 = np.sum(mask_shear_b0) / len(y) * 100 tensile_percentage_b0 = 100 - shear_percentage_b0 print(f"使用拐点k值分界线:") print(f" 剪切裂纹比例: {shear_percentage:.2f}%") print(f" 张拉裂纹比例: {tensile_percentage:.2f}%") print(f"\n使用拐点k值 + 平均临界截距b0分界线:") print(f" 剪切裂纹比例: {shear_percentage_b0:.2f}%") print(f" 张拉裂纹比例: {tensile_percentage_b0:.2f}%") # 10. 创建包含所有计算结果的数据表格 # 添加密度值（原始顺序） cleaned_df['Density'] = z_density # 使用拐点k值进行分类 critical_k = knee_k # 添加分类结果（剪切裂纹或拉伸裂纹） cleaned_df['Crack_Type'] = np.where( cleaned_df['AF'] < critical_k * cleaned_df['RA'], '剪切裂纹', '拉伸裂纹' ) # 添加关键k值列 cleaned_df['Critical_k'] = critical_k # 重命名列以增加可读性 column_mapping = { COLUMN_NAMES['rise_time']: '上升时间(us)', COLUMN_NAMES['amplitude']: '幅值(uV)', COLUMN_NAMES['ring_count']: '振铃计数', COLUMN_NAMES['duration']: '持续时间(us)', 'RA': 'RA(ms/mV)', 'AF': 'AF(kHz)', 'Density': '点密度', 'Crack_Type': '裂纹类型', 'Critical_k': '关键k值' } cleaned_df = cleaned_df.rename(columns=column_mapping) # 选择要保存的列 output_columns = [ '上升时间(us)', '幅值(uV)', '振铃计数', '持续时间(us)', 'RA(ms/mV)', 'AF(kHz)', '点密度', '裂纹类型', '关键k值' ] result_df = cleaned_df[output_columns] # 保存结果表格 result_table_path = os.path.join(output_dir, f'{base_name}_results.xlsx') with pd.ExcelWriter(result_table_path) as writer: # 保存完整结果 result_df.to_excel(writer, sheet_name='计算结果', index=False) # 保存关键参数汇总 shear_count = result_df[result_df['裂纹类型'] == '剪切裂纹'].shape[0] tensile_count = result_df[result_df['裂纹类型'] == '拉伸裂纹'].shape[0] total_count = len(result_df) summary_data = { '参数': [ '关键k值', '平均临界截距b0', '使用k值分界线的剪切裂纹比例(%)', '使用k值分界线的张拉裂纹比例(%)', '使用k+b0分界线的剪切裂纹比例(%)', '使用k+b0分界线的张拉裂纹比例(%)', '总数据点数', '剪切裂纹点数', '拉伸裂纹点数' ], '值': [ critical_k, avg_b0, shear_percentage, tensile_percentage, shear_percentage_b0, tensile_percentage_b0, total_count, shear_count, tensile_count ] } summary_df = pd.DataFrame(summary_data) summary_df.to_excel(writer, sheet_name='参数汇总', index=False) # 添加临界截距b0数据到汇总表 with pd.ExcelWriter(result_table_path, mode='a', if_sheet_exists='replace') as writer: b0_df.to_excel(writer, sheet_name='临界截距b0数据', index=False) # 11. 打印结果摘要 print("\n" + "=" * 50) print(f"分析完成: {base_name}") print("=" * 50) print(f"原始数据点数量: {len(df)}") print(f"有效数据点数量(剔除NaN后): {len(ra_af_df)}") print(f"剔除异常值后数据点数量: {len(cleaned_df)}") print(f"拐点k值: {knee_k:.2f} (剪切裂纹占比 = {knee_percentage:.2f}%)") print(f"当斜率 k > {knee_k:.2f} 时，裂纹分类结果基本稳定") print(f"平均临界截距 b0: {avg_b0:.2f}") print(f"临界截距b0值范围: {min(b0_values):.2f} 到 {max(b0_values):.2f}") print(f"使用k值分界线的剪切裂纹比例: {shear_percentage:.2f}%") print(f"使用k值分界线的张拉裂纹比例: {tensile_percentage:.2f}%") print(f"使用k+b0分界线的剪切裂纹比例: {shear_percentage_b0:.2f}%") print(f"使用k+b0分界线的张拉裂纹比例: {tensile_percentage_b0:.2f}%") print(f"散点图保存至: {scatter_plot_path}") print(f"带分界线的散点图保存至: {scatter_with_lines_path}") print(f"k值百分比曲线图保存至: {percentage_plot_path}") print(f"k值百分比数据保存至: {csv_filename}") print(f"所有斜率b值曲线图（含临界截距b0）保存至: {all_k_b_path}") print(f"临界截距b0数据保存至: {b0_csv}") print(f"拐点k值对应的b曲线图保存至: {target_k_b_path}") print(f"结果表格保存至: {result_table_path}") # 返回关键结果 return { 'knee_k': knee_k, 'knee_percentage': knee_percentage, 'critical_k': critical_k, 'avg_b0': avg_b0, 'b0_values': b0_values, 'shear_percentage': shear_percentage, 'tensile_percentage': tensile_percentage, 'shear_percentage_b0': shear_percentage_b0, 'tensile_percentage_b0': tensile_percentage_b0, 'cleaned_data': cleaned_df, 'scatter_path': scatter_plot_path, 'scatter_with_lines_path': scatter_with_lines_path, 'curve_path': percentage_plot_path, 'all_k_b_path': all_k_b_path, 'target_k_b_path': target_k_b_path, 'table_path': result_table_path } # ========== 使用示例 ========== if __name__ == "__main__": # 配置参数 input_file = r'J:\321\Kneedle值\AE-S-2.xlsx' output_dir = r'J:\321\Kneedle值\results' # 执行分析 results = analyze_ra_af( input_file=input_file, output_dir=output_dir, remove_method='zscore', # 可选 'zscore' 或 'IQR' k_range=(0, 500), # k值范围 k_step=0.1, # k值步长 b_analysis_k_range=(0, 2000, 10), # 用于b值分析的k范围 (start, stop, step) b_analysis_b_range=(-600, 600, 1), # 用于b值分析的b范围 (start, stop, step) density_cmap='Spectral_r' # 密度图颜色 ) # 访问结果 print(f"拐点k值: {results['knee_k']:.2f}") print(f"平均临界截距 b0: {results['avg_b0']:.2f}") print(f"使用k值分界线的剪切裂纹比例: {results['shear_percentage']:.2f}%") print(f"使用k+b0分界线的剪切裂纹比例: {results['shear_percentage_b0']:.2f}%") print(f"结果表格路径: {results['table_path']}") # ---------- 计算拉伸/剪切占比 ---------- total_count = len(df) shear_count = np.sum(df['Type'] == 'Shear Crack') tensile_count = np.sum(df['Type'] == 'Tensile Crack') shear_ratio = shear_count / total_count * 100 tensile_ratio = tensile_count / total_count * 100 # ---------- 绘制饼状图 ---------- fig, ax = plt.subplots(figsize=(6, 6)) labels = [f'剪切裂缝 ({shear_ratio:.1f}%)', f'拉伸裂缝 ({tensile_ratio:.1f}%)'] sizes = [shear_count, tensile_count] colors = ['#ff9999', '#66b3ff'] # 可以改颜色 explode = (0.05, 0.05) # 两块都稍微突出 wedges, texts, autotexts = ax.pie( sizes, labels=labels, colors=colors, explode=explode, autopct='%1.1f%%', textprops={'fontsize': FONT_SIZE, 'fontname': FONT_NAME} ) for autotext in autotexts: autotext.set_fontsize(FONT_SIZE) autotext.set_fontname(FONT_NAME) ax.set_title('拉伸与剪切裂缝占比', fontsize=FONT_SIZE + 2, fontname=FONT_NAME) plt.tight_layout() plt.show()
08-09
plt.title('Shear Crack Ratio vs. Slope k') plt.legend() # 8. 在散点图上绘制最终分界线 plt.figure(figsize=(10, 6)) sc = plt.scatter(ra, af, c=z, s=10, cmap='Spectral_r') plt.colorbar(sc, label='...
--------------------------------------------------------------------------- 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&sup2; 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&sup2;") logger.info(f"平均日灌溉需求: {avg_demand_per_m2:.2f} L/m&sup2;") 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&sup2;)') 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
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/...
import pandas as pd import numpy as np import matplotlib.pyplot as plt from itertools import combinations import random import time import matplotlib as mpl from tqdm import tqdm # 设置中文字体支持 plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'WenQuanYi Micro Hei'] plt.rcParams['axes.unicode_minus'] = False # ==================================================== # 任务1: 订单需求分析 # ==================================================== print("正在执行任务1: 订单需求分析...") # 加载订单数据 demand_df = pd.read_csv('OrderDemand.csv') # 计算每个区域的总需求和日均需求 region_demand = demand_df.groupby('Region')['OrderVolume'].sum().reset_index() region_demand['AvgDailyDemand'] = region_demand['OrderVolume'] / 3 region_demand_sorted = region_demand.sort_values(by='AvgDailyDemand', ascending=False) # 计算每个小时的总需求 hourly_demand = demand_df.groupby('Hour')['OrderVolume'].sum().reset_index() hourly_demand['AvgHourlyDemand'] = hourly_demand['OrderVolume'] / 3 # 识别高峰时段 peak_hours = hourly_demand.nlargest(3, 'AvgHourlyDemand')['Hour'].tolist() # 输出分析结果 print("\n=== 高需求区域Top 5 ===") print(region_demand_sorted.head(5)) print("\n=== 高峰时段 ===") print(f"高峰时段(小时): {peak_hours}") # 可视化需求分析 plt.figure(figsize=(14, 7)) bars = plt.bar(region_demand_sorted['Region'], region_demand_sorted['AvgDailyDemand'], color='skyblue') plt.title('区域日均需求分析', fontsize=16) plt.xlabel('区域', fontsize=12) plt.ylabel('日均需求', fontsize=12) plt.xticks(rotation=45) plt.grid(axis='y', linestyle='--', alpha=0.7) # 添加数据标签 for bar in bars: height = bar.get_height() plt.annotate(f'{height:.1f}', xy=(bar.get_x() + bar.get_width() / 2, height), xytext=(0, 3), textcoords="offset points", ha='center', va='bottom', fontsize=9) plt.tight_layout() plt.savefig('region_demand.png', dpi=300) plt.close() plt.figure(figsize=(14, 7)) plt.plot(hourly_demand['Hour'], hourly_demand['AvgHourlyDemand'], marker='o', linestyle='-', color='royalblue', linewidth=2.5) plt.title('小时平均需求分布', fontsize=16) plt.xlabel('小时', fontsize=12) plt.ylabel('平均需求', fontsize=12) plt.grid(True, linestyle='--', alpha=0.7) plt.xticks(range(0, 24)) # 标记高峰时段 for i, hour in enumerate(peak_hours): plt.axvline(x=hour, color='r', linestyle='--', alpha=0.8) plt.text(hour, max(hourly_demand['AvgHourlyDemand']) * (0.85 - i*0.1), f'高峰 {hour}:00', rotation=90, ha='right', va='top', color='r', fontsize=11, bbox=dict(facecolor='white', alpha=0.8, edgecolor='none')) plt.legend(['平均需求', '高峰时段'], loc='upper left') plt.tight_layout() plt.savefig('hourly_demand.png', dpi=300) plt.close() # ==================================================== # 任务2: 网络优化模型 # ==================================================== print("\n正在执行任务2: 网络优化模型...") # 加载距离矩阵 dist_matrix = pd.read_csv('DistanceMatrix.csv', index_col=0) # 分配区域到最近的DC dc_allocation = {} for region in dist_matrix.columns: closest_dc = dist_matrix[region].idxmin() dc_allocation[region] = closest_dc # 为每个DC创建服务区域列表 dc_regions = {dc: [] for dc in dist_matrix.index} for region, dc in dc_allocation.items(): dc_regions = dc_regions[dc] if region not in dc_regions: dc_regions[dc].append(region) # 生成区域坐标（模拟位置） np.random.seed(42) positions = {} for dc in dist_matrix.index: positions[dc] = (random.uniform(0, 10), random.uniform(0, 10)) for region in dc_regions[dc]: # 区域在DC附近随机分布 positions[region] = ( positions[dc][0] + random.uniform(-2, 2), positions[dc][1] + random.uniform(-2, 2) ) # 计算区域间距离 def calculate_distance(pos1, pos2): return np.sqrt((pos1[0]-pos2[0])**2 + (pos1[1]-pos2[1])**2) # 最近邻算法求解TSP def nearest_neighbor_tsp(nodes, positions): if not nodes: return [] unvisited = nodes.copy() current = unvisited.pop(0) tour = [current] while unvisited: # 找到最近邻居 nearest = min(unvisited, key=lambda x: calculate_distance(positions[current], positions[x])) tour.append(nearest) unvisited.remove(nearest) current = nearest return tour # 为每个DC规划路径 dc_paths = {} for dc, regions in dc_regions.items(): if regions: # 包含DC的节点列表 nodes = [dc] + regions path = nearest_neighbor_tsp(nodes, positions) dc_paths[dc] = path else: dc_paths[dc] = [dc] # 计算路径总距离 def calculate_path_distance(path, positions): if len(path) < 2: return 0 distance = 0 for i in range(len(path)-1): distance += calculate_distance(positions[path[i]], positions[path[i+1]]) # 添加返回起点的距离 distance += calculate_distance(positions[path[-1]], positions[path[0]]) return distance # 输出任务2结果 print("\n=== DC分配结果 ===") for dc, regions in dc_regions.items(): print(f"{dc} 服务区域: {regions}") print("\n=== 配送路径 ===") total_distance = 0 for dc, path in dc_paths.items(): path_distance = calculate_path_distance(path, positions) total_distance += path_distance path_str = ' → '.join(path) + f" → {path[0]}" if len(path) > 1 else path[0] print(f"{dc} 路径: {path_str} | 距离: {path_distance:.2f} km") print(f"总配送距离: {total_distance:.2f} km") # ==================================================== # 任务3: 多车辆路径调度模型 (修复版) # ==================================================== print("\n正在执行任务3: 多车辆路径调度模型...") # 加载车辆信息 vehicle_df = pd.read_csv('VehicleInfo.csv') vehicles = vehicle_df.to_dict('records') # 区域需求数据（来自任务1） region_demands = region_demand.set_index('Region')['AvgDailyDemand'].to_dict() # 增强版节约算法 def enhanced_savings_algorithm(depot, nodes, demands, positions, vehicles, service_time_per_demand=5): # 如果没有节点，返回空列表 if not nodes: return [] # 创建可用车辆列表的深拷贝 available_vehicles = deepcopy(vehicles) routes = [] assigned_nodes = set() # 步骤1: 创建初始路线（每个节点一条路线） for node in nodes: # 计算单点路线的需求和时间 total_demand = demands[node] travel_time = (calculate_distance(positions[depot], positions[node]) * 2) * 1.5 # 往返时间 service_time = total_demand * service_time_per_demand total_time = travel_time + service_time # 寻找合适的车辆 vehicle_found = False for i, vehicle in enumerate(available_vehicles): if vehicle['Capacity'] >= total_demand and vehicle['MaxServiceTime'] >= total_time: routes.append({ 'nodes': [node], 'demand': total_demand, 'time': total_time, 'vehicle_id': vehicle['VehicleID'] }) assigned_nodes.add(node) # 移除已使用的车辆 del available_vehicles[i] vehicle_found = True break # 如果没有找到合适的车辆，使用最小容量的车辆（即使超载） if not vehicle_found and available_vehicles: min_capacity_vehicle = min(available_vehicles, key=lambda x: x['Capacity']) routes.append({ 'nodes': [node], 'demand': total_demand, 'time': total_time, 'vehicle_id': min_capacity_vehicle['VehicleID'] }) assigned_nodes.add(node) available_vehicles.remove(min_capacity_vehicle) # 步骤2: 计算节约值矩阵 savings = [] for i, j in combinations(nodes, 2): # 如果两个点都已被分配，计算节约值 if i in assigned_nodes and j in assigned_nodes: saving = calculate_distance(positions[depot], positions[i]) + \ calculate_distance(positions[depot], positions[j]) - \ calculate_distance(positions[i], positions[j]) savings.append((i, j, saving)) # 按节约值降序排序 savings.sort(key=lambda x: x[2], reverse=True) # 步骤3: 合并路线 for i, j, saving in savings: # 找到包含i和j的路线 route_i = next((r for r in routes if i in r['nodes']), None) route_j = next((r for r in routes if j in r['nodes']), None) # 确保是两条不同的路线 if route_i is None or route_j is None or route_i == route_j: continue # 尝试不同的合并方式 merge_success = False # 方式1: i_end -> j_start if route_i['nodes'][-1] == i and route_j['nodes'][0] == j: new_nodes = route_i['nodes'] + route_j['nodes'] merge_success = True # 方式2: j_end -> i_start elif route_j['nodes'][-1] == j and route_i['nodes'][0] == i: new_nodes = route_j['nodes'] + route_i['nodes'] merge_success = True # 方式3: i_end -> j_end (需要反转j路线) elif route_i['nodes'][-1] == i and route_j['nodes'][-1] == j: new_nodes = route_i['nodes'] + list(reversed(route_j['nodes'])) merge_success = True # 方式4: i_start -> j_start (需要反转i路线) elif route_i['nodes'][0] == i and route_j['nodes'][0] == j: new_nodes = list(reversed(route_j['nodes'])) + route_i['nodes'] merge_success = True if merge_success: # 计算新路线的总需求 total_demand = sum(demands[node] for node in new_nodes) # 计算新路线的总时间 path_distance = calculate_path_distance([depot] + new_nodes + [depot], positions) travel_time = path_distance * 1.5 # 速度40km/h，转换为分钟 service_time = total_demand * service_time_per_demand total_time = travel_time + service_time # 检查车辆是否满足新路线的需求 current_vehicle_id = route_i['vehicle_id'] current_vehicle = next((v for v in vehicles if v['VehicleID'] == current_vehicle_id), None) if current_vehicle and current_vehicle['Capacity'] >= total_demand and current_vehicle['MaxServiceTime'] >= total_time: # 合并路线 new_route = { 'nodes': new_nodes, 'demand': total_demand, 'time': total_time, 'vehicle_id': current_vehicle_id } # 移除旧路线，添加新路线 routes.remove(route_i) routes.remove(route_j) routes.append(new_route) # 释放route_j的车辆 released_vehicle_id = route_j['vehicle_id'] if released_vehicle_id != current_vehicle_id: released_vehicle = next((v for v in vehicles if v['VehicleID'] == released_vehicle_id), None) if released_vehicle: available_vehicles.append(released_vehicle) # 返回格式化结果 result = [] for route in routes: result.append(( route['vehicle_id'], route['nodes'], route['demand'], route['time'] )) return result # 为每个DC规划多车辆路径 dc_vehicle_routes = {} start_time = time.time() for dc, regions in tqdm(dc_regions.items(), desc="规划车辆路径"): if regions: routes = enhanced_savings_algorithm(dc, regions, region_demands, positions, vehicles) dc_vehicle_routes[dc] = routes else: dc_vehicle_routes[dc] = [] # 输出任务3结果 print("\n=== 多车辆路径调度结果 ===") total_vehicles_used = 0 total_vehicles_available = len(vehicles) for dc, routes in dc_vehicle_routes.items(): if not routes: print(f"\n{dc} 车辆调度方案: 无可用路径") continue print(f"\n{dc} 车辆调度方案:") for vehicle_id, route_nodes, demand, time_used in routes: total_vehicles_used += 1 # 获取车辆信息 vehicle = next((v for v in vehicles if v['VehicleID'] == vehicle_id), None) if not vehicle: vehicle_capacity = "未知" vehicle_maxtime = "未知" else: vehicle_capacity = vehicle['Capacity'] vehicle_maxtime = vehicle['MaxServiceTime'] # 计算利用率 capacity_util = (demand / vehicle_capacity * 100) if vehicle_capacity > 0 else 0 time_util = (time_used / vehicle_maxtime * 100) if vehicle_maxtime > 0 else 0 # 打印结果 route_str = f"{dc} → {' → '.join(route_nodes)} → {dc}" print(f" 车辆 {vehicle_id}:") print(f" 路径: {route_str}") print(f" 需求: {demand:.2f} / {vehicle_capacity} ({capacity_util:.1f}%)") print(f" 时间: {time_used:.2f} min / {vehicle_maxtime} min ({time_util:.1f}%)") # 计算总体利用率 utilization_rate = (total_vehicles_used / total_vehicles_available * 100) if total_vehicles_available > 0 else 0 print(f"\n总使用车辆数: {total_vehicles_used} / {total_vehicles_available}") print(f"车辆利用率: {utilization_rate:.1f}%") print(f"路径规划耗时: {time.time()-start_time:.2f}秒") # ==================================================== # 保存结果 # ==================================================== # 保存任务1结果 region_demand_sorted.to_csv('region_demand_analysis.csv', index=False) hourly_demand.to_csv('hourly_demand_analysis.csv', index=False) # 保存任务2结果 dc_allocation_df = pd.DataFrame({ 'Region': list(dc_allocation.keys()), 'AssignedDC': list(dc_allocation.values()) }) dc_allocation_df.to_csv('dc_allocation.csv', index=False) # 保存任务3结果 vehicle_routes_data = [] for dc, routes in dc_vehicle_routes.items(): for vehicle_id, route_nodes, demand, time_used in routes: vehicle = next((v for v in vehicles if v['VehicleID'] == vehicle_id), None) vehicle_capacity = vehicle['Capacity'] if vehicle else 0 vehicle_routes_data.append({ 'DC': dc, 'VehicleID': vehicle_id, 'Route': ' → '.join([dc] + route_nodes + [dc]), 'TotalDemand': demand, 'VehicleCapacity': vehicle_capacity, 'CapacityUtilization(%)': (demand / vehicle_capacity * 100) if vehicle_capacity > 0 else 0, 'TotalTime(min)': time_used, 'MaxServiceTime(min)': vehicle['MaxServiceTime'] if vehicle else 0, 'TimeUtilization(%)': (time_used / vehicle['MaxServiceTime'] * 100) if vehicle else 0 }) vehicle_routes_df = pd.DataFrame(vehicle_routes_data) vehicle_routes_df.to_csv('vehicle_routes.csv', index=False) # 生成配送网络图 plt.figure(figsize=(16, 12)) # 创建颜色映射 dc_colors = {} color_palette = plt.cm.tab10.colors for idx, dc in enumerate(dist_matrix.index): dc_colors[dc] = color_palette[idx % len(color_palette)] # 绘制DC位置 for dc in dist_matrix.index: pos = positions[dc] plt.scatter(pos[0], pos[1], s=400, c=[dc_colors[dc]], marker='s', edgecolor='k', linewidth=2, zorder=10) plt.text(pos[0], pos[1], dc, fontsize=14, ha='center', va='center', fontweight='bold', zorder=11) # 绘制区域位置和路径 for dc, routes in dc_vehicle_routes.items(): color = dc_colors[dc] # 绘制该DC下的所有区域 for region in dc_regions[dc]: pos = positions[region] plt.scatter(pos[0], pos[1], s=120, c=[color], marker='o', alpha=0.8, zorder=9) plt.text(pos[0], pos[1], region, fontsize=10, ha='center', va='bottom', zorder=9) # 绘制该DC的所有车辆路径 for route in routes: _, route_nodes, _, _ = route path = [dc] + route_nodes + [dc] # 绘制路径线 for i in range(len(path)-1): start = positions[path[i]] end = positions[path[i+1]] plt.plot([start[0], end[0]], [start[1], end[1]], c=color, linestyle='-', linewidth=1.5, alpha=0.7, zorder=5) # 标记路径起点 plt.scatter(positions[dc][0], positions[dc][1], s=150, c='red', marker='*', zorder=12) # 添加图例和标题 plt.title('物流配送网络与车辆路径规划', fontsize=18) plt.xlabel('X 坐标', fontsize=14) plt.ylabel('Y 坐标', fontsize=14) plt.grid(True, linestyle='--', alpha=0.5) # 创建自定义图例 from matplotlib.lines import Line2D legend_elements = [ Line2D([0], [0], marker='s', color='w', label='配送中心(DC)', markerfacecolor='red', markersize=12), Line2D([0], [0], marker='o', color='w', label='配送区域', markerfacecolor='blue', markersize=10), Line2D([0], [0], marker='*', color='w', label='路径起点', markerfacecolor='red', markersize=12), Line2D([0], [0], color='green', lw=2, label='配送路径') ] plt.legend(handles=legend_elements, loc='upper left', fontsize=12) plt.tight_layout() plt.savefig('delivery_network.png', dpi=300) plt.close() print("\n所有结果已保存到CSV文件和图表") TypeError: list indices must be integers or slices, not str给我解决完这个问题后修正后的完整的代码
06-03
$$ \text{查询效率} = O(\log n) \quad \text{vs} \quad O(n) \quad (\text{原始线性搜索}) $$ ```python from scipy.spatial import cKDTree def build_spatial_index(positions): points = [(x, y) for (x, y) ...
closest
MountGao的博客
08-10 1515
closest 简介 该closest方法返回当前元素（或当前元素本身）的最接近的祖先，该祖先与参数中给出的选择器匹配。如果没有这样的祖先，它就会返回null。 语法 var elt = element .closest（selectors）; 示例 &lt;article&gt; &lt;div id="div-01"&gt;Here is div...
有限自动机总结
binling的专栏
09-26 6613
有限自动机A用来识别字符串，它由5部分组成： 1）alphabet，字符集 2）states，状态集合 3）init，初始状态 4）trans(s, ch)，状态转移函数 5）end 可接受state 集合 A(str) == true的意思是，A可以接受字符串str，即从初始状态init读入str所有字符之后所达到的状态属于集合end， tran(init, str) 属于 end
一个完整的语法分析、词法分析例子——Universal Pasrser
binling的专栏
11-12 4595
需求：用户用formal notation指定语法、词法，然后可以匹配相应的文本。用法类似正则表达式，只需给出formal notation，不需要为每一种格式的文本单独写匹配器。 formal notation主要是3个部分： 1）BNF 列表 table：给出上下文无关文法的产生式规则，以及所有的符号（终端和非终端） 2）起始符号 start 3）终端符号Regular Express
双向BFS搜索和A*算法
binling的专栏
12-03 4224
适合给出起点和终点，求最短路径的问题 分别从起点和终点扩展，每次选择待扩展节点少的那个方向进行扩展，扩展一层。然后下一轮。 int doubleBFS(int start, int end, vector> &G) { vector marked[2] {vector(G.size()), vector(G.size())}; //including dead nodes and to-be
爬虫中网页分析的几种技术
binling的专栏
10-20 4045
一般来说我们只抓取网页中的特定数据，比如抓取某人所有的blog，我们就只关心list 页面中文章列表那部分的链接和title 有3中可以选技术来分析网页 1）正则匹配 2）一般字符串匹配content.substring(pattern, startIndex)，一般是带有startIndex的substring，而不是每次都是从头匹配 3) 基于sax的事件 4）DOM + XPath
彻底解决数据库去中心化问题
binling的专栏
08-12 4041
1）首先垂直切分，按业务模块把一些紧密联系的表分在一起（按表切分)，形成一个抽象的vertical partition 2）水平切分。垂直切分之后的单个vertical partition（甚至其中单个表）仍然可能太大，一个server放不下，进一步按行切分，这时候需要选择一个聚合根作为记录的horizontal partiotion的标志，一般是userId，对于Forum模块，forumId
最少砝码问题（用一部分数的和/差表示区间上所有的整数）
binling的专栏
08-07 2927
问题1, 需要表示[1, N] 的所有重量，最少需要多少砝码？ 答案：需要 1，2， 3,  ... , logN，也就是 logN个。每个砝码代表二进制数的一位，N有logN个二进制位。 问题2， 需要表示[1, N] 的所有重量，手头已有一些砝码，问：怎样添加新的砝码，使得需要添加的砝码数最少？ 不变式：curMax为当前可以表示的重量的最大值，也就是[1, curMax] 都可
关于线段的算法汇总
binling的专栏
05-15 2858
基本知识： 线段树：存储线段的二叉排序树，以start作为key, 同时维护树中最大的右端点值 max class Node { int start; int end; int max; }经典应用：query any interval that intersect(overlap) with a given interval. 如果是查询所有与给定线段重合的线段，则查
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
表情包 代码片
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值