一级倒立摆系统仿真及分析

1一级倒立摆系统仿真及分析1.摘要本次课程设计，我们小组选择一级倒立摆系统作为物理模型，首先通过物理分析建立数学模型，得到系统的传递函数，通过对传递函数的极点，根轨迹，单位阶跃响应来分析系统稳定性。建立状态空间模型，利用matlab进行能控能观性分析，输入阶跃信号，分析系统输出响应。通过设定初始条件，查看系统稳定性，利用simulink绘制系统状态图。再对系统进行极点配置，进行状态反馈，使得系统在初始状态下处于稳定状态，并绘制系统状态图。2.课程设计目的倒立摆系统是一个经典的快速、多变量、非线性、绝对不稳定系统，是用来检验某种控制理论或方法的典型方案。倒立摆控制理论产生的方法和技术在半导体及精密仪器加工、机器人技术、导弹拦截控制系统和航空器对接控制技术等方面具有广阔的开发利用前景。因此研究倒立摆系统具有重要的实践意义。3.课程设计题目描述和要求本次课程设计我们小组选择环节项目三：系统状态响应、输出响应的测量。环节目的：1.利用MATLAB分析线性定常系统。2.利用SIMULINK进行系统状态空间控制模型仿真，求取系统的状态响应及输出响应。环节内容、方法：1.给定系统状态空间方程，对系统进行可控性、可观性分析。并利用SIMULINK绘制系统的状态图，求取给定系统输入信号和初始状态时的状态响应及输出响应。2.给定两个系统的状态空间模型，分别求两个系统的特征值；将两个系统的系统矩阵化为标准型；求出给定系统初始状态时，状态的零输入响应；求两个系统的传递函数并分析仿真结果。24.课程设计报告内容4.1数学模型的建立及分析对于倒立摆系统，由于其本身是自不稳定的系统，实验建模存在一定的困难。但是经过小心的假设忽略掉一些次要的因素后，倒立摆系统就是一个典型的运动的刚体系统，可以在惯性坐标系内应用经典力学理论建立系统的动力学方程。下面我们采用其中的牛顿－欧拉方法建立直线型一级倒立摆系统的数学模型。在忽略了空气阻力，各种摩擦之后，可将直线一级倒立摆系统抽象成小车和匀质杆组成的系统，如下图1所示图l直线一级倒立摆系统我们不妨做以下假设：M小车质量、m摆杆质量、b小车摩擦系数、l摆杆转动轴心到杆质心的长度、I摆杆惯、F加在小车上的力、x小车位置、φ摆杆与垂直向上方向的夹角、θ摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下）。图2，图3分别是系统中小车和摆杆的受力分析图。其中，N和P为小车与摆杆相互作用力的水平和垂直方向的分量。注意：在实际倒立摆系统中检测和执行装置的正负方向已经完全确定，因而矢量方向定义如图所示，图示方向为矢量正方向：3图2小车隔离受力图图3摆杆隔离受力图分析小车水平方向所受的合力，可以得到以下方程：NxbFxM（1）由摆杆水平方向的受力进行分析可以得到下面等式：sin122lkdmN即：sincos2mlmlxmN（2）把这个等式代入上式中，就得到系统的第一个运动方程：FmlmxbxMmsincos2（3）为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，4可以得到下面方程：cos22ldtdmmgP即：cossin2mlmlmgP(4)力矩平衡方程如下：INlPlcossin(5)方程中力矩的方向，由于sinsin,coscos,，故等式前面有负号。合并这两个方程，约去P和N，得到第二个运动方程：cossin2xmlmglmlI(6)设（φ是摆杆与垂直向上方向之间的夹角），假设φ与l（单位是弧度）相比很小，即φ《1，则可以进行近似处理：0,sin,1cos2dtd用u代表被控对象的输入力F，线性化后两个运动方程如下：umlxbxmMxmlmglmlI2(7)4.2传递函数建立及分析对方程组(7)进行拉普拉斯变换，得到：sUssmlssbXssXmMssmlXsmglssmlI22222(8)注意：推导传递函数时假设初始条件为0。由于输出为角度φ，求解方程组(3.8)的第一个方程，可以得到：ssgmlmlIsX22)((9)把上式代入方程组(8)的第二个方程，得到：5sUssmlsssgmlmlIbsssgmlmlImM22222(10)整理后得到传递函数：sqbmglsqmglmMsqmlIbssqmlsUs23242(11)其中：22mlmlIMmq实际的系统模型如下：M小车质量1.096Kgm摆杆质量0.109Kgb小车摩擦系数0.1N/m/secl摆杆转动轴心到杆质心的长度0.25mI摆杆惯量0.0034kg*m*mT采样频率0.005秒将参数代入（11）中，求得q=0.0115635，ssssssUs31.283.270883167.035655.2342利用MATLAB求出该传递函数的极点为:s1=5.2728s2=-5.2781S3=-0.0830s4=0由此可以看出存在正极点，故系统不稳定。对该系统输入单位阶跃信号，其输出信号为：clearnum=[2.35655300]den=[10.0883167-27.83-2.310]g=tf(num,den)[y,t,x]=step(g)6plot(t,y)再利用用MATLAB绘制出根轨迹：clearnum=[2.35655300]den=[10.0883167-27.83-2.310]rlocus(num,den)[k,poles]=rlocfind(num,den);根轨迹如下图所示：通过对传递函数的分析，可知该系统是不稳定的，需要添加控制器来控制系统使其稳定。74.3状态空间结构方程建立及分析系统状态空间方程为DnCXyBuAXx(12)方程组（3.12）对,x解代数方程，得到解如下：uMmlmMImlMmlmMImMmglxMmlmMImlbuMmlmMImlIMmlmMIglmxMmlmMIbmlIxxx2222222222)()()((13)整理后得到系统状态空间方程：222222222201000000000100ImlbxxImlmglIMmMmlxxIMmMmlIMmMmlumlmglMmmlbIMmMmlIMmMmlIMmMml（14）uxxxy0000000001（15）由公式(3.7)的第一个方程为：xmlmglmlI2(16)对于质量均匀分布的摆杆有：231mlI(17)于是可以得到：xmlmglmlml2231(18)8化简得到：xllg4343(19)设xuxxX,,,,，则有：uxxxyulxxlgxx00010000014301004300100000000010(20)把上述参数代入，可以得到系统的实际模型摆杆角度和小车位移的传递函数：26705.00102125.00275.022sssXs(21)摆杆角度和小车加速度之间的传递函数为：26705.00102125.002725.02ssVs(22)摆杆角度和小车所受外界作用力的传递函数：ssssssUs31.283.270883167.035655.2342(23)以外界作用力作为输入的系统状态方程：uxxxyuxxxx000100000135655.20883167.0008285.27225655.0010000629317.00883167.000010(24)以小车加速度作为输入的系统状态方程：9uxxxx301004.2900100000000010uxxxy0001000001（25）首先我们对以小车加速度作为输入的系统状态方程（25）式进行能控能观性分析。04.2900100000000010A3010B01000001C00D利用matlab对系统进行仿真求系统能观能控性代码：A=[0,1,0,0,;0,0,0,0;0,0,0,1;0,0,29.4,0];B=[0;1;0;3];C=[1,0,0,0;0,0,1,0];D=[0;0];Qc=ctrb(A,B);n=rank(Qc);if(n==4),disp('系统能控');else,disp('系统不能控');endQo=obsv(A,C);m=rank(Qo);if(m==4),disp('系统能观');else,disp('系统不能观');End运行结果如下：利用matlab中的simulink绘制系统状态图10现在我们向系统输入单位阶跃信号，利用matlab仿真观察系统单位阶跃响应。代码：A=[0,1,0,0,;0,0,0,0;0,0,0,1;0,0,29.4,0];B=[0;1;0;3];C=[1,0,0,0;0,0,1,0];D=[0;0];step(A,B,C,D)运行结果如下：由此可以看出输出的x与φ均偏向于无穷大，故可以看出系统是不稳定的。;11给系统一个初始状态x0=[3;0;2;0],利用matlab仿真求出系统输出响应代码：A=[0,1,0,0,;0,0,0,0;0,0,0,1;0,0,29.4,0];B=[0;1;0;3];C=[1,0,0,0;0,0,1,0];D=[0;0];x0=[3;0;2;0];sys=ss(A,B,C,D);[y,t,x]=initial(sys,x0,0:0.01:5);plot(t,x);运行结果如下：由上文得出该系统的极点为:s1=5.2728s2=-5.2781S3=-0.0830由于存在正极点，故引起该系统的不稳定，现利用极点配置法来重修配置系统极点，使系统能够处于稳定状态。对于直线一级倒立摆的极点配置转化来说：要按上述系统设计控制器。则要求具有较短，约3s的调整时间和合适的阻尼比ξ=0.5.要使系统具备能控，能观且易验证。步骤为：计算特征值。根据要求，设调整时间为3s，并留有一定的余量，选择期望的闭环极点：s=μi(i=1，2，3，4）,其中：101u,102u,ju3223，ju3224，其中μ3，μ4是一对具有ξ=0.5，ωn=4的主导闭环极点。μ1，μ2位于主导闭环极点的左边，其影响较小，因此期望的特征跟方程为：12s4+24s3+196s2+720s+1600=