【MATLAB源码-第363期】基于matlab的马尔可夫链Markov Chain蒙特卡罗算法的倒立摆模型参数估计。

操作环境：

MATLAB 2024a

1、算法描述

该系统是一个结合数值仿真与贝叶斯推断的完整闭环流程，用于对旋转倒立摆模型中的若干关键参数进行估计与识别。首先，用户需要给定一个旋转倒立摆模型的动力学方程表达式，在本例中已经由若干MATLAB函数完成了该动力学模型的定义。系统工作正式从加载已知的仿真数据开始，即通过之前使用某种输入电压激励（如线性调频信号）对旋转倒立摆进行仿真，并将该仿真产生出的时间序列与对应的状态变量（包含臂角、臂角速度、摆角、摆角速度等）保存到一个MAT文件中。这里的仿真数据可以看作是“真实测量”或参考数据，后续参数估计的核心目标便是让估计出的模型在相同输入条件下，其输出轨迹尽可能地逼近这组参考数据。

在系统启动时，将事先准备好的MAT文件载入内存，从而得到时间向量与参考状态矩阵。这些数据将作为后续误差计算与似然函数评估的基准。在本系统中，动力学方程以函数形式存在，即rip_dynamics1和rip_dynamics这两个函数分别定义了带固定参数和带待估参数的非线性动力学，具体包括了摆与臂的质量、长度、转动惯量、阻尼系数，以及电机与齿轮的传动比、效率、电阻和电机常数等物理量。这样一来，仿真部分就拥有灵活的接口：对固定参数，可以直接调用rip_dynamics1进行数值积分；对待估计的参数，rip_dynamics函数能够在运行时将这些参数传递给动力学方程，从而根据给定的参数向量计算系统的响应。

仿真核心依赖MATLAB自带的高精度常微分方程求解器ode45。该求解器基于经典四阶龙格－库塔算法，能够在保证数值稳定性的前提下，自适应地调整步长，以满足用户指定的相对与绝对误差容限。具体来说，在仿真过程中，每当需要得知系统在某个时刻的状态时，ode45就会根据当前状态信息、传入的函数句柄以及仿真的时间步长范围自动计算下一个步骤的近似解。这种数值积分方式不仅能够有效处理无刚性冲击但存在非线性耦合项的旋转倒立摆系统，同时也便于后续在MCMC采样循环中多次调用。每次参数提案生成后，都会调用simulate_system函数，该函数内部使用ode45根据给定参数进行一次完整的仿真，并返回整个时间序列上各状态量的数值解。这里需要注意的是，MATLAB的ode45在内部通过求解器对象存储中间计算信息，然后使用deval函数对给定的时间向量进行插值，从而在用户指定的离散时刻点上一并输出完整的状态轨迹。

有了仿真功能之后，系统接下来要做的便是比较仿真输出与参考数据之间的差异，并据此构造一个基于误差的度量函数。在本例中，选择了加权均方误差作为误差指标，具体流程是先将仿真输出的状态矩阵与参考状态矩阵在相同时间点处进行逐元素相减，得到误差矩阵；然后，为了反映摆角和摆角速度等关键物理变量对整体拟合质量的重要性，系统为各个状态量分配了不同的权重，这里假定摆角和摆角速度的权重较高，而臂角与臂角速度的权重可以设为零或很小。这样做的目的在于，如果在参数估计过程中让摆的运动与参考数据保持高度一致，对于控制器设计或后续实验验证会更具实用价值。权重矩阵与误差矩阵结合后，通过对每个时间点的加权平方误差求和，再对所有时间点求平均，就得到一个标量形式的加权均方误差值。这个误差值既可以反映模型输出与真实数据整体趋势的偏离情况，也能量化评估某一组参数在匹配参考数据方面的好坏。

在确定了误差度量之后，就进入到贝叶斯推断阶段。系统设计了一套先验分布，通常先验信息来自于经验或者同类电机、齿轮系统的参数统计数据。为了简化，先验分布在本例中全部选用高斯分布，并为每个待估参数指定了一个先验均值与先验标准差，意味着在未观察到数据之前，我们相信这些参数最可能聚集于某个区间内，但也允许一定范围的随机波动。具体参数包括臂阻尼系数、摆阻尼系数、转矩常数、反电动势常数以及电机和齿轮效率等。设定好先验分布之后，就可以将先验概率评价转化为对数先验，形式上会对参数偏离先验均值的程度给予惩罚，即偏离越大则对数先验值越低。

为了将误差与概率空间联系起来，需要构造似然函数。在本系统中，我们假设测量噪声服从高斯分布，均值为零且方差可以根据实验情况预先设定。在给定某一组参数下进行仿真后，可以得到仿真输出与参考数据之间的偏差曲线，将该偏差除以噪声标准差并做平方和，就可以得到一个拟合度量——误差越小意味着仿真输出更加接近参考数据。通过将该度量映射到对数似然空间，实际就是在计算负半乘法将误差平方和缩放后取负值，从而保证在贝叶斯框架下似然函数符合高斯噪声假设：拟合效果越好，对数似然就越大。

有了对数先验与对数似然之后，就可以得到对数后验，这部分工作由compute_log_posterior函数完成。它会先调用simulate_system_mcmc函数进行一次完整的仿真，得到当前参数下的摆角轨迹；然后将参考摆角轨迹与当前仿真摆角轨迹之间的残差通过高斯噪声假设计算对数似然值，再将当前参数与先验分布的偏离程度计算对数先验值，两者相加便得到对数后验值。这样一来，我们就有了一个衡量当前参数优劣的标量指标：对数后验值越大，说明在满足先验分布的前提下更能匹配观测数据。

核心算法部分是马尔可夫链蒙特卡罗采样，具体使用Metropolis-Hastings算法。在开始正式迭代之前，用户需设定链长度和弃置期等超参数，例如链上采样迭代次数可能设为一万步，而弃置期则取前一半作为“热身”阶段。链的初始位置通常选取先验分布的均值作为初始参数向量，然后利用compute_error函数得到对应的初始误差，再根据噪声假设计算出初始的对数似然，最后与对数先验相加得出初始对数后验值。完成初始化后，进入迭代循环。在每一次迭代中，系统首先在参数空间中按照高斯随机游走的方式生成一个新的提案点，这里的高斯随机游走简单地是将当前参数加上一个零均值、方差为先验部分标准差乘以步长系数得到的随机扰动。这样产生的提案参数向量在数值上与当前参数向量存在一定偏差，偏差大小由步长参数决定，步长过小会导致采样链缓慢探索参数空间；步长过大则可能导致提案大幅偏离先验区间，反而被拒绝概率增加影响采样效率。为了保证物理意义，系统对提案参数设定了硬性边界约束，比如阻尼系数必须为正且在合理范围内，电机效率必须在零到一之间等。如果提案参数落在这些边界之外，则直接拒绝此次提案并保持当前参数不变，将当前参数直接记录到采样链里，同时误差值也保持不变，进入下一次迭代。

如果提案参数在允许范围内，则系统会重新调用simulate_system函数，通过ODE求解器在给定参数下对系统进行完整仿真，并基于仿真输出与参考数据计算新的加权均方误差，进而计算对数似然与对数先验，再相加得到提案的对数后验值。下一步便是计算接受率的核心逻辑：以对数后验值差为基准，如果提案对数后验大于当前值，则无条件接受提案；如果提案对数后验不大于当前，则以二者差值的指数值作为接受概率，与一个[0,1]之间的随机数进行比较，如果随机数小于该概率，则仍然接受提案，否则拒绝提案。通过这样的接受/拒绝机制，采样链能够在保证遍历参数空间同时向高后验区域集中，使得最终保留下来的样本集能够近似地服从后验分布。

在采样循环过程中，为了让用户直观了解采样进展与模型误差变化，系统在初始化阶段就创建了两个实时绘图窗口：一个用于呈现随迭代次数变化的误差曲线，另一个用于呈现各参数随迭代次数的轨迹曲线。在每次迭代结束时，系统都会将更新后的误差值以及当前链中所有样本的各参数取值传递给绘图句柄，MATLAB绘图引擎就会更新曲线，使得用户能够看到误差曲线如何逐步下降以及参数轨迹如何在目标区间附近上下波动并最终收敛。为了避免因频繁绘图而导致性能瓶颈，系统使用了“限制绘图更新频率”的方式（调整drawnow的调用方式），确保在每次迭代中只有一个绘图动作而非针对每个参数多次更新。这样既能保证可视化的实时性，也能最大程度减少绘图开销对仿真速度产生的影响。

当全部迭代结束后，系统会进入后处理阶段。首先，根据用户预设的弃置期长度，将采样链中前期未收敛的部分舍弃，仅保留后半段或其他设定部分作为有效样本集合。这些样本被认为已经在高后验区域进行充分采样，能够真实反映参数后验分布的形状。然后，系统使用这些有效样本计算各参数的均值，将其作为最终的参数点估计结果并打印到命令行或对话框中。此外，为了让用户直观了解估计结果的分布特性，系统还绘制了每个参数的后验分布直方图。具体做法是在同一图形窗口中，以若干子图的形式分别绘制各参数的直方图，并在每个子图上标出先验分布的中心位置以及后验样本均值所在位置，用不同颜色的竖线加以区分。通过这样的可视化呈现，用户可以一方面直观对比先验与后验的差异，另一方面观察后验分布的集中程度与模式，比如是否呈单峰分布或多峰分布，以及标准差有多小等。

除了后验直方图之外，系统还计算了样本集的相关矩阵，用以衡量不同参数之间在后验分布下的线性相关程度。相关矩阵的每个元素表示两个参数的样本联合变化趋势是否正相关、负相关或基本不相关，这对后续要做的模型简化或联合辨识有参考价值。例如，如果两个参数高度正相关，说明在数据提供的信息下，它们不能同时独立估计，可能需要重新考虑参数化方式或二次设计实验。系统会将相关矩阵以表格形式打印出来，让用户清晰地看到各参数之间的相关性系数。

在进行完样本统计分析之后，系统还会基于最优参数点估计（即后验均值）与参考数据重新进行一次仿真，以便绘制“真实摆角轨迹”与“估计摆角轨迹”的对比图。该图同样使用时间向量作为横坐标，以角度（以度为单位）作为纵坐标，通过不同颜色或线型线条展现真实与估计之间的重叠情况。如果估计结果较好，两条曲线将基本重合；如果存在偏差，则能够从该图形上直观地观察出建模误差或参数估计的不准确性，以便用户进一步调整模型结构或采样设置。

最后，系统基于保留下来的后验样本计算赤池信息量准则（AIC），这是一种在给定数据下衡量模型优劣的综合指标，既考虑了模型拟合误差大小（用对数似然或误差度量体现），也对模型复杂度（参数维数）给予惩罚，从而帮助用户在多个候选模型之间进行比较。AIC的计算方法是将模型参数个数乘以二，再加上数据点数乘以误差对数值，并进行适当缩放。系统会将最终的AIC数值打印出来，以便用户判断模型是否过于复杂或是否需要简化。

总体而言，该系统将旋转倒立摆动力学建模、数值积分、误差度量、贝叶斯框架、MCMC采样、实时可视化和后验分布分析多个技术环节无缝衔接，使得用户只需提供基础的初始状态、时间序列和先验信息，就能自动完成从数据驱动到参数估计再到结果可视化的闭环工作流程。系统中的每一项技术都各司其职：数值积分保证了对非线性动力学方程的准确仿真；加权均方误差使得误差度量更贴合物理需求；高斯先验与对数变换为后验构造提供了理论依据；Metropolis-Hastings算法确保了在高维参数空间中的随机游走能够向高后验区域聚集并获得收敛分布；实时绘图让用户随时掌握迭代进度与模型收敛情况；后验直方图和相关矩阵则为参数不确定性与相互关系提供了完整的统计视角，而AIC则帮助用户做出模型选择决策。

2、仿真结果演示

3、关键代码展示

略

4、MATLAB 源码获取

  V

点击下方名片关注公众号获取