倒立摆系统实验设计报告

1、学生实验报告课程名称：倒立摆系统课程设计组号：7姓名：学号：邮箱：2010年11月11 1日目录倒立摆系统的构成 3单级倒立摆数学模型的建立 3传递函数 6状态空间方程 6系统 MATLAB 仿真和开环响应 7稳定性与可控性分析 11控制器设计 12基于状态反馈的控制算法设计与仿真 LQR12极点配置法 16PID控制算法19实验结果及与仿真结果的对比分析 29感想和建议 30倒立摆系统的构成图1倒立摆系统的组成框图如图1所示为倒立摆的结构图。系统包括计算机、 运动控制卡、伺服机构、 倒立摆本体和光电码盘几大部分， 组成了一个闭环系统。 光电码盘1将小车的位移、速度信号反馈给伺 服驱动器和运动

2、控制卡，摆杆的位置、速度信号由光电码盘 2反馈回控制卡。计算机从运动 控制卡中读取实时数据，确定控制决策（小车向哪个方向移动、移动速度、加速度等），并由运动控制卡来实现该控制决策，产生相应的控制量，使电机转动，带动小车运动，保持摆 杆平衡。单级倒立摆数学模型的建立在忽略了空气流动，各种摩擦之后，可将倒立摆系统抽象成小车和匀质杆组成的系统, 如下图2所示图2单级倒立摆模型示意图 那我们在本实验中定义如下变量：（本实验系统 （本实验系统 （本实验系统小车质量摆杆质量小车摩擦系数0.5 Kg)0.2 Kg)0.1 N/m/sec)0.3 m)x00摆杆转动轴心到杆质心的长度（摆杆惯量 （0.006

3、kg*m*m加在小车上的力小车位置摆杆与垂直向上方向的夹角摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下） 图3是系统中小车和摆杆的受力分析图。其中,下面我们对这个系统作一下受力分析。N和P为小车与摆杆相互作用力的水平和垂直方向的分量。注意：在实际倒立摆系统中检测和执行装置的正负方向已经完全确定, 义如图，图示方向为矢量正方向。因而矢量方向定-blRicticr图3倒立摆模型受力分析分析小车水平方向所受的合力，可以得到等式：应用Newt on方法来建立系统的动力学方程过程如下:分析小车水平方向所受的合力，可以得到以下方程：Mx丄 F - bx - N由摆杆水平方向的受力进行分析可以得到

4、下面等式:N =mdt2(x - l sin T)即N =mx 亠 ml vcosv _ ml/sin 71把这个等式代入上式中，就得到系统的第一个运动方程:n aaQ(M m)x bx cos - mlsin - F(1)为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得到下面方程:(l COST)即：P - mg - -mb sin - mH2 cost力矩平衡方程如下:-Pl sin - Nl cos)- I注意：此方程中力矩的方向，由于j -二 ，cos - -cosj,s in- - si nr，故等式前面有负号。合并这两个方程，约去 P和N，由I二吉ml2得到

5、第二个运动方程:4ml mgl sin)- -mlxcos(2)设- r： ( 是摆杆与垂直向上方向之间的夹角)，假设与1 (单位是弧度)相比d0很小，即1,则可以进行近似处理：cos = ，sin - - ，()2 =0。用U来代dt表被控对象的输入力 F，线性化后两个运动方程如下：心x( 3)i(M +m)x+bx ml = u传递函数对方程组（3）进行拉普拉斯变换，得到14(M +m)X(s)s2+bX(s)s-ml (s)s2 =U(s)注意：推导传递函数时假设初始条件为0。由于输出为角度，求解方程组（4）的第一个方程，可以得到4 g _ X（s） l2G（s）3 s把上式代入方程组（

6、4）的第二个方程，得到,(l +ml2) g |乱 2,(l +ml2) g |民,2(M m)(s)s2 b2 .：(s)s_ml(s)s =U(s)| ml sml s整理后得到传递函数：ml 2 sU(s)二q41 bml2 3 (M m)mgl 2 bmglssssqqq其中q=(M m)(l ml2)(ml)2状态空间方程系统状态空间方程为X 二 AX Buy 二 CX Du方程组（3）对x, 解代数方程，得到解如下:-4b(4M m)3mg(4M m)4(4M +m)-3bx 3g(M m)(4M m)lu(4M m)l (4M m)l整理后得到系统状态空间方程:00L*1-4b(

7、4M m)0-3b(4M m)l03mg(4M m)03g(M m)(4M m)l0 x010u0系统Matlab仿真和开环响应实际系统参数如下，求系统的传递函数、Matlab 仿真。M小车质量m摆杆质量b 小车摩擦系数I摆杆转动轴心到杆质心的长度I摆杆惯量F加在小车上的力x小车位置0摆杆与垂直方向的夹角T采样频率状态空间万程，并进行脉冲响应和阶跃响应的1.096 Kg0.109 Kg0 .1N/m/sec0.2 5m0.0034 kg*m*m0.005 秒传递函数：在 Matlab 中，拉普拉斯变换后得到的传递函数可以通过计算并输入分子和分母矩阵来实现。求系统传递函数的m-文件内容如下：M

8、= 1.096;m = 0.109;b = 0.1;I= 0.0034;g = 9.8;l = 0.25;q = (M+m)*(l+m*L2)-(m*l)A2;%simplifies in putnum = m*l/q 0den = 1 b*(I+m*lA2)/q -(M+m)*m*g*l/q -b*m*g*l/qt = 0 : 0.01 : 5;impulse ( num , den , t )axis ( 0 1.1 0 70 )执行上面的文件， 就可以求出系统传递函数的分子与分母多项式的 Matlab 表示：num =2.3566 0den =1.0000 0.0883 -27.8285

9、 -2.3094可以 得 到 系 统 开 环 脉冲 响 应的 曲 线 如 下 ：70Impulse Resp onse6050403020OedLIP m6.0图4系统开环脉冲响应曲线状态空间法：状态空间法可以进行单输入多输出系统设计，（从实验二开始，我们将尝试同时对摆杆角度和小车位置进行控制）。为了更具挑战性，给小车一个阶跃输入信号，设计指标如下：小车位置x和摆杆角度0的稳定时间小于5秒；位置x的上升时间小于 0.5秒；摆杆角度的超调量小于 20度（0.35弧度）。下面，我们用 Matlab求出系统的状态空间方程各矩阵，并仿真系统的开环阶跃响应。在这里同样给出了一个 m-文件，执行这个文件，

10、Matlab将会给出系统状态空间方程的A,B, C和D矩阵，并可以绘出在给定输入为一个0.2 m的阶跃信号时系统的响应曲线。M = 1.096;m = 0.109; b = 0.1;I= 0.0034;g = 9.8;l = 0.25;p = l*(M+m)+M*m*L2;%de nomi nator for the A and B matriciesA = 0100;0 -(I+m*lA2)*b/p (mA2*g*lA2)/p0;0 0 0 1;0 -(m*l*b)/pm*g*l*(M+m)/p 0B = 0;(I+m*lA2)/p;0;m*l/pC = 1 0 0 0;0 0 1 0D =

11、 0;0T=0:0.005:10;U=0.2*ones(size(T);Y,X=lsim(A,B,C,D,U,T);plot(T,Y)axis(0 2.5 0 100)执行该 m 文件，可以求出系统的状态空间A、B、C、D 矩阵 ,得到开环系统阶跃响应的曲线。A =0 1.0000 000 -0.0883 0.629300 0 01.00000 -0.2357 27.82850B =00.883202.3566C =100908070605040302010 -0 L00.511.522.5图5倒立摆状态空间开环系统阶跃响应曲线图中，实线是摆杆角度响应曲线，虚线是小车位置响应曲线。稳定性与可控

12、性分析我们先看一看系统的稳定性，将数据代入状态方程中，利用 的零极点。源代码如下：M = 1.096;m = 0.109;b = 0.1;I= 0.0034;l = 0.25;a b c d=wer_ss(M,m,b,l);%自己编写的函数，建立模型之用,matlab程序可以求出系统具体程序见下面sysc=ss(a,b,c,d);sysd=c2d(sysc,0.005);da db dc dd=ssdata(sysd);z p gain=ss2zp(da,db,dc,dd,1)-0.99991.02750.9733-0.99991.0000 + 0.0000i1.0000 - 0.0000i1

13、.00000.99961.02850.9723gain =1.0e-004 *0.11130.3338wer_ss源程序：fun cti ona b c d=wer_ss(M,m,b,l)a=01 0 0;0 -4*b/(4*M+m)3*m*9.8/(4*M+m)0;00 0 1;0-3*b/(4*M+m)*l)3*9.8*(M+m)/(4*M+m)*l) 0;b=0;4/(4*M+m);0;3/(4*M+m)*l);c=1 0 0 0;0 0 1 0;d=0;0由得到的p (极点)可知，有的极点在单位圆外，所以可知原系统是不稳定的。同样，我们可以利用matlab来得到系统的能控性，源代码如下

14、：ud=ctrb(da,db);ran k(ud)ans =4由得到的rank( ud)的值可知，原系统的能控性矩阵为4，所以我们可知原系统是能控的。控制器设计基于状态反馈的控制算法设计与仿真 LQR由理论分析知，可以设计基于最优控制的状态调节器，使系统闭环稳定。设状态反馈调节律的形式为u(k)二-Krx(k)通过使性能指标函数1 1 N JLJ Xt(N)PX(N) XT(k)Qr(k)X(k) uT(k)Rr(k)u(k)2 2y为最小，根据在附录1中我们所介绍的求得T1 TKr 工Rr H P(K 1)H H P(k 1)G其中p由下列黎卡提方程获得P(k) =Qr GTP(k 1)G

15、_GTP(k 1)HRr H Tp(k 1)H H Tp(k 1)G其中Qr，Rr分别用来对状态向量 x(k)，控制向量U(k)引起的性能度量的相对重要性进 行加权。在实际运算中我们运用 Matlab控制系统工具箱中的dlqr”函数直接进行运算。利用dlqr函数，我们需要提供两个权值矩阵：Q、R。通常我们取R=1，而对于Q我们只能通过不断的凑取来得到。源代码如下： R=1; Q=10 0 0 0;0 0 0 0;0 0 1 0;0 0 0 0Q =10 0000 0000 0100 000 T=0.005; syse k=wer_lqr(da,db,dc,dd,Q,R,T);%wer_lqr

16、是自己定义的函数，具体见下面程序 x0=0.05;0;0.0175;0; t=0:0.005:10; y x1=in itial(syse,x0,t);plot(t,y(:,1),red,t,y(:,2),blue)wer_lqr源程序：fun cti onsysresult k=wer_lqr(da,db,dc,dd,Q,R,T);%sysresult k=wer_lqr(da,db,dc,dd,Q,R,T);k S e=dlqr(da,db,Q,R);G=da-db*k;sysresult=ss(G,db,dc,dd,T);我们已开始的 Q为：Q1=10 0 0 0;0 0 0 0;0 0

17、 1 0;0 0 0 0；结果得到的图为图6；Q2取为：Q=100 0 0 0;0 0 0 0;0 0 10 0;0 0 0 0;结果得至曲勺图为图 7 ；通过比较，我们发现当Q11、Q33比值一定时，取大的值时系统的响应速度加快，但是超调加大；反之则响应变慢但是超调减小。0.03 -小车位移-0.01 Lc1J012345678910图6 Q1响应图0 12345678910图7 Q2响应图 在左右权衡之间，我们最终选取了：Q=300 0 0 0;0 0 0 0;0 0 30 0;0 0 0 0此时的响应曲线如图8，而k值为：-16.6147-12.422656.590910.24440.0

18、50.040.030.020.01-0.01 匚0小车位移8910图8最优的响应曲线此时的单位阶跃响应曲线为图9:0.020.010-0.01-0.02-0.03-0.04-0.05-0.06-0.0710摆角小车位移图9单位阶跃响应曲线从仿真效果来看，零状态响应和单位阶跃响应都符合要求。极点配置法采用极点配置法设计多输出的倒立摆系统的控制方案。 制摆杆和小车的位置。图10是控制系统的示意图。可以用完全状态反馈来解决,图10控制系统框图假定所有的状态变量都可以测量和反馈， 可以证明：若所研究的系统是状态完全可控的， 那么，利用状态反馈的方法， 经过适当的状态反馈增益矩阵， 就可以把闭环系统的极

19、点配置 到任何期望的位置。设开环控制系统的离散状态方程为：x(k+1) = Gx(k) + Hu(k)其中，假设系统是状态完全可控的x(k)为在第k次采样时刻的状态矢量(n维矢量)u(k)为在第k次采样时刻的控制信号(标量)G = n n矩阵H = n 1矩阵设极点配置的控制律形式为u(k) 一KpX(k)K式中Kp是状态反馈增益矩阵(1 n矩阵)，于是该系统就成为一个闭环控制系统。其闭环状态方程为x(k 1) =(G HKp)x(k)G HK注意，p的特征值就是所要求的闭环极点我们希望利用状态反馈u(k) = -K pX(k)把闭环极点布置在，z八2 ,z = %。即要求特征方程为：zl -

20、G+HKp =(z-气)(z-巴)(z 叫)_ n -1 .n-2 -a1z a2zanz a 0根据Cayley_hamiton定理，经过推导(此略)可以得到Kp = 0 00 1H GHGn4H (G)其中(G) =Gn a1Gn a2GnanG a.IKk上式给出所要求的状态反馈增益矩阵K P。矩阵K P的这种特殊表达式就是常说的阿克曼公式。状态反馈增益矩阵按这样的方法确定，即可使误差（由扰动所引起的）以足够快的速度K降到零。注意，对于一个给定的系统，矩阵p并不是唯一的，而是取决于所期望的闭环极点位置（它决定响应速度）的选择。选择期望的闭环极点或期望的特征方程是在误差矢量响应的快速性与对

21、扰动和测量噪声敏感型之间的一个折衷方案。 应的速度提高，那么扰动和测量噪声的有害影响往往也会增强。也就是说，如果我们使误差响 在确定给定系统的状态反馈Kk增益矩阵p时，通常是通过比较按不同的期望闭环极点或期望特征方程得到的矩阵p，K并从中选出使整个系统达到的特性最好的那个矩阵p。在实际设计时，我们运用Matlab控制系统工具箱中的“ place ”函数直接进行仿真和运算。先在连续域中进行计算，然后再转到离散域。根据系统的性能要求，我们可取 =0.707，Wn =3，这时完全满足题中的性能要求。那么我们就可以进行离散域设计了： z仁exp(-3*2A0.5)/2+(-3*2A0.5)*j/2)*

22、0.005) z1 =0.9894 - 0.0105i z2=exp(-3*2A0.5)/2-(-3*2A0.5)*j/2)*0.005) z2 =0.9894 + 0.0105i z3=exp(-10*0.005)z3 =0.9512 z4=exp(-12*0.005)z4 =0.9418 p=z1 z2 z3 z4;K=place(da,db,p)K =-38.6579-25.5096103.324717.9041 G=da-db*K; syse=ss(G,db,dc,dd,0.005); t=0:0.005:10; x0=0.05; 0; 0.0175; 0; y1,x=in itial

23、(syse,x0,t); plot(t,y1(:,1),red,t,y1(:,2),blue)0.060.050.040.03小车位移0.020.01-0.01 L0910图11极点配置图零输入响应 而它的单位阶跃响应是：u=on es(1,le ngth(t); y,x=lsim(syse,u,t)0.0150.010.005摆角0-0.005-0.01-0.015-0.02小车位移-0.025-0.03图12极点配置单位阶跃响应 从仿真效果来看，也是基本上达到了系统的要求。PID控制算法目的：设计PID控制器，使得当在小车上施加1N的脉冲信号时，闭环系统的响应指标为:1、稳定时间小于5秒2

24、、 稳态时摆杆与垂直方向的夹角变化小于0.1弧度分析：系统输出量为摆杆的位置，它的初始位置为垂直向上，我们给系统施加一个扰动，观察摆杆的响应。系统框图如下：图13 PID系统框图图中KD (s)是控制器传递函数，G(s)是被控对象传递函数。考虑到输入r(s) =0，结构图可以很容易的变换成CoritrollRf图14 PID系统反馈控制框图该系统的输出为y(s)G(s)1 KD(s)G(s)F(s)num门而臥而F(s)(denPID)(den)num(de nPID)(den PID)(den) (nu mPID)( num)其中，num 被控对象传递函数的分子项被控对象传递函数的分母项nu

25、mPID PID控制器传递函数的分子项denPID PID控制器传递函数的分母项被控对象的传递函数是mls(s)qnumU (s)3 售 bml2 2 (M m)mgl bmgl dens +ssqqq其中2 2q 二(Mm)(I ml ) -(ml)PID控制器的传递函数为KD(s)二 KDs KP Ks2KDs2 KPs KIsnu mPID den PID调节PID控制器的各个参数，以得到满意的控制效果。前面讨论的输出量只考虑了摆杆角度，那么，在我们施加扰动的过程中，小车位置如何变化？考虑小车位置，得到改进的系统框图如下:ControllerControllerr(j) = 04-图15

26、改进的PID系统控制框图由于输入信号r(s)=O，所以可以把结构图转换成:PgConWllerConWllerKD&)KD&)图16等价PID系统控制框图其中，反馈环代表我们前面设计的摆杆的控制器。注：从此框图我们可以看出此处只对摆杆角度进行了控制，并没有对小车位置进行控制。小车位置输出为:num2X(s)二G2(s)1 KD (s)G1 (s)F(s)二den21 + (numPID Xnurni)(de nPID)(denJF(s)(num2)(denPID)(dem)(de nPID)(de nj(de n2) (n umPID)( num/de n2)F(s)其中，num1，den1，

27、 num2， den2分别代表被控对象 1和被控对象2传递函数的分子和分母。numPID和denPID代表PID控制器传递函数的分子和分母。下面我们来求g2(s)，根据前面实验二的推导，有X(s)二白-耳：(s)3 s可以推出小车位置的传递函数为G2 (s)X(s)U(s)mglq423 +Ebml 2 (M +m)mgls bmgl sssqqq其中q 二(Mm)(I ml2) -(ml)2可以看出，dem = de n2=de n，小车的算式可以简化成:F(s)X(s)血碣(知叫)(denPID)(den ) + k(numPID)(numi)Matlab仿真实际系统参数如下:M小车质里1

28、.096 Kgm摆杆质量0.109 Kgb小车摩擦系数0 .1N/m/secl摆杆转动轴心到杆质心的长度0.2 5mI摆杆惯量0.0034 kg*m*mF加在小车上的力x小车位置0摆杆与垂直方向的夹角T采样频率0.005 秒注意：在进行实际系统的Matlab仿真时，请将采样频率改为实际系统的采样频率1. 参照上例的处理方法，把实际系统参数代入，进行PID控制仿真，并找到合适的参数；2. 将小车推到导轨正中间位置，并且使摆杆处于自由下垂的静止状态；3. 给计算机和电控箱通电；4. 打开计算机，在 DOS操作系统下，键入 Pend”启动程序，并按 S”使系统处于准备状 态；5. 按“f”键起摆，等

29、摆杆立起来并稳定下来之后（此时为LQR控制），选择”控制器”菜单中的“ PID”，输入参数，观察小车和摆杆的运动；（注意由于控制器只对摆杆进行了控制，所以在 PID中输入参数后小车可能向一个方向运动，此时需用手轻轻扶一下摆杆）6. 按“T”停止摆杆，当摆杆处于自由下垂的静止状态，并在”控制器”菜单中选择“PID”，并参照仿真结果，输入PID控制器参数，观察小车和摆杆的运动；（通过调整参数可以控制摆杆摆起并能够竖直向上，此时可能需用手轻轻扶一下摆杆，以避免小车“撞墙”)。7. 如果控制效果不理想，调整控制器参数，直到获得较好的控制效果；&运动曲线文件保存在当前路径下，其中pos_x.dat为小车

30、位置，pos_t.dat为摆杆角度，con_u.dat为控制量；9. 到MatLab中运行如下指令，MatLab会将保存的曲线重新绘制出来。S=load(路径文件名)plot(S)1摆杆角度控制首先，需要把输出为摆杆角度勺寸系统的传递函数用Matlab表示出来，建立一个m-文件，将下面几行表示传递函数的语句拷贝进去，其中k代表比例系数 Kp :M = 1.096;m = 0.109;b = 0.1;I = 0.0034;g = 9.8;l = 0.25;q =(M+m)*(l+m*L2) -(m*l)A2;%simplifies in putnum = m*l/q 0 0den = 1 b*(

31、l+m*L2)/q -(M+m)*m*g*l/q -b*m*g*l/q 0kd=1k=1ki=1n umPID= kd k ki ;den PID= 1 0 ;n umc= conv ( n um, den PID )denc= polyadd ( con v(de nPID, den ), conv( n umPID, num )t = 0 : 0.05 : 5;impulse ( n umc , dene , t )其中函数polyadd是求两个多项式之和的函数，它不是Matlab工具，因此必须把它拷贝到polyadd.m文件中，并把该文件的目录用 addpath命令加到路径中。Polya

32、dd函数内容如 下：fun cti on poly=polyadd(poly1,poly2)if len gth(poly1)0poly=zeros(1,mz),short+long;elsepoly=long+short;end在这里我们假定比例、积分和微分控制都是必需的。m-文件中加入下面的语现在，就可以进行系统脉冲响应的 PID 控制仿真了。在前面的 句就可以得到系统的脉冲响应仿真结果：运行结果及响应曲线如下：num =2.3566 0 0den =1.0000 0.0883 -27.8285 -2.3094 0kd =1k =1ki =1numc =2.3566 0 0 0denc =

33、1.0000 2.4449 -25.4720 0.0471 0 0x 10Impulse Resp onsew ith P ID control:k=1,kd=1,ki=1987654321000.511.522.53Time (sec)3.544.55Kp，Kd和Ki，直到获得满意图17初始PID参数摆角状态图系统响应是不稳定的，不能满足要求，需要调整参数的控制结果。首先增加比例系数Kp，观察它对响应的影响，取KP=100，kd=1.系统响应如下：0 00.0883-27.8285-2.30940num =2.3566 den =1.0000 kd =1k =100ki =1numc =2.

34、3566dene =1.00002.4449207.82680.04715ToImpulse Responsew ith PID control k=1OO,ki=1,kd=1LLIiI II i0505o-0-0.15-0.200.511.522.533.544.55Time (sec)图18调节PID参数后摆角输出图系统稳定时间约为4秒，满足要求。由于此时稳态误差为 0,所以不需要改变积分环节(你可以改变积分系数，观察系统响应如何变化)；系统响应的超调量比较大，为了减小超调，增加微分系数Kd，取Kd=2o，响应结果和响应曲线如下：num =2.3566 00den =1.00000.088

35、3 -27.8285 -2.30940kd =20k =100ki =1numc =2.3566dene =1.0000 47.2194 207.82680.0471000.10.090.080.070.060.050.040.030.020.010Impulse Responsew ith PID control k=100,ki=20,kd=100.511.522.533.544.55Time (sec)图19微调PID参数后摆角状态曲线系统稳定时间约为1秒，超调约为0.04，响应满足指标要求。2 小车位置变化仿真仿真小车位置变化的m-文件内容如下：M = 1.096;m = 0.109;

36、b = 0.1;I = 0.0034;g = 9.8;l = 0.25;q =(M+m)*(I+m*lA2) -(m*l)A2;%simplifies in putnum1 = m*l/q 0 0den1 = 1 b*(l+m*L2)/q -(M+m)*m*g*l/q -b*m*g*l/q 0num2 = -(I+m*lA2)/q 0 m*g*l/qden2 = den1kd = 20k = 100ki = 1numPID = kd k ki;den PID = 1 0;numc = conv(num2,denPID)denc = polyadd(conv(denPID,den2),conv(

37、numPID,num1) t=0:0.05:5;impulse(numc,denc,t) 仿真结果如下：num1 =2.3566 0 0den1 =1.0000 0.0883 -27.8285 -2.3094 0num2 =-0.8832 0 23.0942den2 =1.0000 0.0883 -27.8285 -2.3094 0kd =20k =100ki =1numc =-0.8832 0 23.0942 0denc =1.0000 47.2194 207.8268 0.0471 0 00.60.5 _” _0.4 _e 0.3pm 0.2J0.1 L.J *0干-0.1 II :00.511.522.533.544.55Time (sec)图20小车位置曲线实验结果及与仿真结果的对比分析我们组于2010年10月前往倒立摆实验室做实验，我们共做了PID、LQR和极点配置法，由于事先已经将参数设计好并经过精心的仿真修改，所以实