《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解决实际问题，如果你觉得这个工具好用，欢迎关注长安牧笛！
