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 非線性控制系統的階躍響應