直线一级倒立摆系统的建模及仿真

计算机控制技术课程设计实验：直线一级倒立摆系统的建模及仿真、已知条件:摆杆角度为输出，小车的位移为输入。导轨中点为坐标轴的中心即零点，右向为坐标值增加的方向，杆偏移其瞬时平衡位置右侧的角度为正值。二、任务要求：总体任务通过调节PID参数，设计PID控制器实现摆杆在受到干扰的情况下，依然能恢复平衡。具体包括以下几部分：理论推导包括倒立摆系统的动力学模型，传递函数，离散传递函数，状态空间或差分方程，稳定性分析，PID控制器设计程序实现实现内容：倒立摆系统模型，控制器以及仿真结果的显示。开发语言和工具：Matlabm文件或C/C++（工具：VC++或其它）PID控制参数设定及仿真结果。分别列出不同杆长的仿真结果（例如：L=0.25和L=0.5）。将理论推导、程序实现、仿真结果写成实验报告。具体求解过程如下：一，倒立摆系统动力学模型的建立1）图1摆杆的受力分析图以摆杆为研究对象，对其进行受力分析，如图11）FxCP+mgxCP-M=0IC I式中，F为杆的惯性力，表达式为F=ma=mCa+at+an)=F+Ft+FnIC IC C PCPCPIPicpicp，m为杆的质量，g为重力加速度，M为杆的惯性力偶。惯性力及惯性力偶的大小分别为d2x d20 1 d202）F=ma=m,Ft=ma=m,M=Ja=_mL22）IPPdt2ICP dt2Ic3 dt2式中，a为杆的角加速度，a为小车的加速度，2L为杆的长度，0为杆偏离中P心位置的角度，X偏离轨道中心的位移。对（2）式代入（1）式，并整理可得TOC\o"1-5"\h\zL^0-gsin0=-cos0 （3）3dt2 dt2当摆动较小时，可以进行近似处理sin0q0,cos0沁1。故（3）式可化为4d20 d2xL -g0=- （4）3dt2 dt2对（4）式进行拉普拉斯变换得3Ls20（s）-g0（s）=-s2X（s） （5）

则系统的开环传递函数为()0(s) -s2 -as2GW_ _ _—X(s)4Ls2-g s2-ag3式中，a_-4L二，离散传递函数，差分方程、状态空间及稳定性分析1，离散传递函数将连续系统离散化，根据连续传递函数G(s)可求得相应的脉冲传递函数为1-bz1-bz-i+z-27)1-e-TsG(s)_(1-z-1)Z-asss2-agG(z)=Z式中，b=eagT+e—agT。将参数L二1,g二10,T二0.01代入(7)式得该参数下的脉冲传递函数为G(zG(z)_—0.75+1.50z-i-0.7503z-1-2.001z-i+z-28)2，根据离散传递函数求系统差分方程由阳)®(z)-a+(a+ab/2)z-1-abz-2/2Glz丿_ _1-bz-1+z-1-bz-1+z-2可得z-1x(z)-吳-2x(z)弘)-bz啟)+z-1x(z)-吳-2x(z)进行反变换即可得到对应的差分方程6(k)6(k)-b6(k-1)+e(k-2)=一ax(k)'ab'

+a+——xI2丿(k-])一#x(k-2)10)将参数L二1,g二10,T二0.01代入(10)式得该参数下的差分方程为9(k)—2.0010(k-l)+e(k-2)=—0.75x(k)+1.50x(k-1)—0.7503x(k-2)(11)3，根据离散传递函数求系统状态空间表达式根据G(z)=®(zG(z)=®(z)X(z)(a+ab/2)z-1-abz-2/21-bz-1+z-2=—a+'ab) 7a———z—1—bz—2I2丿1—bz—1+z—2设W(z)= (12)1—bz—1+z—2W(z)=bz-1W(z)-z-2W(z)+X(z)®(z)=-aX(z)+a—岁丿z-1W(z)+”—字丿丿—2®(z)=-aX(z)+a—岁丿z-1W(z)+”—字丿丿—2W(z)13)选取状态变量X(z)=z-1W(z),X(z)=z-2W(z)=z-1X(z)12114)将W(z)代入X(z),X(z)，取z反变换，可得状态方程12x(k+1)x(k+1)2rx(汀Tx1(k)2+_0_-10x(k)15)(k)=-ax(k)x(k)1x(k)2将参数L=1,g=10,T=0.01代入(15)(k+1)「=(k+1)_=0(k)=-0.75x(k)-0.0003x(k)-0.0003x(k)124，稳定性分析对于开环系统，由传递函数可得系统的特征方程为A(z)=z2-2.001z+1式得该参数下的状态方程为x1x22.00081-10x(k)

x1(k)2x(k)16)17)特征方程的根为z=1.032,z=0.969。12由于特征根中有一个大于1，位于单位圆外，故系统是不稳定的。四，PID控制器设计由以上分析可知，系统是不稳定的，为使倒立摆在受到干扰时能保持稳定必须对系统进行PID控制器的设计。

倒立摆计算机控制系统的框图如下倒立摆控制系统结构图模拟PID控制器的基本算式为x(k)=Ke(k)+K丈e(j)+Ke(k)-e(k-1) (18)PIDj=0根据式(16)和(18)可求得在不同K,K,K下的摆角值。具体的仿真过

PID程如下：1,取K=-2.0,K=-0.3固定不变，K分别取0.5,1.0,1.5,2.0,倒立摆的输出变I D P化曲线如下：KP=0.5,KI=-2.0,KD=-0.3(a1)KP=1.0,KI=-2.0,KD=-0.3t/s(b1)KP=1.5,KI=-2.0,KD=-0.3t/s(c1)由图（al）、（bl）、（cl）、（dl）可知，当K=0.5时，倒立摆在受到干扰时可以P达到稳定；继续增大K到1.0,系统的调节时间变短；当K=1.5时，虽然倒立PP摆仍能稳定，但调节时间变长；当K=2.0时，倒立摆已无法恢复稳定。P2,取K=1.0,K=—0.3固定不变，K分别取-0.5,-1.5,-2.5,-3.5，倒立摆的输出变PDI化曲线如下：KP=1.0,KI=-0.5,KD=-0.3(a2)KP=1.0,KI=-1.5,KD=-0.3t/s(b2)KP=1.0,KI=-2.5,KD=-0.3t/s(c2)（d2）由图（a2）、（b2）、（c2）、（d2）可知，当K=-0.5时，倒立摆在受到干扰时开始I可以达到稳定，但当时间变长后无法达到稳定；继续减小K到-1.5，系统可以达I到稳定；当K=-2.5时，系统的调节时间变短；当K减小到-3.5时，倒立摆已II无法恢复稳定。3,取K=1.0,K=-2.0固定不变，K分别取-0丄-0.3,-0.5,-0.7，倒立摆的输出变PID化,曲线如下：KP=1.0,KI=-2.0,KD=-0.1(a3)KP=1.0,KI=-2.0,KD=-0.3t/s(b3)KP=1.0,KI=-2.0,KD=-0.5t/s(c3)（d3）由图（a3）、（b3）、（c3）、（d3）可知，当K从-0.1减小到-0.7时，倒立摆在受到I干扰时达到稳定的调节时间先减小后增大，到-0.7时系统已无法恢复稳定。五，改变杆长后的仿真结果将杆长由L=1改为L=0.5，系统的仿真结果如下：1,取K=—1.2,K=—0.1固定不变，K分别取020.4,0.8,1.0，倒立摆的输出变IDP化曲线如下：KP=0.2,KI=-1.2,KD=-0.1(a11)KP=0.4,KI=-1.2,KD=-0.1(b11)KP=0.8,KI=-1.2,KD=-0.1(c11)012345678910t/s（d11）由图（all）、（bll）、（ell）、（dll）可知，当K=0.2时，倒立摆在受到干扰时可P以达到稳定；继续增大K到0.4,系统的调节时间变短；当K=0.8时，虽然倒PP立摆仍能稳定，但调节时间变长；当K=l.0时，倒立摆已无法恢复稳定。P2,取K=0.4,K=-0.l固定不变，K分别取-l.0,-l.3,-l.6,-l.9,倒立摆的输出变PDI化曲线如下：KP=0.4,KI=-1.0,KD=-0.1(a12)KP=0.4,KI=-1.3,KD=-0.1(b12)KP=0.4,KI=-1.6,KD=-0.1(c12)（d12）由图（al2）、（bl2）、（cl2）、（dl2）可知，当K=—1.0时，倒立摆在受到干扰时I开始可以达到稳定，但当时间变长时会失去稳定；继续减小K到-1.3，系统可以I达到稳定；当K=—1.6时系统可以稳定但调节时间变长；当K减小到-1.9时，II倒立摆已无法恢复稳定。3，取K=0.4,K=—1.2固定不变，K分别取-0.01,-0.1,-0.2,-0.24，倒立摆的输出PID变化,曲线如下：KP=0.4,KI=-1.2,KD=-0.01(a13)KP=0.4,KI=-1.2,KD=-0.1(b13)KP=0.4,KI=-1.2,KD=-0.2(c13)(d13)由图(a13)、(bl3)、(cl3)、(dl3)可知，当k【从-0.01减小到-0.24时，倒立摆在受到干扰时达到稳定的调节时间先减小后增大，到-0.24时系统已无法恢复稳定。附录：1，PID设计过程的Matlab程序如下l=1;g0=l0;%重力加速度tt=0.0l;%采样周期a=3/(4*1);%代数anumGs=[-a00];%连续传递函数的分子denGs=[l0-a*g0];%连续传递函数的分母sysGs=tf(numGs,denGs);%系统连续传递函数[F,G,C,D]=ssdata(c2d(sysGs,tt,'zoh'));%求得ZOH下以tt为采样周期的状态空间(离散)t(l)=0;x=zeros(l,l000);%一行1000列0矩阵theta(l)=0;e=zeros(l,l000);%—行1000列0矩阵s=0;KP=0.5;KI=-2.0;KD=-0.3;%设置PID参数y=[0;0];fork=2:l000;t(k)=t(k-1)+tt;ifk==20delta=0.1;elsedelta=0;endtheta(k-1)=theta(k-1)+delta;e(k)=-theta(k-l);%PID公式的第一项s=s+e(k);%PID公式的第二项x(k)=KP*e(k)+KI*s+KD*(e(k)-e(k-1));%PID公式theta(k)=C*y+D*x(k);%初始输出y为零y=F*y+G*x(k);endfigure(1)plot(t,theta),gridtitle('KP=0.5,KI=-2.0,KD=-0.3')xlabel('t/s'),ylabel('\theta')2,为便于改变参数，可利用Matlab的GUI界面来方便的进行PID参数的设计，所设计的界面如下图所示：此时需将程序命令添加到Plot按钮下，即functionplot_push_Callback(hObject,eventdata,handles)l=str2double(get(handles.L_edit,'String'));tt=str2double(get(handles.T_edit,'String'));g0=str2double(get(handles.G_edit,'String'));KP=str2double(get(handles.KP_edit,'String'));KI=str2double(get(handles.KI_edit,'String'));KD=str2double(get(handles.KD_edit,'String'));a=3/(4*l);numGs=[-a00];denGs=[10-a*g0];sysGs=tf(numGs,denGs);[F,G,C,D]=ssdata(c2d(sysGs,tt,'zoh'));t(1)=0;x=zeros(1,1000);theta(1)=0;e=zeros(1,1000);s=0;y=[0;0];fork=2:1000;t(k)=t(k-1)+tt;ifk==20delta=0.1