核心算法解析

1. 算法框架與初始化

class EnhancedCFO: def __init__(self, objective_func, dim=10, pop_size=30, max_iter=200, lb=-10, ub=10):

• ??改進點??：針對傳統優化算法后期易停滯的問題，結合了精英策略、多樣性控制和自適應參數
• ??關鍵特性??：
  • latin_hypercube_init()：拉丁超立方抽樣確保種群初始分布均勻
  • 動態參數：diversity_threshold,?crossover_rate,?scaling_factor?隨迭代自適應調整
  • 三重監控：stagnation_count,?diversity_history,?improvement_rate

2. 創新機制分析

a. ??強化的停滯檢測?? (update_best):

if len(self.best_history) > 5: best_range = abs(max(self.best_history[-5:]) - min(...)) tolerance = min(best_range*0.1, abs(self.best_fit)*5e-3)

• 使用滑動窗口計算動態容差閾值
• 相對改進率計算：(best_fit - current_best)/|best_fit|
• 優于傳統的固定閾值法

b. ??三階段優化策略?? (adaptive_params):

if progress < 0.3: # 探索階段 phase_shift = 1 cr = 0.9; sf = 0.8 elif progress < 0.7: # 平衡階段 phase_shift = 2 cr = 0.85; sf = 0.7 else: # 開發階段 phase_shift = 3 cr = 0.75; sf = 0.6

• 不同階段啟用不同算子組合
• 參數隨迭代進度平滑過渡

c. ??精英協作機制?? (elite_mutation):

# 精英間維度交換 swap_dims = np.random.choice(self.dim, size=int(self.dim*0.3)) self.positions[elite1, dim] = self.positions[elite2, dim]

• 前10%精英交換部分維度信息
• 突破局部最優的協同進化策略

d. ??方向性多樣化?? (diversify_population):

if np.random.rand() < 0.7: direction = position - mean_position # 遠離中心 else: direction = elite_ref - position # 趨向精英

• 70%概率遠離種群中心，30%概率趨向精英
• 平衡探索與開發的創新策略

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
import warnings
from scipy.spatial.distance import pdist
import matplotlib.cm as cm
from matplotlib.animation import FuncAnimation
import matplotlib.gridspec as gridspec
import matplotlib.animation as animation
from sklearn.decomposition import PCA
import matplotlib.colors as mcolors
from matplotlib.animation import FuncAnimation
from sklearn.decomposition import PCAwarnings.filterwarnings("ignore")class EnhancedCFO:def __init__(self, objective_func, dim=10, pop_size=30, max_iter=200, lb=-10, ub=10):"""增強版冬蟲夏草優化算法 - 針對后期停滯問題優化主要改進:1. 強化的停滯檢測機制 - 加入適應度改善率和多樣性指標2. 增強的多樣化策略 - 多樣化時引入方向性搜索3. 精英種群協作機制 - 防止局部最優4. 自適應參數調整優化 - 更平滑的參數變化"""self.obj_func = objective_funcself.dim = dimself.pop_size = pop_sizeself.max_iter = max_iterself.lb = lbself.ub = ubself.range = ub - lb# 算法參數self.p = 0.6  # 精英寄生概率self.diversity_threshold = 0.15  # 多樣性閾值初始值# 初始化種群self.positions = self.latin_hypercube_init()self.fitness = np.array([objective_func(x) for x in self.positions])# 初始化最優解min_idx = np.argmin(self.fitness)self.best_pos = self.positions[min_idx].copy()self.best_fit = self.fitness[min_idx]self.second_best = self.best_pos.copy()self.second_fit = self.best_fitself.best_history = [self.best_fit]  # 存儲歷史最優值# 記錄優化過程self.convergence_curve = [self.best_fit]self.execution_time = 0self.stagnation_count = 0self.stagnation_threshold = 25  # 初始停滯閾值self.diversity_history = [self.calculate_diversity()]self.improvement_rate = []  # 記錄改進率歷史# 自適應參數self.crossover_rate = 0.8self.scaling_factor = 0.7self.adaptive_step_size = 0.3# 優化階段標志self.phase_shift = 0self.diversify_count = 0self.last_improvement = 0self.runtime_start = time.time()# 新增可視化數據記錄self.pop_history = [self.positions.copy()]  # 記錄每代種群位置self.iteration_markers = []  # 記錄關鍵迭代點self.pca_model = None  # PCA模型用于高維可視化# 修復缺失的屬性初始化self.adaptive_step_size = 0.3self.crossover_rate = 0.8self.scaling_factor = 0.7# 修復參數歷史記錄初始化self.adaptive_history = {'crossover_rate': [self.crossover_rate],'scaling_factor': [self.scaling_factor],'step_size': [self.adaptive_step_size]}def latin_hypercube_init(self):"""改進的拉丁超立方抽樣初始化 - 增強邊界多樣性"""samples = np.zeros((self.pop_size, self.dim))for i in range(self.dim):segment_size = (self.ub - self.lb) / self.pop_sizestart = self.lb + segment_size * np.random.rand()samples[:, i] = np.linspace(start, self.ub, self.pop_size)np.random.shuffle(samples[:, i])return samplesdef calculate_diversity(self):"""計算種群多樣性指數 - 使用平均成對距離"""distances = pdist(self.positions, 'euclidean')return np.mean(distances) / (self.range * np.sqrt(self.dim))def update_best(self, t):"""更新全局最優解并檢測停滯狀態 - 增強停滯檢測"""min_idx = np.argmin(self.fitness)current_best = self.fitness[min_idx]# 基于歷史變化的動態容差if len(self.best_history) > 5:best_range = max(1e-10, abs(max(self.best_history[-5:]) - min(self.best_history[-5:])))tolerance = min(best_range * 0.1, abs(self.best_fit) * 5e-3)else:tolerance = max(1e-10, abs(self.best_fit) * (1e-3 + 1e-4 * (t / self.max_iter)))improvement = self.best_fit - current_best# 計算改進率（相對改進量）if improvement > tolerance:if self.best_fit != 0:rate = improvement / abs(self.best_fit)else:rate = improvementself.improvement_rate.append(rate)self.second_best = self.best_pos.copy()self.second_fit = self.best_fitself.best_fit = current_bestself.best_pos = self.positions[min_idx].copy()self.stagnation_count = 0self.last_improvement = tself.best_history.append(self.best_fit)return Trueelse:self.stagnation_count += 1self.improvement_rate.append(0)  # 沒有改進return Falsedef adaptive_params(self, t):"""動態調整算法參數 - 平滑過渡"""progress = t / self.max_iter# 基于改進率的自適應參數if len(self.improvement_rate) > 10:recent_improvement = np.mean(self.improvement_rate[-10:])step_size = max(0.1, min(0.5, 0.4 * np.exp(-10 * recent_improvement)))self.adaptive_step_size = step_sizeelse:step_size = self.adaptive_step_size# 基于種群多樣性的參數調整diversity = self.diversity_history[-1]self.diversity_threshold = max(0.05, 0.2 * np.exp(-5 * diversity))# 三階段策略if progress < 0.3:self.phase_shift = 1self.stagnation_threshold = max(20, 25 * (1 - progress))cr = 0.9sf = 0.8elif progress < 0.7:self.phase_shift = 2self.stagnation_threshold = max(25, 30 * (1 - progress))cr = 0.85sf = 0.7else:self.phase_shift = 3self.stagnation_threshold = max(15, 20 * (1 - progress))  # 降低后期閾值cr = 0.75sf = 0.6# 基于停滯計數的參數調整stagnation_factor = min(1.0, self.stagnation_count / self.stagnation_threshold)cr *= 1 - stagnation_factor * 0.3sf *= 1 + stagnation_factor * 0.2  # 增加變異強度self.crossover_rate = crself.scaling_factor = sf# 記錄參數歷史self.adaptive_history['crossover_rate'].append(cr)self.adaptive_history['scaling_factor'].append(sf)self.adaptive_history['step_size'].append(self.adaptive_step_size)def hybrid_exploration(self, t):"""增強混合探索策略 - 引入方向導數"""improved = Falsefor i in range(self.pop_size):if np.random.rand() < 0.7:# 基于方向的探索direction = self.best_pos - self.positions[i]direction /= np.linalg.norm(direction) + 1e-10# 自適應步長step_size = self.adaptive_step_size * (1.0 + np.random.randn() * 0.5)step_vector = step_size * direction + np.random.randn(self.dim) * 0.3new_position = np.clip(self.positions[i] + step_vector, self.lb, self.ub)else:# 差分進化 - 引入最優向量idxs = np.random.choice(self.pop_size, 3, replace=False)r1, r2, r3 = self.positions[idxs]if np.random.rand() < 0.5:  # 50%概率使用最優向量best_weight = np.random.rand() * 0.7mutation = (r1 + self.best_pos) * best_weight + self.scaling_factor * (r2 - r3)else:mutation = r1 + self.scaling_factor * (r2 - r3)cross_points = np.random.rand(self.dim) < self.crossover_rateif not np.any(cross_points):cross_points[np.random.randint(self.dim)] = Truenew_position = np.where(cross_points, mutation, self.positions[i])new_position = np.clip(new_position, self.lb, self.ub)# 評估并更新new_fitness = self.obj_func(new_position)if new_fitness < self.fitness[i]:self.positions[i] = new_positionself.fitness[i] = new_fitnessif new_fitness < self.best_fit:improved = True# 更新全局最優improved = self.update_best(t) or improvedreturn improveddef elite_mutation(self, intensity=1.0):"""精英突變算子 - 自適應突變強度"""elite_count = max(3, int(self.pop_size * 0.15))elite_indices = np.argsort(self.fitness)[:elite_count]for idx in elite_indices:# 基于停滯計數的自適應突變強度mutation_intensity = (self.ub - self.lb) * 0.08 * (intensity + self.stagnation_count / self.stagnation_threshold)mutation = mutation_intensity * np.random.randn(self.dim)new_position = np.clip(self.positions[idx] + mutation, self.lb, self.ub)new_fitness = self.obj_func(new_position)if new_fitness < self.fitness[idx]:self.positions[idx] = new_positionself.fitness[idx] = new_fitnessif new_fitness < self.best_fit:self.update_best(self.last_improvement)# 精英協作機制 - 前10%精英交換信息for i in range(1, len(elite_indices)):# 精英個體間交換部分維度信息swap_dims = np.random.choice(self.dim, size=max(1, int(self.dim * 0.3)), replace=False)for dim_idx in swap_dims:temp = self.positions[elite_indices[i - 1], dim_idx]self.positions[elite_indices[i - 1], dim_idx] = self.positions[elite_indices[i], dim_idx]self.positions[elite_indices[i], dim_idx] = temp# 重新評估適應度self.fitness[elite_indices[i - 1]] = self.obj_func(self.positions[elite_indices[i - 1]])self.fitness[elite_indices[i]] = self.obj_func(self.positions[elite_indices[i]])def spiral_rising(self, t):"""增強螺旋上升算子 - 非線性適應性螺旋"""improved = False# 基于進化階段的自適應螺旋參數a = 0.4 - 0.3 * (t / self.max_iter)b = 0.5 + 1.0 * (t / self.max_iter)# 三引導源（動態權重）weights = np.random.dirichlet(np.ones(3), 1)[0]target = weights[0] * self.best_pos + weights[1] * self.second_best + weights[2] * self.positions[np.random.randint(self.pop_size)]for i in range(self.pop_size):# 自適應螺旋參數r_val = a * np.exp(b * 2 * np.pi * np.random.rand())# 方向向量direction = target - self.positions[i]direction /= np.linalg.norm(direction) + 1e-10# 螺旋運動加上隨機擾動random_component = 0.1 * np.random.randn(self.dim)new_position = self.positions[i] + r_val * direction + random_componentnew_position = np.clip(new_position, self.lb, self.ub)# 評估并更新new_fitness = self.obj_func(new_position)if new_fitness < self.fitness[i]:self.positions[i] = new_positionself.fitness[i] = new_fitnessif new_fitness < self.best_fit:improved = True# 更新全局最優improved = self.update_best(t) or improvedreturn improveddef cooperative_parasitism(self, t):"""協同寄生機制 - 多宿主策略"""improved = Falsefitness_order = np.argsort(self.fitness)elite_count = max(5, int(self.pop_size * 0.3))elite_pool = fitness_order[:elite_count]ordinary_pool = fitness_order[elite_count:]for i in range(self.pop_size):# 精英宿主概率自適應調整elite_prob = max(0.4, min(0.8, self.p * (1 - t / self.max_iter)))# 多宿主策略（隨機選擇1-3個宿主）host_count = np.random.randint(1, 4)displacement = np.zeros(self.dim)for _ in range(host_count):if np.random.rand() < elite_prob and elite_pool.size > 0:host_idx = elite_pool[np.random.randint(len(elite_pool))]else:host_idx = ordinary_pool[np.random.randint(len(ordinary_pool))] if ordinary_pool.size > 0 else i# 自適應步長控制step_size = self.adaptive_step_size * (1.0 + np.random.randn() * 0.3)displacement += step_size * (self.positions[host_idx] - self.positions[i])displacement /= host_count  # 平均位移new_position = np.clip(self.positions[i] + displacement, self.lb, self.ub)# 評估并更新new_fitness = self.obj_func(new_position)if new_fitness < self.fitness[i]:self.positions[i] = new_positionself.fitness[i] = new_fitnessif new_fitness < self.best_fit:improved = True# 更新全局最優improved = self.update_best(t) or improvedreturn improveddef diversify_population(self, t):"""增強多樣化策略 - 方向性多樣化"""self.diversify_count += 1elite_count = max(5, int(self.pop_size * 0.25))elite_indices = np.argsort(self.fitness)[:elite_count]elite_positions = self.positions[elite_indices].copy()# 方向性多樣化mean_position = np.mean(self.positions, axis=0)for i in range(self.pop_size):if i not in elite_indices:# 選擇一個精英作為方向參考elite_ref = elite_positions[np.random.randint(elite_count)]# 多樣化方向：遠離種群中心或向精英方向探索if np.random.rand() < 0.7:# 遠離種群中心direction = self.positions[i] - mean_positionelse:# 向精英方向探索direction = elite_ref - self.positions[i]direction /= np.linalg.norm(direction) + 1e-10# 自適應多樣化步長mutate_factor = 0.4 + 0.3 * self.stagnation_count / self.stagnation_thresholdmutation = mutate_factor * (self.ub - self.lb) * direction * np.random.randn()self.positions[i] = np.clip(self.positions[i] + mutation, self.lb, self.ub)self.fitness[i] = self.obj_func(self.positions[i])# 更新精英位置 - 小幅擾動for i, idx in enumerate(elite_indices):mutation = 0.1 * (self.ub - self.lb) * np.random.randn(self.dim)elite_positions[i] = np.clip(elite_positions[i] + mutation, self.lb, self.ub)self.positions[idx] = elite_positions[i]self.fitness[idx] = self.obj_func(elite_positions[i])# 更新全局最優min_idx = np.argmin(self.fitness)if self.fitness[min_idx] < self.best_fit:self.best_pos = self.positions[min_idx].copy()self.best_fit = self.fitness[min_idx]self.best_history.append(self.best_fit)self.stagnation_count = 0diversity = self.calculate_diversity()self.diversity_history.append(diversity)self.iteration_markers.append(("diversify", t, self.best_fit))def optimize(self):"""主優化流程 - 增強停滯檢測與處理"""print(f"Starting Enhanced CFO v2 | DIM: {self.dim} | POP: {self.pop_size}")print("Iter    Best Fitness    Diversity    Stag    ImpRate")print("----------------------------------------------------")self.runtime_start = time.time()for t in range(self.max_iter):self.adaptive_params(t)# 基于改進率的停滯檢測stagnation_condition = Falseif len(self.improvement_rate) > 5:# 計算最近n次迭代的平均改進率recent_improvement = np.mean(self.improvement_rate[-5:])# 檢測改進停滯：連續5次平均改進率小于閾值stagnation_condition = (recent_improvement < 1e-5 andself.stagnation_count > self.stagnation_threshold)# 執行多樣化if stagnation_condition:self.diversify_population(t)# 多樣化后進行一次精英突變增強self.elite_mutation(intensity=1.5)# 重置無重大改進的計數器current_improvement = False# 階段優化策略if self.phase_shift == 1:current_improvement = self.hybrid_exploration(t) or current_improvementcurrent_improvement = self.cooperative_parasitism(t) or current_improvementelif self.phase_shift == 2:current_improvement = self.hybrid_exploration(t) or current_improvementcurrent_improvement = self.cooperative_parasitism(t) or current_improvementcurrent_improvement = self.spiral_rising(t) or current_improvementelse:current_improvement = self.spiral_rising(t) or current_improvementcurrent_improvement = self.cooperative_parasitism(t) or current_improvementself.elite_mutation()  # 后期增強精英突變# 記錄當前種群狀態self.pop_history.append(self.positions.copy())diversity = self.calculate_diversity()self.diversity_history.append(diversity)# 記錄收斂曲線self.convergence_curve.append(self.best_fit)# 進度報告if t % 20 == 0 or t == self.max_iter - 1 or current_improvement:imp_rate = self.improvement_rate[-1] if self.improvement_rate else 0print(f"{t:4d}    {self.best_fit:.6e}    {diversity:.4f}    {self.stagnation_count:2d}    {imp_rate:.3e}")# 完成優化self.execution_time = time.time() - self.runtime_startprint("\nOptimization Completed!")print(f"Execution time: {self.execution_time:.2f} seconds")print(f"Best Fitness: {self.best_fit:.8e}")print(f"Diversifications: {self.diversify_count}")# 訓練PCA模型用于高維可視化if self.dim > 3:all_positions = np.vstack(self.pop_history)self.pca_model = PCA(n_components=2)self.pca_model.fit(all_positions)return self.best_pos, self.best_fit, self.convergence_curvedef enhanced_visualization(self, save_path=None):"""增強版可視化：多維分析+動畫+熱力圖"""fig = plt.figure(figsize=(18, 12), constrained_layout=True)fig.suptitle(f'Enhanced CFO Optimization Visualization | Dim={self.dim}, Pop={self.pop_size}, Iter={self.max_iter}',fontsize=16, fontweight='bold')gs = gridspec.GridSpec(2, 2, figure=fig, width_ratios=[1, 1], height_ratios=[1.5, 1])# 1. 搜索空間可視化面板ax1 = fig.add_subplot(gs[0, 0], projection='3d' if self.dim >= 3 else None)# 為低維問題創建搜索網格if self.dim == 2:x = np.linspace(self.lb, self.ub, 100)y = np.linspace(self.lb, self.ub, 100)X, Y = np.meshgrid(x, y)Z = np.zeros_like(X)for i in range(100):for j in range(100):Z[i, j] = self.obj_func(np.array([X[i, j], Y[i, j]]))# 繪制等高線contour = ax1.contourf(X, Y, Z, levels=50, cmap='viridis', alpha=0.7)fig.colorbar(contour, ax=ax1, label='Fitness Value')# 繪制初始和最終種群ax1.scatter(self.pop_history[0][:, 0], self.pop_history[0][:, 1],c='blue', s=40, marker='o', alpha=0.6, label='Initial Population')ax1.scatter(self.pop_history[-1][:, 0], self.pop_history[-1][:, 1],c='green', s=40, marker='s', alpha=0.6, label='Final Population')# 繪制最佳解軌跡best_traj = [pop[np.argmin([self.obj_func(x) for x in pop])] for pop in self.pop_history]best_traj = np.array(best_traj)ax1.plot(best_traj[:, 0], best_traj[:, 1], 'r-', linewidth=2, label='Best Solution Path')# 標記多樣化位置for marker in self.iteration_markers:iter_idx = marker[1]if len(self.pop_history) > iter_idx:diver_pop = self.pop_history[iter_idx]ax1.scatter(diver_pop[:, 0], diver_pop[:, 1], c='red', s=60,marker='*', alpha=0.8,label='Diversification' if marker == self.iteration_markers[0] else '')ax1.set_title('Search Space & Population Trajectory (2D)')ax1.set_xlabel('Dimension 1')ax1.set_ylabel('Dimension 2')ax1.legend()elif self.dim >= 3:if self.dim > 3:# 繪制PCA投影pca_history = [self.pca_model.transform(pop) for pop in self.pop_history]# 繪制所有種群的PCA投影colors = plt.cm.jet(np.linspace(0, 1, len(pca_history)))for i, proj in enumerate(pca_history):ax1.scatter(proj[:, 0], proj[:, 1], color=colors[i], s=20,alpha=0.2 + 0.8 * i / len(pca_history))# 繪制最佳解的軌跡best_proj_history = []for i, pop in enumerate(self.pop_history):best_idx = np.argmin([self.obj_func(x) for x in pop])best_proj_history.append(pca_history[i][best_idx])best_proj = np.array(best_proj_history)ax1.plot(best_proj[:, 0], best_proj[:, 1], 'r-', linewidth=2, alpha=0.7)ax1.scatter(best_proj[-1, 0], best_proj[-1, 1], s=100, c='gold',marker='*', edgecolor='k', label='Best Solution')# 標記多樣化位置for marker in self.iteration_markers:iter_idx = marker[1]if len(pca_history) > iter_idx:ax1.scatter(pca_history[iter_idx][:, 0], pca_history[iter_idx][:, 1],s=80, facecolors='none', edgecolors='r',linewidth=1.5,label='Diversification' if marker == self.iteration_markers[0] else '')ax1.set_xlabel('Principal Component 1')ax1.set_ylabel('Principal Component 2')ax1.set_title('Population Dynamics Projection (PCA)')# 添加顏色條表示迭代進度sm = plt.cm.ScalarMappable(cmap=plt.cm.jet, norm=plt.Normalize(vmin=0, vmax=self.max_iter))sm.set_array([])cbar = fig.colorbar(sm, ax=ax1, label='Iteration')else:# 3D空間的軌跡可視化ax1 = fig.add_subplot(gs[0, 0], projection='3d')# 繪制初始和最終種群ax1.scatter(self.pop_history[0][:, 0], self.pop_history[0][:, 1],self.pop_history[0][:, 2], c='blue', s=30, label='Initial Population')ax1.scatter(self.pop_history[-1][:, 0], self.pop_history[-1][:, 1],self.pop_history[-1][:, 2], c='green', s=50, marker='s', label='Final Population')# 繪制最佳解軌跡best_traj = [pop[np.argmin([self.obj_func(x) for x in pop])] for pop in self.pop_history]best_traj = np.array(best_traj)ax1.plot(best_traj[:, 0], best_traj[:, 1], best_traj[:, 2],'r-', linewidth=2, label='Best Solution Path')# 標記多樣化位置for marker in self.iteration_markers:iter_idx = marker[1]if len(self.pop_history) > iter_idx:diver_pop = self.pop_history[iter_idx]ax1.scatter(diver_pop[:, 0], diver_pop[:, 1], diver_pop[:, 2],c='red', s=60, marker='*', alpha=0.8,label='Diversification' if marker == self.iteration_markers[0] else '')ax1.set_title('Population Trajectory in 3D Search Space')ax1.set_xlabel('Dimension 1')ax1.set_ylabel('Dimension 2')ax1.set_zlabel('Dimension 3')ax1.legend()# 2. 收斂分析面板ax2 = fig.add_subplot(gs[0, 1])# 繪制收斂曲線convergence_line = ax2.semilogy(self.convergence_curve, 'b-', linewidth=1.5, label='Best Fitness')[0]ax2.set_xlabel('Iteration')ax2.set_ylabel('Fitness Value (log scale)', color='b')ax2.tick_params(axis='y', labelcolor='b')ax2.grid(True, linestyle='--', alpha=0.7)# 添加多樣性曲線ax3 = ax2.twinx()diversity_line = ax3.plot(self.diversity_history, 'g-', linewidth=1.5, label='Diversity')[0]ax3.set_ylabel('Diversity Index', color='g')ax3.tick_params(axis='y', labelcolor='g')# 標記多樣化操作for marker in self.iteration_markers:t, fitness = marker[1], marker[2]ax2.axvline(x=t, color='r', linestyle='--', alpha=0.6)ax2.text(t, self.convergence_curve[t] * 10, 'Diversification',rotation=90, fontsize=9, ha='right', va='bottom')# 添加圖例lines = [convergence_line, diversity_line]ax2.legend(lines, [l.get_label() for l in lines], loc='upper right')ax2.set_title('Convergence & Diversity Analysis')# 3. 參數動態變化面板ax4 = fig.add_subplot(gs[1, 0])# 修復：確保所有參數歷史長度一致min_length = min(len(self.adaptive_history['crossover_rate']),len(self.adaptive_history['scaling_factor']),len(self.adaptive_history['step_size']))# 提取參數歷史iterations = np.arange(min_length)# 繪制參數變化曲線ax4.plot(iterations, self.adaptive_history['crossover_rate'][:min_length],'o-', markersize=4, linewidth=1.5, color='blue', label='Crossover Rate')ax4.plot(iterations, self.adaptive_history['scaling_factor'][:min_length],'s-', markersize=4, linewidth=1.5, color='green', label='Scaling Factor')ax4.plot(iterations, self.adaptive_history['step_size'][:min_length],'^-', markersize=4, linewidth=1.5, color='purple', label='Step Size')ax4.set_xlabel('Iteration')ax4.set_ylabel('Parameter Value')ax4.legend()ax4.grid(True, linestyle='--', alpha=0.7)ax4.set_title('Algorithm Parameter Adaptation')# 標記參數變化點for marker in self.iteration_markers:if marker[1] < min_length:ax4.axvline(x=marker[1], color='r', linestyle='--', alpha=0.3)# 4. 種群密度與位置熱圖面板ax5 = fig.add_subplot(gs[1, 1])if self.dim == 2:# 創建密度網格x_bins = np.linspace(self.lb, self.ub, 50)y_bins = np.linspace(self.lb, self.ub, 50)density = np.zeros((len(y_bins) - 1, len(x_bins) - 1))# 統計所有迭代中的個體位置for pop in self.pop_history:for ind in pop:x_idx = np.digitize(ind[0], x_bins) - 1y_idx = np.digitize(ind[1], y_bins) - 1if 0 <= x_idx < density.shape[1] and 0 <= y_idx < density.shape[0]:density[y_idx, x_idx] += 1# 應用對數縮放density = np.log1p(density)# 繪制熱力圖im = ax5.imshow(density, cmap='hot', origin='lower',extent=[self.lb, self.ub, self.lb, self.ub],aspect='auto', alpha=0.8)# 添加等高線x = np.linspace(self.lb, self.ub, 100)y = np.linspace(self.lb, self.ub, 100)X, Y = np.meshgrid(x, y)Z = np.zeros_like(X)for i in range(100):for j in range(100):Z[i, j] = self.obj_func(np.array([X[i, j], Y[i, j]]))contours = ax5.contour(X, Y, Z, levels=15, colors='cyan', alpha=0.5)ax5.clabel(contours, inline=True, fontsize=8, fmt='%1.1f')# 繪制最佳解位置ax5.plot(self.best_pos[0], self.best_pos[1], '*', markersize=15,markeredgecolor='black', markerfacecolor='gold',label=f'Best Solution: {self.best_fit:.2e}')# 添加顏色條fig.colorbar(im, ax=ax5, label='Log Population Density')ax5.set_title('Population Distribution Heatmap & Fitness Landscape')ax5.legend()else:ax5.text(0.5, 0.5, 'Population heatmap only available for 2D problems',ha='center', va='center', fontsize=12)ax5.set_title('Population Heatmap Not Available for High Dimensions')# 保存結果if save_path:plt.savefig(save_path, dpi=150, bbox_inches='tight')print(f"Enhanced visualization saved to {save_path}")plt.show()# 如果維度為2，創建進化過程動畫if self.dim == 2 and len(self.pop_history) > 1:print("\nGenerating optimization animation...")self.create_animation()return figdef create_animation(self, filename='cfo_animation.gif'):"""創建優化過程動畫 (僅適用于2D問題)"""if self.dim != 2:print("Animation only available for 2D problems")return# 創建網格用于可視化x = np.linspace(self.lb, self.ub, 100)y = np.linspace(self.lb, self.ub, 100)X, Y = np.meshgrid(x, y)Z = np.zeros_like(X)for i in range(100):for j in range(100):Z[i, j] = self.obj_func(np.array([X[i, j], Y[i, j]]))# 減少幀數以加速動畫step = max(1, len(self.pop_history) // 100)frames = range(0, len(self.pop_history), step)fig, ax = plt.subplots(figsize=(10, 8))contour = ax.contourf(X, Y, Z, levels=50, cmap='viridis', alpha=0.7)fig.colorbar(contour, ax=ax, label='Fitness Value')ax.set_xlim(self.lb, self.ub)ax.set_ylim(self.lb, self.ub)ax.set_xlabel('Dimension 1')ax.set_ylabel('Dimension 2')ax.set_title('Enhanced CFO Optimization Process')# 初始散點圖scat = ax.scatter([], [], s=50, c='white', ec='black', alpha=0.8, zorder=3)best_point = ax.scatter([], [], s=100, marker='*', c='gold', ec='black', zorder=4)best_traj, = ax.plot([], [], 'r-', linewidth=1.5, zorder=2)best_positions = []iteration_text = ax.text(0.05, 0.95, '', transform=ax.transAxes,fontsize=12, bbox=dict(facecolor='white', alpha=0.8))fitness_text = ax.text(0.05, 0.90, '', transform=ax.transAxes,fontsize=10, bbox=dict(facecolor='white', alpha=0.8))# 動畫更新函數def update(frame):ax.collections = [ax.collections[0]]  # 只保留背景pop = self.pop_history[frame]scat.set_offsets(pop)# 更新顏色表示適應度colors = np.array([self.obj_func(x) for x in pop])norm = mcolors.Normalize(vmin=colors.min(), vmax=colors.max())cmap = cm.plasmascat.set_color(cmap(norm(colors)))# 更新最佳解best_idx = np.argmin(colors)current_best = pop[best_idx]best_point.set_offsets([current_best])# 記錄最佳解軌跡best_positions.append(current_best)if len(best_positions) > 1:best_traj.set_data(*zip(*best_positions))# 更新文本iteration_text.set_text(f'Iteration: {frame}/{len(self.pop_history) - 1}')fitness_text.set_text(f'Best Fitness: {colors.min():.4e}\nPopulation Diversity: {self.diversity_history[frame]:.4f}')return scat, best_point, best_traj, iteration_text, fitness_text# 創建動畫ani = FuncAnimation(fig, update, frames=frames, interval=100, blit=True)# 保存動畫print(f"Saving animation to {filename}...")ani.save(filename, writer='pillow', fps=10)print("Animation saved successfully!")plt.close(fig)return ani# 測試函數定義
def test_function(func_name, dim=10):"""測試函數工廠"""def sphere(x):return np.sum(x ** 2)def rastrigin(x):A = 10return A * dim + np.sum(x ** 2 - A * np.cos(2 * np.pi * x))def rosenbrock(x):return np.sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)def ackley(x):a = 20b = 0.2c = 2 * np.pin = 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 griewank(x):sum1 = np.sum(x ** 2) / 4000prod1 = np.prod(np.cos(x / np.sqrt(np.arange(1, len(x) + 1))))return sum1 - prod1 + 1def schwefel(x):return 418.9829 * dim - np.sum(x * np.sin(np.sqrt(np.abs(x))))func_map = {'sphere': sphere,'rastrigin': rastrigin,'rosenbrock': rosenbrock,'ackley': ackley,'griewank': griewank,'schwefel': schwefel}return func_map.get(func_name.lower(), sphere)def benchmark(test_funcs=None, dim=10, runs=5):"""基準測試套件"""if test_funcs is None:test_funcs = ['sphere', 'rastrigin', 'rosenbrock', 'ackley', 'griewank']results = {}for func_name in test_funcs:print(f"\n=== Benchmarking {func_name} function (Dim={dim}) ===")best_fitness = []exec_times = []for i in range(runs):print(f"\nRun {i + 1}/{runs}")optimizer = EnhancedCFO(objective_func=test_function(func_name, dim),dim=dim,pop_size=50,max_iter=300)_, fitness, _ = optimizer.optimize()best_fitness.append(fitness)exec_times.append(optimizer.execution_time)# 每次運行后調用增強可視化if dim in [2, 3]:optimizer.enhanced_visualization(save_path=f"{func_name}_result_{i + 1}.png")avg_fitness = np.mean(best_fitness)avg_time = np.mean(exec_times)std_fitness = np.std(best_fitness)results[func_name] = {'best_fitness': best_fitness,'average_fitness': avg_fitness,'std_fitness': std_fitness,'execution_times': exec_times,'average_time': avg_time}print(f"\n{func_name} Results:")print(f"  Average Fitness: {avg_fitness:.6e}")print(f"  Fitness Std Dev: {std_fitness:.6e}")print(f"  Average Execution Time: {avg_time:.2f}s")# 打印總結報告print("\n=== Benchmark Summary ===")for func, res in results.items():print(f"{func:15} | Best: {min(res['best_fitness']):.2e} | Avg: {res['average_fitness']:.2e} | Std: {res['std_fitness']:.2e}")return resultsif __name__ == "__main__":# 單個函數優化示例（添加可視化）optimizer = EnhancedCFO(objective_func=test_function("rastrigin"),dim=2,  # 改為2維以啟用完整可視化pop_size=40,max_iter=100,lb=-5.12,ub=5.12)best_solution, best_fitness, convergence = optimizer.optimize()optimizer.enhanced_visualization(save_path="enhanced_cfo_visualization.png")print("\n=== Solution Validation ===")print(f"Best Fitness: {best_fitness:.10f}")print("Best Solution Vector:")print(np.around(best_solution, 4))# 使用gif動畫需要啟用以下行# optimizer.create_animation()# 完整基準測試（可選）# benchmark(dim=2, runs=1)  # 二維測試以便可視化