Stability analysis of multi-serial-link mechanism driven by antagonistic multiarticular artificial muscles

Artificial multiarticular musculoskeletal systems consisting of serially connected links driven by monoarticular and multiarticular muscles, which are often inspired by vertebrates, enable robots to elicit dynamic, elegant, and flexible movements. However, serial links driven by multiarticular muscles can cause unstable motion (e.g., buckling). The stability of musculoskeletal mechanisms driven by antagonistic multiarticular muscles depends on the muscle configuration, origin/insertion of muscles, spring constants of muscles, contracting force of muscles, and other factors. We analyze the stability of a multi-serial-link mechanism driven by antagonistic multiarticular muscles aiming to avoid buckling and other undesired motions. We theoretically derive the potential energy of the system and the stable condition at the target point, and validate the results through dynamic simulations and experiments. This paper presents the static stability criteria of serially linked robots, which are redundantly driven by monoarticular and multiarticular muscles, resulting in the design and control guidelines for those robots.


Background
Vertebrates, including humans, have multiarticulate musculoskeletal systems to voluntarily exert large forces and perform flexible movements with adaptability to the environment. A musculoskeletal system has two types of muscles: (1) multiarticular muscles, which span two or more joints and act on the joints simultaneously, and (2) monoarticular muscles, which act on a single joint. Multiarticular muscles are considered to enable large movements involving large parts of the musculoskeletal system, while monoarticular muscles are mainly used for small movements and joint support.
We have developed an artificial multiarticulate musculoskeletal robot that mimics the structure of necks of giraffes and humans [1,2]. During the development, we faced a considerable technical challenge in attaining a stable stability. The serial link mechanisms that are redundantly driven by monoarticular and multiarticular muscles easily cause unintended postures. For example, when the monoarticular muscles are not strong enough, the robot often deforms to a bow-shape (as shown in the Fig. 1 and the left in Fig. 2) or a zigzag-shape, even if the muscles on the right side and the left side of the robot contract at an equal rate (as shown in the right in Fig. 2). We refer to these instability postures as "buckling" since a bow-shape corresponds to the first order buckling and zigzag-shape corresponds to a higher order buckling in beam structure. The buckling of such serial-link mechanisms is well described in [3]. In this study, we aim to clarify the stability criteria.
stiffness. The author referred the stable and unstable configurations, as the minimum and maximum points of the potential energy, respectively. In this study, we provide a similar stability analysis by considering the musculoskeletal robot with multi-articular muscles.
Various types of robots driven by multiarticular muscles mimicking upper and lower limbs, vertebrae, and other anatomical structures have been proposed. The endpoint stiffness control of robots mimicking human upper limbs has been studied for a three-bar serial-link model driven by monoarticular and biarticular muscles [4,5]. In addition, a stable muscle arrangement has been studied using similar models [6]. However, since only two joints and biarticular muscles have been considered, no zigzag buckling can occur. Moreover, although the stability from the viewpoint of potential energy has been discussed in [6], the tensile force in muscles were assumed to be constant relative to its length. Such as assumptions fails to describe the pneumatic artificial muscles, as the tensile force and stiffness change with the length and the driving force, respectively.  For multiarticulate musculoskeletal lower-limb robots, the reaction forces from the ground and the jumping force were controlled by considering the muscle elasticity [7][8][9]. Moreover, two joints with biarticular muscles were also considered, as in the case of the upper-limb robots. The stiffness of the biarticular muscle and the muscle model of the lower limb considering it as a passive spring were studied in [7]. Subsequently, antagonistic muscles spanning different joints were studied in [8], while a model with antagonistic biarticular muscles spanning the same joints was presented in [9]. However, none of these studies have addressed instability in the musculoskeletal system.
Studies on multiarticulate musculoskeletal robots mimicking the upper and lower limbs have dealt with systems driven by multiarticular muscles spanning three or more joints. However, the instability conditions that we address have been neglected. Multiarticulate musculoskeletal robots consisting of three or more joints and multiarticular muscles spanning the joints have been inspired by the spine of animals such as humans [10,11], lancelets [12], and snakes [13,14]. The robots introduced in [10,11] consisted of silicone rubber to support each joint, but their instability conditions were neglected. In [12], lancelet swimming was studied using a nine-joint musculoskeletal system with seven pairs of tri-articular muscles, which were arranged antagonistically. However, each antagonistic pair was driven by one motor, preventing the analysis of simultaneous contraction of antagonistic muscles and system buckling. In [13], snake locomotion was simulated, and the efficiency of multiarticular muscles was studied by changing the number of joints spanned by monoarticular up to tri-articular muscles. In [14], the superiority of multiarticular muscles over monoarticular muscles was verified with respect to force per cross-sectional area and energy efficiency. However, instability was not addressed.
In [15], the stability of multiarticular tendon-driven robots was studied in general considering the elastic energy of the tendons, and the minimum length up to which the tendons are required to be stretched was theoretically derived. However, in the corresponding wire rope-pulley system, the moment arm around each joint exerted by tension of each wire was assumed to be constant. Consequently, the analysis of this system cannot be applied to musculoskeletal systems, in which the moment arm changes nonlinearly with respect to the joint angles.

Purpose
Although various studies discussing the multiarticulate musculoskeletal systems are available, the literature dealing with multiarticular muscles spanning three or more joints is very limited. Moreover, most of the studies neglect buckling that may occur when activating multiarticular muscles. Additionally, the posture stability analysis by considering simultaneous contraction of antagonistic multiarticular muscles and dependence of muscle stiffness on the driving force has not been studied. Therefore, studies focused on the determination of appropriate force and stiffness of monoarticular and multiarticular muscles for stable operation of a musculoskeletal system provide good research opportunities.
Since thin and flexible artificial muscles were not previously available for practical use, this type of analysis has not been a significant issue in robotics. However, with the recent availability of such thin artificial muscles, the postural stability of musculoskeletal robots becomes significantly challenging.
Therefore, in this study, considering that the stiffness and the maximum contraction ratio of artificial muscles physically change with the driving force, we aimed to clarify the required parameters of monoarticular and multiarticular muscles to avoid buckling during multiarticular muscle activation. The two important contributions of this paper are as follows: (1) a clarity on the important parameters that are required for the static stability of the multiarticular musculoskeletal serial-link robot without buckling will be provided; (2) the analysis would be validated with experiments using thin McKibben artificial muscles. The static stability analysis of a multiarticulate musculoskeletal system using the potential energy approach is carried out for a multi-serial-link mechanism. The analysis consists of modeling the artificial muscles as spring elements. In addition, we considered static stability as the system convergence to a target angle and buckling instability if the system buckles in zigzag or bends uniformly like a bow. In addition to the model analysis, we verified the stability formulation through dynamic simulations implemented in Math-Works MATLAB and experiments using the McKibben artificial muscles.

Modeling and stability analysis
Our model of the musculoskeletal system is a planar multi-serial-link mechanism that is antagonistically driven by the artificial muscles, as shown in Fig. 3a. Figure 3b illustrates the fixed position of the multiarticular muscles spanning from the joint p to the joint q. Figure 3 shows the m-th pair of antagonistic muscles attached to the vertical fixtures at a distance a m from the axis along each link and at distances e 0,m and e 1,m from each joint along each link. Thus, e 0,m and e 1,m represent the distance of a joint with a number and that of a joint with the subsequent number, respectively. In addition, the stiffness of the left and right sides of the m-th pair of antagonistic muscles, the link lengths, and the joint angles are denoted as K L,m , K R,m , l 0 . . . l N , and θ 1 . . . θ N , respectively. It is noted that both sides of each antagonistic muscle are driven with the same contraction force, as detailed below.

McKibben artificial muscle
In this study, we used McKibben pneumatic artificial muscles. For the experiments, we used McKibben muscles with a diameter of 2 mm. This type of artificial muscle mainly consists of an inner rubber tube to which air pressure is applied, and an outer sleeve of a woven fiber intersects at a certain angle. The sleeve contracts through expansion in the radial direction of the inner tube.
Among the various formulations, the simplest equation to describe the McKibben muscle operation has been proposed by Schulte [16]: where F , D, P, θ 0 , and ε are the contraction force, diameter of the inner tube, air pressure applied to the tube, woven angle of the sleeve, and contraction ratio of the artificial muscle, respectively. Equation (1) is based on the virtual work principle. The work done by the compressed air in the rubber tube and that done by the artificial muscle are assumed to be equal. Thus, the corresponding geometric calculation neglects the elasticity of the rubber tube and sleeve strands. The maximum contraction force, which occurs when the contraction rate is 0, appears to agree with experimental data, but the maximum contraction ratio may be inconsistent.
In the next section, the McKibben artificial muscle is modelled as a spring with stiffness K and natural length s 0 , as expressed in Eq. (2). The maximum contraction force, F Max , is derived from Eq. (1), and the maximum contraction ratio, ε Max , is obtained experimentally. Note that S represents the default length of the muscle, and muscle stiffness K varies with the driving pressure P.

Potential energy of multi-serial-link mechanism driven by antagonistic artificial muscles
We also consider the stability of the multiarticulate musculoskeletal system driven by antagonistic muscles. We focus on a straight posture, which is common for the spine. If the total potential energy, U , which is produced by the muscles considered as springs, reaches a local minimum value at the target angles, the multiarticular musculoskeletal system is said to be stable. A similar stability analysis by considering the potential energy is conducted by Ochi et al. [6], and Ozawa et al. [15]. Note that, unlike the analysis of the latter researchers, in this paper, it is not necessary to consider the gravitational potential energy and kinematic energy terms because we focused on the static stability of the planar serial link. Therefore, considering an N-joint serial-link mechanism, a sufficient condition for stability is given by the entire potential energy, U (θ) ( θ ∈ R N ), which satisfies the following conditions: where H ∈ R N ×N is the Hessian matrix describing the second-order partial derivative of U . Specifically, the element at row i and column j of matrix H is given by Equation (4) is equivalent to the following equation: The potential energy of the entire system, U , can be expressed as the sum of the potential energy produced by the antagonistic muscles: where M is the number of pairs of the antagonistic muscles, s L,m and s R,m are the endpoint vectors of the left and right sides of antagonistic muscle m, respectively, and |s| is the norm of vector s . We set the natural length of both muscle sides to s 0,m . As shown in Fig. 3b, the endpoints of the multiarticular muscles that span from joint p to joint q have distances e 0,m and e 1,m from each joint. As the muscle endpoints are attached to the next links, e 0,m and e 1,m must satisfy the following constraints: By using the angles of the mechanism, s L,m and s R,m can be described as follows: where φ n = n k=p θ k . The only difference between s L,m and s R,m is the change in the positive and negative signs of a m given the symmetry of the antagonistic muscles.
For an antagonistic muscle m, the first-and secondorder partial derivatives of its potential energy U m are given by: where U m is assumed to be a continuous and secondorder differentiable function in R N , and the order of differentiation can be ignored.
As the target angles are set to θ = − → 0 ∈ R N , we have 0 ≤ e 0,m ≤ l p−1 , 0 ≤ e 1,m ≤ l q s L,m = e 0,m + q−1 n=p l n cosφ n + e 1,m cosφ q + a m sinφ q a m + q−1 n l n sinφ n + e 1,m sinφ q − a m cosφ q s R,m = e 0,m + q−1 n=p l n cosφ n + e 1,m cosφ q − a m sinφ q −a m + q−1 n l n sinφ n + e 1,m cosφ q + a m cosφ q where ε Max is the maximum contraction ratio of the muscle, and the muscle length, S m , can be expressed as given below: We consider contraction ratio ε Max to be always positive because the McKibben artificial muscles are assumed to be of the contraction type. Therefore, the second term related to ε Max in Eq. (14) is always negative. Thus, a larger a m leads to more positive elements, H ij , in the Hessian matrix, whereas large values of e 0,m , e 1,m , and link lengths lead to more negative elements. From Eq. (14), we can check the stability of the system by determining whether H ij satisfies the Eq. (4) or Eq.

Stability of monoarticular and multiarticular muscles
For a negative diagonal element of the Hessian matrix of the potential energy, Eq. (4) is not satisfied. This is because any negative H ii results in the following left-hand side of Eq. (4) obtained by setting x as The abovementioned condition is related to the stability of joint i but not to that of the entire system. Therefore, for system stability, all the diagonal elements of the Hessian matrix should be positive.
When joints q and p are the same, that is, when considering the antagonistic monoarticular muscles that span joint p, the corresponding Hessian matrix of potential energy is given by Equation (15) shows that the monoarticular muscles affect only the p-th diagonal element of the Hessian matrix, H pp . Since the diagonal elements must be positive for stability, the antagonistic pairs of the monoarticular muscles stabilize the joint p when H m,pp is positive. This can be expressed as follows: Equation (16) indicates that e 0,m and/or e 1,m must be sufficiently smaller than a m for the monoarticular muscles to stabilize the joint.
When the joints p and q are different, the diagonal elements of the Hessian matrix of potential energy produced by antagonistic multiarticular muscles can expressed as In the Eq. (17), the p-th and q-th diagonal elements become positive when e 0,m and e 1,m are sufficiently close to 0. In contrast to this, the diagonal elements for the joints between p and q are usually negative because a m is not much larger than the sum of the link lengths. Thus, multiarticulate antagonistic muscles cause instability, especially in the joints between them.

Simulations
We first validated our formulation through calculations and simulations implemented in MathWorks MATLAB Simulink.  For the simulations and the experiments reported below, we considered a three-joint multi-serial-link mechanism with four pairs of antagonistic muscles comprising three pairs of monoarticular muscles and one pair of tri-articular muscles. The parameters of the McKibben artificial muscles used for the simulations and the experiments are as listed in Table 1. A mechanism considering the parameters listed In Table 1 are illustrated in the Fig. 4.
It is to be noted that, the heights of both the muscle endpoints are different in the experimental setup. Fig. 5 Results of dynamic simulations. One spring representing a tri-articular muscle (TM) is attached to both ends of the system, and three springs representing monoarticular muscles (MMs) are attached around each joint. Note that all the MMs are driven by the same pressure at each condition Table 3 Determinants of A 1 , A 2 , and A 3 indicating stability in different cases  Therefore, a height gap of h = 9 mm is considered in the simulation. Moreover, the vertical distance of all muscles from links a 1 , . . . , a 4 was set to 19.5 mm. Assuming that the antagonistic artificial muscles are tilted by angle ψ in the height direction from the horizontal line, the default muscle length, S m , and the vectors of the antagonistic muscles vary from those mentioned above. As a result, the Hessian matrix of the potential energy produced by the muscles that is derived using Eq. (13) can be rewritten as follows: with angle ψ expressed as follows: Table 2 lists the experimental parameters associated to the monoarticular and the tri-articular muscles of the McKibben artificial muscles. The stiffness of the artificial muscles, K m , was derived from the Eq. F Max that is obtained using the Eq. (1). In addition, ε Max was determined from measurements on a tri-articular muscle at varying applied pressure P . In this case, the ε Max is measured as the contraction ratio of the muscle when no tension was applied to it. Specifically, the monoarticular and the tri-articular muscles were driven by air pressures of 0.0, 0.3, 0.4, and 0.5 MPa. We measured 15 trials per air pressure except for the pressure of 0.0 MPa, which is equivalent to the muscles being not attached to the joint because no interference occurs to the muscles and links in the simulation. The simulation results are shown in the Fig. 5.
Using the values in Tables 1 and 2, the determinants of A 1 , A 2 , and A 3 as defined in Eq. (5) are obtained and are listed in the Table 3. When no pressure (0.0 MPa) is applied to the monoarticular muscles and a pressure of 0.4 or 0.5 MPa is applied to the tri-articular muscles, at least one of the determinants from A 1 , A 2 , and A 3 becomes negative. Thus, instability occurs as confirmed by the simulation results shown in Fig. 5. When the determinant of A 1 is negative, instability appears as zigzag buckling.

Experiments
We also verified our formulation experimentally using a three-joint multi-serial-link mechanism and the McKibben artificial muscles. A diagram and a photograph of the experimental setup are shown in Fig. 6. An air compressor was used to provide the compressed air. The air pressure was controlled using a pressure valve and applied to the monoarticular and the tri-articular muscles. Bearings were used at all joints to reduce friction. Whereas the friction between the link and the ground was reduced by using light weight links which are fabricated using a 3D printer. The parameter of the McKibben muscles used in the experiment are as listed in the Table 2.
Experiments are performed for the cases with the same conditions that are considered during the simulation. The obtained experimental results are as shown in Fig. 7. When the pressure applied to the monoarticular muscles is 0.0 MPa, each joint angle changes substantially, clearly destabilizing the system at the target posture. When monoarticular muscles are driven at any pressure except for 0.5 MPa and tri-articular muscles are driven at 0.4 or 0.5 MPa, the joint angles change in a smaller proportion than without using the monoarticular muscles. When the monoarticular muscles are driven at 0.5 MPa, it is difficult to recognize the changes. When the monoarticular muscles are driven at 0.3 MPa, no remarkable changes occur regardless of the pressure applied to the tri-articular muscles. Overall, the experimental results are consistent with those obtained from our formulation.

Discussion
The results presented in the Table 3 show that the system becomes theoretically unstable at the target posture θ = 0 for the monoarticular muscles driven at 0.0 MPa or the tri-articular muscles driven at 0.4 or 0.5 MPa. In particular, when the monoarticular muscles are inactive (at 0.0 MPa pressure), the determinant of A 1 becomes negative, and the serial-link mechanism buckles to a zigzag shape. This phenomenon was successfully predicted by the simulations and was obtained using the experiments also (refer to the Figs. 5 and 7). For other theoretically unstable conditions, the simulation predicted the system to bend like a bow as shown in the Fig. 5. On the other hand, the experimental results shown in the Fig. 7 show apparent instability of the serial-link mechanism because changes in each joint angle can be observed, but the mechanism does not bend with the extent shown in the simulations. In the theoretically stable conditions, wherein the tri-articular muscles are driven at 0.0 or 0.3 MPa, the mechanism does not bend and remains stable. This observation is verified using the simulation and the experimental results.
When the tri-articular muscles are driven at 0.5 MPa, the mechanism tends to stabilize as the monoarticular muscles are driven at higher pressures in the experiments, contradicting the simulation results. This inconsistent result is attributable to the effect of joint friction, which increases with the muscle contraction forces.
In addition to the effect of joint friction, the friction caused by the interference between the artificial muscles and the links tends to stabilize the posture of the mechanism. Therefore, the serial-link mechanism can remain stable at a posture around the target, at least when the theoretically stable conditions are satisfied.

Conclusion
Considering monoarticular and multiarticular muscles in a multiarticulate musculoskeletal system, we analyze the conditions for system stability while driving the multiarticular muscles. Theoretically, we determine stability assuming that it can be reached when the potential energy produced by each articular muscle is locally minimal at the target point. In addition, to stabilize each joint, the distances from the fixed point of the monoarticular muscles to the joint must be sufficiently small. We analyzed the stability of a three-joint multiarticulate musculoskeletal system and validated our formulation through dynamic simulations and experiments on a mechanism driven by McKibben artificial muscles. We confirmed that sufficiently high stiffness of monoarticular muscles and locally minimal potential energy of the system at the target posture lead to stability in the musculoskeletal mechanism consisting of serial links and the multiarticulate McKibben muscles.