(完整版)一级倒立摆系统分析

一级倒立摆的系统分析一、倒立摆系统的模型建立如图1-1所示为一级倒立摆的物理模型图1-1一级倒立摆物理模型对于上图的物理模型我们做以下假设：M：小车质量m：摆杆质量b：小车摩擦系数l：摆杆转动轴心到杆质心的长度I：摆杆惯量F：加在小车上的力x：小车位置ɸ：摆杆与垂直向上方向的夹角θ：摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下）图1-2是系统中小车和摆杆的受力分析图。其中，N和P为小车与摆小车导轨xlF摆杆ɸɸ杆相互作用力的水平和垂直方向的分量。注意：实际倒立摆系统中的检测和执行装置的正负方向已经完全确定，因而矢量方向定义如图所示，图示方向为矢量正方向。图1-2小车及摆杆受力分析分析小车水平方向受力，可以得到以下方程：Mẍ＝Ｆ-ｂẋ-Ｎ(1-1)由摆杆水平方向的受力进行分析可以得到以下方程：N=md2dt2(x+lsinθ)(1-2)即：N=mẍ+mlθ̈cosθ−mlθ̇2sinθ(1-3)将这个等式代入式（1-1）中，可以得到系统的第一个运动方程：(M+m)ẍ+bẋ+mlθ̈cosθ−mlθ̇2sinθ=F(1-4)为推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得出以下方程：P−mg=md2dt2(lcosθ)(1-5)P−mg=−mlθ̈sinθ−mlθ̇2cosθ(1-6)利用力矩平衡方程可以有：mgNPθIθ’’FNPMXBx’X’’−Plsinθ−Nlcosθ=Iθ̈(1-7)注意：此方程中的力矩方向，由于θ=π+ɸ，cosɸ=−cosθ，sinɸ=−sinθ，所以等式前面含有负号。合并两个方程，约去P和N可以得到第二个运动方程：(I+ml2)θ̈+mglsinθ=−mlẍcosθ(1-8)设θ=π+ɸ，假设ɸ与1（单位是弧度）相比很小，即ɸ＜＜1，则可以进行近似处理：cosθ=−1，sinθ=−ɸ，（dθdt）2=0。用u来代表被控对象的输入力F，线性化后的两个运动方程如下：{(I+ml2)ɸ̈−mglɸ=mlẍ(M+m)ẍ+bẋ−mlɸ̈=u(1-9)假设初始条件为0，则对式（1-9）进行拉普拉斯变换，可以得到：{(I+ml2)Φ(s)s2−mglΦ(s)=mlX(s)s2(M+m)X(s)s2+bX(s)s−mlΦ(s)s2=U(s)(1-10)由于输出为角度ɸ，求解方程组的第一个方程，可以得到：X(s)=[(I+ml2)ml−gs2]Φ(s)(1-11)或改写为：Φ(s)X(s)=mls2(I+ml2)s2−mgl(1-12)如果令v=ẍ,则有：Φ(s)V（s）=ml(I+ml2)s2−mgl(1-13)如果将上式代入方程组的第二个方程，可以得到：(M+m)[(I+ml2)ml−gs]Φ(s)s2+b[(I+ml2)ml+gs2]Φ(s)s−mlΦ(s)s2=U(s)(1-14)整理后可得传递函数：Φ(s)U(s)=mlqs2s4+b(I+ml2)qs3−(M+m)mglqs2−bmglqs(1-15)其中q=[(M+m)(I+ml2)−(ml)2]假设系统状态空间方程为：Ẋ=AX+Buy=CX+Du(1-16)方程组对ẍ,ɸ̈解代数方程，可以得到解如下：{ẋ=ẋẍ=−(I+ml2)bI(M+m)+Mml2ẋ+m2gl2I(M+m)+Mml2ɸ+(I+ml2)I(M+m)+Mml2uɸ̇=ɸ̇ɸ̈=−mlbI(M+m)+Mml2ẋ+mgl(M+m)I(M+m)+Mml2ɸ+mlI(M+m)+Mml2u(1-17)整理后可以得到系统状态空间方程：[ẋẍɸ̇ɸ̈]=[01000−(I+ml2)bI(M+m)+Mml2m2gl2I(M+m)+Mml2000010−mlbI(M+m)+Mml2mgl(M+m)I(M+m)+Mml20][xẋɸɸ̇]+[0(I+ml2)I(M+m)+Mml20mlI(M+m)+Mml2]uy=[xɸ]=[10000010][xẋɸɸ̇]+[00]u(1-18)由（1-9）的第一个方程为：(I+ml2)ɸ̈−mglɸ=mlẍ对于质量均匀分布的摆杆可以有：I=13ml2于是可以得到：(13ml2+ml2)ɸ̈−mglɸ=mlẍ化简可以得到：ɸ̈=3g4lɸ+34lẍ(1-19)设X={x,ẋ,ɸ,ɸ̇}，u̇=ẍ则有：[ẋẍɸ̇ɸ̈]=[010000000001003g4l0][xẋɸɸ̇]+[01034l]u̇y=[xɸ]=[10000010][xẋɸɸ̇]+[00]u̇(1-20)以上公式推理是根据牛顿力学的微分方程验证的。在实际系统中模型参数如下：M小车质量1.096Kgm摆杆质量0.109Kgb小车摩擦系数0.1N/m/secl摆杆转动轴心到杆质心的长度0.25mI摆杆惯量0.0034kg*m*m将上述参数代入，就可以得到系统的实际模型。摆杆角度和小车位移的传递函数：Φ(s)X(s)=0.02725s20.0102125s2−0.26705(1-21)摆杆角度和小车加速度之间的传递函数为：Φ(s)V（s）=0.027250.0102125s2−0.26705(1-22)摆杆角度和小车所受外界作用力的传递函数：Φ(s)U（s）=2.35655ss3+0.0883167s2−27.9169s−2.30942(1-23)以外界作用力作为输入的系统状态方程：[ẋẍɸ̇ɸ̈]=[01000−0.08831670.629317000010−0.23565527.82850][xẋɸɸ̇]+[00.88316702.35655]uy=[xɸ]=[10000010][xẋɸɸ̇]+[00]u(1-24)以小车加速度作为输入的系统状态方程：[ẋẍɸ̇ɸ̈]=[0100000000010029.40][xẋɸɸ̇]+[0103]u̇y=[xɸ]=[10000010][xẋɸɸ̇]+[00]u̇（1-25）综述可知以上就是一级倒立摆系统的模型建立过程，最终得出了实际模型的传递函数和状态空间方程。二、系统模型的转换以小车加速度作为输入的系统状态方程为例，将系统状态方程转化为能控标准型，能观标准型和约当标准型。由系统状态方程可知：A=[0100000000010029.40]B=[0103]C=[10000010]D=[00]1、转化为能控标准型定出系统特征多项式：a=poly(A)a=1.0000-0.0000-29.400000由此可知a0=0,a1=0,a2=-29.4,a3=0。b3=C*Bb3=00b2=C*A*B+a3*C*Bb2=13b1=C*A^2*B+a3*C*A*B+a2*C*Bb1=00b0=C*A^3*B+a3*C*A^2*B+a2*C*A*B+a1*C*Bb0=-29.40000所以系统的能控标准型为：[ẋ1ẍ1ɸ1̇ɸ1̈]=[0100001000010029.40][x1x1̇ɸ1ɸ1̇]+[0001]u̇y=[−29.40100030][x1x1̇ɸ1ɸ1̇]+[00]u̇2、转化为能观标准型利用对偶性求出能观标准型为：[ẋ1ẍ1ɸ1̇ɸ1̈]=[0000100001029.40010][x1x1̇ɸ1ɸ1̇]+[0−29.4003300]u̇y=[0001][x1x1̇ɸ1ɸ1̇]3、转化为约当标准型首先求出系统的特征值以及相应的特征向量：A=[0100;0000;0001;0029.40]A=01.00000000000001.00000029.40000[V,D]=eig(A)V=001.0000-1.00000000.00000.1814-0.1814000.98340.983400D=5.42220000-5.42220000000000其中D表示A全部特征值构成的对角阵，V表示相对应的特征向量。求出变换矩阵V的逆：V1=inv(V)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=1.720635e-292.V1=1.0e+291*000.00000.000000-0.00000.00000.00002.49480002.494800计算变换后的系数矩阵：A1=V1*A*VA1=5.42220000.0000-5.4222000000.00000000B1=V1*BB1=1.0e+291*0.00000.00002.49482.4948所以系统的约当标准型为：[ẋ1ẍ1ɸ1̇ɸ1̈]=[5.42220000−5.42220000000000][x1x1̇ɸ1ɸ1̇]+1.0e+291∗[002.49482.4948]u̇三、开环阶跃响应曲线及分析利用已知的状态空间方程来进行阶跃响应分析，在MATLAB中可以写入以下命令：A=[0100;0000;0001;0029.40]；B=[0;1;0;3]；C=[1000;0100]；D=[0;0]；step(A,B,C,D)可以看出，在单位阶跃响应作用下，小车位置和摆杆角度都是发散的。四、判断系统稳定性判断系统的稳定性可以利用根轨迹来判断，已知实际系统的开环传递函数为：Φ(s)V（s）=0.027250.0102125s2−0.26705，则其根轨迹图形可以利用MATLAB键入如下命令来完成。num=[0.02725];den=[0.01021250-0.26705];z=roots(num)z=Emptymatrix:0-by-1p=roots(den)p=5.1136-5.1136rlocus(num,den)可以看出系统没有零点，有两个极点，并且有一个极点为正。由画出的根轨迹图形可以看出闭环传递函数的一个极点位于复平面的右半平面，这就意味着系统是一个不稳定的系统。五、能控性和能观性分析对于系统的能控性和能观性分析，可以利用能控性秩判据和能观性秩判据。能控性秩判据：对于n维连续时间线性时不变系统，构成能控性判别矩阵：Qc=[B⋮AB⋮⋯⋮An−1B]，则系统完全能控的充要条件为：rankQc=rank[B⋮AB⋮⋯⋮An−1B]=n能观性秩判据：对于n维连续时间线性时不变系统，构成能观性判别矩阵：Qo=[CCA⋮CAn−1]，则系统完全能观的充要条件为：rankQo=rank[CCA⋮CAn−1]=n利用MATLAB键入以下命令来进行判断：A=[0100;0000;0001;0029.40];B=[0;1;0;3];C=[1000;0100];D=[0;0];Qc=[BA*BA^2*BA^3*B]Qc=01.0000001.000000003.0000088.20003.0000088.20000R1=rank(Qc)R1=4Qo=[C;C*A;C*A^2;C*A^3]Qo=10000100010000000000000000000000R2=rank(Qo)R2=2可以看出，系统的完全能控矩阵的秩等于系统的状态变量维数，系统的输出完全能观测矩阵的秩等于系统输出向量y的维数，所以系统是可以完全能控完全能观测的系统。六、根轨迹校正以及仿真已知系统的传递函数：G(s)=0.027250.0102125s2−0.26705设计控制器使得调整时间ts=0.5s(2%)；最大超调Mp≤10%。计算整理可得超前校正装置的零点和极点分别为：zc=−6.92214；zp=−26.4568，由此可得校正后的传递函数：Q=G(s)K(s)=K(s+6.92214)s+26.45680.027250.0102125s2−0.26705利用MATLAB命令观察校正后的根轨迹图形：图中可以看出，系统三条根轨迹都位于左半平面，只需要选定适当的K就可以使得系统稳定。现在利用MATLAN仿真，选择不同的K观察系统的阶跃响应曲线。K=141.137时，画出闭环控制系统：响应曲线图：K=120.115时，响应曲线：K=50