传递函数离散化方法：差分，双线性Tustin，预矫正双线性，Matlab实现

传递函数严格离散化

  傅里叶变换、拉普拉斯变换，以及$Z$变换的定义如下所示：
$$傅里叶变换:\{\begin{array}{l}F\left(\omega \right)={\int }_{-\infty }^{\infty }f\left(t\right)e^{-j\omega t}dt\\ f\left(t\right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }F\left(\omega \right)e^{j\omega t}d\omega \end{array}$$

$$拉普拉斯变换:\{\begin{array}{l}F\left(s\right)={\int }_{-\infty }^{\infty }f\left(t\right)e^{-st}dt\\ f\left(t\right)=\frac{1}{2\pi j}{\int }_{\sigma -j\infty }^{\sigma +j\infty }F\left(s\right)e^{st}ds\end{array}$$

$$Z变换:\{\begin{array}{l}F\left(z\right)=\sum f\left[n\right]z^{-n}\\ f\left[n\right]=\frac{1}{2\pi j}\oint F\left(z\right)z^{n-1}dz\end{array}$$

  可以利用如下映射关系，实现传递函数$F(s)$、$F(\omega )$与$F(z)$之间的等价转换，详见《信号与系统》的“离散系统的Z域分析”章节：

$$s=j\omega =\frac{1}{T}\ln z，或，z=e^{sT}=e^{j\omega T}$$

  稳定性分析： 连续型传递函数稳定的充要条件为，所有极点均位于$s$域左半平面；离散型传递函数稳定的充要条件为，所有极点均位于$z$域原点为圆心的单位圆内。对于连续传函离散化问题，$s$域左半平面所有极点，经$z=e^{sT}$映射后，均位于$z$域原点为圆心的单位圆内。因此，严格离散化前后，连续系统稳定，则离散系统稳定。

离散化方法1，前向欧拉，前向差分

  对于一个连续函数$y(t)$在$kT$时刻的微分可以做以下近似：
$$y^{\prime }(k)\approx \frac{y(k+1)-y(k)}{T}$$

  对上式的左右两端分别进行拉氏变换和Z变换，可以得到如下映射关系：
$$sY(s)\leftrightarrow \frac{z-1}{T}Y(z)$$

  因此，可以令$s=\frac{z-1}{T}$实现S域传递函数到Z域传递函数的转换。值得注意的是，这种传递函数离散化是近似的，$s=\frac{1}{T}\ln z$才是严格的，但并不方便数字实现。实际上，$s=\frac{z-1}{T}$正是$s=\frac{1}{T}\ln z$，在$s=0$即$z=1$时的近似表示。

  稳定性分析： 假设连续传递函数是稳定的，即传递函数极点位于$s$域左半平面，分析离散化后极点在$z$域平面位置，
$$s=\frac{z-1}{T}\Rightarrow z=sT+1$$
  前向欧拉离散化方法将$s$域左半平面映射到$z$域$z<1$区域，显然包含不稳定区域，稳定极点可能变为不稳定极点。因此，稳定的连续系统变换所得离散系统可能不稳定，这种方法很少用。

离散化方法2，后向欧拉，后向差分

  同理还有后向差分，将$y(t)$在$kT$时刻的微分做以下近似：
$$y^{\prime }(k)\approx \frac{y(k)-y(k-1)}{T}$$
  对上式的左右两端分别进行拉氏变换和Z变换，可以得到如下映射关系：
$$sY(s)\leftrightarrow \frac{1-z^{-1}}{T}Y(z)$$
  因此，可以令$s=\frac{1-z^{-1}}{T}$实现S域传递函数到Z域传递函数的转换。实际上，$s=\frac{1-z^{-1}}{T}$也是$s=\frac{1}{T}\ln z$，在$s=0$即$z=1$时的近似表示。

  稳定性分析： 假设连续传递函数是稳定的，即连续传函极点位于$s$域左半平面，分析离散化后极点在$z$域平面位置，
$$s=\frac{1-z^{-1}}{T}\Rightarrow z=\frac{1}{1-sT}$$
  后向欧拉离散化方法将$s$域左半平面映射到$z$域以$(0.5,0)$为圆心，以$0.5$为半径的圆内（具体分析见附录）。也就是说，原来的连续系统稳定，则离散系统稳定。

离散化方法3，双线性变换，Tustin

  由于$z=e^{sT}=\frac{e^{sT/2}}{e^{-sT/2}}\approx \frac{1+sT/2}{1-sT/2}$，分子分母进行一阶泰勒线性展开，故得名双线性变换。进一步可以得到：
$$s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$$

  稳定性分析： 假设连续传递函数是稳定的，分析离散化后极点在$z$域平面位置，
$$s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}\Rightarrow z=\frac{1+sT/2}{1-sT/2}$$
  双线性变换离散化方法将$s$域左半平面映射到$z$域以$(0,0)$为圆心，以$1$为半径的单位圆内（具体分析见附录），这正是我们很希望看到的结果。因此，原来的连续传函稳定，则离散传函稳定。

离散化方法4，预矫正双线性变换，Tustin with Pre-Warping

问题描述

  上述方法均可实现连续系统的离散化，但因为离散化过程中采用了近似化手段，不可避免地引入了误差，使得离散化后的频域特性并不完全等价于连续系统，以下以双线性变换为例，进行说明。设计电网系统的五次谐波消抑制传递函数：
$$C_a(s)=\frac{K_{i}s}{s^2+{(2\pi \ast 250)}^2}$$
  取$K_i=1$，采用双线性变换$s=\frac{2}{T}\frac{z-1}{z+1}$进行离散化，其中T取0.001，得到离散化后的传递函数：
$$C_d(z)=C_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}=\frac{K_i(\frac{2}{T}\frac{z-1}{z+1})}{{\left(\frac{2}{T}\frac{z-1}{z+1}\right)}^2+{(2\pi \ast 250)}^2}=\frac{2000z^2-2000}{\left[{2000}^2+{(2\pi \ast 250)}^2\right]z^2+2\left[{(2\pi \ast 250)}^2-{2000}^2\right]z+\left[{(2\pi \ast 250)}^2+{2000}^2\right]}$$

  利用如下代码，对比连续传递函数和离散传递函数的波特图：

a=(2*pi*250)^2;
b=2000^2;
Ca = tf([1,0],[1 0 a]);
Cd = tf([2000 0 -2000],[b+a 2*a-2*b b+a],0.001);
bode(Ca, Cd)
legend('Ca','Cd')

  当然，也可以利用matlab实现双线性变换，两者效果是一致的（更多函数介绍请参考matlab函数用法：传递函数tf，波特图bode，离散化c2d）：

Ca = tf([1,0],[1 0 (2*pi*250)^2]);
Cd = c2d(Ca,0.001,'tustin');
bode(Ca, Cd);
legend('Ca','Cd')

  连续传递函数和离散传递函数的波特图对比结果如下所示。可以看出，近似的离散化过程使得离散传递函数特性发生了改变，离散化后传递函数的特性与所设计连续传递函数的特性不一致。

解决方案

  对离散化所引入的误差进行分析，$C_d(z)$代表数字控制器，$C_a(s)$代表模拟控制器，两者之间的转换关系如下：
$$C_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}=C_d(z)$$

  下面探究一下$C_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}$所对应的真实连续系统，注：$\sin \alpha =-j0.5\left(e^{j\alpha }-e^{-j\alpha }\right),\cos \alpha =0.5\left(e^{j\alpha }+e^{-j\alpha }\right)$:
$$\{\begin{array}{rl}{C}_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}& =C_a(\frac{2}{T}\frac{e^{sT}-1}{e^{sT}+1})\\ & =C_a\left(\frac{2}{T}\frac{e^{j\omega T}-1}{e^{j\omega T}+1}\right)\\ & =C_a\left[\frac{2}{T}\frac{e^{j\frac{\omega T}{2}}\left(e^{j\frac{\omega T}{2}}-e^{-j\frac{\omega T}{2}}\right)}{e^{j\frac{\omega T}{2}}\left(e^{j\frac{\omega T}{2}}+e^{-j\frac{\omega T}{2}}\right)}\right]\\ & =C_a\left[j\ast \frac{2}{T}\tan \left(\frac{\omega T}{2}\right)\right]\end{array}$$

  $C_d(z)$并不对应$C_a(j\omega )$，而对应$C_a\left[j\frac{2}{T}\tan \left(\frac{\omega }{2}T\right)\right]$，离散系统的$\omega$会映射到连续系统的$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)$，这就是频率畸变现象。也就是说，分别以$\omega$为横轴画出$C_d(z)$的波特图，和以$\frac{2}{T}\tan \left(\frac{\omega }{2}T\right)$为横轴画出$C_a(s)$的波特图，两个波形图应该是完全对应的。当$\frac{\omega T}{2}$足够小时，$j\frac{2}{T}\tan \left(\frac{\omega }{2}T\right)$是近似等于$j\omega$的。但对于高频谐波，或者特定频率的分析，就不得不考虑这部分影响了。

  回顾上述所设计$C_a(s)$时，对于$2\pi \ast 250$ rad/s的角速度变化，传递函数取得最大值。由于离散系统的$\omega$会映射到连续系统的$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)$，利用$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)=2\pi \ast 250$求解离散系统的原始$\omega$，解得$\omega =1331.5$ rad/s，这与上述波特图结果是一致的。

  解决方案：离散系统的${\omega }_d$与连续系统的${\omega }_a=\frac{2}{T}\tan \left(\frac{{\omega }_{d}T}{2}\right)$之间存在映射关系，根据目标${\omega }_d$计算${\omega }_a$，基于${\omega }_a$完成连续传递函数的设计，采用双线性变换得到离散传递函数。正确的设计思路为：（1）首先，明确所需要的离散传递函数的特性，比如，目标谐振频率${\omega }_d=2\pi \ast 250$ rad/s；（2）计算上述离散频率所映射的连续频率，${\omega }_a=\frac{2}{T}\tan \left(\frac{{\omega }_{d}T}{2}\right)=2000$ rad/s = 318 Hz；(3) 基于映射的连续角频率${\omega }_a$设计连续传递函数，$C_a(s)=\frac{K_{i}s}{s^2+2000^2}$，并完成离散化，$C_d(z)=C_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}=\frac{K_i(\frac{2}{T}\frac{z-1}{z+1})}{{\left(\frac{2}{T}\frac{z-1}{z+1}\right)}^2+2000^2}$。

  基于如下代码重绘波特图。可以看出，尽管连续传递函数在2000 rad/s (318Hz)处实现最大增益，离散传递函数在目标频率点1570 rad/s （250Hz）处取得最大增益，与预期一致。

Ca = tf([1,0],[1 0 (2000)^2]);
Cd = c2d(Ca,0.001,'tustin');
bode(Ca, Cd);
legend('Ca','Cd')

总结

  对离散化的简单场景，优先考虑后向差分离散化方法；对于频率有较高准确度的场景，如谐振控制器，建议使用预矫正双线性变换方法。

附录

附录1：后向欧拉相关证明

  对于$s$域的左半平面，其边界为$s$平面虚轴。证明：该边界在$z$域的映射为，以$(0.5,0)$为圆心，以$0.5$为半径的圆，如下所示:
$$\{\begin{array}{rl}z& =\frac{1}{1-sT}\\ & =\frac{1}{1-j\omega }\\ & =\frac{1+j\omega }{\left(1-j\omega \right)\left(1+j\omega \right)}\\ & =\frac{1}{1+{\omega }^2}+j\frac{\omega }{1+{\omega }^2}\end{array}\to \{\begin{array}{rl}{\Vert z-0.5\Vert }^2& ={\left(\frac{1}{1+{\omega }^2}-0.5\right)}^2+{\left(\frac{\omega }{1+{\omega }^2}\right)}^2\\ & =\frac{1}{4}\left[{\left(\frac{2}{1+{\omega }^2}-1\right)}^2+{\left(\frac{2\omega }{1+{\omega }^2}\right)}^2\right]\\ & =\frac{1}{4}\left[{\left(\frac{1-{\omega }^2}{1+{\omega }^2}\right)}^2+{\left(\frac{2\omega }{1+{\omega }^2}\right)}^2\right]\\ & =\frac{1}{4}，证毕\end{array}$$

  进一步证明，$s$域内任何左半平面的点，映射到上述圆内：
$$\{\begin{array}{rl}z& =\frac{1}{1-sT},Re(s)<0\\ & =\frac{1}{a-bj},a>1\\ & =\frac{a+bj}{\left(a-bj\right)\left(a+bj\right)}\\ & =\frac{a}{a^2+b^2}+j\frac{b}{a^2+b^2}\end{array}\to \{\begin{array}{rl}{\Vert z-0.5\Vert }^2& ={\left(\frac{a}{a^2+b^2}-0.5\right)}^2+{\left(\frac{b}{a^2+b^2}\right)}^2\\ & =\frac{1}{4}\left[{\left(\frac{2a}{a^2+b^2}-1\right)}^2+{\left(\frac{2b}{a^2+b^2}\right)}^2\right]\\ & <\frac{1}{4}\left[{\left(\frac{2a}{a^2+b^2}-1\right)}^2+{\left(\frac{2ab}{a^2+b^2}\right)}^2\right]\\ & <\frac{1}{4}\left[{\left(\frac{a^2-b^2}{a^2+b^2}\right)}^2+{\left(\frac{2ab}{a^2+b^2}\right)}^2\right]\\ & =\frac{1}{4}，证毕\end{array}，其中\{\begin{array}{rl}& \frac{2a}{a^2+b^2}-1<\frac{a^2-b^2}{a^2+b^2}\\ & \leftarrow 2a-a^2-b^2<a^2-b^2\\ & \leftarrow 2a<2a^2\\ & \leftarrow 1<a\end{array}$$

  Matlab绘图验证：

  Matlab代码如下：

clc; clear; close all;

% 定义参数
T = 1;

% 定义s=0（虚部从-无穷到正无穷）
omega = linspace(-1000, 1000, 20001);  % 使用大的范围代替无穷
s1 = 0.0 + 1j * omega;
s2 = -0.1 + 1j * omega;
s3 = -1.0 + 1j * omega;
s4 = -10.0 + 1j * omega;

% 映射到z平面
z1 = 1 ./ (1 - s1 * T);
z2 = 1 ./ (1 - s2 * T);
z3 = 1 ./ (1 - s3 * T);
z4 = 1 ./ (1 - s4 * T);

% 绘制z平面的值
figure;
hold on;

% 绘制s=0的映射
plot(real(z1), imag(z1), 'b', 'DisplayName', 's = 0 (Im: -\infty to +\infty)');

% 绘制s=-0.1的映射
plot(real(z2), imag(z2), 'r', 'DisplayName', 's = -0.1 + j\omega');

% 绘制s=-1.0的映射
plot(real(z3), imag(z3), 'g', 'DisplayName', 's = -1.0 + j\omega');

% 绘制s=-10.0的映射
plot(real(z4), imag(z4), 'm', 'DisplayName', 's = -10.0 + j\omega');

% 设置图形属性
xlabel('Real Part');
ylabel('Imaginary Part');
title('Mapping of s-plane lines to z-plane');
legend show;
grid on;
axis equal;
hold off;

附录2：双线性变换相关证明

  对于$s$域的左半平面，其边界为$s$平面虚轴。证明：该边界在$z$域的映射为，以$(0,0)$为圆心，以$1$为半径的圆，如下所示:
$$\{\begin{array}{rl}z& =\frac{1+sT/2}{1-sT/2}\\ & =\frac{1+j\omega }{1-j\omega }\\ & =\frac{1-{\omega }^2}{1+{\omega }^2}+j\frac{2\omega }{1+{\omega }^2}\end{array}\to \{\begin{array}{rl}{\Vert z\Vert }^2& ={\left(\frac{1-{\omega }^2}{1+{\omega }^2}\right)}^2+{\left(\frac{2\omega }{1+{\omega }^2}\right)}^2\\ & =1\end{array}$$
  对于边界内的点，类比后向欧拉方法，可证明其在$z$平面的映射在单位圆内。给出Matlab绘图结果：

  Matlab代码如下：

clc; clear; close all;
% 定义参数
T = 1;

% 定义s=0（虚部从-无穷到正无穷）
omega = linspace(-1000, 1000, 20001);  % 使用大的范围代替无穷
s1 = 0.0 + 1j * omega;
s2 = -0.1 + 1j * omega;
s3 = -1.0 + 1j * omega;
s4 = -10.0 + 1j * omega;

% 映射到z平面
z1 = (1 + s1 * T) ./ (1 - s1 * T);
z2 = (1 + s2 * T) ./ (1 - s2 * T);
z3 = (1 + s3 * T) ./ (1 - s3 * T);
z4 = (1 + s4 * T) ./ (1 - s4 * T);

% 绘制z平面的值
figure;
hold on;

% 绘制s=0的映射
plot(real(z1), imag(z1), 'b', 'DisplayName', 's = 0 (Im: -\infty to +\infty)');

% 绘制s=-0.1的映射
plot(real(z2), imag(z2), 'r', 'DisplayName', 's = -0.1 + j\omega');

% 绘制s=-1.0的映射
plot(real(z3), imag(z3), 'g', 'DisplayName', 's = -1.0 + j\omega');

% 绘制s=-10.0的映射
plot(real(z4), imag(z4), 'm', 'DisplayName', 's = -10.0 + j\omega');

% 设置图形属性
xlabel('Real Part');
ylabel('Imaginary Part');
title('Mapping of s-plane lines to z-plane');
legend show;
grid on;
axis equal;
hold off;