非連續(xù)系統(tǒng)的Simulink仿真方法研究

1、第18卷第7期 系統(tǒng) 仿 真 學(xué) 報© V ol. 18 No. 72006年7月 Journal of System Simulation July, 2006非連續(xù)系統(tǒng)的Simulink 仿真方法研究向 博1,高丙團1,張曉華1,胡廣大2(1.哈爾濱工業(yè)大學(xué)電氣工程學(xué)院,哈爾濱 150001;2.哈爾濱工業(yè)大學(xué)航天學(xué)院,哈爾濱 150001 摘 要:針對滑模變結(jié)構(gòu)控制、最優(yōu)控制等控制系統(tǒng)中出現(xiàn)的由于不連續(xù)性問題而直接影響Simulink 仿真速度和仿真效果的情況,分析了該現(xiàn)象產(chǎn)生的原因,并提出了相應(yīng)的解決策略。結(jié)合吊車防擺變結(jié)構(gòu)控制系統(tǒng),采用提出的解決方案進行了Simulink

2、仿真實驗研究,仿真實驗的結(jié)果證明所提出的策略對于解決這類不連續(xù)性仿真問題具有較好的效果。最后,給出了解決和避免這種問題的幾點建議。 關(guān)鍵詞:非連續(xù)系統(tǒng);Simulink 仿真;滑模變結(jié)構(gòu)控制;吊車中圖分類號:TP391.9 文獻標識碼:A 文章編號:1004-731X (2006 07-1750-05Simulink Simulation Research of Uncontinuous SystemXIANG Bo, GAO Bing-tuan, ZHANG Xiao-hua, HU Guang-da(1.School of Electrical engineering, Harbin In

3、stitute of Technology, Harbin 150001, China;2.School of Astronautics, Harbin Institute of Technology, Harbin 150001, ChinaAbstract: The Simulink simulation speed and effect are decreased because of the uncontinuous functions in the uncontinuous systems, such as sliding mode variable structure contro

4、l systems, optimal control systems and so on. The reasons caused the problem were analyzed in detail and corresponding resolving methods were also developed. Based on the presented methods, the Simulink simulation on variable structure control of an overhead crane system was studied, and the simulat

5、ion results demonstrate the effectiveness of the prompted Simulink simulation schemes. Finally, some suggestions on the problem were summed up.Key words: uncontinuous system; Simulink simulation; sliding mode variable structure control; overhead crane引 言 Simulink 是MathWorks 公司于1990年推出的產(chǎn)品,是用于MATLAB 下

6、建立系統(tǒng)框圖和仿真的環(huán)境。該環(huán)境剛推出時的名字叫Simulab ,由于其名字很類似于當時的一個很著名的語言Simula 語言,所以次年更名為Simulink 。從名字上看,立即就能看出該程序有兩層含義1,首先,“Simu ”一詞表示它可以用于計算機仿真,而“Link ”一詞表明它能進行系統(tǒng)連接,即把一系列模塊連接起來,構(gòu)成復(fù)雜的系統(tǒng)模型。正是由于它的這兩大功能和特色,使得它成為仿真領(lǐng)域首選的計算機環(huán)境。Simulink 是一個用來對動態(tài)系統(tǒng)進行建模、仿真分析的軟件包,它支持連續(xù)、離散及兩者混合的線性和非線性系統(tǒng),也支持具有多種采樣速率的多速率系統(tǒng)。其利用數(shù)學(xué)模型來獲取系統(tǒng)的一些重要特性參數(shù),這

7、些數(shù)學(xué)模型通常是以時間為變量的常微分方程來描述,并用數(shù)值方法進行計算機仿真求解。所以,Simulink 軟件包常被用于運動學(xué)仿真來反復(fù)求解 機構(gòu)運動約束方程2,通過積分獲得最終的速度(或加速度,收稿日期:2005-05-25 修回日期:2006-04-03作者簡介:向博(1980-,男,四川廣元人,碩士生,研究方向為基于滑模變結(jié)構(gòu)控制的吊車防擺系統(tǒng);張曉華(1961-,男,黑龍江哈爾濱人,教授,博導(dǎo),哈爾濱工業(yè)大學(xué)電氣工程學(xué)院常務(wù)副院長,研究方向為管道機器人,MATLAB 仿真技術(shù)。從而確定機構(gòu)運動的位置(和速度。圖1為運動學(xué)仿真的流程圖,它描述了諸如MATLAB 仿真工具箱Simulink

8、仿真軟件包的仿真步驟和流程2。圖1 運動學(xué)仿真流程框圖利用Simulink 仿真軟件包來求解機構(gòu)運動,比通常求解機構(gòu)運動的方法具有許多優(yōu)點。最顯著的優(yōu)點是位置求解問題可以用軟件包隱式求解,使用者只需為仿真提供適當?shù)某跏紬l件即可以求解機構(gòu)在任意時刻的位置問題。但是,如果在仿真模型中一旦出現(xiàn)了非連續(xù)模塊,則可能出現(xiàn)系統(tǒng)仿真停滯的情形。本文針對這一現(xiàn)象,以“基于滑模變結(jié)構(gòu)控制的橋式吊車防擺系統(tǒng)”為例,詳細分析了該現(xiàn)象產(chǎn)生的原因,并提出了相應(yīng)的解決策略。1 問題的提出本文以“基于滑模變結(jié)構(gòu)控制的吊車防擺系統(tǒng)”為例,來研究非連續(xù)系統(tǒng)的Simulink 仿真方法。2006年7月 向 博,等:非連續(xù)系統(tǒng)的S

9、imulink 仿真方法研究 July, 2006根據(jù)文獻3,吊車的物理模型如圖2所示:x 圖2 吊車的物理模型重物通過繩索與小車相連,小車在行走電機的水平驅(qū)動力f 的作用下在水平軌道上運動。小車的質(zhì)量為M ,重物的質(zhì)量為m ,繩索的長度為l ,小車與水平軌道的摩擦阻尼系數(shù)為D ,g 為重力加速度時間常數(shù)。利用拉格朗日方程,我們可以得到吊車的數(shù)學(xué)模型為: f Mx Dx mg x l g =+=+(1 1.1 滑模變結(jié)構(gòu)控制器的設(shè)計設(shè)1x =,2x = ,1l m =,29.8/g m s =,并取控制:a x =則根據(jù)式(1,可寫出如下二階系統(tǒng):12219.8xx x x a =+ (2 取

10、滑模面函數(shù)為:112s c x x =+ 10c > (3采用等速趨近率控制策略:sgn(dss dt= 0> (4 設(shè)李亞甫諾夫函數(shù):212V s =(5 則可得:sgn(0Vs s s s s =< (6 由于在切換面鄰域內(nèi)函數(shù)式(5是正定義的,而式(6中2s 的導(dǎo)數(shù)是負半定義的。所以,在s =0附近V 是一個非增函數(shù),因此,如果滿足條件式(6,則定義函數(shù)式(5是系統(tǒng)的一個條件李亞普諾夫函數(shù)。系統(tǒng)本身也就穩(wěn)定于條件0s =。系統(tǒng)將沿滑模面在t 趨于無窮大時歸于平衡點。由(2、(3、(4式變換可得:1121129.8sgn(a x c x c x x =+(7 式(2、(7

11、即構(gòu)成了具有自動防擺功能的二維滑模變結(jié)構(gòu)控制器。1.2 系統(tǒng)仿真針對前面所設(shè)計的滑模變結(jié)構(gòu)控制系統(tǒng)進行仿真,搭建的仿真模型如圖3所示。在此模型中,我們?nèi)坎捎肧imulink 工具箱中的模塊搭建模型。圖3 滑模變結(jié)構(gòu)控制器仿真模型取11c =,1=,1(01x =,2(00x =,仿真時間T=10s , 采用Variable-step 的ODE45算法。仿真結(jié)果如圖4所示:圖4 系統(tǒng)時域仿真結(jié)果由上圖可以看出,當系統(tǒng)仿到1s 左右時,速度非常慢,停滯不前,出現(xiàn)了圖形未完全仿真完成的情形,這顯然是我們不希望看到的。2 問題的分析為了更好的分析該問題產(chǎn)生的原因以及提出相應(yīng)的解決策略,首先,我們認為

12、有必要從Simulink 仿真運行原理開始進行分析。2.1 Simulink 仿真運行原理對Simulink 系統(tǒng)模型的仿真主要包含兩個階段4:初始化階段和運行階段1.初始化階段在初始化階段主要要完成以下工作:(1 每個模塊的所有參數(shù)都傳遞給MATLAB 進行求值,得到的數(shù)值作為實際的參數(shù)使用;(2 展開模型的層次結(jié)構(gòu),每個子系統(tǒng)被他們所包含的模塊替代,帶有觸發(fā)和使能模塊的子系統(tǒng)被視為原子單元進行處理;2006年7月系統(tǒng)仿真學(xué)報 July, 2006(3 檢查信號的寬度和模塊的連接情況,提取狀態(tài)和輸入/輸出依賴關(guān)系方面的信息,確定模塊的更新順序;(4 確定狀態(tài)的初值和采樣時間;完成這些工作后就

13、可以進行仿真了。2.運行階段初始化之后,仿真進入運行階段。仿真是由求解器控制的,它計算模塊的輸出,更新離散狀態(tài),計算連續(xù)狀態(tài)。在采用變步長求解器時,求解器還會確定時間步長。計算連續(xù)狀態(tài)包含下面兩個步驟:(1 首先,求解器為待更新的系統(tǒng)提供當前狀態(tài)、時間和輸入值。反過來,求解器需要狀態(tài)導(dǎo)數(shù)的值;(2 然后,求解器對狀態(tài)的導(dǎo)數(shù)進行積分,計算新的狀態(tài)的值。狀態(tài)計算完成后,再進行一次模塊的輸出更新,這時,一些模塊可能會發(fā)出過零的警告,促使求解器探測出發(fā)生過零的準確時間。Simulink的仿真過程是在Simulink求解器和系統(tǒng)相互作用之下完成的,系統(tǒng)和求解器在仿真過程中的對話作用如圖5所示。 圖5 系

14、統(tǒng)與求解器之間的對話在圖5中,求解器的作用是傳遞模塊的輸出,對狀態(tài)導(dǎo)數(shù)進行積分,并確定采樣時間。系統(tǒng)的作用是計算模塊的輸出,狀態(tài)進行更新,計算狀態(tài)的導(dǎo)數(shù),生成過零事件、從求解器傳遞給系統(tǒng)的信息包括時間、輸入和當前狀態(tài);反過來,系統(tǒng)為求解器提供模塊的輸出、狀態(tài)的更新和狀態(tài)的導(dǎo)數(shù)。2.2 過零檢測過零檢測5是通過在系統(tǒng)和求解器之間建立對話的方式工作,對話包含的一個內(nèi)容是事件通知,即系統(tǒng)告知求解器在前一時間步發(fā)生了個事件。事件由過零表示,過零在下列兩個條件下產(chǎn)生:(1 信號在上一個時間步改變了符號(含變?yōu)镺和離開O;(2 模塊在上一個時間步改變了模式(例如積分器進入了飽和區(qū)段。過零是一個重要的事件,

15、表征系統(tǒng)中的不連續(xù)性。例如,響應(yīng)中的跳變。如果仿真中不對過零進行檢測,可能會導(dǎo)致不準確的仿真結(jié)果。當采用變步長求解器時,Simulink 能夠檢測到過零(使用固定步長的求解器,Simulink不檢測過零。當一個模塊通知系統(tǒng)前一時間步發(fā)生了過零,變步長求解器就會縮小步長,即便絕對誤差和相對誤差是可接受的?？s小步長的目的是判定事件發(fā)生的準確時間。當然,這樣會降低仿真的速度,但這樣做對有些模塊來講是至關(guān)重要和必要的,因為這些模塊的輸出可能表示了一個物理值,它的零值有著重要的意義。事實上,只有少量的模塊能夠發(fā)出過零事件通知。每個模塊發(fā)出專屬于自己的事件通知,而且可能與不止一個類型的事件發(fā)生關(guān)聯(lián)。能夠產(chǎn)

16、生過零通知的Simulink模塊如圖6所示:圖6 能夠產(chǎn)生過零通知的SIMULINK模塊各個模塊的過零的類型是有差異的。例如,abs模塊在輸入改變符號時產(chǎn)生個事件,而saturation模塊能夠生成兩個不同的過零,一個用于下飽和,一個用于上飽和。signals&system庫中的HitCrossing模塊輸入穿過零點時產(chǎn)生一個過零,可以用來為不帶過零能力的模塊提供過零檢測的能力。一些過零只是用來通知求解器模式已經(jīng)發(fā)生了改變,另外一些則與信號相關(guān),用于觸發(fā)其他模塊。觸發(fā)包含三種類型:(1 上升沿信號上升到或穿過零或者信號離開零變正;(2 下降沿信號下降到或穿過零或者信號離開零變負;(3

17、雙邊沿上升或下降兩者之一發(fā)生。注意,Fcn模塊不支持過零,結(jié)果一些拐角點被漏掉了。2.3 算法的選擇Simulink提供了許多微分方程的解法6,如Runge-Kutta 法、Adams法、Gear法、Euler法、Linsim法等等。需要說明的是,對于不同的仿真問題,需要按照系統(tǒng)的特性,選擇不同的仿真方法。一般來說:Runge-Kutta方法(包含rk45、rk23變步長方法適合于高度非線性或不連續(xù)的系統(tǒng),不適合于剛性系統(tǒng)(即有快變特性又有慢變特性的系統(tǒng);Adams方法適合于非線性小、時間常數(shù)變化小的系統(tǒng);Gear方法是專門用于剛性系統(tǒng)仿真的方法,但對非剛性系統(tǒng)較差;Euler方法比較差,盡量

18、避免使用;Linsim方法適合于近似線性的系統(tǒng),對線性剛性系統(tǒng)有很大的優(yōu)越性。一般來說,使用變步長的自適應(yīng)算法是較好的選擇。這類算法會依照給定的精確度在各積分區(qū)內(nèi)自適應(yīng)地尋找各自的最大步長進行積分,從而使得速率最高。Simulink的變步2006年7月向博,等:非連續(xù)系統(tǒng)的Simulink仿真方法研究July, 2006長解法能夠把積分段分到足夠細以至到滿足精度要求的解。2.4 系統(tǒng)Simulink仿真停滯原因分析綜上所述,我們可以發(fā)現(xiàn):由于系統(tǒng)存在不連續(xù)模塊signum,所以,當系統(tǒng)于1s左右到達滑模面(s=0時,signum 模塊向系統(tǒng)發(fā)出過零通知。而當采用變步長求解器時, Simulin

19、k能夠檢測到過零。所以,當Simulink檢測到過零時,便自動縮小步長,可是,在下一仿真步里,系統(tǒng)繼續(xù)過零。如圖7和如圖8所示,因為滑模面在1s處不能正常歸零,所以signum模塊就反復(fù)過零,同時一直向求解器發(fā)出過零通知。求解器便相應(yīng)的一直不停的縮小步長。 圖7 滑膜面和signum模塊的時域響應(yīng) 圖8 1s左右時的滑模面和signum響應(yīng)的局部放大圖如圖9所示,系統(tǒng)大約經(jīng)過12個仿真步到達1s處時,步長急劇縮小至接近于零。這樣,由于仿真步長太小,系統(tǒng)便在不連續(xù)處形成了過多的點,超出了系統(tǒng)可用的內(nèi)存和資源,使得系統(tǒng)進展緩慢,仿真停滯不前。 圖9 系統(tǒng)仿真步長3 問題的對策通過以上分析,我們總結(jié)

20、、概括得出了以下系統(tǒng)仿真速度慢的原因:(1 系統(tǒng)方程中存在不連續(xù)函數(shù)sgn(s;(2 Simulink仿真模型中存在能夠產(chǎn)生過零通知的signum模塊;(3 采用的變步長求解器具有過零檢測并自動調(diào)整步長的功能。3.1 解決策略因此,基于以上原因,我們提出了以下4種解決策略:(1 取消Zero crossing detection功能;(2 采用fixed-step求解器;(3 采用不能夠產(chǎn)生過零通知的Fcn函數(shù)模塊;(4 柔化sgn(s函數(shù),使其連續(xù)化。以上所提出的4種解決策略中,采用單獨任何一種或幾種都可行。以下,我們嘗試采用第3種策略來觀察系統(tǒng)的仿真效果。即采用Fcn函數(shù)模塊來代替圖3中的

21、其他模塊所表示的數(shù)學(xué)關(guān)系,搭建的仿真模型如圖10所示: 圖10 采用Fcn函數(shù)模塊后的系統(tǒng)仿真模型對于函數(shù)fun,我們?nèi)?x,1x,a分別作為輸入u1, u2, u3。根據(jù)式(2可得函數(shù)fun的表達式為:-9.8*u2+u3對于函數(shù)fun1,我們?nèi)?x,1x,1c,分別作為輸入u1, u2, u3, u4。根據(jù)式(7可得函數(shù)fun1的表達式為: 9.8*u2-u3*u1-u4*sgn(u3*u2+u1取11c=,1=,系統(tǒng)仿真結(jié)果如圖11所示: 圖11 系統(tǒng)時域圖2006年7月 系 統(tǒng) 仿 真 學(xué) 報 July, 2006由于Fcn 模塊不支持過零,所以,系統(tǒng)在不連續(xù)的情況下仍然能迅速完成仿真

22、。除此之外,我們還采用了fixed-step 求解器和變步長下置Zero crossing detection 為off 的仿真方法,其仿真結(jié)果亦如圖11所示。可見,這前三種方法的仿真結(jié)果一樣。由上面的仿真結(jié)果可以看出:(1 取消系統(tǒng)的過零檢測功能之后,仿真速度快; (2 但是由于仍然存在不連續(xù)模塊signum ,所以加速度存在較大的抖振。3.2 sgn(s函數(shù)連續(xù)化 由于加速度存在劇烈的抖振,使得控制系統(tǒng)難于工程實現(xiàn),為了能消除系統(tǒng)因為不連續(xù)性而存在的抖振現(xiàn)象,同時加快仿真速度,且不影響系統(tǒng)的仿真效果,我們設(shè)計了第四種控制策略-sgn(s函數(shù)連續(xù)化。連續(xù)化前后的函數(shù)sgn(s函數(shù)和sat(s

23、函數(shù)如圖12所示: 圖12 sgn(s和sat(s函數(shù)其中,sat(s函數(shù)的表達式為:1(1s sat s s = <(8 為大于零的正數(shù),且。當取無窮小時,sat(s 函數(shù)便非常逼近sgn(s函數(shù)。采用sat(s函數(shù)來代替sgn(s函數(shù)后的系統(tǒng)仿真模型如圖13所示: 圖13 柔化sgn(s函數(shù)后的系統(tǒng)仿真模型根據(jù)式(8可得函數(shù)fun 的表達式為:sin(u/0.1414。其表示,s 時sat(s的值。仿真結(jié)果為:(I 當0=時,連續(xù)函數(shù)sat(s即轉(zhuǎn)化為符號函數(shù)sgn(s,仿真時將出現(xiàn)如前面所述的問題(仿真速度慢、抖振等。(II 當0.1=時,其它參數(shù)與圖3中相同,可得系統(tǒng)仿真結(jié)果如圖

24、14所示:圖14 采用連續(xù)函數(shù)(0.1=后的系統(tǒng)時域圖仿真結(jié)果說明:采用連續(xù)函數(shù)sat(s代替不連續(xù)函數(shù)sgn(s后,系統(tǒng)仿真速度加快,不影響系統(tǒng)控制效果,且加速度不存在抖振現(xiàn)象,易于工程實現(xiàn),滿足我們所期望的結(jié)果。(III 當0.0001=時,系統(tǒng)的仿真結(jié)果如圖15所示:圖15 0.0001=時的系統(tǒng)仿真結(jié)果仿真結(jié)果說明:當0.0001=時,加速度出現(xiàn)劇烈的抖振,性能接近理想的滑模變結(jié)構(gòu)控制。通過以上的仿真結(jié)果,我們可以總結(jié)如下:(1 采用sgn(s函數(shù)連續(xù)化的方法,系統(tǒng)仿真速度快; (2 越小,系統(tǒng)仿真結(jié)果越逼近原理想的滑模變結(jié)構(gòu)控制系統(tǒng),但抖振也越大,相應(yīng)的仿真速度亦越慢;反之,越大,仿

25、真結(jié)果失真也越大,但系統(tǒng)抖振越小,相應(yīng)的仿真速度亦越快。因此,我們建議:(1 如果希望既不影響系統(tǒng)的仿真速度,又不損失精度,0.0001為宜;第 18 卷第 7 期 2006 年 7 月 Vol. 18 No. 7 系 統(tǒng) 仿 真 學(xué) 報 表4 步長 h=0.1 求解例 2 在 x July, 2006 4 數(shù)值試驗 例 1（參見5） 常數(shù)延遲微分方程初值問題 = 10 處的整體誤差. x = 10 處整體誤差 1.6235864E-4 6.70327E-4 y '(x = 1000y(x + 997e y(x 1 +1000 997e 3x y(x =1+ e 3 3 0 x4 x0

26、 h=0.2 上面兩個例子驗證了我們所構(gòu)造的方法的收斂階， 并且 驗證了方法在求解剛性問題時是有效的。 問題的精確解是 y ( x = 1 + e 的整體誤差列在表 3 中。 表3 步長 h=0.1 h=0.01 3 x 。 用表 2 列出三階方法求解例 1，在積分區(qū)間0,4內(nèi)得到 參考文獻: 1 H Podhaisky, B A Schmitt, R Weiner. Design, analysis and testing of some parallel two-step W-methods for stiff system J. Applied Numerical Mathematics

27、(S0168-9274, 2002, 42: 381-395. 2 3 H. Podhaisky, R. Weiner, B.A. Schmitt, Two-step W-methods for stiff ODE systems J. Vietnam J. Math. (S0866-7179, 2002, 30: 591-603. Bernhard A. Schmitt , Rüdiger Weiner, Parallel Two-Step W-Methods with 4 Peer Variables J. SIAM Journal on Numerical Analysis(S

28、0036-1429, 2004, 42: 265-282. K J int Hout, M N Spijker. Stability analysis of numerical methods for delay differential equations J. Numer. Math. (S0006-3835, 1991, 59: 807-814. 5 6 曹學(xué)年，劉德貴，李壽佛. 求解延遲微分方程的 ROSENBROCK 方法的漸近穩(wěn)定性J. 2002，14(3：290-292. E Hairer, S P NØrsett, G.Wanner. Solving Ordinary

29、 Differential Equations I. Nonstiff problems M. Berlin: Springer-Verlag , 1993. 求解例 1 在區(qū)間0,4上的整體誤差 區(qū)間0,4上的整體誤差 1.1762236E-05 1.5437829E-08 例 2（參見6） y1' ( x = y1 ( x y2 ( x 1 + y2 ( x 10 ' y2 ( x = y1 ( x y2 ( x 1 y2 ( x ' y3 ( x = y2 ( x y2 ( x 10 100 y3 ( x 當 x 0 時， y1(x = 5 ， y2 ( x = 0.1 ， y3 ( x = 1 ，積分區(qū)間 是 0,10 ，我們用表 1 列出 GP-穩(wěn)定的二階方法求解例 2。 例 2 是文獻6的 348 頁中傳染病模型的例子， 我們在這 里增加了-100y3(x這一項，使問題成為剛性的。表 4 列出了 例 2 的計算結(jié)果。 （上接第 1754 頁） (