GSA 2009 算法 Python 实现：5个基准函数测试与收敛曲线可视化

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 / mass

Python实现完整架构设计

1. 算法参数初始化

GSA实现需要合理设置以下关键参数：

class GSA: def __init__(self, dim, pop_size=30, max_iter=100, G0=100, alpha=20): 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. 核心迭代流程实现

算法主循环包含四个关键阶段：

1. 适应度评估与质量计算
2. 引力与加速度计算
3. 速度与位置更新
4. 引力常数衰减

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_fitness

5大经典基准函数测试

我们选取以下具有不同特性的测试函数验证算法性能：

1. Sphere函数：单峰对称函数，测试基础收敛性能

   f_1(x) = \sum_{i=1}^n x_i^2

2. Rastrigin函数：多峰函数，测试跳出局部最优能力

   f_2(x) = 10n + \sum_{i=1}^n [x_i^2 - 10\cos(2\pi x_i)]

3. Ackley函数：非线性多峰函数，具有陡峭谷地

   f_3(x) = -20\exp(-0.2\sqrt{\frac{1}{n}\sum_{i=1}^n x_i^2}) - \exp(\frac{1}{n}\sum_{i=1}^n \cos(2\pi x_i)) + 20 + e

4. Rosenbrock函数：非凸病态函数，测试长距离探索能力

   f_4(x) = \sum_{i=1}^{n-1} [100(x_{i+1} - x_i^2)^2 + (1-x_i)^2]

5. Griewank函数：高度多模态但全局最优易找，测试局部开发能力

   f_5(x) = 1 + \frac{1}{4000}\sum_{i=1}^n x_i^2 - \prod_{i=1}^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(f'Convergence Curve on {func_name} Function') ax.set_xlabel('Iteration') ax.set_ylabel('Best Fitness') ax.set_yscale('log') line, = ax.plot([], [], 'b-', lw=2) current_iter = ax.text(0.02, 0.95, '', transform=ax.transAxes) def init(): line.set_data([], []) return line, current_iter def update(frame): xdata = range(frame+1) ydata = best_fitness[:frame+1] line.set_data(xdata, ydata) ax.relim() ax.autoscale_view() current_iter.set_text(f'Iteration: {frame+1}/{len(best_fitness)}') return line, current_iter ani = FuncAnimation(fig, update, frames=len(best_fitness), init_func=init, blit=True, interval=100) plt.show() return ani

2. 参数敏感性分析

考察引力常数G0和alpha对算法性能的影响：

def parameter_sensitivity_analysis(func, dim=2): G0_values = [50, 100, 150] alpha_values = [10, 20, 30] fig, axes = plt.subplots(len(G0_values), len(alpha_values), figsize=(15,10), sharey=True) for i, G0 in enumerate(G0_values): for j, alpha in enumerate(alpha_values): gsa = GSA(dim=dim, G0=G0, alpha=alpha) fitness = gsa.run(func) ax = axes[i,j] ax.plot(fitness) ax.set_title(f'G0={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_G

2. 精英保留策略

防止优秀个体在迭代过程中被破坏：

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_agent

3. 混合局部搜索

在后期引入局部搜索增强开发能力：

def local_search(self, agent, func, radius=0.1, attempts=10): 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

完整工程实现建议

1. 模块化设计：

   /gsa_optimizer │── core/ │ ├── gsa.py # 主算法实现 │ ├── benchmark.py # 测试函数集合 │── visualization/ │ ├── convergence.py │ ├── surface_plot.py │── examples/ # 使用示例 │── tests/ # 单元测试

2. 性能优化技巧：

   • 使用Numpy向量化运算替代循环
   • 对大规模问题实现并行化评估
   • 添加边界处理机制防止解越界

3. 扩展接口设计：

   class BaseOptimizer: def __init__(self, dim, **kwargs): pass def optimize(self, func): raise NotImplementedError class GSA(BaseOptimizer): # 实现具体算法

典型问题解决方案

当遇到算法早熟收敛时，可尝试以下调整：

1. 增加种群多样性：

   def add_diversity(self, threshold=0.1): std = np.std(self.agents, axis=0) if np.mean(std) < threshold: # 随机重置部分个体 idx = np.random.randint(0, self.pop_size) self.agents[idx] = np.random.uniform(-1, 1, self.dim)

2. 动态调整探索开发平衡：

   def dynamic_exploration(self, t): # 早期侧重探索，后期侧重开发 exploration_factor = 1 - (t/self.max_iter)**2 return exploration_factor

3. 混合其他优化算法优点：

   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)