Stability Analysis of a Simple Walking Model Driven by a Rhythmic ...
Proceedings of 2004 IEEE/RSJ International Conference on Intelligent Robots and Systems September 28 - October 2, 2004, Sendai, Japan
Stability Analysis of a Simple Walking Model Driven by a Rhythmic Signal Shinya Aoi and Kazuo Tsuchiya Dept. of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University Yoshida-honmachi, Sakyo-ku, Kyoto, 606-8501, Japan Email: {shinya aoi, tsuchiya}@kuaero.kyoto-u.ac.jp
Abstract— Human walk is achieved by the interaction between the dynamics of the human mechanical system and the rhythmic signals of the central pattern generator (CPG). In this paper, we analyze dynamic properties of a simple walking model of a biped robot driven by a rhythmic signal of an oscillator. This model consists of a hip and two legs which are connected at the hip. The swing leg is controlled by the rhythmic signal which is open loop. We analyze the stability of the periodic walking motion using a Poincar´e map. As a result, it is revealed that the simple walking model has the self-stabilization property, that is, the walking motion of the simple walking model converges to a periodic walking motion.
M g
u l θ2 θ1
y m
m x
I. I NTRODUCTION Walk of a legged robot can be considered to be an organized motion of a mechanical system which has many degrees of freedom. To reveal the laws of the dynamics which govern the walking motion, simple models are used in the literature. In human walk, it is revealed that the swing leg of human can be modeled as a double pendulum [8]. McGeer [7] developed a biped robot which has two legs, showed that the biped robot can walk stably down a slight slope without other energy input or control, and named the walking motion passive dynamic walking. The passive dynamic walking was analyzed using a 2-D model [4]. As a result, it revealed that the passive dynamic walking is a natural motion of this type of mechanism under gravity and has the self-stabilization property. That is, the passive dynamic walking of a biped robot has the property that the walking motion converges to a stable limit cycle due to the dynamic property of the mechanical system under gravity. The basin of attraction of the limit cycle has been researched [9] and the analysis based on a 3-D model has been performed [1]. On the other hand, a rhythmic motion such as a walking motion is achieved by the interaction between the dynamics of the mechanical system and the rhythmic signals of the central pattern generator (CPG) [10]. As a simple model of running motion of a legged robot driven by a rhythmic signal, the spring-loaded inverted pendulum (SLIP) is proposed [3]. Using SLIP, the stability and the basin of attraction of the running motion of biped robots has been analyzed [5]. In this paper, we propose a simple walking model of a biped robot driven by a rhythmic signal of an oscillator which is open loop. The simple walking model has a hip and two legs which are connected at the hip. We
0-7803-8463-6/04/$20.00 ©2004 IEEE
Fig. 1.
Schematic of the simple walking model of a biped robot
analyze the dynamic properties of the walking motion using the simple walking model. First, we consider the walking motion in which the tipping motion of the robot from the normal to the ground is small and obtain approximate periodic solutions of the walking motion. Next, we use a Poincar´e map and analyze the stability of the periodic motion. As a result, it is revealed that the simple walking model has the self-stabilization property, that is, the walking motion of the simple walking model converges to a periodic walking motion. II. A S IMPLE M ODEL A. Equations of motion for the swing phase In this paper, we consider a simple model of a biped robot which is shown in Fig. 1. The robot has a hip and two legs which are connected at the hip. We assume that the robot is constrained on the x-y plane and walks to the x-direction. Masses are concentrated at the hip and the tips of the legs. The hip mass is M , the leg mass is m, and the leg length is l. It has two degrees of freedom, θ1 and θ2 , which are functions of time t. Specifically, θ1 is the angle of the stance leg relative to the normal to the ground and is not directly controlled. On the other hand, θ2 is the angle of the swing leg relative to the stance leg and is directly controlled by the torque u. g is the acceleration due to gravity. The dimensionless equations of motion for the swing phase are given by 1 + 2β(1 − cos θ2 ) −β(1 − cos θ2 ) θ¨1 −β(1 − cos θ2 ) β θ¨2
1365
+
−β sin θ2 (θ˙22 − 2θ˙1 θ˙2 ) −βθ˙2 sin θ2
1
β sin(θ1 − θ2 ) − β sin θ1 − sin θ1 + −β sin(θ1 − θ2 ) 0 0 + = (1) u 0 where β ≡ m/M > 0, τ ≡ g/l t, u ≡ u/mgl, and ∗˙ indicates the derivative with respect to τ . In order to control the angle θ2 , we design a desired angle θ2d . First, we introduce an oscillator (Rejφ ), where R is the amplitude which is constant and φ is the phase whose angular velocity, ω > 0, is constant, that is,
R˙ = 0,
φ˙ = ω
(2)
Second, we design the desired angle θ2d which is driven by a rhythmic signal of the oscillator. In particular, θ2d is given by the simple function of the phase of the oscillator θ2d = θ2d (φ) = S cos φ
(3)
where S indicates the parameter which determines the stride. The torque u to the controlled angle θ2 is given so as to track the desired angle θ2d using the high-gain local feedback u = −Kp {θ2 − θ2d (φ)} − Kd {θ˙2 − θ˙2d (φ)}
(4)
where Kp and Kd are gain constants. In light of the above description, the state variables are defined as
III. A PPROXIMATE A NALYSIS A. Assumption for the controlled angle The angular velocity of the controlled angle θ˙2 becomes discontinuous according to the transition rules (6). This influence and the difference between the initial conditions of the angle θ2 and the desired angle θ2d which is given in (3) result in an error between the angle θ2 and the desired angle θ2d . But, as long as the angle θ2 is controlled by the sufficiently high-gain feedback control torque, we can assume that the angle θ2 can catch up with the desired angle θ2d in a sufficiently small period of time. Therefore, we assume that the angle θ2 is identical to the desired angle θ2d which is driven by a rhythmic signal of the oscillator, that is, θ2 (τ ) = θ2d (φ(τ ))
where Sc := { q | r(q) = 0 } and −θ1 θ˙1 2 cos 2θ1 θ˙1 f(q) = f2 (q) , h(q) = 2 + β(1 − cos 4θ1 ) ω φ−π f2 (q) = {β(1 − cos θ2d )θ¨2d + β sin θ2d (θ˙2 − 2θ˙1 θ˙2d )
2d
q T = [ θ1 , θ˙1 , θ2 , θ˙2 , φ ]
−β sin(θ1 − θ2d ) + β sin θ1 + sin θ1 }/ {1 + 2β(1 − cos θ2d )}
B. Collision model A collision occurs when the swing leg touches the ground. The geometric condition of the collision, the collision condition, is given by r(q) ≡ 2θ1 − θ2 = 0
(7)
Then, the state variables are redefined as q = [ θ1 , θ˙1 , φ ] and the equations of motion, the collision condition, and the transition rules are summarized by
/ Sc q˙ = f(q), q− ∈ (8) + − q = h(q ), q − ∈ Sc T
(5)
At the collision, there is an impulse at the tip of the swing leg. We assume the followings: 1) the swing leg neither slips and rebounds; 2) when the stance leg leaves the ground, there is no interaction with the ground; 3) the control torque is not impulsive; and 4) the swing leg instantaneously becomes the stance leg and vice versa. Then, the relations between the states just before the collision and the states just after the collision, the transition rules, are given as follows [6]: + −θ1− θ1 + 2 cos 2θ1− ˙− θ˙ 1 2 + β(1 − cos 4θ− ) θ1 1 + − θ = (6) −θ 2 2 + 2 cos 2θ− (1 − cos 2θ− ) 1 1 ˙− θ˙ 2 2 + β(1 − cos 4θ− ) θ1 1 φ+ φ− − π where ∗− indicates the state just before the collision, ∗+ indicates the state just after the collision, and the relationship between the states φ+ and φ− comes from + − the condition θ2d = −θ2d .
B. Approximate periodic solutions We analyze the periodic motion of the simple walking model under the assumption that the controlled angle θ2 is identical to the desired angle θ2d which is driven by a rhythmic signal of the oscillator. Since the equation (8) can not be solved directly due to the strong nonlinearity, we attempt to solve (8) approximately and obtain the approximate periodic solutions of the motion from just after a collision to just before the next collision. We consider the walking motion in which the tipping motion of the robot from the normal to the ground is small. In this case, the angles θ1 and θ2d are small and hence we define S ≡ , where is a small parameter. Then, the scaled variables Θ1 and Θ2d are defined as θ1 (τ ) = Θ1 (τ ), θ2d (φ(τ )) = Θ2d (φ(τ ))
(9)
The substitution of (9) into the second equation of the equations of motion (8) gives ¨ 1 − βΘ2d − Θ1 Θ
˙ 1Θ ˙ 2 Θ2d ¨ 1 − β Θ2 Θ ¨ 2d + 2β Θ ˙ 2d Θ2d − β Θ +3 βΘ22d Θ 2d 2 2d β β β 1 + Θ21 Θ2d − Θ1 Θ22d + Θ32d + Θ31 + O(5 ) = 0 (10) 2 2 6 6 where Θ1 = Θ1 (τ ) and Θ2d = Θ2d (φ(τ )). The step period is defined as T (= π/ω = O(0 )) which is the time interval
1366
from just after a collision to just before the next collision and the substitution of T and (9) into the collision condition and the transition rules (8) gives the boundary conditions for the approximate periodic solutions as follows: 1 Θ1 (T ) = − Θ2d (φ(0)) 2 1 Θ1 (0) = Θ2d (φ(0)) (11) 2 ˙ [2 + β {1 − cos(4Θ1 (T ))}] Θ1 (0) ˙ 1 (T ) = 2 cos(2Θ1 (T ))Θ where τ = 0 indicates the time just after a collision and τ = T indicates the time just before the next collision. Next, we consider power-series expansions for the states Θ1 and φ using the parameter µ ≡ 2 Θ1 (τ ) = X0 (τ ) + µX1 (τ ) + µ2 X2 (τ ) + · · · φ(τ ) = (ωτ + Φ0 ) + µΦ1 + µ2 Φ2 + · · ·
(12)
Moreover from (3) and (12) it follows that Θ2d (φ(τ )) = Y0 (τ ) + µY1 (τ ) + µ2 Y2 (τ ) + · · ·
The solutions to (16) are β e−τ 1 eτ + + X0 (τ ) = cos Φ0 ω2 + 1 2 1 − eT 1 − e−T β (17) − 2 cos(ωτ + Φ0 ) ω +1 Φ0 = 0 (Forward), π (Backward) From now on, we use Φ0 = 0 which indicates that the simple walking model walks forward. Next, from (14) and (15) the equations for X1 (τ ) and Φ1 are given by ¨0 − β Y 2 Y¨0 + 2β X˙ 0 Y˙ 0 Y0 ¨1 − βY1 − X1 + βY 2 X X 0 2 0 β β β 1 −βY˙02 Y0 + X02 Y0 − X0 Y02 + Y03 + X03 = 0 2 2 6 6 1 (18) X1 (T ) = − Y1 (0) 2 1 X1 (0) = Y1 (0) 2 ˙ X1 (0) + 4βX02 (T )X˙ 0 (0) = X˙ 1 (T ) − 2X02 (T )X˙ 0 (T ) The solution Φ1 to (18) is
(13) Φ1 =
where Y0 (τ ) Y1 (τ )
= cos(ωτ + Φ0 ) = − sin(ωτ + Φ0 )Φ1 .. .
The substitution of (12) and (13) into (10) gives
¨ 0 − βY0 − X0 + 3 X ¨1 − βY1 − X1 + βY 2 X ¨0 X 0 β β − Y02 Y¨0 + 2β X˙ 0 Y˙ 0 Y0 − β Y˙ 02 Y0 + X02 Y0 2 2 β β 3 1 3 2 − X0 Y0 + Y0 + X0 + O(5 ) = 0 (14) 2 6 6 where Xi = Xi (τ ), Yi = Yi (τ ) (i = 0, 1, · · ·). The substitution of (12) and (13) into (11) gives 1 (T ) + (0) Y X 0 0
2 1 3 X (T ) + (0) + O(5 ) = 0 Y + 1 1 2
1 X0 (0) − Y0 (0) (15)
2 1 5 3 + X1 (0) − Y1 (0) + O( ) = 0 2 X˙ 0 (0) − X˙ 0 (T ) + 3 X˙ 1 (0) + 4βX02 (T )X˙ 0 (0) −X˙ 1 (T ) + 2X 2 (T )X˙ 0 (T ) + O(5 ) = 0 0
From (14) and (15), the equations for X0 (τ ) and Φ0 are given by ¨0 (τ ) − βY0 (τ ) − X0 (τ ) = 0 X 1 X0 (T ) = − Y0 (0) 2 (16) 1 (0) = X Y 0 0 (0) 2 ˙ X0 (0) = X˙ 0 (T )
eT + 1 (1 + 2β)(1 + 2β + ω2 ) eT − 1 8βω
(19)
The solution X1 (τ ) is so complicated that we do not display the details here. It concludes that the approximate periodic solutions of the states θ1 and φ become β e−τ 1 eτ θ1 (τ ) = + + ω2 + 1 2 1 − eT 1 − e−T β − 2 cos(ωτ ) + · · · (20) ω +1 T e + 1 (1 + 2β)(1 + 2β + ω2 ) + · · · (21) φ(τ ) = ωτ + 2 T e −1 8βω Moreover, the desired angle θ2d becomes eT + 1 eT − 1 (1 + 2β)(1 + 2β + ω2 ) × sin(ωτ ) + · · · (22) 8βω
θ2d (φ(τ )) = cos(ωτ ) − 3
First, we show the comparisons of the approximate analysis (Approx. Analy.) and the numerical simulation (Num. Sim. (Fix)) in which we assume that the controlled angle θ2 is identical to the desired angle θ2d . The parameters are set as follows: g = 9.80665 [m/s2 ] and l = 0.1 [m]. Fig. 2a shows the periodic solutions with T = 0.1 g/l, β = 0.2, and = 6 [deg]. Fig. 2b and c show the results of θ1+ and φ+ , respectively, versus the parameter . Moreover, we show the result of the numerical simulation (Num. Sim. (F.B.)) in which we do not assume that the controlled angle θ2 is identical to the desired angle θ2d and control the angle θ2 by the high-gain feedback torque √ using the parameters Kp = 100/(lg) and Kd = 20/(l gl) (Fig. 2). From these results, it is verified that even if we assume that the controlled angle θ2 is identical to the desired angle θ2d , the approximate analysis is valid.
1367
0.15
To analyze stability with numerical simulations, we use a Newton-Raphson method to find the fixed point at least to O(10−7 ) and add a perturbation to O(10−4 ) so as to find the eigenvalues of the Jacobian matrix J(q ∗ ).
Approx. Analy. Num. Sim. (Fix) Num. Sim. (F.B.)
0.1 Angles [rad]
0.05
θ1
0 -0.05
θ2
-0.1 -0.15 0
0.02
0.04 0.06 Time [s]
0.08
0.1
(a) Angles versus time ( = 6 [deg]) 0.18
1
Approx. Analy. Num. Sim. (Fix) Num. Sim. (F.B.)
0.16 0.14
Approx. Analy. Num. Sim. (Fix) Num. Sim. (F.B.)
0.8
φ+ [rad]
θ+1 [rad]
0.12 0.1 0.08
0.6 0.4
0.06 0.04
0.2
0.02 0
0 0
5
10
15
20
0
5
[deg]
(b) θ1+ versus Fig. 2.
10
15
20
A. Jacobian matrix of the Poincar´e map In order to find the Jacobian matrix J(q ∗ ) of the Poincar´e map at the fixed point q ∗ of the periodic walking motion whose states include discontinuity due to a collision, we need to consider the influences of the collision. Coleman et al. show in [2] that the Jacobian matrix J(q ∗ ) is equivalent to the product of the three matrices B, D, and E which are shown below. First, we define the actual periodic solution as q ∗ (τ ), the actual step period as τ ∗ , and the perturbed state from the actual periodic solution q ∗ (τ ) from just after the collision to the next collision as q ∗ (τ ) + δ qˆi (τ ), where qˆi (τ ) is subject to the initial condition qˆi (0) = qˆi+ . Then, the matrices B and D are given by B = Dq h(q ∗ (τ ∗ )) q˙∗ (τ ∗ )Dq r(q ∗ (τ ∗ )) D=I− Dq r(q ∗ (τ ∗ ))q˙∗ (τ ∗ )
[deg]
(c) φ+ versus
Periodic solutions ( T = 0.1(g/l) 1/2, β = 0.2 )
IV. S TABILITY A NALYSIS One of the advantages to obtain the approximate periodic solutions as in the last section is to be able to approximately analyze the stability of the periodic solutions. In particular, in this paper we use a Poincar´e section in order to analyze stability. On a Poincar´e section, a periodic solution appears as a fixed point and the stability of the periodic solution can be determined by the eigenvalues of the Jacobian matrix of the return map around the fixed point. In this paper, the Poincar´e section is taken when a collision occurs and the state just after the collision is the state on the Poincar´e section. The return map from one point to the next point (Poincar´e map) is denoted as q → p(q), then + qi+1 = p(qi+ )
(23)
where qi+ is the state just after the ith collision. Note that the fixed point q ∗ on the Poincar´e section satisfies q ∗ = p(q ∗ )
(24)
Next, we add a perturbation δ qˆi+ just after a collision, where δ is small and ˆ ∗ indicates the perturbation. The state at the instant is defined as qi+ (= q ∗ + δ qˆi+ ) and the next + state on the Poincar´e section is defined as qi+1 (= p(qi+ ) = + ∗ q + δ qˆi+1). Linearizing the Poincar´e map p at the fixed point q ∗ , the Jacobian matrix J(q ∗ ) satisfies + qˆi+1 = J(q ∗ )ˆ qi+
(25)
Therefore, we can conclude that the periodic walk is asymptotically (resp., marginally) stable if all the eigenvalues of the Jacobian matrix J(q ∗ ) are inside (resp., inside or on) the unit circle on the complex plane, that is, all the magnitudes of the eigenvalues are less than (resp., less than or equal to) 1, and otherwise the periodic walk is unstable.
(26)
∂ where Dq ≡ ∂q . And the evolved perturbation after the + next collision qˆi+1 has the following relation: + = BDˆ qi (τ ∗ ) qˆi+1
(27) ∗
The substitution of the perturbed state q (τ ) + δ qˆi (τ ) into the equations of motion (8) gives qˆ˙ (τ ) Dq f(q ∗ (τ ))ˆ qi (τ ) (28) i
The matrix E is derived by integrating (28) as follows: qˆi (τ ∗ ) E qˆi+
(29)
From the results of the last section, the states θ1 and θ˙1 are O(1 ) and φ is O(0 ). Therefore, we define a ˙ perturbation as qˆT = [ θˆ1 , θˆ1 , φˆ ], then it is written by 0 0 θˆ1 θˆ1 qˆ = θˆ˙ 1 = 0 0 θˆ˙ 1 ≡ I˜q˜ (30) ˆ ˆ 0 0 1 φ φ Thus, (27) is modified by + ˜ D˜ ˜ qi (τ ∗ ) q˜i+1 =B
= where q˜i (τ ) = I˜−1 qˆi (τ ), ˜ D ˜ =I ˜ = I˜−1 B I, B q˜i+
(31)
I˜−1 qˆi+ , and ˜−1 ˜ DI
Also (28) is modified by q˜˙ i (τ ) I˜−1 Dq f(q ∗ (τ ))I˜q˜i (τ ) ∗ ≡ D qi (τ ) q f (q (τ ))˜
(32)
By integrating (32), we can obtain ˜ q˜+ q˜i (τ ∗ ) E
(33)
i
It follows from (31) and (33) that the eigenvalues of the Jacobian matrix J(q ∗ ) are given by ˜D ˜ E) ˜ Λ(J(q ∗ )) = Λ(B (34) where Λ(∗) indicates the spectrum of the matrix ∗.
1368
2
B. Approximate stability analysis of periodic motion In the previous section, we obtain the approximate periodic solutions (20) and (21). The substitution of the approximate solutions (20) and (21) to (31) and (33) gives ˜D ˜ E˜ by the power series of the parameter 2 the matrix B
where
0
0 = BDE ∗
0
0
2
1 − βω d
∗ ω(1 + ω2 )d
(35)
λ2 − 2(1 − βω2 d)λ + 1 = 0 Thus, the eigenvalues λ1,2 are given by λ1,2 = 1 − βω2 d ± (1 − βω2 d)2 − 1
(36)
(37)
Case 1. λ1,2 are complex conjugate In this case, it follows that (1 − βω2 d)2 − 1 < 0 ⇐⇒ 0 < βω2 d < 2 and the periodic solutions are marginally stable from the following: ¯ = |λ|2 = 1 λ1 λ2 = λλ
(38)
Case 2. λ1,2 are real and distinct In this case, it follows that (1 − βω2 d)2 − 1 > 0 ⇐⇒ βω2 d > 2 and the periodic solutions are unstable since there are stable and unstable eigenvalues from the following: λ2 − 1 = 2(1 − βω2 d)λ − 2 = 2 (1 − βω2 d)2 − 1 (1 − βω2 d)2 − 1 ± (1 − βω2 d) ≷ 0 (39) × Case 3. λ1,2 are degenerate In this case, it follows that (1 − βω2 d)2 − 1 = 0 ⇐⇒ βω2 d = 0, 2 and the periodic solutions are marginally stable from the following: λ1,2 = 1, −1
(40)
It concludes that the periodic solutions are marginally stable for βω2 d ≤ 2 and otherwise the periodic solutions are unstable. In order to analyze the relationship of the stability between the parameters ω and β, we define g0 (ω, β) ≡ βω2 d − 2. Then, it follows that ∂g0 β(1 − e−π/ω ) π(1 + 2β + ω2 ) (ω, β) = − ∂ω (1 + 2β + ω2 )2 ×(1 + eπ/ω ) + 2ω(1 + 2β)(1 − eπ/ω ) < 0 (41) ∂g0 eπ/ω + e−π/ω − 2 > 0 (42) (ω, β) = ω2 (1 + ω2 ) ∂β (1 + 2β + ω2 )2
1.2 1 Unstable
0.6 0.4 0
βω − (2 − βω2 d) 1 + ω2 1 − βω2 d
eT + e−T − 2 1 is so complicated that > 0, and BDE 1 + 2β + ω2 we do not display the details here. First, we analyze the stability of the periodic solutions 0 has one zero to O(0 ). From (35) the matrix BDE eigenvalue (λ3 = 0) and the other two eigenvalues λ1,2 are obtained from the equation
1.4
0.8
d≡
Marginally Stable
1.6
ω0(β) [rad]
˜D ˜ E˜ = BDE 1 + O(4 ) 0 + 2 BDE B
1.8
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 β
Fig. 3.
1
ω0 (β) versus β
Moreover, it follows that lim g0 (ω, β) = ∞ ω→+0 lim g0 (ω, β) = −2 ω→+∞
(43) lim g0 (ω, β) = −2 ∞
β→+0 2 2n 2ω 1 π >0 π2 − 4 + lim g0 (ω, β) = β→+∞ 2 (2n)! ω n=2 From (41) and (42), g0 (ω, β) is monotonically decreasing with respect to ω and monotonically increasing with respect to β. Therefore, from (43), ω > 0 which satisfies g0 (ω, β) = 0 exists and is unique with respect to each β > 0 and also β > 0 which satisfies g0 (ω, β) = 0 exists and is unique with respect to each ω > 0. From the implicit function theorem, we can define the functions ω0 (β) and ω0−1 (ω) such that g0 (ω0 (β), β) = g0 (ω, ω0−1 (ω)) = 0. Thus, it concludes that the periodic solutions are marginally stable for ω ≥ ω0 (β) or, equivalently, β ≤ ω0−1 (ω). Fig. 3 shows the plot of ω0 (β). Next, we analyze the stability of the periodic solutions to O(2 ) since in the analysis of the stability to O(0 ) it is only revealed that the periodic solutions are ˜D ˜E ˜ to O(2 ) marginally stable or unstable. The matrix B is so complicated that we do not display the details here. Fig. 4a shows the contour of the maximum magnitude ˜D ˜ E˜ to of the three eigenvalues λ1,2,3 of the matrix B O(2 ) with respect to ω and β for = 6 [deg] and the thick line indicates the function ω0 (β). From this figure, it reveals that there is a region where the periodic solutions are asymptotically stable while from the stability to O(0 ) it only reveals that the periodic solutions are marginally stable or unstable. Moreover, Fig. 4b and c show the comparisons of analytic and numerical results of the magnitudes of the eigenvalues with respect to ω (with β = 0.2) and β (with ω = 1.5 [rad]), respectively. Next, we analyze the stability of the periodic solutions with respect to the parameter . The stability of the periodic solutions to O(0 ) changes over the boundary ω = ω0 (β). On the contrary, the stability of the periodic solutions to O(2 ) changes over the boundary on which the maximum ˜D ˜E ˜ to O(2 ) magnitude of the eigenvalues of the matrix B is equal to 1 as shown in Fig. 4a. In general, we can not define the boundary over which the stability of the periodic solutions to O(2 ) changes as a function in the β-ω plane with respect to . Here, we assume that the
1369
2.2
1.8 0.93
1.8
0.95
0.96
1.6
ω [rad]
1.6
Asymptotically Stable 0.97
1.4
ω0
1.2 1
3.0
0.8
8.0
0.6
30.0
0.4 0.1
0.2
0.3
0.4
1.4
1.0
ω [rad]
2
1.2 Approx. Analy. = 0 [deg] = 6 [deg] (ω(β)) =12 [deg] =18 [deg] Num.Sim.(Fix) = 6 [deg]
1
Unstable
0.8 0.5
0.6
0.7
0.8
0.9
0.1
1
0.2
0.3
0.4
β
2.5
2 1.5 1 0.5
Fig. 5.
Approx. Analy. Num. Sim. (Fix)
1.5 1 0.5 0
0 0.9
1
1.1
1.2 1.3 ω [rad]
1.4
1.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
β
(b) Magnitudes of the eigenvalues versus ω (with β = 0.2) Fig. 4.
1.6
0.7
0.8
0.9
1
Analytic result, ω (β) versus β, and numerical result of the boundary of the stability region 2
2
Magnitudes of Eigenvalues
Approx. Analy. Num. Sim. (Fix)
Magnitudes of Eigenvalues
Magnitudes of Eigenvalues
3
0.6
β
(a) Contour of the maximum magnitude of the eigenvalues
2.5
0.5
(c) Magnitudes of the eigenvalues versus β (with ω = 1.5 [rad])
Approx. Analy. Num. Sim. (Fix)
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
Eigenvalues ( = 6 [deg])
0
5
10
15
20
25
30
[deg]
influences of the parameter upon the boundary are small since the parameter is small. Therefore, we assume that the difference between the boundaries to O(0 ) and O(2 ) is small as shown in Fig. 4a and hence we define the function ω (β) with respect to each > 0 such that max |λi | = 1 for ω = ω(β) i=1,2,3 max |λi | < 1 for ω > ω(β) (44) i=1,2,3 max |λi | > 1 for ω < ω(β) i=1,2,3
Fig. 5 shows the comparison of the analytic results of the function ω (β) and the numerical results of the boundary of the stability region for some values of . Moreover, Fig. 6 shows the comparison of the results of the approximate analysis and the numerical simulation of the magnitudes of the eigenvalues with respect to with β = 0.2 and ω = π/0.1 l/g [rad]. From these results, it is revealed that the asymptotic stability region is enlarged due to the parameter and it is verified that the stability analysis of the periodic motion is valid especially in the range where the tipping motion of the robot from the normal to the ground is small. V. C ONCLUSION In this paper, we analyzed the walking motion of a biped robot driven by a rhythmic signal of the oscillator based on the simple walking model. As a result, it is revealed that the simple walking model has the self-stabilization property. But this model has several limitations. One is that the rhythmic signal of the oscillator has no feedback from the dynamics of the biped robot. Another is that the pitch motion of the biped robot of this model has no influences from the roll motion of the biped robot since the biped robot of this model is constrained on a plane. Therefore,
Fig. 6.
Magnitudes of the eigenvalues versus (β = 0.2, ω = π/0.1(l/g)1/2 [rad])
we will consider to construct a model which includes these influences and provide the analysis of these influences. ACKNOWLEDGMENT This paper is supported in part by Center of Excellence for Research and Education on Complex Functional Mechanical Systems (COE program of the Ministry of Education, Culture, Sports, Science and Technology, Japan). R EFERENCES [1] J. Adolfsson, H. Dankowicz, and A. Nordmark, 3D passive walkers: Finding periodic gaits in the presence of discontinuities, Nonlinear Dynamics 24, pp. 205-229, 2001. [2] M. Coleman, A. Chatterjee, and A. Ruina, Motions of a rimless spoked wheel: A simple three-dimensional system with impacts, Dynamics and Stability of Systems, Vol 12, No. 3, pp. 139-160, 1997. [3] R.J. Full and D.E. Koditschek, Templates and anchors: Neuromechanical hypotheses of legged locomotion on land, J. Exp. Biol. 202, pp. 3325-3332, 1999. [4] M. Garcia, A. Chatterjee, A. Ruina, and M. Coleman, The simplest walking model: Stability, complexity, and scaling, ASME J. of Biomechanical Eng., Vol. 120, No. 2, pp. 281-288, 1998. [5] R. Ghigliazza, R. Altendorfer, P. Holmes, and D.E. Koditschek, A simply stabilized running model, SIAM J. on Applied Dynamical Systems, Vol. 2, No. 2, pp. 187-218, 2003. [6] J.W. Grizzle, G. Abba, and F. Plestan, Asymptotically stable walking for biped robots: Analysis via systems with impulse effects, IEEE Trans. Automatic Control, Vol. 46, No. 1, pp. 51-64, 2001. [7] T. McGeer, Passive dynamic walking, Int. J. of Robotics Res., Vol. 9, No. 2, pp. 62-82, 1990. [8] S. Mochon and T. McMahon, Ballistic walking: An improved model, Math. biosc. 52, pp. 241-260, 1980. [9] A.L. Schwab and M. Wisse, Basin of attraction of the simplest walking model, Proc. of ASME Int. Conf. on Noise and Vibration, 2001. [10] G. Taga, A model of the neuro-musculo-skeletal system for human locomotion II. - Real-time adaptability under various constraints, Biol. Cybern. 73, pp. 113-121, 1995.
1370
Stability Analysis of a Simple Walking Model Driven by a Rhythmic Signal Shinya Aoi and Kazuo Tsuchiya Dept. of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University Yoshida-honmachi, Sakyo-ku, Kyoto, 606-8501, Japan Email: {shinya aoi, tsuchiya}@kuaero.kyoto-u.ac.jp
Abstract— Human walk is achieved by the interaction between the dynamics of the human mechanical system and the rhythmic signals of the central pattern generator (CPG). In this paper, we analyze dynamic properties of a simple walking model of a biped robot driven by a rhythmic signal of an oscillator. This model consists of a hip and two legs which are connected at the hip. The swing leg is controlled by the rhythmic signal which is open loop. We analyze the stability of the periodic walking motion using a Poincar´e map. As a result, it is revealed that the simple walking model has the self-stabilization property, that is, the walking motion of the simple walking model converges to a periodic walking motion.
M g
u l θ2 θ1
y m
m x
I. I NTRODUCTION Walk of a legged robot can be considered to be an organized motion of a mechanical system which has many degrees of freedom. To reveal the laws of the dynamics which govern the walking motion, simple models are used in the literature. In human walk, it is revealed that the swing leg of human can be modeled as a double pendulum [8]. McGeer [7] developed a biped robot which has two legs, showed that the biped robot can walk stably down a slight slope without other energy input or control, and named the walking motion passive dynamic walking. The passive dynamic walking was analyzed using a 2-D model [4]. As a result, it revealed that the passive dynamic walking is a natural motion of this type of mechanism under gravity and has the self-stabilization property. That is, the passive dynamic walking of a biped robot has the property that the walking motion converges to a stable limit cycle due to the dynamic property of the mechanical system under gravity. The basin of attraction of the limit cycle has been researched [9] and the analysis based on a 3-D model has been performed [1]. On the other hand, a rhythmic motion such as a walking motion is achieved by the interaction between the dynamics of the mechanical system and the rhythmic signals of the central pattern generator (CPG) [10]. As a simple model of running motion of a legged robot driven by a rhythmic signal, the spring-loaded inverted pendulum (SLIP) is proposed [3]. Using SLIP, the stability and the basin of attraction of the running motion of biped robots has been analyzed [5]. In this paper, we propose a simple walking model of a biped robot driven by a rhythmic signal of an oscillator which is open loop. The simple walking model has a hip and two legs which are connected at the hip. We
0-7803-8463-6/04/$20.00 ©2004 IEEE
Fig. 1.
Schematic of the simple walking model of a biped robot
analyze the dynamic properties of the walking motion using the simple walking model. First, we consider the walking motion in which the tipping motion of the robot from the normal to the ground is small and obtain approximate periodic solutions of the walking motion. Next, we use a Poincar´e map and analyze the stability of the periodic motion. As a result, it is revealed that the simple walking model has the self-stabilization property, that is, the walking motion of the simple walking model converges to a periodic walking motion. II. A S IMPLE M ODEL A. Equations of motion for the swing phase In this paper, we consider a simple model of a biped robot which is shown in Fig. 1. The robot has a hip and two legs which are connected at the hip. We assume that the robot is constrained on the x-y plane and walks to the x-direction. Masses are concentrated at the hip and the tips of the legs. The hip mass is M , the leg mass is m, and the leg length is l. It has two degrees of freedom, θ1 and θ2 , which are functions of time t. Specifically, θ1 is the angle of the stance leg relative to the normal to the ground and is not directly controlled. On the other hand, θ2 is the angle of the swing leg relative to the stance leg and is directly controlled by the torque u. g is the acceleration due to gravity. The dimensionless equations of motion for the swing phase are given by 1 + 2β(1 − cos θ2 ) −β(1 − cos θ2 ) θ¨1 −β(1 − cos θ2 ) β θ¨2
1365
+
−β sin θ2 (θ˙22 − 2θ˙1 θ˙2 ) −βθ˙2 sin θ2
1
β sin(θ1 − θ2 ) − β sin θ1 − sin θ1 + −β sin(θ1 − θ2 ) 0 0 + = (1) u 0 where β ≡ m/M > 0, τ ≡ g/l t, u ≡ u/mgl, and ∗˙ indicates the derivative with respect to τ . In order to control the angle θ2 , we design a desired angle θ2d . First, we introduce an oscillator (Rejφ ), where R is the amplitude which is constant and φ is the phase whose angular velocity, ω > 0, is constant, that is,
R˙ = 0,
φ˙ = ω
(2)
Second, we design the desired angle θ2d which is driven by a rhythmic signal of the oscillator. In particular, θ2d is given by the simple function of the phase of the oscillator θ2d = θ2d (φ) = S cos φ
(3)
where S indicates the parameter which determines the stride. The torque u to the controlled angle θ2 is given so as to track the desired angle θ2d using the high-gain local feedback u = −Kp {θ2 − θ2d (φ)} − Kd {θ˙2 − θ˙2d (φ)}
(4)
where Kp and Kd are gain constants. In light of the above description, the state variables are defined as
III. A PPROXIMATE A NALYSIS A. Assumption for the controlled angle The angular velocity of the controlled angle θ˙2 becomes discontinuous according to the transition rules (6). This influence and the difference between the initial conditions of the angle θ2 and the desired angle θ2d which is given in (3) result in an error between the angle θ2 and the desired angle θ2d . But, as long as the angle θ2 is controlled by the sufficiently high-gain feedback control torque, we can assume that the angle θ2 can catch up with the desired angle θ2d in a sufficiently small period of time. Therefore, we assume that the angle θ2 is identical to the desired angle θ2d which is driven by a rhythmic signal of the oscillator, that is, θ2 (τ ) = θ2d (φ(τ ))
where Sc := { q | r(q) = 0 } and −θ1 θ˙1 2 cos 2θ1 θ˙1 f(q) = f2 (q) , h(q) = 2 + β(1 − cos 4θ1 ) ω φ−π f2 (q) = {β(1 − cos θ2d )θ¨2d + β sin θ2d (θ˙2 − 2θ˙1 θ˙2d )
2d
q T = [ θ1 , θ˙1 , θ2 , θ˙2 , φ ]
−β sin(θ1 − θ2d ) + β sin θ1 + sin θ1 }/ {1 + 2β(1 − cos θ2d )}
B. Collision model A collision occurs when the swing leg touches the ground. The geometric condition of the collision, the collision condition, is given by r(q) ≡ 2θ1 − θ2 = 0
(7)
Then, the state variables are redefined as q = [ θ1 , θ˙1 , φ ] and the equations of motion, the collision condition, and the transition rules are summarized by
/ Sc q˙ = f(q), q− ∈ (8) + − q = h(q ), q − ∈ Sc T
(5)
At the collision, there is an impulse at the tip of the swing leg. We assume the followings: 1) the swing leg neither slips and rebounds; 2) when the stance leg leaves the ground, there is no interaction with the ground; 3) the control torque is not impulsive; and 4) the swing leg instantaneously becomes the stance leg and vice versa. Then, the relations between the states just before the collision and the states just after the collision, the transition rules, are given as follows [6]: + −θ1− θ1 + 2 cos 2θ1− ˙− θ˙ 1 2 + β(1 − cos 4θ− ) θ1 1 + − θ = (6) −θ 2 2 + 2 cos 2θ− (1 − cos 2θ− ) 1 1 ˙− θ˙ 2 2 + β(1 − cos 4θ− ) θ1 1 φ+ φ− − π where ∗− indicates the state just before the collision, ∗+ indicates the state just after the collision, and the relationship between the states φ+ and φ− comes from + − the condition θ2d = −θ2d .
B. Approximate periodic solutions We analyze the periodic motion of the simple walking model under the assumption that the controlled angle θ2 is identical to the desired angle θ2d which is driven by a rhythmic signal of the oscillator. Since the equation (8) can not be solved directly due to the strong nonlinearity, we attempt to solve (8) approximately and obtain the approximate periodic solutions of the motion from just after a collision to just before the next collision. We consider the walking motion in which the tipping motion of the robot from the normal to the ground is small. In this case, the angles θ1 and θ2d are small and hence we define S ≡ , where is a small parameter. Then, the scaled variables Θ1 and Θ2d are defined as θ1 (τ ) = Θ1 (τ ), θ2d (φ(τ )) = Θ2d (φ(τ ))
(9)
The substitution of (9) into the second equation of the equations of motion (8) gives ¨ 1 − βΘ2d − Θ1 Θ
˙ 1Θ ˙ 2 Θ2d ¨ 1 − β Θ2 Θ ¨ 2d + 2β Θ ˙ 2d Θ2d − β Θ +3 βΘ22d Θ 2d 2 2d β β β 1 + Θ21 Θ2d − Θ1 Θ22d + Θ32d + Θ31 + O(5 ) = 0 (10) 2 2 6 6 where Θ1 = Θ1 (τ ) and Θ2d = Θ2d (φ(τ )). The step period is defined as T (= π/ω = O(0 )) which is the time interval
1366
from just after a collision to just before the next collision and the substitution of T and (9) into the collision condition and the transition rules (8) gives the boundary conditions for the approximate periodic solutions as follows: 1 Θ1 (T ) = − Θ2d (φ(0)) 2 1 Θ1 (0) = Θ2d (φ(0)) (11) 2 ˙ [2 + β {1 − cos(4Θ1 (T ))}] Θ1 (0) ˙ 1 (T ) = 2 cos(2Θ1 (T ))Θ where τ = 0 indicates the time just after a collision and τ = T indicates the time just before the next collision. Next, we consider power-series expansions for the states Θ1 and φ using the parameter µ ≡ 2 Θ1 (τ ) = X0 (τ ) + µX1 (τ ) + µ2 X2 (τ ) + · · · φ(τ ) = (ωτ + Φ0 ) + µΦ1 + µ2 Φ2 + · · ·
(12)
Moreover from (3) and (12) it follows that Θ2d (φ(τ )) = Y0 (τ ) + µY1 (τ ) + µ2 Y2 (τ ) + · · ·
The solutions to (16) are β e−τ 1 eτ + + X0 (τ ) = cos Φ0 ω2 + 1 2 1 − eT 1 − e−T β (17) − 2 cos(ωτ + Φ0 ) ω +1 Φ0 = 0 (Forward), π (Backward) From now on, we use Φ0 = 0 which indicates that the simple walking model walks forward. Next, from (14) and (15) the equations for X1 (τ ) and Φ1 are given by ¨0 − β Y 2 Y¨0 + 2β X˙ 0 Y˙ 0 Y0 ¨1 − βY1 − X1 + βY 2 X X 0 2 0 β β β 1 −βY˙02 Y0 + X02 Y0 − X0 Y02 + Y03 + X03 = 0 2 2 6 6 1 (18) X1 (T ) = − Y1 (0) 2 1 X1 (0) = Y1 (0) 2 ˙ X1 (0) + 4βX02 (T )X˙ 0 (0) = X˙ 1 (T ) − 2X02 (T )X˙ 0 (T ) The solution Φ1 to (18) is
(13) Φ1 =
where Y0 (τ ) Y1 (τ )
= cos(ωτ + Φ0 ) = − sin(ωτ + Φ0 )Φ1 .. .
The substitution of (12) and (13) into (10) gives
¨ 0 − βY0 − X0 + 3 X ¨1 − βY1 − X1 + βY 2 X ¨0 X 0 β β − Y02 Y¨0 + 2β X˙ 0 Y˙ 0 Y0 − β Y˙ 02 Y0 + X02 Y0 2 2 β β 3 1 3 2 − X0 Y0 + Y0 + X0 + O(5 ) = 0 (14) 2 6 6 where Xi = Xi (τ ), Yi = Yi (τ ) (i = 0, 1, · · ·). The substitution of (12) and (13) into (11) gives 1 (T ) + (0) Y X 0 0
2 1 3 X (T ) + (0) + O(5 ) = 0 Y + 1 1 2
1 X0 (0) − Y0 (0) (15)
2 1 5 3 + X1 (0) − Y1 (0) + O( ) = 0 2 X˙ 0 (0) − X˙ 0 (T ) + 3 X˙ 1 (0) + 4βX02 (T )X˙ 0 (0) −X˙ 1 (T ) + 2X 2 (T )X˙ 0 (T ) + O(5 ) = 0 0
From (14) and (15), the equations for X0 (τ ) and Φ0 are given by ¨0 (τ ) − βY0 (τ ) − X0 (τ ) = 0 X 1 X0 (T ) = − Y0 (0) 2 (16) 1 (0) = X Y 0 0 (0) 2 ˙ X0 (0) = X˙ 0 (T )
eT + 1 (1 + 2β)(1 + 2β + ω2 ) eT − 1 8βω
(19)
The solution X1 (τ ) is so complicated that we do not display the details here. It concludes that the approximate periodic solutions of the states θ1 and φ become β e−τ 1 eτ θ1 (τ ) = + + ω2 + 1 2 1 − eT 1 − e−T β − 2 cos(ωτ ) + · · · (20) ω +1 T e + 1 (1 + 2β)(1 + 2β + ω2 ) + · · · (21) φ(τ ) = ωτ + 2 T e −1 8βω Moreover, the desired angle θ2d becomes eT + 1 eT − 1 (1 + 2β)(1 + 2β + ω2 ) × sin(ωτ ) + · · · (22) 8βω
θ2d (φ(τ )) = cos(ωτ ) − 3
First, we show the comparisons of the approximate analysis (Approx. Analy.) and the numerical simulation (Num. Sim. (Fix)) in which we assume that the controlled angle θ2 is identical to the desired angle θ2d . The parameters are set as follows: g = 9.80665 [m/s2 ] and l = 0.1 [m]. Fig. 2a shows the periodic solutions with T = 0.1 g/l, β = 0.2, and = 6 [deg]. Fig. 2b and c show the results of θ1+ and φ+ , respectively, versus the parameter . Moreover, we show the result of the numerical simulation (Num. Sim. (F.B.)) in which we do not assume that the controlled angle θ2 is identical to the desired angle θ2d and control the angle θ2 by the high-gain feedback torque √ using the parameters Kp = 100/(lg) and Kd = 20/(l gl) (Fig. 2). From these results, it is verified that even if we assume that the controlled angle θ2 is identical to the desired angle θ2d , the approximate analysis is valid.
1367
0.15
To analyze stability with numerical simulations, we use a Newton-Raphson method to find the fixed point at least to O(10−7 ) and add a perturbation to O(10−4 ) so as to find the eigenvalues of the Jacobian matrix J(q ∗ ).
Approx. Analy. Num. Sim. (Fix) Num. Sim. (F.B.)
0.1 Angles [rad]
0.05
θ1
0 -0.05
θ2
-0.1 -0.15 0
0.02
0.04 0.06 Time [s]
0.08
0.1
(a) Angles versus time ( = 6 [deg]) 0.18
1
Approx. Analy. Num. Sim. (Fix) Num. Sim. (F.B.)
0.16 0.14
Approx. Analy. Num. Sim. (Fix) Num. Sim. (F.B.)
0.8
φ+ [rad]
θ+1 [rad]
0.12 0.1 0.08
0.6 0.4
0.06 0.04
0.2
0.02 0
0 0
5
10
15
20
0
5
[deg]
(b) θ1+ versus Fig. 2.
10
15
20
A. Jacobian matrix of the Poincar´e map In order to find the Jacobian matrix J(q ∗ ) of the Poincar´e map at the fixed point q ∗ of the periodic walking motion whose states include discontinuity due to a collision, we need to consider the influences of the collision. Coleman et al. show in [2] that the Jacobian matrix J(q ∗ ) is equivalent to the product of the three matrices B, D, and E which are shown below. First, we define the actual periodic solution as q ∗ (τ ), the actual step period as τ ∗ , and the perturbed state from the actual periodic solution q ∗ (τ ) from just after the collision to the next collision as q ∗ (τ ) + δ qˆi (τ ), where qˆi (τ ) is subject to the initial condition qˆi (0) = qˆi+ . Then, the matrices B and D are given by B = Dq h(q ∗ (τ ∗ )) q˙∗ (τ ∗ )Dq r(q ∗ (τ ∗ )) D=I− Dq r(q ∗ (τ ∗ ))q˙∗ (τ ∗ )
[deg]
(c) φ+ versus
Periodic solutions ( T = 0.1(g/l) 1/2, β = 0.2 )
IV. S TABILITY A NALYSIS One of the advantages to obtain the approximate periodic solutions as in the last section is to be able to approximately analyze the stability of the periodic solutions. In particular, in this paper we use a Poincar´e section in order to analyze stability. On a Poincar´e section, a periodic solution appears as a fixed point and the stability of the periodic solution can be determined by the eigenvalues of the Jacobian matrix of the return map around the fixed point. In this paper, the Poincar´e section is taken when a collision occurs and the state just after the collision is the state on the Poincar´e section. The return map from one point to the next point (Poincar´e map) is denoted as q → p(q), then + qi+1 = p(qi+ )
(23)
where qi+ is the state just after the ith collision. Note that the fixed point q ∗ on the Poincar´e section satisfies q ∗ = p(q ∗ )
(24)
Next, we add a perturbation δ qˆi+ just after a collision, where δ is small and ˆ ∗ indicates the perturbation. The state at the instant is defined as qi+ (= q ∗ + δ qˆi+ ) and the next + state on the Poincar´e section is defined as qi+1 (= p(qi+ ) = + ∗ q + δ qˆi+1). Linearizing the Poincar´e map p at the fixed point q ∗ , the Jacobian matrix J(q ∗ ) satisfies + qˆi+1 = J(q ∗ )ˆ qi+
(25)
Therefore, we can conclude that the periodic walk is asymptotically (resp., marginally) stable if all the eigenvalues of the Jacobian matrix J(q ∗ ) are inside (resp., inside or on) the unit circle on the complex plane, that is, all the magnitudes of the eigenvalues are less than (resp., less than or equal to) 1, and otherwise the periodic walk is unstable.
(26)
∂ where Dq ≡ ∂q . And the evolved perturbation after the + next collision qˆi+1 has the following relation: + = BDˆ qi (τ ∗ ) qˆi+1
(27) ∗
The substitution of the perturbed state q (τ ) + δ qˆi (τ ) into the equations of motion (8) gives qˆ˙ (τ ) Dq f(q ∗ (τ ))ˆ qi (τ ) (28) i
The matrix E is derived by integrating (28) as follows: qˆi (τ ∗ ) E qˆi+
(29)
From the results of the last section, the states θ1 and θ˙1 are O(1 ) and φ is O(0 ). Therefore, we define a ˙ perturbation as qˆT = [ θˆ1 , θˆ1 , φˆ ], then it is written by 0 0 θˆ1 θˆ1 qˆ = θˆ˙ 1 = 0 0 θˆ˙ 1 ≡ I˜q˜ (30) ˆ ˆ 0 0 1 φ φ Thus, (27) is modified by + ˜ D˜ ˜ qi (τ ∗ ) q˜i+1 =B
= where q˜i (τ ) = I˜−1 qˆi (τ ), ˜ D ˜ =I ˜ = I˜−1 B I, B q˜i+
(31)
I˜−1 qˆi+ , and ˜−1 ˜ DI
Also (28) is modified by q˜˙ i (τ ) I˜−1 Dq f(q ∗ (τ ))I˜q˜i (τ ) ∗ ≡ D qi (τ ) q f (q (τ ))˜
(32)
By integrating (32), we can obtain ˜ q˜+ q˜i (τ ∗ ) E
(33)
i
It follows from (31) and (33) that the eigenvalues of the Jacobian matrix J(q ∗ ) are given by ˜D ˜ E) ˜ Λ(J(q ∗ )) = Λ(B (34) where Λ(∗) indicates the spectrum of the matrix ∗.
1368
2
B. Approximate stability analysis of periodic motion In the previous section, we obtain the approximate periodic solutions (20) and (21). The substitution of the approximate solutions (20) and (21) to (31) and (33) gives ˜D ˜ E˜ by the power series of the parameter 2 the matrix B
where
0
0 = BDE ∗
0
0
2
1 − βω d
∗ ω(1 + ω2 )d
(35)
λ2 − 2(1 − βω2 d)λ + 1 = 0 Thus, the eigenvalues λ1,2 are given by λ1,2 = 1 − βω2 d ± (1 − βω2 d)2 − 1
(36)
(37)
Case 1. λ1,2 are complex conjugate In this case, it follows that (1 − βω2 d)2 − 1 < 0 ⇐⇒ 0 < βω2 d < 2 and the periodic solutions are marginally stable from the following: ¯ = |λ|2 = 1 λ1 λ2 = λλ
(38)
Case 2. λ1,2 are real and distinct In this case, it follows that (1 − βω2 d)2 − 1 > 0 ⇐⇒ βω2 d > 2 and the periodic solutions are unstable since there are stable and unstable eigenvalues from the following: λ2 − 1 = 2(1 − βω2 d)λ − 2 = 2 (1 − βω2 d)2 − 1 (1 − βω2 d)2 − 1 ± (1 − βω2 d) ≷ 0 (39) × Case 3. λ1,2 are degenerate In this case, it follows that (1 − βω2 d)2 − 1 = 0 ⇐⇒ βω2 d = 0, 2 and the periodic solutions are marginally stable from the following: λ1,2 = 1, −1
(40)
It concludes that the periodic solutions are marginally stable for βω2 d ≤ 2 and otherwise the periodic solutions are unstable. In order to analyze the relationship of the stability between the parameters ω and β, we define g0 (ω, β) ≡ βω2 d − 2. Then, it follows that ∂g0 β(1 − e−π/ω ) π(1 + 2β + ω2 ) (ω, β) = − ∂ω (1 + 2β + ω2 )2 ×(1 + eπ/ω ) + 2ω(1 + 2β)(1 − eπ/ω ) < 0 (41) ∂g0 eπ/ω + e−π/ω − 2 > 0 (42) (ω, β) = ω2 (1 + ω2 ) ∂β (1 + 2β + ω2 )2
1.2 1 Unstable
0.6 0.4 0
βω − (2 − βω2 d) 1 + ω2 1 − βω2 d
eT + e−T − 2 1 is so complicated that > 0, and BDE 1 + 2β + ω2 we do not display the details here. First, we analyze the stability of the periodic solutions 0 has one zero to O(0 ). From (35) the matrix BDE eigenvalue (λ3 = 0) and the other two eigenvalues λ1,2 are obtained from the equation
1.4
0.8
d≡
Marginally Stable
1.6
ω0(β) [rad]
˜D ˜ E˜ = BDE 1 + O(4 ) 0 + 2 BDE B
1.8
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 β
Fig. 3.
1
ω0 (β) versus β
Moreover, it follows that lim g0 (ω, β) = ∞ ω→+0 lim g0 (ω, β) = −2 ω→+∞
(43) lim g0 (ω, β) = −2 ∞
β→+0 2 2n 2ω 1 π >0 π2 − 4 + lim g0 (ω, β) = β→+∞ 2 (2n)! ω n=2 From (41) and (42), g0 (ω, β) is monotonically decreasing with respect to ω and monotonically increasing with respect to β. Therefore, from (43), ω > 0 which satisfies g0 (ω, β) = 0 exists and is unique with respect to each β > 0 and also β > 0 which satisfies g0 (ω, β) = 0 exists and is unique with respect to each ω > 0. From the implicit function theorem, we can define the functions ω0 (β) and ω0−1 (ω) such that g0 (ω0 (β), β) = g0 (ω, ω0−1 (ω)) = 0. Thus, it concludes that the periodic solutions are marginally stable for ω ≥ ω0 (β) or, equivalently, β ≤ ω0−1 (ω). Fig. 3 shows the plot of ω0 (β). Next, we analyze the stability of the periodic solutions to O(2 ) since in the analysis of the stability to O(0 ) it is only revealed that the periodic solutions are ˜D ˜E ˜ to O(2 ) marginally stable or unstable. The matrix B is so complicated that we do not display the details here. Fig. 4a shows the contour of the maximum magnitude ˜D ˜ E˜ to of the three eigenvalues λ1,2,3 of the matrix B O(2 ) with respect to ω and β for = 6 [deg] and the thick line indicates the function ω0 (β). From this figure, it reveals that there is a region where the periodic solutions are asymptotically stable while from the stability to O(0 ) it only reveals that the periodic solutions are marginally stable or unstable. Moreover, Fig. 4b and c show the comparisons of analytic and numerical results of the magnitudes of the eigenvalues with respect to ω (with β = 0.2) and β (with ω = 1.5 [rad]), respectively. Next, we analyze the stability of the periodic solutions with respect to the parameter . The stability of the periodic solutions to O(0 ) changes over the boundary ω = ω0 (β). On the contrary, the stability of the periodic solutions to O(2 ) changes over the boundary on which the maximum ˜D ˜E ˜ to O(2 ) magnitude of the eigenvalues of the matrix B is equal to 1 as shown in Fig. 4a. In general, we can not define the boundary over which the stability of the periodic solutions to O(2 ) changes as a function in the β-ω plane with respect to . Here, we assume that the
1369
2.2
1.8 0.93
1.8
0.95
0.96
1.6
ω [rad]
1.6
Asymptotically Stable 0.97
1.4
ω0
1.2 1
3.0
0.8
8.0
0.6
30.0
0.4 0.1
0.2
0.3
0.4
1.4
1.0
ω [rad]
2
1.2 Approx. Analy. = 0 [deg] = 6 [deg] (ω(β)) =12 [deg] =18 [deg] Num.Sim.(Fix) = 6 [deg]
1
Unstable
0.8 0.5
0.6
0.7
0.8
0.9
0.1
1
0.2
0.3
0.4
β
2.5
2 1.5 1 0.5
Fig. 5.
Approx. Analy. Num. Sim. (Fix)
1.5 1 0.5 0
0 0.9
1
1.1
1.2 1.3 ω [rad]
1.4
1.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
β
(b) Magnitudes of the eigenvalues versus ω (with β = 0.2) Fig. 4.
1.6
0.7
0.8
0.9
1
Analytic result, ω (β) versus β, and numerical result of the boundary of the stability region 2
2
Magnitudes of Eigenvalues
Approx. Analy. Num. Sim. (Fix)
Magnitudes of Eigenvalues
Magnitudes of Eigenvalues
3
0.6
β
(a) Contour of the maximum magnitude of the eigenvalues
2.5
0.5
(c) Magnitudes of the eigenvalues versus β (with ω = 1.5 [rad])
Approx. Analy. Num. Sim. (Fix)
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
Eigenvalues ( = 6 [deg])
0
5
10
15
20
25
30
[deg]
influences of the parameter upon the boundary are small since the parameter is small. Therefore, we assume that the difference between the boundaries to O(0 ) and O(2 ) is small as shown in Fig. 4a and hence we define the function ω (β) with respect to each > 0 such that max |λi | = 1 for ω = ω(β) i=1,2,3 max |λi | < 1 for ω > ω(β) (44) i=1,2,3 max |λi | > 1 for ω < ω(β) i=1,2,3
Fig. 5 shows the comparison of the analytic results of the function ω (β) and the numerical results of the boundary of the stability region for some values of . Moreover, Fig. 6 shows the comparison of the results of the approximate analysis and the numerical simulation of the magnitudes of the eigenvalues with respect to with β = 0.2 and ω = π/0.1 l/g [rad]. From these results, it is revealed that the asymptotic stability region is enlarged due to the parameter and it is verified that the stability analysis of the periodic motion is valid especially in the range where the tipping motion of the robot from the normal to the ground is small. V. C ONCLUSION In this paper, we analyzed the walking motion of a biped robot driven by a rhythmic signal of the oscillator based on the simple walking model. As a result, it is revealed that the simple walking model has the self-stabilization property. But this model has several limitations. One is that the rhythmic signal of the oscillator has no feedback from the dynamics of the biped robot. Another is that the pitch motion of the biped robot of this model has no influences from the roll motion of the biped robot since the biped robot of this model is constrained on a plane. Therefore,
Fig. 6.
Magnitudes of the eigenvalues versus (β = 0.2, ω = π/0.1(l/g)1/2 [rad])
we will consider to construct a model which includes these influences and provide the analysis of these influences. ACKNOWLEDGMENT This paper is supported in part by Center of Excellence for Research and Education on Complex Functional Mechanical Systems (COE program of the Ministry of Education, Culture, Sports, Science and Technology, Japan). R EFERENCES [1] J. Adolfsson, H. Dankowicz, and A. Nordmark, 3D passive walkers: Finding periodic gaits in the presence of discontinuities, Nonlinear Dynamics 24, pp. 205-229, 2001. [2] M. Coleman, A. Chatterjee, and A. Ruina, Motions of a rimless spoked wheel: A simple three-dimensional system with impacts, Dynamics and Stability of Systems, Vol 12, No. 3, pp. 139-160, 1997. [3] R.J. Full and D.E. Koditschek, Templates and anchors: Neuromechanical hypotheses of legged locomotion on land, J. Exp. Biol. 202, pp. 3325-3332, 1999. [4] M. Garcia, A. Chatterjee, A. Ruina, and M. Coleman, The simplest walking model: Stability, complexity, and scaling, ASME J. of Biomechanical Eng., Vol. 120, No. 2, pp. 281-288, 1998. [5] R. Ghigliazza, R. Altendorfer, P. Holmes, and D.E. Koditschek, A simply stabilized running model, SIAM J. on Applied Dynamical Systems, Vol. 2, No. 2, pp. 187-218, 2003. [6] J.W. Grizzle, G. Abba, and F. Plestan, Asymptotically stable walking for biped robots: Analysis via systems with impulse effects, IEEE Trans. Automatic Control, Vol. 46, No. 1, pp. 51-64, 2001. [7] T. McGeer, Passive dynamic walking, Int. J. of Robotics Res., Vol. 9, No. 2, pp. 62-82, 1990. [8] S. Mochon and T. McMahon, Ballistic walking: An improved model, Math. biosc. 52, pp. 241-260, 1980. [9] A.L. Schwab and M. Wisse, Basin of attraction of the simplest walking model, Proc. of ASME Int. Conf. on Noise and Vibration, 2001. [10] G. Taga, A model of the neuro-musculo-skeletal system for human locomotion II. - Real-time adaptability under various constraints, Biol. Cybern. 73, pp. 113-121, 1995.
1370
