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

连续传递函数

严格建模：指数形式

根据拉普拉斯变换的性质，
$\left[ {f\left( t \right) \leftrightarrow F\left( s \right)} \right] \Leftrightarrow \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代码如下：

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实现。
$\left\{ \begin{aligned} 原始函数&: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} + 12{x^2} + 60x + 120}}{{{x^3} - 12{x^2} + 60x - 120}}\\ \end{aligned} \right. \to \left\{ \begin{aligned} 延时传函&: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} - 6{t_0}s + 12}}{{{t_0}^2{s^2} + 6{t_0}s + 12}} \\ [3,3]近似&: \frac{{ - {t_0}^3{s^3} + 12{t_0}^2{s^2} - 60{t_0}s + 120}}{{ {t_0}^3{s^3}+ 12{t_0}^2{s^2} + 60{t_0}s + 120}} \\ \end{aligned} \right.$

不同阶次近似的幅值响应无差别（均为0dB），对比相位响应如下：
在这里插入图片描述
对应的Matlab代码如下所示：

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] \Leftrightarrow \left[ {{z^{ - 1}} \leftrightarrow {e^{ - sT}}} \right] \Leftrightarrow \left[ {{z^{ - k}} \leftrightarrow {e^{ - skT}}} \right]$
以延时 $t_0=2T=0.002$ 为例，对比波特图：

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

对应的Matlab代码为：

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代码如下：

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代码如下：

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代码如下：

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.