一種適用于傳遞函數模型離散化的M-函數設計

尹湛華

（廣東松山職業技術學院電氣工程系，廣東韶關512126）

一種適用于傳遞函數模型離散化的M-函數設計

尹湛華

（廣東松山職業技術學院電氣工程系，廣東韶關512126）

針對單輸入單輸出線性時不變連續系統傳遞函數模型離散化問題，提出了一種獨立于MATLAB控制系統工具箱c2d函數的M-函數設計方法.根據各種常用離散化轉換方法的定義及其數學算式，在適當處理傳遞函數模型中純滯后部分的基礎上，著重運用MATLAB軟件的符號計算工具箱函數，編寫了具有多分支主控程序結構的實用M-函數離散化程序.最后通過一個具有代表性的算例，驗證了該M-函數程序的正確性和有效性，可借以提高設計結果的準確性和置信度.

傳遞函數模型；離散化方法；MATLAB；程序設計

連續系統離散相似法應用廣泛，如何將連續模型轉換為離散模型是其首要問題.MATLAB軟件在控制工程領域影響甚巨，舉凡作者教學過程中所用教材和其他參考文獻，只要涉及控制系統計算機輔助設計和仿真分析的連續模型（無論是控制對象還是控制器）離散化問題，幾乎都會用到“MATLAB軟件控制系統工具箱的c2d函數”.但作者在使用過程中發現該函數并不完善.迄今尚未檢索到任何一篇公開討論該函數應用局限性并提出相應解決辦法的文獻，因此本文設計了完全獨立于c2d函數且能適用于SISO系統連續模型各種常見離散化方法的MATLAB離散化程序.本文設計M-函數程序立足于符號計算功能相對較強的MATLAB7.6.0（R2008a），并在較晚近版本MATLAB8.3.0（R2014a）中測得一致結論.

1 問題描述

MATLAB控制系統工具箱函數c2d的應用局限性主要表現在兩個方面.

1.1 離散化方法不全

c2d函數的基本調用格式為sysd=c2d（sysc，Ts，method）.其中method用于說明離散化方法，其功能如表1所示（晚近版本MATLAB中method并不包含prewarp選項，而是改用opts選項并借助函數c2dOptions實現）.然而，各種差分法（含一階后向差分、一階前向差分和中心差分［1］等）也是連續系統模型離散化設計中常用的方法，c2d函數并不包含這些轉換方法，這在多種MATLAB論壇中都有反映，說明該函數所提供的離散化方法不盡全面.

1.2 純滯后環節的近似處理有歧義

在使用c2d函數進行離散化時，若純滯后時間常數是采樣周期的整數倍（可稱之為“整步采樣”），屏幕輸出離散化結果中純滯后因子指數值即為該整數，各種方法都不會出錯.反之，若非整步采樣，一方面，采用雙線性變換等方法時，程序運行后首先警告稱延遲時間被舍入至最接近采樣周期Ts，并建議使用ZOH法或FOH法，否則離散化結果是不精確的；另一方面，該函數代碼程序中調用了通過Pad é逼近近似處理延遲時間的矩陣指數函數expm，離散化時先將純滯后時間常數τ分解為整步部分NgTs和余量θ，形如：

表1 c2d函數method的功能選項

式中，Ng為純滯后因子的指數值（以下簡稱為“整步系數”），Ts為采樣周期.根據MATLAB幫助文檔的解釋，余量θ將按Pad é逼近近似處理后直接被對象傳遞函數中的非純滯后部分吸收，這意味著只要τ＜Ts，則必有整步系數Ng=0，但實驗證明并非如此.筆者就MATLAB幫助文檔中提供的一個實例進行了極限測試，該例所予連續模型的不含純滯后部分的傳遞函數表為：

為方便測試，采樣周期Ts取1 s，若離散化結果記作Hd，則c2d函數將以參數Hd.inputdelay存儲整步系數Ng，編寫測試程序代碼如下：

clc；close all；clear all；Ts=1；m=0；

Hd=c2d（tf（[1-1]，[1 4 5]，'inputdelay'，0.3），Ts，'zoh'）

while（Hd.inputdelay==0）

m=m+1；tol=1-10^（-m）；

H=tf（[1-1]，[1 4 5]，'inputdelay'，tol）；

Hd=c2d（H，Ts，'zoh'）

end

m，tol=vpa（tol），Hd.inputdelay

e=（1-vpa（tol，12））/1

運行這段程序，結果表明在采用ZOH方法進行離散化時，盡管Hd.inputdelay在連續對象的純滯后時間常數τ（即程序中的tol）達到10-11截斷誤差容限之前一直保持為0而未被近似為1，但只要0＜τ＜Ts，則屏幕返回結果總有整步系數Ng=1.這意味著，只要純滯后時間常數與采樣周期不完全整步，c2d函數ZOH方法的屏幕返回結果就會一律將整步系數Ng=τ/Ts沿正無窮方向取整，這顯然與幫助文檔的敘述不符；而用類似程序測試FOH方法，則結果與幫助文檔的描述基本一致.故這種情況下c2d函數的離散化結果也未足置信.

c2d函數還存在零極點匹配法并非全匹配（即未能實現無窮遠處零點的匹配），以及純滯后部分計算意義不明顯等問題，也會給工程應用帶來困擾［2］.鑒于這些問題，本文設計了一種獨立于控制系統工具箱c2d函數，適用于SISO連續系統傳遞函數模型各種常見離散化方法的M-函數.

2 離散化程序設計

以M-函數文件形式、多分支程序結構，運用MATLAB符號運算工具，在適當處理連續模型純滯后環節的基礎上，從定義出發，遍歷各種常用的離散化方法，通過總成得到最終的單輸入單輸出線性時不變連續系統傳遞函數模型離散化結果.

2.1 程序首段

采用M-函數文件形式編程，以C2D_for_SISO_LTI_Tf_Model為名定義函數，并用MATLAB的保留字varargin限定可變輸入變量使其參數數目介于2個到4個之間，否則報錯；輸出變量為Sys_d.該函數的調用格式與c2d函數基本一致，默認轉換方法也是ZOH.不同之處在于：（1）頻率預修正雙線性變換法的轉折頻率ωc既不出現在輸入變量列表中，也不是通過其他函數來定義，而是在程序執行過程中從鍵盤輸入.（2）純滯后環節的近似處理方式不同，第4個輸入變量appx用于選擇純滯后環節的近似處理方法.

2.2 純滯后環節預處理程序段

設連續系統傳遞函數模型為G（s）.程序開始時，先將G（s）分解成不含純滯后部分Gp（s）和純滯后環節e-τs兩部分，再將純滯后環節時間常數τ分解為整步部分NgTs和余量θ，即：

進而可得純滯后環節時間常數的整步系數Ng=floor（τ/Ts）和余量θ=τ-NgTs，其中floor（●）表示沿負無窮方向取整，借助MATLAB函數floor實現.在程序總成段，整步部分z-Ng將作為因子直接代入最終的離散化結果中；而剩余部分e-θs則根據實際情況，采用Pad é逼近（缺省時，利用MATLAB函數pade實現）、1階Taylor逼近（當θ＜＜Ts時）或全極點逼近（當θ接近Ts時）進行近似處理得Gd（s）［3-5］，然后與不含純滯后部分Gp（s）相乘將其吸收，得：

離散化方法為零極點匹配法時采用1階Pad é逼近，否則缺省情況下一概采用3階Padé逼近，如表2所示.該段程序核心代碼如下：

tol=sys_c.inputdelay；

[num_gps，den_gps]=tfdata（sys_c，'v'）；

syms s z t k；

Ng=floor（tol/Ts）；%整步系數

theta=tol-Ng*Ts；%非整步部分時間余量

if theta～=0

switch appx

case{'taylor'，'taylo'，'tayl'，'tay'，'ta'，'t'}

gds=vpa（simplify（1-theta*s））；%1階Taylor逼近

case{'allp'，'all'，'al'，'a'}

gds=expand（simplify（1/（1+theta*s+0.5*theta^2*s^2）））；%全極點逼近

otherwise

switch method

case{'matched'，'matche'，'match'，'matc'，'mat'，'ma'，'m'}

[num_gd，den_gd]=pade（theta，1）；%1階Pade逼近

otherwise

[num_gd，den_gd]=pade（theta，3）；%3階Pade逼近

end

gds=expand（simplify（poly2sym（num_gd，'s'）/poly2sym（den_gd，'s'）））；

end

else

gds=1；

end

g0s=simplify（（poly2sym（num_gps，'s'）/poly2sym（den_gps，'s'））*gds）；

2.3 公式法直接轉換程序段

對于前向差分法、后向差分法、中心差分法［1,6］、雙線性變換法和頻率預修正的改進雙線性變換法，可直接套用表3中連續域算子s到離散域算子z的轉換公式，利用符號計算工具箱復合函數compose直接代入G0（s）進行符號運算即可得到G0（z）.該程序段核心代碼如下：

switch method

表2 純滯后環節時間常數的非整步余量θ近似處理方法

case{'bwd'，'bw'，'bd'，'b'}

s=（（1-z^（-1））/Ts）；%后向差分法

case{'fwd'，'fw'，'fd'}

s=（1-z^（-1））/Ts/（z^（-1））；%前向差分法

case{'ctd'，'ct'，'cd'，'c'}

s=（1-z^（-2））/Ts/2/z^（-1）；%中心差分法

case{'tustin'，'t'，'tu'，'tus'，'tust'，'tusti'，'tst'，'tt'}

s=2*（1-z^（-1））/Ts/（1+z^（-1））；%雙線性變換法

case{'prewarp'，'pre'，'prew'，'prewa'，'prewar'，'pp'，'pw'，'pwp'，'prwp'，'p'}

wc=input（'Please input the critical frequency Wc（in rad/sec） wc='）；

s=wc*（1-z^（-1））/（1+z^（-1））/tan（Ts*wc/2）；%改進的雙線性變換法

end

gpz=compose（g0s，s）；%用符號函數復合工具compose代入

2.4 零極點匹配法程序段

根據定義，零極點匹配法不但要分別進行零點和極點全匹配，還要進行穩態增益匹配.因此，實際編程時還包括用多項式求根函數roots解得連續系統傳遞函數的零點和極點、通過公式實現從s域到z域的映射并用for循環實現分子因式和分母因式連乘，以及匹配穩態增益等3個步驟，實際設計程序時，零點和極點全匹配應使用：

表3 公式法直接轉換的部分常用離散化方法

其中zi（i=1，2，...，m）表示零點，pj（j=1，2，...，n）表示極點；而穩態增益匹配用：

實驗表明，MATLAB控制系統工具箱c2d函數的'matched'方法不是全匹配，未能實現無窮遠處零點的匹配，離散化結果中并不包含（z+1）n-m因子，進而穩態增益匹配結果也將有所不同.該程序段核心代碼如下：

[num，den]=numden（g0s）；

g0s_num=sym2poly（num）；

g0s_den=sym2poly（den）；

zero=roots（g0s_num）；%求連續傳函的零點

pole=roots（g0s_den）；%求連續傳函的極點

%符號計算進行零極點匹配

gps0=compose（g0s，0）；

m=length（zero）；

n=length（pole）；

gpz_num=1；gpz_den=1；

if m＞=1

for i=1∶m

gpz_num=gpz_num*（z-exp（zero（i）*Ts））；

end

end

if n＞=1

for j=1∶n

gpz_den=gpz_den*（z-exp（pole（j）*Ts））；

end

end

gpz=vpa（（z+1）^（n-m）*gpz_num/gpz_den）；%包含無窮遠處零點的完全匹配

gpz0=compose（gpz，1）；

kz=gps0/gpz0；%匹配穩態增益gpz=kz*gpz；

2.5 常用加保持器離散化方法程序段

對于脈沖響應不變法、加零階保持器法、加一階保持器法以及加三角形保持器法［7-8］，需根據變換方法的不同做一些特殊處理：因為e-sTs=z-1，故先可將各種保持器模型中有關e-sTs部分分離出來得Gm（z），然后將其余部分直接與G0（s）合并為Gp（s）再做處理，方法見表4.具體編程步驟可歸納為：（1）通過拉普拉斯反變換函數ilaplace將表4中的Gp（s）轉換為連續時域信號gp（t）.（2）用符號工具箱復合函數compose將gp（t）離散化為脈沖序列gp（k）.（3）用Z變換函數ztrans將gp（k）變換為gp（z），進而得Gp（z）=gp（z）Gm（z）.

MATLAB控制系統工具箱函數c2d并不包含嚴格意義上的“加一階保持器法”，其FOH方法實際上對應于表4中的“加三角形保持器法”［9］.該程序段核心代碼略.

2.6 離散化結果總成程序段

在求得Gp（z）的基礎上，可得離散化結果的符號表達式為：

再經MATLAB函數sym2poly將符號表達式轉換為最終的離散模型脈沖傳遞函數.

表4 各種加保持器離散化方法的特殊處理

3 算例

設某帶時間延遲一階模型（First-Order Lag Plus Delay，FOLPD）的傳遞函數為G（s）=e-1.5s/（4s+1）［10］，采樣周期分別選取Ts=0.5 s和Ts=1 s，依次采用MATLAB固有控制系統工具箱函數c2d和按本文方法設計的離散化M-函數C2D_for_SISO_LTI_Tf_Model遍歷各種離散化方法進行轉換.為便于比較，后者關于純滯后環節的余量近似處理方法均為缺省，即全部按Padé逼近方法處理.算例檢驗結果如表5所示.

表5 屏幕輸出算例離散化結果

當采樣周期選取Ts=0.5 s時，除零極點匹配法之外，MATLAB軟件固有函數c2d和本文M-函數C2D_for_SISO_LTI_Tf_Model對應離散化方法的執行結果并無不同，表明二者在整步采樣的情況下是等效的；但當采樣周期選取Ts=1秒時，結果則出現明顯差異.由于本文所采用的離散化轉換方法都是從嚴格的定義出發的，結果更可信.另外，前向差分法、后向差分法、中心差分法和加一階保持器法亦已在C2D_for_SISO_LTI_Tf_Model中全部呈現，而這些方法無法通過c2d函數實現.

4 結語

面向控制系統工程設計和仿真分析中常見的連續模型轉換為離散模型問題，本文設計了一種獨立于MATLAB控制系統工具箱c2d函數的M-函數離散化程序.該函數在采用適當方法近似處理單輸入單輸出連續定常系統傳遞函數模型純滯后環節時間常數的基礎上，借助MATLAB軟件的符號計算工具，從定義出發以多分支程序結構涵蓋多種常用離散化轉換方法，然后通過總成從計算機屏幕返回最終離散化結果.典型算例表明，使用該M-函數可以有效解決c2d函數所固有的方法不足以及算法歧義等問題，能夠提高離散化結果的準確性和置信度，具有實用意義.

［1］張東峰，惠晶，許勝.基于中心差分法的光伏陣列最大功率點跟蹤[J].電力電子技術，2012，46（9）:10-12.

［2］尹湛華.Dahlin數字控制器的兩種MATLAB輔助解法[J].南方金屬，2013，（4）:53-55，60.

［3］薛定宇.控制系統計算機輔助設計——MATLAB語言與應用[M].2版.北京：清華大學出版社，2006.

［4］尹先斌，周有訓.基于Taylor和Padé逼近的滯后系統IMC-PID研究[J].昆明理工大學學報（理工版），2006，31（2）:76-81.

［5］靳其兵，劉明鑫，馮春蕾.大滯后特性處理的研究和比較[J].北京化工大學學報（自然科學版），2008，35（6）:98-101.

［6］王倩穎，吳斌，歐進萍.考慮作動器時滯及其補償的實時子結構實驗穩定性分析[J].工程力學，2007，24（2）:9-14，8.

［7］喻婭婭，王印松，許鵬春.數字控制中的保持器與系統穩定性的關系[J].自動化技術與應用，2004，23（11）:17-19.

［8］盧志強.連續時間系統課程的離散化方法綜述[J].中國科教創新導刊，2010（34）:85-86.

［9］張紹寧，戴紅纓.MATLAB在數字控制系統設計中的兩則應用問題[J].計算機仿真，2003（Z1）:203-205.

［10］陳雪梅，王萍.達林算法在z域和δ域的仿真[J].天津工業大學學報，2004，23（1）:73-75.

The Design of A M-function Applying to Transfer Function Model Discretization

YIN Zhan-hua
（Department of Electrical Engineering，Guangdong Songshan Polytechnic College，Shaoguan 512126，Guangdong，China）

For the discretization of transfer function model of SISO-LTI continuous systems，the design method is proposed about a M-function which is independent of the C2D function of MATLAB control system Toolbox. According to the definition and mathematical formulas of a variety of common discretization methods，the practical discretization program of the M-function with multi-branch master control structure is coded mostly by the symbolic Toolbox functions of MATLAB software and properly based on preconditioning the time delay part of the transfer function model.Finally，the correctness and validity of the M-function is verified by a typical example which thereby can improve the accuracy and confidence relating to the designing results.

transfer function model；discretization methods；MATLAB；program designing

N945.12

A

1007-5348（2017）06-0004-07

（責任編輯：邵曉軍）

2017-02-19

尹湛華（1969-），男，內蒙古赤峰人，廣東松山職業技術學院電氣工程系副教授，碩士；研究方向：計算機測控技術.