【质量-弹簧-阻尼系统】基于脉冲响应约束的子空间辨识研究（Matlab代码实现）

 💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码、数据、文章

---

💥1 概述

文献来源：

摘要
当只有少量嘈杂的测量数据可用时，子空间辨识方法可能会产生不可靠的模型估计。在这种情况下，通过利用关于系统的先验知识可以提高估计参数的准确性。本文考虑的先验知识是对脉冲响应的约束。这是受到有关系统稳态增益、超调和上升时间的信息的影响，这些信息可以表达为对脉冲响应的约束。所提出的方法包括两个步骤：(1) 带有线性等式和不等式约束的脉冲响应估计，以及 (2) 估计脉冲响应的实现。第一步的问题被证明是一个凸二次规划问题。在将先验知识表达为等式约束的情况下，第一步的问题有一个封闭形式的解决方案。在一般情况下，包括等式和不等式约束，解决方案通过标准的数值优化方法计算。我们通过一个质量-弹簧-阻尼系统来说明该方法的性能。

本文的主要目标是在用户对待识别的过程具有先验信息时，改进标准子空间算法的效率。这些信息可能来自控制系统的物理定律、初步实验（如阶跃响应或正弦输入信号的响应），或专家知识。例如，用户可能了解系统的稳态增益、调节时间或主导时间常数。开发的识别方法是一种通用建模工具，并不局限于特定的应用领域。实际上，任何应用都可以从在识别过程中利用先验知识中受益，只要（1）存在这样的先验知识，以及（2）存在可以利用它的方法。

对于系统的先验知识可以自然地表达为其行为的约束，例如，超调和上升时间是用阶跃响应来定义的（Merćère, 2016）。然而，在参数辨识中，模型由参数向量表示 - 传递函数或状态空间表示的系数。然后，识别问题就变成了参数估计问题，而先验知识的包含需要将其重新表述为对参数向量的约束。这可能是不平凡的，并导致更复杂的优化问题。事实上，对系统行为的线性约束通常会导致对参数向量的非线性约束（Rothenberg, 1973）。

当我们处理子空间辨识时，很难直接将这样的先验知识引入模型结构中。子空间辨识不依赖于显式的成本函数，并使用系统的状态空间表示，该表示已知直到相似变换。因此，将具有物理意义的先验信息引入到由子空间辨识算法估计的状态空间模型中似乎是一个具有挑战性的问题。

在本文中，我们通过以下两步方法绕过了与将先验知识包含在参数估计中相关的困难：
（1）脉冲响应的估计，以及
（2）估计脉冲响应的实现。
在第一步中，先验知识被强加在估计的脉冲响应上。该方法基于Markovsky和Rapisarda（2008）的一个结果，其中它表明，对于精确的数据，线性时不变系统的脉冲响应可以通过解一个超定的线性方程组直接从数据中计算出来。

📚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 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

🌈4 Matlab代码、数据、文章