譚程 梁志珊
(中國石油大學(北京),地球物理與信息工程學院,北京 102249)
雖然分數(shù)階微積分理論擁有與整數(shù)階微積分幾乎一樣長的研究歷史,但由于缺乏明顯的幾何意義,使其應(yīng)用一時受到了限制.直到近年來,在機械、物理、工程、信息科學、材料科學等領(lǐng)域發(fā)現(xiàn)存在分數(shù)階現(xiàn)象,使得分數(shù)階微積分理論有了實際的應(yīng)用背景,從而成為了物理學和工程學的研究熱點[1?3].整數(shù)階微積分是相應(yīng)分數(shù)階微積分的特例情況[4].已有的研究表明,相比于整數(shù)階的系統(tǒng)模型,其分數(shù)階模型能更透徹、更準確的反映系統(tǒng)的物理現(xiàn)象[5].
近年來,對電感和電容數(shù)學建模的研究結(jié)果表明:電感和電容本質(zhì)上都是分數(shù)階的[6],整數(shù)階的電感和電容在實際中并不存在,基于分數(shù)階微積分理論建立的電感和電容模型更能反映其電特性[7,8],以往用來描述電感和電容電特性的整數(shù)階模型是不夠準確的甚至可能是錯誤的.然而,電感和電容又是開關(guān)功率變換器電路中不可或缺的電子器件.以往對開關(guān)功率變換器的模型研究都是建立在電感和電容是整數(shù)階基礎(chǔ)上的,顯然這與其分數(shù)階的本質(zhì)是相違背的,是不科學的,這不能準確的反映開關(guān)功率變換器的動力學特性甚至可能會得出錯誤的結(jié)論.在已有的研究中,文獻[9]分析了分數(shù)階電容對功率因數(shù)校正變換器的影響,卻沒有考慮電感也是分數(shù)階的.文獻[10]所建立的分數(shù)階Buck-Boost變換器的模型也僅考慮了電容是分數(shù)階的.鑒于電感和電容本質(zhì)上都是分數(shù)階的事實,文獻[11]對工作于電感電流連續(xù)模式(continuous conduction mode,CCM)下的Boost變換器進行了分數(shù)階區(qū)間數(shù)學模型和分數(shù)階狀態(tài)平均模型的建立和分析,但沒有進行電路模型的仿真驗證且分數(shù)階CCM Boost變換器的控制輸出傳遞函數(shù)具有右半平面零點(right half plane,RHP)問題;文獻[12]對工作于電感電流斷續(xù)模式(discontinuous conduction mode,DCM)下的Boost變換器進行了分數(shù)階區(qū)間數(shù)學模型和分數(shù)階狀態(tài)平均模型的建立和分析,但也沒有進行電路模型的仿真驗證且分數(shù)階DCM Boost變換器有電感電流紋波大、帶負載能力弱等問題.而工作于電感電流偽連續(xù)模式(pseudo continuous conduction mode,PCCM)下的Boost變換器(又稱三態(tài)Boost變換器)具有比CCM模式和DCM模式變換器更優(yōu)良的工作性能[13].PCCM是一種介于CCM和DCM之間,Boost變換器的第三種工作模式[14].相對于CCM Boost變換器,PCCM Boost變換器的控制輸出傳遞函數(shù)不存在RHP問題,提高了系統(tǒng)的閉環(huán)穩(wěn)定性和動態(tài)響應(yīng)性能;相對于DCM Boost變換器,PCCM Boost變換器具有電感電流紋波小、帶負載能力強的優(yōu)點.因此,對PCCM Boost變換器的分數(shù)階模型的研究是一項具有重要理論意義和實際應(yīng)用價值的課題.
本文以PCCM Boost變換器為研究對象,推導出了PCCM Boost變換器的分數(shù)階區(qū)間數(shù)學模型和分數(shù)階狀態(tài)平均模型,對其電感電流和輸出電壓進行了理論分析以及傳遞函數(shù)的推導,并基于Matlab/Simulink的仿真環(huán)境,通過對推導的數(shù)學模型和電路模型進行仿真,分析了模型誤差產(chǎn)生的原因,驗證了分數(shù)階建模與理論分析的正確性.最后,指出了分數(shù)階Boost變換器工作于電感電流偽連續(xù)模式與連續(xù)模式、斷續(xù)模式的區(qū)別與聯(lián)系.
根據(jù)文獻[8]可知,分數(shù)階電感和分數(shù)階電容的數(shù)學模型為
其中,iL為流過分數(shù)階電感的電流,vL為分數(shù)階電感兩端的電壓,L為分數(shù)階電感值,iC為流經(jīng)分數(shù)階電容的電流,vo為分數(shù)階電容兩端的電壓,C為分數(shù)階電容值,α和β分別為電感和電容的分數(shù)階階數(shù)且0<α,β<1.
Boost變換器工作于電感電流偽連續(xù)模式下的電路原理圖和時序脈沖及電感電流波形,如圖1所示.其中,vin為輸入電壓,R為負載電阻,PS1為控制開關(guān)管S1通斷的周期性脈沖信號,PS2為控制開關(guān)管S2通斷的周期性脈沖信號,周期均為T.則PCCM Boost變換器的工作原理為
1)工作模態(tài)1(0<t<d1T):周期性脈沖信號PS1為高電平,PS2為低電平時,開關(guān)管S1導通、S2關(guān)斷,二極管Di1承受反向電壓而關(guān)斷,持續(xù)時間為d1T.
2)工作模態(tài)2(d1T<t<(d1+d2)T):周期性脈沖信號PS1為低電平,PS2為低電平時,開關(guān)管S1關(guān)斷、S2關(guān)斷,二極管Di1承受正向電壓而導通,持續(xù)時間為d2T.
3)工作模態(tài)3((d1+d2)T<t<T):周期性脈沖信號PS1為低電平,PS2為高電平時,開關(guān)管S1關(guān)斷、S2導通,二極管Di1承受反向電壓而關(guān)斷,持續(xù)時間為d3T,其中d1+d2+d3=1.
圖1 PCCM Boost變換器 (a)電路原理圖;(b)時序脈沖及電感電流波形圖
根據(jù)狀態(tài)平均法、PCCM Boost變換器三個工作模態(tài)的特點及分數(shù)階微積分的性質(zhì)[15],對(2),(3),(4)式在一個工作周期T內(nèi)求平均,可推導得到工作于電感電流偽連續(xù)模式下的狀態(tài)平均模型為
其中,〈vin〉,〈vo〉,〈iL〉分別為vin,vo,iL在一個周期內(nèi)的平均值,令Vin,Vo,IL,D1,D2分別為vin,vo,iL,d1,d2的直流分量,?vin,?vo,?iL,?d1,?d2分別為vin,vo,iL,d1,d2的交流分量.因此,可對?vin,?vo,?iL,?d1,?d2作如下分解:
且
將(8)式中的直流分量分離出來,可得
根據(jù)Caputo分數(shù)階導數(shù)定義可知[15],常數(shù)的任意分數(shù)階導數(shù)等于零,則由(9)式可得系統(tǒng)處于穩(wěn)態(tài)時的工作點為
此外,根據(jù)Caputo分數(shù)階導數(shù)定義,由(2)式可求得電感電流iL在(0,d1T)時間內(nèi)的增加量,即電感電流紋波?iL為
將(6)式代入(11)式且忽略高階小量,并將其直流分量分離后,可得
其中,Γ(·)為伽馬函數(shù)[15].可見,電感電流紋波?iL不僅與電感L、輸入電壓直流分量Vin、占空比直流分量D1以及開關(guān)周期T有關(guān),而且還與電感的階數(shù)α有關(guān),并且與電感的階數(shù)α成反比例關(guān)系.與此同時,當α=1時,(12)式與用整數(shù)階模型所求的結(jié)果一致.
根據(jù)(10)式和(12)式可求得電感電流峰值iL?max和谷值iL?min的表達式分別為
由于PCCM Boost變換器獨特的工作模式,使其在工作模式1和工作模式3時,輸出電壓vo均處于下降過程,因此可知輸出電壓紋波?vo為
其中,?v1為處于工作模式1時的電壓減少量,?v2為處于工作模式3時的電壓減少量.
同理,根據(jù)(2)式和(4)式可求得輸出電壓vo在(0,d1T)和((d1+d2)T,T)時間內(nèi)的減少量,即輸出電壓紋波?vo為
其中,Eβ(·)為Mittag-Leラe函數(shù)[15],Vo?max為輸出電壓峰值,其表達式為
將(10)式和(16)式代入(17)式可得
將(18)式代入(16)式可得電壓輸出紋波表達式為
根據(jù)電路模型和數(shù)學模型的電感電流iL和輸出電壓vo,可得所建數(shù)學模型的誤差百分比為
其中,bfbi為電感電流模型誤差百分比,bfbv為輸出電壓模型誤差百分比,iLd和vod分別為電路模型的電感電流和輸出電壓,iLs和vos分別為數(shù)學模型的電感電流和輸出電壓.
可見,輸出電壓紋波?vo不僅與電容C、輸入電壓直流分量Vin、占空比直流分量D1,D2,D3、負載電阻R以及開關(guān)周期T有關(guān),而且還與電容的階數(shù)β有關(guān),并且與電感的階數(shù)β成反比例關(guān)系.與此同時,當β=1時,(19)式與用整數(shù)階模型所求的結(jié)果一致.
將(8)式的交流分量分離出來,可得
(22)式表明當占空比d1,d2的擾動變量==0時,輸入電壓vin的變化對輸出電壓vo的影響.
(23)式表明當輸入電壓vin的擾動變量=0,占空比d2的擾動變量=0時,占空比d1的變化對輸出電壓vo的影響.
輸出電壓?vo對占空比的傳遞函數(shù)Gvd2(s)為
(24)式表明當輸入電壓vin的擾動變量=0,占空比d1的擾動變量=0時,占空比d2的變化對輸出電壓vo的影響.
(25)式表明當占空比d1,d2的擾動變量==0時,輸入電壓vin的變化對電感電流iL的影響.
(26)式表明當輸入電壓vin的擾動變量=0,占空比d2的擾動變量=0時,占空比d1的變化對電感電流iL的影響.
輸出電流?iL對占空比的傳遞函數(shù)Gid2(s)為
(27)式表明當輸入電壓vin的擾動變量=0,占空比d1的擾動變量=0時,占空比d2的變化對電感電流iL的影響.
根據(jù)(23)式可知,分數(shù)階PCCM Boost變換器的控制輸出傳遞函數(shù)不存在RHP問題.當電感和電容的分數(shù)階階數(shù)α,β都等于1時,(22)–(27)式所示的分數(shù)階傳遞函數(shù)與文獻[14,16]所描述的整數(shù)階傳遞函數(shù)相一致,進一步說明了整數(shù)階系統(tǒng)是其分數(shù)階系統(tǒng)的特例情況,且在頻域和時域響應(yīng)中,α,β這兩個參數(shù)對PCCM Boost變換器系統(tǒng)的動力學特性也產(chǎn)生了極大的影響.
根據(jù)文獻[17],基于分抗鏈[15]和改進的Oustaloup濾波器的分數(shù)階微積分算法[18],可以得到分數(shù)階電感和分數(shù)階電容的等效電路模型,如圖2和圖3所示.當Lα=3 mH,α=0.8時,圖2中各電阻值分別為RL1=7.16 k?,RL2=340.84 ?,RL3=34.25 ?,RL4=3.54 ?,RL5=367 m?,RL6=38 m?,RL7=4 m?,RL8=0.4 m?,RL9=42μ?,RL10=5μ?;各電感值分別為L1=95μH,L2=77μH,L3=131.6μH,L4=231.6 μH,L5=408 μH,L6=719.4 μH,L7=1.268 mH,L8=2.235 mH,L9=3.934 mH.當Cβ=100μH,β=0.8時,圖3中各電阻值分別為RC1=20 m?,RC2=160 m?,RC3=1.5 ?,RC4=14.6 ?,RC5=141 ?,RC6=1.36 k?,RC7=13.131 k?,RC8=126.742 k?,RC9=1.222 M?,RC10=102.85 M?;各電容值分別為C1=6.5μF,C2=13.98μF,C3=24.5μF,C4=43.2 μF,C5=76.2 μF,C6=134.2 μF,C7=236.6 μF,C8=417 μF,C9=736 μF,C10=560μF.
圖2 分數(shù)階電感等效電路模型
圖3 分數(shù)階電容等效電路模型
文獻[11]和文獻[12]雖然分別建立了分數(shù)階CCM Boost變換器的數(shù)學模型和分數(shù)階DCM Boost變換器的數(shù)學模型,并進行了仿真,但由于PCCM Boost變換器獨特的工作模式,使得文獻[11]和文獻[12]所建立的仿真數(shù)學模型并不能適用于PCCM Boost變換器的數(shù)學模型仿真.因此,必須根據(jù)PCCM Boost變換器自身的特點,并根據(jù)文獻[18]所提出的改進的Oustaloup濾波器的分數(shù)階微積分算法和(5)式重新構(gòu)建Matlab/Simulink數(shù)學模型,如圖4所示.其中,Fractional Ints?α為分數(shù)階積分單元,其內(nèi)部結(jié)構(gòu)如文獻[11]所示.在改進的Oustaloup濾波器的分數(shù)階微積分算法中,存在三個關(guān)鍵參數(shù):擬合頻率下限ωb、擬合頻率上限ωh、濾波器階數(shù)2N+1.而在對實際分數(shù)階系統(tǒng)進行數(shù)值仿真時,需根據(jù)系統(tǒng)的頻率范圍選擇擬合頻率段(ωb,ωh)和N 值,一般取ωb·ωh=1. 選取電路參數(shù)為vin=24 V,L=3 mH,C=100μF,d1=0.4,d2=0.2,f=50 kHz.由于開關(guān)頻率f=50 kHz,即ω =2πf=3.14× 105rad/s,考慮還有高于開關(guān)頻率的高頻諧波存在,因而擬合頻率的選取需滿足條件ωh>3.14×105rad/s.因此,選取ωh=1×106rad/s,ωb=1×10?6rad/s,N=10.根據(jù)分數(shù)階電感和電容的等效電路模型所建立的電路仿真模型,如圖5所示.
圖4 分數(shù)階PCCM Boost變換器Simulink數(shù)學仿真模型
圖5 分數(shù)階PCCM Boost變換器Simulink電路仿真模型
圖6 分數(shù)階PCCM Boost變換器電路模型仿真波形(a)電感電流iL;(b)輸出電壓vo
圖7 分數(shù)階PCCM Boost變換器電路模型與數(shù)學模型仿真波形比較 (a)電感電流iL;(b)輸出電壓vo
當取α=0.8,β=0.8時,根據(jù)文獻[11],可求得分數(shù)階Boost變換器工作于臨界狀態(tài)時的負載電阻值為R=1671.687 ?,則當R < 1671.687 ?時,可保證系統(tǒng)工作于偽連續(xù)模式.因此,選取R=50 ?.當PCCM Boost變換器處于穩(wěn)定運行狀態(tài)時,其電路模型的電感電流iL和輸出電壓vo的波形分別如圖6(a)和(b)所示.顯然,此種情況下Boost變換器工作于電感電流偽連續(xù)模式.對其電路模型的電感電流iL和輸出電壓vo在一個開關(guān)周期T內(nèi)進行平均,并與其數(shù)學模型的電感電流iL和輸出電壓vo進行比較,如圖7(a)和(b)所示.由圖7結(jié)果可知,對PCCM Boost變換器所建的分數(shù)階數(shù)學模型是正確的.根據(jù)圖6(a)可測量出?iL=0.842 A,iL?min=6.800 A,iL?max=7.642 A,IL=7.211 A.根據(jù)圖6(b)可測量出?vo=3.200 V,vo?min=68.610 V,vo?max=71.810 V,Vo=70.210 V. 根據(jù)(10)式、(12)式、(13)式和(14)式可分別計算出?iL=0.720 A,iL?min=6.840 A,iL?max=7.560 A,IL=7.200 A;根據(jù)(10)式、(18)式和(19)式可分別計算出?vo=2.288 V,vo?min=70.856 V,vo?max=73.144 V,Vo=72.000 V.可知,對數(shù)學模型的理論分析和電路模型的仿真結(jié)果基本一致,從而表明對工作于偽連續(xù)模式下Boost變換器理論分析的正確性.
圖8 整數(shù)階PCCM Boost變換器Simulink數(shù)學仿真模型
圖9 整數(shù)階PCCM Boost變換器Simulink電路仿真模型
當把圖4中的分數(shù)階積分單元換為整數(shù)階積分單元時,即采用了整數(shù)階數(shù)學模型來描述PCCM Boost變換器,則其整數(shù)階Matlab/Simulink數(shù)學仿真模型如圖8所示;當用整數(shù)階電感和電容取代圖5中的分數(shù)階電感和電容時,則其整數(shù)階電路模型如圖9所示,其電路模型的電感電流iL和輸出電壓vo的波形分別如圖10(a)和(b)所示.可見,系統(tǒng)工作于電感電流偽連續(xù)模式.對其電路模型的電感電流iL和輸出電壓vo在一個開關(guān)周期T內(nèi)進行平均,并和其數(shù)學模型的電感電流iL和輸出電壓vo進行比較,如圖11(a)和(b)所示.根據(jù)圖10(a)可測量出?iL=0.063 A,iL?min=7.092 A,iL?max=7.155 A,IL=7.124 A;根據(jù)圖10(b)可測量出?vo=0.230 V,vo?min=70.042 V,vo?max=70.650 V,Vo=70.540 V. 根據(jù)(10)式、(12)式、(13)式和(14)式可分別計算出?iL=0.064 A,iL?min=7.168 A,iL?max=7.232 A,IL=7.200 A;根據(jù)(10)式、(18)式和(19)式可分別計算出?vo=0.240 V,vo?min=71.880 V,vo?max=72.120 V,Vo=72.000 V.
圖10 整數(shù)階PCCM Boost變換器電路模型仿真波形(a)電感電流iL;(b)輸出電壓vo
對比整數(shù)階α=1,β=1和分數(shù)階α=0.8,β=0.8PCCM Boost變換器的數(shù)值仿真結(jié)果,可知電感電流直流分量IL和輸出電壓直流分量Vo沒有發(fā)生變化;而電感電流紋波?iL、電感電流峰值iL?max和谷值iL?min、輸出電壓紋波?vo、輸出電壓的峰值vo?max和谷值vo?min以及動態(tài)響應(yīng)過程中的上升時間、延遲時間、調(diào)節(jié)時間、峰值時間、超調(diào)量都發(fā)生了很大的變化.這就表明用整數(shù)階模型描述本應(yīng)該用分數(shù)階模型描述的PCCM Boost變換器,將會在電感電流紋波?iL、電感電流峰值iL?max和谷值iL?min、輸出電壓紋波?vo、輸出電壓的峰值vo?max和谷值vo?min以及動態(tài)響應(yīng)過程中的上升時間、延遲時間、調(diào)節(jié)時間、峰值時間、超調(diào)量等方面得到錯誤的結(jié)果.因此,基于實際電感和電容本質(zhì)上都是分數(shù)階的事實,為了能夠更好的描述工作于電感電流偽連續(xù)模式下Boost變換器的動力學特性,須采用其分數(shù)階形式的數(shù)學模型.
圖11 整數(shù)階PCCM Boost變換器電路模型與數(shù)學模型仿真波形比較 (a)電感電流iL;(b)輸出電壓vo
圖12 模型誤差百分比 (a)電感電流iL;(b)輸出電壓vo
根據(jù)(20)式可得模型誤差百分比曲線,如圖12所示.由文獻[19]可知,狀態(tài)空間平均法僅是平均法的一階近似,所建的狀態(tài)空間模型僅能近似的表示電路模型.因此,圖12所示的模型誤差百分比是合理的,所建立的PCCM Boost變換器的分數(shù)階數(shù)學模型是正確的.根據(jù)圖12可知,模型誤差與電感階數(shù)α和電容階數(shù)β有關(guān).因此,須采用分數(shù)階數(shù)學模型來描述PCCM Boost變換器的動力學特性.
本文基于分數(shù)階微積分理論,建立了PCCM Boost變換器的分數(shù)階數(shù)學模型并進行了相應(yīng)的理論分析;通過對其數(shù)學模型和電路模型的仿真對比,可得出如下結(jié)論
1)所建的PCCM Boost變換器的分數(shù)階數(shù)學模型可以準確的描述其電路模型.
2)對于分數(shù)階PCCM Boost變換器,所建分數(shù)階數(shù)學模型的誤差是由于狀態(tài)平均法只是平均法的近似而產(chǎn)生的,且其模型誤差百分比與電感L的分數(shù)階階數(shù)α和電容C的分數(shù)階階數(shù)β有關(guān),即在其他參數(shù)不變的情況下,隨著電感L的分數(shù)階階數(shù)α和電容C的分數(shù)階階數(shù)β的增大而減小.
3)分數(shù)階PCCM Boost變換器的數(shù)學模型形式上雖然和分數(shù)階DCM Boost變換器一樣,但由于PCCM Boost變換器不同于DCM Boost變換器的工作特點,使得DCM Boost變換器分數(shù)階模型得出的結(jié)論并不能直接應(yīng)用于PCCM Boost變換器的分數(shù)階模型理論分析.
4)分數(shù)階PCCM Boost變換器的數(shù)學模型形式上雖然和分數(shù)階CCM Boost變換器不一樣,但通過理論分析及仿真驗證,可知CCM Boost變換器分數(shù)階模型得出的結(jié)論能夠直接應(yīng)用于PCCM Boost變換器的分數(shù)階模型理論分析.
5)對于分數(shù)階PCCM Boost變換器,在其他參數(shù)不變的情況下,其動態(tài)響應(yīng)過程隨著電感L的分數(shù)階階數(shù)α和電容C的分數(shù)階階數(shù)β的增大而增大,即其階躍響應(yīng)的上升時間、延遲時間、調(diào)節(jié)時間、峰值時間、超調(diào)量都將增大.
綜上所述,基于電感和電容本質(zhì)上是分數(shù)階的事實,本文所建的PCCM Boost變換器的分數(shù)階數(shù)學模型是正確的,能夠真實的反映PCCM Boost變換器的動力學特性.
[1]Yang S P,Zhang R X 2008 Acta Phys.Sin.57 6837(in Chinese)[楊世平,張若洵 2008物理學報 57 6837]
[2]Zhang C F,Gao J F,Xu L 2007 Acta Phys.Sin.56 5124(in Chinese)[張成芬,高金峰,徐磊2007物理學報56 5124]
[3]Li C L,Yu S M,Luo X S 2012 Chin.Phys.B 21 172
[4]Kenneth S M,Bertram R 1993 An Introduction to the Fractional Calculus and Fractiona Differential Equations(New Jersey:John Wiley&Sons)p21
[5]Shockooh A,Suarez L 1999 Journal of Viberation and Control.5 331
[6]Bohannan G W 2002 Proceedings of the 41st IEEE International Conference on Decision and Control,Tutorial Workshop 2:Fractional Calculus Applications in Automatic Control and Robotics Las Vegas,USA,December 10–13,2002 p1
[7]Westerlund S,Ekstam L 1994 IEEE Trans.Dielectr.Electr.Insulat.1 826
[8]Westerlund S 2002 Dead matter has memory(Kalmar,Sweden:Causal Consulting)chapt.7
[9]Ahmad W 2003 Proceedings of the 2003 International Symposium on Circuits and Systems Bangkok,Thailand,May 25–28,2003 3 p5
[10]Martinez R,Bolea Y,Grau A,Martinez H 2009 IEEE Conference on Emerging Technologies&Factory Automation Palma de Mallorca,Spain,September 22–25,2009 p1
[11]Wang F Q,Ma X K 2011 Acta Phys.Sin.60 070506(in Chinese)[王發(fā)強,馬西奎 2011物理學報 60 070506]
[12]Wang F Q,Ma X K 2013 Scientia Sinica Technological.43 368(in Chinese)[王發(fā)強,馬西奎 2013中國科學:43 368]
[13]Ma D S,Ki W H 2007 IEEE Trans.Circuit and Systtems II:Express Briefs.54 825
[14]Kanakasabai V,Ramesh O,Dipti S 2002 IEEE Trans.Power Electronics.17 677
[15]Podlubny I 1999 Fractional differential equations(New York:Academic Press)chapt 1–2,4
[16]Yu H K 2010 M.S.Thesis.(Sichuan:Southeast Jiaotong University)(in Chinese)[于海坤2010碩士學位論文(四川:西南交通大學)]
[17]Wang F Q,Ma X K 2013 Chin.Phys.B 22 236
[18]Xue D Y,Chen Y Q 2007 MATLAB Solutions to Mathematical Problems in Control(Beijing:Tsinghua University Press)p435(in Chinese)[薛定宇,陳陽泉2007控制數(shù)學問題的MATLAB求解(北京:清華大學出版社)第435頁]
[19]Cao W S,Yang Y X 2007 Journal of System Simulation.19 1329(in Chinese)[曹文思,楊育霞2007系統(tǒng)仿真學報19 1329]