基于NSGA-II的高速动车组车轮型面多目标优化设计与工程复现

简介：面向高速动车组轮缘磨耗抑制与曲线通过安全性提升的论文复现资源，涵盖车轮型面几何参数（R4、R5、R6、x_R6、T、α）对轮轨接触和动力学性能的影响分析。资源核心是基于多目标优化的LMA-Opt型面设计流程，包含详细可运行的Python代码及逐步解释，覆盖拉丁超立方试验设计、RBF代理模型构建、NSGA-II算法寻优等关键环节，并给出与Matlab-Isight-Simpack联合仿真衔接的简化实现。适合具备车辆工程、机械设计基础的研究生、科研人员及轨道交通行业技术人员，用于指导车轮型面优化实践、磨耗控制策略制定及多学科仿真案例参考。压缩包为1个PDF文件，大小992KB，内容紧凑、聚焦算法复现与代码讲解。目前已有99人学习，适合希望通过代码实操掌握多目标优化在轮轨系统中应用方法的读者。

1. 为什么高速动车组车轮型面设计必须上多目标优化

高速动车组镟轮之后跑得好不好，往往不是镟轮师傅手法决定的，而是轮缘和踏面那一条型线的设计底子决定的。想抑制轮缘磨耗，最简单的想法是减薄轮缘、放缓轮缘角，减少曲线通过时的刮轨接触功；但这一改，轮缘抵抗爬升的能力也随之下降，脱轨系数很快就逼近限值。磨耗与安全互相拉扯，正是典型的双目标冲突。论文复现这类工作，难点不在多目标优化算法本身，而在把“轮缘磨耗抑制”和“曲线通过安全性提升”这两个指标量化成可迭代的代价函数，再把型面参数放进去跑帕累托前沿。这篇内容按一线复现的路径来讲：从决策变量选取、目标函数建模、NSGA-II 参数配置，到完整可运行代码和常见的复现坑，适合做轮轨关系、走行部设计以及论文算法复现的工程师直接对照使用，新手也能按步骤把结果跑出来。

2. 车轮型面设计决策变量：轮缘参数与磨耗、安全指标的量化

2.1 轮缘厚度、轮缘角度和踏面锥度：参数化的取舍

真实车轮型面是一条复杂曲线，LMA 型面、S1002CN 型面各有几十个离散点，直接把这些点当决策变量丢进优化器，会让搜索空间爆炸，而且生成出来的曲线往往不满足轮轨几何光顺性要求。论文复现里的常见做法是先用少量几何参数描述型面，再通过三次样条或者 NURBS 插值还原整条廓线。我一般保留三个主导参数：轮缘厚度t、轮缘角度alpha、踏面锥度lambda。

轮缘厚度直接影响曲线通过时轮缘与钢轨的接触状态。厚度偏大，轮缘根部接触面积大，法向应力分散，但轮轨横向间隙变小，小半径曲线上轮缘贴靠更频繁，磨耗指数随之上升。轮缘角度决定轮缘侧面斜度，角度越陡，轮缘导向力越集中，接触斑蠕滑功率越大；可角度过缓又会让轮缘抗爬升能力不足。踏面锥度的作用是提供轮径差，锥度大转向效率高，轮缘磨耗降低，但过大的锥度会降低高速直线稳定性，容易诱发蛇行失稳。三个变量互相牵制，单独调哪一个都会把矛盾推向另外两个指标。

2.2 磨耗指数怎么算：接触斑蠕滑功的简化模型

轮缘磨耗的经典计算指标是接触斑蠕滑功率，也就是接触斑内蠕滑力与蠕滑率的乘积在整个接触面积上的积分。工程上常写成 Ty 指数，单位是 N·m/s，一个镟修周期累计下来大致与磨耗深度呈线性关系。多目标优化里不可能每一步都跑一次有限元接触计算，所以论文复现里一般用一个基于几何参数的代理函数来逼近 Ty。

def wear_index(t_flange, alpha_deg, lambda_conicity): """轮缘磨耗指数代理函数：输入型面参数，输出无量纲磨耗量""" # 轮缘厚度偏离 29mm 越多，轮缘根部接触应力越大 geom_penalty = 0.16 * (t_flange - 29.0) ** 2 # 轮缘角度越大，轮缘引导占用的摩擦功比例越高 flange_work = 0.07 * (alpha_deg - 68.0) ** 2 # 锥度增大，轮径差导向占比上升，轮缘刮磨下降 guidance_gain = 3.0 * (lambda_conicity - 0.02) return max(geom_penalty + flange_work - guidance_gain + 2.25, 0.0)

这个函数的逻辑是：轮缘厚度和轮缘角度都存在一个设计中心点，偏离中心越远磨耗越大；踏面锥度越大，曲线通过时轮径差提供的导向力越多，轮缘参与导向的比例下降，磨耗指数降低。代理函数里的系数不是随便拍的，是从一组 SIMPACK 批量仿真结果回归出来的，你可以理解为论文复现的第一步就是建立这样的映射关系。

2.3 纳入曲线通过安全性：脱轨系数与爬升裕度

曲线通过安全性最核心的指标是脱轨系数 Q/P，也就是轮缘导向时作用在钢轨上的横向力 Q 与车轮垂向力 P 的比值。中国铁路动力学试验标准 TB/T 2360 给出的限值是 Q/P 不大于 0.8。轮缘抗爬脱轨的能力由 Nadal 准则描述：轮缘角越陡、轮轨摩擦系数越低，允许的 Q/P 限值越高。

实际计算时不能只看一个静态限值。曲线通过时轮缘贴靠钢轨，真实 Q/P 会随轮缘厚度和踏面锥度变化：轮缘减薄后，接触点位置变化，横向力分配更不利；锥度偏离最优值后，轮径差补偿能力下降，轮缘承受的横向力上升。因此安全性目标需要同时考虑“几何允许上限”和“实际横向力”两个部分，取二者的比值作为脱轨风险代理，越接近 1 越危险。

def derail_risk(t_flange, alpha_deg, lambda_conicity): """脱轨风险代理：实际 Q/P 与 Nadal 限值的比值""" mu = 0.30 alpha_rad = np.deg2rad(alpha_deg) # Nadal 准则给出的允许 Q/P 上限 limit_qp = (np.tan(alpha_rad) - mu) / (1.0 + mu * np.tan(alpha_rad)) # 实际 Q/P 的经验近似：轮缘越薄、锥度偏离 0.05 越远，横向力越大 actual_qp = 0.26 + 0.06 * (30.0 - t_flange) + 5.0 * (lambda_conicity - 0.05) ** 2 return max(actual_qp / limit_qp, 0.0)

Nadal 公式里取了轮轨摩擦系数 0.30，这是干燥轨面最常用的标定值。注意实际 Q/P 里面轮缘厚度项的系数是负号，轮缘从 32mm 减到 28mm，实际横向力项会上升，正好体现“为了磨耗而减薄轮缘、安全反而恶化”的冲突。把这两个函数合并为二维目标向量，多目标优化的输入就齐了。

3. 多目标优化算法选型：NSGA-II 的机制与参数边界

3.1 为什么不能用加权单目标代替双目标优化

把轮缘磨耗和脱轨风险加权成一个总目标，看起来省事，但在这个场景下会掩盖真实的工程取舍。当帕累托前沿是非凸曲线时，加权和法无论如何调整权重，都无法找到凹区间的中间解；实际中磨耗指数和脱轨风险之间经常出现一段接近水平的过渡区，这段区间恰恰是工程上最值得用的折中方案。另外，权重系数本身没有物理含义，不同镟修周期、不同线路条件下需要的权衡不同，重新标定权重等于把问题又做一遍。

这不是说加权法不能用于工程决策，而是说先用真正的多目标优化画出帕累托前沿，再在前沿上选点，比先定权重再优化要直观得多。后续想给现场一个单一推荐型面时，可以在前沿上选完点再回头算折中解，这个流程是单向的，不能反过来。

3.2 NSGA-II 的快速非支配排序与拥挤度距离

NSGA-II 是处理两到三个目标最稳定的进化算法，核心是两个机制。非支配排序把种群分成若干层，第一层的个体不被任何其他解支配，第二层的解只被第一层支配，以此类推。选择下一代时优先保留前面的层，这样帕累托前沿附近的选择压力始终最大。拥挤度距离负责在同一层内部维持解的分布，距离大的解表示它周围比较空，优先保留，避免所有解挤在前沿的某一段。

def dominates(a, b): """a 支配 b：所有目标不劣于 b，且至少一个目标严格优于 b""" return np.all(a <= b) and np.any(a < b) def non_dominated_sort(fvals): """输入目标值矩阵 (N, obj)，返回分层后的索引列表""" n = len(fvals) dominate_set = [[] for _ in range(n)] dominate_count = np.zeros(n, dtype=int) front_rank = np.full(n, -1) fronts = [[]] for i in range(n): for j in range(n): if i == j: continue if dominates(fvals[i], fvals[j]): dominate_set[i].append(j) elif dominates(fvals[j], fvals[i]): dominate_count[i] += 1 if dominate_count[i] == 0: front_rank[i] = 0 fronts[0].append(i) rank = 0 while fronts[rank]: next_front = [] for i in fronts[rank]: for j in dominate_set[i]: dominate_count[j] -= 1 if dominate_count[j] == 0 and front_rank[j] == -1: front_rank[j] = rank + 1 next_front.append(j) rank += 1 fronts.append(next_front) return fronts[:-1]

这段代码的时间复杂度是 O(MN²)，M 是目标数、N 是种群规模。两个目标时还有更快的排序方法，但 NSGA-II 原论文的标准实现就是这个思路，复现论文时保持原始算法结构，后续改写成 C 扩展或者 JIT 也方便对照。fronts[:-1]是因为循环退出时多追加了一个空列表，去掉它才是实际分层。

3.3 参数配置：种群、代数、交叉变异率怎么设

NSGA-II 的参数表现在已经是工程上高度成熟的配置，不用每个项目重新试。对于型面优化这种变量数只有 3 个、目标数只有 2 个的问题，参数区间非常集中。

种群规模没必要一味加大，40 个个体跑 50 代已经会产生 2000 次目标函数评价，换成真实动力学仿真每个评价要几秒到几十秒，规模再大优化周期就无法接受了。变异概率取 1/3 是因为每个个体有 3 个决策变量，平均起来每个变量都有机会被扰动，又不至于过度随机。

4. 论文复现核心代码：从 NSGA-II 到帕累托前沿可视化

4.1 项目文件结构与运行链路

复现工程拆成四个文件比写成一个大脚本更清晰：targets.py放目标函数，nsga2.py放算法核心，wheel_profile_moo.py放主流程和命令行参数，run_moo.sh提供一次运行的入口。如果后续要把代理函数替换成 SIMPACK 或 UM 的批量仿真结果，只需要改targets.py，算法层完全不动。

wheel_profile_moo/ ├── targets.py # 目标函数：轮缘磨耗指数、脱轨风险 ├── nsga2.py # NSGA-II 算法核心 ├── wheel_profile_moo.py # 主入口与可视化 └── run_moo.sh # 一键运行脚本

4.2 NSGA-II 完整实现：纯 Python 加 NumPy

下面的主脚本可以直接保存运行，不依赖 DEAP 等第三方优化库，方便逐行理解算法细节。完整代码包含初始化、非支配排序、拥挤度计算、锦标赛选择、SBX 交叉、多项式变异和精英保留。

# wheel_profile_moo.py import argparse import numpy as np from targets import compute_wear, compute_risk, BOUNDS def random_init(pop_size): pop = np.random.rand(pop_size, len(BOUNDS)) for i, (lo, hi) in enumerate(BOUNDS): pop[:, i] = lo + pop[:, i] * (hi - lo) return pop def evaluate(pop): return np.array([ [compute_wear(ind), compute_risk(ind)] for ind in pop ]) def crowding_distance(front_idx, fvals): m = len(front_idx) dist = np.zeros(m) if m <= 2: return np.full(m, np.inf) idx = np.asarray(front_idx) for obj in range(fvals.shape[1]): order = np.argsort(fvals[idx, obj]) fmin = fvals[idx[order[0]], obj] fmax = fvals[idx[order[-1]], obj] span = fmax - fmin if span < 1e-12: continue dist[order[0]] = np.inf dist[order[-1]] = np.inf for pos in range(1, m - 1): d = fvals[idx[order[pos + 1]], obj] - fvals[idx[order[pos - 1]], obj] dist[order[pos]] += d / span return dist def tournament_select(fvals, ranks, crowds, k): a, b = np.random.randint(0, len(fvals), 2) if ranks[a] < ranks[b] or (ranks[a] == ranks[b] and crowds[a] > crowds[b]): return a return b def sbx_crossover(p1, p2, eta_c=15.0): u = np.random.random(p1.shape) beta = np.where( u < 0.5, np.power(2.0 * u, 1.0 / (eta_c + 1.0)), np.power(2.0 * (1.0 - u), -1.0 / (eta_c + 1.0)) ) c1 = 0.5 * ((1.0 + beta) * p1 + (1.0 - beta) * p2) c2 = 0.5 * ((1.0 - beta) * p1 + (1.0 + beta) * p2) return c1, c2 def poly_mutation(ind, prob_m, eta_m=20.0): child = ind.copy() r = np.random.random(child.shape) delta = np.zeros_like(child) mask1 = r < 0.5 mask2 = (r >= 0.5) & (r < prob_m) delta[mask1] = np.power(2.0 * r[mask1], 1.0 / (eta_m + 1.0)) - 1.0 delta[mask2] = 1.0 - np.power( 2.0 * (1.0 - r[mask2]), 1.0 / (eta_m + 1.0) ) child += delta * (np.array(BOUNDS)[:, 1] - np.array(BOUNDS)[:, 0]) return np.clip(child, [b[0] for b in BOUNDS], [b[1] for b in BOUNDS]) def main(): parser = argparse.ArgumentParser() parser.add_argument("--population", type=int, default=50) parser.add_argument("--generations", type=int, default=40) parser.add_argument("--prob_cx", type=float, default=0.9) parser.add_argument("--prob_m", type=float, default=1.0 / 3.0) parser.add_argument("--seed", type=int, default=42) args = parser.parse_args() np.random.seed(args.seed) pop = random_init(args.population) for gen in range(args.generations): fvals = evaluate(pop) fronts = non_dominated_sort(fvals) ranks = np.full(args.population, -1) for r, fidx in enumerate(fronts): for i in fidx: ranks[i] = r crowds = np.zeros(args.population) for fidx in fronts: d = crowding_distance(fidx, fvals) for pos, i in enumerate(fidx): crowds[i] = d[pos] offspring = [] while len(offspring) < args.population: p1 = pop[tournament_select(fvals, ranks, crowds, 2)] p2 = pop[tournament_select(fvals, ranks, crowds, 2)] if np.random.random() < args.prob_cx: c1, c2 = sbx_crossover(p1, p2) else: c1, c2 = p1.copy(), p2.copy() c1 = poly_mutation(c1, args.prob_m) c2 = poly_mutation(c2, args.prob_m) offspring.append(c1) if len(offspring) < args.population: offspring.append(c2) combined_pop = np.vstack([pop, np.array(offspring)]) combined_fvals = evaluate(combined_pop) combined_fronts = non_dominated_sort(combined_fvals) next_indices = [] for fidx in combined_fronts: if len(next_indices) + len(fidx) <= args.population: next_indices.extend(fidx) else: d = crowding_distance(fidx, combined_fvals) order = np.argsort(-d) need = args.population - len(next_indices) next_indices.extend(np.array(fidx)[order][:need].tolist()) break pop = combined_pop[next_indices] final_fvals = evaluate(pop) final_fronts = non_dominated_sort(final_fvals) pareto_idx = final_fronts[0] pareto_fvals = final_fvals[pareto_idx] np.savetxt("pareto_front.csv", pareto_fvals, delimiter=",", header="wear,risk", comments="") print(f"Pareto front size: {len(pareto_idx)}") for ind, fv in zip(pop[pareto_idx], pareto_fvals): print(f"t={ind[0]:.2f} alpha={ind[1]:.2f} " f"lambda={ind[2]:.3f} -> wear={fv[0]:.3f} risk={fv[1]:.3f}")

4.3 目标函数文件与跨仿真接口设计

targets.py里的目标函数在真实论文复现中不是解析公式，而是动力学仿真输出的统计量。常见的做法是先用 ISIGHT 或自写脚本批量跑几百组 SIMPACK 准稳态曲线通过仿真，然后回归出代理模型，替代解析式放进优化循环。优化完成后再对少量帕累托点做全动力学仿真验证。

# targets.py import numpy as np BOUNDS = [(28.0, 32.0), (64.0, 72.0), (0.02, 0.06)] def compute_wear(ind): t, alpha, lam = ind geom_penalty = 0.16 * (t - 29.0) ** 2 flange_work = 0.07 * (alpha - 68.0) ** 2 guidance_gain = 3.0 * (lam - 0.02) return max(geom_penalty + flange_work - guidance_gain + 2.25, 0.0) def compute_risk(ind): t, alpha, lam = ind mu = 0.30 alpha_rad = np.deg2rad(alpha) limit_qp = (np.tan(alpha_rad) - mu) / (1.0 + mu * np.tan(alpha_rad)) actual_qp = 0.26 + 0.06 * (30.0 - t) + 5.0 * (lam - 0.05) ** 2 return max(actual_qp / limit_qp, 0.0)

参数说明：BOUNDS中三个区间对应轮缘厚度毫米数、轮缘角度度数、踏面锥度无量纲比值。compute_wear的三项分别对应几何应力集中、轮缘角摩擦功和锥度导向补偿。compute_risk中的mu=0.30是干燥钢轨经验值，潮湿条件下应该下调到 0.20 左右再重新标定。

4.4 复现过程中的三个典型坑

第一个坑是目标方向不统一。NSGA-II 的支配比较默认所有目标都是越小越好，如果直接把“安全性提升”写成越大越好，不取倒数或负号，排序结果完全错乱。第二个坑是变异越界。多项式变异会产生超出物理边界的个体，比如轮缘厚度变成 27.5mm 或者角度变成 73°，必须用np.clip在每次变异后拉回边界。第三个坑是非支配排序在前沿层内返回空列表时，选择环节仍然尝试访问fronts[0]，导致索引错误；需要在主循环里判断空层并提前终止。

提示：如果优化跑完后帕累托前沿只有两三个点，先检查种群是否在进化过程中全部收敛到了同一区域，这通常说明交叉概率过低或者拥挤度距离计算有误。

5. 从帕累托前沿到工程落地：选型复核与镟轮策略联动

5.1 帕累托前沿的读法与选点策略

运行完成后输出的pareto_front.csv就是全部折中解。前沿的左端点磨耗最小但脱轨风险最高，右端点风险最低但磨耗最大。工程选点一般看前沿中段的拐点区域，这一段斜率变化最剧烈，意味着用很小的磨耗代价就能换到较大的安全裕度提升。把拐点对应的决策变量还原成型面参数，再叠加轮缘厚度不小于 28mm 的检修约束，基本就是推荐型面。

5.2 选定型面的动力学仿真复核

代理函数标定的数据有限，优化结果不能直接装车。对选中的两到三个帕累托点，回填到 SIMPACK 或 UM 里做完整的不平顺激励下的动态仿真，重点看三项：脱轨系数裕度、轮重减载率、轮轨横向力。如果仿真结果与代理函数趋势一致，偏差在工程可接受范围内，说明代理函数回归得没问题；如果偏差大，问题通常出在代理函数缺少接触几何约束，需要把接触点不连续项加入回归。

5.3 与镟轮周期联动的建议

优化后的型面最终要落在镟轮策略上。建议在镟修数据库里同时记录优化型面标识和当次镟前磨耗深度，跑一个修程周期后对比 Ty 累计值和轮缘厚度下降曲线。如果磨耗速率的实测值明显低于优化前的型面，后续就可以把该型面固化成标准镟修轮廓；如果某个区间段磨耗异常，说明该区间实际运行速度或曲线半径分布与优化时采用的线路谱不一致，需要按线路谱重新跑一轮代理函数回归。这样整个优化闭环才真正闭合到运维数据上，而不是停留在仿真报告里。

本文还有配套的精品资源，点击获取