錢 承 胡紅生
嘉興學(xué)院機(jī)電工程學(xué)院,嘉興,314001
壓電作動器以壓電智能材料為主要材料,具有結(jié)構(gòu)緊湊、作動力大、剛性高、位移分辨率高、頻率響應(yīng)快、控制驅(qū)動簡單等優(yōu)點(diǎn),目前壓電微位移器件已得到廣泛的研究和應(yīng)用[1-3]。但壓電材料同時也存在著固有缺陷,如壓電陶瓷材料的遲滯特性、輸入輸出的非線性特性、蠕變特性等[4-5]。這些固有特性的存在使得壓電作動器在實(shí)際應(yīng)用中重復(fù)性和控制精度降低,實(shí)時瞬態(tài)響應(yīng)也相應(yīng)變慢,這會使得壓電作動器在工業(yè)應(yīng)用中的推廣具有一定難度。特別是壓電材料的遲滯特性,該特性使得控制變得十分困難。
目前國內(nèi)外學(xué)者已對壓電材料的遲滯特性開展了廣泛研究,取得了一定的成果并進(jìn)行了相關(guān)的應(yīng)用,主要涉及壓電材料遲滯模型的理論研究或采用智能算法來建立壓電非線性仿真模型以指導(dǎo)應(yīng)用。尤其在遲滯模型方面,截至目前已形成了多種常用模型,分別為Preisach模型[6-7]、Maxwell模型[8]、Dahl模型[9]、Bouc-Wen模型[10-11]等,其中,Preisach模型及其改進(jìn)形式已被廣泛應(yīng)用。文獻(xiàn)[12]基于徑向基神經(jīng)網(wǎng)絡(luò)法來構(gòu)建遲滯非線性模型,并進(jìn)行了補(bǔ)償控制,得到了很好的跟蹤控制性能。文獻(xiàn)[13]提出了基于模糊控制系統(tǒng)的建模方法,并通過所建立的逆模型進(jìn)行前饋遲滯補(bǔ)償控制,在頻率50 Hz和100 Hz下明顯減小了遲滯特性,且該建模方法簡單,可應(yīng)用于實(shí)時在線建模。文獻(xiàn)[14]為取得壓電作動器的跟蹤精度,在一階回轉(zhuǎn)曲線的基礎(chǔ)上采用了雙輸入的Preisach遲滯模型來預(yù)測耦合遲滯特性,通過實(shí)驗(yàn)對比單輸入的Preisach遲滯模型的控制情況,結(jié)果表明雙輸入的Preisach遲滯模型具有更好的性能。文獻(xiàn)[15]將壓電作動器應(yīng)用于顯微操縱器,采用Preisach遲滯理論建立壓電作動器的遲滯模型,并提出利用前饋非線性PID控制法結(jié)合Preisach遲滯補(bǔ)償法來控制定位精度,最終通過實(shí)驗(yàn)驗(yàn)證了該機(jī)構(gòu)控制性能滿足設(shè)計要求。
對于Preisach遲滯理論的建模,其計算的難點(diǎn)為求取各個遲滯單元的加權(quán)系數(shù)。通過一系列的一階回轉(zhuǎn)曲線來求取遲滯單元加權(quán)系數(shù)的方法被證實(shí)為一種求解精度較高的方法,但該方法的缺點(diǎn)是需要通過實(shí)驗(yàn)獲取大量的一階回轉(zhuǎn)曲線,要使所有的實(shí)驗(yàn)均達(dá)到很高的精度,則具有一定的難度。針對上述問題,研究人員采用人工智能技術(shù)來計算遲滯單元加權(quán)系數(shù)[12,16-18],通過測試一定量的一階回轉(zhuǎn)曲線,根據(jù)一階回轉(zhuǎn)曲線建立人工智能模型從而建立更詳細(xì)的遲滯模型,這將減少大量測試帶來的時間及誤差,但辨識精度、效率仍依賴于實(shí)驗(yàn)樣本,且模型參數(shù)隨外部輸入信號變化的自適應(yīng)性較差,擬合的數(shù)學(xué)曲線在實(shí)際工程應(yīng)用中無法精確模擬壓電遲滯特性,故誤差難以控制。
最小二乘支持向量機(jī)(least squares support vector machines,LS-SVM)將傳統(tǒng)支持向量機(jī)(SVM)中的不等式約束改成等式約束,是SVM的一種改進(jìn)形式。LS-SVM將二次規(guī)劃問題轉(zhuǎn)化為線性方程組求解問題,提高了求解速度及收斂精度[19-20]。考慮到LS-SVM的諸多優(yōu)點(diǎn),本文基于Preisach離散遲滯模型,在實(shí)驗(yàn)的基礎(chǔ)上得到壓電作動器的一階回轉(zhuǎn)曲線,通過一階回轉(zhuǎn)曲線求取遲滯單元加權(quán)值并應(yīng)用LS-SVM法建立所需的遲滯模型,引入遺傳算法對模型參數(shù)進(jìn)行尋優(yōu)得到最精確的模型,最后通過仿真對比分析,驗(yàn)證了本文方法的可行性。
典型的Preisach遲滯模型[21-22]通過對遲滯因子的雙重積分來計算模型的輸出,其計算表達(dá)式如下:
f(t)=?α≤βμ(α,β)γα,β(u(t))dαdβ
(1)
式中,α、β分別為遲滯單元的下閾值和上閾值;μ(α,β)為遲滯單元的加權(quán)函數(shù);u(t)為遲滯模型的輸入;γα, β(u(t))為遲滯單元的值。
遲滯單元計算原理如圖1所示。對于圖1所示的遲滯單元模型,也可將其表述在Preisach平面內(nèi),該平面由α、β為橫縱坐標(biāo)且β≥α的三角形區(qū)域構(gòu)成,如圖2所示。
圖1 遲滯單元計算原理Fig.1 Calculation principle of hysteresis unit
圖2 遲滯計算模型的α-β平面Fig.2 α-β plane of hysteresis calculated model
由圖2可以看出,在三角形區(qū)域內(nèi),加權(quán)函數(shù)μ(α,β)不為零。當(dāng)輸入u(t)增大或減小時會激活不同區(qū)域的遲滯單元或使激活的單元變?yōu)榉羌せ顮顟B(tài),這就需要計算各個激活單元的加權(quán)函數(shù)μ(α,β),然后根據(jù)這些加權(quán)函數(shù)求取最終的輸出。常用的方法為基于一階回轉(zhuǎn)曲線來確定加權(quán)函數(shù),且一階回轉(zhuǎn)曲線越多,模型的準(zhǔn)確度就越高。該方法建立的遲滯模型為有限個遲滯單元的并聯(lián)形式,因此其離散遲滯模型的表達(dá)式如下:
(2)
其中,μn為n×n的矩陣,矩陣元素為分割后的壓電遲滯單元加權(quán)量,可通過一階回轉(zhuǎn)曲線實(shí)驗(yàn)獲取的數(shù)據(jù)得到;γi,j(t)為判斷矩陣,可判斷分割后的壓電遲滯單元是否激活,矩陣中元素是隨時間變化的。
相應(yīng)的離散遲滯模型如圖3所示。
圖3 離散遲滯模型并聯(lián)圖Fig.3 Parallel figure of discrete hysteresis model
為此,本文的研究內(nèi)容主要圍繞如何利用有限的一階回轉(zhuǎn)曲線快速有效地確定離散模型單個遲滯單元的加權(quán)量。
支持向量機(jī)的核心理論是通過核函數(shù)定義的非線性變換將n維樣本空間(x1,y1),(x2,y2), …(xl,yl)∈Rn,yi∈{+1,-1}(i=1,2,…,n)映射到一個高維特征空間,在此高維空間中尋找輸入量與輸出量之間的一種非線性關(guān)系。假設(shè)非線性映射為φ(xi)(i=1,2,…,n),將n維樣本數(shù)據(jù)(x1,y1),(x2,y2),…,(xl,yl) ∈Rn映射到高維特征空間,根據(jù)結(jié)構(gòu)風(fēng)險最小化原則,求解如下最優(yōu)化問題:
(3)
式中,Re為經(jīng)驗(yàn)風(fēng)險值;c為懲罰因子;εi為樣本誤差,也稱松弛變量;b為閾值;w為超平面的權(quán)值向量。
將式(3)的標(biāo)準(zhǔn)支持向量機(jī)最優(yōu)化問題轉(zhuǎn)化成二次規(guī)劃問題,利用核函數(shù)取代高維空間中的點(diǎn)積運(yùn)算,則可把式(3)轉(zhuǎn)化成LS-SVM最優(yōu)化問題,優(yōu)化目標(biāo)損失函數(shù)為誤差εi的二次項,其關(guān)系表達(dá)式如下:
(4)
為了求解該優(yōu)化問題,將約束條件改為等式約束,引入拉格朗日算子ai(i=1,2,…,n),構(gòu)造拉格朗日方程以求解此優(yōu)化問題,即
(5)
要使目標(biāo)函數(shù)取得最小值,則要使拉格朗日方程中的變量w、εi、b、ai的偏導(dǎo)數(shù)均為0,即
(6)
通過求解該二次規(guī)劃問題,可構(gòu)造函數(shù)如下:
(7)
(8)
式中,K(x,xi)為核函數(shù);σ為核參數(shù)。
基于該理論可以發(fā)現(xiàn),LS-SVM法是處理非線性、小樣本的回歸預(yù)測計算方法。針對壓電遲滯模型建立的難點(diǎn)正好符合非線性、小樣本的特性,即可通過較少的一階回轉(zhuǎn)曲線來求解本文所需的遲滯單元加權(quán)系數(shù)。
為了將上述理論應(yīng)用于壓電遲滯模型的參數(shù)辨識,必須通過實(shí)驗(yàn)的方式獲取壓電作動器一定量的一階回轉(zhuǎn)曲線。壓電作動器選用的型號為P-843.20(PI公司),其主要參數(shù)如下:開環(huán)輸出位移30 μm;閉環(huán)輸出位移30 μm;集成反饋傳感器為電阻應(yīng)變片式傳感器(strain gauge sensor,SGS);開環(huán)/閉環(huán)精度0.6/0.3 μm;靜態(tài)剛度27 N/μm;推力/抗拉能力800/300 N;電容3.0 μF。本文將壓電驅(qū)動電壓范圍平均分成五等分,以20 V為間隔作為一階回轉(zhuǎn)曲線的各回轉(zhuǎn)電壓值(即對20 V、40 V、60 V、80 V、100 V這5個電壓值進(jìn)行回轉(zhuǎn)),如此便將離散遲滯模型式(2)中的n設(shè)為5,則相應(yīng)的遲滯單元有n(n+1)/2個,通過這些回轉(zhuǎn)曲線便能求出各個遲滯單元的加權(quán)量。
為使得到的一階回轉(zhuǎn)曲線數(shù)據(jù)準(zhǔn)確可靠,實(shí)驗(yàn)進(jìn)行了10個循環(huán)測量,電壓加載步長為1 V,為使得到的數(shù)據(jù)是壓電作動器運(yùn)行在穩(wěn)定狀態(tài)下摒棄了前3次循環(huán)的數(shù)據(jù),將后7個循環(huán)的數(shù)據(jù)進(jìn)行平均計算,得到一階回轉(zhuǎn)曲線如圖4所示。
圖4 一階回轉(zhuǎn)曲線Fig.4 First order reversed curves
根據(jù)該一階回轉(zhuǎn)曲線計算本文所需的初始加權(quán)值μ6:
(9)
該初始加權(quán)值是將遲滯模型在α-β平面內(nèi)劃分為15個遲滯單元計算所得,為了消除Preisach遲滯模型固有的遲滯特性,將該初始加權(quán)值分別除以相應(yīng)遲滯單元的面積,作為該遲滯單元的平均加權(quán)值:
(10)
將式(1)轉(zhuǎn)換為平均加權(quán)值和遲滯單元面積的乘積形式,即
(11)
將LS-SVM理論應(yīng)用于遲滯模型,將上述劃分后相應(yīng)遲滯區(qū)域的形心在α-β平面內(nèi)的坐標(biāo)作為模型的輸入,計算所得的平均加權(quán)系數(shù)則為模型的輸出進(jìn)行訓(xùn)練。對于訓(xùn)練模型,選擇核函數(shù)為高斯徑向基核函數(shù),而懲罰因子c和核參數(shù)σ的選取對訓(xùn)練模型準(zhǔn)確率具有較大的影響,因此必須在模型訓(xùn)練的過程中進(jìn)行尋優(yōu)。本文采用遺傳算法對參數(shù)進(jìn)行優(yōu)化,其步驟如下:①選定訓(xùn)練樣本和校驗(yàn)樣本,設(shè)定懲罰因子c和核參數(shù)σ的區(qū)間(0,100)和(0,10),從而產(chǎn)生LS-SVM參數(shù)初始群體;②設(shè)定交叉概率為0.6,變異概率為0.2,群體規(guī)模為30,進(jìn)化代數(shù)為300;③進(jìn)行訓(xùn)練得出優(yōu)化后的懲罰因子c和核參數(shù)σ。
選定遺傳算法的適應(yīng)度函數(shù):
(12)
式中,yi為期望輸出;f(xi)為實(shí)際輸出;k為一很小的正數(shù),其作用是防止分母為零的情況出現(xiàn),此處取值為10-3。
為了最優(yōu)化模型,選定評價指標(biāo),定義誤差函數(shù)為
(13)
初始交叉概率和初始變異概率分別由下式確定:
(14)
(15)
式中,f′為交叉兩個體較大的適應(yīng)度函數(shù)值;f為個體對應(yīng)的適應(yīng)度函數(shù)大小;favg為樣本的平均適應(yīng)度;fmax為樣本個體的最大適應(yīng)度。
仿真建模和計算均采用MATLAB軟件。圖5所示為初始參數(shù)下模型輸出與實(shí)驗(yàn)數(shù)據(jù)的誤差,可以看出,評價指標(biāo)函數(shù)的值較大,不能滿足模型計算要求。圖6所示為引入遺傳算法后經(jīng)過200步迭代后的誤差,可以看出,評價指標(biāo)函數(shù)值大幅度減小。圖7為在模型迭代穩(wěn)定后評價指標(biāo)函數(shù)圖,可以看出,優(yōu)化前后評價指標(biāo)函數(shù)的值大幅度減小。
圖5 初始參數(shù)樣本誤差Fig.5 Error under initial parameters
圖6 迭代200步后樣本誤差Fig.6 Error after 200 iterations
圖7 迭代穩(wěn)定后樣本誤差Fig.7 Error after stable iterations
經(jīng)過計算,最終確定了懲罰因子c=13.6和核參數(shù)σ=1.4。確定了遲滯訓(xùn)練模型的參數(shù),并建立了最終的計算預(yù)測模型。
將α-β平面劃分為更多遲滯單元的組成形式,并將式(2)中的n設(shè)為20,即根據(jù)第一步建立的模型求取210個平均加權(quán)系數(shù),即可通過較少的一階回轉(zhuǎn)曲線求取α-β平面內(nèi)較多的平均加權(quán)系數(shù),并根據(jù)這些平均加權(quán)系數(shù)求取更精確的壓電作動器的輸出。
通過上述兩個步驟的計算,可得到平均加權(quán)系數(shù)的三維系數(shù)圖,見圖8。
圖8 平均加權(quán)系數(shù)計算值Fig.8 Calculated value of average weighted coefficients
通過平均加權(quán)系數(shù)矩陣和被激活的遲滯單元的激活面積,并利用式(11)即可求得當(dāng)前壓電作動器的輸出值。
為驗(yàn)證遲滯模型的準(zhǔn)確性,本文搭建了壓電作動器遲滯實(shí)驗(yàn)平臺,其主要組成部分包括壓電作動器、壓電驅(qū)動模塊、傳感檢測模塊及計算機(jī)監(jiān)測軟件,如圖9所示。其工作原理為:控制器通過驅(qū)動模塊輸出電壓驅(qū)動壓電作動器,再由傳感模塊對傳感器采集信號并進(jìn)行檢測處理,處理結(jié)果通過控制總線提供給主控模塊進(jìn)行計算、分析并將結(jié)果顯示于電腦顯示屏上,如圖10所示。
1.壓電驅(qū)動模塊 2.微動測量臺架 3.離線傳感模塊4.參數(shù)調(diào)整及顯示模塊 5.計算機(jī) 6.光學(xué)隔振平臺 7.壓電作動器圖9 壓電遲滯實(shí)驗(yàn)平臺Fig.9 Experimental platform of piezoelectric hysteresis
圖10 實(shí)驗(yàn)平臺工作原理Fig.10 Principle of the experimental platform
為了獲得壓電作動器的遲滯特性,控制器的控制模式采用開環(huán)控制。根據(jù)圖10的實(shí)驗(yàn)平臺測試原理,系統(tǒng)利用采集的位移信號作為輸入,利用建立的遲滯模型進(jìn)行計算以獲取輸出電壓,從而完成自動加載實(shí)驗(yàn)。
仿真和實(shí)驗(yàn)驗(yàn)證分別以40 V、60 V、80 V、100 V為回轉(zhuǎn)電壓點(diǎn)獲取一階回轉(zhuǎn)曲線,測試電壓加載步長1 V,一階回轉(zhuǎn)曲線對比結(jié)果如圖11所示。并以每個回轉(zhuǎn)曲線的仿真結(jié)果輸出值與實(shí)驗(yàn)結(jié)果輸出值的平均相對誤差作為評判的參考,結(jié)果見表1。
由圖11和表1可以看出,由本文建立的遲滯模型的仿真結(jié)果和實(shí)驗(yàn)所得的一階回轉(zhuǎn)曲線吻合度較高、誤差較小,因此可得出本文所用的遲滯模型的建立方法為一種切實(shí)可行的方法。
(1)詳細(xì)闡述了壓電作動器遲滯模型的建立理論,以實(shí)驗(yàn)的方式獲得最小二乘支持向量機(jī)模型所需的一階回轉(zhuǎn)曲線。
(2)根據(jù)一階回轉(zhuǎn)曲線計算遲滯模型初始平均加權(quán)系數(shù),并確定訓(xùn)練模型的輸入輸出,導(dǎo)入最小二乘支持向量機(jī)模型進(jìn)行訓(xùn)練。將α-β平面劃分為更多遲滯單元的組成形式,確定預(yù)測計算的輸入值并將其導(dǎo)入訓(xùn)練后的模型得到最終的遲滯單元的平均加權(quán)系數(shù),從而得到壓電作動器的遲滯模型。
(3)通過仿真結(jié)果和實(shí)驗(yàn)結(jié)果的對比,得出本文基于最小二乘支持向量機(jī)所建立的模型能精確描述壓電作動器的遲滯非線性特性,為一種有效的方法。本文方法加快了建模的速度,免去了需大量實(shí)驗(yàn)來獲得一階回轉(zhuǎn)曲線從而計算加權(quán)系數(shù)的麻煩。
(a)電壓U=100 V
(b)電壓U=80 V
(c)電壓U=60 V
(d)電壓U=40 V圖11 仿真和實(shí)驗(yàn)結(jié)果對比圖Fig.11 Contrast figure of simulation and experiment results
回轉(zhuǎn)電壓(V)最小相對誤差(%)最大相對誤差(%)1000.191.83800.231.52600.121.93401.562.01
參考文獻(xiàn):
[1] GUO Jiang, CHEE S K, YANO T, et al. Micro-vibration Stage Using Piezo Actuators [J]. Sensors and Actuators A: Physical,2013,194:119-127.
[2] CAI Kunhai, TIAN Yanling, WANG Fujun, et al. Development of a Piezo-driven 3-DOF Stage with T-shape Flexible Hinge Mechanism [J]. Robotics and Computer-integrated Manufacturing,2016,37:125-138.
[3] 王常松, 梁森, 韋利明. 智能微位移主動隔振控制系統(tǒng)的研究[J]. 振動與沖擊,2015,34(13):211-216.
WANG Changsong, LIANG Sen, WEI Liming. A Smart Micro-displacementactive Vibration Isolation System [J]. Journal of Vibration and Shock,2015,34(13):211-216.
[4] BIGGIO M, BUTCHER M, GIUSTINIANI A, et al. Memory Characteristics of Hysteresis and Creep in Multi-layer Piezoelectric Actuators: an Experimental Analysis [J]. Physica B: Condensed Matter.,2014,435:40-43.
[5] 秦海辰, 尹周平. 壓電陶瓷晶體遲滯特性的本構(gòu)關(guān)系研究 [J]. 中國機(jī)械工程,2014,25(15):2059-2064.
QIN Haichen, YIN Zhouping. Research on Hysteresis Constitutive Relation in Piezoceramic Crystals [J]. China Mechanical Engineering,2014,25(15):2059-2064.
[6] FELIX W, SUTOR A, RUPITSCH S J, et al. A Generalized Preisach Approach for Piezoceramic Materials Incorporating Uniaxial Compressive stress[J]. Sensors and Actuators A: Physical,2012,186(10):223-229.
[7] LIU L, TAN K K, CHEN S, et al. SVD-based Preisach Hysteresis Identification and Composite Control of Piezo Actuators [J]. ISA Transactions,2015,51(3):430-438.
[8] CHOI G H, OH H J, CHIO S G. Repetitive Tracking Control of a Coarse-fine Actuator[C]// Proceedings of the 1999 IEEE/ASME International Conference on Advanced Intelligent Mechatronics. Atlanta,1999:335-340.
[9] XU Q S, LI Y M. Dahl Model-based Hysteresis Compensation and Precise Positioning Control of an XY Parallel Micromanipulator with Piezoelectric Actuation[J]. Journal of Dynamic Systems, Measurement, and Control,2010,132(4):1-12.
[10] WEI Z, WANG Daihua. Non-symmetrical Bouc-Wen Model for Piezoelectric Ceramic Actuators[J]. Sensors and Actuators A: Physical,2012,181(7):51-60.
[11] ZAMAN M A, SIKDER U. Bouc-Wen Hysteresis Model Identification Using Modified Firefly Algorithm [J]. Journal of Magnetism and Magnetic Materials,2015,395:229-233.
[12] 范家華, 馬磊, 周攀, 等. 基于徑向基神經(jīng)網(wǎng)絡(luò)的壓電作動器建模與控制[J]. 控制理論與應(yīng)用,2016,33(7):856-862.
FAN Jiahua, MA Lei, ZHOU Pan, et al. Modeling and Control of Piezoelectric Actuator Based on Radial Basis Function Neural Network[J]. Control Theory & Applications,2016,33(7):856-862.
[13] LI Pengzhi, YAN Feng, GE Chuan, et al. A Simple Fuzzy System for Modelling of Both Rate-independent and Rate-dependent Hysteresis in Piezoelectric Actuators [J]. Mechanical Systems and Signal Processing,2013,36(1):182-192.
[14] DONG Yangyang, HU Hong, WANG Hongjun. Identification and Experimental Assessment of Two-input Preisach Model for Coupling Hysteresis in Piezoelectric Stack Actuators [J]. Sensors and Actuators A: Physical,2014,220:92-100.
[15] TANG Hui, LI Yangmin. Feedforward Nonlinear PID Control of a Novel Micromanipulator Using Preisach Hysteresis Compensator [J]. Robotics and Computer-integrated Manufacturing,2015,34:124-132.
[16] DLALA E. Efficient Algorithms for the Inclusion of the Preisach Hysteresis Model in Nonlinear Finite-element Methods [J]. IEEE Transactions on Magnetics,2011,47(2):395-408.
[17] 趙新龍, 譚永紅, 董建萍. 基于擴(kuò)展輸入空間法的壓電執(zhí)行器遲滯特性動態(tài)建模[J]. 機(jī)械工程學(xué)報,2010,46(20):169-174.
ZHAO Xinlong, TAN Yonghong, DONG Jianping. Dynamic Modeling of Rate-dependent Hysteresis in Piezoelectric Actuators Based on Expanded Input Space Method [J]. Journal of Mechanical Engineering,2010,46(20):169-174.
[18] MA Yingkun, ZHANG Xinong, XU Minglong, et al. Hybrid Model Based on Preisach and Support Vector Machine for Novel Dual-stack Piezoelectric Actuator [J]. Mechanical Systems and Signal Processing,2013,34(1/2):156-172.
[19] 張新鋒, 趙彥. 基于最小二乘支持向量機(jī)的小樣本威布爾可靠性分析[J]. 中國機(jī)械工程,2012,23(16):1967-1971.
ZHANG Xinfeng, ZHAO Yan. Weibull Reliability Analysis in Small Samples Based on LS-SVM[J]. China Mechanical Engineering,2012,23(16):1967-1971.
[20] CHEN T T, LEE S J. A Weighted LS-SVM Based Learning System for Time Series Forecasting [J]. Information Sciences,2015,299:99-116.
[21] KRASNOSEL’SKII M, POKROVSKII A. Systems with Hysteresis [M]. New York: Springer,1989.
[22] SMITH R. Smart Material Systems: Model Development [M]. Philadelphia: Society for Industrial and Applied Mathematics,2005.