Stability Analysis and Control Design for an Underactuated Walking Robot via Computation of a Transverse Linearization Leonid Freidovich1 , Anton Shiriaev1,2 , Ian Manchester1 1 Department

of Applied Physics and Electronics Umeå University, Sweden (http : //www.control.tfe.umu.se/).

2 Department of Engineering Cybernetics Norwegian University of Science and Technology, Norway (http : //www.itk.ntnu.no/english/).

IFAC’08 WC, July 9, Seoul, South Korea

Dynamics of a Three-Link Walking Robot Design of an Orbitally Stabilizing Controller Analysis Based on Virtual Constraints Approach Summary and General Remarks

Outline

1

Dynamics of a Three-Link Walking Robot

2

Design of an Orbitally Stabilizing Controller

3

Analysis Based on Virtual Constraints Approach

4

Summary and General Remarks

Leonid Freidovich, Anton Shiriaev, Ian Manchester

Stability Analysis and Control Design for an Underactuated Walkin

Dynamics of a Three-Link Walking Robot Design of an Orbitally Stabilizing Controller Analysis Based on Virtual Constraints Approach Summary and General Remarks

Outline

1

Dynamics of a Three-Link Walking Robot

2

Design of an Orbitally Stabilizing Controller

3

Analysis Based on Virtual Constraints Approach

4

Summary and General Remarks

Leonid Freidovich, Anton Shiriaev, Ian Manchester

Stability Analysis and Control Design for an Underactuated Walkin

A Three-Link Walking Robot

The generalized coordinates are the absolute angles: θ1 , θ2 , and θ3 . The picture and the model are taken from [Grizzle, Abba, Plestan’99/01]. r = 1 – length of each leg, m = 5 – mass of each leg, MH = 15 – mass of hips, MT = 10 – mass of torso, L = 0.5 – distance between hips and torso.

Dynamics in the Swinging Phase The Lagrangian dynamics are described by  5 25 ¨ θ − cos(θ1 − θ2 )θ¨2 + 2 cos(θ3 − θ1 )θ¨3 2 2 1

 −2 sin(θ3 − θ1 )θ˙32 − sin(θ1 − θ2 )θ˙22 − 13g sin θ1 = −u1 ,   5 θ¨2 2 + g sin θ ¨ ˙ − cos(θ − θ ) θ + sin(θ − θ ) θ 2 = −u2 , 1 2 1 1 2 1 2 2   ¨ 5 cos(θ3 −θ1 )θ¨1 + θ23 + sin(θ3 −θ1 )θ˙12 − g sin θ3 = u1 +u2 ,

where u1 and u2 are the controlled torques applied between the legs and the torso.

It is assumed that the swinging leg, which corresponds to θ2 , moves to be set in front of the stance leg (described by θ1 ), which acts as a pivot until it reaches a certain constant value.

Dynamics in the Swinging Phase (Cont’d) It is convenient to use simultaneously the following notation ¨ + C(q, q) ˙ q˙ + G(q) = B u, D(q) q using the standard notation: for the inertia matrix D(q), the ˙ the vector of matrix of Coriolis and centrifugal forces C(q, q), generalized gravitationalhforces G(q), and with iT q = [θ1 , θ2 , θ3 ]T , q˙ = θ˙1 , θ˙2 , θ˙3 , u = [u1 , u2 ]T , and


x˙ = f x(t), u , h iT with x = θ1 , θ2 , θ3 , θ˙1 , θ˙2 , θ˙3 .

The swinging (single support) phase is assumed to hold until θ1 riches the value of π/8 ≈ 0.3927.

Dynamics in the Swinging Phase (Cont’d) It is convenient to use simultaneously the following notation ¨ + C(q, q) ˙ q˙ + G(q) = B u, D(q) q using the standard notation: for the inertia matrix D(q), the ˙ the vector of matrix of Coriolis and centrifugal forces C(q, q), generalized gravitationalhforces G(q), and with iT q = [θ1 , θ2 , θ3 ]T , q˙ = θ˙1 , θ˙2 , θ˙3 , u = [u1 , u2 ]T , and


x˙ = f x(t), u , h iT with x = θ1 , θ2 , θ3 , θ˙1 , θ˙2 , θ˙3 .

The swinging (single support) phase is assumed to hold until θ1 riches the value of π/8 ≈ 0.3927.

Dynamics in the Swinging Phase (Cont’d) It is convenient to use simultaneously the following notation ¨ + C(q, q) ˙ q˙ + G(q) = B u, D(q) q using the standard notation: for the inertia matrix D(q), the ˙ the vector of matrix of Coriolis and centrifugal forces C(q, q), generalized gravitationalhforces G(q), and with iT q = [θ1 , θ2 , θ3 ]T , q˙ = θ˙1 , θ˙2 , θ˙3 , u = [u1 , u2 ]T , and


x˙ = f x(t), u , h iT with x = θ1 , θ2 , θ3 , θ˙1 , θ˙2 , θ˙3 .

The swinging (single support) phase is assumed to hold until θ1 riches the value of π/8 ≈ 0.3927.

The Impact Model An impact with the ground occurs when the solution of the continuous-time dynamics hits the smooth surface o iT n h Γ− = x = θ1 , θ2 , θ3 , θ˙1 , θ˙2 , θ˙3 ∈ R6 : θ1 = π/8 . After that, the state vector is instantaneously mapped by F : Γ− ∋ x(t−) 7−→ x(t+ ) ∈ Γ+ = F (Γ−) defined by   ˙ ˙ θi (t+ ) = θi (t−) and θi (t+ ) = ωi q(t−), q(t−) ,

where i ∈ {1, 2, 3}, the arguments t− and t+ denote the values right before and right after the impact, correspondingly,

The Impact Model An impact with the ground occurs when the solution of the continuous-time dynamics hits the smooth surface o iT n h Γ− = x = θ1 , θ2 , θ3 , θ˙1 , θ˙2 , θ˙3 ∈ R6 : θ1 = π/8 . After that, the state vector is instantaneously mapped by F : Γ− ∋ x(t−) 7−→ x(t+ ) ∈ Γ+ = F (Γ−) defined by   ˙ ˙ θi (t+ ) = θi (t−) and θi (t+ ) = ωi q(t−), q(t−) ,

where i ∈ {1, 2, 3}, the arguments t− and t+ denote the values right before and right after the impact, correspondingly,

The Impact Model (Cont’d)

 ˙ = 5 θ˙1 − 20θ˙1 cos(2θ1 − 2θ2 ) + 4θ˙1 cos(2θ3 − 2θ1 ) ω1 (q, q)  +2θ˙2 cos(2θ1 − 2θ2 ) /∆(q),  ˙ = 10 2θ˙1 cos(2θ3 −θ1 −θ2 ) ω2 (q, q)  −9θ˙1 cos(θ1 −θ2 ) + θ˙2 /∆(q),  ˙ = 5 12θ˙1 cos(θ1 +θ3 −2θ2 ) − 12θ˙1 cos(θ1 − θ3 ) ω3 (q, q) +θ˙1 cos(3θ1 −2θ2 −θ3 ) − θ˙2 cos(θ2 −θ3 ) − 19 θ˙ 2 3 +θ˙3 cos(2θ1 −2θ2 ) + 2θ˙3 cos(2θ2 −2θ3 ) /∆(q), ∆(q) = −95 + 10 cos(2 θ1 − 2 θ2 ) + 20 cos(2 θ2 − 2 θ3 ).

The Total Hybrid Model After the impact the next step is described by a similar continuous-time dynamics, so that we have: ˙ ¯ x= f1 (·)

F1

Γ1+ −→ Γ1− → {z } | step one x(t0+ ) ˙ x(t) t1 x(t1+ ) ˙ x(t) t2 x(t2+ ) t0

∈ =

def

=

= = def

=

= def

=

˙ ¯ x= f2 (·)

F2

Γ2+ −→ Γ2− → Γ3+ ≡ Γ1+ → . . . . | {z } {z } | ... step two

Γ1+ ,
¯ f1 x(t), u

for t0 < t < t1 ,

arg min {t1 > t0 : x(t1−) ∈ 6 Γ1−} ,
F1 x(t1−) ∈ Γ2+ ,
¯ for t1 < t < t2 , f2 x(t), u

arg min {t2 > t1 : x(t2−) 6∈ Γ2−} ,
F2 x(t2−) ∈ Γ1+ , t2

(redefine and restart),

The Total Hybrid Model After the impact the next step is described by a similar continuous-time dynamics, so that we have: ˙ ¯ x= f1 (·)

F1

Γ1+ −→ Γ1− → {z } | step one x(t0+ ) ˙ x(t) t1 x(t1+ ) ˙ x(t) t2 x(t2+ ) t0

∈ =

def

=

= = def

=

= def

=

˙ ¯ x= f2 (·)

F2

Γ2+ −→ Γ2− → Γ3+ ≡ Γ1+ → . . . . | {z } {z } | ... step two

Γ1+ ,
¯ f1 x(t), u

for t0 < t < t1 ,

arg min {t1 > t0 : x(t1−) ∈ 6 Γ1−} ,
F1 x(t1−) ∈ Γ2+ ,
¯ for t1 < t < t2 , f2 x(t), u

arg min {t2 > t1 : x(t2−) 6∈ Γ2−} ,
F2 x(t2−) ∈ Γ1+ , t2

(redefine and restart),

The Total Hybrid Model (Cont’d)

Due to the symmetry of the mass-length distribution we have:  Γ1− = Γ2+ = Γ− = x ∈ R6 : θ1 = π/8 ,  Γ1+ = Γ2− = Γ+ = F (Γ−) = x ∈ R6 : θ2 = π/8 ,



¯ F1 x = F x , f1 x(t), u = f x(t), u ,



¯ F2 x = F P x , f2 x(t), u = P −1 f P x(t), u , where the 6 × 6 permutation matrix P denotes the linear transformation renaming the legs (θ1 ↔ θ2 and θ˙1 ↔ θ˙2 ).

The goal of the control design is to create (plan and orbitally stabilize) a periodic two-step motion, reminiscent to walking.

The Total Hybrid Model (Cont’d)

Due to the symmetry of the mass-length distribution we have:  Γ1− = Γ2+ = Γ− = x ∈ R6 : θ1 = π/8 ,  Γ1+ = Γ2− = Γ+ = F (Γ−) = x ∈ R6 : θ2 = π/8 ,



¯ F1 x = F x , f1 x(t), u = f x(t), u ,



¯ F2 x = F P x , f2 x(t), u = P −1 f P x(t), u , where the 6 × 6 permutation matrix P denotes the linear transformation renaming the legs (θ1 ↔ θ2 and θ˙1 ↔ θ˙2 ).

The goal of the control design is to create (plan and orbitally stabilize) a periodic two-step motion, reminiscent to walking.

Dynamics of a Three-Link Walking Robot Design of an Orbitally Stabilizing Controller Analysis Based on Virtual Constraints Approach Summary and General Remarks

Outline

1

Dynamics of a Three-Link Walking Robot

2

Design of an Orbitally Stabilizing Controller

3

Analysis Based on Virtual Constraints Approach

4

Summary and General Remarks

Leonid Freidovich, Anton Shiriaev, Ian Manchester

Stability Analysis and Control Design for an Underactuated Walkin

Control Design: Steps 1 and 2 [Step 1.] Define the outputs to be regulated
y1c = θ2 −py1c (θ1 ) ≡ θ2 +θ1 − a5 +a6 θ1 +a7 θ12 +a8 θ13 ,

2 y2c = θ3 −py2c (θ1 ) ≡ θ3 − a1 +a2 θ1 +a3 θ12 +a4 θ13 θ12 −θ1d with

a1 = 0.512,

a2 = 0.073,

a3 = 0.035,

a4 = −0.819,

a5 = −2.27,

a6 = 3.26,

a7 = 3.11,

a8 = 1.89,

[Step 2.] Perform the partial linearizing feedback transformation with respect to these outputs −1     0 v 3×2 1c u = Jc (θ1 , θ˙1 ) −1 v2c D (q)B   ˙ q ˙
, −Jc (θ1 , θ1 ) ˙ q˙ − G(q) D −1 (q) −C(q, q)

Control Design: Steps 1 and 2 [Step 1.] Define the outputs to be regulated
y1c = θ2 −py1c (θ1 ) ≡ θ2 +θ1 − a5 +a6 θ1 +a7 θ12 +a8 θ13 ,

2 y2c = θ3 −py2c (θ1 ) ≡ θ3 − a1 +a2 θ1 +a3 θ12 +a4 θ13 θ12 −θ1d with

a1 = 0.512,

a2 = 0.073,

a3 = 0.035,

a4 = −0.819,

a5 = −2.27,

a6 = 3.26,

a7 = 3.11,

a8 = 1.89,

[Step 2.] Perform the partial linearizing feedback transformation with respect to these outputs −1     0 v 3×2 1c u = Jc (θ1 , θ˙1 ) −1 v2c D (q)B   ˙ q ˙
, −Jc (θ1 , θ1 ) ˙ q˙ − G(q) D −1 (q) −C(q, q)

Control Design: Steps 2 (Cont’d) and Step 3 [Step 2 (cont’d).] The matrix Jc (θ1 , θ˙1 ) =

"

′′ (θ ) θ˙ , 0, 0, −p ′ (θ ), 1, 0 −py1c 1 1 y1c 1 ′′ (θ ) θ˙ , 0, 0, −p ′ (θ ), 0, 1 −py2c 1 1 y2c 1 ˙

#

∂[yc , yc ] is the Jacobian matrix ∂[θ for yc = [y1c , y2c ]T and 1 , θ2 , θ3 ] y˙ c = [y˙ 1c , y˙ 2c ]T . With this transformation the continuous-time dynamics can be rewritten as

θ¨1 = . . .

¨2c = v2c , y

¨1c = v1c , y

[Step 3.] Apply the continuous finite-time stabilizing feedback 9

11

9

vic = −10 10 sign(y˙ ic ) |y˙ ic | 10 − 102 sign(φic ) |φic | 11 , 9

φic = yic +

10 10 110

11

sign (y˙ ic ) (|y˙ ic |) 10 ,

i = 1, 2.

Control Design: Steps 2 (Cont’d) and Step 3 [Step 2 (cont’d).] The matrix Jc (θ1 , θ˙1 ) =

"

′′ (θ ) θ˙ , 0, 0, −p ′ (θ ), 1, 0 −py1c 1 1 y1c 1 ′′ (θ ) θ˙ , 0, 0, −p ′ (θ ), 0, 1 −py2c 1 1 y2c 1 ˙

#

∂[yc , yc ] is the Jacobian matrix ∂[θ for yc = [y1c , y2c ]T and 1 , θ2 , θ3 ] y˙ c = [y˙ 1c , y˙ 2c ]T . With this transformation the continuous-time dynamics can be rewritten as

θ¨1 = . . .

¨2c = v2c , y

¨1c = v1c , y

[Step 3.] Apply the continuous finite-time stabilizing feedback 9

11

9

vic = −10 10 sign(y˙ ic ) |y˙ ic | 10 − 102 sign(φic ) |φic | 11 , 9

φic = yic +

10 10 110

11

sign (y˙ ic ) (|y˙ ic |) 10 ,

i = 1, 2.

Stability result [Grizzle, Abba, Plestan’01] It is shown that the controller exponentially orbitally stabilize a periodic motion in the closed-loop system. The achieved motion starts with the continuous-time “step one” dynamics initiated at  x(t0+ ) = x0 ≈ −0.392699, 0.392699, 0.496465, T 0.926751, −0.239267, 1.483783 ∈ Γ1+

and reaches

 x((t0 + T )−) = xT ≈ 0.392699, −0.392699, 0.496466, T 1.549891, −2.031622, −0.431507 ∈ Γ1−

in the half-period T ≈ 1.12163.

After that, the impact occurs followed by an identically controlled continuous-time “step two” dynamics and then by the final jump due to impact.

Are there alternatives? Proof of orbital stability is based on the fact that the controller brings the trajectories to a two-dimensional manifold in a sufficiently short finite time. The region of attraction is quite small. Let us redesign the controller using a simpler to implement linear approximation of the feedback after linearization: vic = −13.023 yic − 106.45 (yic + 0.0699 y˙ ic ) The stability proofs based on technique of Grizzle, Abba, and Plestan or generalizations do not work. Let us analyze the new closed-loop system using analytically computed transverse linearization. Such an analytical construction can be done using knowledge of the virtual constraint for the continuous part of the cycle, which is the geometrical synchronization law.

Are there alternatives? Proof of orbital stability is based on the fact that the controller brings the trajectories to a two-dimensional manifold in a sufficiently short finite time. The region of attraction is quite small. Let us redesign the controller using a simpler to implement linear approximation of the feedback after linearization: vic = −13.023 yic − 106.45 (yic + 0.0699 y˙ ic ) The stability proofs based on technique of Grizzle, Abba, and Plestan or generalizations do not work. Let us analyze the new closed-loop system using analytically computed transverse linearization. Such an analytical construction can be done using knowledge of the virtual constraint for the continuous part of the cycle, which is the geometrical synchronization law.

Poincaré map

S

Figure: A surface S transversal to the flow (locally, around the intersection with the cycle) is called a Poincaré section. Poincaré map P : S → S is defined by the first hit rule; the orbit corresponds to a fixed point of P(·).

Linearization of a Poincaré map

S

TS

Figure: Linearization dP of the Poincaré map P(·) is defined on the tangent plane TS. Eigenvalues of dP say on exponential stability (instability) of the orbit.

Moving Poincaré section and transverse dynamics S(0) TS(0) S(t) TS(t)

Figure: Moving Poincaré section is a family S(t)t∈[0,T ) of transversal surfaces. It allows to introduce transversal coordinates (on each surface) and a system ζ˙ = A(t) ζ , with A(t) = A(t + T ), where ζ(t) are linear increments of the transversal coordinates that leave on TS(t).

Dynamics of a Three-Link Walking Robot Design of an Orbitally Stabilizing Controller Analysis Based on Virtual Constraints Approach Summary and General Remarks

Outline

1

Dynamics of a Three-Link Walking Robot

2

Design of an Orbitally Stabilizing Controller

3

Analysis Based on Virtual Constraints Approach

4

Summary and General Remarks

Leonid Freidovich, Anton Shiriaev, Ian Manchester

Stability Analysis and Control Design for an Underactuated Walkin

Computing the impact-invariant outputs Since in the closed-loop system there exists an exponentially orbitally stable trajectory, θ1 = θ⋆(t),

θ2 = θ2⋆(t),

θ1 = θ3⋆(t)

there must exist two functions: py1 (θ) and py2 (θ), defining an induced virtual holonomic constraint, such that the two outputs y1 = θ2 − py1 (θ1 )

and y2 = θ3 − py2 (θ1 )

are identically equal to zeros along the planned trajectory of the closed-loop hybrid system: θ2⋆(t) ≡ py1 (θ⋆(t)) and θ3⋆(t) ≡ py2 (θ⋆(t)).

Numerically computed virtual constraint 0.4

0.6 0.4 py1(θ)

0.3 0.2

0.2 0 −0.2

θ ⋆ (t)

0.1 −0.4 −0.4

−0.2

0 θ

0.2

0.4

−0.2

0 θ

0.2

0.4

0 0.56 −0.1 0.54 py2(θ)

−0.2 −0.3 −0.4

0.52 0.5

0

0.5 t

1

0.48 −0.4

Figure: (a) θ⋆ (t) – the first component of the desired trajectory; (b) p1y (θ) and p1y (θ) – the approximations for the functions, defining the imposed constraint.

Computing the virtual limit system The function θ⋆(t) is a solution of the reduced dynamics α(θ) θ¨ + β(θ) θ˙ 2 + γ(θ) = 0, where

′ (θ) cos θ − p (θ) α(θ) = − 25 1 + py1 y1

′ (θ) + 1 cos θ − p (θ) +5 py2 y2
′ (θ) + 2p ′ (θ) , + 45 25 + py1 y2

5 ′′ ′′ (θ) cos θ − p (θ) β(θ) = − 2 py1 (θ) cos θ − py1 (θ) + 5py2 y2   2
′ ′′ (θ) py2 (θ) − 1 sin θ − py2 (θ) + 54 py1 +5   2
′′ (θ) + 5 ′ (θ) + 25 py2 1 − p sin θ − py1 (θ) , y1 2  

sin p (θ) − 2 sin p (θ) − 13 sin(θ) . γ(θ) = 5g y1 y2 2

Conserved quantity and its properties The reduced dynamics has the conserved quantity Z θ
2 γ(s) ˙ = θ˙ 2 − ψ θ⋆(0), θ θ˙⋆(0)2 + ds, ψ(s, θ) I(θ, θ) α(s) o θ⋆ (0) n R s β(s) ds ψ(s1 , s2 ) = exp − s12 2α(s) It is equivalence to the r Euclidean distance 2  2  ˙ = k (·) · min I(θ, θ) θ − θ⋆(τ ) + θ˙ − θ˙⋆(τ ) τ

Satisfies the passivity-like relation   d 2 2 β(θ) ˙ I =θ W− I dt α(θ) α(θ) for

α(θ)θ¨ + β(θ)θ˙ 2 + γ(θ) = W

Admits the sensitivity derivatives ∂I ∂I = −2 θ¨⋆(t), ∂θ θ=θ⋆ (t), ∂ θ˙ ˙ θ˙⋆ (t) θ=

θ=θ⋆ (t), ˙ θ˙⋆ (t) θ=

= 2 θ˙⋆(t)

Conserved quantity and its properties The reduced dynamics has the conserved quantity Z θ
2 γ(s) ˙ = θ˙ 2 − ψ θ⋆(0), θ θ˙⋆(0)2 + ds, ψ(s, θ) I(θ, θ) α(s) o θ⋆ (0) n R s β(s) ds ψ(s1 , s2 ) = exp − s12 2α(s) It is equivalence to the r Euclidean distance 2  2  ˙ = k (·) · min I(θ, θ) θ − θ⋆(τ ) + θ˙ − θ˙⋆(τ ) τ

Satisfies the passivity-like relation   d 2 2 β(θ) ˙ I =θ W− I dt α(θ) α(θ) for

α(θ)θ¨ + β(θ)θ˙ 2 + γ(θ) = W

Admits the sensitivity derivatives ∂I ∂I = −2 θ¨⋆(t), ∂θ θ=θ⋆ (t), ∂ θ˙ ˙ θ˙⋆ (t) θ=

θ=θ⋆ (t), ˙ θ˙⋆ (t) θ=

= 2 θ˙⋆(t)

Computing a continuous transverse linearization The five coordinates describing the dynamics transverse to the desired orbit: during the first step (ζ1 ) during the second step (ζ2 )     1 1 I θ1 , θ˙1 , θ⋆(0), θ˙⋆(0) I θ2 , θ˙2 , θ⋆(0), θ˙⋆(0) 2

y1 = θ2 − py1 (θ1 )

2

y1 = θ1 − py1 (θ2 )

3

y2 = θ3 − py2 (θ1 ) ′ (θ ) θ˙ y˙ 1 = θ˙2 − py1 1 1

3

y2 = θ3 − py2 (θ2 ) ′ (θ ) θ˙ y˙ 1 = θ˙1 − py1 2 2

4 5

′ (θ ) θ˙ y˙ 2 = θ˙3 − py2 1 1

4 5

′ (θ ) θ˙ y˙ 2 = θ˙3 − py2 2 2

Each linear system is in the form   ζ˙i = A(t) + B(t) K(t) ζi

i = 1, 2

and with analytically computed functions A(t) and B(t) as functions of θ⋆(t), θ˙⋆(t), and θ¨⋆(t).

Computing linearization for the impact map First of all, we need to compute the Jacobians of the impact map ∂F (·) (dF ) = ˙ [q, q] at the end points of the continuous-time arcs of the cycle. In addition, we need to derive 1

˙ restricted to the The transformations between [q, q] tangent planes to the flow, which coincides with the tangent planes to moving Poincaré sections, and the transversal coordinates ζ.

2

The projection operator along the flow onto the tangent planes to the switching surfaces.

Computing linearization for the impact map First of all, we need to compute the Jacobians of the impact map ∂F (·) (dF ) = ˙ [q, q] at the end points of the continuous-time arcs of the cycle. In addition, we need to derive 1

˙ restricted to the The transformations between [q, q] tangent planes to the flow, which coincides with the tangent planes to moving Poincaré sections, and the transversal coordinates ζ.

2

The projection operator along the flow onto the tangent planes to the switching surfaces.

Computing linearization for the impact map First of all, we need to compute the Jacobians of the impact map ∂F (·) (dF ) = ˙ [q, q] at the end points of the continuous-time arcs of the cycle. In addition, we need to derive 1

˙ restricted to the The transformations between [q, q] tangent planes to the flow, which coincides with the tangent planes to moving Poincaré sections, and the transversal coordinates ζ.

2

The projection operator along the flow onto the tangent planes to the switching surfaces.

Hybrid transverse linearization Finally, for our system ˙ ¯ x= f1 (·)

˙ ¯ x= f2 (·)

F1

Γ1+ −→ Γ1− → {z } | step one

F2

Γ2+ −→ Γ2− → Γ3+ ≡ Γ1+ → . . . . | {z } {z } | ... step two

we have the following linearization ¯1+ Γ |




ζ˙1 = A(t)+B(t) K(t) ζ1

−→

¯2+ Γ |

¯1− Γ

{z step one




ζ˙2 = A(t)+B(t) K(t) ζ2

−→

{z step two

−1 −

MT+ (dF) MT

→

¯2− Γ

}

−1 −

MT+ (dF) MT

→

}

¯3+ ≡ Γ ¯ →.... Γ | 1+ {z } ...

If the zero solution of this linearization is asymptotically stable, then the periodic solution in the original closed-loop system is exponentially orbitally stable.

Analyze of the hybrid transverse linearization The obtained transverse linearization allows us to compute the transition matrix as follows: Solve the initial value problem
˙ Φ(t) = A(t) + B(t) K(t) Φ(t), φ(0) = I5 Compute the matrix   Φ(T ) Ψ(T ) = MT+ (dF ) MT−1 −

The maximal absolute value of the eigenvalues of the transition matrix Ψ(T ) is approximately equal to 0.36377, which is an estimate of the stability degree in a vicinity of the desired trajectory during every half of the period.

Feedback controller based on the linearization

The feedback controller from [Grizzle, Abba, Plestan’01] is locally equivalent to the composition of an appropriate nonlinear feedback transformation followed by a linear stabilizing feedback (for the step one)  h iT τ = K T (θ1 , θ˙1 ) I(θ1 , θ˙1 ), y1 , y2 , y˙ 1 , y˙ 2

where the operator T (·) recovers the time-stamp of the projection onto the cycle, which can be taken as   θ − θ⋆(0+ ) ˙ T (θ, θ) = T mod 1 . θ⋆(T−) − θ⋆(0+ )

Numerically simulations with the linearized controller 0.5 θ1 (t)

2.5

0 −0.5

0

2

4 t, sec

6

8

0

2

4 t, sec

6

8

0

2

4 t, sec

6

8

1

1.5

0.5

θ2 (t)

kζ(t)k = k[I, y1 , y˙ 1 , y2 , y˙ 2 ] k

2

1

0 −0.5 0.8

θ3 (t)

0.5

0

0

5 t, sec

10

0.6 0.4

Figure: Signals in the hybrid closed-loop system with (1) the original controller, (2) the new approximation-based controller, (3) the linearization-based controller for the approximation are impossible to distinguish.

Redesigning controller Since the stability is determined by o n Φ(T )} J (K(t)) = max1≤i≤n eigi {MT+ (dF ) MT−1