💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
💥1 概述
基于脉冲响应约束的质量-弹簧-阻尼系统子空间辨识研究
摘要：当只有少量嘈杂的测量数据可用时，子空间辨识方法可能会产生不可靠的模型估计。在这种情况下，通过利用关于系统的先验知识可以提高估计参数的准确性。本文考虑的先验知识是对脉冲响应的约束。这是受到有关系统稳态增益、超调和上升时间的信息的影响，这些信息可以表达为对脉冲响应的约束。所提出的方法包括两个步骤：(1) 带有线性等式和不等式约束的脉冲响应估计，以及 (2) 估计脉冲响应的实现。第一步的问题被证明是一个凸二次规划问题。在将先验知识表达为等式约束的情况下，第一步的问题有一个封闭形式的解决方案。在一般情况下，包括等式和不等式约束，解决方案通过标准的数值优化方法计算。我们通过一个质量-弹簧-阻尼系统来说明该方法的性能。
本文的主要目标是在用户对待识别的过程具有先验信息时，改进标准子空间算法的效率。这些信息可能来自控制系统的物理定律、初步实验（如阶跃响应或正弦输入信号的响应），或专家知识。例如，用户可能了解系统的稳态增益、调节时间或主导时间常数。开发的识别方法是一种通用建模工具，并不局限于特定的应用领域。实际上，任何应用都可以从在识别过程中利用先验知识中受益，只要（1）存在这样的先验知识，以及（2）存在可以利用它的方法。
对于系统的先验知识可以自然地表达为其行为的约束，例如，超调和上升时间是用阶跃响应来定义的（Merćère, 2016）。然而，在参数辨识中，模型由参数向量表示 - 传递函数或状态空间表示的系数。然后，识别问题就变成了参数估计问题，而先验知识的包含需要将其重新表述为对参数向量的约束。这可能是不平凡的，并导致更复杂的优化问题。事实上，对系统行为的线性约束通常会导致对参数向量的非线性约束（Rothenberg, 1973）。
当我们处理子空间辨识时，很难直接将这样的先验知识引入模型结构中。子空间辨识不依赖于显式的成本函数，并使用系统的状态空间表示，该表示已知直到相似变换。因此，将具有物理意义的先验信息引入到由子空间辨识算法估计的状态空间模型中似乎是一个具有挑战性的问题。
在本文中，我们通过以下两步方法绕过了与将先验知识包含在参数估计中相关的困难：
（1）脉冲响应的估计，以及
（2）估计脉冲响应的实现。
在第一步中，先验知识被强加在估计的脉冲响应上。该方法基于Markovsky和Rapisarda（2008）的一个结果，其中它表明，对于精确的数据，线性时不变系统的脉冲响应可以通过解一个超定的线性方程组直接从数据中计算出来。
本文针对质量-弹簧-阻尼系统（MSDS）在噪声数据下的参数估计问题，提出了一种基于脉冲响应约束的子空间辨识方法。该方法通过两步策略：首先利用脉冲响应约束进行估计，再实现其估计过程，有效解决了传统子空间辨识在噪声环境下的优化复杂性。结合稳态增益和响应特性等先验知识，显著提升了模型估计的可靠性。通过质量-弹簧-阻尼系统的实例验证，证明了该方法在有限嘈杂数据下仍能保持较高的参数估计精度，为复杂系统的建模提供了新的理论工具。
关键词：质量-弹簧-阻尼系统；子空间辨识；脉冲响应约束；先验知识；参数估计
1 引言
质量-弹簧-阻尼系统（Mass-Spring-Damper System, MSDS）作为经典力学模型，广泛应用于机械振动、控制理论及工程抗震等领域。其动态特性由质量、弹簧刚度和阻尼系数共同决定，传统建模方法通常基于整数阶微分方程。然而，实际系统中非理想阻尼特性（如材料黏弹性、流体非牛顿性）常导致整数阶模型精度不足。分数阶微积分虽能更精确描述此类特性，但需引入分数阶导数，增加了模型复杂度。
子空间辨识方法通过状态空间模型直接从输入输出数据中估计系统参数，无需显式成本函数，具有计算效率高的优势。但在噪声数据下，传统子空间辨识易产生不可靠估计。为解决这一问题，本文提出一种基于脉冲响应约束的子空间辨识方法，通过引入稳态增益、超调量等先验知识，将物理约束转化为对脉冲响应的数学约束，从而提升模型估计的鲁棒性。
2 理论基础
2.1 质量-弹簧-阻尼系统动力学模型
MSDS的经典动力学方程为：
编辑
- 欠阻尼（ζ<1）：系统发生衰减振荡，振幅逐渐减小至零。
- 临界阻尼（ζ=1）：系统以最快速度回归平衡位置，无振荡。
- 过阻尼（ζ>1）：系统缓慢回归平衡位置，无振荡。
2.2 子空间辨识方法
子空间辨识通过投影算法从输入输出数据中估计状态空间模型，其核心步骤包括：
- 数据矩阵构造：将输入输出数据排列为Hankel矩阵。
- 投影计算：利用QR分解或SVD分解提取系统阶数和状态序列。
- 参数估计：通过最小二乘法或总体最小二乘法估计状态空间矩阵。
传统子空间辨识未考虑先验知识，在噪声数据下易产生过拟合或欠拟合问题。
2.3 脉冲响应约束的引入
脉冲响应是系统对单位脉冲输入的输出响应，其系数直接反映系统动态特性。通过将先验知识（如稳态增益、超调量、调节时间）转化为对脉冲响应的数学约束，可显著提升估计精度。例如：
- 稳态增益约束：脉冲响应的稳态值等于系统直流增益，即 limt→∞h(t)=G(0)。
- 超调量约束：脉冲响应的最大值与稳态值的比值不超过给定阈值。
- 调节时间约束：脉冲响应进入稳态值±5%范围的时间不超过给定值。
3 基于脉冲响应约束的子空间辨识方法
3.1 方法框架
本文提出的方法分为两步：
- 脉冲响应估计：在先验知识约束下，估计系统的脉冲响应系数。
- 状态空间实现：将估计的脉冲响应转换为状态空间模型。
3.1.1 脉冲响应估计
脉冲响应估计问题可表述为：
编辑
其中，Y为输出数据向量，Φ为输入脉冲响应矩阵，h为脉冲响应系数向量，A和 b为约束矩阵和向量。当约束为线性等式时，问题为凸二次规划，存在封闭解；当包含不等式约束时，需采用数值优化方法（如内点法）求解。
3.1.2 状态空间实现
将估计的脉冲响应 h转换为状态空间模型 (A,B,C,D) 的过程称为实现问题。通过Markovsky和Rapisarda的结果，精确数据下脉冲响应可通过解超定线性方程组直接计算。本文采用基于Ho-Kalman算法的实现方法，确保状态空间模型的最小实现性。
3.2 先验知识的数学表达
先验知识的引入是提升估计精度的关键。本文考虑以下三类约束：
-
编辑
-
编辑
4 实例验证
4.1 系统参数与数据生成
考虑一个质量-弹簧-阻尼系统，参数为 m=1kg，c=0.5Ns/m，k=1N/m。系统阻尼比 ζ=0.25，属于欠阻尼系统。输入信号为高斯白噪声，输出信号添加信噪比为20 dB的高斯噪声。
4.2 辨识结果对比
分别采用传统子空间辨识（N4SID）和本文方法进行辨识，结果如下：
- 脉冲响应估计：
- 传统方法：脉冲响应系数估计误差较大，尤其在高频段。
- 本文方法：通过稳态增益约束和超调量约束，脉冲响应系数估计误差显著降低。
- 状态空间模型：
- 传统方法：状态空间矩阵 A 的特征值偏离真实值，导致动态响应失真。
- 本文方法：状态空间矩阵 A 的特征值与真实值接近，动态响应准确。
- 频率响应：
- 传统方法：频率响应在共振峰处偏差较大。
- 本文方法：频率响应与真实系统高度吻合，尤其在低频段。
4.3 鲁棒性分析
为验证方法的鲁棒性，分别在信噪比为15 dB和10 dB的噪声环境下进行辨识。结果表明，本文方法在低信噪比下仍能保持较高的估计精度，而传统方法性能显著下降。
5 结论
本文提出了一种基于脉冲响应约束的质量-弹簧-阻尼系统子空间辨识方法，通过引入稳态增益、超调量等先验知识，有效提升了模型估计的可靠性。实例验证表明，该方法在噪声数据下仍能保持较高的参数估计精度，为复杂系统的建模提供了新的理论工具。未来研究可进一步探索非线性系统和多输入多输出系统的脉冲响应约束辨识方法。
📚2 运行结果
编辑
编辑
编辑
编辑
编辑
编辑
编辑
编辑
可视化代码：
figure, hold on,
plot(Esys(:, 1), 'b.', 'MarkerSize', 30), h(1) = plot(Esys(:, 1), 'b--', 'linewidth', 2);
plot(Esys(:, 2), 'r.', 'MarkerSize', 30), h(2) = plot(Esys(:, 2), 'r-', 'linewidth', 2);
plot(Esys(:, 3), 'r.', 'MarkerSize', 30), h(3) = plot(Esys(:, 3), 'r-.', 'linewidth', 2);
title(''), xlabel('number of equality constraints'), ylabel('esys')
legend(h, {'uy2ss\_pk', 'uy2ss', 'n4sid'}, 'location', 'SouthWest'),
ax = axis; set(gca, 'xTick', 1:10), set(gca, 'yTick', [0 ax(4)])
ax = axis; axis([1 i ax(3:4)])
set(gca, 'fontsize', 15)
figure, hold on,
plot(Ehh(:, 1), 'b.', 'MarkerSize', 30), h(1) = plot(Ehh(:, 1), 'b--', 'linewidth', 2);
plot(Ehh(:, 2), 'r.', 'MarkerSize', 30), h(2) = plot(Ehh(:, 2), 'r-', 'linewidth', 2);
title(''), xlabel('number of equality constraints'), ylabel('ehh')
legend(h(1:2), {'uy2ss\_pk', 'uy2ss'}, 'location', 'SouthWest'),
ax = axis; set(gca, 'xTick', 1:10), set(gca, 'yTick', [0 ax(4)])
ax = axis; axis([1 i ax(3:4)])
set(gca, 'fontsize', 15)
figure, hold on,
plot(Ehhp(:, 1), 'b.', 'MarkerSize', 30), h(1) = plot(Ehhp(:, 1), 'b--', 'linewidth', 2);
plot(Ehhp(:, 2), 'r.', 'MarkerSize', 30), h(2) = plot(Ehhp(:, 2), 'r-', 'linewidth', 2);
title(''), xlabel('number of equality constraints'), ylabel('ehhp')
legend(h(1:2), {'uy2ss\_pk', 'uy2ss'}, 'location', 'NorthWest'),
ax = axis; set(gca, 'xTick', 1:10), set(gca, 'yTick', [0 ax(4)])
ax = axis; axis([1 i ax(3:4)])
set(gca, 'fontsize', 15)
figure, hold on,
plot(Ehp(:, 1), 'b.', 'MarkerSize', 30), h(1) = plot(Ehp(:, 1), 'b--', 'linewidth', 2);
plot(Ehp(:, 2), 'r.', 'MarkerSize', 30), h(2) = plot(Ehp(:, 2), 'r-', 'linewidth', 2);
plot(Ehp(:, 3), 'r.', 'MarkerSize', 30), h(3) = plot(Ehp(:, 3), 'r-.', 'linewidth', 2);
title(''), xlabel('number of equality constraints'), ylabel('ehp')
legend(h, {'uy2ss\_pk', 'uy2ss', 'n4sid'}, 'location', 'NorthWest'),
ax = axis; set(gca, 'xTick', 1:10), set(gca, 'yTick', [0 ax(4)])
ax = axis; axis([1 i ax(3:4)])
set(gca, 'fontsize', 15)
🎉3 参考文献
文章中一些内容引自网络，会注明出处或
