《Python机器人基础原理》第十三篇：焊接产线仿真——计算雅可比矩阵行列式，依据数值判断当前构型是否靠近奇异位形

 

1. 实际工业应用场景

 

周五上午，某新能源汽车电池托盘焊接产线。一台4-DOF关节型工业机械臂搭载MIG焊枪，对托盘铝合金边框进行连续角焊缝作业。机器人构型：3旋转+1垂直平移（SCARA变体），大臂长500 mm，小臂长400 mm，基座安装高度d₁=350 mm，z轴行程±150 mm。关节硬件限位：θ₁ ±160°、θ₂ ±120°、θ₃ ±180°、d₄ ±150 mm。额定负载8 kg，焊枪及线缆总质量2.3 kg。焊接工艺要求末端线速度恒定在60~80 mm/s，电弧跟踪精度±0.3 mm。

 

上午十点一刻，焊接质量巡检员小刘拿着抽检报告找到控制工程师小王。

 

"第47号到第52号焊缝，边缘有咬边和焊瘤。"小刘指着X光片说，"全都在托盘拐角处，也就是机械臂手臂几乎伸直的位置。"

 

小王调出该工位的轨迹日志，皱眉看了半天："轨迹规划给的末端速度曲线是恒定的，60 mm/s，没问题啊。"

 

小赵从后面走过来，看了眼屏幕："你用的什么方法保证速度恒定？"

 

"逆雅可比算关节速度，然后直接下发给伺服。"

 

"奇异位形附近呢？"

 

小王愣了一下。

 

"你的逆雅可比在接近奇异的时候，关节速度会趋向无穷大。"小赵在纸上画了个构型，"大臂小臂几乎伸直，这时候雅可比矩阵的行列式趋近于零——你有没有做奇异检测？"

 

"没有。我以为轨迹都在工作空间里面就行。"

 

"工作空间内部也有奇异线。你需要实时算雅可比行列式，当行列式绝对值掉到某个阈值以下，就判定当前构型靠近奇异，触发降速或者切换冗余自由度。"

 

本篇解决的问题：给定一组关节角度，计算该位形下雅可比矩阵的行列式值，依据行列式数值是否低于预设阈值，判断当前构型是否靠近奇异位形。对于SCARA构型，末端位置只由θ₁和θ₂两个关节决定（腕部roll关节θ₃只改变姿态、第4关节只改变高度），因此位置奇异的本质是平面2连杆的奇异，用2×2子块行列式 |det| = L₁·L₂·|sin(θ₂)| 判定。通过在工作空间网格上扫描该指标，生成奇异区域热力图，为轨迹规划阶段的奇异规避提供离线校核依据。

 

⚠️ 声明：本篇全部代码为课程离线仿真，运行在普通PC上，用于雅可比行列式计算与奇异位形检测分析。不可直接用于真实机器人控制器、伺服驱动器或焊接产线。真实工业现场需结合实时动力学补偿、关节扭矩限幅与安全监控。代码作者不对直接使用本代码于真实设备造成的任何后果负责。

 

2. 痛点分析

 

课程模块 本篇对应的工程痛点

第3章 连杆参数及运动学求解 → 3.5 速度雅可比 雅可比矩阵行列式趋零是逆运动学数值不稳定的根源

第3章 → 3.6 奇异位形 工作空间内部存在构型奇异线，需要行列式阈值检测

第6章 轨迹规划与控制 → 6.3 笛卡尔空间规划 逆雅可比求关节速度时，奇异附近关节速度发散导致末端失控

第2章 位姿描述与坐标变换 行列式计算依赖雅可比子块列向量的线性独立性

 

本篇代码解决 L3 离线运动学分析层。输入一组关节角度，计算雅可比2×2子块行列式值，与预设阈值比较，输出奇异判定结果。代码不是设备实时控制器，不做伺服闭环，不做在线奇异规避。它做的事情是：在离线环境中扫描工作空间各区域的行列式分布，告诉轨迹规划工程师"哪些区域是危险的奇异区，规划路径要避开"。

 

不覆盖：逆雅可比求解（伪逆/阻尼最小二乘/冗余解析）、实时奇异规避算法、动力学计算、碰撞检测、6轴冗余机械臂、力控制。

 

3. 核心逻辑讲解

 

物理直觉

 

雅可比矩阵的行列式衡量的是"关节速度空间到末端速度空间的映射体积缩放比"。想象一个三维的小立方体，每个顶点代表一组略微不同的关节速度组合。经过雅可比矩阵这个"变换透镜"之后，这个立方体在末端速度空间里变成了一个平行六面体。行列式就是这个平行六面体的体积。

 

当机械臂处于普通构型时，各关节对末端速度的贡献方向互不重叠，行列式值较大，映射是"健康的"。当机械臂接近奇异构型（比如大臂小臂完全伸直），两个关节对末端速度的贡献方向几乎平行——平行六面体被压扁了，体积趋近于零，行列式趋近于零。这时候逆映射（从末端速度反推关节速度）就变得极度敏感：末端速度稍有波动，算出来的关节速度就会爆掉。

 

SCARA构型的特殊性：为什么用2×2子块

 

SCARA构型的前3个旋转关节轴线全部竖直（z轴方向[0,0,1]），关节3是腕部roll回转，a₃=0（不产生位移），第4关节只沿z轴平移。这意味着：

 

- 末端位置只由θ₁和θ₂决定（平面2连杆问题）

- θ₃只改变末端姿态角，不产生线速度

- d₄只改变末端高度z，不影响xy平面位置

 

因此，SCARA的位置奇异只由θ₁和θ₂决定，与θ₃和d₄无关。雅可比线速度子块的前3列中，第3列恒为零向量（因为z₂×竖直向量=0），3×3行列式恒等于零——用3×3子块对本构型毫无意义，全部构型都会被误判为奇异。

 

正确的做法：取J_v的前2行（xy平面线速度）和前2列（θ₁、θ₂），构成2×2子块：

 

|det(J₂×₂)| = L₁·L₂·|sin(θ₂)|

 

θ₂→0°（手臂完全伸直）或θ₂→±180°（手臂完全折叠）时行列式趋零，对应两种经典奇异位形。θ₂=90°时行列式最大，构型最"健康"。

 

业务逻辑 → 代码映射

 

输入：一组关节角度 q = [θ₁, θ₂, θ₃, d₄]

  ↓

步骤1：正运动学计算各连杆变换矩阵 T₀¹, T₁², T₂³, T₃⁴

  ↓

步骤2：提取大臂末端位置 p₁ 和末端位置 p_end

  ↓

步骤3：构建雅可比2×2线速度子块——

   J₂×₂ = [[∂x/∂θ₁, ∂x/∂θ₂],

            [∂y/∂θ₁, ∂y/∂θ₂]]

         = [[-L₁sinθ₁-L₂sin(θ₁+θ₂), -L₂sin(θ₁+θ₂)],

            [ L₁cosθ₁+L₂cos(θ₁+θ₂), L₂cos(θ₁+θ₂)]]

  ↓

步骤4：计算 det(J₂×₂) = L₁·L₂·sin(θ₂)

  ↓

步骤5：与阈值比较 |det| < threshold → 标记奇异

  ↓

步骤6：在工作空间(r, θ₂)网格上解析生成行列式分布热力图

 

关键公式文字翻译

 

2×2子块行列式：|det| = L₁·L₂·|sin(θ₂)|。L₁和L₂是两个臂长（m），sin(θ₂)是肘角正弦。物理含义：θ₂越接近0°或±180°，两个臂的运动方向越趋于平行，末端位置对关节角的敏感度越高，行列式越小。单位是m²（两个长度量纲相乘），不是m³/s³。

 

阈值标定：对L₁=500 mm、L₂=400 mm的构型，行列式量级在0~0.2 m²之间。取阈值1e-4 m²，对应|sin(θ₂)|<5e-4，即θ₂距奇异位置约0.03°以内才触发——这个阈值偏紧，适合高精度的焊接场景。若想更早预警，可放宽到1e-3 m²（对应θ₂距奇异约0.14°）。

 

4. OOP代码实现

 

技术栈

 

严格限定：numpy, scipy, matplotlib。不引入pybullet、roboticstoolbox等任何第三方机器人库。sympy仅用于符号推导验证，本篇核心计算全部使用数值方法。

 

项目结构树

 

robot_singularity_sim/

├── run_singularity_sim.py

├── requirements.txt

├── singularity/

│ ├── __init__.py

│ ├── models.py

│ ├── forward_kinematics.py

│ ├── jacobian_solver.py

│ ├── singularity_detector.py

│ └── visualize.py

├── tests/

│ ├── __init__.py

│ └── test_singularity_sim.py

└── results/

 

核心源码

 

"singularity/models.py" — dataclass数据模型

 

"""

4-DOF工业机械臂数据模型

连杆参数、D-H表、关节限位、奇异检测阈值结构体

"""

 

from dataclasses import dataclass, field

from typing import List, Tuple

import numpy as np

 

 

@dataclass

class JointLimit:

    """

    单关节硬件限位

    min_val: 最小允许值 (rad 或 m)

    max_val: 最大允许值 (rad 或 m)

    """

    min_val: float

    max_val: float

    joint_name: str = ""

 

    def __post_init__(self):

        if self.min_val >= self.max_val:

            raise ValueError(

                f"JointLimit: min_val({self.min_val})必须严格小于"

                f"max_val({self.max_val})")

        if not np.isfinite(self.min_val) or not np.isfinite(self.max_val):

            raise ValueError("JointLimit: 限位值必须为有限数值")

 

    def contains(self, value: float) -> bool:

        """

        判断关节值是否在限位区间内

        浮点比较使用1e-10容差，避免边界舍入误判

        """

        return (self.min_val - 1e-10) <= value <= (self.max_val + 1e-10)

 

 

@dataclass

class JointLimitSet:

    """

    完整关节限位集合（4-DOF工业机械臂）

    """

    limits: List[JointLimit]

 

    def __post_init__(self):

        if len(self.limits) != 4:

            raise ValueError(f"需要4个关节限位, 实际{len(self.limits)}个")

 

    @classmethod

    def default_4dof_welding(cls):

        """

        默认4-DOF焊接机械臂关节限位

        θ1: ±160° θ2: ±120° θ3: ±180° d4: ±150 mm

        """

        pi = 3.141592653589793

        return cls(limits=[

            JointLimit(-160 * pi / 180, 160 * pi / 180, 'joint1_base_yaw'),

            JointLimit(-120 * pi / 180, 120 * pi / 180, 'joint2_elbow_yaw'),

            JointLimit(-180 * pi / 180, 180 * pi / 180, 'joint3_wrist_roll'),

            JointLimit(-0.150, 0.150, 'joint4_z_lift'),

        ])

 

    def check_all_with_indices(self, joint_values: List[float]) -> Tuple[bool, List[int]]:

        """

        检查全部关节是否在限位内

        返回：(全部通过?, 超限关节索引列表)

        使用索引而非字符串匹配，避免子串匹配脆弱性

        """

        if len(joint_values) != 4:

            raise ValueError(f"需要4个关节值, 实际{len(joint_values)}个")

 

        all_pass = True

        violation_indices = []

 

        for i, (val, limit) in enumerate(zip(joint_values, self.limits)):

            if not limit.contains(val):

                all_pass = False

                violation_indices.append(i)

 

        return all_pass, violation_indices

 

 

@dataclass

class DHParam:

    """标准D-H参数"""

    theta: float

    d: float

    a: float

    alpha: float

 

    def __post_init__(self):

        for name, val in [('theta', self.theta), ('d', self.d),

                           ('a', self.a), ('alpha', self.alpha)]:

            if not np.isfinite(val):

                raise ValueError(f"DHParam.{name}必须为有限数值")

 

 

@dataclass

class SCARAParams:

    """

    4-DOF工业机械臂完整参数

    大臂500mm, 小臂400mm, d1=350mm

    """

    links: List[DHParam]

    joint_names: List[str] = field(default_factory=lambda: [

        'joint1_base_yaw', 'joint2_elbow_yaw',

        'joint3_wrist_roll', 'joint4_z_lift'

    ])

 

    @classmethod

    def default_welding(cls):

        pi = 3.141592653589793

        links = [

            DHParam(0.0, 0.350, 0.500, 0.0),

            DHParam(0.0, 0.0, 0.400, 0.0),

            DHParam(0.0, 0.0, 0.0, 0.0),

            DHParam(0.0, 0.0, 0.0, -pi/2),

        ]

        return cls(links=links)

 

 

@dataclass

class SingularityThreshold:

    """

    奇异检测阈值配置

    det_threshold: 2×2子块行列式绝对值低于此值判定为奇异

    单位：m²（两个长度量纲相乘）

    默认1e-4 m²对应|sin(θ₂)|<5e-4，即θ₂距奇异约0.03°

    """

    det_threshold: float = 1e-4

 

    def __post_init__(self):

        if self.det_threshold <= 0:

            raise ValueError("det_threshold必须为正数")

        if not np.isfinite(self.det_threshold):

            raise ValueError("det_threshold必须为有限数值")

 

"singularity/forward_kinematics.py" — 正运动学求解器

 

"""

正运动学求解器：计算各连杆变换矩阵，供雅可比和奇异检测使用

"""

 

import numpy as np

from singularity.models import SCARAParams

 

 

class ForwardKinematicsSCARA:

    """

    4-DOF SCARA正运动学求解器

    输出各连杆在基座标系中的变换矩阵列表

    """

 

    def __init__(self, params: SCARAParams):

        self.params = params

 

    def compute_all_transforms(self, joint_values):

        """

        计算从基座到各连杆的齐次变换矩阵

        返回：T_list[i] = T_0^{i+1}，i=0,1,2,3

        同时返回各连杆坐标系原点在基座标系中的位置 p_list[i]

        以及各连杆z轴方向在基座标系中的单位向量 z_list[i]

 

        注意：z_list[3]是第4连杆系z轴（经α₄=-π/2翻转），

        不是关节4平移方向；平移方向用z_list[2]（第3连杆z轴）

        """

        if len(joint_values) != 4:

            raise ValueError(f"需要4个关节值, 实际{len(joint_values)}个")

 

        T_list = []

        p_list = []

        z_list = []

 

        T = np.eye(4)

        for i in range(3):

            self.params.links[i].theta = joint_values[i]

            T_i = self._dh_matrix(self.params.links[i])

            T = T @ T_i

            T_list.append(T.copy())

            p_list.append(T[0:3, 3].copy())

            z_list.append(T[0:3, 2].copy())

 

        # 第4关节：平移，α₄=-π/2

        self.params.links[3].d = joint_values[3]

        T_3 = self._dh_matrix(self.params.links[3])

        T = T @ T_3

        T_list.append(T.copy())

        p_list.append(T[0:3, 3].copy())

        z_list.append(T[0:3, 2].copy())

 

        return T_list, p_list, z_list

 

    @staticmethod

    def _dh_matrix(param):

        ct, st = np.cos(param.theta), np.sin(param.theta)

        ca, sa = np.cos(param.alpha), np.sin(param.alpha)

        return np.array([

            [ct, -st*ca, st*sa, param.a*ct],

            [st, ct*ca, -ct*sa, param.a*st],

            [0, sa, ca, param.d],

            [0, 0, 0, 1]

        ], dtype=float)

 

"singularity/jacobian_solver.py" — 雅可比矩阵求解器

 

"""

雅可比矩阵求解器：输入关节角度，计算解析雅可比矩阵

供奇异检测器调用

"""

 

import numpy as np

from singularity.models import SCARAParams

from singularity.forward_kinematics import ForwardKinematicsSCARA

 

 

class JacobianSolver:

    """

    4-DOF SCARA雅可比矩阵解析求解器

    """

 

    def __init__(self, params: SCARAParams):

        self.params = params

        self.fk = ForwardKinematicsSCARA(params)

 

    def compute_jacobian(self, joint_values) -> np.ndarray:

        """

        计算给定关节角度下的6×4雅可比矩阵

        """

        if len(joint_values) != 4:

            raise ValueError(f"需要4个关节值, 实际{len(joint_values)}个")

 

        T_list, p_list, z_list = self.fk.compute_all_transforms(joint_values)

 

        p_end = p_list[-1]

        z_base = np.array([0.0, 0.0, 1.0])

 

        J = np.zeros((6, 4))

 

        # 第1列：joint1 (旋转关节，绕基座z轴)

        J[0:3, 0] = np.cross(z_base, p_end - np.zeros(3))

        J[3:6, 0] = z_base

 

        # 第2列：joint2 (旋转关节，绕第1连杆z轴)

        J[0:3, 1] = np.cross(z_list[0], p_end - p_list[0])

        J[3:6, 1] = z_list[0]

 

        # 第3列：joint3 (旋转关节，绕第2连杆z轴)

        J[0:3, 2] = np.cross(z_list[1], p_end - p_list[1])

        J[3:6, 2] = z_list[1]

 

        # 第4列：joint4 (平移关节，沿第3连杆z轴方向)

        # 注意：用z_list[2]（第3连杆z轴），不用z_list[3]（经α₄翻转）

        J[0:3, 3] = z_list[2]

        J[3:6, 3] = np.array([0.0, 0.0, 0.0])

 

        return J

 

    def compute_jacobian_2x2(self, joint_values) -> np.ndarray:

        """

        计算末端位置对θ₁、θ₂的2×2雅可比子块（SCARA位置奇异核心）

        x = L1*cos(θ1) + L2*cos(θ1+θ2)

        y = L1*sin(θ1) + L2*sin(θ1+θ2)

        J₂×₂ = [[-L1*sin(θ1)-L2*sin(θ1+θ2), -L2*sin(θ1+θ2)],

                [ L1*cos(θ1)+L2*cos(θ1+θ2), L2*cos(θ1+θ2)]]

        """

        if len(joint_values) != 4:

            raise ValueError(f"需要4个关节值, 实际{len(joint_values)}个")

 

        theta1 = joint_values[0]

        theta2 = joint_values[1]

        L1 = self.params.links[0].a

        L2 = self.params.links[1].a

 

        J22 = np.array([

            [-L1*np.sin(theta1) - L2*np.sin(theta1 + theta2),

             -L2*np.sin(theta1 + theta2)],

            [L1*np.cos(theta1) + L2*np.cos(theta1 + theta2),

             L2*np.cos(theta1 + theta2)]

        ], dtype=float)

 

        return J22

 

"singularity/singularity_detector.py" — 奇异检测器（核心业务类）

 

"""

奇异检测器：计算雅可比2×2子块行列式，依据阈值判定奇异

核心业务：|det(J₂×₂)| = L₁·L₂·|sin(θ₂)| < threshold → 奇异

 

为什么用2×2而不是3×3：

SCARA前3轴全竖直，θ₃是腕部roll（a₃=0），第4关节只沿z平移

→ p_end - p₂ 是纯竖直向量

→ J_v第3列 = z₂ × (p_end - p₂) = [0,0,1] × [0,0,d₄] = [0,0,0] 恒为零

→ 3×3子块行列式恒等于零，与θ₁θ₂θ₃取值无关

→ 全部构型被误判为奇异，指标失去意义

"""

 

import numpy as np

from singularity.models import SCARAParams, SingularityThreshold

from singularity.jacobian_solver import JacobianSolver

 

 

class SingularityDetector:

    """

    基于雅可比2×2子块行列式绝对值的奇异位形检测器

    |det(J₂×₂)| = L₁·L₂·|sin(θ₂)|

    单位：m²

    """

 

    def __init__(self, params: SCARAParams,

                 threshold: SingularityThreshold = None):

        self.params = params

        self.L1 = params.links[0].a

        self.L2 = params.links[1].a

        self.threshold = threshold or SingularityThreshold()

        self.jacobian_solver = JacobianSolver(params)

 

    def compute_analytical_determinant(self, joint_values) -> float:

        """

        解析计算2×2子块行列式

        det = L₁·L₂·sin(θ₂)

        这是SCARA位置奇异的精确指标，与θ₁无关

        """

        theta2 = joint_values[1]

        return float(self.L1 * self.L2 * np.sin(theta2))

 

    def compute_jacobian_determinant(self, joint_values) -> float:

        """

        数值计算2×2子块行列式（验证用）

        取J_v的前2行（xy平面）和前2列（θ₁、θ₂）

        """

        J22 = self.jacobian_solver.compute_jacobian_2x2(joint_values)

        return float(np.linalg.det(J22))

 

    def is_singular(self, joint_values, threshold_override=None) -> tuple:

        """

        判断当前构型是否靠近奇异

        返回：(is_singular: bool, det_value: float, threshold_used: float)

        """

        thresh = threshold_override if threshold_override is not None \

                 else self.threshold.det_threshold

        # 使用解析公式（更精确，无数值误差）

        det_val = self.compute_analytical_determinant(joint_values)

        is_sing = abs(det_val) < thresh

        return is_sing, det_val, thresh

 

    def get_distance_to_singularity(self, joint_values) -> float:

        """

        计算当前构型到奇异边界的角距离（度）

        用于量化"还有多远到奇异"

        """

        theta2 = joint_values[1]

        # 奇异发生在θ₂=0或θ₂=±180°

        dist_to_zero = abs(theta2)

        dist_to_180 = abs(abs(theta2) - np.pi)

        return float(min(dist_to_zero, dist_to_180) * 180 / np.pi)

 

    def scan_workspace_grid(self, r_range, theta2_range):

        """

        在工作空间网格上解析生成行列式分布

        由于det = L₁·L₂·sin(θ₂)与θ₁无关，也近似与r无关

        （给定θ₂后r由θ₁决定，但det只依赖θ₂）

        返回：grid_data[i][j] = (r, theta2, det_val, is_singular)

        """

        grid_data = []

 

        for r in r_range:

            row = []

            for theta2 in theta2_range:

                # 解析计算行列式（与θ₁无关）

                det_val = self.L1 * self.L2 * np.sin(theta2)

                is_sing = abs(det_val) < self.threshold.det_threshold

                row.append((r, theta2, det_val, is_sing))

            grid_data.append(row)

 

        return grid_data

 

"singularity/visualize.py" — 可视化模块

 

"""

奇异检测仿真结果可视化：4张核心图表

"""

 

import numpy as np

import matplotlib.pyplot as plt

 

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']

plt.rcParams['axes.unicode_minus'] = False

 

 

class Visualizer:

    def __init__(self, results_dir="results"):

        self.results_dir = results_dir

        import os

        os.makedirs(results_dir, exist_ok=True)

 

    def plot_determinant_vs_theta2(self, theta2_values, det_values,

                                    threshold, is_singular_list):

        """

        ① 行列式绝对值随θ₂角度变化曲线

        工程结论：θ₂→0°或±180°时行列式趋零（肘部伸直/折叠奇异）

        """

        fig, ax = plt.subplots(figsize=(10, 5))

 

        ax.plot(np.array(theta2_values) * 180 / np.pi,

                 np.abs(det_values), 'o-', linewidth=1.5,

                 markersize=5, color='steelblue')

 

        # 奇异阈值线

        ax.axhline(y=threshold, color='red', linestyle='--',

                    linewidth=1.5, label=f'Singularity Threshold '

                    f'({threshold} m²)')

 

        # 标注奇异区域

        for i, sing in enumerate(is_singular_list):

            if sing:

                ax.axvline(x=np.array(theta2_values[i]) * 180 / np.pi,

                            color='red', linestyle=':', alpha=0.5)

 

        ax.set_xlabel('Joint2 Angle θ₂ (deg)')

        ax.set_ylabel('|det(J₂×₂)| (m²)')

        ax.set_title('Jacobian 2×2 Sub-block Determinant vs. Elbow Angle θ₂\n'

                      '(|det| = L₁·L₂·|sin(θ₂)|)')

        ax.set_yscale('log')

        ax.grid(True, alpha=0.3)

        ax.legend(fontsize=8)

 

        plt.tight_layout()

        path = f"{self.results_dir}/01_determinant_vs_theta2.png"

        plt.savefig(path, dpi=150)

        plt.close()

        return path

 

    def plot_singularity_heatmap(self, grid_data, r_range, theta2_range,

                                  weld_positions=None):

        """

        ② 工作空间奇异检测热力图

        工程结论：红色区域为奇异区，轨迹规划需避开

        标注叙事中的焊缝位置（第47~52号焊缝在拐角处）

        """

        fig, ax = plt.subplots(figsize=(10, 7))

 

        # 构建行列式绝对值矩阵

        det_matrix = np.zeros((len(r_range), len(theta2_range)))

        sing_matrix = np.zeros((len(r_range), len(theta2_range)))

 

        for i, row in enumerate(grid_data):

            for j, (r, theta2, det_val, is_sing) in enumerate(row):

                det_matrix[i, j] = abs(det_val)

                sing_matrix[i, j] = 1 if is_sing else 0

 

        # 热力图

        R_mesh, T2_mesh = np.meshgrid(

            np.array(r_range) * 1000,

            np.array(theta2_range) * 180 / np.pi)

        c = ax.contourf(R_mesh, T2_mesh, det_matrix.T, levels=20,

                         cmap='viridis')

        # 奇异边界等高线

        ax.contour(R_mesh, T2_mesh, sing_matrix.T, levels=[0.5],

                    colors='red', linewidths=2, alpha=0.8,

                    label='Singularity Boundary')

 

        plt.colorbar(c, ax=ax, label='|det(J₂×₂)| (m²)')

        ax.set_xlabel('Horizontal Reach Distance r (mm)')

        ax.set_ylabel('Joint2 Angle θ₂ (deg)')

        ax.set_title('Singularity Map: |det(J₂×₂)| = L₁·L₂·|sin(θ₂)|\n'

                      '(Red contour: singularity boundary)')

 

        # 标注叙事中的焊缝位置（拐角处，θ₂接近0°）

        if weld_positions:

            for wp in weld_positions:

                ax.plot(wp[0] * 1000, wp[1] * 180 / np.pi, 'r*',

                         markersize=15, label='Weld Seam 47-52' if wp == weld_positions[0] else '')

 

        ax.legend(fontsize=8)

 

        plt.tight_layout()

        path = f"{self.results_dir}/02_singularity_heatmap.png"

        plt.savefig(path, dpi=150)

        plt.close()

        return path

 

    def plot_det_vs_reach(self, reach_distances, det_values,

                           threshold, is_singular_list):

        """

        ③ 行列式绝对值随水平距离变化曲线

        工程结论：由于det与r无关（只依赖θ₂），曲线应为水平线

        若不同r对应的θ₂不同，则展示θ₂变化导致的行列式变化

        """

        fig, ax = plt.subplots(figsize=(10, 5))

 

        ax.plot(np.array(reach_distances) * 1000,

                 np.abs(det_values), 'o-', linewidth=1.5,

                 markersize=5, color='steelblue')

 

        ax.axhline(y=threshold, color='red', linestyle='--',

                    linewidth=1.5, label=f'Threshold ({threshold} m²)')

 

        for i, sing in enumerate(is_singular_list):

            if sing:

                ax.plot(np.array(reach_distances[i]) * 1000,

                         abs(det_values[i]), 'ro', markersize=8,

                         fillstyle='none')

 

        ax.set_xlabel('Horizontal Reach Distance r (mm)')

        ax.set_ylabel('|det(J₂×₂)| (m²)')

        ax.set_title('Jacobian Determinant vs. Horizontal Reach Distance\n'

                      '(det depends only on θ₂, independent of r)')

        ax.set_yscale('log')

        ax.grid(True, alpha=0.3)

        ax.legend(fontsize=8)

 

        plt.tight_layout()

        path = f"{self.results_dir}/03_det_vs_reach.png"

        plt.savefig(path, dpi=150)

        plt.close()

        return path

 

    def plot_analytical_vs_numerical(self, test_configs, analytical_dets,

                                      numerical_dets):

        """

        ④ 解析公式 vs 数值计算对比

        工程结论：两者完全一致，验证解析公式正确性

        """

        fig, ax = plt.subplots(figsize=(10, 6))

 

        n = len(test_configs)

        x = np.arange(n)

        width = 0.35

 

        ax.bar(x - width/2, np.abs(analytical_dets), width,

                label='Analytical: L₁·L₂·sin(θ₂)', alpha=0.8,

                color='steelblue', edgecolor='black')

        ax.bar(x + width/2, np.abs(numerical_dets), width,

                label='Numerical: det(J₂×₂)', alpha=0.8,

                color='orange', edgecolor='black')

 

        ax.set_xlabel('Test Configuration')

        ax.set_ylabel('|det| (m²)')

        ax.set_title('Analytical vs Numerical Jacobian Determinant')

        ax.set_xticks(x)

        ax.set_xticklabels([cfg['label'] for cfg in test_configs],

                           rotation=45, ha='right')

        ax.legend(fontsize=8)

        ax.grid(True, alpha=0.3, axis='y')

 

        plt.tight_layout()

        path = f"{self.results_dir}/04_analytical_vs_numerical.png"

        plt.savefig(path, dpi=150)

        plt.close()

        return path

 

5.

利用AI解决实际问题，如果你觉得这个工具好用，欢迎关注长安牧笛！