系统辨识是连接仿真与现实的关键桥梁。本章深入探讨如何从实验数据中准确估计物理系统参数，包括质量、惯性、摩擦系数等关键物理量。我们将从经典的最小二乘方法出发，逐步深入到贝叶斯框架下的参数估计，并特别关注接触参数的标定问题。这些技术对于构建高保真度的数字孪生系统至关重要，是实现仿真到现实成功迁移的基础。
在机器人系统中，动力学方程往往可以重新整理为参数的线性形式：
\[\mathbf{Y}(\mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}}) \boldsymbol{\theta} = \boldsymbol{\tau}\]其中 $\mathbf{Y}$ 是回归矩阵（regressor matrix），$\boldsymbol{\theta}$ 是待辨识的参数向量，包含：
对于第 $i$ 个连杆，其最小参数集表示为：
\[\boldsymbol{\theta}_i = [m_i, m_i r_{cx,i}, m_i r_{cy,i}, m_i r_{cz,i}, I_{xx,i}, I_{yy,i}, I_{zz,i}, I_{xy,i}, I_{xz,i}, I_{yz,i}]^T\]系统辨识的成功很大程度上取决于激励轨迹的质量。理想的激励轨迹应该：
常用的激励轨迹形式：
\[q_j(t) = \sum_{k=1}^{N_f} \left[ a_{jk} \sin(\omega_k t) + b_{jk} \cos(\omega_k t) \right]\]其中频率 $\omega_k$ 的选择需要避免共振，通常选择为基频的整数倍。
优化激励轨迹的目标函数：
\[\min_{\mathbf{a}, \mathbf{b}} \text{cond}(\mathbf{Y}^T\mathbf{Y})\]其中 $\text{cond}(\cdot)$ 是条件数，用于衡量数值稳定性。
对于在线参数估计，递归最小二乘提供了计算高效的更新方式：
\[\hat{\boldsymbol{\theta}}_{k+1} = \hat{\boldsymbol{\theta}}_k + \mathbf{K}_{k+1} (\boldsymbol{\tau}_{k+1} - \mathbf{Y}_{k+1} \hat{\boldsymbol{\theta}}_k)\]其中增益矩阵：
\[\mathbf{K}_{k+1} = \mathbf{P}_k \mathbf{Y}_{k+1}^T (\lambda + \mathbf{Y}_{k+1} \mathbf{P}_k \mathbf{Y}_{k+1}^T)^{-1}\]协方差矩阵更新：
\[\mathbf{P}_{k+1} = \frac{1}{\lambda} (\mathbf{I} - \mathbf{K}_{k+1} \mathbf{Y}_{k+1}) \mathbf{P}_k\]遗忘因子 $\lambda \in (0, 1]$ 用于处理时变参数。
辨识得到的参数必须满足物理约束：
惯性张量正定性： \(\mathbf{I}_i = \begin{bmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{xy} & I_{yy} & I_{yz} \\ I_{xz} & I_{yz} & I_{zz} \end{bmatrix} \succ 0\)
带约束的最小二乘问题可以表述为半定规划（SDP）：
\[\begin{aligned} \min_{\boldsymbol{\theta}} & \quad \|\mathbf{Y}\boldsymbol{\theta} - \boldsymbol{\tau}\|^2 \\ \text{s.t.} & \quad \mathbf{I}_i(\boldsymbol{\theta}) \succeq \epsilon \mathbf{I}_3, \quad \forall i \\ & \quad m_i \geq \epsilon_m \end{aligned}\]贝叶斯方法将参数估计问题转化为后验分布的推断：
\[p(\boldsymbol{\theta} | \mathcal{D}) = \frac{p(\mathcal{D} | \boldsymbol{\theta}) p(\boldsymbol{\theta})}{p(\mathcal{D})}\]其中：
| $p(\mathcal{D} | \boldsymbol{\theta})$ 是似然函数，描述观测数据的生成过程 |
| $p(\boldsymbol{\theta} | \mathcal{D})$ 是后验分布，包含了参数的不确定性信息 |
对于动力学辨识，似然函数通常假设为高斯噪声模型：
\[p(\boldsymbol{\tau} | \mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}}, \boldsymbol{\theta}) = \mathcal{N}(\mathbf{Y}\boldsymbol{\theta}, \boldsymbol{\Sigma}_{\tau})\]选择共轭先验可以得到解析的后验分布。对于高斯似然，正态-逆威沙特分布是共轭先验：
\[p(\boldsymbol{\theta}, \boldsymbol{\Sigma}) = \text{NIW}(\boldsymbol{\mu}_0, \boldsymbol{\Lambda}_0, \nu_0, \boldsymbol{\Psi}_0)\]后验更新公式：
\[\begin{aligned} \boldsymbol{\mu}_n &= \frac{\boldsymbol{\Lambda}_0 \boldsymbol{\mu}_0 + n\bar{\boldsymbol{\theta}}}{\boldsymbol{\Lambda}_0 + n} \\ \boldsymbol{\Lambda}_n &= \boldsymbol{\Lambda}_0 + n \\ \nu_n &= \nu_0 + n \\ \boldsymbol{\Psi}_n &= \boldsymbol{\Psi}_0 + \mathbf{S} + \frac{n\boldsymbol{\Lambda}_0}{n + \boldsymbol{\Lambda}_0}(\bar{\boldsymbol{\theta}} - \boldsymbol{\mu}_0)(\bar{\boldsymbol{\theta}} - \boldsymbol{\mu}_0)^T \end{aligned}\]对于复杂的后验分布，MCMC方法提供了通用的采样框架。Hamiltonian Monte Carlo (HMC) 特别适合高维参数空间：
动量变量引入： \(p(\boldsymbol{\theta}, \mathbf{p}) = p(\boldsymbol{\theta}) p(\mathbf{p}) = p(\boldsymbol{\theta}) \mathcal{N}(\mathbf{0}, \mathbf{M})\)
哈密顿动力学： \(\begin{aligned} \frac{d\boldsymbol{\theta}}{dt} &= \mathbf{M}^{-1} \mathbf{p} \\ \frac{d\mathbf{p}}{dt} &= -\nabla_{\boldsymbol{\theta}} U(\boldsymbol{\theta}) \end{aligned}\)
| 其中 $U(\boldsymbol{\theta}) = -\log p(\boldsymbol{\theta} | \mathcal{D})$ 是势能函数。 |
变分推断通过优化寻找最接近真实后验的分布：
\[q^*(\boldsymbol{\theta}) = \arg\min_{q \in \mathcal{Q}} \text{KL}(q(\boldsymbol{\theta}) \| p(\boldsymbol{\theta} | \mathcal{D}))\]对于平均场近似： \(q(\boldsymbol{\theta}) = \prod_i q_i(\theta_i)\)
坐标上升更新： \(\log q_j^*(\theta_j) = \mathbb{E}_{q_{-j}}[\log p(\boldsymbol{\theta}, \mathcal{D})] + \text{const}\)
摩擦参数的辨识需要特殊设计的实验。对于库仑摩擦模型：
\[\mathbf{f}_t = \begin{cases} -\mu_s \|\mathbf{f}_n\| \frac{\mathbf{v}_t}{\|\mathbf{v}_t\|} & \text{if } \|\mathbf{f}_t^{\text{trial}}\| > \mu_s \|\mathbf{f}_n\| \\ \mathbf{f}_t^{\text{trial}} & \text{otherwise} \end{cases}\]标定实验设计：
倾斜平面法：测量临界倾角 $\theta_c$ \(\mu_s = \tan(\theta_c)\)
拉力测试：施加水平力直到滑动 \(\mu_s = \frac{F_{\text{critical}}}{mg}\)
振荡运动：通过能量耗散估计动摩擦系数 \(E_{\text{lost}} = \mu_k mg \cdot d\)
恢复系数 $e$ 定义了碰撞前后的相对速度关系：
\[e = -\frac{v_{\text{sep}}}{v_{\text{app}}}\]标准测量方法：
垂直跌落测试： \(e = \sqrt{\frac{h_{\text{bounce}}}{h_{\text{drop}}}}\)
摆锤碰撞：通过角度测量 \(e = \frac{\cos(\theta_{\text{min}}) - \cos(\theta_0)}{\cos(\theta_{\text{max}}) - \cos(\theta_0)}\)
高速摄像分析：直接测量碰撞前后速度
恢复系数的速度依赖性模型： \(e(v) = e_0 \cdot \exp(-\alpha v)\)
柔性接触模型（Hunt-Crossley）的参数辨识：
\[\mathbf{f}_n = k_n \delta^n - c_n \delta^n \dot{\delta}\]其中：
通过力传感器数据拟合：
\[\min_{k_n, c_n, n} \sum_i \|\mathbf{f}_{\text{measured},i} - \mathbf{f}_{\text{model},i}(k_n, c_n, n)\|^2\]使用Levenberg-Marquardt算法求解非线性最小二乘问题。
对于面接触，压力分布的测量对于验证接触模型至关重要：
接触面积与法向力的关系（对于弹性接触）： \(A = \pi \left(\frac{3FR}{4E^*}\right)^{2/3}\)
其中有效弹性模量： \(\frac{1}{E^*} = \frac{1-\nu_1^2}{E_1} + \frac{1-\nu_2^2}{E_2}\)
将数据集分为训练集和验证集，评估模型的泛化能力：
将数据分为k个子集
for i = 1 to k:
使用除第i个子集外的所有数据训练
在第i个子集上验证
计算平均性能指标
留一法（LOO-CV）：k = n的特殊情况
定义归一化预测误差：
\[\text{NPE} = \frac{1}{N} \sum_{i=1}^N \frac{\|\boldsymbol{\tau}_{\text{pred},i} - \boldsymbol{\tau}_{\text{meas},i}\|}{\|\boldsymbol{\tau}_{\text{meas},i}\|}\]频域分析通过传递函数评估：
\[H(\omega) = \frac{\mathcal{F}\{\boldsymbol{\tau}_{\text{pred}}\}}{\mathcal{F}\{\boldsymbol{\tau}_{\text{cmd}}\}}\]相位和幅值裕度提供稳定性指标。
在闭环控制下验证模型：
位置跟踪误差： \(e_p = \sqrt{\frac{1}{T} \int_0^T \|\mathbf{q}_{\text{ref}}(t) - \mathbf{q}(t)\|^2 dt}\)
力矩预测误差： \(e_{\tau} = \sqrt{\frac{1}{T} \int_0^T \|\boldsymbol{\tau}_{\text{model}}(t) - \boldsymbol{\tau}_{\text{real}}(t)\|^2 dt}\)
能量一致性： \(\Delta E = |E_{\text{model}} - E_{\text{real}}|\)
使用统计方法验证模型假设：
残差白度检验（Ljung-Box test）： \(Q = n(n+2) \sum_{k=1}^h \frac{r_k^2}{n-k}\) 其中 $r_k$ 是滞后k的自相关系数
正态性检验（Shapiro-Wilk test）： 验证残差的高斯分布假设
异方差性检验（Breusch-Pagan test）： 检测误差方差的非恒定性
平衡模型复杂度与拟合质量：
赤池信息准则（AIC）： \(\text{AIC} = 2k - 2\ln(L)\)
贝叶斯信息准则（BIC）： \(\text{BIC} = k\ln(n) - 2\ln(L)\)
最小描述长度（MDL）： \(\text{MDL} = -\ln(L) + \frac{k}{2}\ln(n)\)
其中 $k$ 是参数数量，$L$ 是似然值，$n$ 是样本数。
仿真与现实之间的参数差异是策略迁移失败的主要原因之一。典型的参数不确定性来源：
强化学习训练过程中的参数自适应：
初始化: θ_sim 从制造商规格
while 训练进行中:
1. 在真实机器人上执行探索动作
2. 收集状态-动作-力矩数据 (q, q̇, q̈, τ)
3. 使用RLS更新参数估计 θ̂_real
4. 计算参数不确定性范围 Δθ = |θ̂_real - θ_sim|
5. 在仿真中使用 θ ∼ Uniform(θ̂_real - Δθ, θ̂_real + Δθ)
6. 训练策略对参数变化的鲁棒性
使用MAML（Model-Agnostic Meta-Learning）加速参数辨识：
\[\theta^* = \arg\min_\theta \sum_{i=1}^N \mathcal{L}(\theta - \alpha \nabla_\theta \mathcal{L}_i(\theta))\]其中每个任务 $i$ 对应不同的机器人实例或工作条件。
设计信息量最大的探索动作：
\[\mathbf{u}^* = \arg\max_{\mathbf{u}} I(\boldsymbol{\theta}; \mathbf{y} | \mathbf{u}, \mathcal{D})\]其中 $I(\cdot;\cdot)$ 是互信息，衡量动作 $\mathbf{u}$ 对参数 $\boldsymbol{\theta}$ 的信息增益。
实践中使用近似方法：
Karl Johan Åström（1934-）是现代控制理论的奠基人之一，他在1970年代开创性地将统计方法引入控制系统设计。
Åström的工作直接影响了：
他的名言：”All models are wrong, but some are useful”深刻影响了系统辨识领域的哲学思考。
将物理定律作为神经网络训练的约束：
\[\mathcal{L} = \mathcal{L}_{\text{data}} + \lambda \mathcal{L}_{\text{physics}}\]其中物理损失项： \(\mathcal{L}_{\text{physics}} = \|\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) - \boldsymbol{\tau}\|^2\)
将动力学建模为连续时间系统：
\[\frac{d\mathbf{x}}{dt} = f_\theta(\mathbf{x}, t)\]其中 $f_\theta$ 是神经网络参数化的向量场。
优势：
保持能量守恒的神经网络架构：
\[\begin{aligned} \frac{d\mathbf{q}}{dt} &= \frac{\partial H_\theta}{\partial \mathbf{p}} \\ \frac{d\mathbf{p}}{dt} &= -\frac{\partial H_\theta}{\partial \mathbf{q}} \end{aligned}\]其中 $H_\theta(\mathbf{q}, \mathbf{p})$ 是神经网络表示的哈密顿量。
使用高斯过程捕获模型不确定性：
\[\mathbf{f}(\mathbf{x}) \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}'))\]核函数选择：
稀疏高斯过程用于大规模数据： \(q(\mathbf{f}) = \int p(\mathbf{f} | \mathbf{u}) q(\mathbf{u}) d\mathbf{u}\)
其中 $\mathbf{u}$ 是诱导点。
结合物理模型与数据驱动补偿：
\[\boldsymbol{\tau} = \boldsymbol{\tau}_{\text{physics}}(\mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}}, \boldsymbol{\theta}) + \boldsymbol{\tau}_{\text{residual}}(\mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}})\]其中：
这种方法保持了可解释性，同时提高了预测精度。
系统辨识是连接理论模型与实际系统的关键技术。本章介绍了从经典最小二乘到现代贝叶斯方法的参数估计技术，强调了物理一致性约束的重要性。关键要点：
核心公式回顾：
| 贝叶斯后验：$p(\boldsymbol{\theta} | \mathcal{D}) \propto p(\mathcal{D} | \boldsymbol{\theta}) p(\boldsymbol{\theta})$ |
练习 12.1：最小二乘辨识 一个2自由度机械臂的动力学方程可以写成线性参数形式 $\mathbf{Y}\boldsymbol{\theta} = \boldsymbol{\tau}$。给定100组测量数据 $(q_i, \dot{q}_i, \ddot{q}_i, \tau_i)$，如何判断激励轨迹是否充分？
Hint：检查回归矩阵 $\mathbf{Y}$ 的秩和条件数。
练习 12.2：递归最小二乘 推导遗忘因子 $\lambda = 0.95$ 时，历史数据的有效窗口长度（即多少步之前的数据权重降到原来的5%）。
Hint：权重按 $\lambda^k$ 指数衰减。
练习 12.3：摩擦系数测量 设计一个实验测量机器人关节的粘性摩擦系数。假设摩擦模型为 $\tau_f = b\dot{q}$，其中 $b$ 是粘性摩擦系数。
Hint：考虑恒速运动时的力矩平衡。
练习 12.4：贝叶斯参数估计 假设机器人连杆质量的先验分布为 $m \sim \mathcal{N}(1.0, 0.1^2)$ kg。经过n次测量后，后验分布为 $m \sim \mathcal{N}(\mu_n, \sigma_n^2)$。推导需要多少次独立测量才能将不确定性（标准差）降低到0.01 kg？
Hint：使用贝叶斯更新公式，假设测量噪声标准差为 $\sigma_m$。
练习 12.5：物理一致性约束 证明对于任意刚体，其惯性张量必须满足三角不等式。提示：考虑主轴坐标系。
Hint：在主轴坐标系中，惯性张量是对角的。
练习 12.6：主动学习 设计一个基于高斯过程的主动学习策略，用于机器人动力学模型的在线改进。定义获取函数（acquisition function）并说明其物理意义。
Hint：考虑预测不确定性和模型改进的权衡。
练习 12.7：开放性思考题 在仿真到现实迁移中，是否应该让仿真参数完全匹配真实系统？讨论过度拟合真实参数可能带来的问题。
问题：回归矩阵病态导致参数估计不稳定
症状：
解决方案：
# 使用SVD进行稳健求解
U, S, Vt = np.linalg.svd(Y)
# 截断小奇异值
threshold = 1e-6 * S[0]
S_inv = np.where(S > threshold, 1/S, 0)
theta = Vt.T @ np.diag(S_inv) @ U.T @ tau
问题：辨识得到的参数违反物理定律
常见情况：
预防措施：
问题：某些参数无法辨识或辨识精度差
诊断方法：
# 计算Fisher信息矩阵
FIM = Y.T @ Y
# 检查可辨识性
eigenvalues = np.linalg.eigvals(FIM)
if min(eigenvalues) < 1e-6:
print("参数不可辨识")
改进激励轨迹：
问题：使用错误的模型结构
例子：
检测方法：
问题：忽略参数估计的不确定性
后果：
正确做法：