概述
多目标粒子群优化（MOPSO） 是粒子群优化（PSO）的一种扩展，用于解决具有多个目标函数的优化问题。MOPSO的目标是找到一组非支配解（Pareto最优解），这些解在不同目标之间达到平衡。
主要步骤
- 初始化粒子群：
随机生成一组粒子，每个粒子有一个初始位置和速度。 - 评估适应度：
计算每个粒子在所有目标函数上的适应度值。
更新粒子的个人最佳位置和个人最佳适应度。 - 非支配解选择：
使用支配关系判断粒子之间的优劣。
找到一组非支配解，构成当前的Pareto前沿。 - 更新速度和位置：
根据认知项（粒子自身的最佳位置）和社会项（全局最佳位置）更新粒子的速度。
根据更新后的速度调整粒子的位置。 - 迭代优化：
重复执行适应度评估、非支配解选择和速度位置更新，直到达到最大迭代次数或满足其他停止条件。
代码功能
该代码实现了一个多目标粒子群优化（Multi-Objective Particle Swarm Optimization, MOPSO）算法，并使用Python编写。MOPSO算法用于解决具有多个目标函数的优化问题。具体来说，该代码实现了以下功能：
- 粒子群初始化：
创建一组粒子，每个粒子有一个初始位置和速度。
粒子的位置在给定的边界范围内随机生成，初始速度设为零。 - 适应度评估：
对每个粒子计算其适应度值，即每个目标函数在该粒子位置上的值。
更新粒子的个人最佳位置和个人最佳适应度。 - 速度和位置更新：
根据认知项（粒子自身的最佳位置）和社会项（全局最佳位置）更新粒子的速度。
根据新的速度更新粒子的位置，并确保位置在搜索空间内。 - 非支配解选择：
使用支配关系判断粒子之间的优劣，找到一组非支配解。
随机选择一个非支配解作为全局最佳位置。 - 多目标优化：
在指定的迭代次数内重复执行适应度评估、速度和位置更新以及非支配解选择。
最终返回找到的所有非支配解的位置。 - 结果可视化：
计算每个非支配解在两个目标函数上的值。
使用Matplotlib库绘制这些解在两个目标函数值上的分布图，展示Pareto前沿。
代码
import numpy as np
import random
import matplotlib.pyplot as plt
class Particle:
def __init__(self, position, velocity):
self.position = position # 粒子的位置
self.velocity = velocity # 粒子的速度
self.best_position = position.copy() # 粒子的个人最佳位置
self.best_fitness = None # 粒子的个人最佳适应度
self.fitness = None # 粒子的当前适应度
def initialize_population(dimensions, population_size, bounds):
# 初始化粒子群
particles = []
for _ in range(population_size):
position = [random.uniform(bounds[i][0], bounds[i][1]) for i in range(dimensions)] # 随机生成初始位置
velocity = [0.0 for _ in range(dimensions)] # 初始速度为0
particles.append(Particle(position, velocity))
return particles
def evaluate_fitness(particle, objectives):
# 计算粒子的适应度
particle.fitness = [objective(particle.position) for objective in objectives] # 计算当前适应度
if particle.best_fitness is None or dominates(particle.fitness, particle.best_fitness):
# 如果当前适应度优于个人最佳适应度，则更新个人最佳位置和适应度
particle.best_fitness = particle.fitness.copy()
particle.best_position = particle.position.copy()
def update_velocity(particle, global_best_positions, w=0.7, c1=1.5, c2=1.5):
# 更新粒子的速度
for i in range(len(particle.velocity)):
r1 = random.random()
r2 = random.random()
cognitive = c1 * r1 * (particle.best_position[i] - particle.position[i]) # 认知项
social = c2 * r2 * (global_best_positions[i] - particle.position[i]) # 社会项
particle.velocity[i] = w * particle.velocity[i] + cognitive + social # 更新速度
def update_position(particle, bounds):
# 更新粒子的位置
for i in range(len(particle.position)):
particle.position[i] += particle.velocity[i]
# 确保粒子的位置在搜索空间内
if particle.position[i] < bounds[i][0]:
particle.position[i] = bounds[i][0]
elif particle.position[i] > bounds[i][1]:
particle.position[i] = bounds[i][1]
def dominates(p1, p2):
# 检查p1是否支配p2
return all(f1 <= f2 for f1, f2 in zip(p1, p2)) and any(f1 < f2 for f1, f2 in zip(p1, p2))
def find_global_best(particles):
# 找到全局最佳位置
non_dominated_sorting = []
for particle in particles:
is_dominated = False
for other in particles:
if dominates(other.fitness, particle.fitness):
is_dominated = True
break
if not is_dominated:
non_dominated_sorting.append(particle)
if non_dominated_sorting:
# 随机选择一个非支配解作为全局最佳位置
return random.choice(non_dominated_sorting).position
else:
return particles[0].position
def mopso(objectives, dimensions, bounds, population_size, iterations):
# 多目标粒子群优化算法
particles = initialize_population(dimensions, population_size, bounds)
for particle in particles:
evaluate_fitness(particle, objectives)
for _ in range(iterations):
global_best_positions = find_global_best(particles)
for particle in particles:
update_velocity(particle, global_best_positions)
update_position(particle, bounds)
evaluate_fitness(particle, objectives)
# 返回找到的最佳位置
return [p.position for p in particles]
# 示例用法
def obj1(x):
return sum(xi**2 for xi in x) # 目标函数1
def obj2(x):
return sum((xi-2)**2 for xi in x) # 目标函数2
dimensions = 2 # 维度数
bounds = [(-5.12, 5.12)] * dimensions # 每个维度的边界
population_size = 50 # 粒子群大小
iterations = 1000 # 迭代次数
solutions = mopso([obj1, obj2], dimensions, bounds, population_size, iterations)
# 计算最终解的适应度
final_fitness = [(obj1(sol), obj2(sol)) for sol in solutions]
# 绘制结果
plt.scatter([f1 for f1, f2 in final_fitness], [f2 for f1, f2 in final_fitness])
plt.xlabel('Objective 1')
plt.ylabel('Objective 2')
plt.title('Pareto Front')
plt.grid(True)
plt.show()
