基于分數階微積分的油氣懸架建模與試驗分析

孫會來，金 純，張文明，李 昊，田海勇

（1．北京科技大學 機械工程學院，北京 100083；2．中國北車股份有限公司，北京 100083）

油氣懸架具有良好的非線性阻尼特性，結構簡單，性能可靠，其特性對車輛行駛平順性和操縱穩定性有重要影響，在大客車、載重貨車以及工程機械車輛上被廣泛采用，特別是在載重自卸車領域應用更加普遍。針對油氣懸架建模及其運動特性國內外學者做出了大量研究，Concordia大學的研究人員［1］先后建立了單筒單氣室油氣懸架缸的參數化計算模型和非線性模型，采用基于能量的頻率相關等效線性化技術，其等效線性參數化模型被廣泛采用；Gao等［2］對重型車輛主動油氣懸架進行了建模及仿真；孫濤等［3］建立了油氣懸架長通孔紊流阻尼模型；楊波等［4］基于懸架非線性阻尼對雙氣室懸架進行了建模研究。在上述研究中建立的油氣懸架數學模型在考慮油液可壓縮的特性時，多認為油液含氣以后體積彈性模量發生變化。且同時做了氣體的質量保持不變，無泄漏，無形變等多項假設［5－8］。而黃夏旭等［9］研究又表明氣體溶解效應和油液可壓縮性在油氣懸架設計研究過程中不可忽略。

懸架系統的實際工作要更為復雜，不僅與本身參數有關，同時受到各處零部件相對運動位移、溫度、磨損情況、油液黏度及其在油氣懸架內流動及可壓縮特性的影響。以上諸多研究中總是側重于某幾項參數的變化及對懸架系統特性的影響，或者是通過增加計算項來保證計算精度，同時假設過多計算復雜難以從總體實現準確快速描述。

傳統整數階導數模型需要通過引入多個導數項和材料參數，本文提出根據油氣懸架多相介質的力學特點，引入分數階微積分理論對運動微分方程中的阻尼損耗項進行修正建模并分析其等效黏彈性，建立某型工程車輛油氣懸架的振動分數階振動方程并進行數值求解，探討在該模型及不同工況下懸架的振動特性。

1 油氣懸架結構及原理

單氣室油氣懸架主要由缸筒以及活塞桿和活塞組件組成，油氣懸架上端設有氣室或儲能器。油缸氣室與缸筒上端通過浮動活塞隔開，整個懸架缸內形成工作腔和環形腔，活塞桿壁上設有阻尼孔和單向閥。工作腔上部充入惰性氣體，下部及環形腔充入油液。圖1是根據上述描繪簡化后的油氣懸架物理模型。

圖1 單氣室油氣懸架原理Fig．1 Single-chamber pneumatic suspension principle

圖2 油氣懸架實車測試Fig．2 Pneumatic suspension substantial vehicle test

當懸架處在壓縮行程時，活塞桿向上運動，工作腔惰性氣體體積減小，壓力升高，油液受到壓縮經單向閥和阻尼孔流入環形腔，使氣室容積減小，氣壓力升高。在這此過程單向閥及阻尼孔會同時連通工作腔與環形腔，油液流量大，產生的油液阻尼力比較小，主要由氣室內氣體受到壓縮產生彈性作用來抑制活塞桿的向上運動。當懸架處在伸張行程時，活塞桿向下運動，工作腔的體積增大進而壓力減小，環形腔中的油液流向工作腔，此時單向閥關閉，油液只流經阻尼孔，流速較高，產生較大的阻尼力，抑制了活塞桿的向下運動，產生的阻尼力較大，有較強衰減振動的作用。工程車輛油氣懸架安裝如圖2所示，油氣彈簧懸架與輪胎組成振動系統起衰減振動、減緩沖擊的作用。

經過如上分析及在試驗中發現油液黏度及油液在油氣懸架內流動、氣體在油液中的溶解率的隨溫度和壓力而變化等這些都不同程度的影響到了油氣懸架特性。油氣懸架阻尼特性等同于純黏性體的計算精度有待提高，往往需要涉及過多參數。油氣懸架特性和外力大小，溫度的改變、力的作用時間及加載歷史都有關系，這些符合黏彈性物質的力學特點。黏彈性材料的力學性質既表現出彈性，同時又有黏性特性，介于理想彈性體和牛頓流體之間，其應力－應變響應依賴于時間和應變率有記憶性，與荷載和變形歷史有關。而分數階黏彈性微積分模型僅用較少的參數便能準確地描述大量復雜黏彈性材料廣泛頻率范圍內的動力學特性［10］，能較準確地描述黏彈性材料的動力學特性。

基于以上思想本文將油氣懸架等效為黏彈性系統進行建模，引入分數階模型描述阻尼損耗項對油氣懸架系統進行分析，并進一步確定分數階次，觀察在不同激勵下分數階建模與整數階區別。

2 油氣懸架分數階模型

2．1 分數階微積分定義

分數階微積分即非整數階微積分，數學家從不同的角度入手，給出了分數階微積分不同的定義。每種定義有其表達式與性質，其中被多數人所采納的定義形式是 Riemann-Liouville（R-L）分數階微積分如下［11］：

式中：Γ為伽馬函數；n為任意整數；aDqt為分數階微分算子；q為分數階階次；a為積分算子下限；t為積分算子上限。

在工程中應用較多及本文采用的是更為方便的Caputo定義，這種計算形式避免了極限求解，且初始條件以整數階微積分的形式給出，其定義如下：

（R－L）分數階微積分定義與Caputo分數階微積分定義都是在時域范圍內對函數f（x）的定義。（RL）的Laplace變換涉及到分數階積分的初始值與原（R－L）分數微分的初始值，雖可求出含（R－L）分數階次導數的解，但很難解賦予合理的物理意義。而Caputo分數階微積分定義的最大優點是其初始值與整數階微積分一樣，有明確的物理意義，具有超奇異性。

2．2 油氣懸架的分數階模型

如前所述，影響懸架輸出特性的因素眾多，試驗條件也無法做到規避所有干擾項，建立數學模型和模擬分析前，都需根據具體情況做出相應假設。本文研究油氣懸架阻尼力等效分數階黏彈性的表現，可借助筒式液力減振器的研究經驗［12］，對工作參數做了適當的簡化和忽略，做出如下假設：

（1）懸架中活塞、活塞桿與缸筒之間不產生泄漏；

（2）液壓油的重力遠小于其工作情況下的運動力，故不計工作油液的重力勢能；

（3）不計壓力變化而引起系統構件的彈性變形；

（4）被研究封閉區域內同一瞬時壓力處處相等；

（5）不考慮油氣懸架內部溫度變化對性能的影響，即假設其工作工程中是恒溫過程。

油氣懸架的受力分析對象為活塞桿筒組件，活塞桿在內部做往復運動。油氣懸架存在壓縮運動和拉伸運動，這兩種運動中流體的運動方向是相反，因此在對應的兩種運動狀態下對活塞桿筒組件進行力平衡分析。油氣懸架在受到外部激勵時，活塞桿及活塞組件和缸筒之間要產生相對運動。從輸出力的角度出發整體來分析其振動特點。油氣懸架在受到外部激勵時，活塞桿及活塞組件和缸筒之間產生相對運動時，其輸出作用力主要由氣體彈性力、油液阻尼力及活塞和缸筒之間的摩擦力組成。輸出合力公式如下：

式中：F為油氣懸架缸筒輸出力；Fk為氣體彈性力項；Fc為油液阻尼力項；Ff為活塞和缸筒之間密封件的摩擦力項

對于彈性力、阻尼力及摩擦力項的計算，先取基本整數階常用計算方法［13］進行說明，彈性力項為：

式中：P為工作腔氣體壓力 ；A3為活塞桿面積；M為流體密度；γ為氣體多變指數；V0為靜平衡位置氣體體積；x為活塞和缸筒間相對位移。

油氣懸架密封件的摩擦力使得油氣懸架的動靜太特性呈現出遲滯性，物理特性類似阻尼可以借助液壓缸及活塞桿的密封件摩擦阻力的算法來計算：

式中：μ1為密封圈摩擦系數；ΔP為密封圈兩側壓力差；D為液壓缸缸筒內徑d為活塞桿外徑；bD，bd為活塞、活塞桿密封圈寬度；kDkd為密封圈摩擦修正系數。

油氣懸架阻尼力項根據懸架阻尼孔的節流孔口形式，選取薄壁孔口類型的阻尼流量計算方法如下：

式中：Cz，Az為阻尼孔的流量系數、過流面積；Cd，Ad為工作腔氣體壓力；ρ為油液密度。

sign為符號函數，根據輸入信號的速度方向取值，壓縮行程為1，伸張行程－1。

先只考慮壓縮行程即sigx＝1。聯立方程（4）～（6）可以得到懸架輸出力整數階方程，便于計算進行無量 綱 化 處 理，令 a1＝ （A1－A3）3r／［2n2（CzAz＋CdAd）2］，a2＝Mg·Vγ0，a3＝V0，a4＝A3，得到簡化系統方程：

根據圖1運行原理及其幾何關系可知：

式中：z（t）為懸架上端輸出位移；y（t）為懸架下方輸入激勵位移。

上式描繪了三者位移及加速度關系，在油氣懸架的下方施加一個正弦位移激勵 y（t）＝A sinωt作為輸入，相當于路面激勵，A為激勵振幅，ω為激勵頻率。

聯立式（3）～（9），在整數階方程基礎上對于耗散阻尼項引入分數階參數，同時計入自身重力，建立系統等效分數階 Bagley-Torvik阻尼振動方程［14］，進一步得到油氣懸架輸出力的非線性分數階方程：

式中：D為分數階微分算子；α為分數階算子階次。

方程（10）為分數階非線性常微分方程，彈性力項中變量x位于分母上，阻尼力項亦為平方是非線性項，這使得懸架的振動輸出力是非線性的。大部分分數階微分方程難以求得解析解，此類分數階微分方程的解析解也只能由比較特殊的函數來表示，針對工程需求一般求取此類方程的數值解。

目前分數階微積分常用的計算方法包括解析法、數值法和濾波器算法。本文采用改進的Oustaloup濾波器［15］算法來仿真未知信號的分數階微積分，其原理是用高階的傳遞函數近似替代分數階算子。Oustaloup算法就是在頻率段內，用泰勒展開并取一階近似得到傳遞函數。另一方面，考慮到濾波器分子和分母階次的一致性會導致仿真過程出現代數環，故而進行改進，后面再接一個低通濾波器。這樣可編寫濾波器零極點及增益函數，在Simulink中建立起如圖3的分數階微分器模塊。進而搭建懸架系統微分方程，求得數值解。

圖3 分數階微分濾波器Fig．3 Fractional differential filter

3 懸架分數階模型數值求解和試驗

3．1 分數階求解及試驗對比

自行設計油氣懸架試驗臺能獲取油氣懸架輸出力特性的動態數據，方便于對不同工況下的工作特性進行分析和比較，驗證油氣懸掛分數階模型的正確性。文中試驗使用的試驗臺其結構圖如4所示，主要由驅動裝置、懸架系統、傳感器、采集系統和簧上負載組成。設計思路主要是參照1／4車的振動模型，模擬車輛實際運行過程中油氣懸架的工作狀態，測試其受到外部激勵和載荷時的輸出特性。

整個臺架驅動力由下端地基平臺上的變頻三相異步電動機發出，經減速器后將力傳導至與支架接觸的偏心輪上。偏心輪與傳動梁保持接觸，傳動梁通過關節軸承與懸架活塞桿鉸接。偏心輪隨電機勻速旋轉同時驅動傳動梁做上下方向的擺動，擺動傳遞至懸架缸活塞桿與傳動梁鉸接處按照幾何關系即為豎直方向的正弦激勵。為防止試驗時質量塊左右晃動，在質量塊兩側設計了限位支架。本設計激勵頻率可通過調節電機頻率來隨時改變；激勵幅值可通過調節偏心輪的偏心距來方便的調整；懸架簧上質量可通過改變質量塊數量來調整大小。

對于橫梁臂自身質量及質量塊質量需通過杠桿比進行等效，即為簧上質量部分，通過氣室靜態氣壓推算簧上質量予以校正。試驗臺系統構件引起的彈性變形以及系統自身各處阻尼較小可以忽略不計。試驗臺下放為鋼結構試驗臺地基，地基自身振動引發的誤差極小也可忽略不計。每次試驗前需對偏心距進行校正，保證激勵振幅準確性。傳動梁及橫臂梁的杠桿比問題，需首先參照幾何關系進行等效轉換和試驗校正。

圖4 油氣懸架試驗臺Fig．4 Test stand of hydro-pneumatic suspension

為驗證分數階理論的準確性與整數階對比，完整描述試驗過程及結果，本文同時根據試驗臺結構設計構建了臺架三維模型如圖5，進而將模型導入多體動力學軟件ADAMS／View中，通過設置各零部件的材料屬性并添加零部件約束及力關系，與LMSAmesim液壓軟件聯立得到試驗臺多剛體整數階仿真模型，目的是作為整數階結果與論文的分數階模型進行對比。

試驗選用懸架基本參數如表1所示，為某噸位自卸車實車懸架。研究主要對比相同工況下懸架缸的輸出特性并針對不同工況展開試驗，即簧上載荷的振動加速度、懸架缸輸出力幅值和頻率、運動行程等。

圖5 試驗臺架軟件仿真Fig．5 Software simulation of test stand

表1 試驗用懸架基本參數Tab．1 The basic parameters of test suspension

本文針對偏心輪偏心距為5 mm，轉動頻率為1 Hz正弦激勵情況先進行分析，這也是工程車輛常見激勵頻率。取多個不同分數階次于振動穩態后在時域范圍內進行對比分析發現，與整數階次結果的周期和相位基本一致，而幅值差異明顯，可見分數階次的選取對振動結果的描述有很大影響。為了最為吻合試驗數據，分析前述方程，將最優階次的選取轉化為單變量的最優數學問題［16］，單變量為分數階次α，根據計算模型與實測數據間誤差最小原則進行數值求解選取最優階次。定義誤差評價函數E為計算結果與實測數據的差方均值，這樣更為直觀方便，取毫米為單位。如下：

式中：E為誤差評定；n為對比采樣點個數；xli為整數階結果第i個采樣點位移；xi為實測結果第i個采樣點位移。

同樣x2i為分數階結果中第i個采樣點位移，代入后則為分數階計算結果與實際測量誤差。約束條件的選取設在頻段［0．4，1．4］內進行，迭代步長為 0．001，觀察最優階次變化后目標函數的變化如圖7。迭代求解得分數階建模在分數階次取0．752～1的范圍內比整數階結果誤差要小，在1處即為整數階結果同為1．68的差方均值。迭代結果最小誤差即最優值出現在0．912處，其幅值結果最接近實測結果。下文以0．912次分數階結果、整數階次仿真結果、實測數據進一步對比分析。

圖6 差方均值隨迭代階次變化Fig．6 Error changes with the order

實測數據是懸架穩定振動后測得，為了便于說明對比圖的繪制中，實測數據相位滯后或提前90°，取Matlab分數階數值結果、ADAMS整數階結果及實測結果合并對比，如圖7為懸架缸缸筒相對位移對比，同取在4 s時域范圍內對比，后續分析也如是。

分析圖7可知缸筒與活塞桿的相對位移中分數階建模結果更加貼近于實際情況，整數階結果振幅比試驗數據振幅要更大，最大振幅要超過實測值約1．5 mm。在分數階原理上，對于阻尼項階數與運動關系可知，隨著分數階次的增大振動會越來越激烈，振動位移幅度也會越大，即前面已有的運動狀態對下一時刻影響力減弱。反之越小于1，則增加這種相互關系，使得阻尼力變大，振幅減小。明顯分數階建模在幅值方面更為精確，尤其是在位移下止點位置。受阻尼力影響，活塞壓縮和伸張行程是一個不重合過程。在伸張行程，在阻尼的作用下使得其伸張行程要更小，分數階的引入強調了阻尼的非線性，分數階模型結果在伸張后半程和壓縮前半程的優勢更為明顯，最大值最小值的出現時間與實測結果更為相近，改變趨勢更為吻合。

圖7 缸筒相對位移對比Fig．7 data of displacement curve comparison diagram

圖8 缸筒加速度對比Fig．8 Accelerated speed of cylinder comparison

圖9 懸架缸輸出力曲線對比Fig．9 Output force of suspension cylinders comparison

圖8 為懸架缸缸筒加速度曲線，與圖6相似，頻率和幅值都很接近，明顯仿真結果的曲線比試驗結果更為平滑，這是因為在整數階仿真模型及軟件中把試驗臺當作剛體處理，沒有考慮傳動梁等構件自身的阻尼及多處振動，但其對結果影響不大。

圖9為懸架缸輸出力曲線，動平衡后油氣懸架的輸出力呈周期性變化，結果曲線幅值變化頻率與實測數據相近。油氣懸架本身的高度耦合多復雜工況使得黏彈性特性出現且最終趨向穩態。在壓縮行程分數階模型優勢并不明顯，而在伸張行程分數階建模遠比整數階結果更貼近實測數據，在伸張行程主要是阻尼力在起作用，分數階模型的引入是針對阻尼力進行了調整，其優勢在伸張行程體現明顯，阻尼力的調整同樣使得輸出力較之于整數階相位上有延遲，其輸出力最大值出現要晚，這與實測數據更為相近。

3．2 多工況下分數階次與誤差變化

為了進一步驗證分數階建模在油氣懸架多工況中具有廣泛的有效性，可改變激勵振幅和頻率，觀察最優分數階次的變化，同時對比其與整數階建模誤差變化。本文通過控制偏心輪距和激勵電機的轉速改變激勵頻率和振幅，進行多次試驗，并參照分數階結果差方均值最小為目標進行最優階次的確定。在同一振幅5 mm改變激勵頻率，在同一頻率1 Hz下改變振幅分別進行最優化求解和試驗，通過定量計算擬合曲線得到圖10最優階次隨激勵頻率變化以及圖11隨激勵振幅變化。

分析圖10和11發現隨著頻率的改變最優階次下降并不大，高頻激勵下最終穩定趨于0．9，而隨著激勵振幅的變化下降較大趨于穩定0．86，這是因為在高幅快速振動下阻尼力起到了更為明顯的作用，其等效黏彈性明顯，二者最后都趨于穩定，這說明分數階建模的取值在高幅階段可以穩定選取，建模方便直接。

圖10 最優階次隨激勵頻率變化Fig．10 Optimal order changes with the excitation frequency

圖11 最優階次隨激勵振幅變化Fig．11 Optimal order changes with the excitation amplitude

圖12 不同振幅下差方均值對比Fig．12 Error comparison at different amplitude

為了進一步分析和觀察分數階最優階次確定后與整數階建模的差別，可對二者按評價指標進行觀察。取多組不同振幅下數據進行評價然后擬合繪圖得到圖12即不同振幅激勵下差方均值變化，差方的平均值越小說明精確度越高，可見隨著振幅增大二者誤差都有進一步的增大，但是整數階誤差增大更為明顯，進一步驗證了分數階建模的有效性。另一方面，取多組不同頻率下試驗數據擬合可得不同頻率下差方均值對比，因最優階次下降不大，二者誤差并不明顯，但在總體范圍內變頻和變幅值的分數階建模誤差都要小于整數階。結合圖11可知等效黏彈性的分數階建模方式能進行有效修正，得到更準確的結果。

圖13 阻尼力對比Fig．13 Damping force comparison

分析分數階建模的特點可以對阻尼力項進行說明，圖13為初始試驗條件即1 Hz和5 mm激勵下二者數值計算的阻尼力項對比，在伸張行程分數階可以達到穩定的550 N，而整數階結果只有400 N，分數階阻尼對系統消耗更大，而二者彈性力項變化不大，原理上這造成了懸架的相對位移非標準正弦波，在阻尼力的非線性作用下使得向下運動的伸張行程偏大，向上運動的壓縮行程偏小。

4 結 論

（1）根據油氣懸架工作原理，突出油氣懸架的多相介質特點，采用等效黏彈性阻尼的分析方法，引入分數階導數建立了油氣懸架的 Bagley-Torvik方程。利用Oustaloup算法對分數階非線性微分方程進行數值求解得到油氣懸架分數階輸出力，位移，加速度等多項數值結果。

（2）將分數階、整數階仿真與試驗結果進行對比，驗證了分數階模型的正確性。在1 Hz，5 mm激勵下測試表明取分數階次0．912時與實驗結果最為相近，該簡潔模型能較精確地描述油氣懸架的特性。

（3）改變激勵振幅和頻率，在不同工況下進行對比分析最優階次的變化，在高頻激勵下最優階次趨于0．9，在高幅激勵下趨于0．86，在整個測試激勵范圍內分數階模型誤差都要小于整數階，這說明了分數階建模的有效性。

［1］Pazooki A，Cao Dong-pu，Rakheja S．Experimental and analytical evaluations of a torsio-elastic suspension for off-road vehicles［J］．SAE International Journal of Materials and Manufacturing，2010，3（1）：326－338．

［2］Gao B，Darling J，Tilley D，et al．Control of a hydropneumatic active suspension based on a non－linear quarter-car model［J］．Journal of Systems＆Control Engineering，2006，220（1）：15－31．

［3］孫濤，喻凡，張振東．重載車輛油氣懸架不確定分析及H∞控制器設計［J］．振動與沖擊，2007，26（9）：51－54．SUN Tao，YU Fan，ZHANGZhen-dong．Uncertainty analysis for a hydro-pneumatic suspension system of heavy duty vehicles and H∞controller design［J］．Journal of Vibration and Shock，2007，26（9）：51－54．

［4］楊波，陳思忠，王勛，等．雙氣室油氣懸架特性研究［J］．機械工程學報，2009，45（5）：276－280．YANG Bo，CHEN Sizhong，WANG Xun，et al．Research of Twin accumulator Hydro－pneumatic Suspension［J］．Journal of Mechanical Engineerin，2009，45（5）：276－280．

［5］王增全，申焱華，楊玨．連通式油氣懸架數學模型及特性分析［J］．農業工程學報，2012，28（5）：60－65．WANG Zeng-quan，SHEN Yan-hua，YANG Jue．Mathematical model and characteristics analysis of interconnected hydropneumatic suspension［J］．Transactions of the Chinese Society of Agricultural Engineering，2012，28（5）：60－65．

［6］徐道臨，張林，周加喜．重型礦用自卸車油氣懸架參數優化［J］．振動與沖擊，2012，31（24）：98－101．XU Dao-lin， ZHANG Lin， ZHOU Jia-xi． Parametric optimization of hydro pneumatic suspension of a heavy mining dumper［J］．Journal of Vibration and Shock，2012，31（24）：98－101．

［7］黃夏旭，申焱華，楊玨，等．基于集中參數熱模型法的自卸車油氣懸架系統熱分析［J］．農業工程學報，2013，29（10）：64－70．HUANG Xia-xu， SHEN Yan-hua， YANG Jue， et al．Thermal analysis of hydro-pneumatic suspension system for dumper based on a lumped-parameter thermal model［J］．Transactions of the Chinese Society of Agricultural Engineering，2013，29（10）：64－70．

［8］張輝，何景峰，韓俊偉．油氣懸架加載試驗臺正弦波復現控制方法［J］．農業機械學報，2009，40（7）：35－39．ZHANG Hui，HE Jing-feng， HAN Jun-wei． Sine wave replication control methodon hydro pneumatic suspension force loading system［J］．Transactions of the Chinese Society for Agricultural Machinery，2009，40（7）：35－39．

［9］黃夏旭，楊玨，申焱華，等．基于氣體溶解與油液可壓縮性的油氣懸架性能研究［J］．農業機械學報，2013，44（6）：14－18．HANG Xia-xu， YANG Jue， SHEN Yan-hua， et al．Characteristics analysis of hydro pneumatic suspension based on gas dissolution and oil compressibility［J］．Transactions of the Chinese Society for Agricultural Machinery， 2013，44（6）：14－18．

［10］ Song D Y，J TQ．Study on the constitutive equation with fractonal derivative for viscoelastic fluids－modified Jeffreys models and its application［J］．Rhelo Acta 1998，21（37）：512－517．

［11］Connolly J A．The numerical solution of fractional and distributed order differential equations［D］．University of LIverpool，2004．

［12］ Kurfess，Thomas R，Messman，et al．A methodology to integrate a nonlinear shock absorber dynamics into a vehicle model for system identification［J］．SAE International Journal of Materials and Manufacturing，2011，4（1）：527－534．

［13］ Wilkinson，Peter M，Van Dierendonck，et al．Theoretical model for the influence of gas properties and pressure on single-bubble formation at an orifice［J］．Source：Chemical Engineering Science，1994，9（49）：1429－1438．

［14］沈淑君，劉發旺．解分數階Bagley-Torvik方程的一種計算有效的數值方法［J］．廈門大學學報（自然科學版），2004，7（3）：306－311．SHEN Shu-jun， LIU Fa-wang． A computational efficient numerical method of solving fractional Bagley－Torvik equation［J］．Journal of Xiamen University（Natural Science）2004，7（3）：306－311．

［15］薛定宇．高等應用數學問題的MATLAB求解［M］．北京：清華大學出版社，2008：398－399．

［16］劉鯖潔，陳桂明，楊旗．基于Matlab工具的遺傳算法求解有約束最優化問題［J］．兵工自動化，2008，27（11）：43－44．LIU Qing-jie， CHEN Gui-ming， YANG Qi． Solving constrained optimization through genetic algorithm based on matlab toolbox ［J］． Advanced Manufacture and Management，2008，27（11）：43－44．