基于FOLPD模型的PID控制器参数的整定

在PID控制器的诸多经典参数整算法中，绝大多数的算法都是在带有时间延迟的一阶模型（FOPDT）的基础上提出的，模型的一般形式为：
主要是因为大部分过程控制的受控对象模型的响应曲线和一阶模型的响应较为类似，可以直接进行拟合。
一般的过程控制对象的阶跃响应曲线形状如下图所示：

可根据读图、最优降阶方法0/1，最小二乘参数识别的方法求出k,L,T
下面给出一个编写的matlab函数用来求k，L，T

function [K,L,T,G1]=getfopdt(key,G)
switch key
case 1, [y,t]=step(G);
   fun=@(x,t)x(1)*(1-exp(-(t-x(2))/x(3))).*(t>x(2));
   x=lsqcurvefit(fun,[1 1 1],t,y); K=x(1); L=x(2); T=x(3);
case 2, [Kc,Pm,wc,wcp]=margin(G);
   ikey=0; L=1.6*pi/(3*wc); K=dcgain(G); T=0.5*Kc*K*L;
   if isfinite(Kc), x0=[L;T];
      while ikey==0, u=wc*x0(1); v=wc*x0(2);
         FF=[K*Kc*(cos(u)-v*sin(u))+1+v^2; sin(u)+v*cos(u)];
   J=[-K*Kc*wc*sin(u)-K*Kc*wc*v*cos(u),-K*Kc*wc*sin(u)+2*wc*v;
      wc*cos(u)-wc*v*sin(u), wc*cos(u)]; x1=x0-inv(J)*FF;
         if norm(x1-x0)<1e-8, ikey=1; else, x0=x1; end, end
      L=x0(1); T=x0(2);
   end
case 3, [n1,d1]=tfderv(G.num{1},G.den{1});
   [n2,d2]=tfderv(n1,d1); K1=dcgain(n1,d1);
   K2=dcgain(n2,d2); K=dcgain(G); Tar=-K1/K;
   T=sqrt(K2/K-Tar^2); L=Tar-T;
case 4
   Gr=opt_app(G,0,1,1);
   L=Gr.ioDelay;
   T=Gr.den{1}(1)/Gr.den{1}(2); K=Gr.num{1}(end)/Gr.den{1}(2);
end
G1=tf(K,[T 1],'iodelay',L);
function [e,f]=tfderv(b,a)
f=conv(a,a); na=length(a); nb=length(b);
e1=conv((nb-1:-1:1).*b(1:end-1),a);
e2=conv((na-1:-1:1).*a(1:end-1),b); maxL=max(length(e1),length(e2));
e=[zeros(1,maxL-length(e1)) e1]-[zeros(1,maxL-length(e2)) e2];

其中输入参数的key变量表示采用的各种方法，可以取值为1,2,3,4.对于已知的阶跃响应数据，通过响应曲线识别一阶模型，key=1，G为受控对象模型。函数返回值为一阶近似模型参数k，L，T和近似的传递函数模型G1.
输入变量key=2时表示采用的是基于频率响应的近似方法，
输入变量key=3时表示采用的是基于传递函数的辨识方法，
输入变量key=4时表示采用的是最优降阶方法。

Ziegler-Nichols经验公式

有了这些参数，如何来设计PID控制器呢，我们可以采用经验公式来设计PID控制器。由k，L，T可以计算出一个中间变量a，a=KL/T

这样就可以根据上表设计出P、PI、PID控制器。设计方法非常简单直观，根据此算法编写一个matlab函数帮助计算

function [Gc,Kp,Ti,Td]=ziegler(key,vars)
switch length(vars)
   case 3,
      K=vars(1); Tc=vars(2); N=vars(3);
      if key==1, Kp=0.5*K; Ti=inf; Td=0;
      elseif key==2, Kp=0.4*K; Ti=0.8*Tc; Td=0;
      elseif key==3, Kp=0.6*K; Ti=0.5*Tc; Td=0.12*Tc; end
   case 4
      K=vars(1); L=vars(2); T=vars(3); N=vars(4); a=K*L/T;
      if key==1, Kp=1/a; Ti=inf; Td=0;
      elseif key==2, Kp=0.9/a; Ti=3*L; Td=0;
      elseif key==3, Kp=1.2/a; Ti=2*L; Td=L/2; end
   case 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)); Td=0;
      elseif key==3, Ti=Tc*(1+sin(pb))/(pi*cos(pb)); Td=Ti/4;
end, end
Gc=pidstd(Kp,Ti,Td,N);

这里的key=1,2,3分别对应于P、PI、PID控制器，用户可以选择控制器类型。vars=[K,L,T,N]。

Chien-Hrones-Reswick参数整定算法

在实际应用中，传统的Ziegler-Nichols算法有各种各样的变形，其中的Chien-Hrones-Reswick（CHR）就是其中的一种改进。其允许带有较大的阻尼，以确保没有超调量的快速响应。分为有0%超调量和20%超调量。在CHR算法中直接使用了时间常数。

按此算法编写的matlab函数：

 function [Gc,Kp,Ti,Td]=chrpid(key,vars)
K=vars(1); L=vars(2); T=vars(3); N=vars(4); ov=vars(5)+1;
a=K*L/T; KK=[0.3,0.35,1.2,0.6,1,0.5; 0.7,0.6,1,0.95,1.4,0.47];
if key==1, Kp=KK(ov,1)/a; Ti=inf; Td=0;
elseif key==2, Kp=KK(ov,2)/a; Ti=KK(ov,3)*T; Td=0;
else, Kp=KK(ov,4)/a; Ti=KK(ov,5)*T; Td=KK(ov,6)*L; end
Gc=pidstd(Kp,Ti,Td,N);

该函数的调用格式为 [Gc,Kp,Ti,Td]=chrpid(key,vars)，其返回的变量和函数ziegler（）是完全一致的。key=1,2,3分别对应于P、PI、PID控制器，变量vars可以表示为vars=[k,L,T,N,Os]，Os=0对应于没有超调量的控制，Os=1对应于有20%超调量的控制。
假设对象模型为一个六阶的传递函数G（s）=1/(s+1)^6
可以用下门面的matlab语句设计出Ziegler-Nichols PID控制器与两种准则下的CHR控制器。

s=tf('s');
G=1/(s+1)^6;
 [k,L,T]=getfopdt(4,G);
 N=10;
[Gc1,Kp,Ti,Td]=ziegler(3,[k,L,T,N])
[Gc2,Kp,Ti,Td]=chrpid(3,[k,L,T,N,0]);
[Gc3,Kp,Ti,Td]=chrpid(3,[k,L,T,N,1]);

绘制出不同控制器下的闭环系统的阶跃响应曲线：

step(feedback(G*Gc1,1),'-',feedback(G*Gc2,1),'--',feedback(G*Gc3,1),':',50)

从曲线的结果可以看出，虽然没有超调量的设定点控制器响应速度较慢，但控制的效果还比较理想，优于经典的 Ziegler-Nichols算法。