Firefly Algorithm(萤火虫算法,FLA)是一种启发式优化算法,其灵感来源于萤火虫的闪烁行为。FLA算法通过模拟萤火虫群体中的个体之间的相互吸引和相对亮度来搜索解空间,主要应用于全局优化问题,例如路径规划、函数优化等。
Firefly Algorithm(萤火虫算法,FLA)是一种启发式优化算法,其灵感来源于萤火虫的闪烁行为。FLA算法通过模拟萤火虫群体中的个体之间的相互吸引和相对亮度来搜索解空间,主要应用于全局优化问题,例如路径规划、函数优化等。
实现 FLA的基本步骤如下所示。
(1)初始化萤火虫群体:在解空间中随机生成一定数量的萤火虫,每个萤火虫代表一个可能的解。初始亮度可以通过问题的适应性函数计算得到。
(2)定义亮度和吸引力:萤火虫的亮度表示其解的优劣,亮度值越高表示解越好。吸引力表示萤火虫之间的相互影响,通常与亮度和距离有关。
(3)更新萤火虫位置:每个萤火虫根据其亮度和与其他萤火虫的距离来更新其位置。亮度越大的萤火虫具有更强的吸引力,其他萤火虫会朝着更亮的方向移动。
(4)调整亮度和位置:根据问题的适应性函数,重新计算每个萤火虫的亮度。然后,根据吸引力和亮度,调整每个萤火虫的位置。
(5)重复迭代过程:重复执行步骤3和步骤4,直到满足停止条件,例如达到最大迭代次数或适应性函数值足够小。
FLA被广泛应用于各种优化问题,包括路径规划、函数优化、组合优化等。它是一种适应性强、易于实现的优化算法,对于一些复杂的问题能够取得良好的效果。请看下面的实例,演示了使用FLA优化无人机路径规划的过程。本项目使用FLA实现无人机路径规划的优化操作,通过最小化路径长度和能量消耗来优化无人机的轨迹,同时避免环境中的障碍物。FLA算法通过调整多维搜索空间中分子的位置来寻找最优解。为了评估该方法的性能,研究人员在模拟环境中进行了实验,考虑了各种障碍物配置。
实例4-10:使用FLA优化无人机的路径规划(源码路径:codes\4\gui\fla-algorithm.ipynb)
在本项目中,使用Fick's Law Algorithm (FLA)进行无人机路径规划和优化。通过使用FLA算法,本项目能够在考虑路径长度、能耗和避免环境障碍的情况下自动搜索最优路径,实现了路径规划模型、随机解生成、成本函数定义以及FLA算法的主要步骤。通过调用FLA函数,得到最优路径解,并提供了绘制二维和三维路径规划图的功能,使用户能够直观地了解无人机在复杂环境中的优化路径规划效果。项目整体为无人机路径规划提供了一个基于FLA的优化方法,并通过实验验证了其性能。
实例文件fla-algorithm.ipynb的具体实现流程如下所示。
(1)导入需要的库
导入了用于处理数据和进行可视化的必要库,包括NumPy用于线性代数操作,Pandas用于数据处理和CSV文件的输入输出,Matplotlib用于绘图,Scipy中的io模块用于保存.mat文件,以及一些警告处理。其中 %matplotlib inline 是 Jupyter Notebook 中的魔法命令,用于在输出单元格中显示图形。
(2)建模
定义函数 create_model,用于创建一个模型,包括起点 (xs, ys, zs)、终点 (xt, yt, zt)、障碍物的位置 (xobs, yobs, zobs) 以及它们的半径 (robs)。同时,在定义的模型中还包含了路径网格的数量 (n) 和空间边界的范围 (xmin, xmax, ymin, ymax, zmin, zmax),这个模型将用于后续的路径规划和优化工作。
def create_model():
# Source
xs = 0
ys = 0
zs = 0
# Target (Destination)
xt = 4
yt = 5
zt = 5
xobs = [1.5, 4.0, 1.2]
yobs = [4.5, 3.0, 1.5]
zobs = [4.5, 3.0, 1.5]
robs = [0.8, 0.6, 0.8]
n = 10
xmin = -10
xmax = 10
ymin = -10
ymax = 10
zmin = -10
zmax = 10
model = {
"xs": xs,
"ys": ys,
"zs": zs,
"xt": xt,
"yt": yt,
"zt": zt,
"xobs": xobs,
"yobs": yobs,
"zobs": zobs,
"robs": robs,
"n": n,
"xmin": xmin,
"xmax": xmax,
"ymin": ymin,
"ymax": ymax,
"zmin": zmin,
"zmax": zmax,
}
return model
(3)创建随机解决方案
函数 create_random_solution的功能是根据给定的模型生成一个随机的解决方案。其中,解决方案包括在模型定义的路径网格上随机生成的一组 x 和 y 坐标。这将用作路径规划的初始解决方案。
def create_random_solution(model):
n = model["n"]
xmin = model["xmin"]
xmax = model["xmax"]
ymin = model["ymin"]
ymax = model["ymax"]
sol1 = {
"x": np.random.uniform(xmin, xmax, size=n),
"y": np.random.uniform(ymin, ymax, size=n),
}
return sol1
(4)定义成本函数
函数 my_cost的功能是计算给定解决方案 sol1 的成本。在函数内部,首先通过调用 create_model 函数创建了模型,然后通过调用 parse_solution 函数将输入的分子位置解析成路径规划模型的具体形式。最后,通过计算路径长度 sol["L"] 以及一个惩罚项 beta * sol["Violation"] 的和来得到成本 z。其中,beta 是一个调整参数。在路径优化中,将通过这个成本函数评估解决方案的质量。
def my_cost(sol1):
model = create_model()
new_Molecules = {'x': np.sort(sol1[0:10]), 'y': np.sort(sol1[10:20])}
sol = parse_solution(new_Molecules, model)
beta = 100
z = sol["L"] * (1 + beta * sol["Violation"])
return z
(5)解析解决方案
定义函数parse_solution,将输入的分子位置解析成路径规划模型的具体形式。在函数内部,通过使用样条插值 splrep 和 splev 对给定的 x 和 y 坐标进行插值,得到路径的平滑曲线。然后,计算路径的长度 L,并检查路径是否避开了障碍物,通过计算路径与障碍物的距离,生成了违规项 Violation。
from scipy.interpolate import splprep, splev, splrep, interp1d
def parse_solution(sol1, model):
x = sol1["x"]
y = sol1["y"]
xs = model["xs"]
ys = model["ys"]
xt = model["xt"]
yt = model["yt"]
xobs = model["xobs"]
yobs = model["yobs"]
robs = model["robs"]
XS = [xs] + x.tolist() + [xt]
YS = [ys] + y.tolist() + [yt]
k = len(XS)
TS = np.linspace(0, 1, k)
tt = np.linspace(0, 1, 100)
xx = splrep(TS, XS, k=3)
yy = splrep(TS, YS, k=3)
x_interpolated = splev(tt, xx)
y_interpolated = splev(tt, yy)
dx = np.diff(x_interpolated)
dy = np.diff(y_interpolated)
L = np.sqrt(np.sum(np.square(dx) + np.square(dy)))
nobs = len(xobs)
Violation = 0
for k in range(nobs):
d = np.sqrt((x_interpolated - xobs[k])**2 + (y_interpolated - yobs[k])**2)
v = np.maximum(1 - d/(robs[k]+0.3), 0)
Violation += np.mean(v)
sol2 = {
"TS": TS,
"XS": XS,
"YS": YS,
"tt": np.linspace(0, 1, 100),
"xx": x_interpolated,
"yy": y_interpolated,
"dx": dx,
"dy": dy,
"L": L,
"Violation": Violation,
"IsFeasible": (Violation == 0),
}
return sol2
在上面的代码中,函数 parse_solution将返回一个包含路径的各种信息的字典 sol2,包括插值后的路径坐标 (x_interpolated, y_interpolated)、路径长度 (L)、是否违规 (IsFeasible) 等。这个信息将在路径优化过程中用于评估解决方案的质量。
(6)绘制解决方案
- 函数 PlotSolution用于绘制路径规划的解决方案。在函数内部,首先调用 create_model 函数创建模型,然后调用 parse_solution 函数将输入的分子位置解析成具体的路径信息。接着,通过 Matplotlib 绘制了包括障碍物、起点、终点以及路径在内的可视化图形。
def PlotSolution(bst):
model = create_model()
bst_Molecules = {'x': np.sort(bst[0:10]), 'y': np.sort(bst[10:20])}
sol = parse_solution(bst_Molecules, model)
xs = model['xs']
ys = model['ys']
xt = model['xt']
yt = model['yt']
xobs = model['xobs']
yobs = model['yobs']
robs = model['robs']
XS = sol['XS']
YS = sol['YS']
xx = sol['xx']
yy = sol['yy']
theta = np.linspace(0, 2*np.pi, 100)
for k in range(len(xobs)):
plt.fill(xobs[k] + robs[k]*np.cos(theta), yobs[k] + robs[k]*np.sin(theta), [0.5, 0.7, 0.8])
plt.plot(xx, yy, 'k', linewidth=2)
plt.plot(XS, YS, 'ro')
plt.plot(xs, ys, 'bs', markersize=12, markerfacecolor='y')
plt.plot(xt, yt, 'kp', markersize=16, markerfacecolor='g')
plt.grid(True)
plt.axis('equal')
plt.show()
具体来说,函数 PlotSolution绘制了障碍物的轮廓、规划的路径、起点、终点,并使用不同的颜色和形状进行标识。最终,通过调用 plt.show() 显示绘制的图形。这个函数可用于可视化路径规划的结果。
- 函数 PlotSolution3D 用于绘制三维空间中的路径规划解决方案。在函数内部,通过调用 create_model 函数创建模型,然后通过调用 parse_solution 函数将输入的分子位置解析成具体的路径信息。接着,使用 Matplotlib 的 3D 投影 projection='3d' 创建了一个三维坐标轴。
# 绘制三维解决方案
def PlotSolution3D(bst):
# 创建模型
model = create_model()
# 提取最佳分子位置
bst_Molecules = {'x': np.sort(bst[0:10]), 'y': np.sort(bst[10:20])}
# 解析最佳解决方案
sol = parse_solution(bst_Molecules, model)
# 提取模型参数
xs = model['xs']
ys = model['ys']
xt = model['xt']
yt = model['yt']
xobs = model['xobs']
yobs = model['yobs']
robs = model['robs']
# 提取解析后的路径信息
XS = sol['XS']
YS = sol['YS']
xx = sol['xx']
yy = sol['yy']
# 创建三维坐标轴
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# 绘制障碍物
theta = np.linspace(0, 2*np.pi, 100)
for k in range(len(xobs)):
# 创建立方体障碍物的顶点
vertices = np.array([[xobs[k]-robs[k], yobs[k]-robs[k], 0],
[xobs[k]-robs[k], yobs[k]+robs[k], 0],
[xobs[k]+robs[k], yobs[k]+robs[k], 0],
[xobs[k]+robs[k], yobs[k]-robs[k], 0],
[xobs[k]-robs[k], yobs[k]-robs[k], robs[k]*2],
[xobs[k]-robs[k], yobs[k]+robs[k], robs[k]*2],
[xobs[k]+robs[k], yobs[k]+robs[k], robs[k]*2],
[xobs[k]+robs[k], yobs[k]-robs[k], robs[k]*2]])
# 创建立方体障碍物的面
faces = np.array([[0, 1, 2, 3],
[4, 5, 6, 7],
[0, 1, 5, 4],
[1, 2, 6, 5],
[2, 3, 7, 6],
[3, 0, 4, 7]])
# 计算障碍物的高度范围
zmin = np.min(vertices[:, 2])
zmax = np.max(vertices[:, 2])
# 创建 Poly3DCollection 对象并添加到绘图中
obstacle = Poly3DCollection(vertices[faces], alpha=0.4, cmap='plasma')
obstacle.set_array((vertices[:, 2] - zmin) / (zmax - zmin))
obstacle.set_facecolor([0.5, 0.7, 0.8])
ax.add_collection3d(obstacle)
# 绘制规划路径
ax.plot(xx, yy, np.zeros_like(xx), 'k', linewidth=2)
# 标记起点和终点
ax.plot(XS, YS, np.zeros_like(XS), 'ro')
ax.plot(xs, ys, np.zeros_like(xs), 'bs', markersize=12, markerfacecolor='y')
ax.plot(xt, yt, np.zeros_like(xt), 'kp', markersize=16, markerfacecolor='g')
# 设置坐标轴范围
ax.grid(True)
ax.set_xlim3d(np.min(xobs) - np.max(robs), np.max(xobs) + np.max(robs))
ax.set_ylim3d(np.min(yobs) - np.max(robs), np.max(yobs) + np.max(robs))
ax.set_zlim3d(0, np.max(robs) * 2)
在上述代码中,函数 PlotSolution3D绘制了三维空间中的立方体障碍物、规划的路径、起点、终点,并使用不同的颜色和形状进行标识。最终,通过设置坐标轴的范围,使其适应障碍物和路径的范围,并通过 ax.set_xlim3d、ax.set_ylim3d、ax.set_zlim3d 来设置坐标轴的范围。这个函数可用于可视化三维路径规划的结果。
(7)定义FLA优化函数
定义函数FLA,实现了一种基于领域分组和自适应温度的自组织火焰优化算法(FLA),用于求解优化问题。该算法通过模拟分子在温度场中的扩散与聚集行为,采用方向操作符(DO)、方向流动操作符(TDO)和稳态操作符(SSO)来更新分子群体的位置,以寻找全局最优解。算法采用自适应降温策略和分组策略,最终输出优化结果及其对应的最小成本。
def FLA( NoMolecules, T, lb, ub, dim, objfunc):
C1 = 0.5
C2 = 2
c3 = 0.1
c4 = 0.2
c5 = 2
D = 0.01
eps = 2.2204*(10**-16)
X = lb + np.random.rand(NoMolecules, dim) * (ub - lb) # intial positions
FS = np.zeros(NoMolecules)
for i in range(NoMolecules):
FS[i] = CostFunction(X[i])
BestF = np.min(FS)
IndexBestF = np.argmin(FS)
Xss = X[IndexBestF,:]
n1 = round(NoMolecules/2)
n2 = NoMolecules - n1
X1 = X[0:n1,:]
X2 = X[n1:NoMolecules,:]
FS1 = np.zeros(n1)
FS2 = np.zeros(n2)
for i in range(n1):
FS1[i] = CostFunction(X1[i,:])
for i in range(n2):
FS2[i] = CostFunction(X2[i,:])
FSeo1, IndexFSeo1 = np.min(FS1), np.argmin(FS1)
FSeo2, IndexFSeo2 = np.min(FS2), np.argmin(FS2)
Xeo1, Xeo2 = X1[IndexFSeo1,:], X2[IndexFSeo2,:]
vec_flag = [1, -1]
if FSeo1 < FSeo2:
FSss, YSol = FSeo1, Xeo1
else:
FSss, YSol = FSeo2, Xeo2
T = Max_iteration
X1new = np.zeros((30,dim))
X2new = np.zeros((30,dim))
CNVG = []
# Main Loop
for t in range(Max_iteration):
TF = np.sinh(t / T) ** C1
X = np.concatenate((X1, X2), axis=0)
# DO
if TF < 0.9:
DOF = np.exp(-(C2*TF - np.random.rand()))**C2
# direction of flow
TDO = c5*TF - np.random.rand()
if TDO < np.random.rand():
# select no of molecules
M1N = c3*n1;
M2N = c4*n1;
NT12 = np.round((M2N-M1N) * np.random.rand(1, 1) + M1N)[0,0]
for u in range(1, int(NT12)):
flag_index = np.floor(2 * np.random.rand() + 1)
DFg = vec_flag[int(flag_index)-1]
Xm2 = np.mean(X2)
Xm1 = np.mean(X1)
J = -D * (Xm2 - Xm1) / np.linalg.norm(Xeo2 - X1[u, :] + np.finfo(float).eps)
X1new[u,:] = Xeo2 + DFg * DOF * np.random.rand(1, dim) * (J * Xeo2 - X1[u, :])
for u in range(int(NT12)+1, n1):
for tt in range(dim):
p = np.random.rand()
if p < 0.8:
X1new[u,tt] = Xeo1[tt]
elif p < 0.9:
r3 = np.random.rand()
X1new[u,tt] = X1[u,tt] + DOF * ((ub-lb) * r3 + lb)
else:
X1new[u,tt] = X1[u,tt]
for u in range(n2):
r4 = np.random.rand()
X2new[u,:] = Xeo2 + DOF * ((ub-lb) * r4 + lb)
else:
M1N = 0.1 * n2
M2N = 0.2 * n2
Ntransfer = int(round((M2N - M1N) * float(np.random.rand()) + M1N))
for u in range(1, Ntransfer+1):
flag_index = np.floor(2 * np.random.rand() + 1)
DFg = vec_flag[int(flag_index)-1]
R1 = np.random.randint(n1)
Xm1 = np.mean(X1, axis=0)
Xm2 = np.mean(X2, axis=0)
J = -D * (Xm1 - Xm2) / np.linalg.norm(Xeo1 - X2[u-1, :] + np.finfo(float).eps)
X2new[u-1, :] = Xeo1 + DFg * DOF * np.random.rand(1, dim) * (J * Xeo1 - X2[u-1, :])
for u in range(Ntransfer, n2):
for tt in range(dim):
p = np.random.rand()
if p < 0.8:
X2new[u, tt] = Xeo2[tt]
elif p < 0.9:
r3 = np.random.rand()
X2new[u, tt] = X2[u, tt] + DOF * ((ub - lb) * r3 + lb)
else:
X2new[u, tt] = X2[u, tt]
for u in range(n1):
r4 = np.random.rand()
X1new[u, :] = Xeo1 + DOF * ((ub - lb) * r4 + lb)
else:
if TF <= 1:
for u in range(n1):
flag_index = np.floor(2 * np.random.rand() + 1)
DFg = vec_flag[int(flag_index)-1]
Xm1 = np.mean(X1, axis=0)
Xmeo1 = Xeo1
J = -D * (Xmeo1 - Xm1) / (np.linalg.norm(Xeo1 - X1[u, :] + np.finfo(float).eps))
DRF = np.exp(-J / TF)
MS = np.exp(-FSeo1 / (FS1[u] + np.finfo(float).eps))
R1 = np.random.rand(dim)
Qeo = DFg * DRF * R1
X1new[u, :] = Xeo1 + Qeo * X1[u, :] + Qeo * (MS * Xeo1 - X1[u, :])
for u in range(n2):
flag_index = np.floor(2 * np.random.rand() + 1)
DFg = vec_flag[int(flag_index)-1]
Xm2 = np.mean(X2, axis=0)
Xmeo2 = Xeo2
J = -D * (Xmeo2 - Xm2) / (np.linalg.norm(Xeo2 - X2[u, :] + np.finfo(float).eps))
DRF = np.exp(-J / TF)
MS = np.exp(-FSeo2 / (FS2[u] + np.finfo(float).eps))
R1 = np.random.rand(dim)
Qeo = DFg * DRF * R1
X2new[u, :] = Xeo2 + Qeo * X2[u, :] + Qeo * (MS * Xeo1 - X2[u, :])
else:
# Steady state operator (SSO)
for u in range(n1):
flag_index = int(np.floor(2 * np.random.rand() + 1))
DFg = vec_flag[int(flag_index)-1]
Xm1 = np.mean(X1, axis=0)
Xm = np.mean(X, axis=0)
J = -D * (Xm - Xm1) / np.linalg.norm(Xss - X1[u, :] + eps)
DRF = np.exp(-J / TF)
MS = np.exp(-FSss / (FS1[u] + eps))
R1 = np.random.rand(dim)
Qg = DFg * DRF * R1
X1new[u, :] = Xss + Qg * X1[u, :] + Qg * (MS * Xss - X1[u, :])
for u in range(n2):
Xm1 = np.mean(X1, axis=0)
Xm = np.mean(X, axis=0)
J = -D * (Xm1 - Xm) / np.linalg.norm(Xss - X2[u, :] + eps)
for j in range(n1):
FU = X1new[j, :] > ub
FL = X1new[j, :] < lb
X1new[j, :] = X1new[j, :] * ~(FU + FL) + ub * FU + lb * FL
v = CostFunction(X1new[j, :])
if v < FS1[j]:
FS1[j] = v
X1[j, :] = X1new[j, :]
for j in range(n2):
FU = X2new[j, :] > ub
FL = X2new[j, :] < lb
X2new[j, :] = X2new[j, :] * ~(FU + FL) + ub * FU + lb * FL
v = CostFunction(X2new[j, :])
if v < FS2[j]:
FS2[j] = v
X2[j, :] = X2new[j, :]
FSeo1, IndexFSeo1 = np.min(FS1), np.argmin(FS1)
FSeo2, IndexFSeo2 = np.min(FS2), np.argmin(FS2)
Xeo1 = X1[IndexFSeo1, :]
Xeo2 = X2[IndexFSeo2, :]
if FSeo1 < FSeo2:
FSss = FSeo1
YSol = Xeo1
else:
FSss = FSeo2
YSol = Xeo2
CNVG.append(FSss)
if FSss < BestF:
BestF = FSss
Xss = YSol
# print(f"iteration {t}: Best Cost Bw = {BestF}")
return Xss, BestF, CNVG
(8)调用FLA函数实现优化
接下来看最后的代码,使用FLA求解给定的优化问题。首先,通过定义成本函数(my_cost)和设置算法的参数(dim、Max_iteration、NoMolecules、lb、ub等),将成本函数传递给FLA算法。接着,调用FLA函数得到最优解 Xfood、最小成本 Xvalue 以及迭代过程中的成本变化 CNVG。最后,通过绘制迭代过程中成本的变化曲线展示算法的优化性能。
#%%
CostFunction = lambda x: my_cost(x) # Cost Function
dim = 20
Max_iteration = 400
NoMolecules = 30
lb = 0
ub = 4
eps = 2.2204*(10**-16)
[Xfood, Xvalue, CNVG] = FLA(NoMolecules, Max_iteration, lb, ub, dim, CostFunction)
plt.figure()
plt.plot(CNVG, color='r')
plt.xlim([1, Max_iteration])
plt.show()
在上述代码中,使用plt.plot(CNVG, color='r')语句绘制了迭代过程中成本的变化曲线。如图4-9所示,这条曲线反映了FLA算法在每次迭代中优化目标函数的效果。
图4-9 成本的变化曲线图
通过如下代码绘制了FLA算法得到的最优解 Xfood 对应的路径规划图,如图4-10所示。该图包括了无人机路径、起始点、目标点以及环境中的障碍物,展示了FLA算法优化后的路径规划结果。
PlotSolution(Xfood)
图4-10 FLA最优解 Xfood 对应的路径规划图
通过如下代码绘制了FLA算法得到的最优解 Xfood 对应的三维路径规划图,如图4-8所示。这幅图展示了无人机在三维空间中的路径,同时包括了起始点、目标点以及环境中的障碍物。通过这个图,可以更全面地观察无人机路径规划在三维空间中的优化效果。
PlotSolution3D(Xfood)
图4-8 FLA最优解 Xfood 对应的三维路径规划图