求解0.5f == 0.5为true,0.3f == 0.3为false

浮点数精度解析与计算陷阱
本文深入探讨了计算机处理浮点数时的精度问题,通过实例解释了为何0.5f、0.25f能准确判断为true,而0.3f则不能。文章详细阐述了浮点数在内存中的存储方式及其导致的精度损失,揭示了这一现象背后的数学原理。
求解0.5f == 0.5为true,0.3f == 0.3为false,0.25f == 0.25又为true,为何

精度问题。 计算机处理浮点数的时候不是你看到的0.3或者0.5而是一组浮点数存在于内存中。
就拿你举得3个数为例
0.3在计算机中浮点表示为(00111110100110011001100110011010) 因为精度原因后面被截断所以实际换算回十进制的时候约为0.299999
0.5在计算机中浮点表示为00111111000000000000000000000000
0.25在计算机中浮点表示为00111110100000000000000000000000 不存在截断所以数值相等
所以0.5f==0.5,0.25f==0.25但是0.3f!=0.3
实际上只有1除以2的倍数得到的小数才能被判断成true因为在换算成浮点数的时候不存在截断。而其他的小数都会存在精度截断的情况。

补充:
更正一下
实际上只有1除以(2的n次方)得到的小数才能被判断成true因为在换算成浮点数的时候不存在截断。而其他的小数都会存在精度截断的情况。

copy自搜搜问问http://wenwen.soso.com/z/q488392311.htm?ch=izw.d.wt
,,如有问题请反驳。。
-----------小数转二进制
可以这样:首先将一个小数如:235.725的小数部分取出,即:0.725,将其乘以进制数二进制就乘以2后得到1。45,取其整数部分1为二进制小数的第一项(十分位),在将小数部分0。45乘2得0。9,取其整数部分为二进制小数的第二位(百分位)0,在将其小数部分0。9乘2,得1。8,取其整数部分为二进制小数的第三位(千分位)1,取其小数部分0。8再乘2……以此类推,直到值为0或形成循环小数则停止。
import time import math import numpy as np import matplotlib.pyplot as plt from scipy.spatial import ConvexHull from scipy.spatial.distance import cdist from scipy.stats import qmc # 全局字体设置,解决中文和特殊符号显示问题 plt.rcParams["font.family"] = ["SimHei", "Arial Unicode MS"] plt.rcParams[&#39;axes.unicode_minus&#39;] = False # 解决负号显示问题 # 1. 目标函数向量化:支持批量解输入(核心优化点1) def zdt1(x): """ZDT1测试函数(支持2D数组输入:n行×dim列)""" x = np.atleast_2d(x) # 确保输入是2D数组 n, dim = x.shape f1 = x[:, 0] # 所有解的第0维(批量取列) g = 1 + 9 * np.mean(x[:, 1:], axis=1) # 每行(每个解)的均值 h = 1 - np.sqrt(f1 / g) f2 = g * h return np.column_stack((f1, f2)) # 输出n行×2列 def zdt2(x): """ZDT2测试函数(支持批量解输入)""" x = np.atleast_2d(x) n, dim = x.shape f1 = x[:, 0] g = 1 + 9 * np.mean(x[:, 1:], axis=1) h = 1 - (f1 / g) ** 2 f2 = g * h return np.column_stack((f1, f2)) def zdt3(x): """ZDT3测试函数(支持批量解输入)""" x = np.atleast_2d(x) n, dim = x.shape f1 = x[:, 0] g = 1 + 9 * np.mean(x[:, 1:], axis=1) h = 1 - np.sqrt(f1 / g) - (f1 / g) * np.sin(10 * np.pi * f1) return np.column_stack((f1, g * h)) def zdt4(x): """严格修正的ZDT4实现""" x = np.atleast_2d(x) n, dim = x.shape f1 = x[:, 0] # 严格遵循原始论文公式 g = 1 + 10 * (dim - 1) + np.sum( x[:, 1:] ** 2 - 10 * np.cos(4 * np.pi * x[:, 1:]), axis=1 ) h = 1 - np.sqrt(f1 / g) f2 = g * h return np.column_stack((f1, f2)) def zdt6(x): """ZDT6测试函数(支持批量解输入)""" x = np.atleast_2d(x) n, dim = x.shape # ZDT6公式:f1 = 1 - exp(-4x1)sin^6(6πx1) f1 = 1 - np.exp(-4 * x[:, 0]) * (np.sin(6 * np.pi * x[:, 0])) ** 6 # g(x) = 1 + 9*(sum(x2~xd)/(d-1))^0.25 g = 1 + 9 * (np.sum(x[:, 1:], axis=1) / (dim - 1)) ** 0.25 # h = 1 - (f1/g)^2 h = 1 - (f1 / g) ** 2 f2 = g * h return np.column_stack((f1, f2)) # --------------------------- 新增:二维DTLZ1/2/3目标函数 --------------------------- def dtlz1(x): """ 二维DTLZ1测试函数(M=2目标,决策变量n=3) 特性:线性Pareto前沿、多模态 决策变量范围:x ∈ [0,1]^3 真实PF:f1 + f2 = 0.5(f1,f2 ≥ 0) """ x = np.atleast_2d(x) # 确保输入为2D数组(n行×3列) n, dim = x.shape assert dim == 3, "DTLZ1(二维目标)需输入3个决策变量(n=3)" # 向量化计算g(x),确保形状为(n,) g = 100 * ( 2 + # 3-2+1=2(n-M+1) (x[:, 2] - 0.5) ** 2 - np.cos(20 * np.pi * (x[:, 2] - 0.5)) ) # 目标函数计算(向量化操作,避免广播错误) f1 = 0.5 * x[:, 0] * x[:, 1] * (1 + g) f2 = 0.5 * (1 - x[:, 0]) * (1 + g) return np.column_stack((f1, f2)) def dtlz2(x): """ 二维DTLZ2测试函数(M=2目标,决策变量n=3) 特性:凸形Pareto前沿、非线性目标关联 决策变量范围:x ∈ [0,1]^3 真实PF:f1² + f2² = 1(f1,f2 ≥ 0) """ x = np.atleast_2d(x) # 确保输入为2D数组(n行×3列) n, dim = x.shape assert dim == 3, "DTLZ2(二维目标)需输入3个决策变量(n=3)" # 向量化计算g(x)(已正确实现,保持并优化注释) g = np.sum((x[:, 2] - 0.5) ** 2, axis=0) # 修正:对于单个变量无需axis=1,避免降维问题 # 目标函数计算 f1 = (1 + g) * np.cos(x[:, 0] * np.pi / 2) * np.cos(x[:, 1] * np.pi / 2) f2 = (1 + g) * np.cos(x[:, 0] * np.pi / 2) * np.sin(x[:, 1] * np.pi / 2) return np.column_stack((f1, f2)) def dtlz3(x): """ 二维DTLZ3测试函数(M=2目标,决策变量n=3) 特性:凸形Pareto前沿、强多模态 决策变量范围:x ∈ [0,1]^3 真实PF:f1² + f2² = 1(f1,f2 ≥ 0) """ x = np.atleast_2d(x) # 确保输入为2D数组(n行×3列) n, dim = x.shape assert dim == 3, "DTLZ3(二维目标)需输入3个决策变量(n=3)" # 向量化计算g(x),确保形状为(n,) g = 100 * ( 2 + # 3-2+1=2(n-M+1) (x[:, 2] - 0.5) ** 2 - np.cos(20 * np.pi * (x[:, 2] - 0.5)) ) # 目标函数计算(向量化操作) f1 = (1 + g) * np.cos(x[:, 0] * np.pi / 2) * np.cos(x[:, 1] * np.pi / 2) f2 = (1 + g) * np.cos(x[:, 0] * np.pi / 2) * np.sin(x[:, 1] * np.pi / 2) return np.column_stack((f1, f2)) def generate_true_front(func, num_points=1000): """ 生成ZDT测试函数真实Pareto前沿,包含ZDT6的实现 :param zdt_func: 目标函数(zdt1/zdt2/zdt3/zdt4/zdt6) :param num_points: 采样点数(越大越精确) :return: 真实前沿的(f1, f2)数组 """ # 以下为原有其他ZDT函数的处理逻辑,保持不变 f1 = np.linspace(0, 1, num_points) f2 = [] # --------------------------- 新增:DTLZ系列前沿生成 --------------------------- if func == dtlz1: # DTLZ1真实PF:f1 + f2 = 0.5(f1∈[0,0.5],f2=0.5-f1) f1 = np.linspace(0, 0.5, num_points) f2 = 0.5 - f1 elif func == dtlz2 or func == dtlz3: # DTLZ2/3真实PF:f1² + f2² = 1(f1∈[0,1],f2=√(1-f1²)) # 注:两函数PF形状相同,仅决策空间多模态差异 theta = np.linspace(0, np.pi / 2, num_points) # 极角(0~90°) f1 = np.cos(theta) f2 = np.sin(theta) if func == zdt1: for f in f1: f2_val = 1 - np.sqrt(f) f2.append(f2_val) elif func == zdt2: for f in f1: f2_val = 1 - f ** 2 f2.append(f2_val) elif func == zdt3: # ZDT3 真实前沿:g(x)=1(x2~xd=0),逐点计算并保留有效分段 for f in f1: g = 1.0 h = 1 - np.sqrt(f / g) - (f / g) * np.sin(10 * np.pi * f) f2_val = g * h # 仅保留物理意义有效的点(过滤极端异常值) if -1.5 <= f2_val <= 1.5: f2.append(f2_val) else: f2.append(np.nan) # 无效点标记为NaN elif func == zdt4: for f in f1: g_fixed = 1.0 # 固定g(x)=1 h = 1 - np.sqrt(f / g_fixed) f2_val = g_fixed * h f2.append(f2_val) elif func == zdt6: # ZDT6真实前沿生成:x2~xd=0时,g(x)=1 for x1 in np.linspace(0, 1, num_points): f1_val = 1 - np.exp(-4 * x1) * (np.sin(6 * np.pi * x1)) ** 6 g = 1.0 # x2~xd=0时g(x)=1 h = 1 - (f1_val / g) ** 2 f2_val = g * h f2.append(f2_val) f1 = [1 - np.exp(-4 * x1) * (np.sin(6 * np.pi * x1)) ** 6 for x1 in np.linspace(0, 1, num_points)] # 转换为numpy数组并过滤NaN f2 = np.array(f2) front = np.column_stack((f1[:len(f2)], f2)) # ZDT3 额外过滤被支配点,确保严格前沿 if func == zdt3: to_keep = np.ones(len(front), dtype=bool) for i in range(len(front)): f1_i, f2_i = front[i] for j in range(len(front)): if i == j: continue f1_j, f2_j = front[j] # 若点j支配点i,则过滤i if f1_j <= f1_i and f2_j <= f2_i and (f1_j < f1_i or f2_j < f2_i): to_keep[i] = False break front = front[to_keep] return front def calculate_igd(algorithm_front, true_front): """计算反向世代距离(IGD)""" if len(algorithm_front) == 0: return float(&#39;inf&#39;) distances = cdist(true_front, algorithm_front).min(axis=1) return np.mean(distances) def initialization(pop_size, dim, lb, ub): """优化版种群初始化:支持固定维度""" lb = np.array(lb) if not isinstance(lb, np.ndarray) else lb ub = np.array(ub) if not isinstance(ub, np.ndarray) else ub lb = np.repeat(lb, dim) if lb.size == 1 else lb ub = np.repeat(ub, dim) if ub.size == 1 else ub # 分离固定维度(上下界相等)和可采样维度 fixed_dims = lb == ub sample_dims = ~fixed_dims n_sample = np.sum(sample_dims) # 初始化种群(固定维度先设为lb值) population = np.full((pop_size, dim), lb) # 仅对可采样维度做LHS采样 if n_sample > 0: sampler = qmc.LatinHypercube(d=n_sample) samples = sampler.random(n=pop_size) # 映射到可采样维度的上下界 scaled_samples = qmc.scale( samples, lb[sample_dims], ub[sample_dims] ) population[:, sample_dims] = scaled_samples return population def Domination(X, FitX): N = X.shape[0] dominated_count = np.zeros(N, dtype=int) # 生成所有解对的支配关系矩阵(向量化比较) dominate_matrix = np.all(FitX[:, np.newaxis] <= FitX, axis=2) & np.any(FitX[:, np.newaxis] < FitX, axis=2) # 统计每个解被多少个其他解支配 dominated_count = np.sum(dominate_matrix, axis=0) # 筛选非支配解 mask = dominated_count == 0 return X[mask], FitX[mask] def crowding_distance(fitness): """计算每个解的拥挤距离""" n, m = fitness.shape if n <= 2: return np.ones(n) * np.inf # 解数量过少时,距离设为无穷大 crowd_dist = np.zeros(n) for j in range(m): sorted_idx = np.argsort(fitness[:, j]) # 按第j个目标排序 sorted_vals = fitness[sorted_idx, j] # 边界解距离设为无穷大(优先保留) crowd_dist[sorted_idx[0]] = np.inf crowd_dist[sorted_idx[-1]] = np.inf # 中间解计算距离 if sorted_vals[-1] != sorted_vals[0]: # 避免除零 crowd_dist[sorted_idx[1:-1]] += (sorted_vals[2:] - sorted_vals[:-2]) / (sorted_vals[-1] - sorted_vals[0]) return crowd_dist def crowding_distance_crop(X, FitX, L): """使用拥挤距离裁剪种群到大小L""" crowd_dist = crowding_distance(FitX) sorted_indices = np.argsort(-crowd_dist) # 按距离降序 return X[sorted_indices[:L], :], FitX[sorted_indices[:L], :] def UpdateArchive(NumObj, archive_size, A, FitA, X, FitX): """优化的存档更新函数""" A = np.atleast_2d(A) if len(A) > 0 else np.empty((0, X.shape[1])) FitA = np.atleast_2d(FitA) if len(FitA) > 0 else np.empty((0, NumObj)) X = np.atleast_2d(X) FitX = np.atleast_2d(FitX) if len(X) == 0: return A, FitA # 合并存档解和新解 combined_pop = np.vstack([A, X]) combined_fit = np.vstack([FitA, FitX]) total = len(combined_pop) # 向量化计算支配关系矩阵 dominate_matrix = ( np.all(combined_fit[:, np.newaxis] <= combined_fit, axis=2) & np.any(combined_fit[:, np.newaxis] < combined_fit, axis=2) ) # 计算每个解的被支配计数 dominated_count = np.sum(dominate_matrix, axis=0) # 筛选所有非支配解 non_dominated_mask = (dominated_count == 0) non_dominated_pop = combined_pop[non_dominated_mask] non_dominated_fit = combined_fit[non_dominated_mask] # 按拥挤距离裁剪到存档大小 if len(non_dominated_pop) > archive_size: non_dominated_pop, non_dominated_fit = crowding_distance_crop( non_dominated_pop, non_dominated_fit, archive_size ) return non_dominated_pop, non_dominated_fit def soft_rime_search(rimepop, A, FitA, it, max_iter, lb, ub, dim, W=5): """简化版软霜搜索函数""" if max_iter < 1: raise ValueError("max_iter不能小于1") N = rimepop.shape[0] new_rimepop = rimepop.copy() progress = it / max_iter # 迭代进度(0~1) # 软霜因子 rime_factor = (np.random.rand() - 0.5) * 2 * np.cos((np.pi * it / (max_iter / 5))) * ( 1 - np.round(it * W / max_iter) / W) # 阶段式参数 E = np.sqrt(progress) w = 1 / (1 + np.exp(-9 * (progress - 0.5))) # 引导解选择 if len(A) > 0: guides = A[np.random.choice(len(A), size=N)] else: guides = np.mean(rimepop, axis=0, keepdims=True).repeat(N, axis=0) # 位置更新 for i in range(N): global_guide = guides[i] for j in range(dim): if np.random.rand() < E: # 正向引导(向引导解方向) forward_dir = (1 + rime_factor) * (global_guide[j] - new_rimepop[i, j]) * np.random.rand() # 反向引导 if progress < 0.5: reverse_guide = lb[j] + ub[j] - new_rimepop[i, j] rev_strength = 1 + (1 - progress) * (1 + rime_factor) else: arch_mean = np.mean(A[:, j]) if len(A) > 0 else np.mean(rimepop[:, j]) reverse_guide = arch_mean * 2 - new_rimepop[i, j] rev_strength = 1 + (1 - progress) # 反向引导强度 reverse_dir = rev_strength * (reverse_guide - new_rimepop[i, j]) * (np.random.rand() * 2 - 1) new_rimepop[i, j] += (1 - w) * reverse_dir + w * forward_dir return np.clip(new_rimepop, lb, ub) # 边界修正 def hard_rime_puncture(new_rimepop, A, FitA, objective_func, it, max_iter): """优化后的硬霜穿刺函数""" N, dim = new_rimepop.shape # 存档解太少时不穿刺,直接返回 if len(A) < 2: new_FitX = objective_func(new_rimepop) # 使用向量化目标函数 return new_rimepop, new_FitX # 计算新解的目标值(批量计算) new_FitX = objective_func(new_rimepop) # 向量化Pareto等级计算 def pareto_ranking(pop_fitness): n = len(pop_fitness) ranks = np.zeros(n, dtype=int) dominated_count = np.zeros(n, dtype=int) # 向量化构建支配关系矩阵 dominate_matrix = ( np.all(pop_fitness[:, np.newaxis] <= pop_fitness, axis=2) & np.any(pop_fitness[:, np.newaxis] < pop_fitness, axis=2) ) np.fill_diagonal(dominate_matrix, False) dominated_count = np.sum(dominate_matrix, axis=0) # 计算Pareto等级 current_rank = 1 while np.any(ranks == 0): mask = (dominated_count == 0) & (ranks == 0) if not np.any(mask): break ranks[mask] = current_rank # 向量化更新被支配计数(优化点) dominated_count[np.any(dominate_matrix[mask], axis=0)] -= 1 current_rank += 1 return ranks # 计算新解的Pareto等级 ranks = pareto_ranking(new_FitX) max_rank = np.max(ranks) if np.max(ranks) > 0 else 1 # 平衡触发概率 progress = it / max_iter progress_factor = 1 - progress ** 1.5 if max_rank > 1: normalized_ranks = (ranks - 1) / (max_rank - 1) else: normalized_ranks = np.zeros_like(ranks) trigger_prob = 0.4 + (normalized_ranks * 0.3) trigger_prob = trigger_prob * progress_factor trigger_prob = np.clip(trigger_prob, 0.1, 0.9) # 多样性引导解选择 crowd_dists = crowding_distance(FitA) sorted_indices = np.argsort(-crowd_dists) guide_count = max(3, int(10 - 7 * progress)) top_guides = A[sorted_indices[:guide_count]] # 交叉穿刺 for i in range(N): if np.random.rand() < trigger_prob[i]: idx1, idx2 = np.random.choice(len(top_guides), 2, replace=False) guide1, guide2 = top_guides[idx1], top_guides[idx2] for j in range(dim): if np.random.rand() < 0.5: new_rimepop[i, j] = guide1[j] else: new_rimepop[i, j] = guide2[j] # 重新评估穿刺后的解(批量计算) new_FitX = objective_func(new_rimepop) return new_rimepop, new_FitX # 2. 替换低效的双重循环(mask生成优化,核心优化点2) def update_population(combined_pop, combined_fit, N, gen, Maxgen, Xmin, Xmax): """种群更新:筛选非支配解+对被支配解执行云变异后补充""" # 1. 筛选非支配解 non_dominated_pop, non_dominated_fit = Domination(combined_pop, combined_fit) nd_size = len(non_dominated_pop) if nd_size >= N: # 非支配解充足时,按拥挤距离选择 crowd_nd = crowding_distance(non_dominated_fit) sorted_nd_idx = np.argsort(-crowd_nd) return non_dominated_pop[sorted_nd_idx[:N]] else: # 2. 提取被支配解(需补充的解) remaining = N - nd_size mask = np.zeros(len(combined_pop), dtype=bool) for i, pop in enumerate(combined_pop): for nd_pop in non_dominated_pop: if np.allclose(pop, nd_pop, atol=1e-6): mask[i] = True break other_pop = combined_pop[~mask] # 被支配解 other_fit = combined_fit[~mask] # 3. 对被支配解按拥挤距离排序 if len(other_pop) == 0: # 极端情况:重复选择非支配解 selected = np.random.choice(nd_size, remaining, replace=True) supplementary = non_dominated_pop[selected] else: crowd_other = crowding_distance(other_fit) sorted_other_idx = np.argsort(-crowd_other) selected_other = other_pop[sorted_other_idx[:remaining]] # 4. 对选中的被支配解执行云变异 mutated_supplementary = [] for idx in range(len(selected_other)): mutated = cloud_mutation( X=selected_other, n=idx, gen=gen, Maxgen=Maxgen, Xmin=Xmin, Xmax=Xmax, is_non_dominated=False ) mutated_supplementary.append(mutated) supplementary = np.array(mutated_supplementary) # 5. 组合非支配解和变异后的补充解 return np.vstack([non_dominated_pop, supplementary]) def cloud_mutation(X, n, gen, Maxgen, Xmin, Xmax, is_non_dominated): """优化版云变异""" D = X.shape[1] progress = gen / Maxgen # 差异化变异概率 if is_non_dominated: pm = 0.3 - 0.2 * progress else: pm = 0.8 - 0.4 * progress pm = np.clip(pm, 0.05, 0.95) if np.random.rand() >= pm: return X[n, :].copy() # 差异化变异维度数量 if is_non_dominated: dim_ratio = 0.2 - 0.1 * progress dim_count = max(1, int(dim_ratio * D)) else: dim_ratio = 0.6 - 0.3 * progress dim_count = max(1, int(dim_ratio * D)) mutate_dims = np.random.choice(D, dim_count, replace=False) # 动态调整En y = X[n, :].copy() for i in mutate_dims: Ex = X[n, i] dim_std = np.std(X[:, i]) norm_std = dim_std / ((Xmax[i] - Xmin[i]) / 2 + 1e-10) base_En = (Xmax[i] - Xmin[i]) / (10 * np.sqrt(gen + 1)) En = base_En * (1.2 - 0.8 * norm_std) En = max(En, 1e-10) He = En / 10 Enn = np.random.normal(En, He) Enn = max(abs(Enn), 1e-10) CloudDrp_pos = np.random.normal(Ex, abs(Enn)) CloudDrp_certainty = np.exp(-(CloudDrp_pos - Ex) ** 2 / (2 * Enn ** 2)) adjusted_pos = Ex + (CloudDrp_pos - Ex) * CloudDrp_certainty y[i] = np.clip(adjusted_pos, Xmin[i], Xmax[i]) return y def MORIME_zdt(N, max_iter, lb, ub, dim, objective_func, archive_size, NumObj=2): """多目标霜冰优化算法(适配ZDT测试函数)""" start_time = time.time() # 初始化种群 rimepop = initialization(N, dim, lb, ub) # 计算初始种群的目标函数值(批量计算) FitX = objective_func(rimepop) # 初始化存档 A, FitA = Domination(rimepop, FitX) # 防御性检查 if A.size == 0: rimepop = initialization(N, dim, lb, ub) FitX = objective_func(rimepop) A, FitA = Domination(rimepop, FitX) # 主循环 for it in range(max_iter): if it % 10 == 0: print(f"迭代次数: {it + 1}/{max_iter}") # 阶段1:软霜搜索 new_rimepop = soft_rime_search(rimepop, A, FitA, it, max_iter, lb, ub, dim) # 阶段2:硬霜穿刺 new_rimepop, new_FitX = hard_rime_puncture(new_rimepop, A, FitA, objective_func, it, max_iter) # 更新存档 A, FitA = UpdateArchive(NumObj, archive_size, A, FitA, new_rimepop, new_FitX) # 更新种群 combined_pop = np.vstack([new_rimepop, A]) combined_fit = np.vstack([new_FitX, FitA]) # 调用修改后的种群更新函数 rimepop = update_population( combined_pop, combined_fit, N, gen=it, Maxgen=max_iter, Xmin=lb, Xmax=ub ) # 重新计算种群目标值(批量计算) FitX = objective_func(rimepop) run_time = time.time() - start_time return A, FitA, run_time def plot_results(FitA, true_front, name, igd_value, plot_hull=False): """绘制ZDT函数的Pareto前沿对比图""" plt.figure(figsize=(10, 6)) # 绘制真实帕累托前沿 plt.plot(true_front[:, 0], true_front[:, 1], &#39;b-&#39;, lw=2, label=&#39;真实帕累托前沿&#39;) # 绘制算法得到的非支配解 plt.scatter(FitA[:, 0], FitA[:, 1], c=&#39;red&#39;, s=50, edgecolors=&#39;k&#39;, alpha=0.7, label=&#39;算法求得的非支配解&#39;) # 可选绘制凸包 if plot_hull and len(FitA) > 2: try: hull = ConvexHull(FitA) for simplex in hull.simplices: plt.plot(FitA[simplex, 0], FitA[simplex, 1], &#39;g--&#39;, lw=1, alpha=0.5) except: pass # 图表设置(显示IGD值) plt.xlabel(&#39;目标函数 f1&#39;, fontsize=12) plt.ylabel(&#39;目标函数 f2&#39;, fontsize=12) plt.title(f&#39;{name}函数的Pareto最优前沿对比(IGD={igd_value:.6f})&#39;, fontsize=14) plt.grid(alpha=0.3) plt.legend() plt.show() # 主函数:包含ZDT4边界修正(核心优化点3) def main(N=100, max_iter=100, archive_size=200): # 定义要测试的函数,包含ZDT系列和新增的DTLZ系列 functions = [ # ZDT系列函数 (函数, 名称, 决策变量维度) (zdt1, "ZDT1", 10), (zdt2, "ZDT2", 10), (zdt3, "ZDT3", 10), (zdt4, "ZDT4", 10), (zdt6, "ZDT6", 10), # 新增DTLZ系列函数 (函数, 名称, 决策变量维度) (dtlz1, "DTLZ1", 3), (dtlz2, "DTLZ2", 3), (dtlz3, "DTLZ3", 3) ] # 存储结果 results = {} # 循环运行每个测试函数 for func, name, dim in functions: print(f"\n======================================") print(f"========== {name} 优化开始 ==========") print(f"======================================\n") # 设置各函数的边界 if name == "ZDT4": # ZDT4原始定义:x1∈[0,1], x2~xd∈[-5,5] lb = np.zeros(dim) ub = np.ones(dim) lb[1:] = -5.0 # 第2到第d维下界为-5 ub[1:] = 5.0 # 第2到第d维上界为5 elif name.startswith("ZDT"): # 其他ZDT函数:所有维度∈[0,1] lb = np.zeros(dim) ub = np.ones(dim) else: # DTLZ系列 # DTLZ1/2/3所有维度∈[0,1] lb = np.zeros(dim) ub = np.ones(dim) # 运行优化算法 A, FitA, run_time = MORIME_zdt( N, max_iter, lb, ub, dim, func, archive_size ) # 生成真实帕累托前沿 true_front = generate_true_front(func) # 计算IGD值 igd_value = calculate_igd(FitA, true_front) # 保存结果 results[name] = { &#39;solutions&#39;: A, &#39;fitness&#39;: FitA, &#39;run_time&#39;: run_time, &#39;igd&#39;: igd_value, &#39;true_front&#39;: true_front } # 打印结果统计 print(f"\n========== {name} 优化结果 ==========") print(f"算法运行时间: {run_time:.2f}秒") print(f"找到的非支配解数量: {len(FitA)}") print(f"IGD值: {igd_value:.6f}") # 绘制结果 plot_results(FitA, true_front, name, igd_value, plot_hull=False) # 打印所有IGD值对比 print("\n======================================") print("各函数的IGD值对比:") for func_name, res in results.items(): print(f"{func_name}: {res[&#39;igd&#39;]:.6f}") print("======================================") return results if __name__ == "__main__": results = main() 请你检查我的DTLZ函数是否编写错误,为什么这三个函数的图像都很差
08-19
import numpy as np import math import matplotlib.pyplot as plt 修正后的参数设置 R = 300.0 # 基准球面半径 (m) F = 0.466 * R # 半径差 (m) focal_z = 161.2 # 馈源器坐标 p0 = 279.6 # 初始抛物线参数 N_samples = 1000 # 抛物线采样点数量 convergence_threshold = 1e-7 # 收敛阈值 max_iterations = 10 # 最大迭代次数 var_window = 10 # 方差计算窗口大小 def parabola_point(x, p): “”" 计算抛物线上点的坐标 抛物线方程: x² = 2pz 整理得: z = x²/2p “”" z = (x**2) / (2 * p) return np.array([x, z]) def reflection_point(Q, p): “”" 计算抛物线上点Q的反射线在接收平面上的落点横坐标 接收平面: z = focal_z “”" x, z = Q # 顶点特殊处理 if abs(x) < 1e-5: return 0.0 # 计算抛物线在点Q处的导数 (dz/dx) dz_dx = x / p # 计算法向量 (归一化) norm_vec = np.array([-dz_dx , 1.0]) norm_length = np.linalg.norm(norm_vec) # 处理零向量情况 if norm_length < 1e-10: return 0.0 norm = norm_vec / norm_length # 入射方向 (来自正上方的平行电磁波) incident_dir = np.array([0, -1.0]) # 计算反射方向 dot_product = np.dot(incident_dir, norm) reflection_dir = incident_dir - 2 * dot_product * norm # 避免除零错误 if abs(reflection_dir[2]) < 1e-10: return 0.0 # 计算反射线与接收平面 (z = focal_z) 的交点 intersection_x = (focal_z - p/2 ) * reflection_dir[0] / reflection_dir[1] return intersection_x def reflection_constraint_satisfied(p, verbose=False): “”" 检查反射约束是否满足 所有采样点的反射线落点必须在接收圆盘内 (|x| <= 0.5) “”" # 在抛物线口径内均匀采样 (x从-150到150) x_samples = np.linspace(-150, 150, N_samples) max_violation = 0.0 for x in x_samples: Q = parabola_point(x, p) intersect_x = reflection_point(Q, p) # 记录最大偏差 deviation = abs(intersect_x) if deviation > max_violation: max_violation = deviation # 检查落点是否在接收圆盘内 if deviation > 0.5: if verbose: print(f"p={p}, x={x}: 落点 x={intersect_x:.6f} 超出范围 (偏差={deviation-0.5:.6f})") return False, max_violation return True, max_violation def objective_function(p): “”" 计算目标函数M(p) 量化抛物线与基准球面的逼近程度 “”" # 在抛物线口径内均匀采样 x_samples = np.linspace(-150, 150, N_samples) M = 0.0 for x in x_samples: # 抛物线上点 Q_parabola = parabola_point(x, p) # 计算该点到原点的距离 r = np.linalg.norm(Q_parabola) # 径向投影到基准球面上的点 if r > 0: Q_sphere = (R / r) * Q_parabola else: Q_sphere = np.array([0, -R]) # 顶点处理 # 计算两点之间的距离 (形状差异) distance = np.linalg.norm(Q_parabola - Q_sphere) # 计算权重: (x/150)^2 + 1 weight = (x / 150)**2 + 1 # 累加到目标函数 M += weight * (distance ** 2) return M def find_feasible_p_range(): “”" 确定满足反射约束的p的可行范围 使用更合理的搜索范围和策略 “”" # 设置合理的搜索范围 (基于理论计算) low_bound = 100.0 high_bound = 300.0 search_step = 0.1 feasible_points = [] print("搜索满足反射约束的p范围...") print(f"搜索范围: [{low_bound}, {high_bound}], 步长: {search_step}") # 扫描p值 p_values = np.arange(low_bound, high_bound + search_step, search_step) for p in p_values: satisfied, max_dev = reflection_constraint_satisfied(p) if satisfied: feasible_points.append(p) print(f"p = {p:.2f}: 满足约束 (最大偏差={max_dev:.6f})") else: print(f"p = {p:.2f}: 不满足约束 (最大偏差={max_dev:.6f})") if not feasible_points: raise ValueError("未找到满足反射约束的p值") p_min = min(feasible_points) p_max = max(feasible_points) return p_min, p_max def pso_optimization(p_min, p_max): “”" 使用粒子群优化算法在可行范围内搜索最优p 最小化目标函数M(p) “”" # PSO参数 n_particles = 30 w = 0.5 # 惯性权重 c1 = 1.5 # 个体学习因子 c2 = 2.0 # 社会学习因子 # 初始化粒子 particles = np.random.uniform(p_min, p_max, n_particles) velocities = np.zeros(n_particles) personal_best_positions = particles.copy() personal_best_values = np.array([objective_function(p) for p in particles]) # 全局最优 global_best_idx = np.argmin(personal_best_values) global_best_position = personal_best_positions[global_best_idx] global_best_value = personal_best_values[global_best_idx] # 收敛跟踪 global_best_history = [global_best_value] # PSO主循环 for iteration in range(max_iterations): for i in range(n_particles): # 更新粒子速度和位置 r1, r2 = np.random.random(2) velocities[i] = (w * velocities[i] + c1 * r1 * (personal_best_positions[i] - particles[i]) + c2 * r2 * (global_best_position - particles[i])) particles[i] += velocities[i] # 确保粒子在可行范围内 particles[i] = np.clip(particles[i], p_min, p_max) # 计算目标函数 current_value = objective_function(particles[i]) # 更新个体最优 if current_value < personal_best_values[i]: personal_best_values[i] = current_value personal_best_positions[i] = particles[i] # 更新全局最优 if current_value < global_best_value: global_best_value = current_value global_best_position = particles[i] # 记录全局最优值 global_best_history.append(global_best_value) # 输出当前迭代信息 print(f"Iter {iteration+1}: 全局最优 p = {global_best_position:.6f}, M = {global_best_value:.6f}") # 检查收敛性 (最近var_window次迭代的方差) if len(global_best_history) > var_window: last_values = global_best_history[-var_window:] variance = np.var(last_values) if variance < convergence_threshold: print(f"在第{iteration+1}次迭代收敛") break return global_best_position, global_best_value 主程序 def main(): print(“===== FAST主动反射面形状调节 - 问题1求解 =====”) print(f"初始抛物线参数 p0 = {p0}“) print(f"焦点位置: (0, 0, {focal_z})”) # 测试初始值 print("\n测试初始值 p0 =", p0) satisfied, max_dev = reflection_constraint_satisfied(p0, verbose=True) if satisfied: print(f"✅ p0={p0} 满足反射约束 (最大偏差={max_dev:.6f})") else: print(f"❌ p0={p0} 不满足反射约束 (最大偏差={max_dev:.6f})") # 步骤1: 确定满足反射约束的p范围 print("\n步骤1: 确定满足反射约束的p范围...") try: p_min, p_max = find_feasible_p_range() print(f"\n反射约束确定的可行范围: p ∈ [{p_min:.6f}, {p_max:.6f}]") except ValueError as e: print(f"错误: {e}") # 如果自动搜索失败,使用理论值作为备选 p_min, p_max = 139.7, 139.9 print(f"使用理论范围: p ∈ [{p_min:.6f}, {p_max:.6f}]") # 步骤2: 在可行范围内使用PSO优化目标函数 print("\n步骤2: 在可行范围内使用PSO优化...") best_p, best_M = pso_optimization(p_min, p_max) # 计算顶点坐标 vertex_z = focal_z - best_p # 输出结果 print("\n===== 优化结果 =====") print(f"最优抛物面参数 p = {best_p:.6f}") print(f"顶点坐标: (0, 0, {vertex_z:.6f})") print(f"焦点坐标: (0, 0, {focal_z:.6f})") print(f"目标函数值 M(p) = {best_M:.6f}") # 验证反射约束 print("\n===== 反射约束验证 =====") satisfied, max_dev = reflection_constraint_satisfied(best_p, verbose=True) if satisfied: print(f"✅ 所有采样点的反射线均落在接收圆盘内 (最大偏差={max_dev:.6f})") else: print(f"❌ 警告: 部分采样点的反射线未落在接收圆盘内 (最大偏差={max_dev:.6f})") # 可视化结果 visualize_results(best_p) def visualize_results(p): “”“可视化优化结果 (使用英文标签)”“” # 创建抛物线采样点 x_samples = np.linspace(-150, 150, 100) parabola_points = [parabola_point(x, p) for x in x_samples] parabola_z = [point[2] for point in parabola_points] # 创建基准球面采样点 sphere_points = [] for x in x_samples: # 基准球面上对应点 (只考虑下半球) if abs(x) <= R: z_sphere = -np.sqrt(R**2 - x**2) sphere_points.append([x, z_sphere]) sphere_z = [point[2] for point in sphere_points] # 绘制抛物线与基准球面 plt.figure(figsize=(10, 6)) plt.plot(x_samples, parabola_z, &#39;b-&#39;, linewidth=2, label=f&#39;Optimized Parabola (p={p:.2f})&#39;) plt.plot(x_samples[:len(sphere_z)], sphere_z, &#39;r--&#39;, linewidth=2, label=&#39;Reference Sphere&#39;) plt.scatter([0], [focal_z - p], c=&#39;g&#39;, s=100, label=&#39;Parabola Vertex&#39;) plt.scatter([0], [focal_z], c=&#39;m&#39;, s=100, marker=&#39;*&#39;, label=&#39;Focal Point&#39;) # 标注接收圆盘 plt.axhline(y=focal_z, color=&#39;k&#39;, linestyle=&#39;-&#39;, alpha=0.3) plt.fill_between([-0.5, 0.5], focal_z-0.1, focal_z+0.1, color=&#39;c&#39;, alpha=0.3, label=&#39;Receiver Disk&#39;) plt.xlabel(&#39;x (m)&#39;, fontsize=12) plt.ylabel(&#39;z (m)&#39;, fontsize=12) plt.title(&#39;Optimized Parabola vs Reference Sphere&#39;, fontsize=14) plt.legend(loc=&#39;best&#39;) plt.grid(True, linestyle=&#39;--&#39;, alpha=0.7) plt.xlim(-160, 160) plt.ylim(-320, -120) plt.tight_layout() plt.savefig(&#39;parabola_comparison.png&#39;, dpi=300) plt.show() # 绘制反射线落点分布 plt.figure(figsize=(10, 6)) x_samples_detail = np.linspace(-150, 150, N_samples) intersections = [reflection_point(parabola_point(x, p), p) for x in x_samples_detail] plt.plot(x_samples_detail, intersections, &#39;b.&#39;, markersize=3, alpha=0.7) plt.axhline(y=0.5, color=&#39;r&#39;, linestyle=&#39;--&#39;, label=&#39;Receiver Boundary&#39;) plt.axhline(y=-0.5, color=&#39;r&#39;, linestyle=&#39;--&#39;) plt.axhline(y=0, color=&#39;g&#39;, linestyle=&#39;-&#39;, alpha=0.5, label=&#39;Center Line&#39;) plt.fill_between(x_samples_detail, -0.5, 0.5, color=&#39;c&#39;, alpha=0.1, label=&#39;Effective Area&#39;) plt.xlabel(&#39;Parabola Point x-coordinate (m)&#39;, fontsize=12) plt.ylabel(&#39;Intersection x-coordinate (m)&#39;, fontsize=12) plt.title(&#39;Reflection Points on Receiver Plane&#39;, fontsize=14) plt.legend(loc=&#39;best&#39;) plt.grid(True, linestyle=&#39;--&#39;, alpha=0.7) plt.ylim(-1.0, 1.0) plt.tight_layout() plt.savefig(&#39;reflection_points.png&#39;, dpi=300) plt.show() if name == “main”: main() 上面这段python代码有问题,我的所有物理活动是发生在二维平面的,你给我重新输出代码让我能直接在pycharm上正确的运行这段程序,代码中只有p/2是抛物线的焦点,F等与抛物线无关,抛物线的原点与基准球面的原点重合,是(0,0)
08-15
我自己对第二问进行了求解,现在我需要你把我修改一下需求: 确定储水罐位置和容量 初始假设储水罐容量为40,000L,初始成本为0,只有当容量超过40000L时,超出部分才按每升5元的成本计算。 根据农场布局和作物分布,确定储水罐的最佳位置 储水罐应位于能够覆盖尽可能多作物的位置 2. 基于储水罐位置确定喷头布局 确保喷头与储水罐的覆盖范围不重叠 喷头布局应确保全覆盖农场且满足最小间距要求 3. 管道网络优化 从河流引水到储水罐和喷头网络 优化管道布局以最小化成本:import numpy as np import pandas as pd from scipy.optimize import minimize, Bounds from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score import matplotlib.pyplot as plt import warnings import logging import math from matplotlib.patches import Circle, Rectangle from sklearn.preprocessing import StandardScaler import os import networkx as nx from scipy.sparse.csgraph import minimum_spanning_tree # 设置日志记录 logging.basicConfig(level=logging.INFO, format=&#39;%(asctime)s - %(levelname)s - %(message)s&#39;) logger = logging.getLogger(__name__) # 忽略特定警告 warnings.filterwarnings("ignore", category=UserWarning, module="sklearn") warnings.filterwarnings("ignore", category=RuntimeWarning) # 设置中文字体 plt.rcParams[&#39;font.sans-serif&#39;] = [&#39;SimHei&#39;] # 使用黑体 plt.rcParams[&#39;axes.unicode_minus&#39;] = False # 正常显示负号 # ==================== 常量定义 ==================== DEPTH = 0.05 # 5 cm DRY_DENS = 1500 # kg/m&sup3; 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 = { &#39;高粱&#39;: 0.5, &#39;玉米&#39;: 0.3, &#39;大豆&#39;: 0.2 } # 计算垂直分布的高度边界 total_height = FARM_SIZE_Y gaoliang_height = total_height * CROP_AREAS[&#39;高粱&#39;] # 高粱区域高度 yumi_height = total_height * CROP_AREAS[&#39;玉米&#39;] # 玉米区域高度 dadou_height = total_height * CROP_AREAS[&#39;大豆&#39;] # 大豆区域高度 # 作物区域垂直分布 (高粱在下,玉米在中,大豆在上) CROP_REGIONS = { &#39;高粱&#39;: {&#39;x_min&#39;: 0, &#39;x_max&#39;: FARM_SIZE_X, &#39;y_min&#39;: 0, &#39;y_max&#39;: gaoliang_height}, &#39;玉米&#39;: {&#39;x_min&#39;: 0, &#39;x_max&#39;: FARM_SIZE_X, &#39;y_min&#39;: gaoliang_height, &#39;y_max&#39;: gaoliang_height + yumi_height}, &#39;大豆&#39;: {&#39;x_min&#39;: 0, &#39;x_max&#39;: FARM_SIZE_X, &#39;y_min&#39;: gaoliang_height + yumi_height, &#39;y_max&#39;: FARM_SIZE_Y} } SPRINKLER_RADIUS = 15 # 喷头半径 RIVER_POSITION = &#39;south&#39; 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&sup3; # 系统参数 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 = &#39;附件/该地土壤湿度数据.xlsx&#39; if not os.path.exists(file_path): logger.error(f"土壤湿度数据文件不存在: {file_path}") dates = pd.date_range(&#39;2021-07-01&#39;, 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=&#39;JingYueTan&#39;) required_columns = [&#39;DATE&#39;, &#39;5cm_SM&#39;] if not all(col in data.columns for col in required_columns): logger.error(f"Excel文件中缺少必要的列: {required_columns}") dates = pd.date_range(&#39;2021-07-01&#39;, 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[&#39;DATE&#39;] = pd.to_datetime(data[&#39;DATE&#39;]) data.set_index(&#39;DATE&#39;, inplace=True) start_date = pd.Timestamp(&#39;2021-07-01&#39;) end_date = pd.Timestamp(&#39;2021-07-31&#39;) 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[&#39;5cm_SM&#39;].values, &#39;b-&#39;, linewidth=2) plt.axhline(y=MIN_SOIL_MOISTURE, color=&#39;r&#39;, linestyle=&#39;--&#39;, label=&#39;最低土壤湿度阈值&#39;) plt.title(&#39;2021年7月土壤湿度变化&#39;) plt.xlabel(&#39;日期&#39;) plt.ylabel(&#39;土壤湿度 (5cm_SM)&#39;) plt.legend() plt.grid(True) plt.xticks(rotation=45) plt.tight_layout() plt.savefig(&#39;土壤湿度变化图.png&#39;, dpi=300) plt.show() return july_data[&#39;5cm_SM&#39;] except Exception as e: logger.error(f"加载土壤湿度数据时出错: {e}") dates = pd.date_range(&#39;2021-07-01&#39;, 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[&#39;max_demand&#39;] = sprinkler_df[&#39;area&#39;] * 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=&#39;灌溉需求&#39;) plt.title(&#39;2021年7月每日灌溉需求&#39;) plt.xlabel(&#39;日期&#39;) plt.ylabel(&#39;灌溉需求 (L/m²)&#39;) plt.legend() plt.grid(True) plt.xticks(rotation=45) plt.tight_layout() plt.savefig(&#39;灌溉需求变化图.png&#39;, 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[&#39;x_min&#39;] <= x <= region[&#39;x_max&#39;] and region[&#39;y_min&#39;] <= y <= region[&#39;y_max&#39;]: crop_type = crop break sprinklers.append({ &#39;id&#39;: len(sprinklers), &#39;x&#39;: x, &#39;y&#39;: y, &#39;radius&#39;: radius, &#39;crop_type&#39;: 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[&#39;x&#39;], sprinkler[&#39;y&#39;] 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[&#39;area&#39;] = areas return sprinkler_df def validate_sprinkler_spacing(sprinkler_df, min_spacing=15): """验证喷头间距是否≥15m""" points = sprinkler_df[[&#39;x&#39;, &#39;y&#39;]].values num_sprinklers = len(points) min_distance = float(&#39;inf&#39;) 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[&#39;x&#39;], sprinkler_df[&#39;y&#39;], c=&#39;blue&#39;, s=50, label=&#39;喷头&#39;) 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]], &#39;r--&#39;, linewidth=2, label=f&#39;最小间距: {min_distance:.2f}m&#39;) for i, row in sprinkler_df.iterrows(): plt.text(row[&#39;x&#39;], row[&#39;y&#39;], f"S{i}", fontsize=9, ha=&#39;center&#39;, va=&#39;bottom&#39;) plt.text(row[&#39;x&#39;], row[&#39;y&#39;], f"({row[&#39;x&#39;]:.1f},{row[&#39;y&#39;]:.1f})", fontsize=8, ha=&#39;center&#39;, va=&#39;top&#39;) # 绘制垂直分布的作物区域 colors = {&#39;高粱&#39;: &#39;lightgreen&#39;, &#39;玉米&#39;: &#39;lightyellow&#39;, &#39;大豆&#39;: &#39;lightblue&#39;} # 高粱区域(底部) rect_gaoliang = Rectangle( (CROP_REGIONS[&#39;高粱&#39;][&#39;x_min&#39;], CROP_REGIONS[&#39;高粱&#39;][&#39;y_min&#39;]), CROP_REGIONS[&#39;高粱&#39;][&#39;x_max&#39;] - CROP_REGIONS[&#39;高粱&#39;][&#39;x_min&#39;], CROP_REGIONS[&#39;高粱&#39;][&#39;y_max&#39;] - CROP_REGIONS[&#39;高粱&#39;][&#39;y_min&#39;], alpha=0.3, color=colors[&#39;高粱&#39;], label=f&#39;高粱 ({CROP_AREAS["高粱"]}公顷)&#39; ) plt.gca().add_patch(rect_gaoliang) # 玉米区域(中部) rect_yumi = Rectangle( (CROP_REGIONS[&#39;玉米&#39;][&#39;x_min&#39;], CROP_REGIONS[&#39;玉米&#39;][&#39;y_min&#39;]), CROP_REGIONS[&#39;玉米&#39;][&#39;x_max&#39;] - CROP_REGIONS[&#39;玉米&#39;][&#39;x_min&#39;], CROP_REGIONS[&#39;玉米&#39;][&#39;y_max&#39;] - CROP_REGIONS[&#39;玉米&#39;][&#39;y_min&#39;], alpha=0.3, color=colors[&#39;玉米&#39;], label=f&#39;玉米 ({CROP_AREAS["玉米"]}公顷)&#39; ) plt.gca().add_patch(rect_yumi) # 大豆区域(顶部) rect_dadou = Rectangle( (CROP_REGIONS[&#39;大豆&#39;][&#39;x_min&#39;], CROP_REGIONS[&#39;大豆&#39;][&#39;y_min&#39;]), CROP_REGIONS[&#39;大豆&#39;][&#39;x_max&#39;] - CROP_REGIONS[&#39;大豆&#39;][&#39;x_min&#39;], CROP_REGIONS[&#39;大豆&#39;][&#39;y_max&#39;] - CROP_REGIONS[&#39;大豆&#39;][&#39;y_min&#39;], alpha=0.3, color=colors[&#39;大豆&#39;], label=f&#39;大豆 ({CROP_AREAS["大豆"]}公顷)&#39; ) plt.gca().add_patch(rect_dadou) plt.plot([0, FARM_SIZE_X], [0, 0], &#39;b-&#39;, linewidth=4, label=&#39;河流&#39;) plt.title(f&#39;喷头布局图 (最小间距: {min_distance:.2f}m)&#39;) plt.xlabel(&#39;X坐标 (m)&#39;) plt.ylabel(&#39;Y坐标 (m)&#39;) plt.grid(True) plt.legend() plt.tight_layout() plt.savefig(&#39;喷头布局验证图.png&#39;, 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][&#39;max_demand&#39;]) # 创建所有喷头之间的连接(完全网状) for i in range(n): for j in range(i+1, n): x1, y1 = sprinkler_df.iloc[i][[&#39;x&#39;, &#39;y&#39;]] x2, y2 = sprinkler_df.iloc[j][[&#39;x&#39;, &#39;y&#39;]] 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=&#39;length&#39;, data=True)) # 计算每条边的流量 edge_flows = {} visited = set() def dfs(node): visited.add(node) total_flow = sprinkler_df.iloc[node][&#39;max_demand&#39;] 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[[&#39;x&#39;, &#39;y&#39;]].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, &#39;bo-&#39;) plt.xlabel(&#39;聚类数量&#39;) plt.ylabel(&#39;SSE&#39;) plt.title(&#39;肘部法&#39;) plt.grid(True) plt.subplot(1, 2, 2) plt.plot(k_range[1:], silhouette_scores[1:], &#39;ro-&#39;) plt.xlabel(&#39;聚类数量&#39;) plt.ylabel(&#39;轮廓系数&#39;) plt.title(&#39;轮廓系数法&#39;) plt.grid(True) plt.tight_layout() plt.savefig(&#39;聚类数量选择.png&#39;, 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[[&#39;x&#39;, &#39;y&#39;]].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(&#39;inf&#39;) closest_sprinkler_idx = -1 # 找到最近的喷头 for j, sprinkler in sprinkler_df.iterrows(): dist = np.sqrt((tank_pos[0] - sprinkler[&#39;x&#39;])**2 + (tank_pos[1] - sprinkler[&#39;y&#39;])**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[&#39;x&#39;] dy = tank_pos[1] - closest_sprinkler[&#39;y&#39;] angle = np.arctan2(dy, dx) # 将储水罐移动到最小距离之外 new_x = closest_sprinkler[&#39;x&#39;] + np.cos(angle) * MIN_DISTANCE_BETWEEN_SPRINKLER_TANK new_y = closest_sprinkler[&#39;y&#39;] + 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] # 确保储水罐之间也有足够距离 for i in range(len(candidate_positions)): for j in range(i+1, len(candidate_positions)): dist = np.sqrt((candidate_positions[i][0] - candidate_positions[j][0])**2 + (candidate_positions[i][1] - candidate_positions[j][1])**2) if dist < MIN_DISTANCE_BETWEEN_SPRINKLER_TANK: # 调整位置使它们分开 dx = candidate_positions[j][0] - candidate_positions[i][0] dy = candidate_positions[j][1] - candidate_positions[i][1] angle = np.arctan2(dy, dx) # 将第二个储水罐移动到最小距离之外 new_x = candidate_positions[i][0] + np.cos(angle) * MIN_DISTANCE_BETWEEN_SPRINKLER_TANK new_y = candidate_positions[i][1] + 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[j] = [new_x, new_y] # 创建标签(所有喷头都分配到最近的储水罐) labels = [] for i, sprinkler in sprinkler_df.iterrows(): min_dist = float(&#39;inf&#39;) closest_tank_idx = -1 for j, tank_pos in enumerate(candidate_positions): dist = np.sqrt((sprinkler[&#39;x&#39;] - tank_pos[0])**2 + (sprinkler[&#39;y&#39;] - 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[[&#39;x&#39;, &#39;y&#39;]].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][&#39;max_demand&#39;]) # 创建所有喷头之间的连接 for i in range(n_sprinklers): for j in range(i+1, n_sprinklers): x1, y1 = sprinkler_df.iloc[i][[&#39;x&#39;, &#39;y&#39;]] x2, y2 = sprinkler_df.iloc[j][[&#39;x&#39;, &#39;y&#39;]] 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=&#39;length&#39;, data=True)) # 找到离河流最近的喷头作为喷头网络入口 distances_to_river = [] for i in range(n_sprinklers): x, y = sprinkler_df.iloc[i][[&#39;x&#39;, &#39;y&#39;]] 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=&#39;length&#39;, data=True)) # 找到离河流最近的储水罐作为储水罐网络入口 distances_to_river_tank = [np.sqrt((x - RIVER_POINT[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[[&#39;x&#39;, &#39;y&#39;]].values, np.array(tank_positions) ) assignments = np.argmin(distances, axis=1) # 计算每个储水罐的需求 sprinkler_max_demands = sprinkler_df[&#39;max_demand&#39;].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][[&#39;x&#39;, &#39;y&#39;]] x2, y2 = sprinkler_df.iloc[j][[&#39;x&#39;, &#39;y&#39;]] 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][&#39;x&#39;] - tank_positions[root_tank_idx][0])**2 + (sprinkler_df.iloc[root_sprinkler_idx][&#39;y&#39;] - 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[&#39;length&#39;] # 计算两个储水罐之间的流量(取最大值) 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[&#39;max_demand&#39;]) 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[&#39;max_demand&#39;]) 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=&#39;SLSQP&#39;, options={&#39;disp&#39;: True, &#39;ftol&#39;: 1e-6, &#39;maxiter&#39;: 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][[&#39;x&#39;, &#39;y&#39;]] 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][[&#39;x&#39;, &#39;y&#39;]] x2, y2 = sprinkler_df.iloc[j][[&#39;x&#39;, &#39;y&#39;]] 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][&#39;x&#39;] - tank_positions[root_tank_idx][0])**2 + (sprinkler_df.iloc[root_sprinkler_idx][&#39;y&#39;] - 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[&#39;length&#39;] 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 = { &#39;tank_cost&#39;: tank_cost, &#39;main_pipe_cost&#39;: main_pipe_cost, &#39;sprinkler_pipe_cost&#39;: sprinkler_pipe_cost, &#39;tank_pipe_cost&#39;: 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.set_aspect(&#39;equal&#39;) # 绘制农场边界和河流 ax.plot([0, FARM_SIZE_X, FARM_SIZE_X, 0, 0], [0, 0, FARM_SIZE_Y, FARM_SIZE_Y, 0], &#39;k-&#39;, linewidth=2) ax.plot([0, FARM_SIZE_X], [0, 0], &#39;b-&#39;, linewidth=4, label=&#39;河流&#39;) ax.plot(RIVER_POINT[0], RIVER_POINT[1], &#39;bo&#39;, markersize=10, label=&#39;取水点&#39;) # 绘制垂直分布的作物区域 colors = {&#39;高粱&#39;: &#39;lightgreen&#39;, &#39;玉米&#39;: &#39;lightyellow&#39;, &#39;大豆&#39;: &#39;lightblue&#39;} # 高粱区域(底部0-50m) rect_gaoliang = Rectangle( (0, 0), FARM_SIZE_X, 50, alpha=0.2, color=colors[&#39;高粱&#39;], label=&#39;高粱 (0.5公顷)&#39; ) ax.add_patch(rect_gaoliang) # 玉米区域(中部50-80m) rect_yumi = Rectangle( (0, 50), FARM_SIZE_X, 30, alpha=0.2, color=colors[&#39;玉米&#39;], label=&#39;玉米 (0.3公顷)&#39; ) ax.add_patch(rect_yumi) # 大豆区域(顶部80-100m) rect_dadou = Rectangle( (0, 80), FARM_SIZE_X, 20, alpha=0.2, color=colors[&#39;大豆&#39;], label=&#39;大豆 (0.2公顷)&#39; ) ax.add_patch(rect_dadou) # 绘制喷头 for i, sprinkler in sprinkler_df.iterrows(): if i == root_sprinkler_idx: ax.plot(sprinkler[&#39;x&#39;], sprinkler[&#39;y&#39;], &#39;o&#39;, color=&#39;red&#39;, markersize=8, markeredgecolor=&#39;black&#39;, markeredgewidth=1.5, label=&#39;喷头网络入口&#39;) else: ax.plot(sprinkler[&#39;x&#39;], sprinkler[&#39;y&#39;], &#39;o&#39;, color=&#39;blue&#39;, markersize=6) # 绘制喷头覆盖范围(确保是圆形) circle = Circle((sprinkler[&#39;x&#39;], sprinkler[&#39;y&#39;]), SPRINKLER_RADIUS, color=&#39;blue&#39;, alpha=0.1, fill=True) ax.add_patch(circle) # 绘制喷头之间的连接(MST边) for i, j, data in sprinkler_mst: x1, y1 = sprinkler_df.iloc[i][[&#39;x&#39;, &#39;y&#39;]] x2, y2 = sprinkler_df.iloc[j][[&#39;x&#39;, &#39;y&#39;]] ax.plot([x1, x2], [y1, y2], &#39;b-&#39;, linewidth=1.5, alpha=0.7, label=&#39;喷头间管道&#39; 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], &#39;s&#39;, color=&#39;red&#39;, markersize=12, markeredgecolor=&#39;black&#39;, markeredgewidth=2, label=&#39;储水罐网络入口&#39;) else: ax.plot(tank[0], tank[1], &#39;s&#39;, color=&#39;purple&#39;, markersize=10, markeredgecolor=&#39;black&#39;, markeredgewidth=1.5) # 绘制储水罐覆盖范围(确保是圆形) circle = Circle((tank[0], tank[1]), TANK_COVERAGE_RADIUS, color=&#39;purple&#39;, alpha=0.15, fill=True) ax.add_patch(circle) ax.text(tank[0], tank[1], f&#39;T{j+1}&#39;, fontsize=10, ha=&#39;center&#39;, va=&#39;center&#39;, fontweight=&#39;bold&#39;) # 绘制储水罐之间的连接(MST边) for i, j, data in tank_mst: x1, y1 = tank_positions[i] x2, y2 = tank_positions[j] ax.plot([x1, x2], [y1, y2], &#39;g-&#39;, linewidth=2, label=&#39;储水罐间管道&#39; if i == 0 and j == 1 else "") # 标注管道长度 mid_x = (x1 + x2) / 2 mid_y = (y1 + y2) / 2 ax.text(mid_x, mid_y, f&#39;{data["length"]:.1f}m&#39;, fontsize=8, bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.7)) # 绘制总管(河流到分流点) connection_point = sprinkler_df.iloc[root_sprinkler_idx][[&#39;x&#39;, &#39;y&#39;]].values ax.plot([RIVER_POINT[0], connection_point[0]], [RIVER_POINT[1], connection_point[1]], &#39;r-&#39;, linewidth=3, label=&#39;总管&#39;) # 标注总管信息 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&#39;{length:.1f}m\n{Q_total_opt:.0f}L/d\n({Q_m3:.2f}m&sup3;/d)&#39;, fontsize=9, fontweight=&#39;bold&#39;, 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]], &#39;m--&#39;, linewidth=2, label=&#39;分流点到储水罐网络&#39;) # 标注分流点到储水罐管道信息 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&#39;{L:.1f}m\n{Q_total_tank:.0f}L/d\n({Q_m3:.2f}m&sup3;/d)&#39;, fontsize=9, fontweight=&#39;bold&#39;, bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.8)) ax.set_xlabel(&#39;X坐标 (m)&#39;) ax.set_ylabel(&#39;Y坐标 (m)&#39;) ax.set_title(&#39;灌溉系统布线规划图&#39;) # 修改标题 handles, labels = ax.get_legend_handles_labels() by_label = dict(zip(labels, handles)) ax.legend(by_label.values(), by_label.keys(), loc=&#39;best&#39;) ax.grid(True) plt.tight_layout() plt.savefig(&#39;灌溉系统布线规划图.png&#39;, dpi=300) # 修改保存文件名 plt.show() def plot_cost_breakdown(cost_breakdown): """绘制成本分解图""" labels = [&#39;储水罐成本&#39;, &#39;总管成本&#39;, &#39;喷头网络管道成本&#39;, &#39;储水罐管道成本&#39;] sizes = [ cost_breakdown[&#39;tank_cost&#39;], cost_breakdown[&#39;main_pipe_cost&#39;], cost_breakdown[&#39;sprinkler_pipe_cost&#39;], cost_breakdown[&#39;tank_pipe_cost&#39;] ] colors = [&#39;lightblue&#39;, &#39;lightcoral&#39;, &#39;lightgreen&#39;, &#39;lightyellow&#39;] plt.figure(figsize=(10, 8)) plt.pie(sizes, labels=labels, colors=colors, autopct=&#39;%1.1f%%&#39;, startangle=90) plt.axis(&#39;equal&#39;) plt.title(&#39;成本分解&#39;) plt.savefig(&#39;成本分解图.png&#39;, 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({ &#39;储水罐编号&#39;: range(1, num_tanks+1), &#39;X坐标&#39;: [f"{tank[0]:.1f}" for tank in tank_positions], &#39;Y坐标&#39;: [f"{tank[1]:.1f}" for tank in tank_positions], &#39;容量(L)&#39;: [f"{v:.0f}" for v in V_opt], &#39;容量(m&sup3;)&#39;: [f"{v * L_TO_M3:.2f}" for v in V_opt], &#39;需求(L/天)&#39;: [f"{d:.0f}" for d in tank_demands], &#39;需求(m&sup3;/天)&#39;: [f"{d * L_TO_M3:.2f}" for d in tank_demands], &#39;分配到储水罐流量(L/天)&#39;: [f"{q:.0f}" for q in Q_tank_network_opt], &#39;分配到储水罐流量(m&sup3;/天)&#39;: [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[&#39;覆盖喷头数&#39;] = coverage_stats tank_pipe_info = [] for i, j, data in tank_mst: tank_pipe_info.append(f&#39;T{i+1}-T{j+1}: {data["length"]:.1f}m&#39;) results_df[&#39;储水罐间管道&#39;] = [&#39;, &#39;.join(tank_pipe_info)] * num_tanks logger.info("\n优化结果详情:") print(results_df.to_string(index=False)) results_df.to_csv(&#39;灌溉系统优化结果.csv&#39;, index=False, encoding=&#39;utf-8-sig&#39;) logger.info("结果已保存到 &#39;灌溉系统优化结果.csv&#39;") total_demand = np.sum(sprinkler_df[&#39;max_demand&#39;]) 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&sup3;/天)") logger.info(f"总管供水能力: {river_supply:.0f} L/天 ({river_supply * L_TO_M3:.2f} m&sup3;/天)") logger.info(f"分配到喷头网络流量: {Q_sprinkler_opt:.0f} L/天 ({Q_sprinkler_opt * L_TO_M3:.2f} m&sup3;/天)") logger.info(f"分配到储水罐网络流量: {np.sum(Q_tank_network_opt):.0f} L/天 ({np.sum(Q_tank_network_opt) * L_TO_M3:.2f} m&sup3;/天)") logger.info(f"储水罐总容量: {tank_supply:.0f} L ({tank_supply * L_TO_M3:.2f} m&sup3;)") logger.info(f"系统可靠性: {(river_supply + tank_supply)/total_demand*100:.1f}%") logger.info(f"\n成本详情:") logger.info(f"储水罐成本: {cost_breakdown[&#39;tank_cost&#39;]:.2f} 元") logger.info(f"总管成本: {cost_breakdown[&#39;main_pipe_cost&#39;]:.2f} 元") logger.info(f"喷头网络管道成本: {cost_breakdown[&#39;sprinkler_pipe_cost&#39;]:.2f} 元") logger.info(f"储水罐管道成本: {cost_breakdown[&#39;tank_pipe_cost&#39;]:.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. 喷头间距验证: {&#39;通过&#39; if spacing_ok else &#39;失败&#39;}") logger.info(f"2. 系统可靠性: {total_cost:.2f} 元") else: logger.error("优化失败:", result.message) if __name__ == "__main__": main()
最新发布
08-25
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值