趙廣輝，高 鑫，方金福，張青青，于偉民，姚 冰

(1.遼寧科技大學(xué) 機(jī)械工程與自動(dòng)化學(xué)院,遼寧 鞍山 114051；2.營(yíng)口理工學(xué)院，遼寧 營(yíng)口 115000 3.沈陽(yáng)飛機(jī)設(shè)計(jì)研究所 飛行控制部,遼寧 沈陽(yáng) 110035)

在現(xiàn)代工業(yè)生產(chǎn)中，控制系統(tǒng)扮演著越來(lái)越重要的角色.由于信息技術(shù)的快速發(fā)展，控制系統(tǒng)中引進(jìn)離散控制方式執(zhí)行控制算法的做法越來(lái)越普及[1].由于連續(xù)控制系統(tǒng)設(shè)計(jì)有著完善的理論基礎(chǔ)，因而利用連續(xù)控制器理論設(shè)計(jì)離散控制器，使離散控制器能模仿連續(xù)控制器的特性是一種比較方便和流行的做法，這個(gè)過(guò)程也被稱作連續(xù)系統(tǒng)的離散化.離散化從不同控制系統(tǒng)角度來(lái)說(shuō)，分為固定步長(zhǎng)(Fixed Step)和可變步長(zhǎng)(Variable Step)兩種模式.常見(jiàn)的控制系統(tǒng)具有穩(wěn)定的控制計(jì)算周期，因而較多采用固定步長(zhǎng)模式，這也是本文所研究的系統(tǒng)控制模式.

為研究方便，一般的控制系統(tǒng)可以表示為輸出量的微分方程，其中微分方程的階數(shù)也是控制系統(tǒng)的階數(shù).為不失一般性，同時(shí)兼顧推導(dǎo)過(guò)程的簡(jiǎn)潔，本文以一階系統(tǒng)為例，研究不同離散化方法計(jì)算一階系統(tǒng)時(shí)域響應(yīng)的時(shí)間性能和離散化誤差，并進(jìn)行橫向?qū)Ρ龋瑥亩x擇出最優(yōu)方法.一階系統(tǒng)的等效微分方程為：

(1)

假設(shè)系統(tǒng)的s域傳遞函數(shù)為G(s)，這里為分析方便只需研究開(kāi)環(huán)響應(yīng)，那么系統(tǒng)的等效框圖如圖1所示.

圖1 開(kāi)環(huán)系統(tǒng)等效框圖

1 固定步長(zhǎng)離散化方法

1.1 歐拉法(Euler法/前向差分法)

根據(jù)z變換理論[2]，對(duì)于控制系統(tǒng)，從連續(xù)域(s域)到離散域(z域)的映射公式為

z=esT，

(2)

其中，T為采樣周期.

上式可作如下變換：

z=esT≈1+sT，

(3)

即

(4)

式(4)兩邊同時(shí)乘以輸出信號(hào)Y，并分別作各自頻域逆變換得到

(5)

那么

yk+1=yk+Ty′=yk+Tf(xk,yk).

(6)

式(6)即為歐拉(Euler)法近似公式，歐拉法也稱前向差分法.

以積分器為例，其傳遞函數(shù)為[3]

(7)

時(shí)域表達(dá)式為

y′=f(x,y)=x.

(8)

使用式(6)得到的積分器歐拉法離散化計(jì)算公式為

yk+1=yk+Txk.

(9)

圖2 使用歐拉法離散化積分器的面積近似圖

圖2為使用歐拉法離散化積分器的面積近似圖.其中，曲線表示連續(xù)系統(tǒng)中真實(shí)的x(t)，其與橫軸所圍面積代表積分器的實(shí)際輸出.柱狀圖代表x(t)在每個(gè)離散化時(shí)間點(diǎn)的采樣值，所有矩形的面積之和即為使用歐拉法近似的積分器離散化輸出面積.可以看出，歐拉法在每個(gè)離散化時(shí)間點(diǎn)都會(huì)產(chǎn)生誤差，并且隨著時(shí)間的推移誤差會(huì)逐漸累積.一般來(lái)說(shuō)，歐拉法的全局離散誤差[4]為0(T)，即誤差大小和離散化時(shí)間間隔線性相關(guān)，這在工程上是不可接受的.

1.2 雙線性變換法(Tustin法)

式(2)也可作如下變換

(10)

即

(11)

那么對(duì)于傳遞函數(shù)為G(s)的連續(xù)系統(tǒng)，其在z域的表達(dá)式為

(12)

式(12)即為雙線性變換法在頻域上從連續(xù)域到離散域的變換公式.雙線性變換法也稱Tustin法[5].

仍以積分器為例，將式(7)代入式(12)中，可得

(13)

對(duì)式(13)應(yīng)用z逆變換，時(shí)域表達(dá)式為

(14)

圖3為使用雙線性變換法離散化積分器的面積近似圖.折線圖代表x(t)在每個(gè)離散化時(shí)間點(diǎn)的采樣值，所有梯形的面積之和即為使用雙線性變換法近似的積分器離散化輸出面積.可以看出，和歐拉法的矩形近似相比，雙線性變換法的梯形近似離散化誤差更小，精度更高了.

1.3 休恩法(Heun法)

式(11)也可作如下變換

(15)

圖3 使用雙線性變換法離散化積分器的面積近似圖

式(15)兩邊同時(shí)乘以輸出信號(hào)Y，并分別作各自頻域逆變換得到

(16)

那么

(17)

式(17)中等號(hào)右邊的f(xk+1,yk+1)中yk+1為未知數(shù)，要繼續(xù)求解，需要一個(gè)估計(jì)值，歐拉法的解能夠滿足這一目的，將其代入式(17)中后，得到休恩(Heun)法的解：

(18)

對(duì)于式(7)中的積分器，使用休恩法進(jìn)行離散化的過(guò)程和雙線性變換法相同，因而離散化積分器的面積近似圖也和圖3一致.這是因?yàn)榉e分器的時(shí)域表達(dá)式(即式(8))僅與x有關(guān)，那么式(18)中關(guān)于pk+1的計(jì)算就可以忽略，故而使用休恩法對(duì)積分器進(jìn)行離散化得到的時(shí)域表達(dá)式與雙線性變換法相同.但是對(duì)于其它動(dòng)態(tài)環(huán)節(jié)(尤其是時(shí)域表達(dá)式與y有關(guān)的)，由于休恩法的第一步進(jìn)行了歐拉法近似估值，會(huì)導(dǎo)致其離散化誤差通常情況下比雙線性變換法要大一些，離散化精度相應(yīng)低一些，這一點(diǎn)在后面的仿真分析里會(huì)提及.

1.4 泰勒級(jí)數(shù)法

1.4.1 泰勒定理 設(shè)y(t)∈CN+1[t0,b]，且y(t)在不動(dòng)點(diǎn)t=tk∈[t0,b]處有N次泰勒級(jí)數(shù)展開(kāi)[6]：

y(tk+T)=y(tk)+TAN(x(tk),y(tk))+O(TN+1)，

(19)

其中，

(20)

式(20)中的y(j)(t)=f(j-1)(x,y)表示函數(shù)f關(guān)于t的(j-1)次全導(dǎo)數(shù).

N次泰勒方法的一般步驟為

(21)

其中，在每步k=0,1,…,M-1有

dj=y(j)(tk),j=1,2,…,N.

1.4.2 龍格-庫(kù)塔(Runge-Kutta)法 泰勒方法的優(yōu)點(diǎn)是最終全局誤差的階為0(TN)，并且可以通過(guò)選擇較大的N來(lái)得到較小的誤差.然而泰勒方法的缺點(diǎn)是，需要先確定N，并且要計(jì)算高階導(dǎo)數(shù)，它們可能十分復(fù)雜.每個(gè)龍格-庫(kù)塔(Runge-Kutta)方法都由一個(gè)合適的泰勒方法推導(dǎo)而來(lái)，使得其最終全局誤差為0(TN).一種折中方法是每步進(jìn)行若干次函數(shù)求值，從而省去高階導(dǎo)數(shù)計(jì)算.這種方法可構(gòu)造任意N階精度的近似公式，最常用的是N=4階龍格-庫(kù)塔方法(RK4).

標(biāo)準(zhǔn)的N=4階龍格-庫(kù)塔方法計(jì)算過(guò)程如下：

(22)

RK4適用于一般的應(yīng)用，因?yàn)樗浅＞_、穩(wěn)定，且易于編程.許多專家聲稱，沒(méi)有必要使用更高階的方法，因?yàn)樘岣叩木扰c增加的計(jì)算量相抵消.如果需要更高的精度，則應(yīng)該使用更小的步長(zhǎng)或某種自適應(yīng)方法.

1.4.3 Bogacki-Shampine法 更小的N階泰勒方法具有較低的精度，使計(jì)算過(guò)程更加簡(jiǎn)便，也更易于理解.N=2的2階龍格-庫(kù)塔方法(RK2)中的一種情形會(huì)退化為休恩法，計(jì)算過(guò)程如式(18)所示.

N=3的3階龍格-庫(kù)塔方法中最有名的是Bogacki-Shampine法(BS3)，其發(fā)明者Bogacki和Shampine聲稱這種方法是所有3階龍格-庫(kù)塔方法中性能最好的.

Bogacki-Shampine法計(jì)算過(guò)程如下：

(23)

2 性能對(duì)比分析

2.1 時(shí)間性能分析

為分析各離散化方法的計(jì)算效率，以常見(jiàn)的一階系統(tǒng)作為分析對(duì)象，考察各個(gè)方法在對(duì)一階系統(tǒng)進(jìn)行離散化時(shí)的時(shí)間性能，這里以一次迭代需要進(jìn)行的浮點(diǎn)數(shù)加法和乘法計(jì)算次數(shù)來(lái)衡量[7].

典型的一階系統(tǒng)的傳遞函數(shù)為[8]

(24)

時(shí)域表達(dá)式為

(25)

將式(25)代入式(8)、(18)、(22)和(23)分別可得應(yīng)用Euler法、Heun法、BS3法和RK4法離散化一階系統(tǒng)的計(jì)算過(guò)程.而將式(24)代入式(12)再做z逆變換可得應(yīng)用雙線性變換法離散化一階系統(tǒng)的計(jì)算過(guò)程.所有方法的時(shí)間性能對(duì)比如表1所示.

表1 不同離散化方法時(shí)間性能對(duì)比

根據(jù)表1中的時(shí)間性能可知，不同離散化方法的計(jì)算效率差別很大.考慮到進(jìn)行一次浮點(diǎn)數(shù)乘法的時(shí)間遠(yuǎn)大于進(jìn)行一次浮點(diǎn)數(shù)加法的時(shí)間，那么根據(jù)表1可得出下列結(jié)論：

(1)RK4法計(jì)算過(guò)程最復(fù)雜，耗時(shí)最多，但基本和BS3法不相上下；

(2)Heun法比Tustin法耗時(shí)稍多一點(diǎn)，二者所需時(shí)間差不多是RK4法和BS3法的1/5左右；

(3)Euler法耗時(shí)最少，所需時(shí)間是Heun法和Tustin法的一半左右.

2.2 離散化誤差分析

仍以一階系統(tǒng)為分析對(duì)象，下面比較各離散化方法的離散化誤差.式(24)所表示的一階系統(tǒng)的單位階躍響應(yīng)時(shí)域表達(dá)式為

h(t)=1-e-t/τ,t≥0.

(26)

其中輸入為單位階躍函數(shù)[9]

(27)

以式(27)所表示的單位階躍函數(shù)為輸入，分別使用5種離散化方法對(duì)式(24)所表示的一階系統(tǒng)進(jìn)行離散化計(jì)算，得到離散化的輸出，并與式(26)的預(yù)期結(jié)果進(jìn)行作差，即可得到不同離散化方法的離散化誤差.

這里假設(shè)采樣間隔T=0.01 s[11]，仿真時(shí)間為10 s，一階系統(tǒng)參數(shù)τ=2，式(26)的輸入經(jīng)過(guò)離散采樣后xk=u(t)|t=tk=1，初值y0=0.離散化誤差表達(dá)式為

err(tk)=|yk-h(tk)|.

(28)

圖4和圖5所示為不同離散化方法一階系統(tǒng)單位階躍響應(yīng)時(shí)域輸出以及離散化誤差對(duì)比示意圖.5種方法全局離散化誤差為各自誤差曲線與時(shí)間軸所圍面積，在仿真中的結(jié)果分別為0.48、8e-4、4e-4、1e-6和1e-9，從以上仿真可以得出如下結(jié)論：

(1)Euler法、Heun法、Tustin法、BS3法和RK4[12]法離散化誤差分別遞減，單點(diǎn)最大誤差均出現(xiàn)在仿真時(shí)刻第τ=2s處，數(shù)量級(jí)分別為10-3、10-6、10-7、10-9和10-12；

(2)Euler法的全局離散化誤差為0(T)，這在工程應(yīng)用中一般是較大的，因此實(shí)際控制系統(tǒng)已經(jīng)很少采用Euler法作為離散化方法；

(3)Heun法和Tustin法的全局離散化誤差均為0(T2)，只是前者的常系數(shù)較大而已.從仿真結(jié)果可以看出，Tustin法的離散化誤差大約是Heun法的一半左右，這與1.4.3中對(duì)兩種離散化方法的對(duì)比分析結(jié)論一致；

(4)BS3法[13]和RK4法的全局離散化誤差分別為0(T3)和為0(T4)，與理論分析一致；

(5)從單點(diǎn)離散化誤差來(lái)看，Tustin法和RK4法分別為10-7和10-12數(shù)量級(jí).值得一提的是，在實(shí)際控制系統(tǒng)中采用何種離散化方法還與所選用的浮點(diǎn)數(shù)精度有關(guān).如果離散化所選用的為32位浮點(diǎn)數(shù)，根據(jù)IEEE754[14]，其十進(jìn)制有效數(shù)字僅為7～8位，此時(shí)選擇Tustin法完全能夠滿足要求，選取精度更高的離散化方法已經(jīng)沒(méi)有必要；如果所選用的為64位浮點(diǎn)數(shù)，則其十進(jìn)制有效數(shù)字為15～16位，選擇RK4法雖不能完全覆蓋其最高精度，但在實(shí)際工程中10-12數(shù)量級(jí)的誤差基本可以等價(jià)于0，故而是可行的選擇方案.

3 結(jié)束語(yǔ)

本文介紹了連續(xù)控制系統(tǒng)中常用的5種固定步長(zhǎng)離散化方法，從時(shí)域和頻域兩方面對(duì)每種方法的計(jì)算過(guò)程做出推導(dǎo)，并對(duì)各種方法的異同給出了簡(jiǎn)單介紹.以一階系統(tǒng)為分析對(duì)象，從時(shí)間性能和離散化誤差兩方面對(duì)5種方法的性能做出對(duì)比分析.仿真分析結(jié)果表明：

(1)Euler法雖然時(shí)間性能良好，但由于離散化誤差較大而不能適應(yīng)控制系統(tǒng)的精度要求；

(2)Heun法無(wú)論在計(jì)算效率還是離散化誤差性能均不如Tustin法，因而可以完全被Tustin法代替；

(3)BS3法計(jì)算效率和RK4法基本相同，但全局離散化誤差卻高了很多，因而也應(yīng)該被RK4法所代替；

(4)從時(shí)間性能和離散化誤差兩方面來(lái)說(shuō)，在實(shí)際工程應(yīng)用中最廣泛使用的就是Tustin法和RK4法兩種離散化方法.兩種方法的綜合性能明顯優(yōu)于其它方法，但受浮點(diǎn)數(shù)精度制約，應(yīng)用場(chǎng)合應(yīng)有所不同.當(dāng)所使用的浮點(diǎn)數(shù)為二進(jìn)制32位時(shí)，采用Tustin法既能獲得很好的計(jì)算效率，又可完全滿足精度要求；而當(dāng)所使用的浮點(diǎn)數(shù)為二進(jìn)制64位時(shí)，采用RK4法可獲得最佳的綜合性能.

圖4 不同離散化方法一階系統(tǒng)單位階躍響應(yīng)

圖5 不同離散化方法離散化誤差對(duì)比

[1] 趙霞，王祝萍，賈海航.連續(xù)系統(tǒng)離散化方法的比較與解析初探[J].工業(yè)和信息化教育,2015,10:71.

[2] 胡壽松.自動(dòng)控制原理[M].北京:科學(xué)出版社,2011：319-320.

[3] Alan V，Oppenheim Alan S，Willsky with S,et al.Signals & Systems[M].(Second Edition)西安:西安交通大學(xué)出版社,1999：497.

[4] John H Mathews，Kurtis D Fink.數(shù)值方法(MATLAB版)[M].周璐，陳諭，錢方，等譯.北京:電子工業(yè)出版社,2010.

[5] 劉瑞葉，任洪林，李志民.計(jì)算機(jī)仿真技術(shù)基礎(chǔ)[M].北京:電子工業(yè)出版社,2011：62.

[6] 同濟(jì)大學(xué)數(shù)學(xué)系.高等數(shù)學(xué)[M].北京:高等教育出版社,2007：139-141.

[7] Randal E Bryant，David R O’Hallaron.深入理解計(jì)算機(jī)系統(tǒng)[M].北京:機(jī)械工業(yè)出版社,2011：534-536.

[8] Katsuhiko Ogata.現(xiàn)代控制工程(英文版)[M].北京:電子工業(yè)出版社,2014：122.

[9] 王正林，王勝開(kāi)，陳國(guó)順，等.MATLAB/Simulink與控制系統(tǒng)仿真[M].北京:電子工業(yè)出版社,2011：110.

[10] IEEE 754:Standard for Binary Floating-Point Arithmetic[C/OL].http://grouper.ieee.org/groups/754/,2012:118-121.

[11] 王向陽(yáng)，金海波.低比特率Bandelet域圖像壓縮編碼算法研究[J].遼寧師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2012，35(3)：321-323.

[12]白乙拉，李冰，馮景山.以氣溫為邊界條件的水庫(kù)冰蓋厚度的數(shù)值模擬研究[J].遼寧師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2012，35(2)：164-166.

[13] 楊忠志,徐珍珍.應(yīng)用Fukui函數(shù)探討雙苯基-取代的自由基閉環(huán)反應(yīng)的區(qū)位選擇性[J].遼寧師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2012，35(2)：193-195.

[14] 劉小丹，孫曉奇，沈?yàn)I.基于EMD的信號(hào)瞬時(shí)頻率估計(jì)[J].遼寧師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2009，32(1)：58-60.