单阶倒立摆控制系统_Matlab实验报告精品资料.doc

现代系统分析设计工具实验报告 学 院 班 级 学 号 姓 名 指导教师 成 绩 2011年6月目录摘要3第一部分 单阶倒立摆系统建模 4（一） 对象模型4（二）电动机、驱动器及机械传动装置的模型 6第二部分 单阶倒立摆系统分析7第三部分 单阶倒立摆系统控制11（一）内环控制器的设计11（二） 外环控制器的设计 14第四部分 单阶倒立摆系统仿真结果16系统的simulink仿真16摘要：该问题源自对于娱乐型”独轮自行车机器人”的控制，实验中对该系统进行系统仿真，通过对该实物模型的理论分析与实物仿真实验研究，有助于实现对独轮自行车机器人的有效控制。控制理论中把此问题归结为“一阶直线倒立摆控制问题”。另外，诸如机器人行走过程中的平衡控制、火箭发射中的垂直度控制、卫星飞行中的姿态控制、海上钻井平台的稳定控制、飞机安全着陆控制等均涉及到倒立摆的控制问题。实验中通过检测小车位置与摆杆的摆动角，来适当控制驱动电动机拖动力的大小，控制器由一台工业控制计算机（IPC）完成。实验将借助于“Simulink封装技术子系统”，在模型验证的基础上，采用双闭环PID控制方案，实现倒立摆位置伺服控制的数字仿真实验。实验过程涉及对系统的建模、对系统的分析以及对系统的控制等步骤，最终得出实验结果。仿真实验结果不仅证明了PID方案对系统平衡控制的有效性，同时也展示了它们的控制品质和特性。第一部分 单阶倒立摆系统建模（一） 对象模型由于此问题为”单一刚性铰链、两自由度动力学问题”，因此，依据经典力学的牛顿定律即可满足要求。如图1.1所示，设小车的质量为，倒立摆均匀杆的质量为，摆长为，摆的偏角为，小车的位移为，作用在小车上的水平方向上的力为，为摆杆的质心。图1.1 一阶倒立摆的物理模型根据刚体绕定轴转动的动力学微分方程，转动惯量与角加速度乘积等于作用于刚体主动力对该轴力矩的代数和，则1） 摆杆绕其重心的转动方程为 （1-1）2）摆杆重心的水平运动可描述为 （1-2）3） 摆杆重心在垂直方向上的运动可描述为 （1-3）4） 小车水平方向运动可描述为 （1-4）由式（1-2）和式（1-4）得 （1-5） 由式（1-1）、式（1-2）和式（1-3）得 （1-6）整理式（1-5）和式（1-6），得 （1-7）因为摆杆是匀质细杆，所以可求其对于质心的转动惯量。因此设细杆摆长为，单位长度的质量为，取杆上一个微段，其质量为，则此杆对于质心的转动惯量有 杆的质量为 所以此杆对于质心的转动惯量有 由式（2-20）可见，一阶直线倒立摆系统的动力学模型为非线性微分方程组。为了便于应用经典控制理论对该控制系统进行设计，必须将其简化为线性定常的系统模型。若只考虑在其工作点附近（）的细微变化，则可近似认为 在这一简化思想下，系统精确模型式（1-7）可简化为 若给定一阶直线倒立摆系统的参数为：小车的质量；倒摆振子的质量；倒摆长度；重力加速度取，则可得到进一步的简化模型为 (1-8)上式为系统的”微分方程模型”，对其进行拉普拉斯变换可得系统的传递函数模型为 （1-9）（二）电动机、驱动器及机械传动装置的模型假设：选用日本松下电工MSMA021型小惯量交流伺服电动机，其有关参数如下：驱动电压：U=0100V 额定功率：=200W额定转速：n=3000r/min 转动惯量：额定转矩： 最大转矩：电磁时间常数：=0.001s 电机时间常数：=0.003s经传动机构变速后输出的拖动力为：F=016N；与其配套的驱动器为：MSDA021A1A，控制电压： =（0±10）V。 若忽略电动机的空载转矩和系统摩擦，就可以认为驱动器和机械传动装置均为纯比例环节，并假设这两个环节的增益分别为和。对于交流伺服电动机，其传递函数可近似为 由于是小惯性的电动机，其时间常数、相对都很小，这样可以进一步将电动机模型近似等效为一个比例环节。综上所述，电动机。驱动器。机械传动装置三个环节就可以合成为一个比例环节 第二部分 单阶倒立摆系统分析尽管上述数学模型系经 机理建模得出，但其准确性（或正确性）还需运用一定的理论方法加以验证，以保证以其为基础的仿真实验的有效性。采用仿真实验的方法在MATLAB的Simulink图形仿真环境下进行模型验证试验，其原理如图1.2所示。其中，上半部分为精确模型仿真图，下半部分为简化模型仿真图。图1.2 模型验证原理图利用Simulink压缩子系统功能可将验证原理图更加简捷的表示为图1.3所示的形式。其中，由得到的精确模型和简化模型的状态方程，可得到Fcn、Fcn1、Fcn2和Fcn3的函数形式为(0.12*u1+0.036*sin(u3)*power(u2,2)-0.9*sin(u3)*cos(u3)/(0.24-0.09*power(cos(u3),2)(0.3*cos(u3)*u1+0.09*sin(u3)*cos(u3)*power(u2,2)-6*sin(u3)/(0.09*power(cos(u3),2)-0.24) 0.8*u1-0.6*u3 40*u3-2.0*u1图1.3 利用子系统封装后的框图假定使倒立摆在（）初始状态下突加微小冲击力作用，则依据经验知，小车将向前移动，摆杆降倒下。下面利用仿真实验来验证正确数学模型的这一必要性质。编制绘图子程序：% Inverted pendulum% Model test in open loop% Singnals recuperation% 将导入到xy.mat中的仿真实验数据读出load xy.matt=signals(1,:); % 读取时间信号f=signals(2,:); % 读取作用力F信号x=signals(3,:); % 读取精确模型中的小车位置信号q=signals(4,:); % 读取精确模型中的倒立摆摆角信号xx=signals(5,:); % 读取简化模型中的小车位置信号qq=signals(6,:); % 读取简化模型中的倒立摆摆角信号% Drawing control and x (t) response signals% 画出在控制力的作用下的系统响应曲线% 定义曲线的横纵坐标、标题、坐标范围和曲线的颜色等特征figure (1) % 定义第一个图形hf=line (t,f(:); % 连接时间-作用力曲线grid on;xlabel (Time (s) % 定义横坐标ylabel (Force （N）) % 定义纵坐标axis (0 1 0 0.12) % 定义坐标范围axet=axes (position，get (gca,position),XAxisLocation,bottom,YAxisLocation,right,Color,None,XColor,k,YColor,k）; % 定义曲线属性ht=line (t, x,color,r,parent, axet) ; % 连接时间-小车位置曲线ht=line (t, xx,color,r,parent, axet); % 连接时间-小车速度曲线ylabel (Evolution of the x position (m) % 定义坐标名称axis (0 1 0 0.1) % 定义坐标范围title (Response x and x in the meter to a f (t) pulse of 0.1 N) % 定义曲线标题名称gtext (leftarrow f (t), gtext (x (t) rightarrow), gtext(leftarrow x (t)% drawing control and theta (t) response singalsfigure (2)hf=line (t, f (:);grid onxlabel (Time)ylabel (Force in N)axis (0 1 0 0.12)axet=axes (Position, get (gca,Position),XAxisLocation,bottom,YAxisLocation,right,Color,None,XColor,k,YColor,k);ht=line (t, q,color,r,parent,axet);ht=line (t, qq,color,r,parent,axet);ylabel (Angle evolution (red)axis (0 1 -0.3 0)title（response theta(t) and theta(t) in rad to a f(t) pulse of 0.1 N）gtext (leftarrow f(t),gtext (theta（t）rightarrow),gtext (leftarrow theta(t)执行该程序的结果如图1.4所示。从中可见，在0.1N的冲击力作用下，摆杆倒下（由零逐步增大），小车位置逐渐增加，这一结果符合前述实验设计，故可以在一定程度上确认该“一阶倒立摆系统”的数学模型是有效的。同时，由图中也可看出近似模型在0.8s以前与精确模型非常接近，因此，也可以认为近似模型在一定条件下可以表述原系统模型的性质。图1.4 模型验证仿真结果第三部分 单阶倒立摆系统控制由于一阶倒立摆系统位置伺服控制的核心是“在保证摆杆不倒的条件下，使小车位置可控”。因此，依据负反馈闭环控制原理，将系统小车位置作为“外环”，而将摆杆摆角作为“内环”，则摆角作为外环内的一个扰动，能够得到闭环系统的有效抑制（实现其直立不倒的自动控制）。 剩下的问题就是如何确定控制器（校正装置）的结构和参数。（一）内环控制器的设计 图1.5 反馈校正控制的系统内环框图其中，Ks=1.6为伺服电动机与减速机构的等效模型1.控制器的选择内环系统未校正时的传递函数为对于内环反馈控制器可有PD，PI，PID三种可能的结构形式，这里，不妨采用绘制各种控制器结构下“系统根轨迹”的办法加以分析比较，从之选出一种比较适合的控制器结构。各种控制器的开环传函的传递函数分别为：在Matlab下输入以下程序用“凑试”的方法画根轨迹图：num=分子;den=分母;xlabel('Real Axis');ylabel('Imag Axis');axis(横、纵坐标范围);title('Root Locus');grid;rlocus(num,den)下图为各种控制器下的系统根轨迹： （a） PD （b） PD （c）PI d） PID 从根轨迹不难发现，采用PD结构的反馈控制器，结构简单且可保证闭环系统的稳定。所以，选定反馈控制器的结构为PD形式的控制器。2.控制器参数的选定首先暂定。这样可以求出内环的传递函数为： 3.系统内环的simulink仿真及结果仿真结果为：（二） 外环控制器的设计可见，系统开环传递函数可视为一个高阶（4阶）且带有不稳定零点的“非最小相位系统”，为了便于设计，需要首先对系统进行一些简化处理。.系统外环模型的降阶（1）对内环等效闭环传递函数的近似处理 将高次项忽略，有近似条件可由频率特性导出，即由（2）得： （2）对象模型G1(s)的近似处理 由（3）得：由（4）得： ，所以，有取 再由“典型型”系统Bode图特性（ ）知： .用simulink对小车的位置在阶跃信号输入下的响应进行仿真：系统框图为：仿真结果：倒立摆位置在阶跃信号下的响应第四部分 单阶倒立摆系统仿真结果系统的simulink仿真连接图如下：仿真结果为：倒立摆在阶跃信号下摆杆和小车位置的响应从图中容易可以看出建立的一阶倒立摆控制系统在matlab中能够实现倒立摆的要求，能通过电动机牵引机构带动小车的移动来控制摆杆和保持平衡。为了进一步验证在不同摆杆下的，该一阶倒立摆控制系统是否还具有鲁棒特性可以分别取摆杆不同的质量和摆长，进行simulink仿真！由上图可知，建立的一阶倒立摆模型在不同摆长下能实现要求。但摆长不能过长！同理，建立的一阶倒立摆模型在不同质量的摆杆下能也实现要求，但同样不能过重！附录资料：不需要的可以自行删除各类滤波器的MATLAB程序一、 理想低通滤波器IA=imread('lena.bmp');f1,f2=freqspace(size(IA),'meshgrid');Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r>0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 二、理想高通滤波器IA=imread('lena.bmp');f1,f2=freqspace(size(IA),'meshgrid');Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r<0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 三、 Butterworth低通滤波器IA=imread('lena.bmp');f1,f2=freqspace(size(IA),'meshgrid');D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 四、 Butterworth高通滤波器IA=imread('lena.bmp');f1,f2=freqspace(size(IA),'meshgrid');D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=(D*D)/r(i,j); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 五、 高斯低通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');f1,f2=freqspace(size(IA),'meshgrid');D=100/size(IA,1);r=f1.2+f2.2;Hd=ones(size(IA);for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 六、 高斯高通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');f1,f2=freqspace(size(IA),'meshgrid');%D=100/size(IA,1);D=0.3;r=f1.2+f2.2;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1-exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 七、 梯形低通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');f1,f2=freqspace(size(IA),'meshgrid');%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=zeros(size(IA);Hd(r<D0)=1;for i=1:size(IA,1) for j=1:size(IA,2) if r(i,j)>=D0 & r(i,j)<=D1 Hd(i,j)=(D1-r(i,j)/(D1-D0); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 八、 梯形高通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');f1,f2=freqspace(size(IA),'meshgrid');%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=ones(size(IA);Hd(r<D1)=0;for i=1:size(IA,1) for j=1:size(IA,2) if r(i,j)>=D0 & r(i,j)<=D1 Hd(i,j)=(D0-r(i,j)/(D0-D1); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 九、 用其他方法编写的理想低通、理想高通、Butterworth低通、同态滤波程序1、 理想低通i1=imread('lena.bmp');i2=imnoise(i1,'salt & pepper',0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);d0=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); if d<=d0 h=1; else h=0; end y(i,j)=g(i,j)*h; endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1);subplot(2,2,2),imshow(uint8(i2);subplot(2,2,3),imshow(uint8(E2); 2、 理想高通i1=imread('lena.bmp');i2=imnoise(i1,'salt & pepper',0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);n=2;d0=10;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); if d<=d0 h=0; else h=1; end y(i,j)=g(i,j)*h;endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1);subplot(2,2,2),imshow(uint8(i2);subplot(2,2,3),imshow(uint8(E2); 3、 Butterworth低通i1=imread('lena.bmp');i2=imnoise(i1,'salt & pepper',0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);n=2;d0=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); h=1/(1+(d/d0)(2*n); y(i,j)=g(i,j)*h; endendy=ifftshift(y);E1=ifft2(y);E2=real(E1); figuresubplot(2,2,1),imshow(uint8(i1);subplot(2,2,2),imshow(uint8(i2);subplot(2,2,3),imshow(uint8(E2); 4、 同态滤波I=rgb2gray(imread('fabric00.bmp');M,N=size(I);T=double(I);L=log(T);F=fft2(L);A=2;B=0.3;for i=1:M for j=1:N D(i,j)=(i-M/2)2+(j-N/2)2); endendc=1.1;%锐化参数D0=max(M,N);H=(A-B)*(1-exp(c*(-D/(D02)+B;F=F.*H;F=ifft2(F);Y=exp(F);figuresubplot(1,2,1),imshow(I);subplot(1,2,2),imshow(uint8(real(Y); 十、 Gabor滤波器- 30 -