算法概述与环境配置

淘金优化算法（Gold Rush Optimizer, GRO）是一种于2023年提出的新型元启发式算法。该算法通过模拟淘金者的勘探行为，将搜索过程抽象为迁移、开采和协作三个核心阶段，在解决高维复杂优化问题时展现出优异的寻优能力。为了在Python中高效实现该算法，我们需要预先配置以下科学计算与可视化库：

pip install numpy matplotlib tqdm joblib

• NumPy：提供高效的矩阵与向量运算支持，是实现群体智能算法向量化计算的核心。
• Matplotlib：用于构建适应度收敛曲线及搜索空间的动态演化图。
• tqdm：在终端输出迭代进度条，优化长时间运算的交互体验。
• Joblib：用于后续进阶章节中的种群适应度并行评估。

GRO算法核心架构设计

我们首先定义算法的基础类结构。为了提升计算效率，代码设计中将尽量采用NumPy的广播机制和向量化操作，替代低效的Python原生for循环。

import numpy as np
from tqdm import trange

class GROAlgorithm:
    def __init__(self, objective_fn, dimensions=2, swarm_size=40, max_epochs=150, bounds=None):
        self.objective_fn = objective_fn
        self.dimensions = dimensions
        self.swarm_size = swarm_size
        self.max_epochs = max_epochs
        
        # 默认搜索边界设置
        if bounds is None:
            self.lower_bound = np.full(dimensions, -10.0)
            self.upper_bound = np.full(dimensions, 10.0)
        else:
            self.lower_bound = np.array(bounds[0])
            self.upper_bound = np.array(bounds[1])
            
        self.control_param = 1.0  # 核心控制参数 le
        
    def _init_swarm(self):
        """生成初始种群及其适应度"""
        self.positions = np.random.uniform(
            self.lower_bound, self.upper_bound, (self.swarm_size, self.dimensions)
        )
        self.fitness = np.apply_along_axis(self.objective_fn, 1, self.positions)
        self.best_global_idx = np.argmin(self.fitness)
        self.best_global_pos = self.positions[self.best_global_idx].copy()
        self.best_global_fit = self.fitness[self.best_global_idx]

    def _update_control_params(self, epoch):
        """动态调整非线性控制因子"""
        decay_ratio = (self.max_epochs - epoch) / (self.max_epochs - 1)
        self.control_param = (decay_ratio ** 2) * (2 - 1 / self.max_epochs) + (1 / self.max_epochs)
        
    def _enforce_boundaries(self, positions):
        """边界截断处理"""
        return np.clip(positions, self.lower_bound, self.upper_bound)

    def _migration_step(self):
        """阶段一：全局迁移（向最优区域靠拢）"""
        pass

    def _mining_step(self):
        """阶段二：局部开采（精细搜索）"""
        pass

    def _collaboration_step(self):
        """阶段三：信息协作（种群交互）"""
        pass

    def run(self):
        """执行主优化循环"""
        self._init_swarm()
        history = []
        
        for epoch in trange(self.max_epochs, desc="GRO Optimization"):
            self._update_control_params(epoch)
            
            self._migration_step()
            self._mining_step()
            self._collaboration_step()
            
            # 更新全局最优记录
            current_best_idx = np.argmin(self.fitness)
            if self.fitness[current_best_idx] < self.best_global_fit:
                self.best_global_fit = self.fitness[current_best_idx]
                self.best_global_pos = self.positions[current_best_idx].copy()
                
            history.append(self.best_global_fit)
            
        return history, self.best_global_pos

三大核心搜索机制实现

1. 迁移机制 (Migration)

迁移阶段旨在利用当前全局最优解引导种群向高潜力区域移动。通过引入随机扰动因子，平衡算法的探索与开发能力。

def _migration_step(self):
    # 生成随机系数矩阵
    rand_matrix_1 = np.random.rand(self.swarm_size, self.dimensions)
    rand_matrix_2 = np.random.rand(self.swarm_size, self.dimensions)
    
    # 计算移动步长与方向因子
    step_factor = 1.0 + self.control_param * (rand_matrix_1 - 0.5)
    direction_factor = 2.0 * rand_matrix_2
    
    # 计算指向全局最优的位移向量
    displacement = direction_factor * (self.best_global_pos - self.positions)
    new_positions = self.positions + step_factor * displacement
    
    new_positions = self._enforce_boundaries(new_positions)
    new_fitness = np.apply_along_axis(self.objective_fn, 1, new_positions)
    
    # 贪婪选择策略
    mask = new_fitness < self.fitness
    self.positions[mask] = new_positions[mask]
    self.fitness[mask] = new_fitness[mask]

2. 开采机制 (Gold Mining)

开采阶段模拟淘金者在特定矿脉附近的深度挖掘。此过程不依赖全局最优信息，而是通过个体间的随机配对来维持种群多样性，防止算法陷入局部极值。

def _mining_step(self):
    # 为每个个体随机分配一个交互伙伴（排除自身）
    partner_indices = np.random.randint(0, self.swarm_size, self.swarm_size)
    for i in range(self.swarm_size):
        while partner_indices[i] == i:
            partner_indices[i] = np.random.randint(0, self.swarm_size)
            
    partner_positions = self.positions[partner_indices]
    
    # 计算开采扰动
    rand_vec = np.random.rand(self.swarm_size, self.dimensions)
    mining_factor = 2.0 * self.control_param * rand_vec - self.control_param
    
    relative_dist = self.positions - partner_positions
    new_positions = partner_positions + mining_factor * relative_dist
    
    new_positions = self._enforce_boundaries(new_positions)
    new_fitness = np.apply_along_axis(self.objective_fn, 1, new_positions)
    
    mask = new_fitness < self.fitness
    self.positions[mask] = new_positions[mask]
    self.fitness[mask] = new_fitness[mask]

3. 协作机制 (Collaboration)

协作阶段通过随机抽取两个其他个体，利用它们之间的空间差异来更新当前个体的位置。这种无中心化的信息交换模式在处理多峰复杂地形时具有显著优势。

def _collaboration_step(self):
    # 随机抽取两个不同的邻居个体
    neighbor_1 = np.random.randint(0, self.swarm_size, self.swarm_size)
    neighbor_2 = np.random.randint(0, self.swarm_size, self.swarm_size)
    
    # 确保邻居不与当前个体重复，且互不相同
    for i in range(self.swarm_size):
        while neighbor_1[i] == i:
            neighbor_1[i] = np.random.randint(0, self.swarm_size)
        while neighbor_2[i] == i or neighbor_2[i] == neighbor_1[i]:
            neighbor_2[i] = np.random.randint(0, self.swarm_size)
            
    pos_1 = self.positions[neighbor_1]
    pos_2 = self.positions[neighbor_2]
    
    # 基于邻居差异进行位置更新
    random_scale = np.random.rand(self.swarm_size, self.dimensions)
    new_positions = self.positions + random_scale * (pos_2 - pos_1)
    
    new_positions = self._enforce_boundaries(new_positions)
    new_fitness = np.apply_along_axis(self.objective_fn, 1, new_positions)
    
    mask = new_fitness < self.fitness
    self.positions[mask] = new_positions[mask]
    self.fitness[mask] = new_fitness[mask]

算法性能评估与可视化

基准测试函数定义

为了验证算法的鲁棒性，我们选取三种具有不同地形特征的基准函数进行测试：

函数名称	数学表达式	地形特征
Sphere	f(x) = Σx²	单峰、连续凸函数，测试基础收敛速度
Rastrigin	f(x) = 10d + Σ[x² - 10cos(2πx)]	多峰、存在大量局部极小值，测试逃逸能力
Ackley	包含指数与余弦项的复杂组合	全局平坦且中心存在极小值，测试局部开发精度

def sphere_func(x):
    return np.sum(np.square(x))

def rastrigin_func(x):
    return 10 * len(x) + np.sum(np.square(x) - 10 * np.cos(2 * np.pi * x))

def ackley_func(x):
    part1 = -0.2 * np.sqrt(np.mean(np.square(x)))
    part2 = np.mean(np.cos(2 * np.pi * x))
    return -20 * np.exp(part1) - np.exp(part2) + 20 + np.e

# 执行批量测试
test_suite = {
    'Sphere': sphere_func,
    'Rastrigin': rastrigin_func,
    'Ackley': ackley_func
}

convergence_logs = {}
for name, func in test_suite.items():
    optimizer = GROAlgorithm(func, dimensions=2, swarm_size=50, max_epochs=100, bounds=[[-5, -5], [5, 5]])
    history, _ = optimizer.run()
    convergence_logs[name] = history

收敛曲线绘制

通过Matplotlib将各测试函数的适应度下降过程进行对数级可视化：

import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
for name, log in convergence_logs.items():
    plt.plot(log, label=name, linewidth=2)

plt.xlabel('Epoch')
plt.ylabel('Log Fitness')
plt.title('GRO Convergence Analysis on Benchmark Functions')
plt.yscale('log')
plt.legend()
plt.grid(True, which="both", ls="--", alpha=0.5)
plt.tight_layout()
plt.show()

从输出图表中可以观察到：Sphere函数在初期呈现极快的下降斜率；Rastrigin函数由于地形崎岖，收敛曲线存在明显的阶梯状停滞，但最终仍能逼近全局最优；Ackley函数则展现了算法在平坦区域与陡峭极值区之间的良好过渡能力。

进阶优化策略

1. 控制参数非线性衰减

原始GRO的控制参数呈二次曲线衰减。为了在搜索后期提供更平滑的局部挖掘能力，可将其替换为余弦退火机制：

def _update_control_params(self, epoch):
    # 余弦退火策略
    progress = epoch / self.max_epochs
    self.control_param = 0.5 * (1 + np.cos(np.pi * progress))

2. 适应度并行评估

当目标函数计算成本较高（如有限元分析或深度学习模型推理）时，串行评估将成为性能瓶颈。利用joblib可实现种群级别的并行计算：

from joblib import Parallel, delayed

class ParallelGRO(GROAlgorithm):
    def __init__(self, *args, n_cores=-1, **kwargs):
        super().__init__(*args, **kwargs)
        self.n_cores = n_cores

    def _parallel_evaluate(self, positions):
        return np.array(Parallel(n_jobs=self.n_cores)(
            delayed(self.objective_fn)(ind) for ind in positions
        ))

    def _init_swarm(self):
        self.positions = np.random.uniform(
            self.lower_bound, self.upper_bound, (self.swarm_size, self.dimensions)
        )
        self.fitness = self._parallel_evaluate(self.positions)
        self.best_global_idx = np.argmin(self.fitness)
        self.best_global_pos = self.positions[self.best_global_idx].copy()
        self.best_global_fit = self.fitness[self.best_global_idx]

3. 约束条件处理

针对带有不等式或等式约束的实际工程问题，可通过外部罚函数法将约束优化转化为无约束优化：

def apply_penalty(original_fn, constraint_fns, penalty_weight=1e6):
    def penalized_fn(x):
        base_val = original_fn(x)
        penalty = 0.0
        for c_fn in constraint_fns:
            violation = c_fn(x)
            if violation > 0:
                penalty += penalty_weight * (violation ** 2)
        return base_val + penalty
    return penalized_fn

# 使用示例：
# constrained_fn = apply_penalty(my_objective, [constraint_1, constraint_2])
# optimizer = GROAlgorithm(constrained_fn, ...)