Optimal periodic locomotion for a two piece worm with an asymmetric
dry friction model

arXiv:1710.08546v1 [math.OC] 23 Oct 2017

Nak-seung Patrick Hyun1 and Erik Verriest2
Abstract— This paper solves the optimization problem for a
simplified one-dimensional worm model when the friction force
depends on the direction of the motion. The motion of the worm
is controlled by the actuator force f (t) which is assumed to be
piecewise continuous and always generates the same force in the
opposite directions. The paper derives the necessary condition
for the force which maximizes the average velocity or minimizes
the power over a unit distance. The maximum excursion of the
worm body and the force are bounded. A simulation is given
at the end of the paper.

In this paper, we first list the variables that are commonly
used in this paper, and define the model which includes the
friction model and the actuator model. Second, we suggest
two performance indices to be optimized and state the
problem formally. Third, we solve the problem analytically
and find the necessary conditions for the optimal solution.
Lastly, we provide one simulation example to understand the
solution better.
II. M ODELING

I. I NTRODUCTION
The locomotion for legless animals, such as a snake and
a worm, has many interesting features. One of the most
appealing characteristics of these types of animal is that
the periodic change in its shape generates the motion in a
certain direction. Scales on the animal’s skin create different
friction forces depending on the direction of the movement.
Due to this asymmetric friction, these types of animal,
such as snake or worm, can move forward. As a result,
many researchers have focused on analyzing the periodic
kinematics and dynamics of the legless animal locomotion in
an asymmetric friction model, Verriest [1] and Zimmermann
[2].
Even for a simplified one-dimensional worm model it is
complicated to analyze its full dynamics. In the past years, a
lot of research analyzed the one-dimensional toy problem in
depth, Zimmermann [3] and Chernousko [4]. Especially, the
book [3] has successfully characterized the kinematics and
dynamics of the one-dimensional worm problem.
Several researchers, Figurina [5] and Bolotnik [6], have
attempted to solve the optimization problem for a similar
one-dimensional worm problem. Chernousko [4] has solved
the optimal problem when the actuator force is modeled to be
piecewise constant and for the case when the infinite force is
available. He found the optimal control which maximizes the
average velocity. Chernousko [7] also optimizes the energy
consumed over a unit distance for all piecewise constant
forces. Here, we attempt to optimize a similar performance
index for an extended admissible control which includes all
the piecewise continuous forces. In addition, we pose the
realistic constraint that the force and the size of the worm
are bounded.
*This work was supported by a seed grant(RIM 3661304) from the Center
for Robotics and Intelligent Machines at Georgia Institute of Technology
1 School of Electrical and Computer Engineering Georgia Institute of
Technology Atlanta, Georgia 30332–0250 nhyun3@gatech.edu
2 School
of
Electrical
and
Computer
Engineering
Georgia
Institute
of
Technology
Atlanta,
Georgia
30332–0250

erik.verriest@ece.gatech.edu

M otion

M otion

m1

m2

fBW

f

fF W

x(t)
Fig. 1.

One dimensional two piece worm model.

A. Configuration space
Consider a simplified one-dimensional worm consisting
of two parts, head and tail. These parts are connected by
the actuator which exerts the same force f (t) at time t
to both parts (Fig.1). This actuator introduces energy to
system, which converts into kinetic energy for each parts and
compensates the work done by the friction force. Since the
worm contains two independent components, the dynamics
of motion can be written as two uncoupled second order
nonlinear differential equation. The nonlinearity comes from
the differential friction which depends on the direction of the
movement. More detail on the friction will be shown in the
next subsection.
−f (t) + f1 (t) = m1 x¨1 (t)
f (t) + f2 (t) = m2 x¨2 (t)

(1)

In this paper, we define a shape space, R+ , as an extend
of the worm, d(t). The shape can be changed by switching
the sign of the force f as the worm stretches or contracts
its body. With this shape space, the original space can be
prjected to the three dimensional space, where its trajectory
Θ = {(d(t), x˙1 (t), x˙2 (t)) ⊂ R3 |∀t ∈ R+ } lies in R3 ,
assuming that dynamics starts at t = 0 and appropriate initial
conditions are given. Using the change of basis method, we

now define the configuration space in R3 with its homeomorphic trajectory Ω = {(d(t), v(t), u(t)) ⊂ R3 | for ∀t ∈ R+ }.
This paper analyzes the trajectory, Ω of the worm.

B. Friction model
Our friction model is an asymmetric dry friction which
depends on the sign of the velocity but is independent
of the relative speed. The asymmetry in the friction will
result in moving the worm in a certain direction. More
concretely, let f be an asymmetric Coulomb friction which
is an approximate model for dry friction. In addition, let
m = m1 = m2 so that each part experiences the same
friction.

 −ff w
fbw
fi (t) =

f0

if ẋi (t) > 0
if ẋi (t) < 0
if ẋi (t) = 0

(2)


R t +T

v1 + t11 1 |f (t)|dt + α(t − t1 )


R t +T



v2 + t22 2 |f (t)|dt



 v + R t3 +T3 |f (t)|dt − α(t − t )
3
3
Rtt34 +T4
v(t) =

v4 − t4
|f (t)|dt − α(t − t4 )


R t5 +T5



|f (t)|dt
v
−
5

t5


 v − R t6 +T6 |f (t)|dt + α(t − t )
6
6
t6
u
+
β(t
−
t
)
CASE

1
1



u2 − ff w (t − t2 ) CASE



u3 + β(t − t3 ) CASE
u(t) =
u4 + β(t − t4 ) CASE




u

5 − ff w (t − t5 ) CASE


u6 + β(t − t6 ) CASE

CASE
CASE
CASE
CASE
CASE
CASE
1
2
3
4
5
6

1
2
3
(3)
4
5
6

(4)

Similarly, the dynamics of the configuration variable, d(t),
is solved in terms of v(t). With initial conditions d(ti ) = di
and v(ti ) = vi as before for each i case, d(t) can be derived
by integrating the Eqn(3)
Z ti +Ti
v(t)dt + di for any t ∈ [ti , ti + Ti ] (5)
d(t) = 2
ti

where ff w , fbw , f0 > 0 are constant and i = 1 or 2. The
friction is modeled to satisfy fbw ≥ ff w assuming that the
worm moves in forward direction. As a result of asymmetry,
the center of mass will move forward. In other words, the
third component in Ω: u(t), will always be nonnegative for
any t ∈ R+ .

C. Periodic motion in configuration space
There are eight cases of dynamics depending on the sign
of the velocities and the sign of the force acting on each
part. TABLE (I) shows the resulting eight possible cases
TABLE I
E IGHT MODES OF DYNAMICS .
CASE

1

2

3

4

5

6

7

8

f (t)

+

+

+

−

−

−

+

−

x˙1 (t)

+

+

−

−

+

+

−

−

x˙2 (t)

−

+

+

+

+

−

−

−

Duration

T1

T2

T3

T4

T5

T6

T7

T8

of dynamics. The table also contains the duration of each
case indicated by Ti for each case i. Observe that CASE 7
and CASE 8 cannot occur since u(t) is always nonnegative.
It suffices to analyze the motion for the six cases. For
conveninece, mi = 1 for all i. Using the friction model,
Eqn(2), and the dynamics, Eqn(1), we derive Eqn(3) and
Eqn(4). With initial conditions (v(ti ), u(ti )) = (vi , ui ) for
each i case, the Eqn(1) integrates to,

for all i. In order to find a periodic control which generates
a periodic motion in configuration space, it suffices to check
if there exists a periodic f defined as Eqn(1) such that
(d(t), v(t), u(t)) = (d(t + T ), v(t + T ), u(t + T )). If this is
true then the trajectory of (d(t), v(t), u(t)) follows an orbit
in Ω space and so does its projection to a subspace of Ω.
For example, a periodic orbit of the Ω trajectory is shown
in Fig.2, which (a) shows the projection to the (v, u) plane
and (b) shows the projection to the (v, d) plane. Here we
chose |f (t)| = 5 constant, (v(t1 ), u(t1 ), (d1 )) = (−3, 1, 10)
and ρ = 0.5.
D. Periodic actuator force
Most research on the periodic motion of the one dimensional worm model designed an actuator which generates
a simple harmonic force [3], [8] and analyzed its periodic
dynamics. However, in this paper, the worm is assumed to
generate an arbitrary periodic piecewise continuous force.
This actuator force is modeled as follows,
|f (t)| = F (t) t ∈ [t1 , t1 + T ]
f (t + T ) = f (t) for all t ∈ R+

(6)

where T is the period and F (t) is a positive piecewise
continuous function. If F (t) ≥ fbw , the actuator can change
the direction of motion of each part. Previous research [9],
[10] shows a sufficient condition for the existence of periodic
motion in certain nonlinear systems. Since the friction force
and periodic actuator control are bounded, the nonlinearities
in Eqn(1) are bounded by maxt∈[t1 ,t1 +T ] {F (t)} + fbw . By
the construction, both forces are continuous almost everywhere and also we have a freedom in choosing the sign
of the control force. This satisfies the sufficient condition
derived in [10], and so there exists a periodic solution for
Eqn(1).

uHtL
1.2

integral term in Eqn(3) and (5) will cancel out and so
v(t1 + T ) = v(t1 )
d(t1 + T ) = d(t1 ).

1.0
(v_0,u_0)

0.8
0.6
u(t)= -v(t)

u(t)=-v(t)

0.4
0.2

-3

-2

0

-1

1

2

3

vHtL

(a) (v, u) plane

u(t1 + T )

vHtL
3
2

9.0

9.5

10.0

10.5

11.0

11.5

12.0

dHtL

-1
(v_0,u_0)

(b) (d, v) plane
Simulation results for ρ = 0.5, d0 = 10 and F = 5

Here, we list four assumption for the actuator force design.
First, we assume that the worm switches the sign of the
actuator force once within a period. It is required to switch
the sign of the actuator force at least once to pull back to
the original state in the configuration space. This assumption
shows that it tends to only extend its body as much as
possible when f (t) > 0, and contract its body back as
the sign changes to negative, f (t) < 0. In the presence of
inertia, we assume that the switch happens when two parts
move towards or move in the opposite direction in order to
keep the maximum worm size. The physical constraint of
the maximum allowable extension will appear in the next
section.
Second, following the design assumption, we enumerate a
sequence of cases in TABLE (I) which the worm experiences.
Suppose the force is initially set to positive at time t1 as
defined before, then Case1 → Case2 → Case3 will be the
motion sequence for f (t) > 0 in t ∈ [t1 , t1 + T1 + T2 + T3 ].
After the sign switches at t1 + T1 + T2 + T3 , the sequence
changes to Case4 → Case5 → Case6 for f (t) <P
0 in t ∈
i−1
[t1 +T1 +T2 +T3 , t1 +T ]. This shows that ti = t1 + k=1 Tk
for all i = 2, · · · , 6.
Third, there exist a symmetry on the actuator force design
between the positive and negative region. In other words,
the initial state in (v, u)-plane for negative force is equal to
(v4 , u4 ) = (−v1 , u1 ). By using the definition of (v, u), this
shows that (x˙1 (t4 ), x˙2 (t4 )) = (x˙2 (t1 ), x˙1 (t1 )). Therefore,
by designing the negative force f (t) = −f (t − T2 ), the
dynamics in Eqn(1) are the same except for the negative sign
on each equation. This guarantees that such a force design
will generate the periodic motion in (v, d)-space since the

T 1+ρ
2 ( 1−ρ )
T 2ρ
2 ( 1+ρ ) − T1

T2 =
T3 =

-3

Fig. 2.

= u1 + β(T1 + T3 + T4 + T6 ) − ff w (T2 + T5 )
= u1 + 2β(T1 + T3 ) − 2ff (T2 )
= u1 .

In addition, we now have a useful set of equations.

1

-2

Lastly, we define a relation between T1 , T2 and T3 in
order to have u(t1 ) = u(t1 + T ) which finalize the full
periodic motion in the configuration space. From the above
third assumption, we know that T1 = T4 , T2 = T5 , T3 = T6
and so T2 = T1 +T2 +T3 = T4 +T5 +T6 . Here, we claim that
if T2 = 1−ρ
2ρ (T1 + T3 ) then u(t1 + T ) = u(t1 ). By following
the sequence defined above and Eqn(4), we have

(7)
(8)

Hence, for given T and T1 , any piecewise continuous function F (t) for t ∈ [t1 , t1 + T2 ] can generate a periodic motion
with period T in the configuration space.
E. Physical constraints
There are two physical constraints for the worm in this
paper. One is the size of the worm and the other is the
upper limit of the actuator force. The first one states that
the difference between maximum and minimum excursion
of the worm’s body is limited to E(T ) ≤ L and the actuator
was designed to reach this bound, L. The other constraint is
that there exists an upper limit, fu > 0, of the actuator force
F because of the physical limit of the worm’s muscle. With
these two constraints, we can derive F which optimizes given
performances such as the power over unit distance traveled
and the speed of the center of the mass. Here, we define a
set of admissible control forces by
ST = {F ∈ P C[t1 , t1 + T2 ] : F (t) ∈ [fbw , fu ]
for all t ∈ [t1 , t1 + T2 ]}

(9)

where P C is a set of piecewise continuous function.
III. O PTIMIZATION
A. Performance index
Research on the optimization of this one dimensional
worm model was first suggested in [4], and extended on
[7]. The latter paper suggests two performance indices which
we adopt in this paper. One is the power over the unit
distance and the other is the average velocity over one period.
Other papers in the optimization for one-dimensional worm
problem, [5] [6], also focused on maximizing the average
velocity over one period.
Let Eqn(10) be the total distance that the center of mass
traveled and the average power exerted by the actuator,
respectively.
R t +T
X = t11 u(t)dt
R t +T
(10)
P = T1 t11 F (t)v(t)dt

The two performance indeces in this paper are as follows,
V =
Pu =

X
T
P
X.

(11)
(12)

These performance indices are the average velocity over one
period, V , and the average power over unit distance, Pu .
B. Optimal control problem
Past research of the optimal control problem on one
dimensional two body system have modeled the actuator
force to be constant in each dynamic motion, [4] and [7].
In this paper, we not only include the constant model but
also generalize to any piecewise continuous force model in
each dynamics.
There are two optimal control problems in this paper
for the fixed period T . The problem is stated for each
performance index,
1) Find an optimal F ∈ ST for a given T1 with the
maximum excursion constraint, E(T ) = L.
2) Find an optimal F ∈ ST with the maximum excursion
limit, E(T ) = L.
C. Solution to the first problem
Pick T1 such that T1 ∈ [0, T2 − T2 ], then do the following:
1) STEP 1: Compute the boundary condition for each
dynamics and E(T )
2) STEP 2: Compute two performance index, Pu and V
in Eqn(12) and Eqn(11).
3) STEP 3: Solve the optimization problem for a given
Tmin .
4) STEP 4: Find the necessary and sufficient conditions
for the existence of solution.
5) STEP 5: Find the optimal Tmin which minimize the
performance.
Here are the details:
1) STEP 1: First, we compute the boundary condition for
each dynamics. To do this, we need to define three additional
functions, H : [t1 , t1 + T1 ] → R, G : [t2 , t2 + T2 ] → R and
I : [t3 , t3 + T3 ] → R, which satisfy
Z t
H(t) =
F (t)dt
(13)
t1
Z t
G(t) =
F (t)dt
(14)
t2
Z t
I(t) =
F (t)dt.
(15)
t3

The sequence of motion and the TABLE I shows that
x˙2 (t2 ) = 0 and x˙1 (t3 ) = 0. By the dynamics of xi (t)
for i ∈ {1, 2}, we derive the boundary conditions for H
and G. In addition, by using the third assumption of the
actuator force, v4 = −v1 , we derive the boundary condition
for I as well. Here, we use Eqn(3), Eqn(7) and Eqn(8) for

computation.
R t +T
= t11 1 F (t)dt + fbw T1 + (u1 + v1 )
= H(t1 + T1 ) + fbw T1 + (u1 + v1 )
=0
R t +T
x˙1 (t3 ) = t22 2 −F (t)dt − ff w T2 + x˙1 (t2 )
= −G(t2 + T2 ) − ff w T2 + (2u1 + 2βT1 )
=0
R t +T
v̇(t4 ) = t33 3 F (t)dt − αT3 + βT1 − ff w T2 + u1
T ff w
= I(t3 + T3 ) − fbw T1 − 1+ρ
+ u1
= −v1 .
x˙2 (t2 )

By using Eqn(7) and the above equations, the boundary
conditions for H, G and I are
H(t1 + T1 ) = −fbw T1 + c1

(16)

G(t2 + T2 ) = 2βT1 + c2

(17)

I(t3 + T3 ) = −fbw T1 + c3

(18)

T ff w 1−ρ
2 1+ρ

where c1 = −(u1 + v1 ), c2 = 2u1 −
and c3 =
T ff w
−(u1 + v1 ) + 1+ρ
are constant.
Second, we need to compute the difference between maximum and minimum excursion, E(T ). To do this, we need
to find the time when the size of the worm, d(t), are at
minimum and maximum. Since v(t) < 0 for Case1, the
Eqn(5) is decreasing in [t1 , t1 + T1 ]. However, we know that
v(t2 + T2 ) > 0 and so there exist some Tmin ∈ [0, T2 ] such
that v(t2 + Tmin ) = 0. Similarly, v(t) > 0 for Case4 and
so d(t) is increasing in [t4 , t4 + T4 ] but v(t5 + T5 ) < 0.
Therefore, there also exist Tmax ∈ [0, T5 ] = [0, T2 ] such
that v(t5 + Tmax ) = 0. The symmetry assumption for the
actuator force shows that Tmin = Tmax . From this we can
compute E(T ) as follows.
R t +Tmax
v(t)dt
E(T ) = 2 t25+Tmin
R t4
R t +T
= 2( t2 +Tmin v(t)dt + t45 max v(t)dt)
Rt
R t +T
= 2( t24+Tmin v(t)dt − t12 min v(t)dt)
R t2 +Tmin
Rt
= 2(− t1
v(t)dt + t24+Tmin v(t)dt).
By applying Eqn(3) to the above equation, we get
Rt
Rt
E(T ) = 2(−( t12 H(t)dt + t12 (αt + v1 )dt)
R t +T
R t +T
−( t22 min G(t)dt + t22 min (−αT1 − u1 )dt)
R t3
R t3
+( t2 +Tmin G(t)dt + t2 +Tmin (−αT1 − u1 )dt)
Rt
Rt
+( t34 I(t)dt + t34 (−αt + fbw T1 + βT2 + u1 )dt))
= 2h(T1 , u1 , v1 , Tmin )
Rt
Rt
+2( t34 I(t)dt − t12 H(t)dt)
R t +T
Rt
−2( t22 min G(t)dt − t23+Tmin G(t)dt)
where 2h(T1 , u1 , v1 , Tmin ) collects the constant terms of the
equation. By the constraint of this optimization problem, we
want E(T ) = L and so finally we get,
R t4
Rt
I(t)dt − t12 H(t)dt = L2 − h(T1 , u1 , v1 , Tmin )
t3
R t +T
(19)
+ t22 min G(t)dt
Rt
− t23+Tmin G(t)dt

In addition, v(t2 + Tmin ) = 0 shows that
βT1 − u1 = 0 and so we have

R t2 +T
t2

F (t)dt −

G(t2 + Tmin ) = βT1 + u1 .

(20)

2) STEP 2: Now, we compute the performance index.
We start with the average velocity over one period defined
in Eqn(11). To do this, we compute X first. By using the
Eqn(4), Eqn(10), Eqn(7) and Eqn(8), we get

X=

R t1 +T
t1

u(t)dt

Rt
= 2 t14 u(t)dt
= T (u1 + βT1 −

T ff w β
4α ).

T ff w β
4α

(21)

Next, we compute the second performance index, Pu .
Since X was constant for fixed T1 , it is enough to optimize
the power over a unit distance. The symmetry assumption in
the actuator force design shows that the total work, WT,T1 , is
equal to two times of the work done in the half of a period,
2WT /2,T1 . Therefore, it suffices to compute the work done
in the half of a period.

WT /2,T1

R t +T
= t11 2 F (t)v(t)dt
= W1 + W2 + W3

(22)

R t +T
where Wi = tii i F (t)v(t)dt are the work done in the
i−th Case for i ∈ {1, 2, 3}. By using the fundamental
theorem of calculus, we can show that H 0 = F for a.e
t ∈ (t1 , t1 + T1 ), G0 = F for a.e t ∈ (t2 , t2 + T2 ) and
I 0 = F for a.e t ∈ (t3 , t3 + T3 ), where a.e stands for almost
everywhere. This shows that
R t1 +T1
t1

R t2 +T2
t2

F (t)H(t)dt
F (t)G(t)dt

R t3 +T3
t3

F (t)I(t)dt

=

R H(t1 +T1 )

2

1)
(23)
= H(t1 +T
+ H(t1 + T1 )(αT1 + v1 )
R t12+T1
−α t1
H(t)dt
R t +T
= W1const − α t11 1 H(t)dt
R t +T
W2 = t22 2 F (t)v(t)dt
R t +T
2
2)
= G(t2 +T
− (βT1 + u1 ) t22 2 F (t)dt
2
(24)
2
2)
= G(t2 +T
−
(βT
+
u
)G(t
+
T
)
1
1
2
2
2
= W2const
R t +T
W3 = t33 3 F (t)v(t)dt
2

By simply dividing by T , we get the average velocity
performance index,

V = u1 + βT1 −

simplify the equations,
R t +T
W1 = t11 1 F (t)v(t)dt
R t +T
2
1)
= H(t1 +T
+ t11 1 (F (t)αt + v1 F (t))dt
2

HdH

H(0)
2
1)
= H(t1 +T
2
R G(t +T )
= G(0)1 1 GdG
2
2)
= G(t2 +T
2
R I(t +T )
= I(0)1 1 IdI
2
3)
= I(t3 +T
.
2

By using the above equations and Eqn(3), we can compute
each Wi (Ti ) for all i. Here, we use integration by parts to

3)
− (I(t3 + T3 ) + v1 )I(t3 + T3 )
= I(t3 +T
R t32+T3
(25)
+α t3
I(t)dt
R t +T
= W3const + α t33 3 I(t)dt.

Let Wconst (T1 , T, u1 , v1 ) = W1const + W2const + W3const
and by substituting above equations to Eqn(22), we get
WT,T1

= 2W T ,T1
2
= 2Wconst (T1 , T, u1 , v1 )
(26)
R t +T
R t +T
+2α( t33 3 I(t)dt − t11 1 H(t)dt).

By substituting the constraint, Eqn(19), we conclude that
WT,T1

= 2Wconst (T1 , T, u1 , v1 )
+αL − 2αh(T1 , u1 , v1 , Tmin )
R t +T
+2α t22 min G(t)dt
Rt
−2α t23+Tm in G(t)dt.

(27)

3) STEP 3: Since V did not depend on the choice
of F ∈ ST,T1 for fixed T1 , there is nothing to optimize. For the second performance index, we summarize
the problem as follows. For given T1 and L, we want
to find F ∈ ST,T1 such that minimize the total work,
WT,T1 in Eqn(27), and satisfies the boundary conditions in
Eqn(16), Eqn(17), Eqn(18) and the excursion constraint in
Eqn(19). Observe that the only function that changes the
total work is G since other term remain constant for a
given T1 and Tmin . Therefore, it is free to choose H and I
within the admissible control that satisfies all the constraints.
Hence, the optimization problem reducesRto find the optimal
t +T
G which optimizes minG∈P C[t2 ,t2 +T2 ] ( t22 min G(t)dt −
R t3
G(t)dt) with the given boundary constraint.
t2 +Tmin
We regard G(t) as a composition with two different
dynamics

F (t) if t ∈ [t2 , t2 + Tmin ]
G(t) =
−F (t) if t ∈ (t2 + Tmin , t2 + T2 ].
Accordingly, we construct two Hamiltonians,
H 1 (t) = G(t) + λ1 (t)F (t)
H 2 (t) = G(t) − λ2 (t)F (t).
Since F (t) ∈ [fbw , fu ] bounded, by applying the Pontryagain
minimum principle( [11]), we find the optimal control

fbw if λ1 (t) > 0
F1 (t) =
fu if λ1 (t) < 0

for H 1 Hamiltonian. Computing the Euler-Lagrange equation on H 1 (t) gives the differential equation for λ1 (t)
equation.
∂H 1
λ˙1 (t) = −
= −1
∂G
which is solved by λ1 (t) = −t + τ1 for some τ1 constant.
By using the internal point constraint, Eqn(20), we get
R t +T
G(t2 + Tmin ) = t22 min F (t)dt
= fbw τ1 + (Tmin − τ1 )fu
= βT1 + u1
and solving for τ1 gives
τ1

=
=

Tmin fu −(βT1 +u1 )
fu −fbw
(βT +u )
fu − T 1 1
min
Tmin .
fu −fbw

Similarly, we solve for the second Hamiltonian system. By
applying the Pontryagin minimum principle, we find the
optimal control

fu if λ2 (t) > 0
F2 (t) =
fbw if λ2 (t) < 0
and the Euler-Lagrange equation gives
2

∂H
= −1
λ˙2 (t) = −
∂G
which is solved by λ2 (t) = −t + τ2 for some τ2 constant.
Now, by using the final state constraint, Eqn(17), we compute
R t +T
G(t2 + T2 ) = t22 min F (t)dt + G(t2 + Tmin )
= fu τ2 + (T2 − (Tmin + τ2 ))fbw + βT1 + u1
= 2βT1 + c2
and solving for τ2 gives
=
=

Z

t1
t2 +Tmin

t1

Z

t2 +Tmin

fbw dt ≤ G(t2 + Tmin ) ≤

fu dt

t2

t2

Z

t2 +T2

fbw dt ≤ G(t2 + T2 ) − G(t2 + Tmin )
t2 +Tmin
t2 +T2

Z

+u1 )
Since τ1 ∈ [0, Tmin ], if Tmin ≤ (βTf1bw
, the optimal
solution exist and the solution is

fbw if t ∈ [t2 , t2 + τ1 ]
F1 (t) =
fu if t ∈ [t2 + τ1 , t2 + Tmin ].

τ2

4) STEP 4: Here, we list the necessary conditions on
(u1 , v1 , Tmin ) which gives the existence of a solution. Since
F (t) is a positive bounded function and H, G and I are
increasing functions, the following bounds need to hold
Z t1 +T1
Z t1 +T1
fbw dt ≤ H(t1 + T1 ) ≤
fu dt

Tmin fb −2αT2 +βT1 +u1
fu −fbw
−αT2 +βT1 +u1
−fbw
T2 −Tmin
(T2
fu −fbw

− Tmin ).

1 −u1
Since τ2 ∈ [0, T2 − Tmin ], if Tmin ≤ T2 + αT2 −βT
, the
fu
optimal solution exist and the solution is

fu if t ∈ [t2 + Tmin , t2 + τ2 ]
F2 (t) =
fbw if t ∈ [t2 + τ2 , t2 + T2 ]

Hence, by combining the results of F1 and F2 , the optimal
actuator force in [t2 , t2 + T2 ] which minimizes the total
energy consumed for a single period is

 fbw if t ∈ [t2 , t2 + τ1 ]
fu if t ∈ [t2 + τ1 , t2 + Tmin + τ2 ]
F (t) =
(28)

fbw if t ∈ [t2 + τ2 , t2 + T2 ].
In the next section, the necessary condition for existence of
such a solution will be covered.

fu dt ≥ G(t2 + T2 ) − G(t2 + Tmin )
t2 +Tmin
Z t3 +T3

Z

t3 +T3

fbw dt ≤ I(t3 + T3 ) ≤
t3

fu dt.
t3

Using the conditions for Tmin in STEP 3, the above bounds
simplify to
T (η + ρ)
T
− T1 ) ≥ u1 ≥ β( − T1 )
2(1 + ρ)
2
−Kfbw T1 − u1 ≤ v1 ≤ −2fbw T1 − u1
max{γ, T2 (1 + ρ) − γη} ≤ Tmin
Tmin ≤ min{γη, T2 (1 + ηρ ) − γ}
ρ
0 < T1 < T
1+ρ
1≤η
β(

(29)
(30)
(31)
(32)
(33)

1)
where γ = (βT1f+u
, K = min{1 + η, (η − 1) TT31 − 1} and
u
fu
η = fbw ≥ 1. By choosing T1 satisfying Eqn(32), we found
a set of admissible initial conditions, (v1 , u1 ), which satisfy
Eqn(29) and Eqn(30). The bounds for Tmin are then well
defined by the Eqn(29) and Eqn(30).
5) STEP 5: By substituting the optimal G(t) into Eqn(27),
we can compute the power over unit distance from Eqn(12).
A long computation shows that Eqn(12) ends up with a
quadratic equation of Tmin . By using the boundary condition
for Tmin , Eqn(31), we can find the optimal Tmin . Finally, by
choosing F and I that satisfies Eqn(19), the first optimization
problem is solved.

D. Solution to the second problem
1) Maximum average velocity: By using the boundary
condition for T1 in Eqn(32) and the average velocity in
Eqn(21), we can compute T1 which maximizes the average
velocity. Since β > 0 for all admissible T1 , we conclude that
ρ
V reaches its maximum when T1 = T 1+ρ
.
2) Power over unit distance: This problem is solved by
STEP 5 of the first problem. Here, instead of Tmin that
minimizes the performance, we find the optimal (T1 , Tmin )
pair for given boundary conditions. The boundary conditions
for this problem are obtained from Eqn(32) and Eqn(31).
IV. S IMULATIONS
For simulation purpose, we choose the physical constraints
as T = 10, fbw = 1, ff w = 0.1, fu = 5, d1 = 40, t1 = 0 and
L = 32.261. Let uratio and vratio be the relative portion

from the minimum to the maximum of its boundary, Eqn(29)
and Eqn(30), respectively. Assume that uratio = 0.2 and
vratio = 0.5. Now define T1r as the relative portion between
minimum and maximum of the T1 boundary. Similarly,
define Tminr as the relative portion between minimum and
maximum of the Tmin boundary. Fig.3 shows the plot of the
power over unit distance for all possible pairs (T1r , Tminr ).

xHtL coordinate

60

40

20

time
2

Fig. 5.

4

6

8

10

A simulation for the motion in one period.

and after the size of the worm gets close to its minimum or
maximum length.
VI. D ISCUSSION
Fig. 3.

A simulation for the power over an unit distance.

The minimum occurs when T1r = 0.363635 and Tminr =
0.563214. For a given (T1 , Tmin ), we can find the optimal
actuator force model for G by using Eqn(28). Since we have
freedom to choose H and I which satisfies the maximum
excursion constraint, Eqn(19), one of the solutions that meet
the maximum excursion constraint is shown in Fig.4.

Since the set, {t ∈ [t2 , t2 + T2 ] : λi (t) = 0}, only
contains two points, there does not exist a singular solution
which performs better than the local optimal solution. In
addition, all the admissible controls F in Case1 and Case3,
which satisfy the equality constraint Eqn(19), generate the
same power performance. If we consider the problem when
there exist an additional penalty before and after the swtich
happens, then one may find the unique optimal control
among the solutions to Eqn(19).

ControlHf HtLL

R EFERENCES
4
2
time
2

4

6

8

10

-2
-4

Fig. 4. The optimal actuator force model that minimize the power over
unit distance.

The minimum and maximum of excursion occur when t =
2.24 and t = 7.74, respectively. Fig.5 shows the motion of
each part. The blue(top) line is the trajectory of x2 (t) and the
red(bottom) line is the trajectory of x1 (t). The middle line
shows the motion of the center of mass. The optimal force
uses its maximum in an interval about the point of minimum
and maximum excursion.
V. C ONCLUSION
Whereas previous work researched the admissible controls
to be piecewise constant without a bound on the force, we
have shown in this paper that such a shape of control can be
also the locally optimal for the case when we optimize over
the piecewise continuous and bounded controls. However,
due to the freedom in choosing H and I, other solutions
having the same performance exist. The result shows that it
is necessary to allocate the maximum allowable force before

[1] E. I. Verriest and D. Yeung, “Locomotion based on differential
friction,” in Proc. of the American Control Conference, 2008. IEEE,
2008, pp. 862–867.
[2] K. Zimmermann, I. Zeidis, H. Jianjun, and J. Steigenberger, “An
approach to the modelling of worm-like motion systems with finite
degree of freedom–first steps in technical realization,” PAMM, vol. 2,
no. 1, pp. 529–530, 2003.
[3] C. Behn, I. Zeidis, and K. Zimmermann, Mechanics of Terrestrial
Locomotion [electronic resource] : With a Focus on Non-pedal Motion
Systems / by Carsten Behn, Igor Zeidis, Klaus Zimmermann. Berlin,
Heidelberg : Springer-Verlag Berlin Heidelberg, 2009., 2009.
[4] F. Chernous’Ko, “The optimum rectilinear motion of a two-mass
system,” Journal of applied Mathematics and Mechanics, vol. 66,
no. 1, pp. 1–7, 2002.
[5] T. Figurina, “Optimal motion control for a system of two bodies on
a straight line.” Journal of Computer and Systems Sciences International, vol. 46, no. 2, pp. 227–233, 2007.
[6] N. Bolotnik and T. Figurina, “Optimal control of periodic motions
of vibration-driven systems.” in IFAC Proceedings Volumes (IFACPapersOnline), vol. 3, no. PART 1, Institute for Problems in Mechanics, Russian Academy of Sciences, 2007, pp. 142–147.
[7] F. Chernous’Ko, “Analysis and optimization of the rectilinear motion
of a two-body system.” Journal of Applied Mathematics and Mechanics, vol. 75, no. 5, pp. 493–500, 2011.
[8] H.-b. Fang and J. Xu, “Dynamics of a three-module vibration-driven
system with non-symmetric coulombs dry friction,” Multibody System
Dynamics, vol. 27, no. 4, pp. 455–485, 2012.
[9] R. Miller and A. Michel, “On existence of periodic motions in
nonlinear control systems with periodic inputs,” SIAM Journal on
Control and Optimization, vol. 18, no. 5, pp. 585–598, 1980.
[10] P. Nistri, “Periodic control problems for a class of nonlinear periodic
differential systems,” Nonlinear Analysis: Theory, Methods & Applications, vol. 7, no. 1, pp. 79–90, 1983.
[11] A. E. Bryson and Y. C. Ho, Applied Optimal Control. New York:
Blaisdell, 1969.

