圓弧足被動行走機器人動力學(xué)參數(shù)影響研究

張婉婉，高建設(shè)，高家昌，王高峰

（鄭州大學(xué)機械與動力工程學(xué)院，河南 鄭州 450001）

1 引言

由于雙足機器人對復(fù)雜地形具有較強的適應(yīng)性，受到了研究人員廣泛關(guān)注［1-2］。目前，為了實現(xiàn)其穩(wěn)定行走，許多機構(gòu)已經(jīng)開發(fā)了各種具有復(fù)雜控制系統(tǒng)的主動雙足機器人［3-5］。然而，復(fù)雜控制系統(tǒng)的設(shè)計往往會導(dǎo)致機器人不可避免地存在能量損耗大、控制難度高、步態(tài)僵硬等問題。這些問題極大地制約了雙足機器人的發(fā)展和應(yīng)用。

20世紀90年代，文獻［6］提出了被動動力學(xué)理論并實現(xiàn)了二維無膝關(guān)節(jié)和有膝關(guān)節(jié)雙足被動行走機器人在斜坡上的自發(fā)步行。文獻［7］表明被動行走機器人研究不僅有助于理解人類步行的內(nèi)在機理，而且可以指導(dǎo)設(shè)計更高效靈巧的雙足機器人。然而，被動行走機器人的穩(wěn)定行走在很大程度上依賴于初始條件和動力學(xué)參數(shù)［8］。為避免被動行走機器人出現(xiàn)不穩(wěn)定甚至摔倒的情況，動力學(xué)參數(shù)影響研究成為被動行走領(lǐng)域至關(guān)重要的課題。

繼文獻［6］開創(chuàng)性工作之后，許多學(xué)者對雙足被動行走機器人進行了大量相關(guān)性研究。文獻［9］提出了一種最簡模型，研究發(fā)現(xiàn)，當(dāng)斜坡傾角增大時該模型步態(tài)呈現(xiàn)出倍周期分岔現(xiàn)象。文獻［10］提出了一種羅盤模型，除斜坡傾角外，他們還研究了質(zhì)量和腿長對模型步態(tài)的影響。文獻［11］進一步證明，隨著斜坡傾角增大，最簡模型具有周期3 到周期7 步態(tài)。文獻［12］研究表明增大斜坡傾角導(dǎo)致被動雙足機器人穩(wěn)定性下降。文獻［13］詳細分析了羅盤模型的極限環(huán)步態(tài)和能量轉(zhuǎn)化關(guān)系并證明了被動行走過程中機械能守恒。文獻［14］通過計算雅克比矩陣的特征值討論了羅盤模型中機械參數(shù)和斜坡角對穩(wěn)定性的影響。文獻［15］研究了具有點足、圓弧足和平板足的最簡模型，對比分析得出圓弧足可以擴大穩(wěn)定步行范圍而平板足的情況則取決于足形參數(shù)。文獻［16］證明了引入擺臂結(jié)構(gòu)能夠提高被動行走機器人穩(wěn)定性。文獻［17］發(fā)現(xiàn)當(dāng)轉(zhuǎn)動慣量增大、質(zhì)心位置降低或足半徑減小，都會導(dǎo)致圓弧足羅盤模型步態(tài)發(fā)生分岔。文獻［18］證明了具有合理髖關(guān)節(jié)質(zhì)量的圓弧足羅盤模型有更強的行走穩(wěn)定性與魯棒性。

盡管上述幾位學(xué)者針對圓弧足羅盤模型已有了顯著的研究成果，但缺少系統(tǒng)地討論動力學(xué)參數(shù)對圓弧足羅盤模型運動特性的影響，特別是雙參數(shù)變化對機器人穩(wěn)定步態(tài)的影響。另一方面，由于被動行走機器人為非光滑系統(tǒng)，方程中不連續(xù)處雅克比矩陣不可導(dǎo)，采用傳統(tǒng)方法求解李雅普諾夫指數(shù)比較困難。在此首先以圓弧足被動行走機器人為研究對象，對其運動過程進行了分析與建模，著重研究了圓弧足半徑、質(zhì)心位置、轉(zhuǎn)動慣量和斜坡傾角變化時機器人步態(tài)演變規(guī)律，特別是當(dāng)足半徑和質(zhì)心位置同時變化對機器人穩(wěn)定參數(shù)區(qū)間的影響。另外，采用正交擾動向量法求解被動行走機器人的李雅普諾夫指數(shù)，這種方法可以避免直接計算雅可比矩陣。

2 圓弧足被動行走機器人動力學(xué)建模

2.1 模型介紹

首先給出圓弧足被動行走機器人模型，如圖1所示。該模型是由兩個具有完全相同質(zhì)量和長度的剛性桿件在髖關(guān)節(jié)處通過鉸鏈連接而成。由于髖關(guān)節(jié)和圓弧足質(zhì)量可以等效簡化為兩腿質(zhì)量，為簡化模型研究，不考慮髖關(guān)節(jié)和足部質(zhì)量。設(shè)腿質(zhì)量為m，長度為l，腿質(zhì)心與髖關(guān)節(jié)的距離為c，腿相對質(zhì)心的轉(zhuǎn)動慣量為Jc，圓弧足半徑為r且圓心在腿上，斜面傾角為φ，重力加速度為g。

圖1 圓弧足被動行走機器人模型Fig.1 The Model of a Passive Walking Robot with Arc Feet

機器人單步運動過程，如圖2所示。主要是支撐腿和擺動腿之間的轉(zhuǎn)換。機器人運動過程分為兩個階段，第一個階段是擺動階段，第二個階段是碰撞階段。在擺動階段，支撐腿繞其圓弧足做倒立擺運動，擺動腿繞髖關(guān)節(jié)做單擺運動，如圖2（a）、圖2（b）所示。在碰撞階段，擺動腿與斜面發(fā)生碰撞，為清楚了解兩腿角色變換，又將碰撞時刻分為碰撞前和碰撞后，如圖2（c）、圖2（d）所示。碰撞發(fā)生后擺動腿和支撐腿角色互換。

圖2 單步運動過程Fig.2 The Process of Motion in One Step

為了便于進行動力學(xué)建模，參考前人研究［9-10］，對模型做如下假設(shè)：

（1）整個模型為剛體，沒有彈性變形；

（2）髖關(guān)節(jié)與腿之間無阻尼無摩擦；

（3）支撐腿與地面之間為純滾動約束，不產(chǎn)生滑動與變形，且忽略摩擦；

（4）擺動腿與地面碰撞為完全非彈性碰撞且是瞬時發(fā)生。

2.2 擺動階段動力學(xué)方程

該圓弧足被動行走機器人模型為二自由度保守系統(tǒng)，可利用拉格朗日方法進行建模。將系統(tǒng)坐標系建立在斜面上，x軸沿斜面向下，y軸垂直斜面向上，選取支撐腿與斜面垂直時的垂足O為坐標系原點。設(shè)支撐腿與軸的夾角為θs，擺動腿與軸的夾角為θns，取逆時針為正方向，圖2（a）中的θs為正，θns為負。選θs和θns為廣義坐標，則系統(tǒng)微分形式的動力學(xué)方程為：

其中，

2.3 碰撞階段動力學(xué)方程

在碰撞階段，擺動腿足端與斜面發(fā)生碰撞。設(shè)原擺動腿足端與斜面的碰撞點為B，在碰撞前后，擺動腿與斜坡之間的作用力只有沖擊力。故碰撞前后，擺動腿和支撐腿的角度不變，角速度發(fā)生改變。根據(jù)角動量守恒定律，模型整體關(guān)于碰撞點B的角動量守恒，且原支撐腿（即新擺動腿）對于髖關(guān)節(jié)質(zhì)心位置處點的角動量守恒。

設(shè)CH為髖關(guān)節(jié)質(zhì)心位置，Cns為擺動腿質(zhì)心位置，Cs為支撐腿質(zhì)心位置，髖關(guān)節(jié)質(zhì)心速度為v→H，負號和正號的上標分別代表為碰撞前和碰撞后的狀態(tài)符號。碰撞前的坐標系選取為擺動階段建立的坐標系，即以原支撐腿與斜面垂直時的垂足為原點建立坐標系，x軸沿斜面向下，y軸垂直斜面向上。碰撞后，以B點附近新支撐腿與斜面垂直時的垂足為原點建立坐標系，x軸沿斜面向下，y軸垂直斜面向上。

（1）參考點為H時

碰撞前，原支撐腿對參考點H的角動量為：

碰撞后，原支撐腿變?yōu)樾碌臄[動腿，對參考點H的角動量為：

（2）參考點為B時

碰撞前，系統(tǒng)整體對碰撞點B的角動量為：

碰撞后，系統(tǒng)整體對碰撞點B的角動量為：

碰撞前后，擺動腿和支撐腿的角度不變，角速度發(fā)生改變。根據(jù)角動量守恒定律，得L1=L2，L3=L4。由此可得系統(tǒng)碰撞方程為：

2.4 方程無量綱化

為了基于一個貢獻標尺計算不同參數(shù)對機器人動力學(xué)特性的影響，需要對方程進行無量綱化。在此引入無量綱參數(shù)：

其中，

經(jīng)化簡得碰撞方程的無量綱形式為：

其中，

式中：x1—擺動階段支撐腿角位移；x2—擺動階段支撐腿角速度；x3—擺動階段擺動腿角位移；x4—擺動階段擺動腿角速度；y1—碰撞后支撐腿角位移；y2—碰撞后支撐腿角速度；y3—碰撞后擺動腿角位移；y4—碰撞后擺動腿角速度。

取定參數(shù)：

由上面所取定參值可得無量綱常量

β1=0.24，β2=0.5，β3=0.04

3 單參數(shù)變化對機器人步態(tài)的影響

根據(jù)第二節(jié)所述動力學(xué)方程，選擇β1、β2、β3和φ作為目標參數(shù)，它們分別與圓弧足半徑、質(zhì)心位置、轉(zhuǎn)動慣量和斜坡傾角有關(guān)。接下來借助分岔圖進行參數(shù)變化影響研究。

3.1 圓弧足半徑變化的影響

設(shè)置參數(shù)β2=0.5，β3=0.04，φ=0.01，β1從0到0.9變化。機器人支撐腿角位移隨足半徑參數(shù)β1變化的分岔圖，如圖3所示。可以看出，隨著足半徑參數(shù)β1增大，機器人支撐腿角位移未發(fā)生分岔現(xiàn)象，當(dāng)在區(qū)間（0，0.805）變化時，機器人始終具有穩(wěn)定周期一步態(tài)。但當(dāng)足半徑參數(shù)β1在0.805之后增加時，角位移急劇增大并最終發(fā)散，此參數(shù)下的機器人在行走過程中將會極易跌倒。

圖3 β1變化時的分岔圖Fig.3 Bifurcation Diagram When β1 Changes

3.2 質(zhì)心位置變化的影響

設(shè)置參數(shù)β1=0.24，β3=0.04，φ=0.01，β2從（0.5～0.72）變化。機器人支撐腿角位移隨質(zhì)心位置參數(shù)β2變化的分岔圖，如圖4所示。從圖中可以看出，隨著β2增大，機器人步態(tài)出現(xiàn)了倍周期分岔現(xiàn)象。其中，穩(wěn)定的周期一步態(tài)占據(jù)了較大參數(shù)區(qū)間，后在區(qū)間（0.6629，0.7094）內(nèi)表現(xiàn)為周期二步態(tài)且兩個分支存在匯交于一點的現(xiàn)象。當(dāng)β2繼續(xù)增加時，經(jīng)過短暫的周期四步態(tài)后，機器人將失去失穩(wěn)性。

圖4 β2變化時的分岔圖Fig.4 Bifurcation Diagram When β2 Changes

3.3 轉(zhuǎn)動慣量變化的影響

設(shè)置參數(shù)β1=0.24，β2=0.5，φ=0.01，β3從0.04 到0.19 變化。機器人支撐腿角位移隨轉(zhuǎn)動慣量參數(shù)β3變化的分岔圖，如圖5所示。從圖中可以看出，隨著β3增大，機器人步態(tài)同樣出現(xiàn)了倍周期分岔現(xiàn)象。這與質(zhì)心位置增加對機器人步態(tài)的影響相同，只是分岔發(fā)生位置不一致。由圖知，機器人周期步態(tài)分別在0.109和0.184處發(fā)生分岔。以上現(xiàn)象說明轉(zhuǎn)動慣量的增大也會導(dǎo)致機器人周期步態(tài)發(fā)生分岔現(xiàn)象并最終失穩(wěn)。

圖5 β3變化時的分岔圖Fig.5 Bifurcation Diagram When β3 Changes

3.4 斜坡傾角變化的影響

設(shè)置參數(shù)β1=0.24，β2=0.5，β3=0.04，φ從0.01到0.16變化。機器人支撐腿角位移隨斜坡傾角φ變化的分岔圖，如圖6所示。

圖6 φ變化時的分岔圖Fig.6 Bifurcation Diagram When φ Changes

從圖中可以看出，隨著斜坡傾角φ增大，機器人步態(tài)同樣地出現(xiàn)了倍周期分岔現(xiàn)象。與上述倍周期分岔不同的是，此時機器人步態(tài)由倍周期分岔通向混沌。從圖中觀察到，機器人在區(qū)間（0.01，0.138）內(nèi)表現(xiàn)為周期一步態(tài)，區(qū)間（0.138，0.151）內(nèi)為周期二步態(tài)。隨著斜坡傾角φ繼續(xù)增加，機器人經(jīng)過短暫的周期四步態(tài)后進入混沌狀態(tài)。

4 李雅普諾夫指數(shù)計算

上述分岔圖定性地表現(xiàn)出機器人步態(tài)隨參數(shù)變化的演變規(guī)律，而李雅普諾夫指數(shù)能夠定量地確定系統(tǒng)在特定參數(shù)下的步態(tài)特性。由于碰撞發(fā)生，導(dǎo)致被動行走機器人運動方程不連續(xù)。因此，機器人運動軌跡相對于初始條件的導(dǎo)數(shù)即雅克比矩陣變得病態(tài)或根本無法計算，用傳統(tǒng)方法求解李雅普諾夫指數(shù)比較困難。利用初始正交向量的小擾動估計雅可比矩陣的思想求解被動行走機器人李雅普諾夫指數(shù)，這種方法極大地簡化了非光滑系統(tǒng)的李亞普諾夫指數(shù)求解問題，因為這種估計方法可以避免直接計算雅可比矩陣［19］。估計方法過程如下。根據(jù)第二節(jié)被動行走機器人運動方程為：

式中：x=[x1，x2，x3，x4]T；f=[f1，f2，f3，f4]T；Γ—碰撞面。

非線性系統(tǒng)龐加萊截面的一般形式是：

被動行走系統(tǒng)通常選取碰撞面為龐加萊截面。針對圓弧足被動行走機器人，這里選取龐加萊截面為：

對于等式（6）中的離散映射，直接計算雅可比矩陣是不可能的，因為等式（6）的右邊是未知的。此時引入一個擾動向量Δn=[δ1，δ2，δ3，δ4]T，其中，δi(i=1，2，3)是小量值，可以得到：

其中，

Dui(xn)=雅可比矩陣的列向量；

DU(xn)=[Du1(xn)，Du2(xn)，Du3(xn)，Du4(xn)]T—雅可比矩陣。由等式（7）可知，雅可比矩陣中列的近似值可由下列等式求得：

具體求解步驟如下：

（1）設(shè)初始條件x00=[x10，x20，x30，x40]，將x00代入機器人動力學(xué)方程進行求解得到x10；（2）設(shè)擾動初始條件為x01=[x10+δ，x20，x30，x40]，x02=[x10，x20+δ，x30，x40]，x03=[x10，x20，x30+δ，x40]，x04=[x10，x20，x30，x40+δ]，將x01，x02，x03，x04分別代入動力學(xué)方程進行求解得到x11，x12，x13，x14；（3）求解此時雅克比矩陣的估計值為：[(x11-x10)/δ，(x12-x10)/δ，(x13-x10)/δ，(x14-x10)/δ]T；（4）將x10作為初始條件返回第一步，不斷進行迭代即可求解機器人運動軌跡的雅克比矩陣估計值。根據(jù)構(gòu)造的雅克比矩陣利用傳統(tǒng)方法［19］即可求解李雅普諾夫指數(shù)。與圖6相對應(yīng)的李雅普諾夫指數(shù)圖，如圖7所示。從圖中可以清楚地觀察到兩個分岔點和混沌分布位置。在區(qū)間（0.01，0.138）與（0.138，0.151）內(nèi)，最大李雅普諾夫指數(shù)為負，說明機器人此時處于穩(wěn)定運動狀態(tài)，與分岔圖中周期一和周期二參數(shù)區(qū)間相對應(yīng)；在0.138 和0.151位置處，最大李雅普諾夫指數(shù)為零，此時系統(tǒng)處于不穩(wěn)定狀態(tài)，與分岔圖中的分岔點一致。通過對最大李雅普諾夫指數(shù)進行分析，為機器人分岔動力學(xué)提供了強有力的驗證。

圖7 最大李亞普諾夫指數(shù)圖Fig.7 The Largest Lyapunov Exponent Diagram

5 雙參數(shù)變化對機器人步態(tài)的影響

上幾節(jié)討論了單參數(shù)變化對機器人步態(tài)的影響，但不同參數(shù)聯(lián)合作用下的影響仍需探究。圓弧足半徑參數(shù)β1和質(zhì)心位置參數(shù)β2組合變化下對機器人步態(tài)的影響，如圖8所示。與圖8相對應(yīng)的機器人穩(wěn)定行走參數(shù)區(qū)間，如表1所示。由圖8和表1可得，當(dāng)β1取0時，減小質(zhì)心位置參數(shù)β2可增大機器人周期步態(tài)參數(shù)區(qū)間；當(dāng)β1取0.2時，隨著質(zhì)心位置參數(shù)β2的增大，機器人周期步態(tài)參數(shù)區(qū)間先增大后減?。划?dāng)β1取0.4時，增大質(zhì)心位置參數(shù)β2，機器人周期步態(tài)參數(shù)區(qū)間增大；而當(dāng)β1取0.6時，減小質(zhì)心位置參數(shù)β2可以增大機器人周期步態(tài)參數(shù)區(qū)間，且β2取0.6時機器人步態(tài)出現(xiàn)了逆倍周期分岔現(xiàn)象。

表1 機器人穩(wěn)定行走參數(shù)區(qū)間Tab.1 Parameter Ranges of Stable Robot Walking

圖8 β1和β2取值不同時的分岔圖Fig.8 Bifurcation Diagrams of Different Values of β1and β2

另一方面，當(dāng)質(zhì)心位置參數(shù)β2分別取0.5，0.55，0.6，且圓弧足半徑參數(shù)β1在（0，0.4）內(nèi)增大時有利于增加機器人周期步態(tài)參數(shù)區(qū)間，β1取0.6會導(dǎo)致周期參數(shù)區(qū)間減小。以上對圓弧足半徑參數(shù)β1和質(zhì)心位置參數(shù)β2組合變化下機器人穩(wěn)定參數(shù)區(qū)間分析為機器人結(jié)構(gòu)設(shè)計和穩(wěn)定行走提供參考。

6 結(jié)論

（1）詳細介紹了圓弧足被動行走機器人運動過程和動力學(xué)方程的建立。借助分岔圖深入討論了單參數(shù)變化影響，研究結(jié)果表明，當(dāng)圓弧足半徑參數(shù)β1在（0，0.805）內(nèi)增大時機器人仍保持周期一步態(tài)，但當(dāng)足半徑參數(shù)β1大于0.805時會導(dǎo)致機器人無法穩(wěn)定行走，而隨著質(zhì)心位置參數(shù)β2、轉(zhuǎn)動慣量參數(shù)β3和斜坡傾角φ增大，機器人步態(tài)均出現(xiàn)了倍周期分岔現(xiàn)象。

（2）雙參數(shù)變化影響研究表明，當(dāng)β1取0和0.6時，減小質(zhì)心位置參數(shù)β2可增大機器人周期步態(tài)參數(shù)區(qū)間，且β2取0.6時機器人步態(tài)出現(xiàn)了逆倍周期分岔現(xiàn)象。當(dāng)β1取0.4時，增大質(zhì)心位置參數(shù)β2，機器人周期步態(tài)參數(shù)區(qū)間增大；當(dāng)β1取0.2時，隨著質(zhì)心位置參數(shù)β2增大，機器人周期步態(tài)參數(shù)區(qū)間先增大后減??；當(dāng)質(zhì)心位置參數(shù)β2分別取0.5，0.55，0.6，且圓弧足半徑參數(shù)β1在（0，0.4）內(nèi)增大時有利于增加機器人周期步態(tài)參數(shù)區(qū)間，β1取0.6會導(dǎo)致周期參數(shù)區(qū)間減小。上述研究結(jié)果為圓弧足被動行走機器人結(jié)構(gòu)參數(shù)設(shè)計和穩(wěn)定步態(tài)分析與控制提供了重要參考意義。此外，利用正交擾動向量法求解被動行走機器人的李雅普諾夫指數(shù)，為未來研究其他非連續(xù)機械系統(tǒng)提供了思路。