《先进 PID 控制及其 MATLAB 仿真》 刘金琨 著 电 子 工 业 出 版 社 压缩文件里包含了这本书的电子档和全书的程序。一章一个文件夹。

				function dx=Model(t,x)
				global w A alfa J Ce R Km Ku S a1 Fm Fc M
				persistent aa
				dx=zeros(3,1);
				
				a1=1.0;  %Effect on the shape of friction curve
				Fm=50;
				Fc=1.5;
				
				F=-J*x(3);
				if t==0
				   aa=0;
				end
				dF=-J*aa;
				
				if abs(x(2))				   if F>Fm
				      Ff=-Fm;
				      dFf=0;
				   elseif F				      Ff=Fm;
				      dFf=0;
				   else
				      Ff=-F;
				      dFf=dF;
				   end      
				end
				if x(2)>alfa   
				   Ff=Fc+(Fm-Fc)*exp(-a1*x(2));
				   dFf=(Fm-Fc)*exp(-a1*x(2))*(-a1)*x(3);
				elseif x(2)				   Ff=-Fc-(Fm-Fc)*exp(a1*x(2));
				   dFf=-(Fm-Fc)*exp(a1*x(2))*a1*x(3);
				end
				
				if S==1
				   rin=A*sin(w*t);
				   drin=A*w*cos(w*t);
				   ddrin=-A*w*w*sin(w*t);
				   dddrin=-A*w*w*w*cos(w*t);
				end
				if S==2
				   rin=1;
				   drin=0;
				   ddrin=0;
				   dddrin=0;
				end
				if S==3
				   rin=A*sign(sin(0.4*2*pi*t));
				   drin=0;
				   ddrin=0;
				   dddrin=0;
				end
				error=rin-x(1);
				derror=drin-x(2);
				dderror=ddrin-x(3);
				
				u=200*error+40*derror;   %PID
				du=200*derror+40*dderror;
				
				if u>=10   
				   u=10;
				end
				if u				   u=-10;
				end
				
				if M==0
				   Ff=0;dFf=0;   %No Friction
				end
				dx(1)=x(2);
				dx(2)=-Km*Ce/(J*R)*x(2)+Ku*Km*u/(J*R)-Ff/J;
				dx(3)=-Km*Ce/(J*R)*x(3)+Ku*Km*du/(J*R)-dFf/J;
				aa=dx(3);