延时系统建模，整数延时与分数延时，连续传函与离散传函，Pade近似与Thiran近似，Matlab实现

连续传递函数

严格建模：指数形式

  根据拉普拉斯变换的性质，
$$\left[f\left(t\right)\leftrightarrow F\left(s\right)\right]\iff \left[f\left(t-t_0\right)\leftrightarrow e^{-s{t}_0}F\left(s\right)\right]$$

  因此，在连续传递函数中，可以用下式对延时环节进行建模，
$$H\left(s\right)=e^{-s{t}_0}$$
  以$t_0=0.001$为例，绘制波特图如下所示：

对数坐标系（幅值、相位）	自然坐标系（相位响应）

  对应Matlab代码如下：

clc; clear; close all;T = 0.001;w = linspace(1,6280,6280);

% 绘制对数坐标系波特图
s = tf("s");Ca = exp(-s*T);figure;bode(Ca, w);grid on;xlim([1,6280]);

% 绘制相位响应（线性坐标系）
[mag, phase, wout] = bode(Ca, w);figure;plot(wout,squeeze(phase)); 
xlabel('Frequency (rad/s)');ylabel('Phase (degrees)');
title('Phase Response (Linear Scale)');axis([1,6280,-360,0]);grid on;

近似建模：Pade近似

  对于不允许使用指数表达式的场合，采用Pade近似将指数形式表示为分式形式，不同阶次的计算结果如下。Pade近似的原理介绍可以参考Pade近似，原理与Matlab实现。
$$\{\begin{array}{rl}原始函数& :f(x)=e^x\\ [1,1]近似& :-\frac{x+2}{x-2}\\ [2,2]近似& :+\frac{x^2+6x+12}{x^2-6x+12}\\ [3,3]近似& :-\frac{x^3+12x^2+60x+120}{x^3-12x^2+60x-120}\end{array}\to \{\begin{array}{rl}延时传函& :H\left(s\right)=e^{-s{t}_0}\\ [1,1]近似& :\frac{-t_{0}s+2}{t_{0}s+2}\\ [2,2]近似& :\frac{{t_0}^2{s}^2-6t_{0}s+12}{{t_0}^2{s}^2+6t_{0}s+12}\\ [3,3]近似& :\frac{-{t_0}^3{s}^3+12{t_0}^2{s}^2-60t_{0}s+120}{{t_0}^3{s}^3+12{t_0}^2{s}^2+60t_{0}s+120}\end{array}$$

  不同阶次近似的幅值响应无差别（均为0dB），对比相位响应如下：

  对应的Matlab代码如下所示：

clc;clear;close all;T = 0.001;w = linspace(1,6280,6280);
s = tf("s");               Ca = exp(-s*T);         [mag, phase, wout] = bode(Ca,w);
[num1,den1] = pade(T,1);   Ca1 = tf(num1,den1);    [mag1, phase1, wout1] = bode(Ca1,w);
[num2,den2] = pade(T,2);   Ca2 = tf(num2,den2);    [mag2, phase2, wout2] = bode(Ca2,w);
[num3,den3] = pade(T,3);   Ca3 = tf(num3,den3);    [mag3, phase3, wout3] = bode(Ca3,w);

plot(w,squeeze(phase),w,squeeze(phase1)-360,w,squeeze(phase2)-360,w,squeeze(phase3)-720);
legend("指数形式","一阶近似","二阶近似","三阶近似");title('Phase Response (Linear Scale)');
xlabel('Frequency (rad/s)');ylabel('Phase (degrees)');axis([1,6280,-360,0]);grid on;

离散传递函数

整数延时建模

  根据$z$变换的性质，$z$与$s$之间存在如下映射：
$$\left[z\leftrightarrow e^{sT}\right]\iff \left[z^{-1}\leftrightarrow e^{-sT}\right]\iff \left[z^{-k}\leftrightarrow e^{-skT}\right]$$
  以延时$t_0=2T=0.002$为例，对比波特图：

连续传函波特图	离散传函波特图

  对应的Matlab代码为：

clc;clear;close all;T = 0.001;k = 2;w = linspace(1,6280,6280);
s = tf("s");          Ca = exp(-s*k*T);    figure();    bode(Ca,w);grid on
z = tf("z",0.001);    Cd = z^-k;           figure();    bode(Cd,w);grid on

小分数延时建模

  对于延时$t_0<T$的小分数延时情况，有两种近似方案：(1) Pade近似+离散化；(2) Thiran近似离散化方法。 对比结果如下所示，对于小分数延时情况，两种方法所得结果是重叠的。Thiran近似的原理介绍可以查阅期刊原文，具体见参考文献部分。

  对应Matlab代码如下：

clc;clear;close all;T = 0.001;k = 0.4;w = linspace(1,6280,6280);
s = tf("s");    Ca = exp(-s*T*k); 

Cd1 = c2d(pade(Ca,1),T,"tustin")
Cd2 = thiran(k*T,T)

[mag0, phase0, wout0] = bode(Ca,w);    [mag1, phase1, wout1] = bode(Cd1,w);    [mag2, phase2, wout2] = bode(Cd2,w);

plot(w, squeeze(phase0),w,squeeze(phase1),w,squeeze(phase2));
xlim([0,2000]);legend("Continuous","Pade","Thiran"); grid on;

  需要说明的是，提高Pade近似的阶数，并不能改善离散化的近似效果，甚至表现出恶化作用。这是因为误差来源包括两部分，Pade近似带来的误差和Tustin近似带来的误差，一阶Pade与一阶Tustin可以实现最大限度的误差抵消。

  对应Matlab代码如下：

clc;clear;close all;T = 0.001;k = 0.4;w = linspace(1,6280,6280);s = tf("s");   Ca = exp(-s*T*k);

Cd1 = c2d(pade(Ca,1),T,"tustin");    Cd2 = c2d(pade(Ca,2),T,"tustin");    Cd3 = c2d(pade(Ca,3),T,"tustin");
[mag0, phase0, wout0] = bode(Ca,w);     [mag1, phase1, wout1] = bode(Cd1,w);
[mag2, phase2, wout2] = bode(Cd2,w);    [mag3, phase3, wout3] = bode(Cd3,w);
plot(w, squeeze(phase0),w,squeeze(phase1),w,squeeze(phase2)-360,w,squeeze(phase3)-360);grid on;
xlabel('Frequency (rad/s)');ylabel('Phase (degrees)');legend("连续传函","一阶Pade","二阶Pade","三阶Pade");

大分数延时建模

  对于延时$t_0>T$的大分数延时情况，有四种近似方案：(1) 四舍五入取整；(2) Pade近似+离散化；(3) Thiran近似离散化方法；(4) 整数部分严格建模，分数部分Pade近似+离散化；(5) 整数部分严格建模，分数部分Thiran近似离散化。 对比结果如下所示，大分数延时的近似效果为：方案3 > 方案4 = 方案5 > 方案2 > 方案1。

  对应Matlab代码如下：

clc;clear;close all;T = 0.001;k = 1.5;w = linspace(1,6280,6280);s = tf("s");z = tf("z",T);               
CaAll = exp(-s*T*k); CaInt = exp(-s*T*1); CaFra = exp(-s*T*0.5);

Cd0 = z^-2
Cd1 = c2d(pade(CaAll,1),T,"tustin")
Cd2 = thiran(k*T,T)
Cd3 = z^-1*c2d(pade(CaFra,1),T,"tustin")
Cd4 = z^-1*thiran(0.5*T,T)

[mag, phase, wout] = bode(CaAll,w);     [mag0, phase0, wout0] = bode(Cd0,w);    [mag1, phase1, wout1] = bode(Cd1,w);
[mag2, phase2, wout2] = bode(Cd2,w);    [mag3, phase3, wout3] = bode(Cd3,w);    [mag4, phase4, wout4] = bode(Cd4,w);

figure();plot(w,squeeze(phase),'-r','linewidth',2);hold on;
plot(w,squeeze(phase0),'g');            plot(w,squeeze(phase1)-360,'b');
plot(w,squeeze(phase2)-360,'m');        plot(w,squeeze(phase3),'y','linewidth',2);
plot(w,squeeze(phase4),'k');            xlim([0,2000]);
xlabel('Frequency (rad/s)');ylabel('Phase (degrees)');legend("连续模型","方案1","方案2","方案3","方案4","方案5");grid on;

延时系统建模汇总

• 当延时系统建模与分析局限于连续传函，尽量避免近似过程，按照严格指数形式保留准确性；
• 当延时系统建模与分析局限于连续传函，但必须进行有理化处理时，建议采用高阶帕德近似；
• 当连续传函最终需要数字化实现时，整倍延时采用严格离散化方案（一阶Pade+Tustin离散化，等价于，严格离散化）；
• 当连续传函最终需要数字化实现时，小分数延时采用Pade近似+离散化方案，采用一阶近似减小误差（等价于Thiran方案）；
• 当连续传函最终需要数字化实现时，大分数延时采用方案三（Thiran整体离散化，适用于高要求场合）或方案四（整数严格+分数Pade离散，适用于大多数场景）；

参考文献

[1] 余成波. 拉普拉斯变换的性质, 信号与系统.第2版[M].清华大学出版社,2007.
[2] Pade近似，原理与Matlab实现。
[3] Matlab文档：Padé approximation of models with time delay.
[4] Matlab文档：Bode frequency response of dynamic system.
[5] Matlab文档：Generate fractional delay filter based on Thiran approximation.
[6] T. Laakso, V. Valimaki, “Splitting the Unit Delay”, IEEE Signal Processing Magazine, Vol. 13, No. 1, p.30-60, 1996.