# 双线摆转动惯量与阻尼参数辨识（MATLAB/Simulink）

本文件夹包含视频中使用的双线摆实验数据、数据预处理脚本、Simulink 动力学模型，以及 MATLAB Parameter Estimator 的配置代码。它用于从自由衰减角度数据中估计转动惯量 `I`、粘性阻尼系数 `C` 和平方阻尼系数 `KD`。

> 这里提供的是可修改的实验示例，不是对所有双线摆装置都通用的“直接运行即得答案”程序。使用自己的数据时，必须核对质量、悬线几何尺寸、测量轴、安装角和拟合区间。

## 项目如何工作（重要）

`BP_est_SmallAngleApprox.slx` 保存了双线摆系统的 Simulink 动力学模型，它是整个参数辨识过程的核心，并不是只用于显示结果的附件。

MATLAB Simulink Parameter Estimator（Parameter Estimation Tool）以这个 `.slx` 模型为被估计对象。辨识过程中，工具会不断改变模型中的 `I`、`C`、`KD`，反复运行 Simulink 仿真，并比较：

```text
实验测得的角度曲线  <----误差---->  Simulink 模型输出的角度曲线
                                │
                                └── 调整 I、C、KD 后重新仿真
```

直到仿真输出与实验数据之间的误差不再明显下降，工具才给出估计参数。因此，本项目的关系是：

```text
CSV 实验数据
    ↓
estimator.m：读取、坐标修正、去除中心漂移、截取数据
    ↓
init_params.m：提供物理参数以及 I、C、KD 的初值/边界/尺度
    ↓
parameterEstimation_setup.m：把数据、参数和初始状态装入估计会话
    ↓
Parameter Estimator 反复运行 BP_est_SmallAngleApprox.slx
    ↓
得到 I、C、KD 的估计值，并比较实测与仿真曲线
```

如果 `.slx` 模型中的动力学方程、参数名称、输入输出端口或积分器状态被修改，相应的 MATLAB 脚本配置也必须同步修改。反过来，仅修改脚本而不理解模型中的方程，也可能得到数值上收敛但物理上错误的参数。

## 1. 文件说明

| 文件或目录 | 用途 |
|---|---|
| `estimator.m` | 主入口：选择实验、读取 CSV、修正传感器安装角、截取拟合及验证区间并打开 Parameter Estimator。 |
| `init_params.m` | 各类实验的质量、几何尺寸、参数初值、上下界、缩放量和是否参与辨识。 |
| `find_envelope.m` | 检测正负峰值，计算振幅包络和振动中心，用来去除缓慢变化的角度漂移。 |
| `parameterEstimation_setup.m` | 创建 Simulink Design Optimization 参数估计会话；主要由 MATLAB 自动生成。 |
| `BP_est_SmallAngleApprox.slx` | 双线摆系统的 Simulink 动力学模型；Parameter Estimator 以它为基础反复仿真并辨识 `I`、`C`、`KD`。 |
| `BP_est_SmallAngleApprox_R2018b.slx` | 为较早 MATLAB 版本保存的模型副本。 |
| `data/` | 示例实验 CSV 数据。 |
| `slprj/`、`*.slxc` | Simulink 自动生成的缓存，可以删除；再次运行模型时会重新生成。 |

## 2. 软件要求

建议使用 MATLAB，并安装：

- Simulink
- Simulink Design Optimization（Parameter Estimator / `sdo`）
- Signal Processing Toolbox（`findpeaks`）
- Aerospace Toolbox（`angle2dcm`）
- Mapping Toolbox（`wrapToPi`；若没有该工具箱，可用 `mod(x+pi,2*pi)-pi` 替代）

模型的当前版本由较新的 MATLAB 保存。如果无法打开，请尝试 `BP_est_SmallAngleApprox_R2018b.slx`，并相应修改脚本中的模型名称。

## 3. 最快上手

1. 在 MATLAB 中把“当前文件夹”切换到本目录。
2. 打开 `estimator.m`。
3. 设置 `expDataFile`，选择 `data` 文件夹中的 CSV。
4. 设置 `exp_type`，它必须与 `init_params.m` 中的分支名称一致。
5. 在 `init_params.m` 中核对质量 `m`、悬线间距 `D`、悬线长度 `h`，以及 `I`、`C`、`KD` 的初值和边界。
6. 在 `estimator.m` 中根据曲线修改拟合区间 `t1～t2` 和验证区间 `t3～t4`。
7. 运行 `estimator.m`。脚本完成预处理后，会基于 `BP_est_SmallAngleApprox.slx` 创建并打开 MATLAB Parameter Estimator 会话。
8. 在估计器中检查实验信号、Simulink 模型输出和待辨识参数，然后开始 Estimate。工具会反复运行该 `.slx` 模型并调整 `I`、`C`、`KD`；不要只看拟合区间，还要用未参与拟合的数据验证结果。

第一次运行建议直接使用：

```matlab
estimator
```

## 4. CSV 数据格式

示例数据采样率为 100 Hz，CSV 至少应包含下列列名：

```text
sample_index,time_s,roll_rad,pitch_rad,yaw_rad,
gyro_x_rad_s,gyro_y_rad_s,gyro_z_rad_s
```

- 时间单位为秒。
- 姿态角单位为弧度。
- 角速度单位为弧度每秒。
- 当前主脚本以 yaw 方向和修正后的 Z 轴角速度作为示例。若测量其他转轴，需要检查坐标变换和所选角速度分量。
- 如果实际采样率不是 100 Hz，请同时修改 `estimator.m` 中的 `fs`。

## 5. 运行前必须修改或核对的项目

### 实验和数据文件

```matlab
expDataFile = "block_alone.csv";
exp_type = 'aluminum_block';
```

可用实验类型见 `init_params.m`：

- `aluminum_block`：铝块标定实验
- `airplane_x`：航模 X 轴实验
- `airplane_y`：航模 Y 轴实验
- `airplane_z`：航模 Z 轴实验

### 安装角及手动角度偏移

`estimator.m` 根据平均 roll、pitch 估计传感器安装倾角，并用 `manual_offset` 处理剩余的航向安装偏移。传感器安装方式变化后，应重新检查这些处理是否符合实际坐标定义。

### 拟合与验证区间

```matlab
t1 = 25.9;
t2 = 35;
t3 = 15;
t4 = 40;
```

这些数值只适合当前示例数据。选择区间时应避开释放动作、碰撞、明显丢样和传感器饱和。验证区间最好包含未参与拟合的数据；如需让估计器同时载入验证实验，可查看 `parameterEstimation_setup.m` 中 `si.Experiments` 附近的注释。

### 峰值检测条件

`find_envelope.m` 中：

```matlab
"MinPeakDistance", 0.6
"MinPeakProminence", 0.03
```

分别约束峰值的最小时间间隔和最小显著度。振动周期或噪声水平改变时，需要调整，否则可能漏检真实峰值或把噪声当作峰值。

## 6. 参数含义和单位

模型中的阻尼力矩写成：

```text
粘性阻尼力矩：C  · θ̇
平方阻尼力矩：KD · θ̇ · |θ̇|
```

| 参数 | 含义 | 单位 |
|---|---|---|
| `I` | 绕实验转轴的转动惯量 | kg·m² |
| `C` | 粘性转动阻尼系数 | N·m·s/rad |
| `KD` | 平方转动阻尼系数 | N·m·s²/rad² |
| `m` | 被测物及随动物体总质量 | kg |
| `D` | 双线之间的有效间距（以模型定义为准） | m |
| `h` | 悬线有效长度 | m |

弧度在量纲分析中也可视为无量纲。

## 7. 参数估计建议

- 先用周期给 `I` 一个合理初值，再估计阻尼。
- `Value` 必须位于 `Minimum` 和 `Maximum` 之间；分享版本中的数值是实验过程示例，运行自己的数据前务必逐项检查。
- `Scale` 应接近参数的典型数量级，例如 `I≈10^-3` 时可把 Scale 设为 `1e-3`。
- 只想辨识某个参数时，将其他参数的 `Free` 设为 `false`。
- 不要仅凭拟合曲线判断结果；应检查验证区间、多个释放振幅及多次重复实验。
- 如果 `C` 和 `KD` 同时自由时结果随初值大幅变化，说明两种阻尼可能难以从当前数据中独立辨识。可先只辨识 `C`，再比较加入 `KD` 后验证误差是否稳定下降。

## 8. 常见问题

### 找不到 CSV

确认 MATLAB 当前文件夹是本目录。`estimator.m` 会把相对路径 `data` 加入搜索路径。

### `findpeaks` 报错或找不到函数

安装 Signal Processing Toolbox，并检查所选区间内是否同时存在足够的正峰和负峰。

### `angle2dcm` 或 `wrapToPi` 找不到

检查相应工具箱是否安装。`wrapToPi(x)` 可以替换成：

```matlab
x = mod(x + pi, 2*pi) - pi;
```

### 插值结果包含 NaN

拟合或验证区间超出了 CSV 的实际时间范围。确保：

```text
t(1) <= t1 < t2 <= t(end)
t(1) <= t3 < t4 <= t(end)
```

### 参数估计一开始就停止或结果卡在边界

检查参数初值、上下界和 Scale，特别确认初值处于边界以内。然后再检查单位、质量和几何尺寸是否正确。

## 9. 使用自己的数据时的推荐流程

1. 先绘制原始角度和角速度，检查时间戳、丢样和饱和。
2. 确认传感器坐标轴与被测转轴的对应关系。
3. 去除安装倾角和缓慢漂移，但不要滤掉真实振动。
4. 选择一段干净数据拟合，保留另一段或另一组实验验证。
5. 设置具有物理意义的初值、边界和 Scale。
6. 分阶段估计 `I`、`C`、`KD`，最后再考虑小范围联合优化。
7. 报告参数时同时注明实验几何、数据区间、采样率和使用的阻尼模型。

## 10. 免责声明

本项目用于实验演示和参数辨识学习。测量精度会受到悬线长度与间距误差、安装偏心、传感器漂移、释放方式、结构柔性、空气阻力以及模型假设等因素影响。涉及工程设计或安全判断时，请进行独立校准和不确定度分析。
