2011计算机仿真大赛一等奖作品解析

编号：重庆邮电大学第五届计算机仿真大赛一等奖作品题号: C组别：高年级组2011年 4 月 6 日单级移动式倒立摆动态控制仿真摘要单级移动式倒立摆系统是一个典型的单输入双输出的自然不稳定系统，具有非线性、强耦合、多变量、高阶次等特性，作为控制系统的被控对象，它是一个理想的教学实验设备，许多抽象的控制概念都可以通过倒立摆直观地表现出来。

应用上，倒立摆广泛应用于控制理论研究、航空航天控制、机器人、杂技项杆表演等领域，在自动化领域中具有重要的价值。

研究倒立摆的精确控制对工业生产中复杂对象的控制有着重要的应用价值，因此倒立摆仿真或实物控制实验是控制领域中用来检验某种控制理论或方法的典型方案。

许多抽象的控制概念如控制系统的稳定性、可控性、可观性、系统收敛速度和系统抗干扰能力等，都可以通过倒立摆系统直观的表现出来。

在控制过程中，它能有效地反映诸如可镇定性、鲁棒性、随动性以及跟踪等许多控制中的关键问题，是检验各种控制理论的理想模型。

因而其研究具有重大的理论和实践意义。

本文以单级移动式倒立摆为被控对象，参考古典控制理论设计控制器（如PID控制器）的设计方法和现代控制理论设计控制器（最优控制）的设计方法。

针对单级移动式倒立摆系统的单输入双输出、非线性、强耦合的不稳定性，我们采用线性二次最优LQR控制和PID控制，通过MATLAB软件进行实时控制（摆杆的角度和小车的位移）和系统仿真实验实现对小车位置和摆杆偏角的同时闭环控制，达到了预期效果（动态仿真）。

关键词：非线性倒立摆实时控制系统仿真动态1.问题重述单级移动式倒立摆原理如下图1-1 所示，图1-1其中小车质量M = 3KG;摆杆质量m = 0.1KG;摆杆为均匀的，质心l=0.5米;重力加速度g=9.81;设x1为摆杆角度，x2为摆杆角速度，x3为小车位移，x4为小车速度，其中小车水平受力为输入量，建立其状态方程：采用状态反馈进行控制，其控制原理如下图1-2 所示：图1-2计算得k1=-446.7287,k2= -91.0448,k3= -244.6376 ,k4=-110.0896 时系统稳定。

建立倒立摆控制仿真系统，当输入信号为阶跃信号和随机信号时，实时演示倒立摆的工作情况。

2.单级倒立摆模型的建立2.1 倒立摆概述单级倒立摆是一个自然不稳定系统，是处于倒置不稳定状态，人为控制使其处于动态平衡的一种摆，一般是由一个可以在水平轨道上自由移动的小车和倒置摆铰链而成，用一种强有力的控制方法对小车的速度作适当的控制，从而使全部摆杆倒置稳定于正上方，这个系统也叫自动平衡车。

悬挂式倒立摆系统开始工作时，首先使小车按摆杆的自由振荡频率摆动，摆杆随之大幅度摆动，当摆杆接近于倒立摆垂直位置时，自动转换控制方法，使其稳定于倒置状态。

随着倒立摆系统的控制研究的不断深入，倒立摆系统的种类也由简单的单线倒立摆发展为多种形式的倒立摆。

常见的倒立摆系统一般由小车和摆杆两部分构成(另有旋转式倒立摆等形式)，其中摆杆可能是一级、两级甚至多级(级数越多，控制难度越大)，其长度也可能是可变化的。

控制的目标一般都是通过给小车施加一个水平方向的力，使小车在期望位置上稳定，而摆杆达到竖直向上的动态平衡状态。

2.2 数学准备对系统建立数学模型是系统分析、设计的前提，而一个准确又简练的数学模型将大大简化后期的工作。

为了简化系统分析，在实际的模型建立过程中，要忽略空气流动阻力，以及各种次要的摩擦阻力。

这样，可将倒立摆系统抽象成小车和匀质刚性杆组成的系统，如下图所示。

本系统内部各相关参数定义如下：M 小车质量m 摆杆质量b 小车摩擦系数l 摆杆转动轴心到杆质心的长度22I 摆杆惯量f加在小车上的力x小车位置φ摆杆与垂直向上方向的夹角θ 摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下） 下图是系统中小车和摆杆的受力分析图。

其中， N 和 P 为小车与摆杆相互作用力 的水平和垂直方向的分量。

注意：在实际倒立摆系统中检测和执行装置的正负方向已经完全确定，因而矢量 方向定义如图，图示方向为矢量正方向。

应用 Newton 方法来建立系统的动力学方程过程如下： 分析小车水平方向所受的合力，可以得到以下方程：M x= f - bx - N由摆杆水平方向的受力进行分析可以得到下面等式：N = m d( x + l sin θ )dt 2即： N = m x+ ml θ cos θ - ml θ 2 sin θ把这个等式代入上式中，就得到系统的第一个运动方程：(M + m ) x+ bx + ml θ cos θ - ml θ 2 sin θ = F（2-1）为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以 得到下面方程：P - mg = -m d(l cos θ )dt 2即： P - mg = ml θ sin θ + ml θ 2 cos θ力矩平衡方程如下：- Pl sin θ - Nl cos θ = I θ注意：此方程中力矩的方向，由于θ = π + φ , cos φ = - cos θ , s in φ = - s in θ ，故等式前面有负号。

合并这两个方程，约去 P 和 N ，得到第二个运动方程：(I + ml 2 )θ + mgl sin θ = -ml xcos θ （2-2）⎨ 222.3 模型建立2.3.1 微分方程建立设θ = π + φ ，当摆杆与垂直向上方向之间的夹角φ 与 1（单位是弧度）相比很小，即 φ << 1 时，则可以进行近似处理： cos θ = -1 ， sin θ = -φ ， ( d θ ) 2 = 0 。

为了与控dt制理论的表达习惯相统一，即u 一般表示控制量，用u 来代表被控对象的输入力 F ，线 性化后得到该系统数学模型的微分方程表达式：⎪⎧(I + ml 2 )φ- mgl φ = ml x ⎪⎩(M + m )x + bx - ml φ = u (2-3)2.3.2 传递函数建立对方程组（2-3）进行拉普拉斯变换，得到⎪⎧(I + ml 2 )Φ(s )s 2 - mgl Φ(s ) = mlX (s )s 2⎨ ⎪⎩(M + m ) X (s )s 2 + bX (s )s - ml Φ(s )s 2 = U (s )注意：推导传递函数时假设初始条件为 0。

由于输出为角度φ ，求解方程组（2-4）的第一个方程，可以得到（2-4）X (s ) = [ (I + ml ) - ml gΦ(s )s 2把上式代入方程组（2-4）的第二个方程，得到⎡ +2 ⎤ ⎡ 2 ⎤ (I (M + m )⎢ ml ) - g (I + ml ⎥Φ(s )s 2 + b ⎢ ) g - ⎥Φ(s )s - ml Φ(s )s 2 = U (s )⎣ ml s 2 ⎦ ⎣ mls 2 ⎦整理后得到以输入力u 为输入量，以摆杆摆角φ 为输出量的传递函数：ml s2G (s ) = Φ(s ) = qU (s ) s 4 + b ( I + m l ) s 3 - (M + m )mgl s 2 - b mgl sq q q其中q = [(M + m )(I + ml 2 ) - (ml ) 2 ]若取小车位移为输出量，可得传递函数：⎡ 0 1 0 0⎤⎢ M + m ⎢ g 0 0 ⎥ 0⎥0 0 1 0 0 1 00 0T 0]2 2 T( I + m l 2 ) mgls 2 -G (s ) = X (s ) = q q U (s ) s 4 + b ( I + m l ) s 3 - (M + m )mgl s 2 - b mgl sq q q2.4 状态空间模型2.4.1 状态空间模型的建立 状态空间模型的一般形式如下：x= Ax + Bfy = Cx + Df2.4.2 系统向量的设定设 x1 为摆杆角度，x2 为摆杆角速度，x3 为小车位移，x4 为小车速度Tx = ⎡⎣ x 1 x 2 x 3 x 4 ⎤⎦⎡10⎤⎢0 0⎥ A = ⎢Ml ⎥C = ⎢ ⎥ ⎢ 0 0 0 1⎥⎢0 0⎥⎢⎥ ⎢ - m g 0 0 0⎥ ⎢01⎥⎣⎢ M⎥⎦⎣ ⎦B = ⎡0 - 1 01 ⎤ D = [0 0 ⎣⎢ Ml M ⎥⎦2.4.3 系统的状态空间模型将上面的向量带入状态方程的一般形式，得到如下的状态方程：3.单级倒立摆模型性能分析3.1 系统能控性分析系统能控性就是研究系统的全部状态是否都会受到输入的影响，从而实现对系统状态的控制。

能控性体现了系统输入对状态的控制能力。

[1]根据系统能控性定义可知，系统完全能控的充分必要条件是：能控性矩阵W=⎡⎣A AB A n-1 B⎤⎦为行满秩。

c仿真程序：[2]n=length(A);Q=zeros(n);Q(:,1)=B;for i=2:n %按照公式Q=[B AB A^2*B....A^(n-1)*B]Q(:,i)=A*Q(:,i-1);endm=rank(Q);if(m==n)disp('该系统可控！');melsedisp('该系统不可控！');mend仿真结果：T仿真结果分析：图 3-1-1由仿真结果知，可控性矩阵的秩为 4，即为行满秩，故该系统具有可控性。

3.2 系统能观性分析 系统的能观性就是系统状态变量的任何运动完全可以由输出来反映。

能观性表征了输出对状态的反应能力。

[1] 根据系统能观性定义知，能观性完全能观的充分必要条件是：能观性矩阵W 0 = ⎡⎣C CACA n -1 ⎤⎦为列满秩。

仿真程序： V=obsv(A,C)%求能观性矩阵 rank(V)%求能观性矩阵的秩 仿真结果:图 3-2-1 由仿真结果知，可观性矩阵的秩为 4，即为为列满秩，故该系统具有可观性。

3.3 阶跃响应分析3.3.1 开环系统的阶跃响应仿真程序：%求解系统的特征值%p=eig(A)t1=0:0.01:1;step(A,B,C,D,1,t1) %阶跃响应仿真结果：图3-3-1 原系统阶跃响应曲线结果分析：由于特征方程有一个正实部的根，故系统不稳定。

3.3.2 闭环系统的阶跃响应加入反馈K=[-446.7287 -91.0448 -244.6376 -110.0896]后的阶跃响应仿真程序：% 带反馈 T1=0:0.01:5;u1=0.25*ones(size(T1));[Y,X]=lsim(A-B*K,B,C,D,u1,T1);plot(T1,Y(:,2));title('加反馈后系统的阶跃响应曲线');grid on仿真结果：图3-2-2仿真结果分析由图3-2-2 可知，系统的快速性很好，过渡过程时间不超过2.5 秒。