@
目錄
- 前言
- 1.測試函式
- 1.1 針狀函式
- 1.1.1 表示式
- 1.1.2 特徵
- 1.1.3 影像
- 1.2 Brains’s rcos函式
- 1.2.1 表示式
- 1.2.2 特徵
- 1.2.3 影像
- 1.3 Griewank函式
- 1.3.1 表示式
- 1.3.2 特徵
- 1.3.3 影像
- 1.4 Easom’s函式
- 1.4.1 表示式
- 1.4.2 特徵
- 1.4.3 影像
- 1.5 Schwefel’s函式
- 1.5.1 表示式
- 1.5.2 特徵
- 1.5.3 影像
- 1.1 針狀函式
- 2. PSO(Particle swarm optimization)
- 2.1 解決的問題
- 2.2 API
- 2.3 引數
- 2.4 示例
- 2.4.1 不帶約束條件
- 2.4.2 帶約束條件
- 3.Genetic Algorithm(recommended)
- 3.1 解決的問題
- 3.2 API
- 3.2.1 Genetic Algorithm
- 3.2.2 Genetic Algorithm for TSP(Travelling Salesman Problem)
- 3.3 TSP的目標函式定義的示例
- 3.4 示例
- 3.5.1 目標函式最佳化
- 3.5.2 TSP
- 4.Differential Evolution
- 4.1 解決的問題
- 4.2 API
- 4.3 引數
- 4.4 示例
- 特別篇：GA與DE的區別
- 5.SA(Simulated Annealing)
- 5.1 解決的問題
- 5.2 API
- 5.2.1 Simulated Annealing
- 5.2.2 SA for TSP
- 5.3 示例
- 5.3.1 目標函式最佳化
- 5.3.2 TSP
- 5.4 bug
- 6.ACA (Ant Colony Algorithm) for tsp
- 6.1 解決的問題
- 6.2 API
- 6.3 引數
- 6.4 示例
- 7.immune algorithm (IA)
- 7.1 解決的問題
- 7.2 API
- 7.3 引數
- 7.4 示例
- 8.Artificial Fish Swarm Algorithm (AFSA)
- 8.1 解決的問題
- 8.2 API
- 8.3 引數
- 8.4 示例
前言
本文著力於介紹scikit-opt工具包中七大啟發式演算法的API呼叫方法，關於具體的數學原理和推導過程，本文不再介紹，請讀者自行查詢相關文獻。
1.測試函式
為了檢驗這些啟發式演算法的效果，本文使用了以下五種專門用於測試的函式。
1.1 針狀函式
1.1.1 表示式
$$
f(r)=\frac{\sin(r)}{r}+1,r=\sqrt{(x-50){2}+(y-50){2}}+e
$$
$$
0\le x,y\le100
$$
1.1.2 特徵
該函式是一個多峰函式，在(50,50)處取得全域性最大值1.1512，第二最大值在其全域性最大值附近，採用一般的最佳化方法很容易陷入區域性極大值點。
1.1.3 影像
1.2 Brains’s rcos函式
1.2.1 表示式
$$
f(x_{1},x_{2})=a(x_{2}-bx_{1}{2}+cx_{1}-d)+e(1-f)\cos (x_{1})+e
$$
$$
a=1,b=\frac{5.1}{4\pi ^{2}},c=\frac{5}{\pi },d=6,e=10,f=\frac{1}{8 \pi }
$$
1.2.2 特徵
有3個全域性最小值：0.397887，分別位於 (-pi, 12.275), (pi, 2.275), (9.42478, 2.475)。
還有一個區域性最小值為：0.8447
1.2.3 影像
1.3 Griewank函式
1.3.1 表示式
$$
f(x)=\sum_{i=1}{N}\frac{x_{i}{2}}{4000}-\prod_{i=1}^{N}\cos (\frac{x_{i}}{\sqrt[]{i}})+1
$$
$$
-600\le x_{i} \le 600,i=1,2,...,N
$$
1.3.2 特徵
在(0,0,…,0)處取得全域性極小值0
1.3.3 影像
1.4 Easom’s函式
1.4.1 表示式
$$
f(x,y)=-\cos (x) \cos (y)\exp (-(x-\pi)^{2}-(y- \pi)^{2})
$$
$$
-100 \le x,y \le 100
$$
1.4.2 特徵
在(pi,pi)處取得全域性最小值-1。
1.4.3 影像
1.5 Schwefel’s函式
1.5.1 表示式
$$
f(x)=\sum_{i=1}^{N}-x_{i}\sin (\sqrt{\left | x_{i} \right | } )
$$
$$
-500 \le x_{i} \le 500,i=1,2,...,N
$$
1.5.2 特徵
$$
x_{i}=420.9687時有一個全域性最小值-n\times 418.9839
$$
1.5.3 影像
2. PSO(Particle swarm optimization)
2.1 解決的問題
- 目標函式最佳化
2.2 API
class PSO(SkoBase):
def __init__(self, func, n_dim=None, pop=40, max_iter=150, lb=-1e5, ub=1e5, w=0.8, c1=0.5, c2=0.5,
constraint_eq=tuple(), constraint_ueq=tuple(), verbose=False
, dim=None):...
2.3 引數
| 引數 | 含義 |
|---|---|
| func | 欲最佳化的目標函式 |
| n_dim | 所給函式的引數個數 |
| pop | 粒子群規模 |
| max_iter | 最大迭代次數 |
| lb(lower bound) | 函式引數的下限 |
| ub(upper bound) | 函式引數的上限 |
| w | 慣性因子 |
| c1 | 個體學習因子 |
| c2 | 社會學習因子 |
| constraint_eq | 等式約束條件，適用於非線性約束 |
| constraint_ueq | 不等式約束條件，適用於非線性約束 |
引數詳解：
- 定義目標函式func時，引數有且只有一個
- 慣性權重w描述粒子上一代速度對當前代速度的影響。w 值較大，全域性尋優能力強，區域性尋優能力弱；反之，則區域性尋優能力強。 當問題空間較大時，為了在搜尋速度和搜尋精度之間達到平衡，通常做法是使演算法在前期有較高的全域性搜尋能力以得到合適的種子，而在後期有較高的區域性搜尋能力以提高收斂精度。所以w不宜為一個固定的常數。
- 個體學習因子c1 =0稱為無私型粒子群演算法，即“只有社會，沒有自我”，會迅速喪失群體多樣性，容易陷入區域性最優解而無法跳出；社會學習因子c2=0稱為自我認識型粒子群演算法，即“只有自我，沒有社會”，完全沒有資訊的社會共享，導致演算法收斂速度緩慢； c1，c2都不為0，稱為完全型粒子群演算法，完全型粒子群演算法更容易保持收斂速度和搜尋效果的均衡，是較好的選擇。通常設定c1=c2=2，但不一定必須等於2，通常c1=c2∈[0,4]。
- 群體大小pop是一個整數，pop很小時陷入區域性最優解的可能性很大；pop很大時PSO的最佳化能力很好, 但是當群體數目增長至一定水平時，再增長將不再有顯著作用，而且數目越大計算量也越大。群體規模pop一般取20～40，對較難或特定類別的問題 可以取到100～200。
- 定義約束條件時，需要用元組，元組中用匿名函式
2.4 示例
2.4.1 不帶約束條件
# 匯入必要的庫，定義函式
import sko.PSO as PSO
import numpy as np
import matplotlib.pyplot as plt
def Schwefel(x):
'''
Schwefel's 函式，自變數為一個n維向量
該向量的每一個分量 -500<=x(i)<=500
當自變數的每一個分量的值為420.9687，有一個全域性最小值為 -n*418.9839
'''
# 初始化函式值 z 為0
z = 0
# 遍歷每個分量
for i in range(1, x.size+1):
# 計算函式值 z
z = z-x[i-1]*np.sin(np.sqrt(np.abs(x[i-1])))
# 返回目標函式值
return z
# 使用粒子群演算法求解 Schwefel 函式的最小值
# 建立 PSO 物件 ex_pso，維度為 3，種群規模為 400，最大迭代次數為 200
# 變數的上下界為-500到500，慣性權重為1，個體和全域性學習引數 c1 和 c2 都為2
ex_pso=PSO.PSO(func=Schwefel,n_dim=3,pop=400,
max_iter=200,lb=[-500,-500,-500,],
ub=[500,500,500,],w=1,c1=2,c2=2)
# 執行 PSO 演算法
ex_pso.run()
# 輸出找到的最佳解
print('best_x is ', ex_pso.gbest_x, 'best_y is', ex_pso.gbest_y)
# 視覺化迭代過程
# 繪製歷史最佳解的變化曲線
plt.plot(ex_pso.gbest_y_hist)
# 儲存並顯示影像
plt.savefig("PSO for Schwefel.png",dpi=300)
plt.show()
best_x is [425.69669626 421.40732833 425.7049067 ] best_y is [-1251.26934821]
2.4.2 帶約束條件
# 匯入必要的庫，定義目標函式和約束條件
import sko.PSO as PSO
import numpy as np
import matplotlib.pyplot as plt
def Schwefel(x):
'''
Schwefel's 函式，自變數為一個n維向量
該向量的每一個分量 -500=<x(i)<=500
當自變數的每一個分量的值為420.9687
有一個全域性最小值為 -n*418.9839
'''
# 初始化目標函式值 z 為0
z = 0
# 遍歷每個分量
for i in range(1, x.size+1):
# 計算函式值 z
z = z-x[i-1]*np.sin(np.sqrt(np.abs(x[i-1])))
# 返回目標函式值
return z
my_constrain_ueq=(
# 定義不等式約束條件，表示變數之間的關係
lambda x:(x[0]-420)**2 + (x[1]-420)**2+ (x[2]-420)**2 - 500**2 ,
)
# 使用粒子群演算法求解帶有不等式約束的 Schwefel 函式的最小值
# 建立 PSO 物件 ex_pso，維度為 3，種群規模為 400，最大迭代次數為 200
# 變數的上下界為-500到500，慣性權重為1，個體和全域性學習引數 c1 和 c2 都為2
# 定義不等式約束
ex_pso=PSO.PSO(func=Schwefel,n_dim=3,pop=400,
max_iter=200,lb=[-500,-500,-500,],
ub=[500,500,500,],w=1,c1=2,c2=2,
constraint_ueq=my_constrain_ueq)
# 執行 PSO 演算法
ex_pso.run()
# 輸出找到的最佳解
print('best_x is ', ex_pso.gbest_x, 'best_y is', ex_pso.gbest_y)
# 視覺化迭代過程
# 繪製歷史最佳解的變化曲線
plt.plot(ex_pso.gbest_y_hist)
# 儲存並顯示影像
plt.savefig("PSO for Schwefel.png",dpi=300)
plt.show()
best_x is [421.25075001 421.6301641 422.22992202] best_y is [-1256.68262677]
- 可以看到收斂更快了
3.Genetic Algorithm(recommended)
3.1 解決的問題
- 目標函式最佳化
- TSP(Travelling Salesman Problem)問題
3.2 API
3.2.1 Genetic Algorithm
class GA(GeneticAlgorithmBase):
"""genetic algorithm
Parameters
----------------
func : function
The func you want to do optimal
n_dim : int
number of variables of func
lb : array_like
The lower bound of every variables of func
ub : array_like
The upper bound of every variables of func
constraint_eq : tuple
equal constraint
constraint_ueq : tuple
unequal constraint
precision : array_like
The precision of every variables of func
size_pop : int
Size of population
max_iter : int
Max of iter
prob_mut : float between 0 and 1
Probability of mutation
"""
def __init__(self, func, n_dim,
size_pop=50, max_iter=200,
prob_mut=0.001,
lb=-1, ub=1,
constraint_eq=tuple(), constraint_ueq=tuple(),
precision=1e-7, early_stop=None):...
| 引數 | 含義 |
|---|---|
| func | 欲最佳化的目標函式 |
| n_dim | 所給函式的引數個數 |
| size_pop | 種群初始規模 |
| max_iter | 最大迭代次數 |
| prob_mut | 變異機率 |
| lb(lower bound) | 函式引數的下限 |
| ub(upper bound) | 函式引數的上限 |
| constraint_eq | 等式約束條件，適用於非線性約束 |
| constraint_ueq | 不等式約束條件，適用於非線性約束 |
| precision | 精度要求 |
調整這些引數可以對演算法的效能產生影響，具體如下：
- 調整
size_pop：增加種群數量，可以增加搜尋空間，但會導致計算代價增加。 - 調整
max_iter：增加迭代次數可以增加搜尋空間，但會增加執行時間。 - 調整
prob_mut：增加個體基因突變的機率，可以增加種群的多樣性，但會降低穩定性，易陷入區域性最優解。 - 調整
precision：增加精度可以增加結果的精確度，但會增加計算時間和計算量。 - 調整約束條件：新增約束條件可以使得解更符合實際限制，但會增加計算複雜度。
- 調整
early_stop：合理設定提前終止演算法的條件可以減少不必要的執行時間。
3.2.2 Genetic Algorithm for TSP(Travelling Salesman Problem)
class GA_TSP(GeneticAlgorithmBase):
"""
Do genetic algorithm to solve the TSP (Travelling Salesman Problem)
Parameters
----------------
func : function
The func you want to do optimal.
It inputs a candidate solution(a routine), and return the costs of the routine.
size_pop : int
Size of population
max_iter : int
Max of iter
prob_mut : float between 0 and 1
Probability of mutation
"""
def __init__(self, func, n_dim, size_pop=50, max_iter=200, prob_mut=0.001):...
| 引數 | 含義 |
|---|---|
| func | 欲最佳化的目標函式 |
| n_dim | 所給函式的引數個數 |
| size_pop | 種群初始規模 |
| max_iter | 最大迭代次數 |
| prob_mut | 變異機率 |
3.3 TSP的目標函式定義的示例
首先，我們來講解一下TSP問題是什麼：TSP（Traveling Salesman Problem）指的是給定一個包含N個城市的地圖，每對城市之間都有一個距離。現在需要找到一條路徑，使得在該路徑上每個城市僅僅只被經過一次，並且要求走過的總路程最小化。
那麼，為了解決這個問題，需要先生成N個隨機的城市，並計算它們之間的距離。（每個城市用一個座標點表示，計算距離可以使用歐幾里得距離公式）
針對上述問題，程式碼的流程可以描述如下：
- 首先，生成隨機的城市座標。
- 接著，根據這些城市座標計算它們之間的距離。
- 最後，定義一個目標函式，計算按某種順序依次經過這些城市的總路程。
按照這個思路，我們可以寫出如下的程式碼：
import numpy as np
# 隨機生成50個兩維點座標
num_points = 50
points_coordinate = np.random.rand(num_points, 2)
def get_distance(x1, x2):
"""計算兩個點之間的歐幾里得距離"""
return np.sqrt(np.sum(np.square(x1 - x2)))
# 根據城市座標計算各個城市之間的距離
n = points_coordinate.shape[0]
distance_matrix = np.zeros((n, n))
for i in range(n):
for j in range(n):
if i != j:
distance_matrix[i, j] = get_distance(points_coordinate[i], points_coordinate[j])
def cal_total_distance(routine):
'''The objective function. input routine, return total distance.
cal_total_distance(np.arange(num_points))
'''
# 按順序經過所有城市的總路程
return np.sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
# 測試
print(cal_total_distance(np.arange(num_points)))
下面還有一個官方給出的版本，比較簡潔：
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
num_points = 50
points_coordinate = np.random.rand(num_points, 2) # generate coordinate of points
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
def cal_total_distance(routine):
'''The objective function. input routine, return total distance.
cal_total_distance(np.arange(num_points))
'''
num_points, = routine.shape
return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
#測試
print(cal_total_distance(np.arange(num_points)))
3.4 示例
3.5.1 目標函式最佳化
# 匯入需要的庫
import sko
import matplotlib.pyplot as plt
import numpy as np
# 構造目標函式
def pin_func(p):
'''
一個針狀函式，取值範圍為 0 <= x, y <= 100
有一個最大值為1.1512在座標(50,50)
'''
# 解析 p 向量，即 (x, y)
x, y = p
# 計算針狀函式的結果
r = np.square( (x-50)** 2 + (y-50)** 2 )+np.exp(1)
t = np.sin(r) / r + 1
# 轉化為求最小值
return -1 * t
# 使用遺傳演算法求解目標函式的最小值
# 構造遺傳演算法物件 ex_GA，大小為50，迭代次數為800次
ex_GA = sko.GA.GA(func=pin_func, n_dim=2, size_pop=50, max_iter=800, prob_mut=0.001,
lb=[0, 0,], ub=[100, 100,], precision=1e-7)
# 執行遺傳演算法
best_x, best_y = ex_GA.run()
# 輸出找到的最優解
print('best_x:', best_x, '\n', 'best_y:', best_y)
# 視覺化迭代過程
# 獲取歷史搜尋最佳解序列
Y_history = ex_GA.all_history_Y
# 建立畫布物件和子圖物件
fig, ax = plt.subplots(2, 1)
# 繪製搜尋的目標函式值的歷史曲線
ax[0].plot(range(len(Y_history)), Y_history, '.', color='red')
# 計算歷史最小值的累積曲線
min_history = [np.min(Y_history[:i+1]) for i in range(len(Y_history))]
# 繪製歷史最小值的累積曲線
ax[1].plot(range(len(min_history)), min_history, color='blue')
# 儲存並顯示影像
plt.savefig("GA.png",dpi=400)
plt.show()
- 經過實測，GA對目標函式的最佳化不如PSO和DE
3.5.2 TSP
# 匯入必要的庫
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
# 隨機生成點的數目
num_points = 50
# 生成隨機座標陣列
points_coordinate = np.random.rand(num_points, 2)
# 生成距離矩陣，即每對點之間的歐氏距離
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
# 定義目標函式
def cal_total_distance(routine):
'''計算旅行商問題的總距離，輸入路線 routine，返回總距離
範例：cal_total_distance(np.arange(num_points))
'''
# 獲取路線上的點的數目
num_points, = routine.shape
# 計算每對相鄰點之間的距離之和
return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
# 使用遺傳演算法求解旅行商問題
import sko
# 建立 GA_TSP 類物件 ex_ga_tsp
ex_ga_tsp = sko.GA.GA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=500, prob_mut=1)
# 執行 GA 演算法
best_points, best_distance = ex_ga_tsp.run()
# 將結果進行視覺化展示
# 繪製找到的最佳路線
fig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
# 繪製迭代過程中歷代最佳總距離的變化曲線
ax[1].plot(ex_ga_tsp.generation_best_Y)
# 儲存並顯示影像
plt.savefig("GA_TSP.png",dpi=400)
plt.show()
4.Differential Evolution
4.1 解決的問題
- 目標函式最佳化
4.2 API
class DE(GeneticAlgorithmBase):
def __init__(self, func, n_dim, F=0.5,
size_pop=50, max_iter=200, prob_mut=0.3,
lb=-1, ub=1,
constraint_eq=tuple(), constraint_ueq=tuple()):...
4.3 引數
| 引數 | 含義 |
|---|---|
| func | 欲最佳化的目標函式 |
| n_dim | 所給函式的引數個數 |
| size_pop | 種群初始規模 |
| max_iter | 最大迭代次數 |
| prob_mut | 變異機率 |
| lb(lower bound) | 函式引數的下限 |
| ub(upper bound) | 函式引數的上限 |
| constraint_eq | 等式約束條件，適用於非線性約束 |
| constraint_ueq | 不等式約束條件，適用於非線性約束 |
調整這些引數會對演算法的收斂速度、效能和精度產生影響：
F的調整會影響演算法的收斂速度和效能，在實際應用中，通常取值範圍在[0.5, 1]之間，較小的值會降低演算法的收斂速度，較大的值會提高演算法的效能但易於陷入區域性最優。size_pop的調整會影響演算法的搜尋範圍和精度，在實際應用中，通常取值範圍在[10, 100]之間，較小的值會降低演算法的搜尋範圍和精度，較大的值會提高演算法的搜尋範圍但易於陷入區域性最優。max_iter的調整會影響演算法的收斂速度和效能，在實際應用中，通常取值範圍在[50, 1000]之間，較小的值會降低演算法的收斂速度，較大的值會提高演算法的效能但易於陷入區域性最優。prob_mut的調整會影響演算法的搜尋範圍和精度，在實際應用中，通常取值範圍在[0.2, 0.5]之間，較小的值會降低演算法的搜尋範圍和精度，較大的值會提高演算法的搜尋範圍但易於陷入區域性最優。lb和ub的調整會影響演算法的搜尋範圍和精度，較寬的取值範圍會提高演算法的搜尋範圍但易於陷入區域性最優，較窄的取值範圍會提高演算法的精度但易於陷入搜尋空間的侷限性。
4.4 示例
# 匯入必要的模組
import numpy as np
import sko
# 定義 Branin's rcos 函式
# 輸入引數為 p = [x, y]，其中 -5 <= x <= 10, 0 <= y <= 15
# 函式的三個最小值為 (x, y) = (-pi, 12.275), (pi, 2.275), (9.42478, 2.475)，對應的函式值為 0.397887
# 函式還有一個區域性最小值為 0.8447
def Brains_rcos(p):
x, y = p
z = (y - 51*(x**2)/(40*np.pi**2) + 5*x/np.pi - 6)**2 + 10*(1-1/(8*np.pi))*np.cos(x) + 10
return z
# 使用 Differential Evolution（DE）演算法求解 Branin's rcos 函式
# 設定搜尋空間為 2 維，種群大小為 50，最大迭代次數為 800
# 搜尋空間的下界為 [-5, 0]，上界為 [10, 15]
ex_DE = sko.DE.DE(func=Brains_rcos, n_dim=2, size_pop=50, max_iter=800, lb=[-5, 0,], ub=[10, 15,])
# 執行 Differential Evolution 演算法，獲取最優解及其對應的函式值
best_x, best_y = ex_DE.run()
# 列印最優解的取值和相應的函式值
print('best_x:', best_x, '\n', 'best_y:', best_y)
best_x: [9.42477795 2.47499997]
best_y: [0.39788736]
特別篇：GA與DE的區別
差分進化演算法（Differential Evolution，DE）和遺傳演算法（Genetic Algorithm，GA）都是最佳化問題的常用解法，它們本質上都是一種基於種群的進化演算法。但是從演算法的具體實現角度來看，它們之間存在一些區別，包括：
- 變異操作的方式：DE是透過差分操作對整個種群進行變異，而GA是透過單點交叉、多點交叉、變異等操作對個體進行變異。
- 交叉操作的方式：DE不使用交叉操作，而GA是透過交叉操作將兩個個體的染色體交換基因片段。
- 引數的個數：DE只有 3 個引數（個體數、變異機率、縮放因子），而GA有更多的引數，如交叉機率、變異機率、選擇演算法、染色體長度等。
- 演算法的收斂速度：DE透過區域性搜尋和全域性搜尋相結合的方式，通常能夠在較少的迭代次數內得到較為理想的全域性最優解；而GA通常較難避免陷入區域性最優解，收斂速度可能會比DE慢。
總的來說，DE更適合處理具有連續性的目標函式，對噪聲和非線性函式的適應性非常強，而GA更適合解決具有離散性或組合性質的最佳化問題，如旅行商問題等。
5.SA(Simulated Annealing)
5.1 解決的問題
- 目標函式最佳化
- TSP(Travelling Salesman Problem)問題
5.2 API
5.2.1 Simulated Annealing
class SAFast(SimulatedAnnealingValue):
def __init__(self, func, x0, T_max=100, T_min=1e-7, L=300, max_stay_counter=150, **kwargs):...
| 引數 | 含義 |
|---|---|
| func | 欲最佳化的目標函式 |
| $x_{0}$ | 目標函式的初始點 |
| T_max | 初始溫度 |
| T_min | 最小溫度 |
| L | 內迴圈迭代次數 |
| max_stay_counter | 最大停留次數 |
| lb(lower bound) | 函式引數的下限 |
| ub(upper bound) | 函式引數的上限 |
引數詳解：
- func：最佳化器所要最佳化的目標函式。它是一個接受輸入引數x並返回實數值的函式。
- x0：最佳化器的起始點，它是一個長度為n的列表，n為目標函式的自變數數量。
- T_max：初始溫度的上限。在退火過程中，溫度將不斷地減小，直到小於等於T_min。T_max的大小決定了退火過程中的起始溫度，取值較高可以使得更多的移不動的狀態被接受，但同時可能導致計算時間增加和結果不穩定。
- T_min：溫度的下限。當溫度降到T_min以下時，退火過程結束。T_min的大小影響退火過程的終止時間以及最終結果的質量，通常情況下應該設定為一個很小的正數，並與目標函式的值域有關。
- L：每個溫度下的迭代次數。每次迭代中，最佳化器從當前點x嘗試隨機擾動得到一個新的點x_new，如果目標函式值變小，則接受x_new作為新的點；否則以一定機率接受該點以避免陷入區域性最小值。L的大小決定每個溫度下的計算量，取值過小可能導致結果不充分，取值過大可能導致計算時間增加。
- max_stay_counter：接受狀態的最大連續次數。如果連續生成的L個新狀態都沒有比當前狀態更優，則認為演算法已經陷入區域性極小值，此時退出迭代，返回當前狀態作為最優解。max_stay_counter的大小影響演算法的全域性搜尋能力，設定過小可能導致演算法過早退出，在區域性極小值附近震盪，設定過大可能導致演算法無法收斂於正確解。
調整這些引數會對模擬退火演算法的效能和結果產生影響，一般地：
- 如果發現演算法很難找到可接受的狀態，可適當提高T_max並減小L，以增加接受狀態的可能性，並減少在狀態空間中鎖定的時間。
- 如果想縮短演算法時間或者降低計算量、並行化搜尋過程，可以減小L、max_stay_counter或T_max等引數，以減少每個溫度下的迭代次數、搜尋到的狀態數或搜尋的狀態空間範圍。
- 如果想得到更精確的結果，可以增加L、max_stay_counter或T_min等引數，以擴大搜尋範圍或提高演算法在搜尋到區域性最優解後跳出區域性最優解的機率。
5.2.2 SA for TSP
class SimulatedAnnealingBase(SkoBase):
def __init__(self, func, x0, T_max=100, T_min=1e-7, L=300, max_stay_counter=150, **kwargs):...
class SA_TSP(SimulatedAnnealingBase):...
5.3 示例
5.3.1 目標函式最佳化
# 匯入numpy庫，函式部分實現依賴numpy
import numpy as np
# 定義Griewank函式
def Griewank(x):
'''
Griewank函式, 在x=(0, 0, ..., 0)處有全域性極小值0,
-600<=x<=600, x為向量，維數可任意大小
'''
# 計算第一項
y1 = np.sum(x**2 / 4000)
y2 = 1
# 計算第二項
for i in range(1, x.size+1):
y2=y2 * np.cos(x[i-1] / np.sqrt(i))
# 返回Griewank函式值
return y1-y2+1
# 匯入sko庫，這裡是使用SA演算法
import sko
# 初始化SA模擬退火演算法類，併為其提供Griewank函式作為最佳化目標，以及初始值x0、其他配置引數（最大溫度T_max、最小溫度T_min、鏈長L、最大停留次數max_stay_counter、函式引數下界lb、函式引數上界ub）
ex_sa = sko.SA.SA(func=Griewank, x0=[1, 1, 1], T_max=1, T_min=1e-9, L=300, max_stay_counter=150
,lb=[-600,-600,-600,],ub=[600,600,600,])
# 執行模擬退火演算法，並獲取最優的x值和對應的函式值y
best_x, best_y = ex_sa.run()
# 輸出當前模擬得到的最優的x值和最優的函式值y
print('best_x:', best_x, 'best_y', best_y)
# 匯入matplotlib.pyplot庫和pandas庫，用於繪製資料視覺化效果
import matplotlib.pyplot as plt
import pandas as pd
# 繪製最優結果歷史的累計最小值（畫圖並儲存）。累計最小值表示在每個時間點之前，所有的資料中的最小值（而不是在特定時間點的最小值）。
# 透過cummin方法計算，在“axis=0”的方向上進行累計最小值計算。這裡，axis=0 表示跨行（沿著第0軸）進行計算
plt.plot(pd.DataFrame(ex_sa.best_y_history).cummin(axis=0))
# 儲存並顯示影像
plt.savefig("SA.png",dpi=400)
plt.show()
5.3.2 TSP
# 匯入必要的庫
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
# 隨機生成點的數目
num_points = 10
# 生成隨機座標陣列
points_coordinate = np.random.rand(num_points, 2)
# 生成距離矩陣，即每對點之間的歐氏距離
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
# 定義目標函式
def cal_total_distance(routine):
'''
計算旅行商問題的總距離，輸入路線routine，返回總距離
範例：cal_total_distance(np.arange(num_points))
'''
# 獲取路線上的點的數目
num_points, = routine.shape
# 計算每對相鄰點之間的距離之和
return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
# 匯入SA_TSP演算法函式
from sko.SA import SA_TSP
# 初始化SA_TSP演算法類，併為其提供生成的目標函式cal_total_distance、初始路線x0、溫度引數T_max, T_min和其他配置引數（鏈長L）
sa_tsp = SA_TSP(func=cal_total_distance, x0=range(num_points), T_max=100, T_min=1, L=10 * num_points)
# 執行模擬退火演算法，並獲取最優的路線best_points及對應的距離best_distance
best_points, best_distance = sa_tsp.run()
# 輸出執行的結果
print(best_points, best_distance, cal_total_distance(best_points))
# 以下程式碼用於繪製最優路線的視覺化效果
# 匯入FormatStrFormatter庫
from matplotlib.ticker import FormatStrFormatter
# 建立兩張畫布
fig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
# 繪製距離迭代歷史的折線圖
ax[0].plot(sa_tsp.best_y_history)
ax[0].set_xlabel("Iteration")
ax[0].set_ylabel("Distance")
# 繪製最優路線的散點圖
ax[1].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1],
marker='o', markerfacecolor='b', color='c', linestyle='-')
ax[1].xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].yaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].set_xlabel("Longitude")
ax[1].set_ylabel("Latitude")
# 儲存並顯示影像
plt.savefig("SA_TSP.png",dpi=400)
plt.show()
5.4 bug
- 本人發現在使用SA_TSP的途中，只要問題的複雜度上升，SA_TSP幾乎必定會給出錯誤的結果。在上面的例子中，我只給出了10個點位，這時的SA_TSP看起來適應地很好。但是，當我把點位增加到50，SA_TSP給出瞭如下的錯誤結果：
- 如有讀者知道原因，請聯絡本人！
6.ACA (Ant Colony Algorithm) for tsp
6.1 解決的問題
- TSP(Travelling Salesman Problem)問題
6.2 API
class ACA_TSP:
def __init__(self, func, n_dim,
size_pop=10, max_iter=20,
distance_matrix=None,
alpha=1, beta=2, rho=0.1,
):...
6.3 引數
| 引數 | 含義 |
|---|---|
| func | 欲最佳化的目標函式 |
| n_dim | 所給函式的引數個數 |
| size_pop | 種群初始規模 |
| max_iter | 最大迭代次數 |
| distance_matrix | TSP問題的距離矩陣 |
| alpha | 資訊素重要程度因子 |
| beta | 資訊素濃度因子 |
| rho | 資訊素揮發係數 |
size_pop：種群規模對演算法的效率和精度有很大的影響。種群規模越大，演算法在搜尋空間中的覆蓋率越高，但相應的計算開銷也越大，執行時間和記憶體開銷也會增加。如果種群規模過小，演算法可能會陷入區域性最優解，從而影響演算法的質量。max_iter：迭代次數越多，演算法對解的搜尋空間越廣，容易得到更優的解，但是也會消耗更多的計算資源，從而影響演算法的執行效率。distance_matrix：距離矩陣的準確性對演算法的精度和執行效率有很大影響。如果距離矩陣不準確，演算法可能會得到錯誤的結果，從而影響演算法的質量。距離矩陣的計算也需要很大的計算開銷，如果距離矩陣過於複雜，演算法的執行效率也會受到影響。alpha和beta：資訊素重要程度因子和啟發函式重要程度因子都可以影響演算法的全域性搜尋能力和區域性搜尋能力。資訊素重要程度因子越大，螞蟻走過的路徑上資訊素的影響越大，演算法的全域性搜尋能力越強，但是區域性搜尋能力會弱化。啟發函式重要程度因子越大，螞蟻選擇路徑時更加傾向於走近目標城市，因此演算法的區域性搜尋能力會增強，但是全域性搜尋能力會弱化。rho：資訊素揮發程度因子主要用於控制資訊素的更新速度。如果揮發程度因子設定得過大，資訊素的更新速度會很快，導致演算法的不穩定性增加；如果設定得過小，資訊素的更新速度會很慢，影響演算法的收斂速度。
需要根據具體問題的特點和問題需求來選擇和調整引數。
6.4 示例
# 匯入常用庫
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
import pandas as pd
# 設定待確定節點的數量為50
num_points = 50
# 隨機生成50個在[0,1)之間的點
# points_coordinate[i][0]、points_coordinate[i][1]是第i個點的x,y座標
points_coordinate = np.random.rand(num_points, 2)
# 採用歐幾里得距離度量計算每兩個點之間的距離生成一個距離矩陣
# distance_matrix[i][j]是第i個點與第j個點之間的距離
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
# 目標函式: 輸入一個排列向量，返回TSP問題中按該排列順序訪問所有點後的總距離
def cal_total_distance(routine):
'''The objective function. input routine, return total distance.
cal_total_distance(np.arange(num_points))
'''
# 計算訪問所有點所需的總距離
num_points, = routine.shape
return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
# 匯入ACA_TSP蟻群演算法，將目標函式、節點數量、種群數量、迭代次數、距離矩陣等引數傳入該演算法物件中
from sko.ACA import ACA_TSP
aca = ACA_TSP(func=cal_total_distance, n_dim=num_points,
size_pop=50, max_iter=200,
distance_matrix=distance_matrix)
# 執行蟻群演算法，尋找最優解
best_x, best_y = aca.run()
# 輸出結果
print(best_x,best_y)
# 畫出路徑和迭代過程中最優解的變化，儲存圖片並顯示
fig, ax = plt.subplots(1, 2)
#合併最優路徑
best_points_ = np.concatenate([best_x, [best_x[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
# 畫出最優路徑
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
# 畫出迭代過程中最優解的變化
pd.DataFrame(aca.y_best_history).cummin().plot(ax=ax[1])
plt.savefig("ACA_TSP.png",dpi=400)
plt.show()
7.immune algorithm (IA)
7.1 解決的問題
- TSP(Travelling Salesman Problem)問題
7.2 API
class IA_TSP(GA_TSP):
def __init__(self, func, n_dim, size_pop=50, max_iter=200, prob_mut=0.001, T=0.7, alpha=0.95):...
7.3 引數
| 引數 | 含義 |
|---|---|
| func | TSP問題的目標函式 |
| n_dim | 所給函式的引數個數 |
| size_pop | 種群初始規模 |
| max_iter | 最大迭代次數 |
| prob_mut | 個體突變的機率 |
| T | 模擬退火演算法的初始溫度 |
| alpha | 模擬退火演算法中溫度下降的比例 |
- size_pop：增加種群數量可能會導致更多的個體能夠探索更廣闊的解空間，但是也會增加演算法的計算時間和消耗的記憶體。
- max_iter：增加迭代次數可能會導致更好的解或計算時間更長。
- prob_mut：增加突變機率能夠使演算法更有機會在探索過程中跳出區域性最優解，但是也會增加計算複雜度。
- T：增加初始溫度可能會使演算法能夠更快地探索解空間，但是也可能導致演算法探索不夠全面或收斂不夠徹底。
- alpha：增加溫度下降比例可以加速演算法的收斂速度，但是也可能導致演算法過早收斂到次優解。
7.4 示例
# 匯入常用庫
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
# 設定待確定節點的數量為50
num_points = 50
# 隨機生成50個在[0,1)之間的點
# points_coordinate[i][0]、points_coordinate[i][1]是第i個點的x,y座標
points_coordinate = np.random.rand(num_points, 2)
# 採用歐幾里得距離度量計算每兩個點之間的距離生成一個距離矩陣
# distance_matrix[i][j]是第i個點與第j個點之間的距離
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
# 目標函式: 輸入一個排列向量，返回TSP問題中按該排列順序訪問所有點後的總距離
def cal_total_distance(routine):
'''The objective function. input routine, return total distance.
cal_total_distance(np.arange(num_points))
'''
# 計算訪問所有點所需的總距離
num_points, = routine.shape
return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
# 匯入IA_TSP蟻群演算法，將目標函式、節點數量、種群數量、迭代次數、交叉機率、變異機率等引數傳入該演算法物件中
from sko.IA import IA_TSP
ia_tsp = IA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=500, max_iter=800, prob_mut=0.2,
T=0.7, alpha=0.95)
# 執行蟻群演算法，尋找最優解
best_points, best_distance = ia_tsp.run()
# 輸出結果
print('best routine:', best_points, 'best_distance:', best_distance)
# 畫出路徑和迭代過程中最優解的變化，儲存圖片並顯示
fig, ax = plt.subplots(1, 2)
# 合併最優路徑
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
# 畫出最優路徑
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
# 畫出迭代過程中最優解的變化
ax[1].plot(ia_tsp.generation_best_Y)
# 儲存圖片並顯示
plt.savefig("IA_TSP.png",dpi=400)
plt.show()
8.Artificial Fish Swarm Algorithm (AFSA)
8.1 解決的問題
- 目標函式最佳化
8.2 API
class AFSA:
def __init__(self, func, n_dim, size_pop=50, max_iter=300,
max_try_num=100, step=0.5, visual=0.3,
q=0.98, delta=0.5):...
8.3 引數
| 引數 | 含義 |
|---|---|
| func | 欲最佳化的目標函式 |
| n_dim | 所給函式的引數個數 |
| size_pop | 種群初始規模 |
| max_iter | 最大迭代次數 |
| max_try_num | 當某一代的最優解無法更新時，最多重複嘗試的次數 |
| step | 每步的步長 |
| visual | 視覺化介面的比例 |
| q | 溫度調節的引數，用於控制降溫速度 |
| delta | 更新解的相關引數，控制上一次最優解和當前最優解的距離 |
這些引數的調整會影響AFSA演算法的收斂性、全域性最優解的尋優效果以及時間成本。具體來說：
- size_pop：一般來說，種群數量越大，能夠全域性搜尋的能力就越強，但同時也會增加計算時間和記憶體消耗。如果時間和記憶體資源允許，可以考慮適當增加該引數的值。
- max_iter：迭代次數過低可能導致演算法收斂不到最優解，而迭代次數過高可能會浪費計算時間，同時也可能陷入區域性最優解而無法跳出。需要根據實際問題進行調節。
- max_try_num：該引數的增加會增加演算法迭代所需的時間，但同時可以增加演算法跳出區域性最優解的可能性。如果演算法收斂不到最優解，可以適當增加該引數的值。
- step：如果步長設定過大，演算法可能會在搜尋過程中快速跳過最優解附近的板塊，從而無法找到最優解。如果步長過小，演算法搜尋效率較低。需要根據具體問題進行調節。
- visual：視覺化介面可以幫助理解演算法搜尋過程以及尋優效果，但同時也會佔用計算資源。可以根據具體需求進行調節。
- q：較低的q值會減緩溫度（降溫）的下降速度，可能會導致演算法長時間在區域性最優解周圍徘徊。較高的q值將降溫速度加快，可能有助於演算法跳出區域性最優解，但同時也可能會降低演算法全域性最佳化的精度。需要根據實際問題進行調節。
- delta：delta引數越大，跳出區域性最優解的可能性越大，但同時可能會增加演算法的計算時間和空間複雜度。如果演算法很難跳出區域性最優解，可以適當增加該引數的值。
8.4 示例
# 匯入常用庫
import numpy as np
# 目標函式，輸入一個二元向量(p), 返回Easom's函式在該點的函式值
def Easom_s(p):
'''
Easom's函式,-100<=x<=100, -100<=y<=100
有一個全域性最小值-1在座標(pi,pi)
'''
x,y = p
z = -np.cos(x) * np.cos(y) * np.exp(-((x-np.pi)**2 + (y-np.pi)**2))
return z
# 匯入AFSA最佳化演算法，將目標函式、變數維度、種群大小、最大迭代次數、最大嘗試次數、步長、視覺化引數、
# 溫度下降速率、步長下降速率等引數傳入該演算法物件中
from sko.AFSA import AFSA
afsa = AFSA(func=Easom_s, n_dim=2, size_pop=50, max_iter=300,
max_try_num=100, step=5, visual=0.3,
q=0.98, delta=1)
# 執行演算法計算最優解
best_x, best_y = afsa.run()
# 輸出結果
print(best_x, best_y)
[3.14173594 3.14302929] -0.9999968733295336
- 如果把步長縮小成0.5，delta縮小成0.5，那麼這個演算法針對Easom's函式給出了錯誤的答案
[1.30500828 1.30513779] -8.110222525091613e-05
本文由部落格一文多發平臺 OpenWrite 釋出！