SCE-UA算法主要流程步骤
SCE-UA算法（Shuffled Complex Evolution - University of Arizona）是一种高效的全局优化算法，特别适用于求解复杂、非线性、多峰（多极值）且参数较多的优化问题，在水文模型参数率定等领域应用广泛。其核心思想结合了确定性搜索（单纯形法）和随机性搜索（进化算法）的优点，并通过种群划分（复合体）、竞争进化和周期性洗牌重组来平衡全局探索和局部开发能力。
1. 初始化 (Initialization)
-
确定优化问题的维度
n（参数个数）、目标函数f(x)（需要最小化）、每个参数的可行域范围。 -
设定算法参数：
-
p: 复合体的数量。 -
q: 每个复合体包含的点数（通常q = 2 * n + 1）。 -
α: 每个复合体在内部进化过程中，每次调用CCE算法生成的子点数（通常α = 1）。 -
β: 每个复合体在洗牌前内部进化的次数（迭代次数，通常为β = 2 * n + 1）。 -
maxiter: 最大洗牌次数（全局迭代次数）。 -
kstop: 收敛判定所需的连续洗牌次数（未显著改进）。 -
pcento: 收敛判定的目标函数改进百分比阈值。
-
-
随机生成初始种群：在整个可行域内随机生成
s = p * q个样本点x₁, x₂, ..., xₛ。 -
计算目标函数值：评估每个样本点
xᵢ的目标函数值fᵢ = f(xᵢ)。
2. 洗牌与排序 (Shuffling and Sorting)
-
根据目标函数值
fᵢ将所有s个点按升序排序（最小化问题）：f(x₁) ≤ f(x₂) ≤ ... ≤ f(xₛ)，对应的点记为D = {u₁, u₂, ..., uₛ}（u₁是当前最优解）。 -
划分复合体 (Partition into Complexes)：
-
将排序后的点集 D 划分为 p 个复合体 C₁, C₂, …, Cₚ。
-
划分规则：采用“蛇形分配”确保每个复合体包含不同性能的点。
-
第1个点
u₁分配给C₁ -
第2个点
u₂分配给C₂ -
…
-
第
p个点uₚ分配给Cₚ -
第
p+1个点uₚ₊₁分配给Cₚ -
第
p+2个点uₚ₊₂分配给Cₚ₋₁ -
…
-
第
2p个点u₂ₚ分配给C₁ -
重复此模式，直到所有
s个点分配完毕。
-
-
每个复合体
Cₖ包含q个点：Cₖ = {vₖ₁, vₖ₂, ..., vₖ_q}，其中vₖ₁是该复合体内最好的点（函数值最小）。
-
3. 复合体竞争进化 (Competitive Complex Evolution - CCE)
-
对每个复合体
Cₖ (k = 1, 2, ..., p)独立并行地 进行β次进化迭代。 -
在每次进化迭代中，对复合体
Cₖ执行以下操作：-
选择父本子集 (Select Parents)：
- 从复合体
Cₖ的q个点中，根据某种概率分布（通常与目标函数值成反比，即偏好更好的点）随机选择m个不同的点构成一个 单纯形（通常m = n + 1）。这个单纯形作为生成新点的“父本”。
- 从复合体
-
生成子点 (Generate Offspring)：
-
对选出的
m个父本点应用 单纯形进化操作，生成 α 个新的试验点（子点）。主要操作包括：-
反射 (Reflection): 远离最差点（单纯形中函数值最大的点）的方向。
-
收缩 (Contraction): 朝向最好点（单纯形中函数值最小的点）的方向。
-
扩展 (Expansion): 沿反射方向走得更远（如果反射点表现好）。
-
随机化 (Randomization): 若上述操作失败，在父本点附近或整个可行域内随机生成点。
-
-
每次操作后需确保新点满足参数边界约束，若不满足则需处理（如拉回边界、重新生成等）。
-
-
评估与替换 (Evaluate and Replace)：
-
计算每个新子点的目标函数值。
-
用新子点替换掉复合体
Cₖ中的 最差点（函数值最大的点），但需满足 新子点比该最差点更优。如果不满足，则放弃该子点（或尝试其他替换策略）。 -
更新复合体
Cₖ内部点的排序（将新点插入正确位置）。
-
-
4. 重组复合体 (Recombine Complexes)
-
完成所有
p个复合体的β次内部进化迭代后，将 所有复合体中的所有点 重新合并回一个大的点集D'（大小仍为s = p * q）。 -
跳转回 步骤2：洗牌与排序，根据新的目标函数值重新排序
D'并按照蛇形分配规则重新划分成p个新的复合体。这个过程称为 洗牌 (Shuffling)。
5. 检查收敛 (Check Convergence)
-
重复 步骤2 (洗牌与排序) -> 步骤3 (CCE进化) ->步骤4 (重组) 的过程，直到满足以下任一停止条件：
-
达到预设的最大洗牌次数
maxiter。 -
在最近的
kstop次洗牌过程中，种群最优解的目标函数值改进小于pcento%（即连续kstop次洗牌后，最优值的变化百分比小于阈值）。 -
找到满足精度要求的最优解（函数值足够小）。
-
-
如果未收敛，则进行下一次洗牌迭代。
6. 输出结果 (Output Results)
当算法停止时，输出整个进化过程中找到的 最佳点 u₁ 及其对应的目标函数值 f(u₁) 作为优化结果。
关键特点总结
-
** 种群划分与洗牌：** 将大种群划分为多个子群（复合体），子群独立进化（局部开发），然后洗牌重组（全局探索和信息交换），有效平衡探索与开发。
-
** 竞争进化（CCE）：** 每个复合体内部使用基于单纯形的进化操作（反射、收缩等）高效地生成新点并改进解的质量，利用了确定性搜索的快速收敛性。
-
** 复合体多样性：** “蛇形分配”确保每个复合体都包含好、中、差不同质量的点，提高了进化的多样性和鲁棒性。
-
** 全局性：** 周期性洗牌重组机制有效防止了算法过早陷入局部最优，增强了全局搜索能力。
-
** 并行潜力：** 不同复合体之间的进化是独立的，天然适合并行计算。
简而言之，SCE-UA通过不断地将种群打散重组（洗牌）成多个子群（复合体），并在每个子群内部进行高效的局部搜索（CCE进化），最终引导整个种群向全局最优解逼近。这种结构是其处理高维复杂优化问题成功的关键。
具体代码实现
# -*- coding: utf-8 -*-
"""
@file: sceua_v1.py
@author: 温新
@date: 2026/3/11
@description: sceua算法的实现
"""
import numpy as np
class SCEUA:
def __init__(self,
obj_func,
bounds,
ngs=10,
max_call=5000,
max_loop=None,
kstop=100,
pcento=0.1,
peps=1e-6,
seed=None):
"""
SCE-UA算法的实现
:param obj_func: 目标函数，接受 shape=(d,) 的数组，返回标量（越小越好）
:param bounds: [(low1, high1), (low2, high2), ..., (low_d, high_d)]
:param ngs: 子复形的数量, 默认为 10 (典型值：5 ~ 20)
:param max_call: 最大目标函数评估次数
:param max_loop: 最大循环数
:param kstop: 过去进化循环的数量及其各自的目标值，用于评估当前循环的边际改进(百分比)是否小于pcento
:param pcento: 在过去的kstop循环中允许的百分比变化，低于该百分比变化，则认为达到了收敛。
:param peps: 总体中参数的标准化几何范围的值，低于该值则认为达到了收敛。
:param seed: 随机种子
"""
self.obj_func = obj_func
if seed is not None:
np.random.seed(seed)
self.bounds = np.array(bounds)
self.bl = self.bounds[:, 0]
self.bu = self.bounds[:, 1]
self.ngs = ngs # 子复形的数量
self.nopt = self.bounds.shape[0] # 维度，待优化参数的数量
self.npg = 2 * self.nopt + 1 # 每个复形的点数
self.nspl = self.npg # 单个复合形的演化次数
self.nps = self.nopt + 1 # 子复合形(单纯形)的点数
self.npt = self.npg * self.ngs # 总点数
self.max_call = max_call
self.max_loop = max_loop
self.kstop = kstop
self.pcento = pcento
self.peps = peps
self.stochastic_parameters = (self.bu-self.bl) != 0 # 得到一个bool列表，固定参数位置为False
self.x = None # 种群
self.xf = None # 适应度值
self.best_x = None # 当前最优参数
self.best_f = None # 当前最优适应度值
self.history = [] # 记录历史值
self.icall = 0
def _init_population(self):
"""
初始化种群
"""
self.x = np.random.uniform(self.bl, self.bu, size=(self.npt, self.nopt))
self.xf = np.array([self.obj_func(ind) for ind in self.x])
self.icall = self.npt
def _gnrng(self):
""" 计算归一化几何范围 """
if not np.any(self.stochastic_parameters):
return 0.0 # 所有参数固定，视为完全收敛
x_sub = self.x[:, self.stochastic_parameters] # shape: (npt, d_stoch)
bl_sub = self.bl[self.stochastic_parameters]
bu_sub = self.bu[self.stochastic_parameters]
ranges = np.max(x_sub, axis=0) - np.min(x_sub, axis=0) # shape: (d_stoch,)
normalized_ranges = ranges / (bu_sub - bl_sub)
# 防止除零和 log(0)：将极小值截断到 1e-12
normalized_ranges = np.clip(normalized_ranges, 1e-12, 1.0)
log_mean = np.mean(np.log(normalized_ranges))
return np.exp(log_mean)
def simulate(self):
# 初始化种群
self._init_population()
# 1. 按适应度值排序
idx_sorted = np.argsort(self.xf)
self.x = self.x[idx_sorted]
self.xf = self.xf[idx_sorted]
# 记录当前最优
self.best_x = self.x[0]
self.best_f = self.xf[0]
self.history.append(self.best_f)
# 初始化收敛指标
gnrng = self._gnrng()
nloop = 0
criter_change_pcent = 1e-5
while self.icall < self.max_call:
nloop += 1
print(f"复形演化循环 #{nloop} 正在进行中...")
# 2. 划分为 self.ngs 个复合形
cs = []
for i in range(self.ngs):
indices_complex = list(range(i, self.npt, self.ngs))
cx = self.x[indices_complex].copy()
cf = self.xf[indices_complex].copy()
cs.append((cx, cf))
# 3. 对每个复合形独立演化
new_x = []
new_f = []
for cx, cf in cs:
cxnew, cfnew, icalls = self._complex_evolution(cx, cf)
new_x.append(cxnew,)
new_f.append(cfnew)
self.icall += icalls
# 合并所有复合形
self.x = np.vstack(new_x)
self.xf = np.hstack(new_f)
# 4. 洗牌，打乱顺序，打破分组记忆
idx = np.argsort(self.xf)
self.xf = self.xf[idx]
self.x = self.x[idx]
# 记录最好点
self.best_x = self.x[0]
self.best_f = self.xf[0]
self.history.append(self.best_f)
gnrng = self._gnrng() # 计算归一化几何范围
# 检查收敛
if self.icall >= self.max_call:
print(f"*** 优化搜索因为超过最大尝试次数{self.max_call}而中断。")
print(f"搜索在初始循环的尝试次数{self.icall}停止！")
break
elif gnrng < self.peps:
print("种群已经收敛到预先指定的小参数空间")
break
elif nloop > self.kstop:
print("现在正在更新和评估目标函数收敛标准...")
his = np.array(self.history)
absolute_change = (np.abs(his[nloop-1] - his[nloop-self.kstop]))
denominator = np.mean(np.abs(his[(nloop - self.kstop): nloop]))
if denominator == 0:
criter_change_pcent = 0.0
else:
criter_change_pcent = absolute_change / denominator * 100
print(f"更新后的收敛标准：{criter_change_pcent}")
if criter_change_pcent < self.pcento:
print(f"最佳点在最后{self.kstop}个循环中的改进小于用户指定的阈值{self.pcento}")
print("基于目标函数准则实现了收敛！！！")
break
elif self.max_loop and nloop >= self.max_loop:
print(f"已达到每次执行的最大循环数")
break
print(f"搜索在尝试{self.icall}次后停止。" )
print(f"归一化几何范围 = {gnrng}")
print(
f"最佳点在最后{self.kstop}个循环中变化小于{criter_change_pcent}"
)
# 最终排序， 返回最优解
idx_sorted = np.argsort(self.xf)
self.best_x = self.x[idx_sorted[0]]
self.best_f = self.xf[idx_sorted[0]]
return self.best_x, self.best_f, self.history
def _complex_evolution(self, cx, cf):
"""
单个复合形的演化
:param cx: np.array shape：[self.npg, self.nopt] 子复合形的所有点
:param cf: np.array shape: [self.npg] 子复合形点的适应度值
:return: 进化后的子复合形 （cx, cf, icalls）
"""
icalls = 0
for loop in range(self.nspl):
# 取样
lcs = self._sample_simplex_idx()
# 构造单纯形
s = cx[lcs]
sf = cf[lcs]
snew, fnew, n_evals = self._cceua(s, sf)
# 用新点取代单纯形中最差的点
s[-1] = snew
sf[-1] = fnew
# 把单纯形重新放回复合形
cx[lcs] = s
cf[lcs] = sf
# 对复合形进行排序
idx = np.argsort(cf)
cf = cf[idx]
cx = cx[idx]
icalls += n_evals
return cx, cf, icalls
def _sample_simplex_idx(self):
"""
根据线性概率分布对复合体进行采样来选择单纯形
:return: np.array shape:(self.nps,) 索引数组
"""
lcs = np.array([0] * self.nps)
lcs[0] = 0
for k in range(1, self.nps):
for _ in range(1000):
lpos = int(
np.floor(
self.npg + 0.5 - np.sqrt((self.npg + 0.5) ** 2 - self.npg * (self.npg + 1) *np.random.random())
)
)
# 确保lpos在有效范围内
lpos = min(lpos, self.npg-1)
# 检查元素是否已经被选择
if lpos not in lcs[:k]:
lcs[k] = lpos
break
else:
# 如果1000次都失败，随机选择一个未使用的点
available = set(range(self.npg)) - set(lcs[:k])
lcs[k] = np.random.choice(list(available))
lcs.sort()
return lcs
def _cceua(self, s, sf):
"""
这是在单纯形中生成新点的子程序
:param s: the sorted simplex in order of increasing function values
:param sf: function values in increasing order
LIST OF LOCAL VARIABLES
sb(.) = the best point of the simplex
sw(.) = the worst point of the simplex
fw = function value of the worst point
ce(.) = the centroid of the simplex excluding wo
snew(.) = new point generated from the simplex
"""
constant_parameters = np.invert(self.stochastic_parameters)
alpha = 1.0
beta = 0.5
n_evals = 0
# 找最差点(因为已经排序，所以最后一个最差)
sw = s[-1]
fw = sf[-1]
# 计算其余点的质心（除去最差点）
ce = np.mean(s[:-1], axis=0)
# 尝试反射点
snew = ce + alpha * (ce - sw)
snew[constant_parameters] = sw[constant_parameters]
if np.all(snew >= self.bl) and np.all(snew <= self.bu): # 检查是否越界
# 反射点有效
fnew = self.obj_func(snew)
n_evals += 1
if fnew > fw:
# 如果反射失败，则尝试收缩
snew = sw + beta * (ce - sw)
snew[constant_parameters] = sw[constant_parameters]
fnew = self.obj_func(snew)
n_evals += 1
# 如果反射和收缩全都失败，尝试随机点
if fnew > fw:
snew = np.random.uniform(self.bl, self.bu, size=self.nopt)
snew[constant_parameters] = sw[constant_parameters] # 保留固定参数
fnew = self.obj_func(snew)
n_evals += 1
else:
# 反射点无效，尝试收缩
snew = sw + beta * (ce - sw)
snew[constant_parameters] = sw[constant_parameters]
fnew = self.obj_func(snew)
n_evals += 1
# 如果反射和收缩全都失败，尝试随机点
if fnew > fw:
snew = np.random.uniform(self.bl, self.bu, size=self.nopt)
snew[constant_parameters] = sw[constant_parameters] # 保留固定参数
fnew = self.obj_func(snew)
n_evals += 1
return snew, fnew, n_evals
以上代码为作者所写，难免存在问题，欢迎各位讨论交流。
扫一扫
383
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
