2.3 连续系统仿真的离散相似算法

从前面介绍的各种数值积分算法的公式可以看出，它们都是一些差分方程，我们仅仅计算了在各个离散时间点tk处系统的近似解yk。这相当于对连续系统进行了离散化处理，把原来的连续系统模型近似等价为一个离散系统模型。从本质上讲，连续系统的数字仿真就是要找出一个与该系统等价的离散模型。因此，数值积分算法也是离散化算法，只不过在推导过程中是从数值积分的角度出发，没有明确地提出“离散”这个概念。本节将从连续系统离散化的角度出发，建立连续系统模型的等价离散化模型，并用采样系统的理论和方法介绍另一种常用的仿真算法。这种算法使得连续系统在进行（虚拟的）离散化处理后仍保持与原系统“相似”，故称之为离散相似算法。

与数值积分算法相比，离散相似算法的每步计算量要小得多，稳定性也要好得多，因而允许采用较大的计算步长。然而，它通常只适合线性定常系统的仿真，具有一定的局限性。

离散相似算法的物理概念明确，所以特别容易被从事控制系统分析和研究的工程技术人员所掌握。

由于连续系统可以用状态空间模型来表示，也可以用传递函数来描述，因此，与连续系统等价的离散化模型可以通过两种途径来获得。一是对状态空间模型作离散化处理得到离散化状态空间模型，称为时域离散相似算法；另一种是对传递函数进行离散化得到脉冲传递函数，称为z域离散相似算法。考虑到应用范围与实现的便利性，本节仅介绍时域离散相似算法。

2.3.1 时域离散相似算法的基本概念

连续系统的离散化是在原连续系统的输入端、输出端分别加上虚拟的采样开关，使输入、输出信号离散化，此时系统模型即为离散化模型。为了使输入信号u（t）离散化后仍能保持原来的变化规律，在输入采样开关后设置信号保持器（亦称为信号重构器），复现原输入信号u（t），其结构如图2.21所示。事实上，各种保持器要完全不失真地复现输入信号是不可能的，因此，连续系统经过离散化后得到的模型具有一定的近似性，其近似程度与保持器的特点和采样周期有关。

图2.21 连续系统离散化

当连续系统用状态空间模型表示时，在状态方程的输入端、输出端分别加入虚拟的采样开关，并在输入采样开关后加上保持器，如图2.22所示。

图2.22 状态方程的离散化

2.3.2 时域离散化模型的推导

设连续系统的状态空间模型为

对状态方程进行拉氏变换，得

从而，有

对式（2.61）进行拉氏反变换，并利用卷积积分，得

式中

称为系统的状态转移矩阵，它描述了状态向量x（t）由初始状态x（0）向任一时刻t转移的特性。

式（2.62）是连续系统状态方程的解析解。下面对这个解进行离散化，以得到状态方程的等价离散化模型。

设采样周期为T（它相当于数值积分算法的计算步长h），考察kT及（k+1）T两个相邻采样时刻状态向量x（t）的值。将t=（k+1）T代入式（2.62），得

由式（2.63），有

于是，得

这就是两个相邻采样时刻kT和（k+1）T处状态向量x（kT）和x[（k+1）T]之间的基本递推关系。对上式右端积分项中的输入向量u（t）作不同的处理，可以得到不同的时域等价离散化模型。

假定对输入u（t）采用零阶保持器，如图2.23所示，则有

图2.23 零阶、三角形保持器对输入信号的处理

于是，式（2.65）变为

取t=τ-kT进行变量代换，则

式中

这样就得到了采用零阶保持器时的等价离散化模型

如果对输入u（t）采用三角形保持器，即假定在两个相邻采样时间之间输入u（t）按u（kT）与u[（k+1）T]连线斜率变化，如图2.23所示，则有

由于

故

代入式（2.65），整理可得

令

将其代入式（2.74），得到采用三角形保持器时的等价离散化模型

由于在式（2.70）和式（2.76）中，当采样周期T选定之后，Φ（T），Γ（T）及（T）都是常数矩阵，可以一次性计算出来，不再变化。这样，利用式（2.70）或式（2.76），在已知状态变量初值的情况下，可以很容易地求出在各采样时刻状态变量的值，进而求出相应的输出方程的值。

2.3.3 时域离散化模型的特性分析

将式（2.70）或式（2.76）代入测试方程式（2.37），有

而测试方程的解析解有如下关系

由此可见，时域离散相似算法是绝对稳定算法。

由图2.23可知，在输入端采样开关后加入零阶保持器相当于对输入u（t）作了零阶近似；加入三角形保持器相当于对输入u（t）作了一阶近似。因此，式（2.70）和式（2.76）的局部截断误差分别是O（T2）和O（T3），分别相当于欧拉法和RK2法的精度。

对于单输入单输出连续系统的状态方程，式（2.70）和式（2.76）每步的计算量主要是1次矩阵与向量的乘法（即n2次乘法），相当于欧拉法的每步计算量。

2.3.4 时域离散算法仿真实例

MATLAB的控制系统工具箱中提供了实现连续系统和离散系统之间相互转换的函数，最常用的是将连续系统转换成离散系统的c2d函数，其调用格式为：

        sysd=c2d（sysc，T，′method′）

其中，sysc和sysd分别为用结构变量描述的连续系统模型和离散系统模型；T为离散化步长；method为转换时所选用的离散化方法，其选项主要有以下几种：

●zoh采用零阶保持器（默认设置）；

●foh采用一阶保持器；

●imp采用脉冲响应不变法；

●tustin采用双线性变换法（Tustin法）

●prewarp采用改进的Tustin法；

●matched采用SISO系统的零极点根匹配法。

如果仅需对状态方程进行离散化转换，则其调用格式为：

        [Ad，Bd]=c2d（A，B，T）

其中，A和B分别为式（2.59）中状态方程的系统矩阵A和输入矩阵B，Ad和Bd分别为采用零阶保持器所得到式（2.69）中的Φ（T）和Γ（T）。

【例2.9】 考虑如图2.23所示的非线性控制系统，试用时域离散相似算法研究其单位阶跃响应（0≤t≤30s）。

图2.24 非线性控制系统

【解】 图2.24中的传递函数对应的状态空间模型为

在该模型前加上虚拟的采样开关和零阶保持器，如图2.25所示。

对应的离散化状态空间模型为

图2.25 连续状态空间模型的离散化

式中

相应的仿真模型如图2.26所示。

图2.26 离散化仿真模型

编制文件名为exam2-9.m的程序如下：

        % 这是例2.9的仿真程序
        clear
        A=[00；1-1]；B=[1；0]；C=[01]；D=0；  % 连续系统参数矩阵
        sysc=ss（A，B，C，D）；                     % 利用ss函数求取连续系统的模型
        T=0.05；                             % 置离散化步长T=0.05s
        sysd=c2d（sysc，T）；                      %用零阶保持器将连续系统转换成等价
                                                                %的离散模型
        Ad=sysd.a；Bd=sysd.b；Cd=sysd.c；Dd=sysd.d；% 求出离散模型的参数矩阵
        X=[0；0]；                            % 状态向量初始值
        yt=0；tt=0；
        R=1；                               %r（t）=1（t）
        M=30/T；                            % 设定仿真递推计算次数
        K=input（′请输入增益K=′）；
        fork=1：M
            E=R-yt（end）；                     % 求出输入与输出之间的差
          ifE＞=1                          % 处理饱和环节
                E=1；
          elseifE＜=-1
                E=1；
          end
          U=K*E；          % 控制信号
          X=Ad*X+Bd*U；   % 递推计算
          Y=Cd*X+Dd*U；
          yt=[yt，Y]；         % yt为记载各采样（kT）时刻输出响应的行向量
          tt=[tt，k*T]；       % tt为记载各采样（kT）时刻的行向量（与yt对应）
      end
      plot（tt，yt，′k′），grid       % 绘制结果曲线

【说明】在该程序中，结构变量sysc和sysd分别包含了连续系统状态空间模型（A，B，C，D）和离散系统状态空间模型（Ad，Bd，Cd，Dd）的所有信息。对于已知的状态空间模型sys，其参数矩阵A，B，C和D可以分别用指令sys.a，sys.b，sys.c和sys.d求出。

运行exam2_9.m后，并输入各种不同的增益K值，就得到如图2.27所示的仿真结果。

图2.27 非线性控制系统的阶跃响应