import time
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
from scipy.spatial.distance import cdist
from scipy.stats import qmc
# 全局字体设置,解决中文和特殊符号显示问题
plt.rcParams["font.family"] = ["SimHei", "Arial Unicode MS"]
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
# 1. 目标函数向量化:支持批量解输入(核心优化点1)
def zdt1(x):
"""ZDT1测试函数(支持2D数组输入:n行×dim列)"""
x = np.atleast_2d(x) # 确保输入是2D数组
n, dim = x.shape
f1 = x[:, 0] # 所有解的第0维(批量取列)
g = 1 + 9 * np.mean(x[:, 1:], axis=1) # 每行(每个解)的均值
h = 1 - np.sqrt(f1 / g)
f2 = g * h
return np.column_stack((f1, f2)) # 输出n行×2列
def zdt2(x):
"""ZDT2测试函数(支持批量解输入)"""
x = np.atleast_2d(x)
n, dim = x.shape
f1 = x[:, 0]
g = 1 + 9 * np.mean(x[:, 1:], axis=1)
h = 1 - (f1 / g) ** 2
f2 = g * h
return np.column_stack((f1, f2))
def zdt3(x):
"""ZDT3测试函数(支持批量解输入)"""
x = np.atleast_2d(x)
n, dim = x.shape
f1 = x[:, 0]
g = 1 + 9 * np.mean(x[:, 1:], axis=1)
h = 1 - np.sqrt(f1 / g) - (f1 / g) * np.sin(10 * np.pi * f1)
return np.column_stack((f1, g * h))
def zdt4(x):
"""严格修正的ZDT4实现"""
x = np.atleast_2d(x)
n, dim = x.shape
f1 = x[:, 0]
# 严格遵循原始论文公式
g = 1 + 10 * (dim - 1) + np.sum(
x[:, 1:] ** 2 - 10 * np.cos(4 * np.pi * x[:, 1:]),
axis=1
)
h = 1 - np.sqrt(f1 / g)
f2 = g * h
return np.column_stack((f1, f2))
def zdt6(x):
"""ZDT6测试函数(支持批量解输入)"""
x = np.atleast_2d(x)
n, dim = x.shape
# ZDT6公式:f1 = 1 - exp(-4x1)sin^6(6πx1)
f1 = 1 - np.exp(-4 * x[:, 0]) * (np.sin(6 * np.pi * x[:, 0])) ** 6
# g(x) = 1 + 9*(sum(x2~xd)/(d-1))^0.25
g = 1 + 9 * (np.sum(x[:, 1:], axis=1) / (dim - 1)) ** 0.25
# h = 1 - (f1/g)^2
h = 1 - (f1 / g) ** 2
f2 = g * h
return np.column_stack((f1, f2))
# --------------------------- 新增:二维DTLZ1/2/3目标函数 ---------------------------
def dtlz1(x):
"""
二维DTLZ1测试函数(M=2目标,决策变量n=3)
特性:线性Pareto前沿、多模态
决策变量范围:x ∈ [0,1]^3
真实PF:f1 + f2 = 0.5(f1,f2 ≥ 0)
"""
x = np.atleast_2d(x) # 确保输入为2D数组(n行×3列)
n, dim = x.shape
assert dim == 3, "DTLZ1(二维目标)需输入3个决策变量(n=3)"
# 向量化计算g(x),确保形状为(n,)
g = 100 * (
2 + # 3-2+1=2(n-M+1)
(x[:, 2] - 0.5) ** 2 - np.cos(20 * np.pi * (x[:, 2] - 0.5))
)
# 目标函数计算(向量化操作,避免广播错误)
f1 = 0.5 * x[:, 0] * x[:, 1] * (1 + g)
f2 = 0.5 * (1 - x[:, 0]) * (1 + g)
return np.column_stack((f1, f2))
def dtlz2(x):
"""
二维DTLZ2测试函数(M=2目标,决策变量n=3)
特性:凸形Pareto前沿、非线性目标关联
决策变量范围:x ∈ [0,1]^3
真实PF:f1² + f2² = 1(f1,f2 ≥ 0)
"""
x = np.atleast_2d(x) # 确保输入为2D数组(n行×3列)
n, dim = x.shape
assert dim == 3, "DTLZ2(二维目标)需输入3个决策变量(n=3)"
# 向量化计算g(x)(已正确实现,保持并优化注释)
g = np.sum((x[:, 2] - 0.5) ** 2, axis=0) # 修正:对于单个变量无需axis=1,避免降维问题
# 目标函数计算
f1 = (1 + g) * np.cos(x[:, 0] * np.pi / 2) * np.cos(x[:, 1] * np.pi / 2)
f2 = (1 + g) * np.cos(x[:, 0] * np.pi / 2) * np.sin(x[:, 1] * np.pi / 2)
return np.column_stack((f1, f2))
def dtlz3(x):
"""
二维DTLZ3测试函数(M=2目标,决策变量n=3)
特性:凸形Pareto前沿、强多模态
决策变量范围:x ∈ [0,1]^3
真实PF:f1² + f2² = 1(f1,f2 ≥ 0)
"""
x = np.atleast_2d(x) # 确保输入为2D数组(n行×3列)
n, dim = x.shape
assert dim == 3, "DTLZ3(二维目标)需输入3个决策变量(n=3)"
# 向量化计算g(x),确保形状为(n,)
g = 100 * (
2 + # 3-2+1=2(n-M+1)
(x[:, 2] - 0.5) ** 2 - np.cos(20 * np.pi * (x[:, 2] - 0.5))
)
# 目标函数计算(向量化操作)
f1 = (1 + g) * np.cos(x[:, 0] * np.pi / 2) * np.cos(x[:, 1] * np.pi / 2)
f2 = (1 + g) * np.cos(x[:, 0] * np.pi / 2) * np.sin(x[:, 1] * np.pi / 2)
return np.column_stack((f1, f2))
def generate_true_front(func, num_points=1000):
"""
生成ZDT测试函数真实Pareto前沿,包含ZDT6的实现
:param zdt_func: 目标函数(zdt1/zdt2/zdt3/zdt4/zdt6)
:param num_points: 采样点数(越大越精确)
:return: 真实前沿的(f1, f2)数组
"""
# 以下为原有其他ZDT函数的处理逻辑,保持不变
f1 = np.linspace(0, 1, num_points)
f2 = []
# --------------------------- 新增:DTLZ系列前沿生成 ---------------------------
if func == dtlz1:
# DTLZ1真实PF:f1 + f2 = 0.5(f1∈[0,0.5],f2=0.5-f1)
f1 = np.linspace(0, 0.5, num_points)
f2 = 0.5 - f1
elif func == dtlz2 or func == dtlz3:
# DTLZ2/3真实PF:f1² + f2² = 1(f1∈[0,1],f2=√(1-f1²))
# 注:两函数PF形状相同,仅决策空间多模态差异
theta = np.linspace(0, np.pi / 2, num_points) # 极角(0~90°)
f1 = np.cos(theta)
f2 = np.sin(theta)
if func == zdt1:
for f in f1:
f2_val = 1 - np.sqrt(f)
f2.append(f2_val)
elif func == zdt2:
for f in f1:
f2_val = 1 - f ** 2
f2.append(f2_val)
elif func == zdt3:
# ZDT3 真实前沿:g(x)=1(x2~xd=0),逐点计算并保留有效分段
for f in f1:
g = 1.0
h = 1 - np.sqrt(f / g) - (f / g) * np.sin(10 * np.pi * f)
f2_val = g * h
# 仅保留物理意义有效的点(过滤极端异常值)
if -1.5 <= f2_val <= 1.5:
f2.append(f2_val)
else:
f2.append(np.nan) # 无效点标记为NaN
elif func == zdt4:
for f in f1:
g_fixed = 1.0 # 固定g(x)=1
h = 1 - np.sqrt(f / g_fixed)
f2_val = g_fixed * h
f2.append(f2_val)
elif func == zdt6:
# ZDT6真实前沿生成:x2~xd=0时,g(x)=1
for x1 in np.linspace(0, 1, num_points):
f1_val = 1 - np.exp(-4 * x1) * (np.sin(6 * np.pi * x1)) ** 6
g = 1.0 # x2~xd=0时g(x)=1
h = 1 - (f1_val / g) ** 2
f2_val = g * h
f2.append(f2_val)
f1 = [1 - np.exp(-4 * x1) * (np.sin(6 * np.pi * x1)) ** 6 for x1 in np.linspace(0, 1, num_points)]
# 转换为numpy数组并过滤NaN
f2 = np.array(f2)
front = np.column_stack((f1[:len(f2)], f2))
# ZDT3 额外过滤被支配点,确保严格前沿
if func == zdt3:
to_keep = np.ones(len(front), dtype=bool)
for i in range(len(front)):
f1_i, f2_i = front[i]
for j in range(len(front)):
if i == j:
continue
f1_j, f2_j = front[j]
# 若点j支配点i,则过滤i
if f1_j <= f1_i and f2_j <= f2_i and (f1_j < f1_i or f2_j < f2_i):
to_keep[i] = False
break
front = front[to_keep]
return front
def calculate_igd(algorithm_front, true_front):
"""计算反向世代距离(IGD)"""
if len(algorithm_front) == 0:
return float('inf')
distances = cdist(true_front, algorithm_front).min(axis=1)
return np.mean(distances)
def initialization(pop_size, dim, lb, ub):
"""优化版种群初始化:支持固定维度"""
lb = np.array(lb) if not isinstance(lb, np.ndarray) else lb
ub = np.array(ub) if not isinstance(ub, np.ndarray) else ub
lb = np.repeat(lb, dim) if lb.size == 1 else lb
ub = np.repeat(ub, dim) if ub.size == 1 else ub
# 分离固定维度(上下界相等)和可采样维度
fixed_dims = lb == ub
sample_dims = ~fixed_dims
n_sample = np.sum(sample_dims)
# 初始化种群(固定维度先设为lb值)
population = np.full((pop_size, dim), lb)
# 仅对可采样维度做LHS采样
if n_sample > 0:
sampler = qmc.LatinHypercube(d=n_sample)
samples = sampler.random(n=pop_size)
# 映射到可采样维度的上下界
scaled_samples = qmc.scale(
samples,
lb[sample_dims],
ub[sample_dims]
)
population[:, sample_dims] = scaled_samples
return population
def Domination(X, FitX):
N = X.shape[0]
dominated_count = np.zeros(N, dtype=int)
# 生成所有解对的支配关系矩阵(向量化比较)
dominate_matrix = np.all(FitX[:, np.newaxis] <= FitX, axis=2) & np.any(FitX[:, np.newaxis] < FitX, axis=2)
# 统计每个解被多少个其他解支配
dominated_count = np.sum(dominate_matrix, axis=0)
# 筛选非支配解
mask = dominated_count == 0
return X[mask], FitX[mask]
def crowding_distance(fitness):
"""计算每个解的拥挤距离"""
n, m = fitness.shape
if n <= 2:
return np.ones(n) * np.inf # 解数量过少时,距离设为无穷大
crowd_dist = np.zeros(n)
for j in range(m):
sorted_idx = np.argsort(fitness[:, j]) # 按第j个目标排序
sorted_vals = fitness[sorted_idx, j]
# 边界解距离设为无穷大(优先保留)
crowd_dist[sorted_idx[0]] = np.inf
crowd_dist[sorted_idx[-1]] = np.inf
# 中间解计算距离
if sorted_vals[-1] != sorted_vals[0]: # 避免除零
crowd_dist[sorted_idx[1:-1]] += (sorted_vals[2:] - sorted_vals[:-2]) / (sorted_vals[-1] - sorted_vals[0])
return crowd_dist
def crowding_distance_crop(X, FitX, L):
"""使用拥挤距离裁剪种群到大小L"""
crowd_dist = crowding_distance(FitX)
sorted_indices = np.argsort(-crowd_dist) # 按距离降序
return X[sorted_indices[:L], :], FitX[sorted_indices[:L], :]
def UpdateArchive(NumObj, archive_size, A, FitA, X, FitX):
"""优化的存档更新函数"""
A = np.atleast_2d(A) if len(A) > 0 else np.empty((0, X.shape[1]))
FitA = np.atleast_2d(FitA) if len(FitA) > 0 else np.empty((0, NumObj))
X = np.atleast_2d(X)
FitX = np.atleast_2d(FitX)
if len(X) == 0:
return A, FitA
# 合并存档解和新解
combined_pop = np.vstack([A, X])
combined_fit = np.vstack([FitA, FitX])
total = len(combined_pop)
# 向量化计算支配关系矩阵
dominate_matrix = (
np.all(combined_fit[:, np.newaxis] <= combined_fit, axis=2) &
np.any(combined_fit[:, np.newaxis] < combined_fit, axis=2)
)
# 计算每个解的被支配计数
dominated_count = np.sum(dominate_matrix, axis=0)
# 筛选所有非支配解
non_dominated_mask = (dominated_count == 0)
non_dominated_pop = combined_pop[non_dominated_mask]
non_dominated_fit = combined_fit[non_dominated_mask]
# 按拥挤距离裁剪到存档大小
if len(non_dominated_pop) > archive_size:
non_dominated_pop, non_dominated_fit = crowding_distance_crop(
non_dominated_pop, non_dominated_fit, archive_size
)
return non_dominated_pop, non_dominated_fit
def soft_rime_search(rimepop, A, FitA, it, max_iter, lb, ub, dim, W=5):
"""简化版软霜搜索函数"""
if max_iter < 1:
raise ValueError("max_iter不能小于1")
N = rimepop.shape[0]
new_rimepop = rimepop.copy()
progress = it / max_iter # 迭代进度(0~1)
# 软霜因子
rime_factor = (np.random.rand() - 0.5) * 2 * np.cos((np.pi * it / (max_iter / 5))) * (
1 - np.round(it * W / max_iter) / W)
# 阶段式参数
E = np.sqrt(progress)
w = 1 / (1 + np.exp(-9 * (progress - 0.5)))
# 引导解选择
if len(A) > 0:
guides = A[np.random.choice(len(A), size=N)]
else:
guides = np.mean(rimepop, axis=0, keepdims=True).repeat(N, axis=0)
# 位置更新
for i in range(N):
global_guide = guides[i]
for j in range(dim):
if np.random.rand() < E:
# 正向引导(向引导解方向)
forward_dir = (1 + rime_factor) * (global_guide[j] - new_rimepop[i, j]) * np.random.rand()
# 反向引导
if progress < 0.5:
reverse_guide = lb[j] + ub[j] - new_rimepop[i, j]
rev_strength = 1 + (1 - progress) * (1 + rime_factor)
else:
arch_mean = np.mean(A[:, j]) if len(A) > 0 else np.mean(rimepop[:, j])
reverse_guide = arch_mean * 2 - new_rimepop[i, j]
rev_strength = 1 + (1 - progress)
# 反向引导强度
reverse_dir = rev_strength * (reverse_guide - new_rimepop[i, j]) * (np.random.rand() * 2 - 1)
new_rimepop[i, j] += (1 - w) * reverse_dir + w * forward_dir
return np.clip(new_rimepop, lb, ub) # 边界修正
def hard_rime_puncture(new_rimepop, A, FitA, objective_func, it, max_iter):
"""优化后的硬霜穿刺函数"""
N, dim = new_rimepop.shape
# 存档解太少时不穿刺,直接返回
if len(A) < 2:
new_FitX = objective_func(new_rimepop) # 使用向量化目标函数
return new_rimepop, new_FitX
# 计算新解的目标值(批量计算)
new_FitX = objective_func(new_rimepop)
# 向量化Pareto等级计算
def pareto_ranking(pop_fitness):
n = len(pop_fitness)
ranks = np.zeros(n, dtype=int)
dominated_count = np.zeros(n, dtype=int)
# 向量化构建支配关系矩阵
dominate_matrix = (
np.all(pop_fitness[:, np.newaxis] <= pop_fitness, axis=2) &
np.any(pop_fitness[:, np.newaxis] < pop_fitness, axis=2)
)
np.fill_diagonal(dominate_matrix, False)
dominated_count = np.sum(dominate_matrix, axis=0)
# 计算Pareto等级
current_rank = 1
while np.any(ranks == 0):
mask = (dominated_count == 0) & (ranks == 0)
if not np.any(mask):
break
ranks[mask] = current_rank
# 向量化更新被支配计数(优化点)
dominated_count[np.any(dominate_matrix[mask], axis=0)] -= 1
current_rank += 1
return ranks
# 计算新解的Pareto等级
ranks = pareto_ranking(new_FitX)
max_rank = np.max(ranks) if np.max(ranks) > 0 else 1
# 平衡触发概率
progress = it / max_iter
progress_factor = 1 - progress ** 1.5
if max_rank > 1:
normalized_ranks = (ranks - 1) / (max_rank - 1)
else:
normalized_ranks = np.zeros_like(ranks)
trigger_prob = 0.4 + (normalized_ranks * 0.3)
trigger_prob = trigger_prob * progress_factor
trigger_prob = np.clip(trigger_prob, 0.1, 0.9)
# 多样性引导解选择
crowd_dists = crowding_distance(FitA)
sorted_indices = np.argsort(-crowd_dists)
guide_count = max(3, int(10 - 7 * progress))
top_guides = A[sorted_indices[:guide_count]]
# 交叉穿刺
for i in range(N):
if np.random.rand() < trigger_prob[i]:
idx1, idx2 = np.random.choice(len(top_guides), 2, replace=False)
guide1, guide2 = top_guides[idx1], top_guides[idx2]
for j in range(dim):
if np.random.rand() < 0.5:
new_rimepop[i, j] = guide1[j]
else:
new_rimepop[i, j] = guide2[j]
# 重新评估穿刺后的解(批量计算)
new_FitX = objective_func(new_rimepop)
return new_rimepop, new_FitX
# 2. 替换低效的双重循环(mask生成优化,核心优化点2)
def update_population(combined_pop, combined_fit, N, gen, Maxgen, Xmin, Xmax):
"""种群更新:筛选非支配解+对被支配解执行云变异后补充"""
# 1. 筛选非支配解
non_dominated_pop, non_dominated_fit = Domination(combined_pop, combined_fit)
nd_size = len(non_dominated_pop)
if nd_size >= N:
# 非支配解充足时,按拥挤距离选择
crowd_nd = crowding_distance(non_dominated_fit)
sorted_nd_idx = np.argsort(-crowd_nd)
return non_dominated_pop[sorted_nd_idx[:N]]
else:
# 2. 提取被支配解(需补充的解)
remaining = N - nd_size
mask = np.zeros(len(combined_pop), dtype=bool)
for i, pop in enumerate(combined_pop):
for nd_pop in non_dominated_pop:
if np.allclose(pop, nd_pop, atol=1e-6):
mask[i] = True
break
other_pop = combined_pop[~mask] # 被支配解
other_fit = combined_fit[~mask]
# 3. 对被支配解按拥挤距离排序
if len(other_pop) == 0:
# 极端情况:重复选择非支配解
selected = np.random.choice(nd_size, remaining, replace=True)
supplementary = non_dominated_pop[selected]
else:
crowd_other = crowding_distance(other_fit)
sorted_other_idx = np.argsort(-crowd_other)
selected_other = other_pop[sorted_other_idx[:remaining]]
# 4. 对选中的被支配解执行云变异
mutated_supplementary = []
for idx in range(len(selected_other)):
mutated = cloud_mutation(
X=selected_other,
n=idx,
gen=gen,
Maxgen=Maxgen,
Xmin=Xmin,
Xmax=Xmax,
is_non_dominated=False
)
mutated_supplementary.append(mutated)
supplementary = np.array(mutated_supplementary)
# 5. 组合非支配解和变异后的补充解
return np.vstack([non_dominated_pop, supplementary])
def cloud_mutation(X, n, gen, Maxgen, Xmin, Xmax, is_non_dominated):
"""优化版云变异"""
D = X.shape[1]
progress = gen / Maxgen
# 差异化变异概率
if is_non_dominated:
pm = 0.3 - 0.2 * progress
else:
pm = 0.8 - 0.4 * progress
pm = np.clip(pm, 0.05, 0.95)
if np.random.rand() >= pm:
return X[n, :].copy()
# 差异化变异维度数量
if is_non_dominated:
dim_ratio = 0.2 - 0.1 * progress
dim_count = max(1, int(dim_ratio * D))
else:
dim_ratio = 0.6 - 0.3 * progress
dim_count = max(1, int(dim_ratio * D))
mutate_dims = np.random.choice(D, dim_count, replace=False)
# 动态调整En
y = X[n, :].copy()
for i in mutate_dims:
Ex = X[n, i]
dim_std = np.std(X[:, i])
norm_std = dim_std / ((Xmax[i] - Xmin[i]) / 2 + 1e-10)
base_En = (Xmax[i] - Xmin[i]) / (10 * np.sqrt(gen + 1))
En = base_En * (1.2 - 0.8 * norm_std)
En = max(En, 1e-10)
He = En / 10
Enn = np.random.normal(En, He)
Enn = max(abs(Enn), 1e-10)
CloudDrp_pos = np.random.normal(Ex, abs(Enn))
CloudDrp_certainty = np.exp(-(CloudDrp_pos - Ex) ** 2 / (2 * Enn ** 2))
adjusted_pos = Ex + (CloudDrp_pos - Ex) * CloudDrp_certainty
y[i] = np.clip(adjusted_pos, Xmin[i], Xmax[i])
return y
def MORIME_zdt(N, max_iter, lb, ub, dim, objective_func, archive_size, NumObj=2):
"""多目标霜冰优化算法(适配ZDT测试函数)"""
start_time = time.time()
# 初始化种群
rimepop = initialization(N, dim, lb, ub)
# 计算初始种群的目标函数值(批量计算)
FitX = objective_func(rimepop)
# 初始化存档
A, FitA = Domination(rimepop, FitX)
# 防御性检查
if A.size == 0:
rimepop = initialization(N, dim, lb, ub)
FitX = objective_func(rimepop)
A, FitA = Domination(rimepop, FitX)
# 主循环
for it in range(max_iter):
if it % 10 == 0:
print(f"迭代次数: {it + 1}/{max_iter}")
# 阶段1:软霜搜索
new_rimepop = soft_rime_search(rimepop, A, FitA, it, max_iter, lb, ub, dim)
# 阶段2:硬霜穿刺
new_rimepop, new_FitX = hard_rime_puncture(new_rimepop, A, FitA, objective_func, it, max_iter)
# 更新存档
A, FitA = UpdateArchive(NumObj, archive_size, A, FitA, new_rimepop, new_FitX)
# 更新种群
combined_pop = np.vstack([new_rimepop, A])
combined_fit = np.vstack([new_FitX, FitA])
# 调用修改后的种群更新函数
rimepop = update_population(
combined_pop,
combined_fit,
N,
gen=it,
Maxgen=max_iter,
Xmin=lb,
Xmax=ub
)
# 重新计算种群目标值(批量计算)
FitX = objective_func(rimepop)
run_time = time.time() - start_time
return A, FitA, run_time
def plot_results(FitA, true_front, name, igd_value, plot_hull=False):
"""绘制ZDT函数的Pareto前沿对比图"""
plt.figure(figsize=(10, 6))
# 绘制真实帕累托前沿
plt.plot(true_front[:, 0], true_front[:, 1], 'b-', lw=2, label='真实帕累托前沿')
# 绘制算法得到的非支配解
plt.scatter(FitA[:, 0], FitA[:, 1], c='red', s=50, edgecolors='k', alpha=0.7, label='算法求得的非支配解')
# 可选绘制凸包
if plot_hull and len(FitA) > 2:
try:
hull = ConvexHull(FitA)
for simplex in hull.simplices:
plt.plot(FitA[simplex, 0], FitA[simplex, 1], 'g--', lw=1, alpha=0.5)
except:
pass
# 图表设置(显示IGD值)
plt.xlabel('目标函数 f1', fontsize=12)
plt.ylabel('目标函数 f2', fontsize=12)
plt.title(f'{name}函数的Pareto最优前沿对比(IGD={igd_value:.6f})', fontsize=14)
plt.grid(alpha=0.3)
plt.legend()
plt.show()
# 主函数:包含ZDT4边界修正(核心优化点3)
def main(N=100, max_iter=100, archive_size=200):
# 定义要测试的函数,包含ZDT系列和新增的DTLZ系列
functions = [
# ZDT系列函数 (函数, 名称, 决策变量维度)
(zdt1, "ZDT1", 10),
(zdt2, "ZDT2", 10),
(zdt3, "ZDT3", 10),
(zdt4, "ZDT4", 10),
(zdt6, "ZDT6", 10),
# 新增DTLZ系列函数 (函数, 名称, 决策变量维度)
(dtlz1, "DTLZ1", 3),
(dtlz2, "DTLZ2", 3),
(dtlz3, "DTLZ3", 3)
]
# 存储结果
results = {}
# 循环运行每个测试函数
for func, name, dim in functions:
print(f"\n======================================")
print(f"========== {name} 优化开始 ==========")
print(f"======================================\n")
# 设置各函数的边界
if name == "ZDT4":
# ZDT4原始定义:x1∈[0,1], x2~xd∈[-5,5]
lb = np.zeros(dim)
ub = np.ones(dim)
lb[1:] = -5.0 # 第2到第d维下界为-5
ub[1:] = 5.0 # 第2到第d维上界为5
elif name.startswith("ZDT"):
# 其他ZDT函数:所有维度∈[0,1]
lb = np.zeros(dim)
ub = np.ones(dim)
else: # DTLZ系列
# DTLZ1/2/3所有维度∈[0,1]
lb = np.zeros(dim)
ub = np.ones(dim)
# 运行优化算法
A, FitA, run_time = MORIME_zdt(
N, max_iter, lb, ub, dim, func, archive_size
)
# 生成真实帕累托前沿
true_front = generate_true_front(func)
# 计算IGD值
igd_value = calculate_igd(FitA, true_front)
# 保存结果
results[name] = {
'solutions': A,
'fitness': FitA,
'run_time': run_time,
'igd': igd_value,
'true_front': true_front
}
# 打印结果统计
print(f"\n========== {name} 优化结果 ==========")
print(f"算法运行时间: {run_time:.2f}秒")
print(f"找到的非支配解数量: {len(FitA)}")
print(f"IGD值: {igd_value:.6f}")
# 绘制结果
plot_results(FitA, true_front, name, igd_value, plot_hull=False)
# 打印所有IGD值对比
print("\n======================================")
print("各函数的IGD值对比:")
for func_name, res in results.items():
print(f"{func_name}: {res['igd']:.6f}")
print("======================================")
return results
if __name__ == "__main__":
results = main()
请你检查我的DTLZ函数是否编写错误,为什么这三个函数的图像都很差