三级倒立摆建模与控制器设计仿真

1、三级倒立摆摘要：从刚体运动学角度分析了三级倒立摆系统的运动，运用状态空间法对三级倒立摆系统进行数学模型的建立;应用二次型最优控理论实现三级倒立摆控制；最后根据计算结果，通过matlab程序仿真，得出上中下摆的角度和小车的位移曲线图。关键词：倒立摆；最优控；仿真；matlab三级倒立摆系统模型 一个典型的三级直线倒立摆系统主要由机电装置和控制装置两部分组成,机电装置的结构见图!小车图1系统分析示意图 其主要由摆杆、电机、小车、导轨、传动齿轮、齿条及各种传感器等构成假设系统中的每一根摆杆都是匀质刚体，驱动力与放大器的输入成正比且无延迟地直接作用于小车上，并且可以在忽略实验中的库仑摩擦和动摩擦的前提

2、下，图1所示的倒立摆系统的物理分析：yf22三级倒立摆受力分析图：f21nx1f11f12f f11f12 小车受力分析图yf31y下摆受力分析f33f21x23xf33f22f31 中摆受力分析上摆受力分析下摆牛顿定律及动力矩方程：f12-f22=m1d2dt2(r+l1sin1)m1g-f11+f21=m1d2dt2(l1cos1)j1=f11l1sin1+f21l1-l1sin1-f12l1cos1-f22(l1-l1)cos1中摆受力分析： f22-f31=m2d2dt2(r+l1sin1+l2sin2)m2g-f21+f33=m2d2dt2(l1cos1+l2cos2)j2=f21l

3、2sin2+f33l2-l2sin2-f22l2cos2-f31(l2-l2)cos2丄摆受力分析：f31=m3d2dt2(r+l1sin1+l2sin2+l3sin3)m3g-f33=m3d2dt2(l1cos1+l2cos2+l3cos3)j3=f33l3sin3-f31l3cos3 注：l1 l2 l3下，中，上摆质心到铰链处距离； l1 l2 下，中摆的长度； 1 2 3-上，中，下摆偏离竖直方向的角度； m1 m2 m3-下，中，上摆质量;其动力学方程如下：m(1 ,2 ,3)r123=f(1,2 ,3,1,2 ,3)r123+n1 ,2 ,3r123+gu m(1 ,2 ,3)=m

4、0+m1+m2+m3 m1d1cos1+m2d4cos1+m3d4cos1m2d2cos2+m3d5cos2m3d3cos3 m1d1cos1+m2d4cos1+m3d4cos1j1+m1d12+m2d42+m3d42m2d2d4cos(1-1)+m3d4d5cos(2-1)m3d3d4cos(3-1) m2d2cos2+m3d5cos2m2d2d4cos(1-1)+m3d4d5cos(2-1)j2+m2d22+m3d52m3d3d5cos(3-2)m3d3cos3m3d3d4cos(3-1)m3d3d5cos(3-2)j3+m3d32f(1,2 ,3,1,2 ,3)=-f0000 m1d1s

5、in1+m2d4sin1+m3d4sin11-f1+f2f2-m2d2d4sin(2-1)+m3d4d5sin(2-1)-m3d3d4sin(3-1) m2d2sin2+m3d5sin22f2+m2d2d4+m3d4d52sin2-1-f3+f2f3-m3d3d52sin3-2m3d3sin33m3d3d4sin(3-1)f3+m3d3d53sin3-2-f3n(1,2 ,3)=0(m1d1+m2d4+m3d4)sin1(m2d2+m3d5)gsin2m3d3gsin3g=g0000t将此非线性模型在系统的平衡点r=1=2=3=0,r= 1=2=3=0处线性化得到系统的状态方程: x=ax+b

6、uy=cx+du a=04*4i4*4a21a22 , b=04*1b2其中： 将参数带入方程计算：a=00000000000000-3.655145.74-49.782-5.159000.12698-80.216169.9-114.49 010001000.44281-8.3647-4.78172.26900-16.80957.098-62.1436.4401000.016897-0.562720.98145-0.53453 000010-0.00564470.32956-0.659640.5167301-0.0032595-0.0721790.32957-0.76437 b=00008.7

7、2-29.62132.238-3.341tc=10000100001000010000000000000000 d=0000t控制器设计根据由图所示的倒立摆系统模型的推导可知，该系统在一定的前提下，可以转变为线性系统#对于线性系统的控制，需要经过系统的稳定性、能控性和能观性判断后才能选择其最优的控制策略#图!显示了式（$）所示的线性系统外控制器的设计过程状态方程的重构过程，实质上就是本文从式式转化为式式的过程，在matlab中，可以用如下代码表示: 系统状态方程的重构 状态变量 (x1,x2,x3,x4,x5,x6,x7,x8)=(r,zeta1,zeta2,zeta3,r,zeta1,zet

8、a2,zeta3,) 输出变量y=x a=zeros(8,8); a(1:4,5:8)=eye(4,4); a(5:8,1:4)=t0m0 0 0-1n0 0 0 t0-1; a(5:8,5:8)= t0m0 0 0-1-f0 0 0 0 0 0 t0-1; 不用invese而用blackslash是由于blackslash有更好的错误检查机制并且少用时间和内存。 b=zeros(8,1); b(5:8)= t0m0 0 0-1g0 0 0 0t; c=zeros(4,8) c(1:4,1:4)=eye(4,4); d=zeros(4,1);系统的稳定性、能控性及能惯性判断 对于如式（2）式所

9、示的线性定常系统而言，其稳定性可通过计算系统矩阵a 的特征值来确定(李雅普诺夫线性系统稳定判据)，如果特征值均处于s复平面左半平面，则系统稳定。 matlab中，用函数eig（a）来计算系统矩阵的特征值，经过计算，式式所示的线性定常系统的特征值为： r1 =0 13.2715 8.3919 3.1656 -18.2950 -2.9940 -8.5547 -13.7810因此，该系统为不稳定系统，需加外部控制器形如式（2）所示的线性定常系统能控的充分必要条件是： rankb aban-1b=n;matlab中，用“ctrb(a,b)”来求系统能控矩阵pc=b aban-1b的秩，经过计算rank

10、(pc)=8,系统能控。形如式（2）所示的线性定常系统能观的充分必要条件为：rankccacan-1=nmatlab中，用“obsv(a,c)”来求系统能观阵qo=ccacan-1的秩，计算的，rank(qo)=8,系统能观。综上所述，式（2）所示的线性定常系统为不稳定的能控、能观系统，可加外控制器实现系统的稳定。线性最有调节器的设计 定义状态方程的二次性能指标函数为： j=120(xtqx+utru)dt 正定r阵的作用在于限制控制器的输出不会太大而导致难于实现，所以不失一般性，可以将其设为1，lqr方法中，q阵可以选为非负的对角阵，q=diag（ q11，q22，q33，q44，q55，q

11、66，q77，q88）.各主对角系数的相对大小对系统的整体控制性能有很大的影响，目前普遍采用试探实验的方法来确定加权矩阵q的系数大小。观察知：q11大小的变化对小车的位移影响很大，因为小车位移大小受直线导轨长度的物理限制，所以，q11取值可以固定为一个适当大小的数，使小车位移不至过大； q22对系统控制的影响最为明显，q33次之，q44的影响较小。取值时可以去q22>q33>q44,取q=diag(500,500,300,100,0,0,0,0),r=1.求解riccati代数方程 pa+atp-pbr-1btp+q=0 可以得到状态反馈阵 kt=r-1btp=-22.3607-8

12、0.103528.5577-408.4518-30.1003-53.7813-44.7459-53.7370此时各摆角的动态和稳态指标较为理想。初始条件为r0=0，10=0.05>0，20=-0.05<0，30=0.05>0。仿真结果程序：>>a=0,0,0,0,1,0,0,0;0,0,0,0,0,1,0,0;0,0,0,0,0,0,1,0;0,0,0,0,0,0,0,1;0,-3.6551,0.12698,0.44281,-16.809,0.016897,-0.0056447,-0.0032595;0,45.74,-80.216,-8.3647,57.098,-

13、0.56272,0.32956,-0.072179;0,-49.782,169.9,-4.781,-62.143,0.98145,-0.65964,0.32957;0,-5.159,-114.49,72.269,6.4401,-0.53453,0.51673,-0.76437;>>b=0;0;0;0;8.72;-29.621;32.238;-3.341;>> c=1,0,0,0,0,0,0,0;0,1,0,0,0,0,0,0;0,0,1,0,0,0,0,0;0,0,0,1,0,0,0,0;>> d=0;0;0;0;>> r1=eig(a)>

14、> q1=500;q2=500;q3=300;q4=100;q5=0;q6=0;q7=0;q8=0;>> q=q1 0 0 0 0 0 0 0;0 q2 0 0 0 0 0 0;0 0 q3 0 0 0 0 0;0 0 0 q4 0 0 0 0;0 0 0 0 q5 0 0 0;0 0 0 0 0 q6 0 0;0 0 0 0 0 0 q7 0;0 0 0 0 0 0 0 q8;r=1;>> k,p,r2=lqr(a,b,q,r)>> ac=(a-b*k);bc=b;cc=c;dc=d;>> t=0:0.02:20;u=zeros(size

15、(t);>> x0=0;0.05;-0.05;0.05;0;0;0;0;>> y,x=lsim(ac,bc,cc,dc,u,t,x0);>> figure(1);plot(t,y(:,1);>> xlabel('time/sec');ylabel('x(小车)/m'); title('小车位移');grid;>> figure(2);plot(t,y(:,2);>> xlabel('time/sec');ylabel('01/rad');title('下摆角度变化曲线');grid;>> fig