弹簧-质量-阻尼模型

来源：榕意旅游网

弹簧-质量-阻尼系统

1 研究背景及意义

弹簧-质量-阻尼系统是一种比较普遍的机械振动系统，研究这种系统对于我们的生活与科技也是具有意义的，生活中也随处可见这种系统，例如汽车缓冲器就是一种可以耗减运动能量的装置，是保证驾驶员行车安全的必备装置，再者在建筑抗震加固措施中引入阻尼器，改变结构的自振特性，增加结构阻尼，吸收地震能量，降低地震作用对建筑物的影响。因此研究弹簧-质量-阻尼结构是很具有现实意义。

2 弹簧-质量-阻尼模型的建立

数学模型是定量地描述系统的动态特性，揭示系统的结构、参数与动态特性之间关系的数学表达式。其中，微分方程是基本的数学模型，不论是机械的、液压的、电气的或热力学的系统等都可以用微分方程来描述。微分方程的解就是系统在输入作用下的输出响应。所以，建立数学模型是研究系统、预测其动态响应的前提。通常情况下，列写机械振动系统的微分方程都是应用力学中的牛顿定律、质量守恒定律等。

弹簧-质量-阻尼系统是最常见的机械振动系统。机械系统如图2.1所示，

图2.1 弹簧-质量-阻尼系统简图

其中m1，m2表示小车的质量，ci表示缓冲器的粘滞摩擦系数，ki表示弹簧的弹性系数，Fi（t）表示小车所受的外力，是系统的输入即Ui（t）=Fi（t）,Xi(t)表示小车的位移，是系统的输出，即Yi（t）=Xi(t)，i=1,2。设缓冲器的摩擦力与活塞的速度成正比，其中

m1=1kg，m2=2kg，k1=k3=100N/cm，k2=300N/cm，c1=c3=3N•s/cm，c2=6N•s/cm。

由图2.1，根据牛顿第二定律,，建立系统的动力学模型如下：对m1有：

（2-1）

对m2有：

（2-2）

3 建立状态空间表达式

令x3x&1,x4x&2,u1F1,u2F2,则原式可化为：

m1x&3(l1l2)x3l2x4(k1k2)x1k2x2u1(t)m2x&4(l2l3)x4l2x

3(k2k3)x2k2x1u2(t)化简得：

x&u1(t)k2x2(k1k2)x1(l1l2)x3l2x43m 1

x&2(t)k2x1(k1k2)x2(l3l2)x4l2x34um 2

整理得：

0010000001x100x&(k1k2)k2(l1l2)l20m1m1m1mx211x3m1u1u2 k2(k3k2)l2(l)x3l241m2m2m2m02m 2x1y100000x201x3x4

m11,m22,k1k3100,k2300l1l33,l26

2-3）

2-4）

（2-5）（（ 0000代入数据得：A400300150200则系统的状态空间表达式为

100000011000 B C 1096010034.500.5000.0x400300200151000yx010010000001xu9610

34.500.54 化为对角标准型

当系统矩阵A有n个不相等的特征根

(i1,2,3...)时，相应的有n个不相等的特征向量

m(i1,2,3...)，所以有矩阵A的特征矩阵Mmmm...m

ii1234根据矩阵论

线性变换得：TM1zTxxMz

可以使用matlab进行对角标准型的运算，matlab作为一种数算工具，很大程度的方便了了我们的计算，对于这个弹簧-质量-阻尼系统是一个四阶的状态空间表达式，所以可以用matlab简化计算。

（1）求特征值与特征向量

A=[0 0 1 0;0 0 0 1;-400 300 9 6;150 -200 3 -4.5] B=[0 0;0 0;1 0;0 0.5] C=[1 0 0 0;0 1 0 0]

[P,J]=eig(A) 求得结果：

P =

0.0007 - 0.0402i 0.0007 + 0.0402i 0.0401 - 0.0698i 0.0401 + 0.0698i -0.0171 + 0.0157i -0.0171 - 0.0157i 0.0176 - 0.0792i 0.0176 + 0.0792i 0.8650 0.8650 0.6682 + 0.2084i 0.6682 - 0.2084i

-0.3442 - 0.3621i -0.3442 + 0.3621i 0.7050 0.7050 J =

0.3667 +21.5183i 0 0 0 0 0.3667 -21.5183i 0 0 0 0 1.8833 + 8.48i 0 0 0 0 1.8833 - 8.48i （2）P矩阵求逆 PN=inv(P)

求得结果： PN =

3.4167 + 9.7803i -2.1017 - 9.2399i 0.3466 - 0.2323i -0.4703 - 0.10i 3.4167 - 9.7803i -2.1017 + 9.2399i 0.3466 + 0.2323i -0.4703 + 0.10i -3.35 + 3.4224i 3.7199 + 3.2032i 0.2886 - 0.0353i 0.5337 - 0.2409i -3.35 - 3.4224i 3.7199 - 3.2032i 0.2886 + 0.0353i 0.5337 + 0.2409i （3）带入公式BPNB CCP 解得对角标准型为：

0000.3667 21.5183i 0.3466 - 0.2323i-0.2352 - 0.0527i0.3466 + 0.2323i-0.2352 + 0.0527i00.3667 21.5183i00xux0.2886 - 0.0353i0.2669 - 0.1205i00 1.8833  8.48i00001.8833 - 8.48i 0.2886 + 0.0353i0.2669 + 0.1205i

- 0.0698i0.0401 + 0.0698i0.0007 - 0.0402i0.0007 + 0.0402i0.0401yu  + 0.0157i-0.0171 - 0.0157i0.0176 - 0.0792i0.0176 + 0.0792i-0.01715求状态空间表达式的解（1）求状态转移矩阵

T1ATeAtTeTt1

其中，T为特征向量

- 0.0698i0.0401  0.0698ie0.0007 - 0.0402i0.0007  0.0402i0.0401-0.0171  0.0157i-0.0171 - 0.0157i 0.0176 - 0.0792i0.0176  0.0792iAt*e0.86500.86500.6682  0.2084i0.6682 - 0.2084i0.7050 0.7050 -0.3442 - 0.3621i-0.3442  0.3621i0.3667 21.5183it000000*T1000e0.3667 -21.5183it00e1.8833  8.48it0e1.8833 - 8.48it

 3.4167 + 9.7803i -2.1017 - 9.2399i 0.3466 - 0.2323i -0.4703 - 0.10i 3.4167 - 9.7803i -2.1017 + 9.2399i0.3466 + 0.2323i -0.4703 + 0.10i1 T-3.35 + 3.4224i3.7199 + 3.2032i0.2886 - 0.0353i 0.5337 - 0.2409i-3.35 - 3.4224i 3.7199 - 3.2032i 0.2886 + 0.0353i 0.5337 + 0.2409i状态转移矩阵为：

-0.77 + 0.0000i 0.1999 0.4493  -5.5977  -4.4097 - 0.0000i -0.5817 + 0.0000i 0.2247 + 0.0000i 0.5509 + 0.0000iAte -12.5772 + 0.0000i -29.8835 - 0.0000i -2.4502 -1.4700 - 0.0000i -7.2316 + 0.0000i -42.7799 - 0.0000i -0.7350 - 0.0000i -1.7127 - 0.0000i5 可控性与可观性

不同于经典控制理论，能控性和能观性，是一个具有实际意义的概念，经典控制理论中用传

递函数描述系统的输入-输出特性，输出量即被控量，只要系统是因果系统并且稳定，输出量便可以受控，且输出量总是可以被测量的，因而不需要能控能观性的提出。但是现代控制理论是建立在状态空间表达式描述系统的基础上的，状态方程描述输入u（t）引起状态x（t）

t（t）的变化。能控能观便是定性的描述的变化过程，输出方程描述有状态变化引起的输出y

x(t)eAtx(0)eA(t)Bu()d0输入u（t）对状态x(t)的控制能力，输出y（t）对状态x（t）的反应能力，他们分别回答了“输入能否控制状态的变化”------可控性

“状态的变化能否有输出反映出来”----------可观性

另外在工程上常用状态变量作为反馈信息，可是状态x（t）的值通常是难测的，往往需要从

t1At1A(t1)x(t)ex(0)eBu()d0测量到的y（t）中估计出状态，如果输出y（t）不能完全反映出系统的状态x（t），那么就10无法实现对状态的估计。

能控性定义：当系统用状态方程描述时，给定系统的任意初始状态，可以找到允许的输入

量，在有限的时间内使系统的所有状态达到任一终止状态，则称系统是完全可控的。有状态方程x’(t)=Ax(t)+Bu(t) 其解为：

如果有限的时间内0 < t < t1内通过输入量u(t)的作用把系统的所有状态引向状态x(t1)[设x(t1)=0] ，则应有：

n1 t1t1n1t1Aiix(0)eBu()dc()ABu()dABci()u()dix(0-)0即在给定和A、B的条件下求可以使x(t)=x(t1)的u(t)。换言之：上述方程有解则系统能00i0i0控。 n1t1令i(t1)＝ci()u()d  x(0)AiBi(t1) 0i0

0(t) (t)n1写成矩阵形式 x(0 )＝BABAB 1 根据凯莱－哈米尔顿定理，、 eAt可写成有限级数：e-At (t)n-1

如果方程有解，等式右边左侧矩阵应满秩=n，此时系统是可控的。

求可控性：

QBcAB(t)x0(t)1eAtx(00)eA9Bu()d0300000.532.25t AtA(t)A*ABy(t)Cex(0)C[eBu()d]31090301163.500.532.25163.580.875tn=4 满秩所以系统是可控的

可观性定义：当系统用状态方程描述时，给定控制后，如果系统的每一个初始状态x(0-)

都可以在有限的时间内通过系统的输出y(t)唯一确定，则称系统完全可观。若只能确定部分初始状态，则称系统部分可观。有状态方程x’(t)=Ax(t)+Bu(t) y(t)=Cx(t) 其解为 t1t1AtA(t)Aex(0)eBu()d0  x(0)eBu()d 00由于在讨论能观性问题时，输入是给定的，上式右侧第二项是确知的，设u(t)=0。 y(t)=CeAtx(0-)。根据凯莱－哈米尔顿定理， e-At 、eAt可写成有限级数：



如果方程有解，等式右侧中间侧矩阵应满秩。其中，秩=n(系统的阶数)

求可观性：

0101C00CAQo00CA*A400300150200n=4满秩所以系统是可观的

000010 0196366 求系统的输入输出传递函数

Y1(s)/U1(s)Y2(s)/U1(s)

对于两输入两输出的系统求得的传递函数是一个二阶的传递函数阵,其中包含四个传递函数

YY1(s)/U2(s) (s)/(s)U22Transfer function from input 1 to output...

s^2 + 4.5 s + 200 #1: --------------------------------------------

s^4 - 4.5 s^3 + 1.5 s^2 - 1800 s + 3.5e004

n1 At2n1ec0(t)Ic1(t)Ac2(t)Acn1(t)Aci(t)Ai 3 s + 150 i0 #2: --------------------------------------------

s^4 - 4.5 s^3 + 1.5 s^2 - 1800 s + 3.5e004 n12n1y(t)C[c0(t)Ic1(t)Ac2(t)Acn1(t)A]x(0)ci(t)CAix(0)

i0Transfer function from input 2 to output... 3 s + 150

#1: --------------------------------------------

C s^4 - 4.5 s^3 + 1.5 s^2 - 1800 s + 3.5e004 CAx(0)写成矩阵形式 y(t)＝c0(t)c1(t)cn1(t) 

 0.5 s^2 - 4.5 s + 200 n1CA #2: --------------------------------------------

s^4 - 4.5 s^3 + 1.5 s^2 - 1800 s + 3.5e004

矩阵函数阵：

s^2  4.5 s  200 3 s  150  --------------------------------------------  -------------------------------------------- s^4 - 4.5 s^3  1.5 s^2 - 1800 s  3.5e004 s^4 - 4.5 s^3  1.5 s^2 - 1800 s  3.5e004  3 s  150 0.5 s^2 - 4.5 s  20------------------------------------------------------------------------------- -- s^4 - 4.5 s^3  1.5 s^2 - 1800 s  3.5e004 s^4 - 4.5 s^3  1.5 s^2 - 1800 s  3.5e004

7 分析开环稳定性

稳定性定义是系统在受到小的外界扰动后，系统状态方程解的收敛性。系统正常工作要求是系统在受到外界扰动后，虽然其原平衡状态被打破，但在扰动消失后，仍然能恢复到原来的平衡状态，或者趋于另一平衡状态继续工作，且线性系统稳定性与输入作用无关。研究系统的稳定性对于研究系统能否正常工作具有很重要的意义，稳定性是自动控制系统正常工作的必要条件，是系统的重要特征。我们不仅要分析一个系统是否稳定还要解决的问题便是怎样使一个系统稳定。经典控制理论稳定性判别方法有很多，例如代数判据，niquist判据，根轨迹判据等。而现代控制理论经常用李雅普诺夫第二法求稳定性。

（1）利用特征根的方法

根据上述结果求得的特征根为 0.3667 +21.5183i 0.3667 +21.5183i 1.8833 +

8.48i 1.8833 - 8.48i ，四个特征值全部都在坐标轴的右半平面，所以系统是不稳定的。

（2）利用利亚普诺夫第二法求解

v(x)0.5(x22221x2x3x4)v(.....

x)x1x1x2x2x3x3x4x4其中：

xx1.2.3xxxx..34400x1300x29x36x4150x1200x23x34.5x4.22

4v(x)199x2x4399x1x3300x2x39x39x3x4150x1x44.5x40039915000300199可以看出A不是正定的，所以系统不稳定。将其换成矩阵形式A00900(A,B,C)k(ABK,B,C)0004.58 利用状态反馈将系统的闭环极点配置到合适的值

状态反馈是将系统每一个状态变量乘以相应的反馈系数馈送到输入端与参考输入相加，其和

作为受控系统的控制输入。

原受控对象为经过状态反馈后得到的闭环系统为闭环系统期望极点选取原则为以下几点： 1）n维控制系统有n个期望极点；

2）期望极点是物理上可实现的，为实数或共轭复数对；

3）期望极点的位置的选取，需考虑它们对系统品质的影响（离虚轴的位置），及与零点分布状况的关系。

4）离虚轴距离较近的主导极点收敛慢，对系统性能影响最大，远极点收敛快，对系统只有

极小的影响。

闭环极点 0.3667 +21.5183i 0.3667 -21.5183i 1.8833 + 8.48i 1.8833 - 8.48i 配置状态反馈后，系统应稳定，所以期望极点应在虚轴左侧，所以期望闭环极点 -1.1001+.59i -1.1001-.59i -5.99+25.4952 -5.99-25.4952 得到极点配置矩阵

K =

1.0e+003 *

1.2556 -0.0375 0.0157 -0.0332 0.9718 2.69 0.0839 0.004 验证极点配置结果是正确的：ans =

-1.1001 -.59i -1.1001 +.59i -5.99 -25.4952i -5.99 +25.4952i

求得开环传递函数阶跃响应曲线（没有经过状态反馈的）：

没有上升时间

经过状态反馈的传递函数：状态空间表达式为

G(s)C[sI(ABK)]k1B

[(ABK),B,C]

kMatlab解得闭环传递函数：

s^2 + 6.758 s + 18

#1: -------------------------------------------------- s^4 + 13.5 s^3 + 4875 s^2 + 4.86e004 s + 2.843e006

-38.95 s - 335.9

#2: -------------------------------------------------- s^4 + 13.5 s^3 + 4875 s^2 + 4.86e004 s + 2.843e006

Transfer function from input 2 to output... 19.59 s + 168.7

#1: -------------------------------------------------- s^4 + 13.5 s^3 + 4875 s^2 + 4.86e004 s + 2.843e006

0.5 s^2 + 3.371 s + 827.8

#2: -------------------------------------------------- s^4 + 13.5 s^3 + 4875 s^2 + 4.86e004 s + 2.843e006 反馈后的阶跃响应：

u0BxxCyAˆxKegBˆxCˆy

阶跃响应上升时间是0.034s，配置后系统最终稳定。 AKeC9 设计全维状态观测器

不是所有的系统状态物理上都能够直接测量得到。需要从系统的可量测参量，如输入u和输

出y来估计系统状态。状态观测器基于可直接量测的输出变量y和控制变量u来估计状态变量，是一个物理可实现的模拟动力学系统。如果系统状态是完全能观测的，那么根据输出y的测量，可以唯一的确定系统的初始状态，所以只要满足一定的条件，便可以从可测量y和u中把x间接重构出来。全维渐进状态结构图：

系统的原极点为0.3667 +21.5183i 0.3667 +21.5183i 1.8833 + 8.48i 1.8833 - 8.48i 期望极点应是原极点的2-5倍，并不是越快越好。

期望极点为 -1.4668 +86.0732i -1.4668 -86.0732i -7.5332 +33.9456i -7.5332 -33.9456i

51.117.552.65.0 L3199.91608.8122.32246.1相应的全维观测器是：

x(ALC)xBuLy00004003001502000017.51017.551.151.152.610000052.6015.05.0uyx103199.91608.8963199.91608.8010000.5122.32246.134.5122.32246.1001017.551.1000017.551.152.60052.60015.0005.00xuy4003009103199.91608.863199.91608.80015020034.5122.32246.10000.5122.32246.151.5100017.551.117.552.60052.65015.0xuy3599.91308.6103199.91608.82446.1034.500.5122.32246.127.7

. -51.13 1 0 -17.51 52.62 -4.985 0 1 10 带观测器的输出状态空间表达式：A -3600 -1309 9 63 -4.5 27.66 -2446 51.1317.51 -52.62 4.985B 3200 1609 122.3 2246

101 C000000100000 D=0

100010001

分别得到输出和观测状态的传递函数：

Transfer function from input \"y1\" to output...

17.51 s^3 + 59 s^2 + 1.23e005 s + 7.911e006 y1_e: ------------------------------------------------