GSA 2009算法Python实战5大基准函数测试与动态可视化分析引力搜索算法核心思想解析引力搜索算法Gravitational Search Algorithm, GSA是一种受牛顿万有引力定律启发的群体智能优化算法。想象一下宇宙中的天体运动——质量大的物体会对周围物体产生更强的引力引导它们向自己靠近。GSA正是模拟了这一自然现象将优化问题的解空间视为一个引力场每个候选解看作具有特定质量的粒子。算法通过以下物理量建立数学模型质量与解的适应度值正相关优秀解具有更大质量引力粒子间相互吸引的力与距离平方成反比加速度决定粒子位置更新的方向和步长# 基本物理量计算公式示例 def calculate_force(mass1, mass2, distance, G): return G * (mass1 * mass2) / (distance**2 epsilon) # 添加极小值epsilon避免除零 def calculate_acceleration(force, mass): return force / massPython实现完整架构设计1. 算法参数初始化GSA实现需要合理设置以下关键参数参数名称典型值范围作用说明种群规模20-50解的数量影响搜索广度最大迭代100-500控制算法运行时间G0初始值50-100引力常数初始值alpha衰减系数10-20控制引力衰减速度epsilon1e-10防止除零的小常数class GSA: def __init__(self, dim, pop_size30, max_iter100, G0100, alpha20): self.dim dim # 问题维度 self.pop_size pop_size self.max_iter max_iter self.G0 G0 self.alpha alpha self.epsilon 1e-10 self.agents np.random.uniform(-1, 1, (pop_size, dim)) self.velocity np.zeros((pop_size, dim))2. 核心迭代流程实现算法主循环包含四个关键阶段适应度评估与质量计算引力与加速度计算速度与位置更新引力常数衰减def run(self, func): best_fitness [] for t in range(self.max_iter): # 1. 计算适应度和质量 fitness np.array([func(x) for x in self.agents]) mass self._calculate_mass(fitness) # 2. 计算引力和加速度 G self.G0 * np.exp(-self.alpha * t / self.max_iter) acceleration self._calculate_acceleration(mass, G) # 3. 更新速度和位置 self.velocity np.random.rand(self.pop_size, 1) * self.velocity acceleration self.agents self.velocity # 记录最佳适应度 best_fitness.append(np.min(fitness)) return best_fitness5大经典基准函数测试我们选取以下具有不同特性的测试函数验证算法性能Sphere函数单峰对称函数测试基础收敛性能f_1(x) \sum_{i1}^n x_i^2Rastrigin函数多峰函数测试跳出局部最优能力f_2(x) 10n \sum_{i1}^n [x_i^2 - 10\cos(2\pi x_i)]Ackley函数非线性多峰函数具有陡峭谷地f_3(x) -20\exp(-0.2\sqrt{\frac{1}{n}\sum_{i1}^n x_i^2}) - \exp(\frac{1}{n}\sum_{i1}^n \cos(2\pi x_i)) 20 eRosenbrock函数非凸病态函数测试长距离探索能力f_4(x) \sum_{i1}^{n-1} [100(x_{i1} - x_i^2)^2 (1-x_i)^2]Griewank函数高度多模态但全局最优易找测试局部开发能力f_5(x) 1 \frac{1}{4000}\sum_{i1}^n x_i^2 - \prod_{i1}^n \cos(\frac{x_i}{\sqrt{i}})基准函数Python实现def sphere(x): return np.sum(x**2) def rastrigin(x): return 10*len(x) np.sum(x**2 - 10*np.cos(2*np.pi*x)) def ackley(x): part1 -0.2 * np.sqrt(np.mean(x**2)) part2 np.mean(np.cos(2*np.pi*x)) return -20*np.exp(part1) - np.exp(part2) 20 np.e def rosenbrock(x): return sum(100*(x[1:]-x[:-1]**2)**2 (1-x[:-1])**2) def griewank(x): return 1 np.sum(x**2)/4000 - np.prod(np.cos(x/np.sqrt(np.arange(1,len(x)1))))动态可视化分析系统1. 收敛曲线绘制使用Matplotlib创建动态收敛曲线直观展示算法在各测试函数上的表现import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation def plot_convergence(best_fitness, func_name): fig, ax plt.subplots(figsize(10,6)) ax.set_title(fConvergence Curve on {func_name} Function) ax.set_xlabel(Iteration) ax.set_ylabel(Best Fitness) ax.set_yscale(log) line, ax.plot([], [], b-, lw2) current_iter ax.text(0.02, 0.95, , transformax.transAxes) def init(): line.set_data([], []) return line, current_iter def update(frame): xdata range(frame1) ydata best_fitness[:frame1] line.set_data(xdata, ydata) ax.relim() ax.autoscale_view() current_iter.set_text(fIteration: {frame1}/{len(best_fitness)}) return line, current_iter ani FuncAnimation(fig, update, frameslen(best_fitness), init_funcinit, blitTrue, interval100) plt.show() return ani2. 参数敏感性分析考察引力常数G0和alpha对算法性能的影响def parameter_sensitivity_analysis(func, dim2): G0_values [50, 100, 150] alpha_values [10, 20, 30] fig, axes plt.subplots(len(G0_values), len(alpha_values), figsize(15,10), shareyTrue) for i, G0 in enumerate(G0_values): for j, alpha in enumerate(alpha_values): gsa GSA(dimdim, G0G0, alphaalpha) fitness gsa.run(func) ax axes[i,j] ax.plot(fitness) ax.set_title(fG0{G0}, alpha{alpha}) ax.set_yscale(log) plt.tight_layout() plt.show()性能优化关键技巧1. 自适应参数调整策略原始GSA的固定衰减模式可能不适应所有问题改进方案def adaptive_G(t, max_iter): 自适应引力常数调整 initial_G 100 final_G 1e-5 # 非线性衰减 return initial_G * (1 - t/max_iter)**2 final_G2. 精英保留策略防止优秀个体在迭代过程中被破坏def run_with_elitism(self, func): best_agent None best_fitness float(inf) for t in range(self.max_iter): # ...原有计算流程... # 精英保留 current_best_idx np.argmin(fitness) if fitness[current_best_idx] best_fitness: best_fitness fitness[current_best_idx] best_agent self.agents[current_best_idx].copy() # 用精英替换最差个体 worst_idx np.argmax(fitness) self.agents[worst_idx] best_agent3. 混合局部搜索在后期引入局部搜索增强开发能力def local_search(self, agent, func, radius0.1, attempts10): best agent.copy() best_fitness func(best) for _ in range(attempts): candidate best radius * np.random.uniform(-1,1,len(agent)) candidate_fitness func(candidate) if candidate_fitness best_fitness: best candidate best_fitness candidate_fitness return best完整工程实现建议模块化设计/gsa_optimizer │── core/ │ ├── gsa.py # 主算法实现 │ ├── benchmark.py # 测试函数集合 │── visualization/ │ ├── convergence.py │ ├── surface_plot.py │── examples/ # 使用示例 │── tests/ # 单元测试性能优化技巧使用Numpy向量化运算替代循环对大规模问题实现并行化评估添加边界处理机制防止解越界扩展接口设计class BaseOptimizer: def __init__(self, dim, **kwargs): pass def optimize(self, func): raise NotImplementedError class GSA(BaseOptimizer): # 实现具体算法典型问题解决方案当遇到算法早熟收敛时可尝试以下调整增加种群多样性def add_diversity(self, threshold0.1): std np.std(self.agents, axis0) if np.mean(std) threshold: # 随机重置部分个体 idx np.random.randint(0, self.pop_size) self.agents[idx] np.random.uniform(-1, 1, self.dim)动态调整探索开发平衡def dynamic_exploration(self, t): # 早期侧重探索后期侧重开发 exploration_factor 1 - (t/self.max_iter)**2 return exploration_factor混合其他优化算法优点def hybrid_pso_update(self, t): # 结合PSO的速度更新公式 w 0.9 - 0.5*t/self.max_iter # 惯性权重线性递减 c1, c2 1.5, 1.5 # 学习因子 r1, r2 np.random.rand(2) personal_best ... # 个体历史最优 global_best ... # 全局历史最优 self.velocity w*self.velocity c1*r1*(personal_best - self.agents) c2*r2*(global_best - self.agents)