ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

帝国竞争算法ICA Python实现:从原理到可视化调参

帝国竞争算法ICA Python实现:从原理到可视化调参 简介这是一份面向优化算法学习者与Python开发者的帝国竞争算法ICA实现资源将社会政治进化中的殖民扩张与帝国吞并过程抽象为仿生优化机制用于求解多维函数优化问题。资源包共6个文件包含2个Python源码文件、2张结果图片、1份依赖清单与1份说明文档压缩包约159KB体积轻巧便于快速上手。其中源码完整实现帝国初始化、同化、竞争与吞并等核心流程图片分别呈现算法收敛曲线与帝国殖民地分布直观反映迭代过程中的种群演化。运行演示脚本即可复现实验并自动保存收敛曲线与帝国分布图帮助读者理解弱帝国被强帝国吞并、最终仅剩单一帝国的收敛逻辑。目前已有52人学习适合作为智能优化算法的入门实践与可视化参考。1. 帝国竞争算法 ICA 到底能干什么从殖民竞争到函数寻优如果你写过遗传算法或粒子群第一次看帝国竞争算法Imperialist Competitive AlgorithmICA的伪代码大概率会愣一下它把种群拆成若干「帝国」每个帝国由一个殖民国家强势解和一批殖民地弱势解组成殖民地会向宗主国靠拢帝国之间还会互相吞并弱国被强国瓜分最后只剩一个帝国。这套机制听起来像社会演化模拟本质却是一个连续优化算法2007 年由 Atashpaz-Gargari 和 Lucas 提出专门用来啃非线性、多峰、带约束的目标函数。我这次拆的是一份 Python 实现核心卖点有两个一是把 ICA 的完整迭代流程写清楚了二是带可视化能实时看到帝国版图怎么收缩、殖民地怎么迁移。对做优化调度、参数标定、神经网络超参搜索的人来说它比遗传算法多了一层「帝国吞并」的收敛压力前期探索猛后期收敛快。适合谁已经会 Python 基础语法、想找一个能跑通、能改参数、能看过程的启发式算法练手的人。如果你连 numpy 都没装先补环境后面第 2 章会讲。2. ICA 的数学骨架与 Python 落地目标函数、帝国初始化、同化与竞争2.1 为什么选 ICA 而不是 GA/PSO收敛压力与多样性平衡遗传算法靠交叉变异维持多样性粒子群靠个体最优和全局最优牵引两者都容易在后期陷入局部最优。ICA 的独特之处在于「帝国竞争」这一步每次迭代把所有帝国按势力值排序最弱的帝国会被最强帝国吞掉一个殖民地帝国数量逐渐减少。这个机制相当于给算法加了一个外部收敛压力前期帝国多、探索广后期帝国少、收敛快。从参数角度看ICA 需要调的东西比 GA 少帝国数量、殖民地数量、同化系数、革命概率、竞争系数。GA 要调交叉率、变异率、种群规模、选择策略PSO 要调惯性权重、学习因子。ICA 的参数物理意义更直观调起来不容易玄学。我一般会在 5 到 10 个帝国之间试殖民地总数控制在 50 到 200同化系数取 1.5 到 2.0革命概率 0.05 到 0.1。提示帝国数量不是越多越好。帝国太多每个帝国殖民地太少同化步长不够收敛慢帝国太少前期探索不足容易早熟。2.2 目标函数与约束处理把问题写成 ICA 能吃的形式ICA 默认处理无约束连续优化但实际工程问题大多带约束。常见做法是罚函数法把约束违反量乘以一个大系数加到目标函数上。下面是一个带边界约束的测试函数Rastrigin 函数多峰、坑多适合验证 ICA 的全局搜索能力。import numpy as np def rastrigin(x): Rastrigin 测试函数全局最优在 x0 处值为 0。 维度由输入 x 的长度决定这里默认 2 维。 A 10 n len(x) return A * n np.sum(x**2 - A * np.cos(2 * np.pi * x)) def penalty_objective(x, bounds, penalty_coef1e6): 带边界罚函数的目标函数。 x: 决策变量向量 bounds: [(low, high), ...] 每个维度的上下界 penalty_coef: 罚函数系数越大越不允许越界 # 越界惩罚超出边界的部分平方后累加 penalty 0.0 for i, (low, high) in enumerate(bounds): if x[i] low: penalty (low - x[i]) ** 2 elif x[i] high: penalty (x[i] - high) ** 2 return rastrigin(x) penalty_coef * penalty这段代码里rastrigin是标准测试函数penalty_objective把越界量平方后乘系数加到原函数上。罚函数系数不能太小否则算法会故意越界换取更低的目标值也不能太大否则数值尺度失衡同化步长失效。我一般取目标函数量级的 1e3 到 1e6 倍先跑一次看越界情况再调。2.3 帝国初始化与同化代码逐段拆解帝国初始化分两步先随机生成一批国家按目标函数值排序取前 N 个作为宗主国剩下的按势力值比例分配给各帝国。势力值通常用目标函数值的倒数或归一化后的相对值。下面是一个完整的初始化函数。def initialize_empires(n_countries, n_imperials, bounds, objective_func): 初始化帝国结构。 n_countries: 国家总数 n_imperials: 帝国数量宗主国数量 bounds: 每个维度的上下界 objective_func: 目标函数输入向量返回标量 返回: empires 列表每个元素是 dict含 imperialist 和 colonies dim len(bounds) # 随机生成国家每个维度在边界内均匀采样 countries np.random.uniform( low[b[0] for b in bounds], high[b[1] for b in bounds], size(n_countries, dim) ) # 计算每个国家的目标函数值 costs np.array([objective_func(c) for c in countries]) # 按成本升序排序成本越小越强 sorted_idx np.argsort(costs) countries countries[sorted_idx] costs costs[sorted_idx] # 前 n_imperials 个作为宗主国 imperialists countries[:n_imperials] imperialist_costs costs[:n_imperials] # 剩下的作为待分配殖民地 remaining countries[n_imperials:] remaining_costs costs[n_imperials:] # 计算每个帝国的势力值成本越小势力越大 # 用最大成本减去当前成本避免倒数带来的数值问题 max_cost np.max(imperialist_costs) powers max_cost - imperialist_costs 1e-12 powers powers / np.sum(powers) # 按势力值比例分配殖民地 empires [] start 0 for i in range(n_imperials): n_colonies int(round(powers[i] * len(remaining))) # 最后一个帝国拿走剩余所有殖民地避免取整误差 if i n_imperials - 1: n_colonies len(remaining) - start colonies remaining[start:start n_colonies] colony_costs remaining_costs[start:start n_colonies] empires.append({ imperialist: imperialists[i].copy(), imperialist_cost: imperialist_costs[i], colonies: colonies.copy(), colony_costs: colony_costs.copy() }) start n_colonies return empires关键参数说明n_countries是国家总数一般取 100 到 500n_imperials是帝国数量取 5 到 20bounds是每个维度的上下界必须和决策变量维度一致。势力值计算用max_cost - cost而不是1/cost是为了避免成本接近零时数值爆炸。分配殖民地时用round取整最后一个帝国兜底防止殖民地总数对不上。同化步骤是 ICA 的核心每个殖民地沿着指向宗主国的方向移动一个随机步长步长由同化系数控制。常见做法是加一个随机扰动角让殖民地不是直线靠近而是带一定偏转角。def assimilate(empire, assimilation_coef1.8, deviation_anglenp.pi/4): 同化操作殖民地向宗主国移动。 assimilation_coef: 同化系数控制移动步长 deviation_angle: 最大偏转角增加搜索多样性 imperialist empire[imperialist] new_colonies [] for colony in empire[colonies]: # 指向宗主国的方向向量 direction imperialist - colony # 随机偏转角在 [-deviation_angle, deviation_angle] 之间 theta np.random.uniform(-deviation_angle, deviation_angle) # 二维旋转矩阵高维时只对前两维旋转其余维度保持 if len(colony) 2: rot np.array([ [np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)] ]) direction[:2] rot direction[:2] # 移动步长同化系数乘以方向向量再加一点随机扰动 step assimilation_coef * np.random.rand() * direction new_colony colony step new_colonies.append(new_colony) empire[colonies] np.array(new_colonies) return empireassimilation_coef一般取 1.5 到 2.0太小收敛慢太大容易跳过最优解。deviation_angle取 π/4 左右给殖民地一个偏转避免所有殖民地走同一条直线。高维问题时只旋转前两维是一种简化严格做法是对每一维都加独立扰动但计算量会上去。2.4 帝国竞争与吞并最弱帝国怎么被瓜分每次迭代末尾计算每个帝国的总势力值宗主国势力加上殖民地平均势力乘以一个系数。最弱的帝国会被最强帝国吞掉一个殖民地如果最弱帝国只剩宗主国没有殖民地它就被灭国宗主国变成最强帝国的殖民地。def imperialistic_competition(empires, zeta0.02): 帝国竞争最弱帝国被最强帝国吞并一个殖民地。 zeta: 殖民地平均势力在总势力中的权重一般取 0.02 到 0.1 # 计算每个帝国的总势力 total_powers [] for emp in empires: imp_power 1.0 / (emp[imperialist_cost] 1e-12) if len(emp[colonies]) 0: col_power np.mean(1.0 / (emp[colony_costs] 1e-12)) else: col_power 0.0 total_powers.append(imp_power zeta * col_power) total_powers np.array(total_powers) # 最弱和最强帝国索引 weakest_idx np.argmin(total_powers) strongest_idx np.argmax(total_powers) weakest empires[weakest_idx] strongest empires[strongest_idx] if len(weakest[colonies]) 0: # 把最弱帝国的一个殖民地转给最强帝国 colony weakest[colonies][-1] colony_cost weakest[colony_costs][-1] weakest[colonies] weakest[colonies][:-1] weakest[colony_costs] weakest[colony_costs][:-1] strongest[colonies] np.vstack([strongest[colonies], colony]) strongest[colony_costs] np.append(strongest[colony_costs], colony_cost) else: # 最弱帝国没有殖民地灭国宗主国变成最强帝国的殖民地 empires.pop(weakest_idx) strongest[colonies] np.vstack([strongest[colonies], weakest[imperialist]]) strongest[colony_costs] np.append(strongest[colony_costs], weakest[imperialist_cost]) return empireszeta控制殖民地平均势力对帝国总势力的贡献取 0.02 到 0.1 之间。太小帝国只靠宗主国撑门面殖民地质量被忽略太大殖民地多的帝国占优可能拖慢收敛。竞争步骤每轮迭代执行一次帝国数量会逐渐减少直到只剩一个帝国。3. 可视化怎么做收敛曲线、帝国版图与殖民地迁移轨迹3.1 收敛曲线看目标函数值怎么掉下来收敛曲线是最基本的可视化横轴迭代次数纵轴当前最优目标函数值。用 matplotlib 画每轮迭代记录一次全局最优值。import matplotlib.pyplot as plt def run_ica_with_history(n_countries, n_imperials, bounds, objective_func, max_iter200, assimilation_coef1.8, deviation_anglenp.pi/4, zeta0.02): 运行 ICA 并记录每轮最优值返回最优解和历史。 empires initialize_empires(n_countries, n_imperials, bounds, objective_func) history [] best_solution None best_cost np.inf for it in range(max_iter): # 同化每个帝国 for emp in empires: emp assimilate(emp, assimilation_coef, deviation_angle) # 重新计算殖民地成本 emp[colony_costs] np.array([objective_func(c) for c in emp[colonies]]) # 殖民地如果比宗主国强交换位置 if len(emp[colonies]) 0: best_col_idx np.argmin(emp[colony_costs]) if emp[colony_costs][best_col_idx] emp[imperialist_cost]: emp[imperialist], emp[colonies][best_col_idx] \ emp[colonies][best_col_idx].copy(), emp[imperialist].copy() emp[imperialist_cost], emp[colony_costs][best_col_idx] \ emp[colony_costs][best_col_idx], emp[imperialist_cost] # 帝国竞争 empires imperialistic_competition(empires, zeta) # 记录全局最优 current_best min(emp[imperialist_cost] for emp in empires) history.append(current_best) if current_best best_cost: best_cost current_best for emp in empires: if emp[imperialist_cost] current_best: best_solution emp[imperialist].copy() break return best_solution, best_cost, history # 运行并画图 bounds [(-5.12, 5.12), (-5.12, 5.12)] best_sol, best_cost, history run_ica_with_history( n_countries200, n_imperials10, boundsbounds, objective_funclambda x: penalty_objective(x, bounds), max_iter200 ) plt.figure(figsize(8, 5)) plt.plot(history, linewidth2) plt.xlabel(Iteration) plt.ylabel(Best Cost) plt.title(ICA Convergence Curve) plt.grid(True, alpha0.3) plt.tight_layout() plt.show()这段代码把同化、竞争、最优记录串起来。注意每次同化后要重新计算殖民地成本并且检查殖民地是否超过宗主国如果超过就交换位置——这是 ICA 的「宗主国更新」机制保证每个帝国的宗主国始终是当前最强解。收敛曲线一般前期下降快后期平缓如果曲线一直不降检查同化系数和革命概率。3.2 帝国版图可视化用散点图看帝国收缩二维问题时可以把所有国家画在平面上宗主国用大点殖民地用小点不同帝国用不同颜色。每轮迭代画一帧就能看到帝国版图怎么收缩。def plot_empires(empires, ax, bounds): 在二维平面上画出当前帝国分布。 colors plt.cm.tab10(np.linspace(0, 1, len(empires))) for i, emp in enumerate(empires): if len(emp[colonies]) 0: ax.scatter(emp[colonies][:, 0], emp[colonies][:, 1], s15, colorcolors[i], alpha0.6, labelfEmpire {i}) ax.scatter(emp[imperialist][0], emp[imperialist][1], s120, colorcolors[i], marker*, edgecolorsblack) ax.set_xlim(bounds[0]) ax.set_ylim(bounds[1]) ax.set_title(fEmpires: {len(empires)}) ax.grid(True, alpha0.3) # 每隔 20 轮画一次 fig, axes plt.subplots(2, 3, figsize(15, 9)) axes axes.flatten() empires initialize_empires(200, 10, bounds, lambda x: penalty_objective(x, bounds)) for it in range(200): for emp in empires: emp assimilate(emp, 1.8, np.pi/4) emp[colony_costs] np.array([penalty_objective(c, bounds) for c in emp[colonies]]) empires imperialistic_competition(empires, 0.02) if it % 40 0 and it // 40 6: plot_empires(empires, axes[it // 40], bounds) plt.tight_layout() plt.show()散点图能直观看到帝国数量从 10 个逐渐减少到 1 个殖民地不断向宗主国聚拢。如果某个帝国的殖民地一直不收敛检查同化系数是不是太小或者革命概率是不是太低。3.3 殖民地迁移轨迹记录每一步的位置变化想看单个殖民地怎么走可以在同化函数里记录轨迹。下面是一个简化版只跟踪第一个帝国第一个殖民地的路径。def track_colony_path(empires, objective_func, bounds, max_iter100): 跟踪第一个帝国第一个殖民地的移动路径。 path [] emp empires[0] if len(emp[colonies]) 0: return path colony emp[colonies][0].copy() for _ in range(max_iter): path.append(colony.copy()) direction emp[imperialist] - colony theta np.random.uniform(-np.pi/4, np.pi/4) rot np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) direction[:2] rot direction[:2] colony colony 1.8 * np.random.rand() * direction # 越界拉回 for i, (low, high) in enumerate(bounds): colony[i] np.clip(colony[i], low, high) return np.array(path) path track_colony_path(empires, lambda x: penalty_objective(x, bounds), bounds) plt.figure(figsize(6, 6)) plt.plot(path[:, 0], path[:, 1], o-, markersize3, linewidth1) plt.scatter(empires[0][imperialist][0], empires[0][imperialist][1], s150, marker*, colorred, labelImperialist) plt.xlabel(x1) plt.ylabel(x2) plt.title(Colony Migration Path) plt.legend() plt.grid(True, alpha0.3) plt.show()轨迹图能看到殖民地不是直线冲向宗主国而是带随机偏转的折线。偏转角越大探索范围越广但收敛越慢。我一般先用大偏转角跑前期后期把偏转角调小做退火式衰减。4. 避坑与排查ICA 跑不出结果的五个血泪经验4.1 现象收敛曲线一开始就平了目标值几乎不变原因帝国初始化时势力值分配有问题或者同化系数太小殖民地几乎不动。常见做法是检查powers计算如果所有宗主国成本接近max_cost - cost会趋近于零归一化后势力值均匀殖民地分配没区分度。解决把势力值计算改成1 / (cost 1e-12)或者对成本做排序后按排名分配。同化系数从 1.8 起步跑 50 轮看曲线有没有下降。4.2 现象算法早熟所有殖民地挤在一个点原因革命概率太低或者偏转角太小殖民地缺乏多样性。ICA 没有变异操作全靠同化时的随机偏转和革命来维持探索。解决加革命操作每个帝国每轮以一定概率随机重置一个殖民地。革命概率取 0.05 到 0.1太高会破坏收敛太低会早熟。def revolution(empire, bounds, revolution_rate0.1): 革命操作以一定概率随机重置殖民地。 n_colonies len(empire[colonies]) n_revolve int(n_colonies * revolution_rate) if n_revolve 0: return empire idx np.random.choice(n_colonies, n_revolve, replaceFalse) for i in idx: empire[colonies][i] np.random.uniform( low[b[0] for b in bounds], high[b[1] for b in bounds] ) return empire4.3 现象帝国竞争后报错数组维度不匹配原因np.vstack拼接时如果最强帝国没有殖民地strongest[colonies]是空数组直接拼接会出问题。另外empires.pop之后索引会变如果还在循环里用旧索引会越界。解决拼接前判断殖民地是否为空空的话直接赋值而不是 vstack。竞争操作放在每轮迭代末尾操作完重新计算帝国数量不要在遍历过程中 pop。4.4 现象目标函数值出现 NaN 或 inf原因罚函数系数太大越界惩罚把数值撑爆或者目标函数本身有除零、log 负数。ICA 的势力值计算用了倒数成本为零时会出 inf。解决所有倒数计算加1e-12罚函数系数先取 1e3 试跑确认没有数值溢出再加大。目标函数里如果有除法分母加小量。4.5 现象可视化窗口卡死迭代 200 轮跑了十分钟原因每轮都重新画图matplotlib 的plt.show()阻塞主线程。或者目标函数计算太慢每次同化都全量重算。解决可视化每隔 20 到 50 轮画一次用plt.pause(0.01)代替plt.show()。目标函数如果计算量大考虑缓存或向量化numpy 的批量计算比 Python 循环快一个量级。5. 进阶技巧参数退火与多目标扩展的实操细节5.1 同化系数退火前期猛探索后期稳收敛固定同化系数有个矛盾大了前期探索好但后期震荡小了前期太慢。我一般用线性退火从 2.0 降到 0.5迭代 200 轮的话每轮减 0.0075。def run_ica_annealing(n_countries, n_imperials, bounds, objective_func, max_iter200): 带同化系数退火的 ICA。 empires initialize_empires(n_countries, n_imperials, bounds, objective_func) history [] for it in range(max_iter): # 同化系数从 2.0 线性降到 0.5 coef 2.0 - 1.5 * (it / max_iter) # 偏转角从 pi/3 降到 pi/12 angle np.pi/3 - (np.pi/3 - np.pi/12) * (it / max_iter) for emp in empires: emp assimilate(emp, coef, angle) emp[colony_costs] np.array([objective_func(c) for c in emp[colonies]]) if len(emp[colonies]) 0: best_idx np.argmin(emp[colony_costs]) if emp[colony_costs][best_idx] emp[imperialist_cost]: emp[imperialist], emp[colonies][best_idx] \ emp[colonies][best_idx].copy(), emp[imperialist].copy() emp[imperialist_cost], emp[colony_costs][best_idx] \ emp[colony_costs][best_idx], emp[imperialist_cost] empires imperialistic_competition(empires, 0.02) history.append(min(emp[imperialist_cost] for emp in empires)) return history退火的好处是前期殖民地大步跳覆盖更多区域后期小步微调避免在最优点附近震荡。实测在 Rastrigin 函数上退火版比固定系数版收敛精度高一个量级。5.2 多目标扩展用拥挤度距离替代单一势力值ICA 原生是单目标扩展到多目标需要改两处一是帝国势力值用 Pareto 支配关系加拥挤度距离二是殖民地与宗主国比较时用支配关系而不是单值比较。常见做法是 NSGA-II 的非支配排序加拥挤度套到 ICA 的帝国结构上。单目标 ICA多目标 ICA 扩展目标函数值排序非支配排序分层势力值 1/cost势力值 拥挤度距离殖民地与宗主国比大小殖民地与宗主国比支配关系帝国竞争按势力值帝国竞争按 Pareto 前沿质量多目标版的计算量比单目标大不少因为每轮都要做非支配排序。我一般把种群控制在 100 以内迭代 100 到 150 轮再大就跑不动了。5.3 验证方法用标准测试函数对比 GA/PSO写完 ICA 别急着上真实问题先用标准测试函数跑一遍和 GA、PSO 对比。我常用的三个函数Sphere单峰测收敛速度、Rastrigin多峰测全局搜索、Rosenbrock窄谷测方向搜索。每个函数跑 30 次独立实验记录最优值、均值、标准差。def sphere(x): return np.sum(x**2) def rosenbrock(x): return np.sum(100 * (x[1:] - x[:-1]**2)**2 (1 - x[:-1])**2) # 对比实验框架 functions {Sphere: sphere, Rastrigin: rastrigin, Rosenbrock: rosenbrock} bounds_2d [(-5.12, 5.12), (-5.12, 5.12)] for name, func in functions.items(): results [] for run in range(30): _, cost, _ run_ica_with_history( 200, 10, bounds_2d, lambda x: penalty_objective(x, bounds_2d, penalty_coef1e4), max_iter200 ) results.append(cost) print(f{name}: best{np.min(results):.4f}, fmean{np.mean(results):.4f}, std{np.std(results):.4f})跑完对比如果 ICA 在 Rastrigin 上明显优于 GA说明全局搜索能力到位如果在 Rosenbrock 上不如 PSO说明方向搜索精度不够可以调小后期同化系数。从那以后我每次改完 ICA 参数都强制跑一遍这三个函数的 30 次独立实验看均值和标准差有没有退化再上真实问题。希望帮到你。本文还有配套的精品资源点击获取
返回列表