1.符号说明与结构框图

y(k)——系统响应输出的离散值
u(k)——数字PID控制输出的离散值
r(k)——期望输出的离散值（事先已知），在本例中为常数（即阶跃输入）
e(k)——e(k)=r(k)-y(k)，为期望值-实际值，是单位负反馈的误差比较信号
D(z)——数字控制器的脉冲传递函数
G(s)——被控对象的传递函数，G(z)——广义被控对象的脉冲传递函数
$Φ(z)=C(z)/R(z)\Phi(z)=C(z)/R(z)$ ——系统的闭环传递函数， $Φe(z)=E(z)/C(z)\Phi_e(z)=E(z)/C(z)$ ——系统的闭环误差传递函数
系统的结构框图如下：

2.最小拍控制系统构造原则

最少拍系统是指系统对某些典型的输入（阶跃、速度、加速度输入）具有最快的响应特性。具体来说，最少拍系统应满足对典型输入在有限个采样周期内结束过渡过程且稳态误差为零。研究误差序列 $e (k)$ 的特性之前，不妨通过在Z域研究 $Φe(z)\Phi_e(z)$ 的表达式来得出使得 $e (k)$ 最快收敛至0的条件。
最小拍控制系统的设计原理是使得系统的闭环误差传递函数 $Φe(z)=a1z−1+a2z−2+...+anz−n\Phi_e(z)=a_1z^{-1}+a_2z^{-2}+...+a_nz^{-n}$ 具有最简单的表达式。因为 $z^{-k}$ 的反Z变换为 $δ(t−k)\delta(t-k)$ ，因此 $e(k)=a1δ(t−T)+a2δ(t−2T)+...+anδ(t−nT)(an≠0)e(k)=a_1\delta(t-T)+a_2\delta(t-2T)+...+a_n\delta(t-nT)(a_n≠0)$ ，即 $e(1)=a_1,e(2)=a_2,...,e(n)=a_n$ n的值是有限的说明系统能在有限拍内达到最小，n的值越小， $Φe(z)\Phi_e(z)$ 的表达式越简单，能达到稳态值所需的拍数（周期数）越少。而我们的目的就是设计控制器D(z)使其达到稳态误差e(k)（甚至是e(t)）的要求。
设计步骤大体如下：

根据被控对象的数学模型求出广义被控对象的脉冲传递函数G(z)
根据输入信号的类型，确定模型 $Φe(z)=1−Φ(z)\Phi_e(z)=1-\Phi(z)$ ,其对应关系如下表：

典型输入	Z域表达式	$Φe(z)的形式\Phi_e(z)的形式$
$u (t)$	$11−z−1\frac{1}{1-z^{-1}}$	$1-z^{-1})F(z)$
$t u (t)$	$Tz−1(1−z−1)2\frac{Tz^{-1}}{{(1-z^{-1})}^2}$	${(1-z^{-1})}^2F(z)$
$12t2u(t)\frac{1}{2}t^2u(t)$	$T2z−1(1+z−1)2(1−z−1)3\frac{T^2z^{-1}(1+z^{-1})}{2{(1-z^{-1})}^3}$	${(1-z^{-1})}^3F(z)$

当被控对象不含有不稳定零极点，不含有纯滞后环节 $z^{-1}$ 时，F(z)=1，此时的被控对象也称为简单对象。
满足以上条件的任意一个，被控对象称为复杂对象，当G(s)中含有不稳定零点或纯滞后环节时，不能使用D(z)或者 $Φe(z)\Phi_e(z)$ 中增加因式消除，只能保留至 $Φ(z)\Phi(z)$ ，当G(s)中含有不稳定极点时，在 $Φe(z)\Phi_e(z)$ 中增加对应的零点抵消 $Φ(z)\Phi(z)$ 中的不稳定极点。
3. 求控制器的脉冲传递函数D(z)，求取原理如下：
由于 $C (z) = E (z) D (z) G (z)$ ，所以 $C(z)R(z)=D(z)G(z)E(z)R(z)\frac{C(z)}{R(z)}=D(z)G(z)\frac{E(z)}{R(z)}$ ，即 $Φe(z)=D(z)G(z)Φe(z)\Phi_e(z)=D(z)G(z)\Phi_e(z)$ 在单位负反馈系统中 $Φe(z)=1−Φ(z)\Phi_e(z)=1-\Phi(z)$ ，因此 $D(z)=Φ(z)G(z)(1−Φ(z))D(z)=\frac{\Phi(z)}{G(z)(1-\Phi(z))}$
4. 根据D(z)生成控制算法（或者脉冲传递函数）求出输出序列，画出系统的响应曲线。

2.1数字控制器D(z)的构造

我们已经知道 $D(z)=Φ(z)G(z)(1−Φ(z))D(z)=\frac{\Phi(z)}{G(z)(1-\Phi(z))}$ ，因此欲求出D(z)只需知道 $Φ(z)\Phi(z)$ 或 $Φe(z)\Phi_e(z)$ 即可。
$Φ(z)=(a0+a1z−1+a2z−2+...+anz−n)(1−z1z−1)(1−z2z−1)...(1−zkz−1)z−m\Phi(z)=(a_0+a_1z^{-1}+a_2z^{-2}+...+a_nz^{-n})(1-z_1z^{-1})(1-z_2z^{-1})...(1-z_kz^{-1})z^{-m}$
其中 $z_1,z_2,...,z_k$ 为 $Φ(z)\Phi(z)$ 保留的的广义被控对象G(z)不稳定零点， $z^{-m}$ 是闭环传递函数 $Φ(z)\Phi(z)$ 保留的，在G(z)里面出现的纯滞后环节 $z^{-m}$ ，而输入形式决定了前面多项式 $a_0+a_1z^{-1}+a_2z^{-2}+...+a_nz^{-n}$ 的阶数n，阶数n和 $Φe(z)=(1−z−1)pF(z)\Phi_e(z)={(1-z^{-1})}^pF(z)$ 中的p是一致的，即n=p，系数 $a_0、a_1、...、a_n$ 依赖于条件 $Φe(z)=(1−z−1)pF(z)\Phi_e(z)={(1-z^{-1})}^pF(z)$ 而定， $Φe(z)\Phi_e(z)$ 中含有p重零点z=1，因此它的0阶导到p-1阶导在z=1处的值都应该是0，即

${Φe(z)∣z=1=0dΦe(z)dz∣z=1=0d2Φe(z)dz2∣z=1=0...dnΦe(z)dzn∣z=1=0\begin{cases} \left. \Phi_e(z) \right| _{z=1}=0 \\ \left. \frac{{\rm d}\Phi_e(z)}{{\rm d}z} \right| _{z=1}=0 \\ \left. \frac{{\rm d^2}\Phi_e(z)}{{\rm d}z^2} \right| _{z=1}=0 \\ ...\\ \left. \frac{{\rm d^n}\Phi_e(z)}{{\rm d}z^n} \right| _{z=1}=0 \end{cases}$
根据上面的方程组可以解出 $Φe(z)\Phi_e(z)$ 中的未知参数 $a_0、a_1、...、a_n$ ，因此 $Φe(z)\Phi_e(z)$ 已经确定，D(z)即可求得。

3.简单控制对象的最小拍控制器设计

3.1阶跃输入

假设被控对象的传递函数为 $G(s)=9.40.017s+1G(s)=\frac{9.4}{0.017s+1}$ 求阶跃输入下的最小拍控制器D(z)的表达式。在Simulink的仿真中，我们设定了Z变换的采样频率是0.005s，（其他地方默认为-1，即服从系统的采样频率）
Simulink仿真图示意
先求取广义被控对象的传递函数
$=Z[\frac{1-e^{-Ts}}{s}\frac{9.4}{0.017s+1}] =(1-z^{-1})Z[(\frac{9.4}{s(0.017s+1)})] =\frac{2.394z^{-1}}{1-0.745z^{-1}}$
没有不稳定零极点和纯滞后环节，因此选取 $Φe(z)=(1−z−1),Φ(z)=z−1\Phi_e(z)=(1-z^{-1}),\Phi(z)=z^{-1}$
此时 $D(z)=Φ(z)G(z)(1−Φ(z))=0.4174z−0.7452z−1D(z)=\frac{\Phi(z)}{G(z)(1-\Phi(z))}=0.4174\frac{z-0.7452}{z-1}$
我们搭建好Simulink仿真电路如下：
Simulink的仿真图
按照我们的理论，系统的应该在控制算法作用的第一拍（0.005s处）时达到0稳态误差，仿真曲线如下图（给定值为1500）：
最小拍的阶跃响应
可以看到验证结果与我们理论推导相符。注意我们这里示波器取的采样周期为-1（inherited），说明是Matlab内置的采样时间，会比我们定的Ts=1s要短的多，非常接近于连续系统的仿真结果。可以看出响应曲线在上升段有很微小的纹波抖动，在进入稳态值后，纹波消失。产生纹波的原因在于 $Y (z)$

3.2速度输入

假设被控对象的传递函数为 $G(s)=2s(s+1)G(s)=\frac{2}{s(s+1)}$ ，求系统在采样时间1s，输入信号为单位速度输入的条件下的最小拍控制器D(z)。系统的结构框图如下图所示：
系统结构框图
下面我们借助Matlab的命令轻松地解决这个问题。首先按照第二部分最小拍控制系统的设计原则的步骤，我们应该求出带零阶保持器的广义被控对象的Z域脉冲传递函数 $G(z)=Z[1−e−Tss⋅G(s)]G(z)=Z[\frac{1-e^{-Ts}{s}\cdot G(s)}]$ ，在Matlab工作区输入代码并获得输出的脉冲传递函数：

>> syms s;%定义符号变量s
>> s=tf('s');%s定义为tf（transfer function传递函数）类型的结构体
>> Gs=2/(s*(s+1));
>> Ts=1;
>> Gz=c2d(Gs,Ts,'zoh')%带零阶保持器（ZOH）的离散化，采样时间为1sGz =0.7358 z + 0.5285----------------------z^2 - 1.368 z + 0.3679Sample time: 1 seconds
Discrete-time transfer function.>> zpk(Gz)%展示为零极点模式ans =0.73576 (z+0.7183)------------------(z-1) (z-0.3679)Sample time: 1 seconds
Discrete-time zero/pole/gain model.

可以看出脉冲传递函数G(z)中含有1个不稳定极点z=1，但是考虑到此时的 $Φe(z)=(1−z−1)2F(z)\Phi_e(z)={(1-z^{-1})}^2F(z)$ ，具有z=1这个零点，说明D(z)的表达式 $D(z)=Φ(z)G(z)Φe(z)D(z)=\frac{\Phi(z)}{G(z)\Phi_e(z)}$ 中可以将极点z=1约去，不必再增加额外的 $1-z^{-1})$ 因式，因此可以取F(z)=1即 $Φe(z)=(1−z−1)2\Phi_e(z)={(1-z^{-1})}^2$ ，此时 $Φ(z)=1−Φe(z)=1−(1−z−1)2=2z−1−z−2\Phi(z)=1-\Phi_e(z)=1-{(1-z^{-1})}^2=2z^{-1}-z^{-2}$ ，根据 $D(z)=Φ(z)G(z)Φe(z)D(z)=\frac{\Phi(z)}{G(z)\Phi_e(z)}$ 可以求出控制器表达式D(z)。

>> syms z;
>> z=tf('z');
>> Phiez=(1-z^(-1))^2;
>> Phiz=1-Phiez;%得到Φ(z)
>> Dz=Phiz/(Gz*Phiez);%得到数字控制器的脉冲传递函数
>> minreal(zpk(Dz))%得到最简零极点式ans =2.7183 (z-0.5) (z-0.3679)-------------------------(z+0.7183) (z-1)Sample time: 1 seconds
Discrete-time zero/pole/gain model.

我们按照此要去搭建好Simulink仿真图如下

D(z)的Sample time设置为1（因为我们的采样时间是Ts=1s），仿真时间设为10s（差不多可以清楚地看见过渡过程）
速度输入的仿真
黄色为指令信号（单位速度输入），蓝色为响应曲线。我们将示波器的图形（print to figure）打印到图形。利用数据游标工具找到t=1和t=2的点。如果我们的理论正确，那么e(0)=0,e(1)=1,e(2)=0.e(3)=e(4)=…=0;

t	数据游标（采样时刻的值）
1
2
3
4

由于r(k)=k,y(k)的采样值依次为0，2，3.002，3.997，e(k)=r(k)-y(k)基本满足e(0)=0,e(1)=1,e(2)=0.e(3)=e(4)=…=0。
从我们仿真的结果看，采样点处的值只需要两拍就达到了稳态，但是采样点之间有纹波（这是因为 $E(z)=a_1z^{-1}+a_2z^{-2}+...+a_nz^{-n}$ 并不是有限项，U(z)也并不是有限项，即u(k)最终并没有达到恒定，误差也并没有最终严格归0），这种控制属于有纹波的最小拍控制。还可以通过增加积分可牺牲控制拍数的手段，达到无纹波的最小拍控制，这是连续系统所做不到的。
希望本文对您有帮助，感谢大家对本站点的支持。