On Control of Reaching Movements for Musculo-skeletal Redundant Arm Model

This paper focuses on a dynamic sensory-motor control mechanism of reaching movements for a musculo-skeletal redundant arm model. The formulation of a musculo-skeletal redundant arm system, which takes into account non-linear muscle properties obtained by some physiological understandings, is introduced and numerical simulations are perfomed. The non-linear properties of muscle dynamics make it possible to modulate the viscosity of the joints, and the end point of the arm converges to the desired point with a simple task-space feedback when adequate internal forces are chosen, regardless of the redundancy of the joint. Numerical simulations were performed and the effectiveness of our control scheme is discussed through these results. The results suggest that the reaching movements can be achieved using only a simple task-space feedback scheme together with the internal force effect that comes from non-linear properties of skeletal muscles without any complex mathematical computation such as an inverse dynamics or optimal trajectory derivation. In addition, the dynamic damping ellipsoid for evaluating how the internal forces can be determined is introduced. The task-space feedback is extended to the 'virtual spring-damper hypothesis' based on the research by Arimoto et al. (2006) to reduce the muscle output forces and heterogeneity of convergence depending on the initial state and desired position. The research suggests a new direction for studies of brain-motor control mechanism of human movements. 1. Introduction The natural movements of a human seem to be smooth, dexterous and sophisticated compared with today's robots. Until now, these excellent movements have attracted many physiologists, kinesiologists and robotics researchers who hope to understand and extract them. In paticular, obtaining and mimicking these natural movements in robotics has become a hot topic, especially for robots that will be used in close proximity to humans. Obtaining the fundamental elements of how the central nervous system (CNS) generates control strategies to perform these kinds of movements is one of the ultimate objectives and dreams for understanding the nature of human movements. In robotics, many human movements such as object manipulations, walking, etc., have been mimicked. However, achieving highly natural movements, like human behaviour, is still difficult. There are many reasons for these difficulties. One is the existence of ill-posed problems induced by the joint and muscle redundancies that a human body intrinsically possesses (Bernstein 1935). Until now, many solutions for overcoming these ill-posed problems by introducing various optimisation criteria have been proposed not only in physiology but also in …


Introduction
The natural movements of a human seem to be smooth, dexterous and sophisticated compared with today's robots. Until now, these excellent movements have attracted many physiologists, kinesiologists and robotics researchers who hope to understand and extract them. In paticular, obtaining and mimicking these natural movements in robotics has become a hot topic, especially for robots that will be used in close proximity to humans. Obtaining the fundamental elements of how the central nervous system (CNS) generates control strategies to perform these kinds of movements is one of the ultimate objectives and dreams for understanding the nature of human movements. In robotics, many human movements such as object manipulations, walking, etc., have been mimicked. However, achieving highly natural movements, like human behaviour, is still difficult. There are many reasons for these difficulties. One is the existence of ill-posed problems induced by the joint and muscle redundancies that a human body intrinsically possesses (Bernstein 1935). Until now, many solutions for overcoming these ill-posed problems by introducing various optimisation criteria have been proposed not only in physiology but also in robotics (Hogan 1984a Katayama and Kawato 1993;Kawato et al. 1987Kawato et al. , 1990Nakamura 1991). In general, these optimisations require large calculation costs to solve complicated optimal functions as well as an inverse dynamics, which would be very difficult for a human to perform in real time. Indeed, do humans perform actually these complex calculations every time they move?
In this paper, we suggest that the sensory-motor control mechanism, which enables many kinds of dexterous human movements by the CNS, may not be so complicated and may even use simple task-space feedback by utilising intrinsic non-linear mechanical characteristics of skeletal muscles. Morasso (1981) observed human reaching movements and suggested that the central commands in human movements may be coordinated not in joint space but in task space. We focus on a simple reaching task by using a three-link planar arm model driven by nine redundant muscles so as to mimic the construction of a human arm in a simple way. We introduce non-linear muscle characteristics that correspond to the findings of physiological studies (Hill 1938;Mashima et al. 1972) and propose a simple task-space feedback control scheme for the arm model taking into account internal forces induced by redundant muscles. These internal forces can be basically generated by the co-contraction of agonist and antagonist muscles. Arimoto et al. (2015) suggested that the performence of the convergence to the desired position for human-like reaching movements by using a joint redundant robotic arm depends upon the damping shaping of each joint. Unfortunately, humans cannot directly modulate the rotational joint viscosity because each joint is driven by linear movement of antagonistic muscles. How do humans modulate the rotational joint viscosity? Hogan (1984b) remarked and Gribble et al. (2013) observed that the co-contraction of agonist and antagonist muscles contributes to the modulation of the joint impedance. In our control scheme presented here, the co-contraction of antagonistic muscles to invoke the internal forces coming from the redundant muscles can modulate the damping shaping of each joint if we introduce non-linear muscle dynamics based on Mashima et al.'s model (1972).
Basically, a redundant actuation system such as a musculo-skeletal arm can generate internal forces. These internal forces belong to the null-space of the muscle space with respect to the joint space and thereby they are independent of a principal part of the joint torques in relation to generation of movements of the arm end point towards the target. However, it causes an ill-posed problem of how to determine the internal forces uniquely. Until now, many optimisational methods to determine the internal forces uniquely have been proposed by introducing some sort of artificial indices (Delp and Loan 2010;Rasmussen et al. 2013). Unlike these methods, our proposed method does not need any optimisation criteria, because the internal forces can affect only part of the joint torques irrelevant to motion of the arm end point towards the target in the task space. Nevertheless, human-like slightly curved trajectories of the end point can be realised by adequately adjusting the internal forces through the null space with respect to the joint space. This point quite differs from the traditional optimisation methods.
In our previous work (Tahara et al. 2015), we treated the two-joint six-muscle model which has only muscle redundancy, that is, joint redundancy was not taken into consideration. In this paper, we take into account not only the muscle redundancy but also the joint redundancy by introducing a three-joint nine-muscle planar arm model. The crucial difference between these two cases is the convergence analysis that uses Lyapunov's direct method to prove that the closed-loop dynamics of the joint redundant model cannot be applied. To prove the stability and convergence of the closed-loop dynamics under the existences of both muscle and joint redundancies, we introduce the novel concept of 'stability on a manifold' originated by .
In what follows, we formulate the kinematics and dynamics of the three-joint planar arm model driven by six mono-articular muscles and three bi-articular muscles governed by non-linear behaviour. The stability and conver-gence of the overall system is given by introducing a stability theory on a manifold. We report some numerical simulations to demonstrate that the simple sensory-motor control scheme makes it possible to make reaching movements even though both joint and muscle redundancies are taken into consideration. To discuss and evaluate how to choose the internal forces, we introduce a dynamic damping ellipsoid. In addition, the task-space feedback manner is extended to the 'virtual spring-damper hypothesis' proposed by Arimoto et al. (2015) to reduce the muscle output forces and the heterogeneity of convergence, which depends on the initial state and desired position of the arm. Finally, these results suggest a direction for studies of the brainmotor control mechanism of human movements.

Musculo-skeletal redundant system
We introduce a musculo-skeletal redundant arm model that consists of three serial links, six mono-articular muscles, and three bi-articular muscles to imitate the configuration of a human arm. This model is shown in Figure 1. The three bi-articular muscles enhance the end-point stiffness.

Kinematics of muscle-joint space
Assume that the overall movements of the system are limited to a horizontal plane so that gravity does not affect Figure 1. Three-link musculo-skeletal arm model with six mono-articular muscles and three bi-articular muscles. movements at all. Any friction such as Coulomb's friction and static friction is ignored in our model. All skeletal muscles included in the arm can only be linearly contracted and are not curved on the way. Also the mass transfer of the muscles during contraction is omitted because its effect on arm movements is not so crucial. The physical parameters and variables are shown in Figure 1. The lengths In this study, we assume that the Jacobian matrix W is of row full-rank during movement. Therefore, we can take the inverse of Equation (3) as where W = W (W W ) ∈ R signifies the pseudoinverse matrix of W , I − W W k ∈ R stands for the 9 e null-space of W , k ∈ R is an arbitrary vector, and I ∈ R denotes an identity matrix. The physical meaning of e 9 9 9×9 the second term of the right-hand side is the internal force generated by the redundant muscles. This relation does not mean an optimisation, and it expresses only decomposition of the muscle force space with respect to the joint space into two spaces, one is the image space and the other is the kernel space.

Kinematics of joint-task space
The end-point position of the arm model given by x = T (x, y) in Cartesian coordinates is expressed as Taking the time derivative of Equation (5) yieldṡ where J ∈ R stands for the Jacobian matrix from joint τ ∈ R and output forces of the end point in Cartesian coordinates F ∈ R can be represented as 2×3 space to task space. The relation between joint torques Substituting Equation (7) into Equation (4) yields the relation between contractile forces of the muscles F and output forces of the end-point F . It is given by

Non-linear muscle model
It is known that the skeletal muscles are governed by strong non-linearity. Here, we introduce a tangible nonlinear muscle model that is based on several physiological results. In the modeling of the muscles, we assume that the masses of all muscles are included in the masses of each link. The model of muscle is basically according to the two-element model which is composed of contractile element (CE) and series elastic element (SE). In this paper, we want to focus on the dumping effect on each joint generated by CE, and thereby we introduce quite strong assumption that the SE part can be approximated as a rigid element which indicate that its elastic coefficient becomes infinity. Therefore, hereafter, we only consider the CE part in the muscle model. In physiology, the force-velocity relation of the skeletal muscle proposed by Hill (1938), which is expressed as a simple hyperbolic equation, is known as one of the most fundamental muscle features. This equation is˙ṫ he contractile velocity of the muscle, f is the maximum isometric tensile force of the muscle which depends on its length, a is the heat constant and b is the rate constant of energy liberation. The shortening directions oḟ in physiology that a skeletal muscle exerts more tensile force in the lengthening phase than in the shortening phase. Mashima et al. (1972) determined the above feature specifically through some experiments. By considering this feature, we can modify Equation (8) as follows: where f is the output tensile force of the muscle, l is 0 f and l are defined as positive. In addition, it is knowṅ here 0 ≤ α ≤ 1 is the muscle activation level, the coefficients a and b are determined experimentally in Mashima et al. (1972) as a = |0.25 × f |, b = |0.9 × l intrinsic rest length of the muscles. Now, we define the control input to the muscles as¯× α, in order to set the dimension of the input as that of the muscle force. Therefore, Equation (10) can be represented as follows: Since a skeletal muscle intrinsically owns some damping factor, which comes from its component and material, we introduce a muscle intrinsic viscosity c > 0 independently fromt his study can be expressed as The force-velocity curve of the muscle property expressed by Equation (12)  as the contractile velocity of the muscle increases and vice versa. Now, it should be remarked that in Equation (12) thė the overall muscle dynamics in this study can be expressed as follows:⎡ It is remarkable that in Equation (13), the diagonal matrix A( ¯n on-linearity of the muscle model. Also, the skeletal muscle generates only a contractile force, so we define the con-ī s semi-positive definite. Both damping matrices C and C parameter p depending on l satisfies the inequality Equation 0 < p ≤ 1 as long as l is upper bounded. Eventually,¯} is constructed by the elements of the control input vector α (Tahara et al. 2005(Tahara et al. , 2006. This comes from the trol signal α, which is a saturated function so as to keep a semi-positive value, as shown in Figure 3. Therefore, A( α) are also positive definite. In this paper, the positive coefficient c , which is the intrinsic viscosity of muscle, is defined as c = 0.2 by considering physiological aspects.

Control input
In physiological studies of the brain-motor control mechanism in human movements, several control strategies for skeletal muscles have been proposed since Fel'dman first proposed a spring-like model called the λ-model (Fel'dman 1966(Fel'dman , 1986) which exerts muscle force by modulating the threshold of tonic stretch reflex about muscle length. Mclntyre and Bizzi (1993) also proposed the α-model which modulates spring-like characteristics of the muscles by α neuron. These models are based on the equilibrium point (EP) hypothesis proposed by Fel'dman and observed by Bizzi et al. Many studies concerned with multijoint reaching movement, which were based on the EP hypothesis, have been done (Bizzi et al. 1976(Bizzi et al. , 1992Flash 1987;Flash and Hogan 1985). Recently, Arimoto et al. (2015) have proposed a virtual spring hypothesis that a given potential energy in task space can generate a spring-like force on the end point and distribute it to each joint torque T task-space control scheme originated from Takegaki and Arimoto (1981). In developmental psychology, it is said that a potential function that emanates neuro-motor signals emerges from intentional incentive of a voluntary movement (Killeen 1992). The virtual spring hypothesis has emphasised that such a potential-generating neuro-motor signals should be coordinated not in the muscle or joint spaces but in the task space in terms of the distance between a present position of the end point and a desired one such as T 2 T plies the joint configuration of a human body. In the present study, we introduce a simple task-space feedback control scheme with internal forces, which basically follows the virtual spring hypothesis. Now, let us consider the diagonal matrix P expressed by Equation (13). Since the boundaries of its components p are given by 0 < p ≤ 1, P is positive definite. Hence, + long as the muscle Jacobian matrix W does not degeneratē expressed by through the transpose of the Jacobian matrix J . This 1 joint through the transposed Jacobian matrix J which imduring movement. The control input to the muscles α is x K x. This potential is eventually distributed to each where W = W P, and x = x − x is the end point position error between the present position x ∈ R and the desired position x ∈ R in the task space. The positive definite diagonal matrix K = k I ∈ R denotes the po-coefficient of the virtual spring. It should be noted that in Equation (14), the first term of the right-hand side denotes an artificial potential energy in task-space and produces spring-like forces to each muscle through two Ja-+ for the internal forces generated by the redundant muscles, 9 expressed by Equation (14) can be interpreted more specifically as that K x comes from intentional incentive of a T implies the joint configuration of the arm relative to the arm + + 9 e uration of the related muscles relative to the joints.
cobian matrices J and W and the second term stands T where k ∈ R is an arbitrary vector. This control signals movement which is generated in the motor cortex, and J¯ē end-point, W and (I − W W )k come from the config-

Dynamics of three-link planar arm model
The dynamics of a three-link planar arm model can be described by Lagrange's equation of motion as is the inertia matrix, which is sym-τ ∈ R are the input torques to the joints. Substituting Equations (2), (3), (13) and (14) into Equation (15) (16) shows the overall closed-loop dynamics, which is expressed in the joint coordination system.

Stability of the closed-loop dynamics
This section illustrates the convergence of the overall system by introducing the concepts of 'stability on a manifold' and 'transferability to a submanifold'. Now, it is convenient to rewrite Equation (16) as where u = 0 are the input to the closed-loop dynamics of Equation (16). Taking the inner product between the input as Equation (17) and the output as the joint angular velocity Here, we should consider the damping matrix W P( AC + T components of the diagonal matrix A, are defined as a saturated function (see Figure 3) to be semi-positive. Therefore, the matrix ( AC + C ) is positive definite because A ≥ 0, C > 0, C > 0 and P > 0. Hence, the damping T the Jacobian matrix W is of row full-rank. Now, integrating Equation (18)  This inequality means that the closed-loop dynamics of Equation (16) satisfies the passivity condition. However, it should be noted that this system is a joint redundant system because the control input of task-space feedback requires only an end-point position (x, y) although the arm consists of three joints. Then, the scalar function E(t) is non-negative definite, but not positive definite for all statė of this redundant system, we introduce the novel concept of 'stability theory on a manifold' proposed by Arimoto et al. (2015). In this proof, we consider that the one-dimensional variables of (θ, θ ). Therefore, to discuss the convergence submanifold M on the spectrum S ∈ R such that the 3×2 1 1 0 hood of δ stays on N (ε, r ) in the neighbourhood of ε, the reference state (θ , 0) is said to be stable on a manifold.
θ + βk J x yields Definition 2 (Transferability to a submanifold) 0 1 stant ε > 0 and another constant r > 0 (r < r ) such that any solution trajectory starting from an arbitrary initial 6 6 verges asymptotically as t → ∞ to a point on a submani-0 1 1 0 be transferable to a submanifold of M .
2 energy of the system does not increase during movement. In addition, the meaning of Definition 2 is that the desireḋ α Taking the inner product between Equation (16) and transposed Jacobian matrix J ∈ R is of column full-For the reference state (θ = θ , θ = 0) ∈ M , if for an ar-(θ(0), θ (0)) contained in N (δ(ε), r ) in the neighbour- If for the reference state (θ , 0) ∈ M , there exist a con-1 1 1 0 state contained in N (δ(ε), r ) stays on N (ε , r ) and con- The physical meaning of Definition 1 is that the overall state ( x = 0, θ = 0) is guaranteed while the system is stable.
where β is a positive number. We introduce the following equation. dḢ Substituting Equation (22) into Equation (21) yields The second, third and sixth terms of the left-hand side of Equation (21) can be rewritten as follows: T 0 where γ is a positive number. Also, the following inequality can be obtained by considering the positive definiteness Similarly, we can obtain the following inequality: T 0

γβk γ αk x J Hθ ≤ γβk x J H J x + 4
By substituting Equations (26) and (27) into Equation (25), we obtain the following inequality: Here, we introduce the scalar function V (t) such that

Substituting Equations (28) and (29) into Equation (23) yields
f their coefficients is less than or equal to the maximum eigenvalue of H. Also, J consists of a trigonometric function and physical parameters, and they are thereby totally bounded. Thus, as long as x(t) < x(0) , there exists a positive number ζ for h to satisfy the following inequality (Tahara et al. 2015(Tahara et al. , 2016: 2 Now, assume that we can choose the numbers k, β, γ and the internal force vector k to satisfy the following inequalities. In Equation (25), since H and S are in linear homogeneous with respect to the angular velocities θ , and the order 0˙ė By considering these inequalities, we see that the scalar function V (t) expressed by Equation (29) is positive, and the inequality is satisfied. Hence, Equation (23) can be reduced to Therefore, from Equation (37) we obtain and conclude that the scalar function V (t) converges to zero exponentially. The scalar function V (t) plays the role of a modified energy function of E(t) in Equation (19).

Q.E.D.
Therefore, θ → 0 and x → 0 when t → ∞.  Table 1 shows the physical parameters of the arm model, Table 2 shows the insertions of the muscles and Table 3 shows the initial and desired positions of the end point.

Numerical simulation
In this paper, we treat two types of reaching movements: short-range and middle-range reaching. We also determine the internal force vector as k = k e in order to see the influence of k on the dynamics, where k is a positive  and that of the y-position is about 2.0 s when k = 20.0. On the other hand, Figure 10 shows the simulation result of the middle-range reaching movement with k = 10.0 and k = 5.0. The trajectory of the end point is not smooth and describe a more winding route than that of the short-range reaching, and the final posture of the arm seems to be impossible for a real human. In contrast, Figure 11 shows the simulation result of the middle-range reaching movement with k = 10.0 and k = 25.0. The end-point converges to the desired point smoothly, and its trajectory becomes a quasi-straight line even though x ≤ 0.6 m holds for in middle-range reaching. Figure 12 shows trajectories of the end-point position (x, y) of the middle-range reaching movement with k = 10.0 when k is shifted from k = 5.0 to k = 30.0 by 5.0 units. Like the short-range reaching movements shown in Figure 6, the end-point trajectories become smooth and close to straight lines as k increases. Figure  These simulation results suggest that the internal force coefficieint k are chosen according to the virtual spring coefficient k, the end point of the musculo-skeletal redundant arm converges to the desired point smoothly without any complex mathematical optimization. The benefit of our proposed model is that co-contraction of antagonistic muscles can modulate the damping effect in the joint space. Until now, many control schemes for muscle redundant arm model are used by some optimisation method to determine muscle output force uniquely because the internal forces cannot affect to the joint torques. Unlike these methods, in our model, internal forces can affect the joint torques, and thereby the trajectory of the end point is changed by varying the internal force vector k .

Dynamic damping ellipsoid
We introduce a dynamic damping ellipsoid, in order to discuss how the internal force vector k is determined, and to evaluate it. The damping ellipsoid at the end point of the human arm has been proposed by Tsuji et al. (1994). The ellipsoid can express the damping effect at the end point of a human arm when the arm's dynamics can be approximated as being almost governed by only a term related to the velocity. It was derived from the result of measurement and the evaluation of a real human hand's damping effect in a quasi-static situation. Unlike the above ellipsoid, our dynamic damping ellipsoid expresses the damping effect that can be given to the end point of the arm by modulating the internal force vector k included in the control input to the muscles. Figure 16 illustrates the dynamic damping ellipsoid at the end point of a human arm. We consider part of the control input in Equation (14), which can generate the joint torque for only the damping effect as where τ is the vector of the input torques for generating the damping effect. Since the damping force F at the end point in task space can be expressed as  This ellipsoid shows the end-point damping forces generated by the control input during movement. By considering this ellipsoid, we can evaluate the internal forces required to modulate the damping effect. Figure 17 shows the damping ellipsoid at the end point of the middle-range reaching movement with k = 25.0, and Figure 18 shows that with k = 5.0. We see from these figures that in the case of a suitable internal force (k = 25.0), the long axis of  Table 4. Parameters for a virtual spring-damper.
Virtual spring k = 10.0 Virtual damper ζ = 1.0 (k = 5.0) ζ = 0.0 (k = 25.0) e e the damping ellipsoid turns on the shoulder during movements. Therefore, the short axis of the damping ellipsoid turns on the desired point, and the damping effect at the end-point is relatively small in the direction of the desired point during movement. Also at the beginning of the movement, the ellipsoid is relatively large. It becomes gradually small as the end-point velocity increases, and finally the shape of the ellipsoid becomes large again like at the beginning shape. However, in the case of an unsuitable internal force (k = 5.0), the damping ellipsoid is relatively small: the long and short axes of the ellipsoid are nearly equal. By comparing these results we can conclude that the end point can converge to the desired point smoothly by considering the size and direction of the dynamic damping ellipsoid during movement. The size and direction of dynamic damping ellipsoid are governed by the internal force vector k and muscle contractile velocity. However in this paper, we only pointed out a relation between the internal forces and the ellipsoid from the viewpoint of muscle and joint configurations as a posteriori reasoning, and could not disclose any dynamic relation between the internal force vector and the ellipsoid explicitly from the dynamics viewpoint, because an explicit dynamics of co-activation model is not presented here. Therefore, e e a more systematic paradigm for determining k according to the dynamic damping ellipsoid or a totally different modeling of muscle dynamics that may naturally prevent involvement of internal forces is needed in our next work. e

Extension to 'virtual spring-damper hypothesis'
We extend our control stragegy to the 'virtual springdamper hypothesis' proposed by Arimoto et al. (2016). This hypothesis made the point that in reaching movements, if we assume the existence of a virtual spring and a virtual damper between the end point and the desired point, the end point of a human arm can converge to the desired point smoothly according to this virtual spring-damper (see Figure 19). Also in physiology, it is known that the cocontraction of agonist and antagonist muscles is gradually  e reduced as the process of learning a movement progresses (Osu et al. 1997). Now, we hypothesise that the internal forces generated by the co-contraction of agonist and antagonist muscles are reduced, which enables an adequate damping shape in the task space to be gradually obtained as the learning process for reaching movements progresses. We modify the control input to the virtual spring-damper model. It is given bȳ Also, the closed-loop dynamics can be expressed as where the first term of the right-hand side of Equation (43) represents the virtual spring-damper and the second term represents the internal forces generated by redundant muscles. Table 4 shows each parameter for a virtual springdamper k and ζ used in the numerical simulations, and these results are shown below. Figure 20 shows the simulation results with a virtual spring-damper, and Figure 21 shows those results without spring-damper. We see from these figures that the trajectories of the end point in these two cases converge to the desired point smoothly. The final orientation of the wrist joint in the case of using a virtual spring-damper is more natural than that without one. of the muscle's output forces. We see from these figures that the output forces in the case of using a virtual springdamper are obviously reduced compared with those without a virtual spring-damper, even though both end points converge to the desired point smoothly. Figures 24 and 25 show the simulation results for reaching movements of all ranges. We see from these figures that when the desired position is in the extension direction from the initial position (especially F, G, H in Figure 25), there is excessive overshoot, which depends on the configuration of the muscles. However, this excessive overshoot can be reduced by introducing a virtual damper (see Figure 24). Figures 26 and 27 show the end-point velocities of all range-reaching movements with and without a virtual spring-damper. We see from these figures that the velocity profiles become more smooth when a virtual spring-damper is used. In paticular, F, G, H in Figure 26, the oscillation is decreased. These simulation results led us to suppose that reaching movements made before the movements were learned need a certain level of internal force to generate the joint damping effect because there is no adequate virtual damper at the end point. After learning, the internal forces can be reduced by obtaining pro-prioceptive information in the CNS, so an adequate virtual damper can be generated at the end point. However, in order to discuss this more, we need to compare a human's movements and our scheme.

Conclusion
This paper treated reaching movements by using the musculo-skeletal redundant arm model composed of three joints and nine muscles with physiological non-linear muscle properties to mimic a human arm. Asymptotic convergence of the closed-loop dynamic system was obtained based on the concept of stability theory on a manifold. We finally confirm through numerical simulation results that the end point of the musculo-skeletal redundant arm could converge to the desired point using only task-space feedback control with internal force adjustment under the existence of both joint and muscle redundancies, without any complex mathematical computation. In addition, we introduced the dynamic damping elipse to determine the internal force vector and evaluated it. We extended our control scheme to the virtual spring-damper hypothesis proposed by Arimoto et al. (2016). By introducing a virtual damper in addition to the virtual spring at the end point, we can reduce the internal forces. Moreover, the heterogeneity of convergence depending on the initial state and desired position is also reduced. These results lead us to the understanding that learning of human movements might be interpreted as acquisition of adequate coordination of a virtual spring and damper in the task coordinates. However, we have not treated the problem of learning of movements in this paper yet. Therefore, in order to claim that our control scheme is reasonable, we should gain a more physical insight into determination of the internal force vector k . In addition, we have only modeled co-contraction of agonist and antagonist muscles from the viewpoint of their configurations and contractile velocities, but did not present dynamics of muscle co-activation which might be given as a dynamic system of more complex differential equations. By introducing such a dynamic model of coactivation, the structure of the internal force space de-+ 9 e contraction would be disclosed in a more explicit form. Fortunately, since modeling of skeletal muscles is a hot topic in physiology (Houdijk et al. 2016;Perreault et al. 2013) and many new models will be proposed hereafter, we would like to introduce such skeletal muscle models into our control strategy for redundant joint and musculoskeletal systems. On the other hand, one of the advantages of our proposed method is that any planning of end-point trajectories is unnecessary in advance. Similarly, Todorov and Jordan (2012) have also proposed an optimal feedback control theory that does not need end-point trajectories.
Our control method presented here may be interpreted as one of such optimal feedback control schemes by regarding an artificial potential energy produced by the virtual spring as one of the possible cost functions. The equivalence of both methods will be treated in our future works. We would also like to introduce a learning strategy for the internal force to investigate human-like multi-body control mechanisms based on our control strategy. e scribed by (I − W W )k and the mechanisum of co-