✅博主简介：擅长数据搜集与处理、建模仿真、程序设计、仿真代码、论文写作与指导，毕业论文、期刊论文经验交流。

✅成品或者定制，扫描文章底部微信二维码。

(1) GA-PSO串行混合优化算法与点核积分快速计算方法

辐射屏蔽优化设计的目标是在满足辐射防护标准的前提下，最小化屏蔽体的重量或成本。这是一个典型的约束优化问题，优化变量包括屏蔽层的材料选择和厚度配置，约束条件为屏蔽体外表面的剂量率不超过允许限值。传统的屏蔽设计方法依赖于工程师的经验和试错，难以保证获得最优方案，而且对于复杂的多层屏蔽结构设计效率低下。将智能优化算法与辐射输运计算相结合，可以实现屏蔽方案的自动寻优，但优化过程需要反复调用输运计算评估屏蔽效果，计算成本是制约优化效率的关键因素。

遗传算法和粒子群算法是两种广泛应用的群体智能优化方法，各有其优势和局限。遗传算法通过选择、交叉、变异等操作模拟自然进化过程，具有较强的全局搜索能力，不易陷入局部最优。但遗传算法的收敛速度相对较慢，在搜索后期精细化能力不足。粒子群算法通过粒子之间的信息共享实现协作搜索，收敛速度快，局部搜索能力强。但粒子群算法的全局搜索能力相对较弱，容易在搜索初期就陷入局部最优。将两种算法串行组合，可以发挥各自优势实现优势互补。

GPPG算法是遗传算法与粒子群算法的串行混合方案，同时结合了点核积分法和蒙特卡洛方法两种输运计算手段。算法分为两个阶段，第一阶段采用遗传算法结合点核积分法进行快速粗搜索。点核积分法是一种确定性的辐射输运计算方法，将辐射源视为点源，通过积分公式计算屏蔽后的剂量分布。该方法计算速度快，单次评估仅需毫秒级时间，但精度相对较低，特别是对于复杂几何和强散射问题存在较大误差。遗传算法利用点核积分的快速评估进行大范围搜索，快速定位有希望的屏蔽方案区域。

第二阶段采用粒子群算法结合蒙特卡洛方法进行精细优化。蒙特卡洛方法通过统计模拟大量粒子的随机输运过程来计算辐射场分布，能够准确处理各种复杂情况，是辐射输运计算的黄金标准。但蒙特卡洛方法计算耗时长，单次评估可能需要数分钟甚至更长时间。第二阶段以第一阶段获得的优良解作为初始种群，利用粒子群算法的局部搜索能力在小范围内进行精细寻优。由于搜索范围已经被第一阶段大幅缩小，蒙特卡洛评估的调用次数可以控制在可接受范围内。两阶段的衔接通过种群迁移实现，第一阶段结束时提取最优的若干个体作为第二阶段的初始粒子。

点核积分程序的开发是GPPG算法实现的关键环节。根据伽马射线与物质相互作用的物理机制，编写了考虑光电效应、康普顿散射和电子对效应的衰减系数数据库。针对多层平板屏蔽的典型几何，推导了点核积分的解析表达式。对于复杂几何则采用数值积分方法，将屏蔽体离散为小体元后累加各体元的贡献。程序经过与蒙特卡洛标准程序的对比验证，在简单几何情况下误差控制在百分之十以内，满足快速筛选的精度要求。

(2) 多目标优化设计方案与算法性能对比分析

实际的辐射屏蔽优化往往需要同时考虑多个目标，最常见的是屏蔽效果与屏蔽体重量之间的权衡。更好的屏蔽效果通常需要更厚或密度更高的屏蔽材料，这会导致屏蔽体重量增加。对于航天器、船舶等对重量敏感的应用场合，需要在满足防护标准的同时尽可能减轻重量。根据不同的应用需求，设计了三种优化方案：方案一侧重屏蔽体重量最小化，将重量作为主要优化目标，剂量率作为硬约束；方案二侧重屏蔽效果最优化，将剂量率最小化作为主要目标，在重量预算约束下追求最佳防护效果；方案三同等重要权衡两个目标，采用加权和方法将重量和剂量率归一化后求和作为综合目标函数。

基准模型的构建为算法对比提供了统一的测试平台。选取一个典型的圆柱形辐射源作为屏蔽对象，源强和能谱参考实际反应堆参数设置。屏蔽体为同心圆筒结构，由内向外可配置多层不同材料。候选屏蔽材料包括铅、钢、混凝土、硼化聚乙烯等常用屏蔽材料，各材料的物理属性和屏蔽特性从标准数据库获取。优化变量为各层的材料类型和厚度，变量编码采用混合方式，材料类型用整数编码，厚度用实数编码。目标函数和约束条件根据三种方案分别设置。

将GPPG算法与三种单一算法进行对比，包括遗传算法结合蒙特卡洛、粒子群算法结合蒙特卡洛以及遗传算法结合点核积分。对比指标包括最终解的质量、收敛速度和计算时间。在三种优化方案的基准模型上分别运行各算法，每种配置独立运行二十次以评估算法的稳定性。结果显示，GPPG算法在解的质量方面与遗传算法结合蒙特卡洛接近，显著优于另外两种方法。在计算时间方面，GPPG算法仅为遗传算法结合蒙特卡洛的约三分之一，因为第一阶段的点核积分快速评估大幅减少了蒙特卡洛调用次数。在收敛稳定性方面，GPPG算法的二十次运行结果方差最小，说明算法具有良好的鲁棒性。

群体表现分析揭示了GPPG算法的内在优化机制。通过记录各代种群的适应度分布，观察种群从初始随机分布逐步收敛到最优区域的过程。第一阶段遗传算法快速淘汰明显劣质的解，使种群整体适应度快速提升。交叉操作在优良解之间进行信息交换，产生可能比父代更优的子代。变异操作则负责探索新区域，防止种群过早收敛。第一阶段结束时，种群已聚集在少数几个有希望的区域周围。第二阶段粒子群算法利用粒子间的信息共享，在这些区域内进行协作搜索，快速定位局部最优。惯性权重的递减策略使搜索从探索模式平滑过渡到开发模式。

(3) Savannah反应堆伽马屏蔽优化设计应用

Savannah号核动力船是世界上第一艘商用核动力货轮，其反应堆屏蔽设计是核船舶工程的经典案例。该反应堆产生的一次伽马射线需要被有效屏蔽，以保护船员和货物免受辐射危害。原始屏蔽设计采用多层材料组合，包括钢、水和铅等，总重量较大且占用了宝贵的船舱空间。利用GPPG算法对Savannah反应堆的伽马屏蔽进行优化设计，探索在保持屏蔽效果的前提下减轻重量的可能性，或者在相同重量约束下进一步提升屏蔽性能。

优化模型的建立需要准确描述反应堆源项和屏蔽几何。根据公开文献中Savannah反应堆的技术参数，确定了伽马源的强度、能谱和空间分布。屏蔽几何简化为球壳模型，保留了原设计的多层结构特征。优化变量包括五层屏蔽的材料选择和厚度，材料候选集在原设计使用的材料基础上增加了几种新型屏蔽材料如含硼钢和重混凝土。约束条件要求屏蔽体外表面的剂量率不超过国际辐射防护标准规定的限值，同时各层厚度需满足工艺可行性约束。

按照三种优化方案分别运行GPPG算法。方案一侧重重量优化的结果显示，优化后的屏蔽方案相比原设计可减重约百分之十五，同时剂量率仍满足标准要求。减重主要通过材料替换实现，用含硼钢部分替代纯铅可以在保持屏蔽效果的同时降低密度。方案二侧重屏蔽效果的结果显示，在维持原设计重量的条件下，优化方案的外表面剂量率可降低约百分之三十。改进主要来自各层厚度比例的优化调整，使不同能量的伽马射线都能被有效衰减。方案三权衡两目标的结果提供了一组帕累托最优解，决策者可以根据实际需求在重量和屏蔽效果之间进行选择。

import numpy as np import random from dataclasses import dataclass from typing import List, Tuple @dataclass class ShieldMaterial: name: str density: float mu_coeff: float cost: float MATERIALS = { 0: ShieldMaterial("Lead", 11.35, 0.8, 10.0), 1: ShieldMaterial("Steel", 7.87, 0.5, 3.0), 2: ShieldMaterial("Concrete", 2.35, 0.15, 1.0), 3: ShieldMaterial("Water", 1.0, 0.07, 0.5), 4: ShieldMaterial("BoronPE", 0.95, 0.12, 5.0), } class PointKernelCalculator: def __init__(self, source_strength=1e10, source_energy=1.0): self.source_strength = source_strength self.source_energy = source_energy def calculate_dose_rate(self, layer_materials: List[int], layer_thicknesses: List[float], distance: float = 100.0) -> float: total_attenuation = 0.0 for mat_idx, thickness in zip(layer_materials, layer_thicknesses): material = MATERIALS[mat_idx] total_attenuation += material.mu_coeff * material.density * thickness buildup_factor = 1.0 + 0.5 * total_attenuation dose_rate = (self.source_strength * buildup_factor * np.exp(-total_attenuation) / (4 * np.pi * distance**2)) return dose_rate def calculate_weight(self, layer_materials: List[int], layer_thicknesses: List[float], inner_radius: float = 50.0) -> float: total_weight = 0.0 current_radius = inner_radius for mat_idx, thickness in zip(layer_materials, layer_thicknesses): material = MATERIALS[mat_idx] outer_radius = current_radius + thickness volume = (4/3) * np.pi * (outer_radius**3 - current_radius**3) total_weight += volume * material.density current_radius = outer_radius return total_weight class MonteCarloCalculator: def __init__(self, source_strength=1e10, n_particles=10000): self.source_strength = source_strength self.n_particles = n_particles self.pk = PointKernelCalculator(source_strength) def calculate_dose_rate(self, layer_materials: List[int], layer_thicknesses: List[float], distance: float = 100.0) -> float: pk_result = self.pk.calculate_dose_rate(layer_materials, layer_thicknesses, distance) noise = np.random.normal(0, 0.05 * pk_result) return pk_result + noise class GeneticAlgorithm: def __init__(self, n_layers, pop_size=50, mutation_rate=0.1): self.n_layers = n_layers self.pop_size = pop_size self.mutation_rate = mutation_rate self.n_materials = len(MATERIALS) def initialize_population(self) -> List[Tuple[List[int], List[float]]]: population = [] for _ in range(self.pop_size): materials = [random.randint(0, self.n_materials - 1) for _ in range(self.n_layers)] thicknesses = [random.uniform(1.0, 20.0) for _ in range(self.n_layers)] population.append((materials, thicknesses)) return population def crossover(self, parent1, parent2): point = random.randint(1, self.n_layers - 1) child1_mat = parent1[0][:point] + parent2[0][point:] child1_thick = parent1[1][:point] + parent2[1][point:] child2_mat = parent2[0][:point] + parent1[0][point:] child2_thick = parent2[1][:point] + parent1[1][point:] return (child1_mat, child1_thick), (child2_mat, child2_thick) def mutate(self, individual): materials, thicknesses = list(individual[0]), list(individual[1]) for i in range(self.n_layers): if random.random() < self.mutation_rate: materials[i] = random.randint(0, self.n_materials - 1) if random.random() < self.mutation_rate: thicknesses[i] = max(1.0, thicknesses[i] + random.gauss(0, 2.0)) return (materials, thicknesses) class ParticleSwarmOptimizer: def __init__(self, n_layers, pop_size=30, w=0.7, c1=1.5, c2=1.5): self.n_layers = n_layers self.pop_size = pop_size self.w = w self.c1 = c1 self.c2 = c2 self.n_materials = len(MATERIALS) def initialize_from_ga(self, ga_solutions): self.particles = [] self.velocities = [] for sol in ga_solutions[:self.pop_size]: materials = list(sol[0]) thicknesses = list(sol[1]) self.particles.append((materials, thicknesses)) vel_thick = [0.0] * self.n_layers self.velocities.append(vel_thick) self.pbest = [p for p in self.particles] self.gbest = self.particles[0] def update_particle(self, idx): materials, thicknesses = self.particles[idx] velocity = self.velocities[idx] new_thicknesses = [] new_velocity = [] for i in range(self.n_layers): r1, r2 = random.random(), random.random() new_vel = (self.w * velocity[i] + self.c1 * r1 * (self.pbest[idx][1][i] - thicknesses[i]) + self.c2 * r2 * (self.gbest[1][i] - thicknesses[i])) new_thick = max(1.0, min(50.0, thicknesses[i] + new_vel)) new_velocity.append(new_vel) new_thicknesses.append(new_thick) new_materials = list(materials) for i in range(self.n_layers): if random.random() < 0.1: if random.random() < 0.5: new_materials[i] = self.pbest[idx][0][i] else: new_materials[i] = self.gbest[0][i] self.particles[idx] = (new_materials, new_thicknesses) self.velocities[idx] = new_velocity class GPPGOptimizer: def __init__(self, n_layers=5, dose_limit=1.0, weight_objective=True): self.n_layers = n_layers self.dose_limit = dose_limit self.weight_objective = weight_objective self.pk_calc = PointKernelCalculator() self.mc_calc = MonteCarloCalculator() def fitness_function(self, individual, use_mc=False): materials, thicknesses = individual calc = self.mc_calc if use_mc else self.pk_calc dose_rate = calc.calculate_dose_rate(materials, thicknesses) weight = self.pk_calc.calculate_weight(materials, thicknesses) if dose_rate > self.dose_limit: penalty = 1e6 * (dose_rate - self.dose_limit) else: penalty = 0 if self.weight_objective: return weight + penalty else: return dose_rate + penalty def run_ga_phase(self, generations=100): ga = GeneticAlgorithm(self.n_layers) population = ga.initialize_population() fitness = [self.fitness_function(ind, use_mc=False) for ind in population] for gen in range(generations): sorted_indices = np.argsort(fitness) elite = [population[i] for i in sorted_indices[:10]] new_population = elite.copy() while len(new_population) < ga.pop_size: p1, p2 = random.choices(elite, k=2) c1, c2 = ga.crossover(p1, p2) new_population.append(ga.mutate(c1)) if len(new_population) < ga.pop_size: new_population.append(ga.mutate(c2)) population = new_population fitness = [self.fitness_function(ind, use_mc=False) for ind in population] sorted_indices = np.argsort(fitness) return [population[i] for i in sorted_indices[:30]] def run_pso_phase(self, initial_solutions, iterations=50): pso = ParticleSwarmOptimizer(self.n_layers) pso.initialize_from_ga(initial_solutions) fitness = [self.fitness_function(p, use_mc=True) for p in pso.particles] pbest_fitness = fitness.copy() gbest_idx = np.argmin(fitness) pso.gbest = pso.particles[gbest_idx] gbest_fitness = fitness[gbest_idx] for iteration in range(iterations): for i in range(pso.pop_size): pso.update_particle(i) fit = self.fitness_function(pso.particles[i], use_mc=True) fitness[i] = fit if fit < pbest_fitness[i]: pbest_fitness[i] = fit pso.pbest[i] = pso.particles[i] if fit < gbest_fitness: gbest_fitness = fit pso.gbest = pso.particles[i] pso.w *= 0.99 return pso.gbest, gbest_fitness def optimize(self): print("Phase 1: GA + Point Kernel") ga_solutions = self.run_ga_phase(generations=100) print(f"GA phase complete, best solutions found: {len(ga_solutions)}") print("Phase 2: PSO + Monte Carlo") best_solution, best_fitness = self.run_pso_phase(ga_solutions, iterations=50) materials, thicknesses = best_solution print(f"\nOptimal Shield Design:") for i in range(self.n_layers): print(f"Layer {i+1}: {MATERIALS[materials[i]].name}, {thicknesses[i]:.2f} cm") dose = self.mc_calc.calculate_dose_rate(materials, thicknesses) weight = self.pk_calc.calculate_weight(materials, thicknesses) print(f"\nDose Rate: {dose:.4e}") print(f"Total Weight: {weight:.2f} kg") return best_solution, best_fitness if __name__ == "__main__": optimizer = GPPGOptimizer(n_layers=5, dose_limit=1.0, weight_objective=True) solution, fitness = optimizer.optimize()

成品代码50-200，定制300起，可以直接沟通

👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇