主题039：差分进化算法（Differential Evolution）

一、引言

1.1 算法背景

差分进化算法（Differential Evolution, DE）是由Rainer Storn和Kenneth Price于1995年提出的一种基于种群的随机优化算法。最初设计用于解决切比雪夫多项式拟合问题，后来被发现是一种高效的全局优化方法，特别适用于连续优化问题。

1.2 核心思想

DE算法的核心思想可以概括为：

• 差分变异：利用种群中个体间的差分向量进行变异
• 交叉重组：通过交叉操作产生试验个体
• 贪婪选择：一对一竞争，保留更优个体

1.3 算法优势

差分进化算法具有以下优势：

• 简单易实现：只有变异、交叉、选择三个基本操作
• 参数少：主要参数只有F（缩放因子）和CR（交叉概率）
• 全局搜索能力强：差分变异具有良好的探索能力
• 收敛速度快：贪婪选择机制保证了收敛速度
• 适用于连续优化：特别适合处理实数优化问题

---

二、数学模型与理论基础

2.1 基本DE算法

2.1.1 变异操作（Mutation）

对于目标个体$x_i$，变异向量$v_i$通过以下方式产生：

DE/rand/1：
$$v_i=x_{r1}+F\cdot (x_{r2}-x_{r3})$$

DE/best/1：
$$v_i=x_{best}+F\cdot (x_{r1}-x_{r2})$$

DE/current-to-best/1：
$$v_i=x_i+F\cdot (x_{best}-x_i)+F\cdot (x_{r1}-x_{r2})$$

其中：

• $x_{r1},x_{r2},x_{r3}$：从种群中随机选择的三个不同个体
• $x_{best}$：当前最优个体
• $F\in [0,2]$：缩放因子，控制差分向量的幅度

2.1.2 交叉操作（Crossover）

二项式交叉（Binomial Crossover）：

$$u_{i,j}=\{\begin{array}{ll}{v}_{i,j}& \text{if }ran{d}_j\le CR\text{ or }j=j_{rand}\\ x_{i,j}& \text{otherwise}\end{array}$$

其中：

• $CR\in [0,1]$：交叉概率
• $j_{rand}$：随机选择的维度索引，确保至少有一个维度来自变异向量
• $ran{d}_j$：在[0,1]区间均匀分布的随机数

2.1.3 选择操作（Selection）

贪婪选择：

$$x_i^{(t+1)}=\{\begin{array}{ll}{u}_i& \text{if }f(u_i)\le f(x_i^{(t)})\\ x_i^{(t)}& \text{otherwise}\end{array}$$

2.2 变异策略分类

DE算法的变异策略通常表示为：DE/x/y/z

• x：基向量的选择方式

  • rand：随机选择
  • best：选择最优个体
  • current-to-best：当前个体指向最优个体

• y：差分向量的数量（1或2）

• z：交叉方式（bin：二项式，exp：指数）

常用策略：

策略	公式	特点
DE/rand/1	$v=r_1+F\cdot (r_2-r_3)$	探索性强
DE/rand/2	$v=r_1+F\cdot (r_2-r_3+r_4-r_5)$	更强探索性
DE/best/1	$v=best+F\cdot (r_1-r_2)$	开发性强
DE/best/2	$v=best+F\cdot (r_1-r_2+r_3-r_4)$	更强开发性
DE/current-to-best/1	$v=x+F\cdot (best-x)+F\cdot (r_1-r_2)$	平衡探索开发

2.3 参数分析

2.3.1 缩放因子F

• 作用：控制差分向量的缩放幅度
• 推荐范围：$[0.5,1.0]$
• 影响：
  • F太小：搜索步长小，收敛慢，易陷入局部最优
  • F太大：搜索步长大，可能跳过最优解

2.3.2 交叉概率CR

• 作用：控制从变异向量继承的维度比例
• 推荐范围：$[0.8,1.0]$
• 影响：
  • CR太小：继承维度少，收敛慢
  • CR太大：继承维度多，可能丢失种群多样性

2.3.3 种群大小NP

• 推荐值：$NP=5\times D$到$10\times D$，$D$为维度
• 影响：
  • NP太小：种群多样性不足
  • NP太大：计算成本高

---

三、Python实现

3.1 完整代码

"""
主题039：差分进化算法
Differential Evolution (DE) for Structural Optimization
"""

import numpy as np
import matplotlib.pyplot as plt
from typing import Callable, List, Tuple, Optional

# 设置后端
plt.switch_backend('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class DE:
    """差分进化算法"""
    
    def __init__(self,
                 n_vars: int,
                 population_size: int = 50,
                 max_iter: int = 100,
                 F: float = 0.8,  # 缩放因子
                 CR: float = 0.9,  # 交叉概率
                 strategy: str = 'rand/1'):
        """
        初始化差分进化算法
        
        参数:
            n_vars: 设计变量数
            population_size: 种群大小
            max_iter: 最大迭代次数
            F: 缩放因子 (通常0.5-1.0)
            CR: 交叉概率 (通常0.8-1.0)
            strategy: 变异策略
        """
        self.n_vars = n_vars
        self.population_size = population_size
        self.max_iter = max_iter
        self.F = F
        self.CR = CR
        self.strategy = strategy
        
        # 变量边界
        self.xl = None
        self.xu = None
        
        # 目标函数和约束
        self.objective_func = None
        self.constraint_func = None
        
        # 种群
        self.population = None
        self.fitness = None
        self.objectives = None
        
        # 最优解
        self.best_individual = None
        self.best_fitness = float('inf')
        self.best_objective = float('inf')
        
        # 历史记录
        self.history = {
            'best_fitness': [],
            'mean_fitness': [],
            'best_objective': [],
            'diversity': [],
            'F_history': [],
            'CR_history': []
        }
    
    def initialize_population(self):
        """初始化种群"""
        self.population = np.random.uniform(
            self.xl, self.xu, (self.population_size, self.n_vars)
        )
        self.fitness = np.zeros(self.population_size)
        self.objectives = np.zeros(self.population_size)
        
        # 评估初始种群
        for i in range(self.population_size):
            self.objectives[i] = self.objective_func(self.population[i])
            
            if self.constraint_func:
                constraints = self.constraint_func(self.population[i])
                violation = np.sum(np.maximum(0, constraints))
                self.fitness[i] = self.objectives[i] + 1e6 * violation
            else:
                self.fitness[i] = self.objectives[i]
        
        # 找到最优个体
        best_idx = np.argmin(self.fitness)
        self.best_individual = self.population[best_idx].copy()
        self.best_fitness = self.fitness[best_idx]
        self.best_objective = self.objectives[best_idx]
    
    def mutate(self, target_idx: int) -> np.ndarray:
        """
        变异操作
        
        参数:
            target_idx: 目标个体索引
        """
        if self.strategy == 'rand/1':
            # rand/1: v = r1 + F * (r2 - r3)
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2, r3 = np.random.choice(candidates, 3, replace=False)
            mutant = self.population[r1] + self.F * (self.population[r2] - self.population[r3])
            
        elif self.strategy == 'best/1':
            # best/1: v = best + F * (r1 - r2)
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2 = np.random.choice(candidates, 2, replace=False)
            mutant = self.best_individual + self.F * (self.population[r1] - self.population[r2])
            
        elif self.strategy == 'current-to-best/1':
            # current-to-best/1: v = current + F * (best - current) + F * (r1 - r2)
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2 = np.random.choice(candidates, 2, replace=False)
            mutant = (self.population[target_idx] + 
                     self.F * (self.best_individual - self.population[target_idx]) +
                     self.F * (self.population[r1] - self.population[r2]))
        else:
            # 默认 rand/1
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2, r3 = np.random.choice(candidates, 3, replace=False)
            mutant = self.population[r1] + self.F * (self.population[r2] - self.population[r3])
        
        # 边界处理
        mutant = np.clip(mutant, self.xl, self.xu)
        
        return mutant
    
    def crossover(self, target: np.ndarray, mutant: np.ndarray) -> np.ndarray:
        """
        交叉操作（二项式交叉）
        
        参数:
            target: 目标个体
            mutant: 变异个体
        """
        trial = np.zeros(self.n_vars)
        j_rand = np.random.randint(self.n_vars)  # 确保至少有一个维度来自变异个体
        
        for j in range(self.n_vars):
            if np.random.random() < self.CR or j == j_rand:
                trial[j] = mutant[j]
            else:
                trial[j] = target[j]
        
        return trial
    
    def select(self, target_idx: int, trial: np.ndarray):
        """
        选择操作
        
        参数:
            target_idx: 目标个体索引
            trial: 试验个体
        """
        # 评估试验个体
        trial_objective = self.objective_func(trial)
        
        if self.constraint_func:
            constraints = self.constraint_func(trial)
            trial_fitness = trial_objective + 1e6 * np.sum(np.maximum(0, constraints))
        else:
            trial_fitness = trial_objective
        
        # 贪婪选择
        if trial_fitness < self.fitness[target_idx]:
            self.population[target_idx] = trial
            self.fitness[target_idx] = trial_fitness
            self.objectives[target_idx] = trial_objective
            
            # 更新最优
            if trial_fitness < self.best_fitness:
                self.best_individual = trial.copy()
                self.best_fitness = trial_fitness
                self.best_objective = trial_objective
    
    def optimize(self,
                 objective_func: Callable,
                 xl: np.ndarray,
                 xu: np.ndarray,
                 constraint_func: Callable = None) -> np.ndarray:
        """
        执行优化
        
        参数:
            objective_func: 目标函数
            xl, xu: 变量下界和上界
            constraint_func: 约束函数
        """
        self.objective_func = objective_func
        self.constraint_func = constraint_func
        self.xl = xl
        self.xu = xu
        
        # 初始化
        self.initialize_population()
        
        # 迭代优化
        for iteration in range(self.max_iter):
            for i in range(self.population_size):
                # 变异
                mutant = self.mutate(i)
                
                # 交叉
                trial = self.crossover(self.population[i], mutant)
                
                # 选择
                self.select(i, trial)
            
            # 记录历史
            self.history['best_fitness'].append(self.best_fitness)
            self.history['mean_fitness'].append(np.mean(self.fitness))
            self.history['best_objective'].append(self.best_objective)
        
        return self.best_individual

3.2 自适应差分进化（jDE）

class AdaptiveDE(DE):
    """自适应差分进化算法 (jDE)"""
    
    def __init__(self, n_vars: int, population_size: int = 50, max_iter: int = 100,
                 F_init: float = 0.5, CR_init: float = 0.5,
                 tau_F: float = 0.1, tau_CR: float = 0.1,
                 Fl: float = 0.1, Fu: float = 0.9,
                 strategy: str = 'rand/1/bin'):
        """
        初始化自适应DE
        
        参数:
            tau_F, tau_CR: 自适应参数
            Fl, Fu: F的上下界
        """
        super().__init__(n_vars, population_size, max_iter, F_init, CR_init, strategy)
        self.tau_F = tau_F
        self.tau_CR = tau_CR
        self.Fl = Fl
        self.Fu = Fu
        
        # 每个个体的F和CR
        self.F_values = None
        self.CR_values = None
    
    def initialize_population(self):
        """初始化种群和自适应参数"""
        super().initialize_population()
        self.F_values = np.ones(self.population_size) * self.F
        self.CR_values = np.ones(self.population_size) * self.CR
    
    def mutate(self, target_idx: int) -> np.ndarray:
        """自适应变异"""
        # 自适应调整F
        if np.random.random() < self.tau_F:
            self.F_values[target_idx] = self.Fl + np.random.random() * (self.Fu - self.Fl)
        
        self.F = self.F_values[target_idx]
        return super().mutate(target_idx)
    
    def crossover(self, target: np.ndarray, mutant: np.ndarray) -> np.ndarray:
        """自适应交叉"""
        # 自适应调整CR
        target_idx = np.where((self.population == target).all(axis=1))[0][0]
        if np.random.random() < self.tau_CR:
            self.CR_values[target_idx] = np.random.random()
        
        self.CR = self.CR_values[target_idx]
        return super().crossover(target, mutant)

---

四、代码解析

4.1 核心算法流程

差分进化算法流程
│
├─ 1. 初始化
│  ├─ 设置参数（F, CR, NP等）
│  ├─ 随机初始化种群
│  └─ 评估初始种群
│
├─ 2. 迭代优化（直到满足终止条件）
│  ├─ 对种群中每个个体
│  │  ├─ 变异：生成变异向量
│  │  ├─ 交叉：生成试验向量
│  │  └─ 选择：贪婪选择更优个体
│  └─ 更新最优解
│
└─ 3. 输出结果

4.2 关键函数详解

4.2.1 变异操作

def mutate(self, target_idx: int) -> np.ndarray:
    """
    变异操作 - 差分进化的核心
    
    核心逻辑：
    1. 根据策略选择基向量
    2. 添加差分向量（缩放后）
    3. 边界处理
    """
    if self.strategy == 'rand/1':
        # 随机选择3个不同个体
        candidates = [i for i in range(self.population_size) if i != target_idx]
        r1, r2, r3 = np.random.choice(candidates, 3, replace=False)
        
        # 差分变异：基向量 + F * 差分向量
        mutant = self.population[r1] + self.F * (
            self.population[r2] - self.population[r3]
        )
    
    elif self.strategy == 'best/1':
        # 使用最优个体作为基向量
        candidates = [i for i in range(self.population_size) if i != target_idx]
        r1, r2 = np.random.choice(candidates, 2, replace=False)
        
        mutant = self.best_individual + self.F * (
            self.population[r1] - self.population[r2]
        )
    
    # 边界处理：裁剪到可行域
    mutant = np.clip(mutant, self.xl, self.xu)
    
    return mutant

4.2.2 交叉操作

def crossover(self, target: np.ndarray, mutant: np.ndarray) -> np.ndarray:
    """
    二项式交叉操作
    
    核心逻辑：
    1. 对每个维度，以概率CR选择变异向量的值
    2. 确保至少有一个维度来自变异向量
    """
    trial = np.zeros(self.n_vars)
    
    # 随机选择一个维度，强制使用变异向量的值
    j_rand = np.random.randint(self.n_vars)
    
    for j in range(self.n_vars):
        if np.random.random() < self.CR or j == j_rand:
            trial[j] = mutant[j]  # 来自变异向量
        else:
            trial[j] = target[j]  # 来自目标向量
    
    return trial

4.2.3 选择操作

def select(self, target_idx: int, trial: np.ndarray):
    """
    贪婪选择操作
    
    核心逻辑：
    一对一竞争，保留适应度更好的个体
    """
    # 评估试验个体
    trial_objective = self.objective_func(trial)
    trial_fitness = trial_objective  # 简化，实际需考虑约束
    
    # 贪婪选择：如果试验个体更好，则替换
    if trial_fitness < self.fitness[target_idx]:
        self.population[target_idx] = trial
        self.fitness[target_idx] = trial_fitness
        self.objectives[target_idx] = trial_objective
        
        # 更新全局最优
        if trial_fitness < self.best_fitness:
            self.best_individual = trial.copy()
            self.best_fitness = trial_fitness
            self.best_objective = trial_objective

---

五、测试案例与结果分析

5.1 测试函数

5.1.1 Sphere函数

$$f(x)=\sum _{i=1}^n{x}_i^2,x_i\in [-5.12,5.12]$$

• 全局最小值：$f(0)=0$
• 特点：单峰、平滑、易优化

5.1.2 Rastrigin函数

$$f(x)=10n+\sum _{i=1}^n[x_i^2-10\cos (2\pi x_i)]$$

• 全局最小值：$f(0)=0$
• 特点：多峰、具有大量局部最优

5.1.3 Rosenbrock函数

$$f(x)=\sum _{i=1}^{n-1}[100(x_{i+1}-x_i^2{)}^2+(1-x_i{)}^2]$$

• 全局最小值：$f(1)=0$
• 特点：病态、狭窄山谷

5.2 桁架优化案例

问题描述：10杆桁架重量最小化

目标函数：
$$W=\rho \sum _{i=1}^{10}{A}_i{L}_i$$

约束条件：
$$A_i\ge 0.0001{\text{ m}}^2$$

5.3 运行结果

============================================================
差分进化算法 (Differential Evolution)
============================================================

测试1: Sphere函数优化
------------------------------------------------------------
设计变量数: 10
种群大小: 50
迭代次数: 100
缩放因子 F: 0.8
交叉概率 CR: 0.9
变异策略: rand/1
Iteration   0: Best=106.496797, Mean=179.999419, Div=9.4318
Iteration  20: Best=89.834255, Mean=139.772334, Div=8.4454
Iteration  40: Best=85.185135, Mean=126.429153, Div=7.9219
Iteration  60: Best=85.120770, Mean=119.528484, Div=7.6453
Iteration  80: Best=85.120770, Mean=115.833838, Div=7.5245
优化完成!
最优目标值: 1.849224e+00

测试2: Rastrigin函数优化
------------------------------------------------------------
最优目标值: 8.633296

测试3: Rosenbrock函数优化
------------------------------------------------------------
最优目标值: 0.000000

测试5: 桁架重量优化
------------------------------------------------------------
最优重量: 0.0091 吨

5.4 可视化结果

程序生成了以下可视化图片：

1. de_sphere_convergence.png - Sphere函数收敛曲线
2. de_rastrigin_convergence.png - Rastrigin函数收敛曲线
3. de_rosenbrock_convergence.png - Rosenbrock函数收敛曲线
4. de_ackley_convergence.png - Ackley函数收敛曲线
5. de_truss_convergence.png - 桁架优化收敛曲线
6. de_jde_convergence.png - 自适应DE收敛曲线
7. de_strategy_comparison.png - 策略比较图
8. de_parameter_comparison.png - 参数比较图
9. de_comparison.png - 综合对比图

---

六、算法变体与改进

6.1 自适应DE（jDE）

核心思想：每个个体自适应调整F和CR参数

自适应机制：

if rand < τ_F:
    F_i = Fl + rand * (Fu - Fl)
if rand < τ_CR:
    CR_i = rand

优势：

• 无需手动调参
• 适应不同优化阶段
• 提高算法鲁棒性

6.2 自适应DE with External Archive（JADE）

改进点：

• 使用外部存档保存失败解
• 自适应调整F和CR
• 新的变异策略：DE/current-to-pbest/1

6.3 基于成功历史的DE（SHADE）

核心思想：

• 维护F和CR的历史记忆
• 根据成功历史自适应调整参数
• 使用外部存档

---

七、与其他算法的比较

7.1 算法特性对比

特性	DE	GA	PSO	ACO
搜索机制	差分变异	遗传操作	速度更新	信息素引导
参数数量	少（F, CR）	多	中等	多
连续优化	非常适合	适合	非常适合	需离散化
收敛速度	快	中等	快	中等
全局搜索	强	中等	中等	强
实现难度	简单	中等	简单	中等

7.2 适用场景

• DE：连续优化、实数编码问题
• GA：离散优化、组合问题
• PSO：连续优化、神经网络训练
• ACO：组合优化、路径规划

---

八、常见问题与解决方案

8.1 问题1：早熟收敛

现象：算法过早收敛到局部最优。

解决方案：

• 增大F值（增加探索步长）
• 增大种群大小
• 使用rand/1策略（增加探索性）
• 引入重启机制

8.2 问题2：收敛过慢

现象：算法收敛速度太慢。

解决方案：

• 使用best/1策略（增加开发性）
• 减小F值（减小步长）
• 增大CR值（增加继承比例）
• 自适应调整参数

8.3 问题3：参数敏感

现象：算法性能对参数敏感。

解决方案：

• 使用自适应DE（jDE, SHADE）
• 参数敏感性分析
• 元优化方法优化参数

8.4 问题4：边界违反

现象：变异后个体超出边界。

解决方案：

• 裁剪法（Clipping）
• 反射法（Reflection）
• 吸收法（Absorption）
• 重新初始化

---

九、进阶应用

9.1 多目标差分进化

class MultiObjectiveDE(DE):
    """多目标差分进化"""
    
    def __init__(self, n_vars: int, n_objectives: int, ...):
        super().__init__(n_vars, ...)
        self.n_objectives = n_objectives
        self.pareto_front = []
    
    def dominates(self, obj1, obj2):
        """判断obj1是否支配obj2"""
        return all(o1 <= o2 for o1, o2 in zip(obj1, obj2)) and \
               any(o1 < o2 for o1, o2 in zip(obj1, obj2))
    
    def update_pareto_front(self, individual, objectives):
        """更新Pareto前沿"""
        # 实现非支配排序逻辑
        pass

9.2 约束处理

罚函数法：

def constrained_objective(x):
    obj = objective(x)
    violation = constraint_violation(x)
    return obj + penalty * violation

可行性规则：

def constrained_select(fitness_trial, violation_trial, 
                       fitness_target, violation_target):
    if violation_trial == 0 and violation_target == 0:
        return fitness_trial < fitness_target
    elif violation_trial == 0:
        return True
    elif violation_target == 0:
        return False
    else:
        return violation_trial < violation_target

9.3 混合算法

DE + 局部搜索：

def hybrid_de(self):
    # DE全局搜索
    best = self.optimize(...)
    
    # 局部搜索精细优化
    from scipy.optimize import minimize
    result = minimize(
        self.objective_func,
        best,
        method='L-BFGS-B',
        bounds=list(zip(self.xl, self.xu))
    )
    
    return result

---

附录：完整代码文件

• differential_evolution.py - 完整实现代码
• 差分进化算法.md - 本教程文档

"""
主题039：差分进化算法
Differential Evolution (DE) for Structural Optimization
"""

import numpy as np
import matplotlib.pyplot as plt
from typing import Callable, List, Tuple, Optional

# 设置后端
plt.switch_backend('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class DE:
    """差分进化算法"""
    
    def __init__(self,
                 n_vars: int,
                 population_size: int = 50,
                 max_iter: int = 100,
                 F: float = 0.8,  # 缩放因子
                 CR: float = 0.9,  # 交叉概率
                 strategy: str = 'rand/1/bin'):
        """
        初始化差分进化算法
        
        参数:
            n_vars: 设计变量数
            population_size: 种群大小
            max_iter: 最大迭代次数
            F: 缩放因子 (通常0.5-1.0)
            CR: 交叉概率 (通常0.8-1.0)
            strategy: 变异策略
        """
        self.n_vars = n_vars
        self.population_size = population_size
        self.max_iter = max_iter
        self.F = F
        self.CR = CR
        self.strategy = strategy
        
        # 变量边界
        self.xl = None
        self.xu = None
        
        # 目标函数和约束
        self.objective_func = None
        self.constraint_func = None
        
        # 种群
        self.population = None
        self.fitness = None
        self.objectives = None
        
        # 最优解
        self.best_individual = None
        self.best_fitness = float('inf')
        self.best_objective = float('inf')
        
        # 历史记录
        self.history = {
            'best_fitness': [],
            'mean_fitness': [],
            'best_objective': [],
            'diversity': [],
            'F_history': [],
            'CR_history': []
        }
    
    def initialize_population(self):
        """初始化种群"""
        self.population = np.random.uniform(
            self.xl, self.xu, (self.population_size, self.n_vars)
        )
        self.fitness = np.zeros(self.population_size)
        self.objectives = np.zeros(self.population_size)
        
        # 评估初始种群
        for i in range(self.population_size):
            self.objectives[i] = self.objective_func(self.population[i])
            
            if self.constraint_func:
                constraints = self.constraint_func(self.population[i])
                violation = np.sum(np.maximum(0, constraints))
                self.fitness[i] = self.objectives[i] + 1e6 * violation
            else:
                self.fitness[i] = self.objectives[i]
        
        # 找到最优个体
        best_idx = np.argmin(self.fitness)
        self.best_individual = self.population[best_idx].copy()
        self.best_fitness = self.fitness[best_idx]
        self.best_objective = self.objectives[best_idx]
    
    def mutate(self, target_idx: int) -> np.ndarray:
        """
        变异操作
        
        参数:
            target_idx: 目标个体索引
        """
        if self.strategy == 'rand/1':
            # rand/1: v = r1 + F * (r2 - r3)
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2, r3 = np.random.choice(candidates, 3, replace=False)
            mutant = self.population[r1] + self.F * (self.population[r2] - self.population[r3])
            
        elif self.strategy == 'rand/2':
            # rand/2: v = r1 + F * (r2 - r3 + r4 - r5)
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2, r3, r4, r5 = np.random.choice(candidates, 5, replace=False)
            mutant = self.population[r1] + self.F * (
                self.population[r2] - self.population[r3] + 
                self.population[r4] - self.population[r5]
            )
            
        elif self.strategy == 'best/1':
            # best/1: v = best + F * (r1 - r2)
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2 = np.random.choice(candidates, 2, replace=False)
            mutant = self.best_individual + self.F * (self.population[r1] - self.population[r2])
            
        elif self.strategy == 'best/2':
            # best/2: v = best + F * (r1 - r2 + r3 - r4)
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2, r3, r4 = np.random.choice(candidates, 4, replace=False)
            mutant = self.best_individual + self.F * (
                self.population[r1] - self.population[r2] + 
                self.population[r3] - self.population[r4]
            )
            
        elif self.strategy == 'current-to-best/1':
            # current-to-best/1: v = current + F * (best - current) + F * (r1 - r2)
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2 = np.random.choice(candidates, 2, replace=False)
            mutant = (self.population[target_idx] + 
                     self.F * (self.best_individual - self.population[target_idx]) +
                     self.F * (self.population[r1] - self.population[r2]))
            
        elif self.strategy == 'rand-to-best/1':
            # rand-to-best/1: v = r1 + F * (best - r1) + F * (r2 - r3)
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2, r3 = np.random.choice(candidates, 3, replace=False)
            mutant = (self.population[r1] + 
                     self.F * (self.best_individual - self.population[r1]) +
                     self.F * (self.population[r2] - self.population[r3]))
        else:
            # 默认 rand/1
            candidates = [i for i in range(self.population_size) if i != target_idx]
            r1, r2, r3 = np.random.choice(candidates, 3, replace=False)
            mutant = self.population[r1] + self.F * (self.population[r2] - self.population[r3])
        
        # 边界处理
        mutant = np.clip(mutant, self.xl, self.xu)
        
        return mutant
    
    def crossover(self, target: np.ndarray, mutant: np.ndarray) -> np.ndarray:
        """
        交叉操作（二项式交叉）
        
        参数:
            target: 目标个体
            mutant: 变异个体
        """
        trial = np.zeros(self.n_vars)
        j_rand = np.random.randint(self.n_vars)  # 确保至少有一个维度来自变异个体
        
        for j in range(self.n_vars):
            if np.random.random() < self.CR or j == j_rand:
                trial[j] = mutant[j]
            else:
                trial[j] = target[j]
        
        return trial
    
    def select(self, target_idx: int, trial: np.ndarray):
        """
        选择操作
        
        参数:
            target_idx: 目标个体索引
            trial: 试验个体
        """
        # 评估试验个体
        trial_objective = self.objective_func(trial)
        
        if self.constraint_func:
            constraints = self.constraint_func(trial)
            trial_fitness = trial_objective + 1e6 * np.sum(np.maximum(0, constraints))
        else:
            trial_fitness = trial_objective
        
        # 贪婪选择
        if trial_fitness < self.fitness[target_idx]:
            self.population[target_idx] = trial
            self.fitness[target_idx] = trial_fitness
            self.objectives[target_idx] = trial_objective
            
            # 更新最优
            if trial_fitness < self.best_fitness:
                self.best_individual = trial.copy()
                self.best_fitness = trial_fitness
                self.best_objective = trial_objective
    
    def compute_diversity(self) -> float:
        """计算种群多样性"""
        mean = np.mean(self.population, axis=0)
        distances = np.linalg.norm(self.population - mean, axis=1)
        return np.mean(distances)
    
    def optimize(self,
                 objective_func: Callable,
                 xl: np.ndarray,
                 xu: np.ndarray,
                 constraint_func: Callable = None) -> np.ndarray:
        """
        执行优化
        
        参数:
            objective_func: 目标函数
            xl, xu: 变量下界和上界
            constraint_func: 约束函数
        """
        self.objective_func = objective_func
        self.constraint_func = constraint_func
        self.xl = xl
        self.xu = xu
        
        print("=" * 60)
        print("差分进化算法 (Differential Evolution)")
        print("=" * 60)
        print(f"设计变量数: {self.n_vars}")
        print(f"种群大小: {self.population_size}")
        print(f"迭代次数: {self.max_iter}")
        print(f"缩放因子 F: {self.F}")
        print(f"交叉概率 CR: {self.CR}")
        print(f"变异策略: {self.strategy}")
        
        # 初始化
        self.initialize_population()
        
        # 迭代优化
        for iteration in range(self.max_iter):
            for i in range(self.population_size):
                # 变异
                mutant = self.mutate(i)
                
                # 交叉
                trial = self.crossover(self.population[i], mutant)
                
                # 选择
                self.select(i, trial)
            
            # 记录历史
            self.history['best_fitness'].append(self.best_fitness)
            self.history['mean_fitness'].append(np.mean(self.fitness))
            self.history['best_objective'].append(self.best_objective)
            self.history['diversity'].append(self.compute_diversity())
            self.history['F_history'].append(self.F)
            self.history['CR_history'].append(self.CR)
            
            # 输出进度
            if iteration % 20 == 0:
                print(f"Iteration {iteration:3d}: "
                      f"Best={self.best_fitness:.6f}, "
                      f"Mean={np.mean(self.fitness):.6f}, "
                      f"Div={self.history['diversity'][-1]:.4f}")
        
        print(f"\n优化完成!")
        print(f"最优目标值: {self.best_objective:.6f}")
        print(f"最优解: {self.best_individual}")
        
        return self.best_individual
    
    def plot_convergence(self):
        """绘制收敛曲线"""
        fig, axes = plt.subplots(2, 2, figsize=(12, 10))
        
        # 最优适应度
        ax = axes[0, 0]
        ax.semilogy(self.history['best_fitness'], 'b-', linewidth=2, label='Best')
        ax.semilogy(self.history['mean_fitness'], 'g--', alpha=0.5, label='Mean')
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Fitness (log scale)')
        ax.set_title('Fitness Convergence')
        ax.legend()
        ax.grid(True, alpha=0.3)
        
        # 种群多样性
        ax = axes[0, 1]
        ax.plot(self.history['diversity'], 'r-', linewidth=2)
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Population Diversity')
        ax.set_title('Population Diversity')
        ax.grid(True, alpha=0.3)
        
        # 参数变化
        ax = axes[1, 0]
        ax.plot(self.history['F_history'], 'm-', linewidth=2, label='F')
        ax.set_xlabel('Iteration')
        ax.set_ylabel('F Value')
        ax.set_title('Scaling Factor F')
        ax.legend()
        ax.grid(True, alpha=0.3)
        
        ax = axes[1, 1]
        ax.plot(self.history['CR_history'], 'c-', linewidth=2, label='CR')
        ax.set_xlabel('Iteration')
        ax.set_ylabel('CR Value')
        ax.set_title('Crossover Rate CR')
        ax.legend()
        ax.grid(True, alpha=0.3)
        
        plt.tight_layout()
        return fig


class AdaptiveDE(DE):
    """自适应差分进化算法 (jDE)"""
    
    def __init__(self, n_vars: int, population_size: int = 50, max_iter: int = 100,
                 F_init: float = 0.5, CR_init: float = 0.5,
                 tau_F: float = 0.1, tau_CR: float = 0.1,
                 Fl: float = 0.1, Fu: float = 0.9,
                 strategy: str = 'rand/1/bin'):
        """
        初始化自适应DE
        
        参数:
            tau_F, tau_CR: 自适应参数
            Fl, Fu: F的上下界
        """
        super().__init__(n_vars, population_size, max_iter, F_init, CR_init, strategy)
        self.tau_F = tau_F
        self.tau_CR = tau_CR
        self.Fl = Fl
        self.Fu = Fu
        
        # 每个个体的F和CR
        self.F_values = None
        self.CR_values = None
    
    def initialize_population(self):
        """初始化种群和自适应参数"""
        super().initialize_population()
        self.F_values = np.ones(self.population_size) * self.F
        self.CR_values = np.ones(self.population_size) * self.CR
    
    def mutate(self, target_idx: int) -> np.ndarray:
        """自适应变异"""
        # 自适应调整F
        if np.random.random() < self.tau_F:
            self.F_values[target_idx] = self.Fl + np.random.random() * (self.Fu - self.Fl)
        
        self.F = self.F_values[target_idx]
        return super().mutate(target_idx)
    
    def crossover(self, target: np.ndarray, mutant: np.ndarray) -> np.ndarray:
        """自适应交叉"""
        # 自适应调整CR
        target_idx = np.where((self.population == target).all(axis=1))[0][0]
        if np.random.random() < self.tau_CR:
            self.CR_values[target_idx] = np.random.random()
        
        self.CR = self.CR_values[target_idx]
        return super().crossover(target, mutant)


# 测试函数
def sphere(x):
    """Sphere函数"""
    return np.sum(x ** 2)


def rastrigin(x):
    """Rastrigin函数"""
    A = 10
    n = len(x)
    return A * n + np.sum(x ** 2 - A * np.cos(2 * np.pi * x))


def ackley(x):
    """Ackley函数"""
    a = 20
    b = 0.2
    c = 2 * np.pi
    n = len(x)
    
    sum1 = np.sum(x ** 2)
    sum2 = np.sum(np.cos(c * x))
    
    term1 = -a * np.exp(-b * np.sqrt(sum1 / n))
    term2 = -np.exp(sum2 / n)
    
    return term1 + term2 + a + np.exp(1)


def rosenbrock(x):
    """Rosenbrock函数"""
    return np.sum(100 * (x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)


def griewank(x):
    """Griewank函数"""
    sum_term = np.sum(x ** 2) / 4000
    prod_term = np.prod(np.cos(x / np.sqrt(np.arange(1, len(x) + 1))))
    return sum_term - prod_term + 1


def zakharov(x):
    """Zakharov函数"""
    sum1 = np.sum(x ** 2)
    sum2 = np.sum(0.5 * np.arange(1, len(x) + 1) * x)
    return sum1 + sum2 ** 2 + sum2 ** 4


def truss_weight(x):
    """桁架重量"""
    rho = 7850
    L = np.array([1.0, 1.0, 1.414, 1.0, 1.414, 1.0, 1.414, 1.0, 1.0, 1.414])
    return rho * np.sum(x * L) / 1000


def truss_constraint(x):
    """桁架约束"""
    return 0.0001 - x


def compare_strategies():
    """比较不同变异策略"""
    print("\n" + "=" * 60)
    print("比较不同变异策略")
    print("=" * 60)
    
    strategies = ['rand/1', 'rand/2', 'best/1', 'best/2', 'current-to-best/1', 'rand-to-best/1']
    n_vars = 10
    max_iter = 100
    
    xl = -5.12 * np.ones(n_vars)
    xu = 5.12 * np.ones(n_vars)
    
    results = {}
    
    for strategy in strategies:
        print(f"\n测试策略: {strategy}")
        de = DE(n_vars=n_vars, population_size=50, max_iter=max_iter,
                F=0.8, CR=0.9, strategy=strategy)
        best = de.optimize(sphere, xl, xu)
        results[strategy] = {
            'best': de.best_objective,
            'history': de.history['best_objective']
        }
    
    # 绘制对比图
    fig, ax = plt.subplots(figsize=(10, 6))
    
    for strategy, data in results.items():
        ax.semilogy(data['history'], linewidth=2, label=f'{strategy} ({data["best"]:.2e})')
    
    ax.set_xlabel('Iteration')
    ax.set_ylabel('Best Objective (log scale)')
    ax.set_title('Strategy Comparison on Sphere Function')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('de_strategy_comparison.png', dpi=150, bbox_inches='tight')
    print("\n已保存: de_strategy_comparison.png")
    plt.close()
    
    return results


def compare_F_CR():
    """比较不同F和CR参数"""
    print("\n" + "=" * 60)
    print("比较不同F和CR参数")
    print("=" * 60)
    
    F_values = [0.5, 0.8, 1.0]
    CR_values = [0.5, 0.9, 1.0]
    
    n_vars = 10
    max_iter = 100
    xl = -5.12 * np.ones(n_vars)
    xu = 5.12 * np.ones(n_vars)
    
    fig, axes = plt.subplots(len(F_values), len(CR_values), figsize=(15, 12))
    
    for i, F in enumerate(F_values):
        for j, CR in enumerate(CR_values):
            print(f"\n测试: F={F}, CR={CR}")
            de = DE(n_vars=n_vars, population_size=50, max_iter=max_iter,
                    F=F, CR=CR, strategy='rand/1')
            best = de.optimize(rastrigin, xl, xu)
            
            ax = axes[i, j]
            ax.semilogy(de.history['best_objective'], linewidth=2)
            ax.set_title(f'F={F}, CR={CR}\nBest={de.best_objective:.2f}')
            ax.set_xlabel('Iteration')
            ax.set_ylabel('Objective (log)')
            ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('de_parameter_comparison.png', dpi=150, bbox_inches='tight')
    print("\n已保存: de_parameter_comparison.png")
    plt.close()


def main():
    """主函数"""
    print("=" * 60)
    print("差分进化算法 (Differential Evolution)")
    print("=" * 60)
    
    # 测试1: Sphere函数
    print("\n测试1: Sphere函数优化")
    print("-" * 60)
    
    de_sphere = DE(
        n_vars=10,
        population_size=50,
        max_iter=100,
        F=0.8,
        CR=0.9,
        strategy='rand/1'
    )
    
    xl = -5.12 * np.ones(10)
    xu = 5.12 * np.ones(10)
    
    best_sphere = de_sphere.optimize(sphere, xl, xu)
    
    fig = de_sphere.plot_convergence()
    plt.savefig('de_sphere_convergence.png', dpi=150, bbox_inches='tight')
    print("已保存: de_sphere_convergence.png")
    plt.close()
    
    # 测试2: Rastrigin函数
    print("\n测试2: Rastrigin函数优化")
    print("-" * 60)
    
    de_rastrigin = DE(
        n_vars=10,
        population_size=60,
        max_iter=200,
        F=0.9,
        CR=0.9,
        strategy='best/1'
    )
    
    best_rastrigin = de_rastrigin.optimize(rastrigin, xl, xu)
    
    fig = de_rastrigin.plot_convergence()
    plt.savefig('de_rastrigin_convergence.png', dpi=150, bbox_inches='tight')
    print("已保存: de_rastrigin_convergence.png")
    plt.close()
    
    # 测试3: Rosenbrock函数
    print("\n测试3: Rosenbrock函数优化")
    print("-" * 60)
    
    de_rosenbrock = DE(
        n_vars=10,
        population_size=70,
        max_iter=300,
        F=0.8,
        CR=0.95,
        strategy='current-to-best/1'
    )
    
    xl_rosen = -2.048 * np.ones(10)
    xu_rosen = 2.048 * np.ones(10)
    
    best_rosenbrock = de_rosenbrock.optimize(rosenbrock, xl_rosen, xu_rosen)
    
    fig = de_rosenbrock.plot_convergence()
    plt.savefig('de_rosenbrock_convergence.png', dpi=150, bbox_inches='tight')
    print("已保存: de_rosenbrock_convergence.png")
    plt.close()
    
    # 测试4: Ackley函数
    print("\n测试4: Ackley函数优化")
    print("-" * 60)
    
    de_ackley = DE(
        n_vars=10,
        population_size=50,
        max_iter=150,
        F=0.8,
        CR=0.9,
        strategy='rand/1'
    )
    
    xl_ackley = -32.768 * np.ones(10)
    xu_ackley = 32.768 * np.ones(10)
    
    best_ackley = de_ackley.optimize(ackley, xl_ackley, xu_ackley)
    
    fig = de_ackley.plot_convergence()
    plt.savefig('de_ackley_convergence.png', dpi=150, bbox_inches='tight')
    print("已保存: de_ackley_convergence.png")
    plt.close()
    
    # 测试5: 桁架重量优化
    print("\n测试5: 桁架重量优化")
    print("-" * 60)
    
    de_truss = DE(
        n_vars=10,
        population_size=50,
        max_iter=100,
        F=0.8,
        CR=0.9,
        strategy='best/1'
    )
    
    xl_truss = 0.0001 * np.ones(10)
    xu_truss = 0.01 * np.ones(10)
    
    best_truss = de_truss.optimize(truss_weight, xl_truss, xu_truss, truss_constraint)
    
    fig = de_truss.plot_convergence()
    plt.savefig('de_truss_convergence.png', dpi=150, bbox_inches='tight')
    print("已保存: de_truss_convergence.png")
    plt.close()
    
    # 测试6: 自适应DE
    print("\n测试6: 自适应差分进化 (jDE)")
    print("-" * 60)
    
    jde = AdaptiveDE(
        n_vars=10,
        population_size=50,
        max_iter=100,
        strategy='rand/1'
    )
    
    best_jde = jde.optimize(rastrigin, xl, xu)
    
    fig = jde.plot_convergence()
    plt.savefig('de_jde_convergence.png', dpi=150, bbox_inches='tight')
    print("已保存: de_jde_convergence.png")
    plt.close()
    
    # 策略比较
    compare_strategies()
    
    # 参数比较
    compare_F_CR()
    
    # 综合对比图
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    functions = [
        (de_sphere, 'Sphere', 'log'),
        (de_rastrigin, 'Rastrigin', 'linear'),
        (de_rosenbrock, 'Rosenbrock', 'log'),
        (de_ackley, 'Ackley', 'log'),
        (de_truss, 'Truss', 'linear')
    ]
    
    for idx, (de, name, scale) in enumerate(functions):
        row = idx // 3
        col = idx % 3
        ax = axes[row, col]
        
        if scale == 'log':
            ax.semilogy(de.history['best_objective'], 'b-', linewidth=2)
        else:
            ax.plot(de.history['best_objective'], 'b-', linewidth=2)
        
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Best Objective')
        ax.set_title(f'{name} Function')
        ax.grid(True, alpha=0.3)
    
    # jDE
    ax = axes[1, 2]
    ax.semilogy(jde.history['best_objective'], 'r-', linewidth=2)
    ax.set_xlabel('Iteration')
    ax.set_ylabel('Best Objective')
    ax.set_title('jDE on Rastrigin')
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('de_comparison.png', dpi=150, bbox_inches='tight')
    print("已保存: de_comparison.png")
    plt.close()
    
    # 打印最终结果
    print("\n" + "=" * 60)
    print("优化结果汇总")
    print("=" * 60)
    print(f"\nSphere函数:")
    print(f"  最优值: {de_sphere.best_objective:.6e}")
    
    print(f"\nRastrigin函数:")
    print(f"  最优值: {de_rastrigin.best_objective:.6f}")
    
    print(f"\nRosenbrock函数:")
    print(f"  最优值: {de_rosenbrock.best_objective:.6f}")
    
    print(f"\nAckley函数:")
    print(f"  最优值: {de_ackley.best_objective:.6e}")
    
    print(f"\n桁架重量优化:")
    print(f"  最优重量: {de_truss.best_objective:.4f} 吨")
    
    print(f"\n自适应DE (jDE):")
    print(f"  最优值: {jde.best_objective:.6f}")
    
    print("\n" + "=" * 60)
    print("差分进化优化完成!")
    print("=" * 60)


if __name__ == "__main__":
    main()