引力搜索算法

引力搜索算法过程，包括了初始化、适应度评估、质量计算、加速度计算、更新速度和位置的一些步骤。

import numpy as np
import random as rd
from math import exp, sqrt
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation# 定义目标函数
def objFunction(x1, x2):return x1 ** 2 + x2 ** 2# 初始化种群的位置和速度
def init(n):position, velocity = [], []for i in range(n):X1 = rd.uniform(-10, 10)X2 = rd.uniform(-10, 10)V1 = rd.uniform(-3, 3)V2 = rd.uniform(-3, 3)position.append([X1, X2])velocity.append([V1, V2])return position, velocity# 计算适应度值
def fitnessEva(position):fitness = []for i in range(len(position)):fitness.append(objFunction(position[i][0], position[i][1]))return fitness# 找到最佳和最差的适应度值
def findBestAndWorst(position):return min(fitnessEva(position)), max(fitnessEva(position))# 计算每个个体的质量
def calculateMass(fitness):mass = []Mass = []for i in range(len(fitness)):mass.append((fitness[i] - max(fitness)) / (min(fitness) - max(fitness)))for i in range(len(mass)):Mass.append(mass[i] / sum(mass))return Mass# 计算每个个体的加速度
def calculateAcceleration(position, Mass, G, topK):acceleration = []Fi0, Fi1 = 0, 0for i in range(len(position)):for j in range(len(position)):if i != j and j in topK:Fi0 += rd.random() * G * ((Mass[i] * Mass[j]) / (calculateDistance(position[i], position[j]) + r)) * (position[j][0] - position[i][0])Fi1 += rd.random() * G * ((Mass[i] * Mass[j]) / (calculateDistance(position[i], position[j]) + r)) * (position[j][1] - position[i][1])if Mass[i] != 0:acceleration.append([Fi0 / Mass[i] / 10, Fi1 / Mass[i] / 10])  # 这里除10是为了避免粒子的加速度过大else:acceleration.append([10, 10])Fi0 = 0Fi1 = 0return acceleration# 找出适应度更优的前K个粒子
def findTopK(fitness, K):topK = []dic = {}for i in range(len(fitness)):dic[i] = fitness[i]fitness.sort()for i in range(K):topK.append(list(dic.keys())[list(dic.values()).index(fitness[i])])return topK# 更新速度和位置
def updateVelocityAndPosition(acceleration, position, velocity):for i in range(len(velocity)):velocity[i][0] = rd.random() * velocity[i][0] + acceleration[i][0]velocity[i][1] = rd.random() * velocity[i][1] + acceleration[i][1]position[i][0] = position[i][0] + velocity[i][0]position[i][1] = position[i][1] + velocity[i][1]# 计算两个点之间的距离
def calculateDistance(p1, p2):return sqrt((p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2)# 检查位置是否在可行域内
def checkPosition(position):for i in range(len(position)):if position[i][0] < -10:position[i][0] = -10elif position[i][0] > 10:position[i][0] = 10if position[i][1] < -10:position[i][1] = -10elif position[i][1] > 10:position[i][1] = 10if __name__ == '__main__':G = 100r = 1K = 50iterx, maxIterx = 0, 50position, velocity = init(50)fig = plt.figure(figsize=(10, 8))ax = fig.add_subplot(111, projection='3d')scatter = ax.scatter([], [], [], c='pink', marker='o')path_lines = [ax.plot([], [], [], color='lightyellow')[0] for _ in range(50)]ax.set_xlim(-10, 10)ax.set_ylim(-10, 10)ax.set_zlim(0, 200)ax.set_title('Gravitational Search Algorithm Visualization')ax.set_xlabel('X1')ax.set_ylabel('X2')ax.set_zlabel('Fitness')def animate(iteration, position, velocity, G, K, path_history):fitness = fitnessEva(position)  # 适应性评估G = G * exp(-20 * iteration / maxIterx)  # 更新引力常量Mass = calculateMass(fitness)  # 更新粒子质量topK = findTopK(fitness, K)  # 找出适应度更优的前K个粒子acceleration = calculateAcceleration(position, Mass, G, topK)  # 计算粒子加速度updateVelocityAndPosition(acceleration, position, velocity)  # 根据加速度更新速度与位置checkPosition(position)  # 检查粒子是否冲出了解空间K = K - iteration  # 更新K值# Update scatter plotscatter._offsets3d = ([pos[0] for pos in position],[pos[1] for pos in position],fitness)# Update path historyfor i in range(len(position)):path_history[i].append(position[i] + [fitness[i]])path_lines[i].set_data(np.array(path_history[i]).T[:2])path_lines[i].set_3d_properties(np.array(path_history[i]).T[2])best_fitness = min(fitnessEva(position))best_position = position[fitnessEva(position).index(best_fitness)]ax.set_title(f'Iteration {iteration}: Best Fitness = {best_fitness:.4f}')return scatter,path_history = [[] for _ in range(50)]ani = FuncAnimation(fig, lambda frame: animate(frame, position, velocity, G, K, path_history), frames=maxIterx,interval=200, blit=True, repeat=False)plt.show()best_fitness = min(fitnessEva(position))best_position = position[fitnessEva(position).index(best_fitness)]print("最优解:", best_fitness)print("最优解对应的位置:", best_position)