matlab、simulink实现PID设计精品资料.doc

自动化专业07级计算机仿真与MATLAB课程报告题目： 基于MATLAB/Simulink的PID控制系统的设计与仿真班级：姓名： 学号：2010年6月 基于MATLAB/Simulink的PID控制系统的设计与仿真摘要: 介绍了基于Ziegler- Nichols整定方法的PID 控制器设计, 给出了基于MATLAB和Simulink的实现方法和仿真。仿真结果表明, 此算法设计的PID 控制器有良好的性能指标。1 控制对象建模1.1 PID 控制系统的建模PID(Proportional，Integral and Differemial)控制器是一种基于“过去”，“现在”和“未来”信息估计的简单算法。常规PID控制系统原理框图如下图所示，系统主要由PID控制器和被控对象组成。作为一种线性控制器，它根据给定值rin(t)与实际输出值yout(t)构成控制偏差e(t)，将偏差按比例、积分、和微分通过线性组合构成控制量u(t)，对被控对象进行控制。 PID控制系统原理图PID 控制器的数学描述为:其传递函数可表示为:PID控制器各校正环节的作用如下：1比例环节：成比例地反映控制系统的偏差信号e(t)，偏差一旦产生，控制器立即产生控制作用，以减少偏差。2积分环节：主要用于消除静差，提高系统的无差度。积分作用的强弱取决于积分时间常数Ti，Ti越大，积分作用越弱，反之越强。3微分环节：反映偏差信号的变化趋势(变化速率)，并能在偏差信号变得太大之前，在系统中引入一个有效的早期修正信号，从而加快系统的动作速度，减少调节时间。从根本上讲，设计 PID 控制器也就是确定其比例系数K p、积分系数Ti 和微分系数Td , 这三个系数取值的不同，决定了比例、积分和微分作用的强弱。控制系统的整定就是在控制系统的结构已经确定、控制仪表和控制对象等处在正常状态的情况下, 适当选择控制器的参数使控制仪表的特性和控制对象的特性相配合, 从而使控制系统的运行达到最佳状态, 取得最好的控制效果。本文介绍基于Ziegler- Nichols整定方法的 PID 控制器设计。1.2 被控对象的建模在实际的过程控制系统中，有大量的对象模型可以近似地由带有延迟的一阶传递函数模型来表示，该对象的模型可以表示如下：如果不能建立起系统的物理模型，可通过试验测取对象模型的阶跃响应，从而得到模型参数。当然, 我们也可在已知对象模型的情况下, 由 MATLAB 通过 STEP( ) 函数得到对象模型的开环阶跃响应曲线。在被控对象的阶跃响应输出信号图(如图所示)中, 可获取 K、L 和 T 参数。2 PID控制系统的设计Ziegler- Nichols法是一种基于频域设计 PID 控制器的方法。此法首先通过实验获取控制对象单位阶跃响应，获得K、L 和 T 参数。令a=KL/T,我们可以通过下表给出的Ziegler- Nichols经验公式确定P、PI 和 PID 控制器的参数。控制器类型KpTiTdP0PI0PIDZiegler- Nichols法整定控制器参数3 PID 控制系统MATLAB/Simulink仿真分析3.1 在MATLAB 下实现PID 控制器的设计与仿真根据Ziegler- Nichols法，这里编写一个MATLAB函数ziegler，该函数的功能实现由Ziegler- Nichols公式设计PID 控制器，在设计过程中可以直接调用。其源程序如下：function Gc,Kp,Ti,Td,H=ziegler(key,vars)Ti=; Td=; H=1;if length(vars)=4,K=vars(1); L=vars(2); T=vars(3); N=vars(4); a=K*L/T;if key=1, Kp=1/a;elseif key=2, Kp=0.9/a; Ti=3.33*L;elseif key=3 | key=4, Kp=1.2/a; Ti=2.2*L; Td=L/2; endelseif length(vars)=3,K=vars(1); Tc=vars(2); N=vars(3);if key=1, Kp=0.5*K;elseif key=2, Kp=0.4*K; Ti=0.8*Tc;elseif key=3 | key=4, Kp=0.6*K; Ti=0.5*Tc; Td=0.12*Tc; endelseif length(vars)=5,K=vars(1); Tc=vars(2); rb=vars(3); N=vars(5);pb=pi*vars(4)/180; Kp=K*rb*cos(pb);if key=2, Ti=-Tc/(2*pi*tan(pb);elseif key=3|key=4, Ti=Tc*(1+sin(pb)/(pi*cos(pb); Td=Ti/4; endendswitch key case 1, Gc=Kp; case 2, Gc=tf(Kp*Ti,1,Ti,0); case 3, nn=Kp*Ti*Td*(N+1)/N,Kp*(Ti+Td/N),Kp; dd=Ti*Td/N,1,0; Gc=tf(nn,dd);end该函数的 调用格式为：Gc,Kp,Ti,Td=Ziegler(key,vars)其中，key为选择控制器类型的变量：当key=1,2,3时分别表示设计P、PI、PID控制器；若给出的是阶跃响应数据，则变量vars=K,L,T,N;若给出的是频域响应数据，则变量vars=Kc,Tc,N。在实际的过程控制系统中，有大量的对象模型可以近似地由带有延迟的一阶传递函数模型来表示，该对象的模型可以表示如下：这里我们不妨设K=8，T=360,L=180，则对象模型可以表示为：利用ziegler()函数计算系统P、PI、PID控制器的参数，并给出校正后系统阶跃响应曲线。程序如下：K=8;T=360;L=180;num=K;den=T 1;G1=tf(num,den)np,dp=pade(L,2);Gp=tf(np,dp)figure,step(G1*Gp);title('未校正前系统阶跃响应曲线曲线');grid;Gc1,Kp1=ziegler(1,K,L,T,1);Gc1Gc2,Kp2,Ti2=ziegler(2,K,L,T,1);Gc2Gc3,Kp3,Ti3,Td3=ziegler(3,K,L,T,1);Gc3G_c1=feedback(G1*Gc1,Gp);figure,step(G_c1);title('P控制器校正后的系统阶跃响应曲线')grid;G_c2=feedback(G1*Gc2,Gp);figure,step(G_c2);title('PI 控制器校正后的系统阶跃响应曲线');grid;G_c3=feedback(G1*Gc3,Gp);figure,step(G_c3);title('PID控制器校正后的系统阶跃响应曲线');grid;figure,step(G_c1);hold on;step(G_c2);hold on;step(G_c3);title('P、PI、PID控制器校正后的系统阶跃响应曲线');grid;gtext('P')gtext('PI')gtext('PID')运行程序，输出如下：Gc1 = 0.2500 Transfer function:134.9 s + 0.225- 599.4 s （PI控制器的传递函数） Transfer function:2.138e004 s2 + 145.8 s + 0.3- 3.564e004 s2 + 396 s（PID控制器的传递函数）P、PI、PID校正后的系统阶跃响应曲线如下图：由上图可知，用Ziegler- Nichols公式计算P、PI、PID控制器对系统校正后，其阶跃响应曲线中的P、PI控制二者的响应速度基本相同，因为两种控制的比例系数不同，因此系统稳定的输出值不同。PI控制超调量比P控制的要小一些。PID控制比前者的响应速度都快，但超调量最大。 3.2 在Simulink 下实现PID 控制器的设计与仿真这里仍然设被控对象的传递函数是建立Simulink模型：图中，“Integrator”为积分器，“Derivative”为微分器，“Kp”为比例系数，“Ti”为积分时间常数，“Td”为微分时间常数。进行P控制器参数整定时，微分器和积分器的输出不连到系统中，在Simulink中，把微分器和积分器的输出连线断开即可。同理，进行PI控制器参数整定时，微分器的输出连线断开。Ziegler- Nichols整定的第一步是获取开环系统的单位阶跃响应，在Simulink中，把反馈连线、微分器的输出连线、积分器的输出连线都断开，“Kp”的值置为1，连线得：选定仿真时间，仿真运行，运行完毕后，双击“Scope”得到结果：根据Ziegler- Nichols经验公式，可知P控制整定时，比例放大系数Kp=0.25，将“Kp”的值置为0.25，并连上反馈连线，得：选定仿真时间，仿真运行，运行完毕后，双击“Scope”得到结果：上图即为P控制时系统的单位阶跃响应。根据Ziegler- Nichols经验公式，可知PI控制整定时，比例放大系数Kp=0.225，积分时间常数Ti=594,将“Kp”的值置为0.225,“1/Ti”的值为1/594，将积分器的输出连线连上，得：选定仿真时间，仿真运行，运行完毕后，双击“Scope”得到结果：上图即为PI控制时系统的单位阶跃响应。根据Ziegler- Nichols经验公式，可知PID控制整定时，比例放大系数Kp=0.3，积分时间常数Ti=396,微分时间常数Td=90,将“Kp”的值置为0.3,“1/Ti”的值为1/396，“Td”的值置为90，将微分器的输出连线连上，得：选定仿真时间，仿真运行，运行完毕后，双击“Scope”得到结果：由以上三图同样可以看出，P、PI控制二者的响应速度基本相同，但系统稳定的输出值不同。PI控制超调量比P控制的要小一些。PID控制比前者的响应速度都快，但超调量最大。针对该PID 控制器，我们可以通过外加扰动信号来测试其控制效果。如下图，我们在t=4000s时，外加一个幅值为15的扰动信号：将该扰动信号加到系统输入端，如下图：选定仿真时间，仿真运行，运行完毕后，双击“Scope”得到结果：当系统稳定后，若加一个扰动信号，PID控制器可以很快对被控对象的响应进行校正，使其尽快稳定。由上图可以看出，该PID 控制器效果良好。从系统接入PID 控制器前后的阶跃响应曲线中, 我们可以明显地看到系统性能的改善。利用MATLAB/Simulink可以实现PID 控制器的离线设计和整定, 并可实现实验室仿真。但是这种常规的PID 控制不具有自适应性, 在长期工作时对象参数会产生偏移, 系统具有时变不确定性, 也存在非线性, 工况点附近小范围的线性化假设在整个工作范围中不能成立时, 就难以达到理想的控制效果。为此, 我们可以考虑自适应的PID 控制算法。4 SummaryIn the ever-changing 21st century , the past decade has witnessed that automation has been one of the most important technology in engineering science ,and ,as the essential tool for automation , MATLAB has been apparently paying more attention meanwhile by an increasing number of teachers and students whether they major in engineering or not . As an illustration , in our class , students obviously enjoy doing the MATLAB assignment ,combining the practice with theories which we have learned . Personally , I believe this phenomenon should be approved and acclaimed , and it exactly tells the truth that the conventional view that the theory teaching should be one of the very highest priorities for automation is not entirely accurate . Furthermore , as I see, with proper reforms on automation education , programming based on MATLAB will be the latest rage among the passionate young . 附录资料：不需要的可以自行删除各类滤波器的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滤波器