现代控制--倒立摆系统仿真

第1页目录倒立摆系统的构成............................2单级倒立摆数学模型的建立....................2传递函数......................................................................4状态空间方程..............................................................5系统MATLAB仿真和开环响应.....................................5稳定性与可控性分析..........................9控制器设计.................................10基于状态反馈的控制算法设计与仿真LQR...........10极点配置法................................................................14PID控制算法.............................................................17实验结果及与仿真结果的对比分析.............27感想和建议.................................28第2页倒立摆系统的构成图1倒立摆系统的组成框图如图1所示为倒立摆的结构图。系统包括计算机、运动控制卡、伺服机构、倒立摆本体和光电码盘几大部分，组成了一个闭环系统。光电码盘1将小车的位移、速度信号反馈给伺服驱动器和运动控制卡，摆杆的位置、速度信号由光电码盘2反馈回控制卡。计算机从运动控制卡中读取实时数据，确定控制决策（小车向哪个方向移动、移动速度、加速度等），并由运动控制卡来实现该控制决策，产生相应的控制量，使电机转动，带动小车运动，保持摆杆平衡。单级倒立摆数学模型的建立在忽略了空气流动，各种摩擦之后，可将倒立摆系统抽象成小车和匀质杆组成的系统，如下图2所示图2单级倒立摆模型示意图那我们在本实验中定义如下变量：M小车质量（本实验系统0.5Kg）m摆杆质量（本实验系统0.2Kg）b小车摩擦系数（本实验系统0.1N/m/sec）l摆杆转动轴心到杆质心的长度（0.3m）第3页I摆杆惯量（0.006kg*m*m）F加在小车上的力x小车位置φ摆杆与垂直向上方向的夹角θ摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下）下面我们对这个系统作一下受力分析。下图3是系统中小车和摆杆的受力分析图。其中，N和P为小车与摆杆相互作用力的水平和垂直方向的分量。注意：在实际倒立摆系统中检测和执行装置的正负方向已经完全确定，因而矢量方向定义如图，图示方向为矢量正方向。图3倒立摆模型受力分析分析小车水平方向所受的合力，可以得到等式：应用Newton方法来建立系统的动力学方程过程如下：分析小车水平方向所受的合力，可以得到以下方程：NxbFxM由摆杆水平方向的受力进行分析可以得到下面等式：)sin(22lxdtdmN即sincos2mlmlxmN把这个等式代入上式中，就得到系统的第一个运动方程：FmlmlxbxmMsincos)(2（1）为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得到下面方程：第4页cossin)cos(222mlmlmgPldtdmmgP即：力矩平衡方程如下：INlPlcossin注意：此方程中力矩的方向，由于sinsin,coscos,，故等式前面有负号。合并这两个方程，约去P和N，由231mlI得到第二个运动方程：cossin234xmlmglml（2）设（是摆杆与垂直向上方向之间的夹角），假设与1（单位是弧度）相比很小，即《1，则可以进行近似处理：1cos，sin，0)(2dtd。用u来代表被控对象的输入力F，线性化后两个运动方程如下：umlxbxmMxgl)(34（3）传递函数对方程组(3)进行拉普拉斯变换，得到)()()()()()()()(342222sUssmlssbXssXmMssXsgssl（4）注意：推导传递函数时假设初始条件为0。由于输出为角度，求解方程组（4）的第一个方程，可以得到)(]34[)(2ssglsX把上式代入方程组（4）的第二个方程，得到)()()()()()()(22222sUssmlsssgmlmlIbsssgmlmlImM第5页整理后得到传递函数：sqbmglsqmglmMsqbmlssqmlsUs2323442)()()(其中])())([(22mlmlImMq状态空间方程系统状态空间方程为DuCXyBuAXX方程组（3）对,x解代数方程，得到解如下：ulmMlmMmMgxlmMbumMmMmgxmMbxxx)4(3)4()(3)4(3)4(4)4(3)4(4整理后得到系统状态空间方程：ulmMmMxxlmMmMglmMbmMmgmMbxx)4(30)4(400)4()(3)4(3010000)4(3)4(400010uxxxy0001000001系统Matlab仿真和开环响应实际系统参数如下，求系统的传递函数、状态空间方程，并进行脉冲响应和阶跃响应的第6页Matlab仿真。M小车质量1.096Kgm摆杆质量0.109Kgb小车摩擦系数0.1N/m/secl摆杆转动轴心到杆质心的长度0.25mI摆杆惯量0.0034kg*m*mF加在小车上的力x小车位置θ摆杆与垂直方向的夹角T采样频率0.005秒注意：在进行实际系统的Matlab仿真时，请将采样频率改为实际系统的采样频率。传递函数：在Matlab中，拉普拉斯变换后得到的传递函数可以通过计算并输入分子和分母矩阵来实现。求系统传递函数的m-文件内容如下：M=1.096;m=0.109;b=0.1;I=0.0034;g=9.8;l=0.25;q=(M+m)*(I+m*l^2)-(m*l)^2;%simplifiesinputnum=[m*l/q0]den=[1b*(I+m*l^2)/q-(M+m)*m*g*l/q-b*m*g*l/q]t=0:0.01:5;impulse(num,den,t)axis([01.1070])执行上面的文件，就可以求出系统传递函数的分子与分母多项式的Matlab表示：num=2.35660den=1.00000.0883-27.8285-2.3094第7页可以得到系统开环脉冲响应的曲线如下：00.20.40.60.81010203040506070ImpulseResponseTime(sec)Amplitude图4系统开环脉冲响应曲线状态空间法：状态空间法可以进行单输入多输出系统设计，(从实验二开始，我们将尝试同时对摆杆角度和小车位置进行控制)。为了更具挑战性，给小车一个阶跃输入信号，设计指标如下：小车位置x和摆杆角度θ的稳定时间小于5秒；位置x的上升时间小于0.5秒；摆杆角度的超调量小于20度（0.35弧度）。下面，我们用Matlab求出系统的状态空间方程各矩阵，并仿真系统的开环阶跃响应。在这里同样给出了一个m-文件，执行这个文件，Matlab将会给出系统状态空间方程的A，B，C和D矩阵，并可以绘出在给定输入为一个0.2m的阶跃信号时系统的响应曲线。M=1.096;m=0.109;b=0.1;I=0.0034;g=9.8;l=0.25;p=I*(M+m)+M*m*l^2;%denominatorfortheAandBmatricies第8页A=[0100;0-(I+m*l^2)*b/p(m^2*g*l^2)/p0;0001;0-(m*l*b)/pm*g*l*(M+m)/p0]B=[0;(I+m*l^2)/p;0;m*l/p]C=[1000;0010]D=[0;0]T=0:0.005:10;U=0.2*ones(size(T));[Y,X]=lsim(A,B,C,D,U,T);plot(T,Y)axis([02.50100])执行该m文件，可以求出系统的状态空间A、B、C、D矩阵,得到开环系统阶跃响应的曲线。A=01.0000000-0.08830.629300001.00000-0.235727.82850B=00.883202.3566C=10000010D=0第9页0100.511.522.50102030405060708090100图5倒立摆状态空间开环系统阶跃响应曲线图中，实线是摆杆角度响应曲线，虚线是小车位置响应曲线。稳定性与可控性分析我们先看一看系统的稳定性，将数据代入状态方程中，利用matlab程序可以求出系统的零极点。源代码如下：M=1.096;m=0.109;b=0.1;I=0.0034;l=0.25;[abcd]=wer_ss(M,m,b,l);%自己编写的函数，建立模型之用，具体程序见下面sysc=ss(a,b,c,d);sysd=c2d(sysc,0.005);[dadbdcdd]=ssdata(sysd);[zpgain]=ss2zp(da,db,dc,dd,1)z=-0.9999-0.99991.02751.0000+0.0000i0.97331.0000-0.0000ip=第10页1.00000.99961.02850.9723gain=1.0e-004*0.11130.3338wer_ss源程序：function[abcd]=wer_ss(M,m,b,l)a=[0100;0-4*b/(4*M+m)3*m*9.8/(4*M+m)0;0001;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=[1000;0010];d=[0;0]由得到的p（极点）可知，有的极点在单位圆外，所以可知原系统是不稳定的。同样，我们可以利用matlab来得到系统的能控性，源代码如下：ud=ctrb(da,db);rank(ud)ans=4由得到的rank（ud）的值可知，原系统的能控性矩阵为4，所以我们可知原系统是能控的。控制器设计基于状态反馈的控制算法设计与仿真LQR由理论分析知，可以设计基于最优控制的状态调节器，使系统闭环稳定。设状态反馈调节律的形式为)()(kxKkur通过使性能指标函数)]()()()()()([21)()(2110kukRkukXkQkXNPXNXJrTNkrTT为最小，根据在附录1中我们所介绍的求得GkPHHKPHRKTTrr)1(])1([1其中P由下列黎卡提方程获得第11页GkPHHkPHRHkPGGkPGQkPTTrTTr)1(])1([)1()1()(1其中rQ，rR分别用来对状态向量x(k)，控制向量u(k)引起的性能度量的相对重要性进行加权。在实际运算中我们运用Matlab控制系统工具箱中的“dlqr”函数直接进行运算。利用dlqr函数，我们需要提供两个权值矩阵：Q、R。通常我们取R=1，而对于Q我们只能通过不断的凑取来得到。源代码如下：R=1;Q=[10000;0000;0010;0000]Q=10000000000100000