Stability conditions of an ODE arising in human motion and its numerical simulation

TL;DR

This paper analyzes the stability of an ODE model from human musculoskeletal control, providing a sufficient stability condition and validating it through numerical simulations and experiments.

Contribution

It introduces a new stability condition for a musculoskeletal system model and demonstrates its effectiveness with simulations and experimental data.

Findings

A sufficient condition for asymptotic stability is derived.

Numerical simulations confirm the stability condition.

Experimental results support the theoretical findings.

Abstract

This paper discusses the stability of an equilibrium point of an ordinary differential equation (ODE) arising from a feed-forward position control for a musculoskeletal system. The studied system has a link, a joint and two muscles with routing points. The motion convergence of the system strongly depends on the muscular arrangement of the musculoskeletal system. In this paper, a sufficient condition for asymptotic stability is obtained. Furthermore, numerical simulations of the penalized ODE and experimental results are described.

Tables

Table 1. Table 1: Experiment parameters

$L_0$	285.0 [mm]
$L_1$	110.0 [mm]
$b_1$	87.0 [mm]
$b_2$	5.0 [mm]
${\ell }_1$, ${\ell }_2$	99.0 [mm]
$r_1$, $r_2$	35.0 [mm]
$s_1$, $s_2$	35.0 [mm]

Equations

$$I\ddot{\theta }(t)+\mu \dot{\theta }(t)-\tau (\theta (t))=0,t>0,$$

$$\begin{array}{rl}& q_{11}(\theta )={\left\{{\left(L_0-b_1-{\ell }_1\sin \frac{\theta }{2}\right)}^2+{\left(d_1-{\ell }_1\cos \frac{\theta }{2}\right)}^2\right\}}^{1/2},\\ & q_{12}(\theta )={\left\{{\left(r_1-{\ell }_1\sin \frac{\theta }{2}\right)}^2+{\left(s_1-{\ell }_1\cos \frac{\theta }{2}\right)}^2\right\}}^{1/2},\\ & q_{21}(\theta )={\left\{{\left(L_0-b_2+{\ell }_2\cos \left(\pi +\frac{\pi +\theta }{2}\right)\right)}^2+{\left(d_2+{\ell }_2\sin \left(\pi +\frac{\pi +\theta }{2}\right)\right)}^2\right\}}^{1/2},\\ & q_{22}(\theta )={\left\{{\left(r_2+{\ell }_2\cos \frac{\pi -\theta }{2}\right)}^2+{\left(s_2-{\ell }_2\sin \frac{\pi -\theta }{2}\right)}^2\right\}}^{1/2},\end{array}$$

$$\tau (\theta )=-⟨J(\theta ),\left(\begin{array}{c}{F}_1\\ F_2\end{array}\right)⟩.$$

$${\mathbb{R}}^2\ni F:=\left(\begin{array}{c}{F}_1\\ F_2\end{array}\right)=-\tau (\theta )J(\theta {)}^{\ast }+v(\theta ),$$

$$J(\theta {)}^{\ast }:=\{\begin{array}{ll}\frac{J(\theta )}{|J(\theta ){|}^2}& \text{mbox}ifJ(\theta )\ne 0,\\ 0& \text{mbox}ifJ(\theta )=0.\end{array}$$

$$v(\theta )=k\left(\begin{array}{c}\frac{d{q}_2}{d\theta }(\theta )\\ -\frac{d{q}_1}{d\theta }(\theta )\end{array}\right)\text{mbox}foranyk\in \mathbb{R}.$$

$$\tau (\theta )=-⟨\frac{dq}{d\theta }(\theta ),v({\theta }_d)⟩=-⟨\left(\begin{array}{c}\frac{d{q}_1}{d\theta }(\theta )\\ \frac{d{q}_2}{d\theta }(\theta )\end{array}\right),k\left(\begin{array}{c}\frac{d{q}_2}{d\theta }({\theta }_d)\\ -\frac{d{q}_1}{d\theta }({\theta }_d)\end{array}\right)⟩$$

$$I\ddot{\theta }+\mu \dot{\theta }+⟨\frac{dq}{d\theta }(\theta ),v({\theta }_d)⟩=0,t>0.$$

$$\dot{x}(t)=g(x(t)),t>0.$$

$$L_0(x,\xi ):=K(\xi )-P(x).$$

$$\frac{{\partial }^{2}K}{\partial {\xi }^2}(\dot{x}(t))\ddot{x}(t)+\frac{\partial P}{\partial x}(x(t))=Q(t,x(t),\dot{x}(t))\text{mbox}fort>0,$$

$$⟨Q(t,x,\xi ),\xi ⟩\le 0\text{mbox}fort>0,x\in {\mathbb{R}}^n,\xi \in {\mathbb{R}}^n.$$

$$P(x_d)=0\text{mbox}andP(x)>0\text{mbox}forx\in B_R(x_d)\setminus \{x_d\}.$$

$$⟨Q(t,x,\xi ),\xi ⟩<0\text{mbox}fort>0,(x,\xi )\in {\mathbb{R}}^n\times {\mathbb{R}}^n\setminus \{(x_d,0)\},$$

$$V(x,\xi ):=⟨\frac{\partial K}{\partial \xi }(\xi ),\xi ⟩-K(\xi )+P(x).$$

$$\{\begin{array}{ll}& \frac{{\partial }^{2}K}{\partial {\xi }^2}(\xi (t))\dot{\xi }(t)+\frac{\partial P}{\partial x}(x(t))=Q(t,x(t),\xi (t)),\\ & \xi (t)=\dot{x}(t),\end{array}t>0.$$

$$\frac{d}{dt}(V(x(t),\xi (t)))=⟨\frac{{\partial }^{2}K}{\partial {\xi }^2}(\xi )\dot{\xi },x⟩+⟨\frac{\partial P}{\partial x},\xi ⟩=⟨Q,\xi ⟩\le 0.$$

$$\frac{d}{dt}(V(x(t),\xi (t)))<0\text{mbox}fort>0,(x,\xi )\in {\mathbb{R}}^n\times {\mathbb{R}}^n\setminus \{(x_d,0)\},$$

$$\begin{array}{rl}f(\theta )=f(\theta ;i,j)& :={\left\{{\left(a_{ij}+b_{ij}\sin \frac{\theta }{2}\right)}^2+{\left(c_{ij}+b_{ij}\cos \frac{\theta }{2}\right)}^2\right\}}^{1/2}\\ & ={\left\{a_{ij}^2+b_{ij}^2+c_{ij}^2+2b_{ij}\sqrt{a_{ij}^2+c_{ij}^2}\sin \left(\frac{\theta }{2}+{\alpha }_{ij}\right)\right\}}^{1/2}\end{array}$$

$$\sin {\alpha }_{ij}=\frac{c_{ij}}{\sqrt{a_{ij}^2+c_{ij}^2}},\cos {\alpha }_{ij}=\frac{a_{ij}}{\sqrt{a_{ij}^2+c_{ij}^2}}.$$

$$\displaystyle\begin{split}\begin{array}[]{llr}&a_{11}=L_{0}-b_{1},\ b_{11}=-\ell_{1},\ c_{11}=d_{1}&\mbox{for }i=1,\ j=1;\\ &a_{12}=r_{1},\ b_{12}=-\ell_{1},\ c_{12}=s_{1}&\mbox{for }i=1,\ j=2;\\ &a_{21}=-(L_{0}-b_{2}),\ b_{21}=-\ell_{2},\ c_{21}=d_{2}&\mbox{for }i=2,\ j=1;\\ &a_{22}=-r_{2},\ b_{22}=-\ell_{2},\ c_{22}=s_{2}&\mbox{for }i=2,\ j=2.\end{array}\end{split}$$

$$b_{ij}<0\text{mbox}and{c}_{ij}>0\text{mbox}fori,j=1,2.$$

$$\sqrt{a_{ij}^2+c_{ij}^2}\ne |b_{ij}|.$$

$$a_{1j}>0\text{mbox}and{a}_{2j}<0\text{mbox}forj=1,2.$$

$${\alpha }_{11},{\alpha }_{12}\in (0,\pi /2)\text{mbox}and{\alpha }_{21},{\alpha }_{22}\in (\pi /2,\pi ).$$

$$C_{{\theta }_0}=C_{{\theta }_0,ij}:=\sqrt{a_{ij}^2+c_{ij}^2}\sin \left(-\frac{{\theta }_0}{2}+{\alpha }_{ij}\right).$$

$$0<{\theta }_0<\min \{2\min \{{\alpha }_{21},{\alpha }_{22}\}-\pi ,\pi /4\}.$$

$$\frac{c_{2j}}{\sqrt{a_{2j}^2+c_{2j}^2}}=\sin {\alpha }_{2j}<\underset{\theta \in (-{\theta }_0,\pi )}{\sup }\sin \left(\frac{\theta }{2}+{\alpha }_{2j}\right)=\sin \left(-\frac{{\theta }_0}{2}+{\alpha }_{2j}\right)=\frac{C_{{\theta }_0,2j}}{\sqrt{a_{2j}^2+c_{2j}^2}}<1.$$

$$-{\theta }_0<{\theta }_d<\pi -2\max \{{\alpha }_{11},{\alpha }_{12}\}.$$

$$f^{\prime }({\theta }_d)=\frac{b_{ij}\sqrt{a_{ij}^2+c_{ij}^2}}{2}f({\theta }_d{)}^{-1}\cos \left(\frac{{\theta }_d}{2}+{\alpha }_{ij}\right).$$

Full text

Key Words and Phrases: Stability condition, Musculoskeletal system with routing points, Numerical simulation, Experimental result

1. 1 Introduction

2. 2 Preliminaries

3. 3 Stability condition of a 1-link-2-muscle musculoskeletal system with routing points

4. 4 Proof of Theorem 5

5. 5 Numerical and experimental results

6. 5.1 Numerical simulation

7. 5.2 Experimental results

1 Introduction

The mechanism of human motion is expected to be applied to robotics. Some hypotheses such as “the equilibrium point hypothesis” [3, 4] and “the virtual trajectory hypothesis” [5] suggest that human motion generation efficiently utilizes a feed-forward position control. In addition, rapid motion without sensory feedback, such as a finger flick, can be controlled to some extent. Therefore, it is assumed that a human musculoskeletal system satisfies the controllability condition for a feed-forward control. We would like to interpret that assumption in both a mathematical and engineering sense.

The approach of our feed-forward control is to keep muscular tensions balanced at a target position. In 2013, Kino et al. [6] studied the feed-forward control for a 2-link-6-muscle musculoskeletal system, which is modeled after a human arm. They gave an engineering consideration and showed their numerical simulations. We also refer to Kino et al. [7] for a mathematically sufficient condition of the same system as [6] for feed-forward control. However, the system is considered without some characteristics. For instance, each muscle of the system is arranged as straight, even though complicated arrangements actually exist around the joints.

In this paper, we consider a musculoskeletal system (Figure 1) that has a link, a joint and two muscles. This system is modeled after a human finger. The name of each part is defined according to Figure 1. We also use the same symbols for the lengths. Base $L_{0}$ is always fixed. Link $L_{1}$ rotates on the joint in the two-dimensional plane. The red line and the blue line in Figure 1 are muscles and are called Muscle 1 and Muscle 2, respectively. These muscles have routing points $P_{12}$ and $P_{22}$. Generally, human muscles are restricted by the tendons and curve with the rotation of the links. For simplification, in this system, we assume that Muscle 1 and Muscle 2 bend at routing points $P_{12}$ and $P_{22}$ (Figure 3), which are the lengths $\ell_{1}$ and $\ell_{2}$ away from the joint, respectively. We call $\ell_{1}$ and $\ell_{2}$ virtual links (Figure 3), although the system primarily has the unique link, Link $L_{1}$. We let the virtual links $\ell_{1}$ and $\ell_{2}$ rotate by half of the rotation angle of Link $L_{1}$. We suppose that friction of the joint comprises only viscosity friction with a viscosity friction coefficient of $\mu>0$ and that the system ignores a viscoelasticity of muscles. Muscles are imaged to resemble wires. We also ignore gravity effects. In this case, we are concerned with a sufficient condition for the above feed-forward control of the system.

To this end, as a mathematical problem, this paper discusses a sufficient condition for asymptotic stability of a equilibrium point of the dynamics

[TABLE]

where $\theta:(0,\infty)\to(-\theta_{0},\pi)$ is an unknown function. Here, a constant $\theta_{0}>0$ is small enough. A constant $I$, a constant $\mu$ and a function $\tau:(-\theta_{0},\pi)\to\mathbb{R}$ are given. $I$ is the moment of inertia, $\mu$ is the viscosity coefficient of the joint, and $\tau$ is the torque of the link generated by constant muscular tensions balancing at the target position $\theta=\theta_{d}$.

We give a derivation of $\tau$. Let $F_{1}$ and $F_{2}$ be constant muscular tensions of Muscle 1 and Muscle 2, respectively. Let $q_{ij}=q_{ij}(\theta):\mathbb{R}\to\mathbb{R}$ ($i,j=1,2$) be the length of $P_{ij}P_{i(j+1)}$ as shown in Figure 1. They are given by

[TABLE]

where parameters $L_{0},b_{1},b_{2},d_{1},d_{2},\ell_{1},\ell_{2},r_{1},r_{2},s_{1}$ and $s_{2}$ are positive constants. Here, $q_{i}:=q_{i1}+q_{i2}$ denotes the length of Muscle $i$ for $i=1,2$. By the principle of virtual work, we have

[TABLE]

Here, $\langle\cdot,\cdot\rangle:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ denotes the Euclidean inner product. It follows that

[TABLE]

where $q:=(q_{1},q_{2})$, $J(\theta):=\dfrac{dq}{d\theta}(\theta)$ and

[TABLE]

The second term $v(\theta)$ of the right-hand side belongs to $\Set{v\in\mathbb{R}^{2}}{\langle J(\theta),v\rangle=0}$, namely,

[TABLE]

Here, the first term $-\tau(\theta)J(\theta)^{*}$ denotes a force to rotate the link $L_{1}$, that is, a driving force at $\theta$. The second term $v(\theta)$ denotes an orthogonal force to the rotation of $L_{1}$, that is, an inner force at $\theta$. For instance, $F=v(\theta_{d})$ is an internal force vector of Muscle 1 and Muscle 2 balancing at $\theta=\theta_{d}$. Since our feed-forward control is to take $F\equiv v(\theta_{d})$, the torque $\tau(\theta)$ is given by

[TABLE]

for some $k\in\mathbb{R}$. The dynamics of (1.1) becomes the ordinary differential equation

[TABLE]

We note that the solution, $\theta=\theta_{d}$, is an equilibrium point of (1.2).

Our aim in this paper is to give conditions for parameters $L_{0}$, $b_{1}$, $b_{2}$, $d_{1}$, $d_{2}$, $\ell_{1}$, $\ell_{2}$, $r_{1}$, $r_{2}$, $s_{1}$ and $s_{2}$, and an equilibrium point $\theta_{d}$, for the asymptotic stability of $\theta=\theta_{d}$. Since each muscle of the 2-link-6-muscle musculoskeletal system in [7] is straight, they show a sufficient condition by applying the Taylor expansion of the muscular lengths at some angle. In this paper, we use a different method because of a complicated $q_{i}$ by routing points.

This paper is organized in the following way. In Section 2, we recall some known results. Sections 3 and 4 are devoted to showing a sufficient condition. In Section 5, we show the numerical simulation results and experimental results.

2 Preliminaries

We first recall Lyapunov’s stability theorem, which corresponds to that of Theorem 1.30 in [2]. We also refer to [1] for an application to the control of mechanical systems.

Proposition 1** (Theorem 1.30 in [2]).**

Let $x_{0}$ be an equilibrium point of the autonomous ordinary differential equation

[TABLE]

Let continuous function $V:U\to\mathbb{R}$ be a Lyapunov function for (2.1) at $x_{0}$, i.e., $V$ implies the following:

• •

$V(x_{0})=0$;

• •

$V(x)>0$* for $x\in U\setminus\{x_{0}\}$;*

• •

the function $x\mapsto\mathrm{grad}\,V(x)$ is continuous for $x\in U\setminus\{x_{0}\}$, and on this set, $\dfrac{d}{dt}(V(x(t)))\leq 0$,

where $U\subset\mathbb{R}^{n}$ is an open set. Then, $x_{0}$ is Lyapunov stable. In addition, if $V$ is a strict Lyapunov function, i.e., $\dfrac{d}{dt}(V(x(t)))<0$ for $x\in U\setminus\{x_{0}\}$, then $x_{0}$ is asymptotically stable.

Define a Lagrangian $L_{0}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ by

[TABLE]

Here, $K:\mathbb{R}^{n}\to\mathbb{R}$ and $P:\mathbb{R}^{n}\to\mathbb{R}$ are respectively a kinetic energy and a potential energy. We consider the Lagrange equation

[TABLE]

where $Q:(0,\infty)\times\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a generalized force. We always assume that $K$, $P$ and $Q$ are smooth functions.

We next recall an application of Lyapunov’s stability theorem to (2.2). The following proposition provides a sufficient condition for Lyapunov/asymptotic stability of an equilibrium of (2.2). We write the proof for researchers in other fields, although the proof is elementary.

Proposition 2**.**

Let $x=x_{d}\in\mathbb{R}^{n}$ be an equilibrium point of (2.2). Let $K$ be convex such that $K(0)=0$. Let $Q$ satisfy

[TABLE]

Assuming that $P$ is locally positive definite around $x_{d}$, i.e., there exists $R>0$ such that

[TABLE]

Then, $x=x_{d}$ is Lyapunov stable. Furthermore, if

[TABLE]

then $x=x_{d}$ is an asymptotically stable equilibrium point.

Proof 2.1**.**

We set

[TABLE]

We show that $V$ becomes a Lyapunov function for

[TABLE]

at an equilibrium point, $(x,\xi)=(x_{d},0)$. Since $K$ is convex and (2.4) holds, $V$ is locally positive definite around $(x_{d},0)$.

By (2.2) and (2.3), we have

[TABLE]

Thanks to Proposition 1, $x_{d}$ is Lyapunov stable since $V$ is a Lyapunov function for (2.6).

In addition, we assume (2.5). Similarly, it follows that

[TABLE]

which implies that $x_{d}$ is asymptotically stable.

3 Stability condition of a 1-link-2-muscle musculoskeletal system with routing points

In this section, we are concerned with a sufficient condition for the stability of the equilibrium point, $\theta=\theta_{d}$, of the equation (1.2).

Define a function $f=f(\cdot;i,j):(-\theta_{0},\pi)\to\mathbb{R}$ by

[TABLE]

for constants $a_{ij},b_{ij},c_{ij}\in\mathbb{R}$, where $\alpha_{ij}\in[0,2\pi)$ satisfies

[TABLE]

In what follows, the constants $a_{ij},b_{ij},c_{ij}$ are always given by

[TABLE]

We note that $f(\theta;i,j)$ becomes $q_{ij}$ ($i,j=1,2$). We compute $f=f(\cdot;i,j)$ instead of $q_{ij}$ to obtain a sufficient condition for the stability. We note that

[TABLE]

We present a list of hypotheses on $a_{ij},b_{ij},c_{ij}$ and $\theta_{0}$. We assume the following for $f>0$:

[TABLE]

This assumption means that $P_{ij}\neq P_{i(j+1)}$ ($i,j=1,2$) as shown in Figure 1.

Next, we suppose that

[TABLE]

Thus, we have

[TABLE]

Finally, we give an assumption of $\theta_{0}>0$. We set

[TABLE]

To have a comparison between $C_{\theta_{0}}$, $a$ and $c$, we assume that

[TABLE]

In fact, (3.7) implies

[TABLE]

We assume that $k\dfrac{dq_{2}}{d\theta}(\theta_{d})$ and $-k\dfrac{dq_{1}}{d\theta}(\theta_{d})$ in $v(\theta_{d})$ are muscular tensions. Therefore, we require that they are positive. A sufficient condition for that requirement is as follows.

Proposition 3**.**

Assume that (3.4), (3.5) and (3.7) hold. We have $\dfrac{dq_{2}}{d\theta}(\theta_{d})>0$ for $\theta_{d}\in(-\theta_{0},\pi)$. In addition, we assume

[TABLE]

Then, $k\dfrac{dq_{2}}{d\theta}(\theta_{d})$ and $-k\dfrac{dq_{1}}{d\theta}(\theta_{d})$ are positive for any $k>0$.

Remark 4**.**

For a typical 1-link-2-muscle muscloskeletal system without routing points, balanced muscular tensions are always positive at any target angle.

Proof 3.1** (Proof of Proposition 3).**

From (3.1), we have

[TABLE]

By (3.6), it is obvious that

[TABLE]

which implies that $f^{\prime}(\theta_{d};2,j)>0$. Thus, $\dfrac{dq_{2}}{d\theta}(\theta_{d})=\sum_{j=1}^{2}f^{\prime}(\theta_{d};2,j)>0$.

Let us prove that $-\dfrac{dq_{1}}{d\theta}(\theta_{d})>0$ under (3.8). By (3.6) and (3.7), it follows that

[TABLE]

which yields

[TABLE]

Thus, $f^{\prime}(\theta_{d};1,j)<0$ for $j=1,2$, which implies $-\dfrac{dq_{1}}{d\theta}(\theta_{d})=-\sum_{j=1}^{2}f^{\prime}(\theta_{d};1,j)>0$.

The following theorem asserts that there exists a sufficient condition for the asymptotic stability of $\theta=\theta_{d}$.

Theorem 5**.**

Assume that (3.4), (3.5) and (3.7) hold. Then, for $a_{ij},b_{ij},c_{ij}$ given by (3.2), the equilibrium point $\theta_{d}\in\bigcup_{i,j=1,2}\Theta_{ij}\cap\Theta_{0}$ of (1.2) is asymptotically stable, where $\Theta_{ij}$ and $\Theta_{0}$ are defined by

[TABLE]

4 Proof of Theorem 5

By letting

[TABLE]

(2.2) becomes (1.2). From Proposition 2, we see that $\theta_{d}$ is an asymptotically stable point if $P$ takes strictly the minimum at $\theta_{d}$. By the definition of $v(\theta_{d})$,

[TABLE]

By Proposition 3, under (3.8), if

[TABLE]

then

[TABLE]

for $k>0$. Hence, let us discuss the existence of $(\theta,a_{ij},b_{ij},c_{ij})$ that $f^{\prime\prime}(\theta)>0$ holds.

Differentiating (3.1), we have

[TABLE]

where

[TABLE]

By (3.3) and (3.4), the discriminant of a quadratic function implies $f^{\prime\prime}(\theta)>0$ for $\theta\in(-\theta_{0},\pi)$, satisfying

[TABLE]

We observe that $\sqrt{B^{2}-4A^{2}}=|a_{ij}^{2}-b_{ij}^{2}+c_{ij}^{2}|$, which implies that

[TABLE]

in the case when $\sqrt{a_{ij}^{2}+c_{ij}^{2}}>|b_{ij}|$, and

[TABLE]

in the case when $\sqrt{a_{ij}^{2}+c_{ij}^{2}}<|b_{ij}|$.

The following lemma shows a sufficient condition for $f^{\prime\prime}(\theta;1,j)>0$.

Lemma 6**.**

Let $i=1$. Assume that

[TABLE]

where $\sin^{-1}:[-1,1]\to[-\pi/2,\pi/2]$ is the inverse sine function,

[TABLE]

Then, $f^{\prime\prime}(\theta;1,j)>0$ for $j=1,2$.

Proof 4.1**.**

We drop the indexes “$1j$” for simplicity. Since $\alpha\in(0,\pi/2)$,

[TABLE]

We first suppose that $\sqrt{a^{2}+c^{2}}>|b|$. In view of (4.1) and (4.2), it follows that $f^{\prime\prime}(\theta)>0$ for $\theta\in(-\theta_{0},\pi)$, satisfying

[TABLE]

By the definition of $\gamma_{1}$, (4.3) implies (4.4).

In the same manner, we can see that $f^{\prime\prime}>0$ in the case when $\sqrt{a^{2}+c^{2}}<|b|$, which completes the proof.

We next show a sufficient condition for $f(\theta;2,j)>0$.

Lemma 7**.**

Let $i=2$. Let (3.7) hold. Assume that

[TABLE]

where

[TABLE]

Then, $f^{\prime\prime}(\theta;2,j)>0$ for $j=1,2$.

Proof 4.2**.**

We drop the indexes “$2j$” for simplicity.

We first assume that $|b|<\min\left\{\sqrt{a^{2}+c^{2}},\ C_{\theta_{0}}\right\}$. By (3.7) and $\alpha\in(\pi/2,\pi)$, we have

[TABLE]

In view of (4.1), it follows that $f^{\prime\prime}(\theta)>0$ for $\theta\in(-\theta_{0},\pi)$, satisfying

[TABLE]

which holds by (4.5).

We similarly conclude that $f^{\prime\prime}(\theta)>0$ for

[TABLE]

in the case when $|b|>\sqrt{a^{2}+c^{2}}\max\left\{C_{\theta_{0}}^{-1}\sqrt{a^{2}+c^{2}},\ 1\right\}$.

Thanks to Lemma 6 and 7, it follows that $\theta_{d}\in\bigcup_{i,j=1,2}\Theta_{ij}\cap\Theta_{0}$ is an asymptotic stable point of (1.2). In fact, the following example is contained in the above.

Example 4.3**.**

We take $L_{0}-b_{1}=2\kappa$, $\ell_{1}=\kappa$, $d_{1}=2\kappa$, $r_{1}=\dfrac{\kappa}{2\sqrt{2}}$, $s_{1}=\dfrac{\kappa}{2\sqrt{2}}$, $L_{0}-b_{2}=2\sqrt{2}\kappa$, $\ell_{2}=\kappa$, $d_{2}=2\sqrt{2}\kappa$, $r_{2}=\dfrac{\kappa}{2\sqrt{2}}$ and $s_{2}=\dfrac{\kappa}{2\sqrt{2}}$ for any $\kappa\in\mathbb{R}$. Then, we have

[TABLE]

Thanks to Proposition 3, inner forces $v_{d1},v_{d2}>0$. Theorem 5 shows that $\theta=\theta_{d}\in(0,\pi/6)$ is asymptotically stable.

5 Numerical and experimental results

This section shows the numerical simulation results and experimental results in a stable case based on Example 4.3. We also give an unstable case to show that the motion convergence of the system in Figure 1 is not always stable, although we did not discuss an unstable condition above.

5.1 Numerical simulation

Since the solution may be unstable, we use a penalized equation for (1.2) in order to simulate.

[TABLE]

where $\theta_{\text{max}}=41\pi/180$, $\theta_{\text{min}}=-\pi/180$, $\varepsilon=1/1000$ and

[TABLE]

Let $\theta^{\varepsilon}$ be a solution of (5.1). We note that $\theta^{\varepsilon}$ coincides with a solution of (1.2) if it is always in $(\theta_{\text{min}},\theta_{\text{max}})$. We also note that $\theta^{\varepsilon}$ is an approximation of a solution of (1.2) with $\theta_{\text{min}}\leq\theta\leq\theta_{\text{max}}$ for small $\varepsilon>0$. More precisely, $\theta^{\varepsilon}$ converges to a solution of the bilateral obstacle problem,

[TABLE]

as $\varepsilon\to 0$. In fact, it holds true in at least the viscosity solution sense with an appropriate initial condition. We refer to [8] for this fact.

We first show a numerical simulation result in Figure 4 in the case of Example 4.3. We take the same parameters as Example 4.3 (unit: [mm]) for $\kappa=30$, and $L_{0}=70$ [mm], $L_{1}=15$ [mm], $I=4.2\times 10^{-3}$ [kg $\mathrm{m}^{2}$], $\mu=0.1$ [-]. The following figures are the shape of the potential energy $P(\theta)$ and a motion behavior of the joint angle $\theta$ in the case when $\theta_{d}=\pi/12$.

We next show an unstable case in Figure 5, where the parameters are too far from a condition in Theorem 5. We take the same values of $L_{0}$, $L_{1}$, $I$, $\mu$ and $\theta_{d}$ as the above simulation. We also choose $b_{i}=20$ [mm], $d_{i}=30$ [mm], $\ell_{i}=30$, $r_{i}=15$ [mm], $s_{i}=25$ [mm] for $i=1,2$.

5.2 Experimental results

We verify a sufficient condition obtained in Theorem 5 through an experiment with a real 1-link-2-muscle musculoskeletal system with routing points that is mechanically made (Figure 6).

The real system has a joint, and two fishing lines instead of muscles. The link moves in a two-dimensional plane. Thus, we can ignore gravity effects. A rotary encoder is installed to measure the joint angle in the joint. For the target system shown in Figure 1, the rotation angle of the virtual links $\ell_{1}$ and $\ell_{2}$ depends on the joint angle. To accomplish this for the real system, we use several gears and replicate the virtual links. The fishing lines pass to contact points $P_{11},P_{21}$ with the swing guide pulleys from the endpoints $P_{13},P_{23}$ of the link through the routing points $P_{12},P_{22}$, respectively. The small pulleys of routing points $P_{12},P_{22}$ make the fishing lines smooth. A weight at an endpoint of a fishing line generates the fishing line’s tension.

We observe the motion behavior of the joint angle when we provide the real system the step input $v(\theta_{d})=(v_{d1},v_{d2})$ balancing at a target angle $\theta_{d}$. By changing the placement of the swing pulleys, we can replace parameters $(d_{1},d_{2},b_{1},b_{2})$ by another. However, $(s_{1},s_{2},r_{1},r_{2},\ell_{1},\ell_{2})$ cannot be changed. We arrange the real systems in a stable case and an unstable case by replacing $(d_{1},d_{2})$.

Finally, we give motion behaviors of the joint angle through an experiment in a stable case based on Example 4.3 and an unstable case, as shown in Figure 7. We take the initial angle $\theta(0)=-5\pi/18$ [rad], the initial angular velocity $\dot{\theta}(0)=0$ [rad/s] and the target angle $\theta_{d}=\pi/12$ [rad]. Letting $\kappa\fallingdotseq 140/\sqrt{2}$ in Example 4.3, we choose the parameters (Table 1).