Python实现周期性两组三维坐标的最小距离最优匹配

循环球体堆叠仿真的周期衔接平滑问题本质是带权指派问题，有成熟可直接落地的解决方案，无需自研算法。
对应仿真逻辑为：单周期内完成N步排料仿真，球体持续从底部移出、顶部随机插入，周期末状态需要和下一周期初的状态做球体编号重排，保证衔接处同编号球体坐标差最小，且底部球体匹配精度优先级高于顶部球体。

这类问题属于加权二分图最小权完美匹配（即带权指派问题），可直接用匈牙利算法求解，针对百级球体规模计算耗时在毫秒级，完全满足使用要求：

• 将周期末（t=Nsteps）的所有球体作为二分图左侧节点，周期初（t=0）的所有球体作为二分图右侧节点
• 左右节点间的边权为「两球三维欧氏距离 × 对应位置的精度权重」
• 求解总边权最小的完美匹配，得到的映射关系就是需要的初始坐标索引重排规则
• 精度权重按z轴位置设置：z值越小（越靠近底部）权重越大，z值越大（越靠近顶部）权重越小，可自然实现「底部匹配优先、顶部精度放宽」的要求，不需要额外添加复杂约束。

直接使用scipy库内置的linear_sum_assignment函数即可，这是经过工业界验证的优化版匈牙利算法实现，不需要安装额外复杂依赖。
完整可运行代码如下：

import numpy as np
from scipy.optimize import linear_sum_assignment

def adjust_initial_positions(initial_positions: np.ndarray, 
                             final_positions: np.ndarray,
                             z_axis_index: int = 2,
                             z_min: float = 0,
                             z_max: float = 100,
                             bottom_priority: float = 10) -> np.ndarray:
    """
    求解初始位置的最优重排索引，保证周期衔接平滑
    参数说明：
    initial_positions: 周期初始位置数组，shape=(Nspheres, 3)
    final_positions: 周期结束位置数组，shape=(Nspheres, 3)
    z_axis_index: 竖直向上的轴对应的索引，默认2对应z轴
    z_min: 箱体底部z坐标
    z_max: 箱体顶部z坐标
    bottom_priority: 底部球体的匹配优先级权重，值越大底部匹配精度越高，顶部权重固定为1
    """
    n_spheres = initial_positions.shape[0]
    # 计算两两球体之间的欧氏距离矩阵，shape=(n_spheres, n_spheres)
    # dist[i,j] = 周期末第i个球和周期初第j个球的三维距离
    dist = np.linalg.norm(
        final_positions[:, np.newaxis, :] - initial_positions[np.newaxis, :, :],
        axis=2
    )
    # 计算每个周期末位置球体的精度权重：越靠底部权重越大
    z_final = final_positions[:, z_axis_index]
    # z坐标归一化到0-1区间，0对应顶部，1对应底部
    z_norm = np.clip((z_max - z_final) / (z_max - z_min), 0, 1)
    weight = 1 + (bottom_priority - 1) * z_norm
    # 构造带权代价矩阵
    cost_matrix = dist * weight[:, np.newaxis]
    # 匈牙利算法求解最小权完美匹配
    row_ind, col_ind = linear_sum_assignment(cost_matrix)
    # 返回重排索引：col_ind[i]表示周期末第i个球匹配周期初第col_ind[i]个球
    return col_ind

# 测试逻辑
if __name__ == "__main__":
    Nspheres = 100
    Nsteps = 5
    coordinates = np.random.uniform(0,100, (Nsteps, Nspheres, 3))
    initial_positions = coordinates[0]
    final_positions = coordinates[Nsteps-1]
    
    indices_adjust_initial_positions = adjust_initial_positions(initial_positions, final_positions)
    adjusted_initial_positions = initial_positions[indices_adjust_initial_positions] 

    mean_error = np.mean(np.abs(final_positions-adjusted_initial_positions))
    max_error = np.max(np.abs(final_positions-adjusted_initial_positions))
    print(f"平均匹配误差: {mean_error:.3f}, 最大匹配误差: {max_error:.3f}")

    Ncycles = 5
    simulation_coordinates = np.empty((Nsteps*Ncycles, Nspheres, 3))
    simulation_coordinates[:Nsteps] = np.array(coordinates)
    for n in range(1, Ncycles):
        new_cycle_coordinates = simulation_coordinates[Nsteps*(n-1):Nsteps*(n), indices_adjust_initial_positions, :]
        simulation_coordinates[Nsteps*n:Nsteps*(n+1)] = new_cycle_coordinates
    print(f"循环仿真序列维度: {simulation_coordinates.shape}")

• 如果需要进一步强化底部匹配精度，只要调大bottom_priority参数即可，比如设置为100时，求解器会优先保证底部球体的匹配距离最小，顶部球体即使匹配距离稍大也不会显著影响总代价
• 如果需要完全忽略某一高度以上的球体匹配精度，可以直接将对应高度以上的权重设为固定极小值（比如0.01）
• 该方法得到的是全局最优解，不存在局部最优问题，针对千级以内球体规模的计算速度都能满足实时仿真要求。

内容的提问来源于stack exchange，提问作者yvrob