微分方程的离散化

当下的控制器通常为数字控制器而非模拟控制器，其控制输入是离散化的，因此，我们需要考虑将现实的连续系统进行离散化。

理论基础

实践中，对连续信号进行离散化的方法是 采样 和 保持 。

采样

考虑连续信号 \(x(t)\) ，我们可以将针对连续信号的采样电路近似为一个采样开关，开关闭合时，采样结果等于原始信号；开关打开时，采样结果为 0 .

若开关闭合时间无限趋近于 0 ，则为理想采样。称采样开关闭合的周期为采样周期 \(T_s\) ，其倒数为采样频率 \(f\)，再乘以 \(2\pi\) 就是采样角频率 \(\omega_s\) 。根据采样定理，采样频率至少需要大于原始信号最高频率的两倍，才能从采样结果中恢复原始信号。这是理论下限，实践下限一般为 5-10 倍。

\[x^{*}(t) = x(t)\sum_{k=-\infty}^{+\infty}\delta(t-kT_s)\]

保持

由于采样结果 \(x^{*}(t)\) 几乎处处为 0 ，我们不能直接将其用于后续控制系统，需要增加一个保持器。保持器的作用是将瞬时的采样信号按某种方式保持到下一采样时刻。最常用的是零阶保持器（ZOH，zero-order hold），即将 \(kT_s\) 时刻的采样值 \(x(kT_s)\) 一直保持到下一采样时刻 \((k+1)T_s\) .

\[x_h(t) = \sum_{k = 0}^{\infty}x(kT_s)\left[u(t-kT_s)-u(t-kT_s-T_s)\right]\]

保持器本质上是一种低通滤波器，并且具有一定的滞后性。

ZOH 的频率特性

我们假设采样信号是一个理想单位脉冲 \(x^{*}(t) = \delta(t)\) ，其拉氏变换结果为 1 . 经过 ZOH 后，输出信号为 \(x_{h}(t) = u(t) - u(t-T_s)\) ，其拉氏变换结果为 \(\dfrac{1}{s}\left(1-e^{-T_ss}\right)\) ，由此可写出 ZOH 的传递函数

\[H_o(s) = \frac{X_{h}(s)}{X^{*}(s)} = \frac{1-e^{-T_s s}}{s}\]

令 \(s = j\omega\) ，可得 ZOH 的频率特性

\[\begin{aligned} H_o(j\omega) &= \frac{e^{-j\omega T_s/2}\left(e^{j\omega T_s/2}-e^{-j\omega T_s/2}\right)}{j\omega} = \frac{2j\sin(\omega T_s/2)e^{-j\omega T_s/2}}{j\omega}\\ &= \frac{2\sin(\omega T_s/2)}{\omega}e^{-j\omega T_s/2}\\ |H_o(j\omega)| &= \frac{2|\sin(\omega T_s/2)|}{\omega} = T_s\frac{|\sin(\omega T_s/2)|}{\omega T_s/2}\\ \angle H_o(j\omega) &= \begin{cases}-\omega T_s/2\ ,& \sin(\omega T_s/2) > 0\\ \pi+2k\pi - \omega T_s/2\ ,& \sin(\omega T_s/2) < 0 \end{cases} \end{aligned}\]

可见，频率越高，幅度越小，滞后越大。

无误差离散化

Warning

无误差离散化不代表离散系统与原始系统性质完全相同！离散系统是原始系统采样+保持（ZOH）的结果。

SISO 系统

考虑一阶惯性环节系统

\[T\frac{\mathrm{d}y(t)}{\mathrm{d}t}+y(t) = Kx(t) \]

在输入端引入虚拟采样开关和零阶保持器，用以表示实践中对输入信号的采样和保持。则实际系统的输入信号为

\[x_o(t) = x(kT_s)\ , kT_s\le t < (k+1)T_s\ ,k=0,1,2\cdots\]

仅考虑 \(kT_s\le t < (k+1)T_s\) 这一时间范围，输入信号恒为 \(x(kT_s)\) ，求解输出信号，若假设 \(t' = t-kT_s\) 和 \(K_m = Kx(kT_s)\) ，则待求解的微分方程为：

\[T\dot{y} + y = K_m\]

先令 \(K_m\) 为 0 ，求通解（零输入分量）。该形式的微分方程解的形式为 \(Ce^{at}\) ，代入可求得 \(a = -\dfrac{1}{T}\) ，因此零输入分量为 \(y_{zi}(t) = Ce^{-\frac{t-kT_s}{T}}\) . 再求特殊解（零状态分量）。观察可知 \(y_{zs}(t) = K_m = Kx(kT_s)\) 即为一特殊解。因此，解为

\[y(t) = y_{zi}(t)+y_{zs}(t) = Ce^{-\frac{t-kT_s}{T}}+Kx(kT_s)\]

考虑 \(t=kT_s\) 时，输出信号为 \(y(kT_s)\) ，代入可解出 \(C\) .

\[y(t) = y_{zi}(t)+y_{zs}(t) = y(kT_s)e^{-\frac{t-kT_s}{T}}+Kx(kT_s)\left(1-e^{-\frac{t-kT_s}{T}}\right)\]

输出信号是连续的，我们可以用上式求 \(t=(k+1)T_s\) 时的输出信号

\[y[(k+1)T_s] = y(kT_s)e^{-\frac{T_s}{T}}+Kx(kT_s)\left(1-e^{-\frac{T_s}{T}}\right)\]

上式已经是一个离散的差分方程，等价于

\[y[k+1] = e^{-\frac{T_s}{T}}y[k] + K\left(1-e^{-\frac{T_s}{T}}\right)x[k]\]

如果 \(K = 1\) （当 \(x(t)\) 稳定时， \(y(t)\) 将收敛到 \(x(t)\) ），上式具有很漂亮的形式：

\[y[k+1] = \alpha y[k] + (1-\alpha)x[k]\]

MIMO 系统

考虑状态空间方程

\[\dot{\mathbf{x}} = A\mathbf{x} + B\mathbf{u}\]

参见 线性定常系统状态空间方程的解 ，其解为

\[\mathbf{x}(t) = e^{A(t-t_0)}\mathbf{x}_0 + \int_{t_0}^{t}e^{A(t-\tau)}B\mathbf{u}\mathrm{d}\tau\]

同样地，仅考虑 \(kT_s\le t < (k+1)T_s\) 这一时间范围，输入信号恒为 \(\mathbf{u}(kT_s)\) ，求解状态向量：

\[\begin{aligned} \mathbf{x}((k+1)T_s) &= e^{AT_s}\mathbf{x}(kT_s) + \int_{kT_s}^{(k+1)T_s}e^{A((k+1)T_s-\tau)}B\mathbf{u}(kT_s)\mathrm{d}\tau\\ &= e^{AT_s}\mathbf{x}(kT_s) + \left[\int_{T_s}^{0} e^{As}\mathrm{d}(-s)\right]B\mathbf{u}(kT_s)\\ &= e^{AT_s}\mathbf{x}(kT_s) + \left(\int_{0}^{T_s} e^{As}\mathrm{d}s\right)B\mathbf{u}(kT_s)\\ \Rightarrow \mathbf{x}[k+1] &= e^{AT_s}\mathbf{x}[k] + \left(\int_{0}^{T_s} e^{As}\mathrm{d}s\right)B \mathbf{u}[k] \end{aligned}\]

实例

考虑系统状态空间方程为

\[\dot{\mathbf{x}} = \begin{bmatrix}0 & 1\\ 0 & 0\end{bmatrix}\mathbf{x}+\begin{bmatrix}0\\1\end{bmatrix}\mathbf{u}\]

易得

\[\begin{aligned} e^{AT_s} &= I+AT_s = \begin{bmatrix}1 & T_s\\ 0 & 1\end{bmatrix}\\ \left(\int_{0}^{T_s} e^{As}\mathrm{d}s\right)B &= \begin{bmatrix}T_s & \frac{1}{2}T_s^2\\ 0 & T_s\end{bmatrix}\begin{bmatrix}0\\1\end{bmatrix}=\begin{bmatrix}\frac{1}{2}T_s^2\\T_s\end{bmatrix} \end{aligned}\]

所以离散化后的结果为

\[\mathbf{x}[k+1] = \begin{bmatrix}1 & T_s\\ 0 & 1\end{bmatrix}\mathbf{x}[k]+\begin{bmatrix}\frac{1}{2}T_s^2\\T_s\end{bmatrix}\mathbf{u}[k]\]

近似差分离散化

针对计算精度要求不高的场合，可以考虑采用近似方法得到系统的差分方程。原理是泰勒展开一阶近似表示一阶导数。

考虑系统微分方程为

\[\dot{y} = f(x(t), y(t))\]

前向差分/前向欧拉积分法/显式欧拉法：

\[\begin{aligned}y(t+T_s) &= y(t) + y'(t)T_s + y''(t)\frac{T_s^2}{2} + y^{(3)}(t)\frac{T_s^{3}}{6}+\cdots\\ y'(t) &\approx \frac{y(t+T_s)-y(t)}{T_s}\\ \Rightarrow y[k+1] &= y[k] + f(x[k], y[k])T_s \end{aligned}\]

这是一个显式方程，已知 \(y[k]\) 的情况下，可以直接计算出 \(y[k+1]\) .

后向差分/后向欧拉积分法/隐式欧拉法：

\[\begin{aligned}y(t-T_s) &= y(t) - y'(t)T_s+y''(t)\frac{T_s^2}{2}-y^{(3)}(t)\frac{T_s^3}{6}+\cdots\\ y'(t) &\approx \frac{y(t)-y(t-T_s)}{T_s}\\ \Rightarrow y[k-1] &= y[k]-f(x[k],y[k])T_s\\ y[k] &= y[k-1] + f(x[k],y[k])T_s \end{aligned}\]

这是一个隐式方程，已知 \(y[k-1]\) 的情况下，需要解该方程才能计算出 \(y[k]\) .

中心差分：

\[\begin{aligned}y(t+T_s)-y(t-T_s) &= 2 y'(t)T_s+y^{(3)}(t)\frac{T_s^3}{3}+\cdots\\ y'(t) &\approx \frac{y(t+T_s)-y(t-T_s)}{2T_s}\\ \Rightarrow y[k+1]-y[k-1] &= 2f(y[k],x[k])T_s\\ y[k+1] &= 2f(y[k],x[k])T_s+y[k-1] \end{aligned}\]

中心差分应该属于显式欧拉法。

误差分析

显然，前向/后向差分的误差均为 \(\dfrac{O\left(T_s^2\right)}{T_s} = O(T_s)\) 级别的一阶误差，而中心差分的误差为 \(\dfrac{O\left(T_s^3\right)}{T_s} = O\left(T_s^2\right)\) 级别的二阶误差。因此，中心差分的误差最小，前向/后向差分误差级别相同。

频率特性分析

考虑一个频率为 \(\omega\) 的信号 \(\sin(\omega t)\) ，求导后变为 \(\omega\sin\left(\omega t+\frac{\pi}{2}\right)\) . 这说明，求导算子的幅频特性是 \(\omega\) ，相频特性是相位超前 \(\frac{\pi}{2}\) .

下面为了方便分析，我们考察频率为 \(\omega\) 的信号 \(x(t) = e^{j\omega t}\) ，其求导结果为 \(\omega e^{j(\omega t + \frac{\pi}{2})}\) ，与上一段的分析结果相同。

在 \(t=kT_s\) 处进行前向差分，有 \(x[k] = e^{j\omega kT_s}\)

\[\begin{aligned} \frac{x[k+1]-x[k]}{T_s} &= \frac{e^{j\omega(k+1)T_s}-e^{j\omega kT_s}}{T_s}\\ &= \frac{e^{j\omega T_s}-1}{T_s}e^{j\omega k T_s} \end{aligned}\]

可见，前向差分算子的频率响应是 \(\dfrac{e^{j\omega T_s}-1}{T_s} = \dfrac{2\sin\frac{\omega T_s}{2}}{T_s}e^{j\left(\frac{\omega T_s}{2}+\frac{\pi}{2}\right)}\) . 其幅频特性是 \(\left|\dfrac{2\sin\frac{\omega T_s}{2}}{T_s}\right|\) ，相频特性是相位超前 \(\dfrac{\omega T_s}{2}+\dfrac{\pi}{2}\) . 相对于求导算子，相位超前 \(\dfrac{\omega T_s}{2}\) .

同样的，在 \(t=kT_s\) 处进行后向差分

\[\begin{aligned} \frac{x[k]-x[k-1]}{T_s} &= \frac{e^{j\omega kT_s}-e^{j\omega (k-1)T_s}}{T_s}\\ &= \frac{1-e^{-j\omega T_s}}{T_s}e^{j\omega k T_s} \end{aligned}\]

可见，后向差分算子的频率响应是 \(\dfrac{1-e^{-j\omega T_s}}{T_s} = \dfrac{2\sin\frac{\omega T_s}{2}}{T_s}e^{j\left(\frac{\pi}{2}-\frac{\omega T_s}{2}\right)}\) . 其幅频特性是 \(\left|\dfrac{2\sin\frac{\omega T_s}{2}}{T_s}\right|\) ，相频特性是相位超前 \(\dfrac{\pi}{2}-\dfrac{\omega T_s}{2}\) . 相对于求导算子，相位滞后 \(\dfrac{\omega T_s}{2}\) .

稳定性分析

考虑 \(y' = \lambda y, \lambda \in \mathbb{C}\) ，其解为 \(y = y(0)e^{\lambda t}\) ，当且仅当 \(y(0)\) 不为 0 且 \(\mathbf{Re}(\lambda) < 0\) 时稳定。

若采用前向差分离散化，则有 \(y[k+1]-y[k] = \lambda T_s y[k]\) ，即

\[y[k+1] = (1+\lambda T_s)y[k]\]

当且仅当 \(|1+\lambda T_s| < 1\) 时，系统才能稳定。这代表 \(\lambda T_s\) 必须位于复平面上，以 \((-1,0)\) 为圆心，以 1 为半径的一个圆盘内。

1. 原系统稳定，前向差分后不一定稳定。
2. 若 \(\lambda = j\omega\) （ \(y\) 为纯周期信号），前向差分系统会发散。

若采用后向差分离散化，则有 \(y[k]-y[k-1] = \lambda T_s y[k]\) ，即

\[y[k] = \frac{1}{1-\lambda T_s}y[k-1]\]

当且仅当 \(\left|\dfrac{1}{1-\lambda T_s}\right| < 1\) 时，系统才能稳定。这代表 \(\lambda T_s\) 必须位于复平面上，以 \((1,0)\) 为圆心，以 1 为半径的一个圆盘外。

1. 原系统稳定，后向差分后必然稳定。
2. 若 \(\lambda = j\omega\) （ \(y\) 为纯周期信号），后向差分系统会收敛。

Z 变换分析

对于一阶导数 \(y'\) ，其复频域表示为 \(sY(s)\)

前向差分：

\[\dfrac{y[k+1]-y[k]}{T_s} \overset{Z 变换}\rightarrow \dfrac{z-1}{T_s}Y(z) \Rightarrow s = \frac{z-1}{T_s}\]

后向差分：

\[\dfrac{y[k]-y[k-1]}{T_s} \overset{Z 变换}\rightarrow \dfrac{1-z^{-1}}{T_s}Y(z) \Rightarrow s = \frac{1-z^{-1}}{T_s}\]

双线性变换（Tustin 变换）：

\[s = \frac{2}{T_s}\frac{1-z^{-1}}{1+z^{-1}}\]