数学建模学习-多目标优化(Multi-objective Optimization)教程(35)
文章目录
写在最前
注意本文的相关代码及例子为同学们提供参考，借鉴相关结构，在这里举一些通俗易懂的例子，方便同学们根据实际情况修改代码，很多同学私信反映能否添加一些可视化，这里每篇教程都尽可能增加一些可视化方便同学理解，但具体使用时，同学们要根据实际情况选择是否在论文中添加可视化图片。
系列教程计划持续更新，同学们可以免费订阅专栏，内容充足后专栏可能付费，提前订阅的同学可以免费阅读，同时相关代码获取可以关注博主评论或私信。
算法简介
多目标优化(Multi-objective Optimization)是数学建模中一个重要的研究方向，它处理的是同时优化多个目标函数的问题。在实际工程中，我们经常需要在多个相互冲突的目标之间寻找平衡，例如:
- 产品设计中的成本与性能
- 投资组合中的收益与风险
- 机器学习中的模型复杂度与预测准确率
本教程介绍一种经典的多目标优化算法 - NSGA-II(Non-dominated Sorting Genetic Algorithm II，非支配排序遗传算法II)。该算法由Deb等人于2002年提出，是目前最流行的多目标优化算法之一。
算法特点
NSGA-II算法具有以下特点：
-
快速非支配排序
- 使用O(MN²)的时间复杂度进行非支配排序，其中M是目标函数个数，N是种群大小
- 能够有效处理大规模种群
-
精英保留策略
- 采用父代和子代种群合并的方式
- 确保优秀个体不会在进化过程中丢失
-
拥挤度距离
- 引入拥挤度距离来保持种群的多样性
- 避免解集中在Pareto前沿的某个局部区域
-
无需用户定义参数
- 不需要指定权重或其他用户定义的参数
- 能够自动找到均匀分布的Pareto最优解集
算法原理
NSGA-II算法的核心思想是通过非支配排序和拥挤度距离计算来评价个体的优劣。其主要步骤如下：
-
非支配排序
- 对种群中的所有个体进行非支配排序
- 将个体分成不同的等级(前沿)
- 第一等级的个体不被任何其他个体支配
- 第二等级的个体仅被第一等级的个体支配
- 以此类推
-
拥挤度距离计算
- 对同一等级的个体计算拥挤度距离
- 拥挤度距离表示个体与其邻近个体的距离
- 距离越大表示该区域的解越稀疏
-
选择操作
- 基于锦标赛选择
- 优先选择等级高的个体
- 当等级相同时，选择拥挤度距离大的个体
-
遗传操作
- 使用模拟二进制交叉(SBX)
- 使用多项式变异
- 生成新的子代种群
环境准备
本教程使用Python实现NSGA-II算法，需要以下库：
import numpy as np
import matplotlib.pyplot as plt
代码实现
1. 问题定义
首先，我们定义一个简单的双目标优化问题：
def objectives(self, x):
"""计算目标函数值"""
f1 = x**2 # 第一个目标：最小化x^2
f2 = (x-2)**2 # 第二个目标：最小化(x-2)^2
return np.array([f1, f2])
2. 种群初始化
def initialize_population(self):
"""初始化种群"""
pop = np.random.uniform(
self.bounds[0][0],
self.bounds[0][1],
(self.pop_size, self.n_variables)
)
return pop
3. 非支配排序
def non_dominated_sort(self, population):
"""非支配排序"""
n_pop = len(population)
domination_count = np.zeros(n_pop) # 支配个体i的个体数量
dominated_solutions = [[] for _ in range(n_pop)] # 个体i支配的个体集合
fronts = [[]] # 非支配等级
# 计算支配关系
for i in range(n_pop):
for j in range(i+1, n_pop):
obj_i = self.objectives(population[i])
obj_j = self.objectives(population[j])
if all(obj_i <= obj_j) and any(obj_i < obj_j): # i支配j
dominated_solutions[i].append(j)
domination_count[j] += 1
elif all(obj_j <= obj_i) and any(obj_j < obj_i): # j支配i
dominated_solutions[j].append(i)
domination_count[i] += 1
# 找出第一个非支配前沿
for i in range(n_pop):
if domination_count[i] == 0:
fronts[0].append(i)
# 生成所有非支配前沿
i = 0
while fronts[i]:
next_front = []
for j in fronts[i]:
for k in dominated_solutions[j]:
domination_count[k] -= 1
if domination_count[k] == 0:
next_front.append(k)
i += 1
if next_front:
fronts.append(next_front)
return fronts
4. 拥挤度距离计算
def crowding_distance(self, population, front):
"""计算拥挤度距离"""
if len(front) <= 2:
return np.full(len(front), np.inf)
distances = np.zeros(len(front))
for m in range(self.n_objectives):
# 获取该前沿个体的第m个目标函数值
values = np.array([self.objectives(population[i])[m] for i in front])
sorted_indices = np.argsort(values)
# 边界点设为无穷大
distances[sorted_indices[0]] = np.inf
distances[sorted_indices[-1]] = np.inf
# 计算中间点的拥挤度距离
for i in range(1, len(front)-1):
distances[sorted_indices[i]] += (
values[sorted_indices[i+1]] - values[sorted_indices[i-1]]
) / (values.max() - values.min())
return distances
5. 选择操作
def tournament_selection(self, population, fronts):
"""锦标赛选择"""
selected = np.zeros((self.pop_size, self.n_variables))
for i in range(self.pop_size):
# 随机选择两个个体
a, b = np.random.randint(0, self.pop_size, 2)
# 找到两个个体所在的前沿
a_front = next(j for j, front in enumerate(fronts) if a in front)
b_front = next(j for j, front in enumerate(fronts) if b in front)
# 根据非支配等级和拥挤度距离选择更好的个体
if a_front < b_front:
selected[i] = population[a]
elif b_front < a_front:
selected[i] = population[b]
else:
# 如果在同一前沿，选择拥挤度距离大的
a_crowd = self.crowding_distance(population, fronts[a_front])[list(fronts[a_front]).index(a)]
b_crowd = self.crowding_distance(population, fronts[b_front])[list(fronts[b_front]).index(b)]
if a_crowd > b_crowd:
selected[i] = population[a]
else:
selected[i] = population[b]
return selected
6. 遗传操作
def crossover(self, parent1, parent2):
"""模拟二进制交叉"""
nc = 20 # 交叉分布指数
u = np.random.random()
if u <= 0.5:
beta = (2*u)**(1/(nc+1))
else:
beta = (1/(2*(1-u)))**(1/(nc+1))
child1 = 0.5*((1+beta)*parent1 + (1-beta)*parent2)
child2 = 0.5*((1-beta)*parent1 + (1+beta)*parent2)
# 确保在边界内
child1 = np.clip(child1, self.bounds[0][0], self.bounds[0][1])
child2 = np.clip(child2, self.bounds[0][0], self.bounds[0][1])
return child1, child2
def mutation(self, individual):
"""多项式变异"""
nm = 20 # 变异分布指数
r = np.random.random()
if r < 0.5:
delta = (2*r)**(1/(nm+1)) - 1
else:
delta = 1 - (2*(1-r))**(1/(nm+1))
mutated = individual + delta*(self.bounds[0][1] - self.bounds[0][0])
mutated = np.clip(mutated, self.bounds[0][0], self.bounds[0][1])
return mutated
算法运行结果
我们以一个简单的双目标优化问题为例，目标是同时最小化f₁(x)=x²和f₂(x)=(x-2)²。这两个目标是相互冲突的，因为f₁的最优解是x=0，而f₂的最优解是x=2。
运行算法200代后，我们可以看到Pareto前沿的演化过程：
-
初始种群(第0代)：
[外链图片转存中…(img-vmYXJh5O-1737349350276)]
-
第50代：
[外链图片转存中…(img-QDLvAMV9-1737349359032)]
-
第100代：
[外链图片转存中…(img-JYZjKRWx-1737349367777)]
-
第150代：
-
最终结果(第199代)：
[外链图片转存中…(img-yaVcmBU6-1737349385423)]
从结果可以看出：
6. 算法成功找到了一系列非支配解，形成了清晰的Pareto前沿
7. 解集在目标空间中分布均匀，说明算法保持了良好的多样性
8. 前沿的形状符合理论预期，验证了算法的正确性
全部代码
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import matplotlib as mpl
import os
# 设置中文字体
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
mpl.rcParams['axes.unicode_minus'] = False
# 获取当前文件所在目录的绝对路径
CURRENT_DIR = os.path.dirname(os.path.abspath(__file__))
IMAGES_DIR = os.path.join(CURRENT_DIR, 'images')
# 确保images目录存在
if not os.path.exists(IMAGES_DIR):
os.makedirs(IMAGES_DIR)
class NSGA2:
def __init__(self, pop_size=100, n_generations=200):
self.pop_size = pop_size # 种群大小
self.n_generations = n_generations # 迭代次数
self.n_objectives = 2 # 目标函数个数
self.n_variables = 1 # 决策变量个数
self.bounds = [(0, 1)] # 决策变量范围
def objectives(self, x):
"""计算目标函数值"""
x = x.flatten()[0] # 确保x是标量
f1 = x**2 # 第一个目标：最小化x^2
f2 = (x-2)**2 # 第二个目标：最小化(x-2)^2
return np.array([f1, f2])
def initialize_population(self):
"""初始化种群"""
pop = np.random.uniform(
self.bounds[0][0],
self.bounds[0][1],
(self.pop_size, self.n_variables)
)
return pop
def non_dominated_sort(self, population):
"""非支配排序"""
n_pop = len(population)
domination_count = np.zeros(n_pop) # 支配个体i的个体数量
dominated_solutions = [[] for _ in range(n_pop)] # 个体i支配的个体集合
fronts = [[]] # 非支配等级
# 计算支配关系
for i in range(n_pop):
for j in range(i+1, n_pop):
obj_i = self.objectives(population[i])
obj_j = self.objectives(population[j])
if all(obj_i <= obj_j) and any(obj_i < obj_j): # i支配j
dominated_solutions[i].append(j)
domination_count[j] += 1
elif all(obj_j <= obj_i) and any(obj_j < obj_i): # j支配i
dominated_solutions[j].append(i)
domination_count[i] += 1
# 找出第一个非支配前沿
for i in range(n_pop):
if domination_count[i] == 0:
fronts[0].append(i)
# 生成所有非支配前沿
i = 0
while i < len(fronts) and fronts[i]: # 修改这里的条件
next_front = []
for j in fronts[i]:
for k in dominated_solutions[j]:
domination_count[k] -= 1
if domination_count[k] == 0:
next_front.append(k)
i += 1
if next_front:
fronts.append(next_front)
return fronts
def crowding_distance(self, population, front):
"""计算拥挤度距离"""
if len(front) <= 2:
return np.full(len(front), np.inf)
distances = np.zeros(len(front))
for m in range(self.n_objectives):
# 获取该前沿个体的第m个目标函数值
values = np.array([self.objectives(population[i])[m] for i in front])
sorted_indices = np.argsort(values)
# 边界点设为无穷大
distances[sorted_indices[0]] = np.inf
distances[sorted_indices[-1]] = np.inf
# 计算中间点的拥挤度距离
for i in range(1, len(front)-1):
distances[sorted_indices[i]] += (
values[sorted_indices[i+1]] - values[sorted_indices[i-1]]
) / (values.max() - values.min())
return distances
def tournament_selection(self, population, fronts):
"""锦标赛选择"""
selected = np.zeros((self.pop_size, self.n_variables))
for i in range(self.pop_size):
# 随机选择两个个体
a, b = np.random.randint(0, self.pop_size, 2)
# 找到两个个体所在的前沿
a_front = next(j for j, front in enumerate(fronts) if a in front)
b_front = next(j for j, front in enumerate(fronts) if b in front)
# 根据非支配等级和拥挤度距离选择更好的个体
if a_front < b_front:
selected[i] = population[a]
elif b_front < a_front:
selected[i] = population[b]
else:
# 如果在同一前沿，选择拥挤度距离大的
a_crowd = self.crowding_distance(population, fronts[a_front])[list(fronts[a_front]).index(a)]
b_crowd = self.crowding_distance(population, fronts[b_front])[list(fronts[b_front]).index(b)]
if a_crowd > b_crowd:
selected[i] = population[a]
else:
selected[i] = population[b]
return selected
def crossover(self, parent1, parent2):
"""模拟二进制交叉"""
nc = 20 # 交叉分布指数
u = np.random.random()
if u <= 0.5:
beta = (2*u)**(1/(nc+1))
else:
beta = (1/(2*(1-u)))**(1/(nc+1))
child1 = 0.5*((1+beta)*parent1 + (1-beta)*parent2)
child2 = 0.5*((1-beta)*parent1 + (1+beta)*parent2)
# 确保在边界内
child1 = np.clip(child1, self.bounds[0][0], self.bounds[0][1])
child2 = np.clip(child2, self.bounds[0][0], self.bounds[0][1])
return child1, child2
def mutation(self, individual):
"""多项式变异"""
nm = 20 # 变异分布指数
r = np.random.random()
if r < 0.5:
delta = (2*r)**(1/(nm+1)) - 1
else:
delta = 1 - (2*(1-r))**(1/(nm+1))
mutated = individual + delta*(self.bounds[0][1] - self.bounds[0][0])
mutated = np.clip(mutated, self.bounds[0][0], self.bounds[0][1])
return mutated
def evolve(self):
"""进化过程"""
# 初始化种群
population = self.initialize_population()
# 记录每代的非支配解
all_fronts = []
# 主循环
for gen in range(self.n_generations):
# 非支配排序
fronts = self.non_dominated_sort(population)
all_fronts.append(fronts[0]) # 记录当前代的非支配解
# 选择
selected = self.tournament_selection(population, fronts)
# 交叉和变异
offspring = np.zeros_like(selected)
for i in range(0, self.pop_size, 2):
if np.random.random() < 0.9: # 交叉概率
offspring[i], offspring[i+1] = self.crossover(
selected[i], selected[i+1]
)
else:
offspring[i] = selected[i]
offspring[i+1] = selected[i+1]
if np.random.random() < 0.1: # 变异概率
offspring[i] = self.mutation(offspring[i])
if np.random.random() < 0.1:
offspring[i+1] = self.mutation(offspring[i+1])
# 合并父代和子代
population = offspring
# 每隔50代绘制一次当前的Pareto前沿
if gen % 50 == 0 or gen == self.n_generations - 1:
self.plot_pareto_front(population, gen)
print(f"Generation {gen}: Plotting Pareto front...")
return population, all_fronts[-1]
def plot_pareto_front(self, population, generation):
"""绘制Pareto前沿"""
plt.figure(figsize=(10, 8))
# 绘制所有解
f1 = [self.objectives(x)[0] for x in population]
f2 = [self.objectives(x)[1] for x in population]
plt.scatter(f1, f2, c='blue', label='所有解')
# 绘制非支配解
front = self.non_dominated_sort(population)[0]
f1_pareto = [self.objectives(population[i])[0] for i in front]
f2_pareto = [self.objectives(population[i])[1] for i in front]
plt.scatter(f1_pareto, f2_pareto, c='red', label='非支配解')
plt.xlabel('目标函数1 (f1)', fontproperties=font)
plt.ylabel('目标函数2 (f2)', fontproperties=font)
plt.title(f'第{generation}代Pareto前沿', fontproperties=font)
plt.legend(prop=font)
plt.grid(True)
# 保存图片
save_path = os.path.join(IMAGES_DIR, f'pareto_front_gen_{generation}.png')
plt.savefig(save_path)
plt.close()
def main():
# 创建NSGA-II实例
nsga2 = NSGA2(pop_size=100, n_generations=200)
# 运行算法
final_pop, final_front = nsga2.evolve()
# 输出最终的非支配解
print("\n最终的非支配解:")
for i in final_front:
x = final_pop[i]
objectives = nsga2.objectives(x)
print(f"x = {x[0]:.4f}, f1 = {objectives[0]:.4f}, f2 = {objectives[1]:.4f}")
if __name__ == "__main__":
main()
应用场景
NSGA-II算法在实际工程中有广泛的应用，例如：
-
工程设计优化
- 结构优化：强度vs重量
- 电路设计：功耗vs性能
- 控制系统：响应速度vs稳定性
-
金融投资
- 投资组合优化：收益vs风险
- 资产配置：流动性vs收益率
-
机器学习
4 - 神经网络结构优化：准确率vs计算复杂度- 特征选择：特征数量vs模型性能
-
供应链管理
- 库存优化：库存成本vs服务水平
- 运输规划：运输成本vs交付时间
算法优化建议
-
参数调优
- 种群大小：根据问题规模和计算资源调整
- 迭代次数：根据收敛情况确定
- 交叉变异概率：影响种群多样性
-
目标函数设计
- 合理设计目标函数的量纲和尺度
- 考虑目标之间的权衡关系
-
约束处理
- 添加惩罚项处理约束
- 使用修复策略保持解的可行性
-
性能提升
- 使用并行计算加速种群评价
- 采用自适应参数调整策略
- 结合问题特点设计专门的遗传算子
总结
NSGA-II算法是一种强大的多目标优化工具，它能够：
- 同时处理多个目标函数
- 自动保持种群多样性
- 不需要预先设定目标权重
- 得到一系列均匀分布的非支配解
在实际应用中，需要根据具体问题特点进行适当的调整和优化。同时，理解算法的核心思想和原理对于正确使用和改进算法至关重要。
同学们如果有疑问可以私信答疑，如果有讲的不好的地方或可以改善的地方可以一起交流，谢谢大家。
扫一扫
1625
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
