一、原子軌道搜尋算法
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）