尹湛華

（廣東松山職業(yè)技術(shù)學(xué)院電氣工程系，廣東韶關(guān)512126）

一種適用于傳遞函數(shù)模型離散化的M-函數(shù)設(shè)計

尹湛華

（廣東松山職業(yè)技術(shù)學(xué)院電氣工程系，廣東韶關(guān)512126）

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

傳遞函數(shù)模型；離散化方法；MATLAB；程序設(shè)計

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

1 問題描述

MATLAB控制系統(tǒng)工具箱函數(shù)c2d的應(yīng)用局限性主要表現(xiàn)在兩個方面.

1.1 離散化方法不全

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

1.2 純滯后環(huán)節(jié)的近似處理有歧義

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

表1 c2d函數(shù)method的功能選項

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

為方便測試，采樣周期Ts取1 s，若離散化結(jié)果記作Hd，則c2d函數(shù)將以參數(shù)Hd.inputdelay存儲整步系數(shù)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

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

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

2 離散化程序設(shè)計

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

2.1 程序首段

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

2.2 純滯后環(huán)節(jié)預(yù)處理程序段

設(shè)連續(xù)系統(tǒng)傳遞函數(shù)模型為G（s）.程序開始時，先將G（s）分解成不含純滯后部分Gp（s）和純滯后環(huán)節(jié)e-τs兩部分，再將純滯后環(huán)節(jié)時間常數(shù)τ分解為整步部分NgTs和余量θ，即：

進而可得純滯后環(huán)節(jié)時間常數(shù)的整步系數(shù)Ng=floor（τ/Ts）和余量θ=τ-NgTs，其中floor（●）表示沿負無窮方向取整，借助MATLAB函數(shù)floor實現(xiàn).在程序總成段，整步部分z-Ng將作為因子直接代入最終的離散化結(jié)果中；而剩余部分e-θs則根據(jù)實際情況，采用Pad é逼近（缺省時，利用MATLAB函數(shù)pade實現(xiàn)）、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）；%整步系數(shù)

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 公式法直接轉(zhuǎn)換程序段

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

switch method

表2 純滯后環(huán)節(jié)時間常數(shù)的非整步余量θ近似處理方法

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）；%用符號函數(shù)復(fù)合工具compose代入

2.4 零極點匹配法程序段

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

表3 公式法直接轉(zhuǎn)換的部分常用離散化方法

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

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

[num，den]=numden（g0s）；

g0s_num=sym2poly（num）；

g0s_den=sym2poly（den）；

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

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

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

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；%匹配穩(wěn)態(tài)增益gpz=kz*gpz；

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

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

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

2.6 離散化結(jié)果總成程序段

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

再經(jīng)MATLAB函數(shù)sym2poly將符號表達式轉(zhuǎn)換為最終的離散模型脈沖傳遞函數(shù).

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

3 算例

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

表5 屏幕輸出算例離散化結(jié)果

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

4 結(jié)語

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

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

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

［3］薛定宇.控制系統(tǒng)計算機輔助設(shè)計——MATLAB語言與應(yīng)用[M].2版.北京：清華大學(xué)出版社，2006.

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

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

［6］王倩穎，吳斌，歐進萍.考慮作動器時滯及其補償?shù)膶崟r子結(jié)構(gòu)實驗穩(wěn)定性分析[J].工程力學(xué)，2007，24（2）:9-14，8.

［7］喻婭婭，王印松，許鵬春.數(shù)字控制中的保持器與系統(tǒng)穩(wěn)定性的關(guān)系[J].自動化技術(shù)與應(yīng)用，2004，23（11）:17-19.

［8］盧志強.連續(xù)時間系統(tǒng)課程的離散化方法綜述[J].中國科教創(chuàng)新導(dǎo)刊，2010（34）:85-86.

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

［10］陳雪梅，王萍.達林算法在z域和δ域的仿真[J].天津工業(yè)大學(xué)學(xué)報，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-），男，內(nèi)蒙古赤峰人，廣東松山職業(yè)技術(shù)學(xué)院電氣工程系副教授，碩士；研究方向：計算機測控技術(shù).