一、原子轨道搜索算法
1、原子轨道搜索算法简介
原子轨道搜索算法（AOS），是Azizi M于2021年提出的一种新型智能优化算法，被收录在Applied Mathematical Modelling中。该算法主要是基于量子力学的一些原理以及量子的原子模型行为，考虑了原子核周围电子的特性。
（原文链接：https://www.sciencedirect.com/science/article/pii/S0307904X20307198；
原文PDF，无需积分：https://download.csdn.net/download/weixin_58427214/87360792）
2、原子轨道搜索算法的特点
① 寻优能力强。论文中处理不同数学测试函数的全局最优方面优于其他替代元启发式算法。
② 收敛速度快，能够以较低的所需函数评估收敛到更好的结果。
3、原子轨道搜索算法的主要步骤
原子轨道搜索算法实质上是三层循环。第一层循环是记录迭代的次数，第二层循环是对每一个虚拟层进行操作，第三层循环是根据光子速率（PR）对电子进行一定的扰动。
① 初始化候选解集
② 确定能量函数
③ 按照高斯分布分虚拟层
④ 计算虚拟层中候选解的结合态
、结合能
和最低能量原子的能量
⑤ 光子作用以及粒子或磁场的相互作用
⑥ 计算所有候选解的能量Es、结合态BS、结合能BE和最低能量原子的能量LE
原子轨道搜索的伪代码如下图所示。
二、原子轨道搜索算法的基本实现
（以下图函数为测试函数）
1、初始化候选解集
① 候选解的定义
② 初始解生成公式
③ 代码实现
根据候选解的数量、维度以及可行域的上下限随机生成初始的候选解集。
def initialization(m_candidates, p_dimensions, min_value, max_value):
"""
随机生成初始化候选解集
:param m_candidates: 候选解的数量
:param p_dimensions: 解的维度
:param min_value: 可行域的下限
:param max_value: 可行域的上限
:return: 初始化的候选解集
"""
candidates = []
for i in range(m_candidates):
Xi = []
for j in range(0, p_dimensions):
random_value = random.uniform(min_value, max_value)
Xi.append(random_value)
candidates.append(Xi)
return candidates
2、确定能量函数
① 能量函数的定义
基于量子的原子模型，每个电子都有一个能量状态，在数学模型中被认为是候选解的目标函数值。具有更好目标函数值的候选解代表具有较低能级的电子，而具有较高能级的电子代表具有较差目标函数值的候选解。
② 代码实现
经测试，将
为能量函数得到的结果较好，因此能量函数设置为
。
# 确定能量函数
def E(X):
"""
求一个候选解的能量
:param X: 候选解
:return: 该候选解的能量
"""
abs_sum = 0
abs_prod = 1
for x in X:
abs_sum += abs(x)
abs_prod *= abs(x)
result = abs_sum + abs_prod
return result / pow(10, 12)
3、按照高斯分布分虚拟层
① 虚拟层划分的原理
基于量子的原子模型，电子在原子核周围的位置由电子概率密度图确定，该电子概率密度图在数学模型中由概率密度函数(PDF)决定。根据概率论，变量的概率密度函数是表示该变量在特定范围内的可能性的函数。通过考虑在核周围假想创建的层，PDF用于确定候选解在这些层中的位置。候选解决方案按升序或降序排序（基于最小化或最大化优化问题），其中具有更好目标函数值的候选被认为具有更高的排名。
下图描绘了利用基于正态高斯分布的概率密度函数确定虚拟层中候选解（电子）位置的示意图。应该注意的是，在第二层中找到电子的总概率高于第一层，因此第二层（L1和L2之间）的概率密度值高于第一层（L0和L1之间）。
② 第k层虚拟层的候选解和能量的定义
③ 代码实现
通过高斯分布生成定义域内的数字，以确定候选解所在的虚拟层次。
# 通过高斯分布生成在定义域内的数字
def gauss(u, sigma, min_value, max_value):
"""
通过高斯分布生成在定义域内的数字
:param u: 均值
:param sigma: 标准差
:param min_value: 定义域的下限
:param max_value: 定义域的上限
:return: 高斯分布随机生成的数字
"""
random_value = random.gauss(u, sigma)
while random_value < min_value or random_value >= max_value:
random_value = random.gauss(u, sigma)
return random_value
按照高斯分布分虚拟层
# 按照高斯分布分虚拟层
def create_layers(Es, m_candidates, n):
"""
按照高斯分布分虚拟层
:param Es: 所有候选解的能量
:param m_candidates: 候选解的数量
:param n: 虚拟层的数量
:return: 所有候选解（原始下标）在虚拟层的分布
"""
# 将候选解的能量排序
Es_sorted = sorted(Es, reverse=True)
# 计算每个候选解所在层
count = [0] * n
for i in range(m_candidates):
# 根据电子的概率密度函数，且因为第二虚拟层的概率密度大于第一虚拟层，所以以1.5为均值
random_value = gauss(1.5, (n - 1.5) / 3, 0, n)
count[math.floor(random_value)] += 1
layers = {}
index = 0
for i in range(n):
layers[i] = []
for j in range(count[i]):
# 用在原始候选集中的下标标识
layers[i].append(list.index(Es, Es_sorted[index]))
index += 1
return layers
4、计算候选解的结合态、结合能和最低能量原子的能量
每个虚构层中具有最佳目标函数值的候选解是每个层中具有最低能级的电子(
)。此外，所有候选者中具有最佳目标函数值的候选解被认为是原子中具有最低能级（LE）的电子。结合能表示从其壳中去除电子所需的能量，在数学模型中根据每一层中候选解的位置和目标函数值确定。
① 第k层虚拟层的结合态和结合能的定义
② 所有候选解的结合态和结合能的定义
③ 代码实现
第k层虚拟层中结合态、结合能以及最低能量原子的能量在原子轨道搜索算法主函数中实现，未单独实现函数。
计算所有候选解的能量Es、结合态BS、结合能BE和最低能量原子的能量LE代码如下。
# 计算所有候选解的能量Es、结合态BS、结合能BE和最低能量原子的能量LE
def calculate_Es_BS_BE_LE(Xs, m_candidates, p_dimensions):
"""
计算所有候选解的能量Es、结合态BS、结合能BE和最低能量原子的能量LE
:param Xs: 候选解集
:param m_candidates: 候选解的数量
:param p_dimensions: 解的维度
:return: 所有候选解的能量Es、结合态BS、结合能BE和最低能量原子的能量LE
"""
# 计算候选解的能量
Es = []
for Xi in Xs:
Ei = E(Xi)
Es.append(Ei)
# 计算结合态
BS = np.zeros(p_dimensions)
for i in range(m_candidates):
BS = BS + np.array(Xs[i])
BS /= m_candidates
# 计算结合能
BE = sum(Es) / m_candidates
# 计算最低能量
LE = min(Es)
return Es, BS, BE, LE
5、光子作用以及粒子或磁场的相互作用
光子速率(PR)被确定为一个参数，表示考虑光子对电子的作用的概率。每个电子在(0,1)范围内生成一个均匀分布的随机数(φ)。如果每个电子的随机生成数(φ)大于PR(φ≥PR)，则光子对电子的作用是可能的，因此电子在原子核周围不同层之间的运动基于发射和吸收来考虑光子。如果每个电子的随机生成数 (φ)小于PR(φ<PR)，则光子对电子的作用是不可能的，因此电子在原子核周围不同层之间的运动是基于其他一些作用来考虑的，例如作为与粒子或磁场的相互作用，也会导致能量的吸收或发射。
① 发射光子
当每个电子的随机生成数(φ)大于PR(φ≥PR)时，电子在原子核周围不同层之间的运动基于发射和吸收来考虑光子。如果候选解的能量大于该虚拟层的结合能，则发射光子。
② 吸收光子
如果候选解的能量小于该虚拟层的结合能，则吸收光子。
③ 粒子或磁场的相互作用
当每个电子的随机生成数(φ)小于PR(φ<PR)时，与粒子或磁场的相互作用，也会导致能量的吸收或发射。
④ 代码实现
光子作用以及粒子或磁场的相互作用在原子轨道搜索算法主函数中实现，未单独实现函数。
6、原子轨道搜索算法主函数
# 原子轨道搜索算法
def AOS(m_candidates, p_dimensions, min_value, max_value, max_iterations, PR):
"""
AOS算法
:param m_candidates: 候选解的数量
:param p_dimensions: 候选解的维度
:param min_value: 可行域的下限
:param max_value: 可行域的上限
:param max_iterations: 最大迭代次数
:param PR: 光子速率
:return: 全局最低能量和全局最优候选解以及最低能量变化路径
"""
# 随机生成初始候选解集
Xs = initialization(m_candidates, p_dimensions, min_value, max_value)
# 计算所有候选解的能量Es、结合态BS、结合能BE和最低能量原子的能量LE
Es, BS, BE, LE = calculate_Es_BS_BE_LE(Xs, m_candidates, p_dimensions)
# 迭代过程中全局的最低能量、最优解和最低能量的变化
global_LE = LE
global_X = Xs[list.index(Es, LE)]
all_LE = []
all_LE.append(LE)
# 开始AOS迭代
count = 0
while count < max_iterations:
# 随机生成虚拟层的层数n
n = random.randint(2, 10)
# 将候选解分层
layers = create_layers(Es, m_candidates, n)
# 遍历每一层
for k in range(n):
# 计算第k虚拟层的结合态BSk、结合能BEk和最小能量LEk
p = len(layers[k])
if p <= 0:
continue
BSk = np.zeros(p_dimensions)
BEk = 0
LEk = 0xfffffff
for i in range(p):
BSk = BSk + np.array(Xs[layers[k][i]])
BEk += Es[layers[k][i]]
if Es[layers[k][i]] < LEk:
LEk = Es[layers[k][i]]
BSk /= p
BEk /= p
# 光子对电子进行作用
for i in range(p):
# 生成随机参数
phi, alpha, beta, gamma = random.random(), random.random(), random.random(), random.random()
# 若phi大于等于PR，则光子进行作用
if phi >= PR:
# 若大于该层结合能，则释放光子
if Es[layers[k][i]] >= BEk:
tmp = Xs[layers[k][i]] + \
(alpha * (beta * np.array([LE] * p_dimensions) - gamma * BS)) / (k + 1)
# 若小于该层结合能，则吸收光子
else:
tmp = Xs[layers[k][i]] + alpha * (beta * LEk - gamma * BSk)
# 若phi小于PR，则光子不进行作用
else:
random_vector = []
for j in range(p_dimensions):
random_vector.append(random.random())
tmp = Xs[layers[k][i]] + np.array(random_vector)
if in_range(tmp, min_value, max_value):
Xs[layers[k][i]] = tmp
# 计算当前候选解集的能量Es、结合态BS、结合能BE和最低能量原子的能量LE
Es, BS, BE, LE = calculate_Es_BS_BE_LE(Xs, m_candidates, p_dimensions)
count += 1
all_LE.append(LE)
# 判断是否为全局最优解
if LE < global_LE:
global_LE = LE
global_X = Xs[list.index(Es, LE)]
return global_LE, global_X, all_LE
7、问题求解
①参数设置
原子轨道搜索算法参数设置：循环迭代次数500，候选解集中候选解的数量100，光子速率0.1。
② 代码实现
# 原子轨道搜索算法求解
if __name__ == '__main__':
start_time = time.process_time()
best_e, best_x, all_ex = AOS(100, 30, -10, 10, 500, 0.1)
end_time = time.process_time()
print("函数的最小值为:", best_e * pow(10, 12))
print("优化求得的解为:", best_x)
print("总耗时:", end_time - start_time)
iteration = len(all_ex)
plt.xlabel("Iteration")
plt.ylabel("f2(x)")
plt.plot(range(iteration), all_ex)
plt.show()
③ 输出结果
函数的最小值为: 9.82747835036586
优化求得的解为: [ 0.26662175 0.42657184 0.51771326 0.25205068 -0.3727787 1.29653549
0.44019179 0.32447502 0.10145573 0.20620475 0.27570473 0.11991767
-0.71445652 -0.69569835 0.21877323 -0.15953558 0.15230652 0.17432502
-0.49591347 -0.20148722 0.24911851 0.20729576 -0.16116592 -0.2007519
0.63004922 -0.18744668 0.03416354 -0.29581042 -0.03657444 0.41238463]
总耗时: 1.0
（在课程大作业中，还在该测试函数上实现了模拟退火算法，并与原子轨道搜索算法进行测试对比。课程大作业报告及代码：https://download.csdn.net/download/weixin_58427214/87360935）