Muscle Coordination Control for an Asymmetrically Antagonistic-Driven Musculoskeletal Robot Using Attractor Selection

Recently, numerous musculoskeletal robots have been developed to realize the flexibility and dexterity analogous to human beings and animals. However, because the arrangement of many actuators is complex, the design of the control system for the robot is difficult and challenging. We believe that control methods inspired by living things are important in the development of the control systems for musculoskeletal robots. In this study, we propose a muscle coordination control method using attractor selection, a biologically inspired search method, for an antagonistic-driven musculoskeletal robot in which various muscles (monoarticular muscles and a polyarticular muscle) are arranged asymmetrically. First, muscle coordination control models for the musculoskeletal robot are built using virtual antagonistic muscle structures with a virtually symmetric muscle arrangement. Next, the attractor selection is applied to the control model and subsequently applied to the previous control model without muscle coordination to compare the control model's performance. Finally, position control experiments are conducted, and the effectiveness of the proposed muscle coordination control and the virtual antagonistic muscle structure is evaluated.


Introduction
Human beings and animals move flexibly and dexterously by controlling their musculoskeletal system with the brain. To understand and imitate the flexible and dexterous motion of human beings and animals, several musculoskeletal robots have recently been developed. The primary driving mechanism in musculoskeletal robots is a tendon-driven assembly using motors [1,2] or pneumatic artificial muscles (PAMs). In particular, using PAMs as actuators enables flexible movement of the musculoskeletal robot compared with conventional robots that use motors because the compressibility and low viscosity of air provide compliance and rapid contraction. The extremely high power-to-weight ratio of the PAM is also good for flexible and dynamic motion.
The musculoskeletal robot comprises antagonisticdriven systems. An antagonistic-driven system includes two or more actuators for joint movement. The actuators are antagonistically arranged around one link, and their output characteristics and arrangements are, for the most part, symmetrical.
Many studies have proposed the musculoskeletal robots comprising this simple antagonistic system, and various control methods (e.g., PID control, neural network, and fuzzy logic) have also been proposed for musculoskeletal robots [3][4][5][6][7][8]. However, the drive system of our musculoskeletal robot [9] differs from a simple antagonistic-driven system because the output characteristics and the arrangements of the actuators of the robot are not symmetrical.
This robot has two kinds of PAM actuators, a monoarticular muscle that drives one joint and a polyarticular muscle that drives multiple joints consecutively. Therefore, the mechanism that drives each joint is not symmetrical. Furthermore, the actuators are not arranged symmetrically, although each actuator is antagonistically arranged around each link. Since each actuator is not arranged symmetrically and the sensors for measuring each joint angle are not arranged, design of the control system for the robot is more difficult and challenging.
Honda et al. [9,10] proposed a biologically inspired control method using muscle coordination for the robot. They hypothesized that human beings use the synergies between antagonistic muscle pairs when joints move, and they defined two parameters, the antagonistic muscle ratio and the antagonistic muscle activity, as the key parameters in human muscle coordination. The parameters are computed by a PID controller [10] and are implemented to control the angle of the robot's joints.
Although good control performance was obtained for one joint, the controller did not work well for multiple joints of the robot [9]. An adaptive method that dynamically and adaptively searches the parameters for muscle coordination was required because PAMs have time variance, compliance, high hysteresis, and nonlinearity. In general, this search problem can be formulated as a combinatorial optimization problem of minimizing an object function subject to the search variables, which requires a precise model (object function) in advance, but system identification against an asymmetrically antagonistic-driven PAM system having nonlinear dynamics is difficult.
As an alternative to such a model-based theoretical approach, we believe that heuristic control methods inspired from living things are important for the control of the musculoskeletal robot. The musculoskeletal structures of human beings and animals are antagonistic-driven systems. They are also asymmetrically antagonistic-driven systems because monoarticular muscles and polyarticular muscles are arranged around various joints. Human beings and animals move flexibly by controlling their various muscles dexterously.
Recent research about the mechanisms of living things indicates that a biological system behaves flexibly using noise [11]. Escherichia coli (E. coli) cells usually prefer to switch to an adaptive attractor using noise to survive better in a new external environment after the environmental conditions have been changed. This adaptive behavior of the E. coli cells is known as attractor selection [12]. The novel control method based on the attractor selection has been proposed and applied to a signal-control method for traffic networks [13], network management [14,15], android motion generation [16], robot navigation and locomotion [17][18][19], robotic arm control [20][21][22], and endoscopic surgery [23].
The control method was described by a stochastic differential equation, input variables for the network systems or robots were computed by solving the equation, and the systems accomplished the tasks without the dynamics and model of the systems and environments. Since the attractor selection is conducted adaptively using noise in the systems or environments, the control method is robust for changes of tasks and environments. Attractor selection was applied to an asymmetrically antagonistic-driven musculoskeletal robot ( Figure 1) with muscles arranged asymmetrically [24]. From the control experiment, the position of the tip of the robot was moved to the desired position by searching pressure for four PAMs individually using the attractor selection. The control time had to be more than 100 s to accomplish tasks. Therefore, modification of the control method is required to accomplish tasks quickly.
In this study, we propose a novel muscle coordination control method for the asymmetrically antagonistic-driven musculoskeletal robot using the attractor selection. The primary difference between the previous method [24] and the proposed one is that the proposed method introduced a virtual antagonistic muscle structure as a muscle coordination control model. Instead of individually and directly searching the PAM pressure in the actual asymmetric antagonistic muscle structure, the new method indirectly searches pressure for actual PAMs via a virtual symmetric antagonistic muscle structure.
First, muscle coordination control models of the musculoskeletal robot were built using virtual antagonistic muscle structures with a virtually symmetric arrangement of muscles. Next, the attractor selection was applied to the control model and also applied to the previous control model without the muscle coordination to compare control performance. Finally, position control experiments were conducted, and the effectiveness of the proposed muscle coordination control applied attractor selection and the virtual antagonistic muscle structure was evaluated.

Purpose: Control of the Position of an Asymmetrically
Antagonistic-Driven Musculoskeletal Robot. In this study, a five-fingered robot hand (Figure 1(a)) inspired by the right hand of a human being was used, and the position of the tip of the index finger of the robot hand was controlled as a musculoskeletal robot. Figure 1(b) shows the structure of the index finger. It has a three-degrees-of-freedom mechanism that is formed by four links (a distal phalange, an intermediate phalange, a proximal phalange, and a metacarpal) and three joints (distal interphalangeal (DIP), proximal interphalangeal (PIP), and metacarpophalangeal (MP)). Four PAMs are used: one actuator on the intermediate phalange, one on the proximal phalange, and two on the metacarpal for flexion and extension. The actuator for extension (referred to as the "extensor") is a polyarticular muscle in which wire covers all joints. Therefore, the number of actuators for flexion and extension is not symmetrical.
Let P e , P MPf , P PIPf , and P DIPf be the pressure supplied to the muscles A, B, C, and D, respectively, in Figure 1. The subscript "e" indicates the extensor, "MPf" indicates the flexor (the actuator for flexion) for driving the MP joint, "PIPf" indicates the flexor for driving the PIP joint, and "DIPf" indicates the flexor for driving the DIP joint, and they are calculated as follows: where x e ∈ 0, 1 , x MPf ∈ 0, 1 , x PIPf ∈ 0, 1 , and x DIPf ∈ 0, 1 are the normalized search variables and P max is the maximum pressure supplied to the actuator. Since the actuators of the robot hand break if a pressure of more than 0.2 MPa is supplied, the value of P max is set to 0.19 MPa.

Muscle Coordination Hypothesis.
Honda et al. suggested the hypothesis that muscle coordination of human beings is a coordination of antagonistic muscles (the extensor and flexor) [9,10]. They express the ratio of the coordination using two parameters, the antagonistic muscle ratio (Ar) and antagonistic muscle activity (Ac). The Ar is the value that regulates the ratio of pressure (P e and P f ) between the antagonistic muscle e (the extensor) and f (the flexor). Ar is calculated between 0 and 1 and is defined as follows: Ac is the sum of pressure for the extensor e and the flexor f and is calculated by In this study, Ac is always set to the maximum pressure for driving joints sufficiently. That is, Pressures P e and P f are calculated from (5), (6), and (7).
Equations (8) and (9) show that the pressures of the extensor and the flexor are determined by searching for the normalized variable Ar ∈ 0, 1 . Figure 2(a) shows the principle of muscle coordination using the antagonistic muscle ratio Ar and antagonistic muscle activity Ac.

Two Types of Virtual Antagonistic Muscle Structures.
Four actuators were used in the musculoskeletal robot: the polyarticular muscle was used as the extensor and three monoarticular muscles were used as the flexor (Figure 1 To calculate the pressure for each actuator based on the muscle coordination hypothesis, two methods that make the total number of extensors and flexors match virtually are proposed. The first method is composed as follows. The virtual antagonistic muscle structure (Figure 2(b)) is composed so that the number of extensors is increased from one to three extensors. Next, Ar is applied to each antagonistic muscle to drive each joint. Hence, Ar is described as Ar MP , Ar PIP , and Ar DIP , and the pressures for the extensor and the flexor at each antagonistic muscle for driving each joint are denoted as P MPe , P MPf , P PIPe , P PIPf , P DIPe , and P DIPf . The pressures are calculated as follows: P MPe = P max · Ar MP , 10 P PIPe = P max · Ar PIP , 11 P DIPe = P max · Ar DIP , 12 Finally, P e that applies to the real extensor on the real form ( Figure 2(d)) is calculated using The second method is distinguishable in the following ways. First, the other virtual antagonistic muscle structure ( Figure 2(c)) reduces the number of flexors to one. Next, Ar is applied to the virtual antagonistic muscle structure, and the pressures P e for the virtual extensor and P f for the virtual flexor are calculated using (8) and (9). Here, P f is the total pressure for the three real flexors (P MPf , P PIPf , and P DIPf ) and is defined as P f ≤ P max . Finally, P f is distributed to the three real flexors on the real form ( Figure 2(d)). Two distribution ratios (Dr MP and Dr PIP ) are used, and the pressure for the three real flexors is calculated by Pressure for the real extensor is the same as that for the virtual extensor and is calculated using (8).

Attractor Selection
Model. Let us consider the combinatorial optimization problem of minimizing an object function U x subject to the search variable x ∈ X 1 , X 2 , … , X N , where X i is a feasible solution (an attractor) and N means the number of attractors. Rather than seek an optimal solution, we try to quickly find good approximate solutions. To accomplish this, our study uses the revised attractor selection model, which has been generalized as a stochastic differential equation [24]: where t is time, x is the search variable or state ( ∈ 0, 1 for our case), the value Activity ∈ 0, 1 is the degree of accomplishment of the task, η is assumed as noise, and f x is the function that makes x converge to a suitable attractor. Typically, the function f can be represented as f x = −∂U x /∂x if the objective function U x is known precisely in advance. This model searches for a solution (the attractor) that successfully accomplishes the task using noise, and the Acti vity makes the behavior of the total system change. Notice that as Activity increases, the term f x · Activity becomes more dominant in (20) and the state transition becomes more deterministic. Consequently, state x tends to be entrained into a suitable attractor, where it remains despite the persistent noise. By contrast, decreasing the Activity increases the dominance of the noise η, thereby flattening the potential landscape. In this scenario, the transition becomes more probabilistic, like a random walk, and x is driven away from the attractor.
The function f x can be designed freely, even if the objective function U x is unknown or not precisely described. Two elements are required. The value of x must converge to the attractors, and the x must remain at a suitable attractor. To satisfy these elements, f x is defined as follows: where X i is the ith attractor, N is the number of attractors, k d is the power that attracts x, and k w is the range in which the attractor's power k d is effective.
2.5. Employing the Attractor Selection to Determine the Pressure Supplied to Each Actuator. The selection of the attractor is employed to determine the pressure supplied to each actuator. The three methods used to determine the pressure using attractor selection are presented as follows.
The first method uses the real musculoskeletal structure (Figure 2(d)) and directly calculates the four pressures by searching for the four variables (x e , x MPf , x PIPf , and x DIPf ) using the attractor selection. We refer to this as a pressure search-type controller. The pressure supplied to each actuator is calculated by (1), (2), (3), and (4). This method is the same as reported in our previous method in [24].
The second method uses the first virtual antagonistic muscle structure (Figure 2(b)). It searches for the three ratios (Ar MP , Ar PIP , and Ar DIP ) using the attractor selection and supplies the pressure to each actuator using (10), (11), (12), (13), (14), (15), and (16). We refer to this as an Ar searchtype controller.
The third method uses the second virtual antagonistic muscle structure (Figure 2(c)). It searches for the Ar value and the two distribution ratios (Dr MP and Dr PIP ) using the attractor selection and then supplies the pressure to each actuator using (8), (9), (17), (18), and (19). This is referred to as an Ar and Dr search-type controller. Ac. When Ar is 0 (left picture), the pressure for the flexor is high (P f = Ac) and the extensor is low (P e = 0). When Ar is 0.5 (middle picture), the pressure for the flexor and the extensor is identical (P f = P e = Ac/2). When Ar is 1 (right picture), the pressure for the flexor is low (P f = 0) and the extensor is high (P e = Ac). Therefore, the pressures of the extensor and the flexor are determined by searching for the normalized variable Ar ∈ 0, 1 . (b) One of the virtual antagonistic structures. Two muscles are added virtually to the musculoskeletal system (shown in (d)) to form a symmetrically antagonistic arrangement on the intermediate phalange and the proximal phalange. Therefore, the virtual antagonistic structure has three simple antagonistic muscle structures (inset picture in (b)). (c) The other virtual antagonistic muscle structure. Two muscles on the intermediate phalange and the proximal phalange are decreased virtually from the musculoskeletal robot (shown in (d)), and the muscles are symmetrically arranged on the metacarpal. Therefore, the virtual antagonistic structure has one simple antagonistic muscle structure (inset picture in (c)). (d) The real form of the musculoskeletal robot that is controlled. The value of Ar (shown in (a)) is applied to the virtual antagonistic structure (shown in (b) or (c)) virtually transformed from the musculoskeletal robot (shown in (d)).
Notice that the search space for all the variables (x e , x MPf , x PIPf , and x DIPf ; Ar MP , Ar PIP , and Ar DIP ; and Ar, Dr MP , and Dr PIP ) ranges from 0 to 1.

Experimental Procedures
A control experiment that makes the tip of the index finger of the musculoskeletal robot (see the bottom-right inset in Figure 3) move to the desired position was conducted using the above three controllers. Figure 3 depicts the experiment setup, which comprises the musculoskeletal robot (SQUSE hand G type, SQUSE Inc.), the control PC (MDV ADVANCE ST 6300B (MouseComputer Co. Ltd.), Windows XP, and an Intel Core i7 920 (2.67 GHz)), an A/D converter (AI-1664L-LPE, CONTEC Co. Ltd.), two D/A converters (AO-1616L-LPE, CONTEC Co. Ltd.), a digital output board (RRY-32-PE, CONTEC Co. Ltd.), a motion capture system (Nobby Tech. Ltd.), regulators (ITV0030, SMC Corporation), solenoid valves (S070-5DCO-32, SMC Corporation), and an air compressor (DPP-AYAD, Koganei Corporation). The sampling frequency was set at 100 Hz. The input signals were voltages generated by the control PC and were converted to pressures by the regulators. The output signals were the coordinates of the tip position [E X , E Y , E Z ] sensed by the motion capture system. The Euclidean norm (the distance between the desired position [E Xd , E Yd , E Zd ] and the tip position) was used as the evaluation index of the controller. The Euclidean norm is described by l and computed using And Activity of the attractor selection model is calculated from 0 to 1 using Activity = − l l max + 1 23 Here, l max is the maximum value of a norm computed from the desired position and the tip position on either the maximum flexion, which is a steady state in which maximum pressure (0.19 MPa) is supplied to all flexors, or the maximum extension, which is a steady state in which maximum pressure is supplied to only the extensor of the robot.
Each norm was obtained in advance, and the larger norm was selected as l max . Two tasks were conducted in the experiment. The flexion task involved making the robot flex to a desired position from an extended state, the extension task involved extending the robot to a desired position from a flexed state, and the tasks were changed after a constant time. First, the flexion task was conducted. After a constant time, the desired position was changed and the extension task was conducted. The control time was set at 60 s, and the task was changed 30 s after the control was started. The first position and the desired position were defined when the robot was in a steady state after constant pressure was applied to each actuator. Each position was captured in advance.
In this experiment, the pressure value, represented by P e , P MPf , P PIPf , and P DIPf (0.05, 0.15, 0.05, and 0.05, resp.), was applied to each actuator to determine the desired position for the flexion task, and 0.15, 0.05, 0.05, and 0.05, respectively, were applied to each actuator to determine the desired position for the extension task. The initial values of x e , x MPf , x PIPf , and x DIPf ; Ar MP , Ar PIP , and Ar DIP ; and Ar, Dr MP , and Dr PIP were set to 0.9, 0.1, 0.1, and 0.1; 0.9, 1.0, and 0.0; and 0.9, 1.0, and 0.0, respectively, and the noise η was generated between −10 and 10. The parameter values in (21) were set as follows: N = 11, k d = 0 01, k w = 0 01, and X i = 0 1 × i.  show the results using the pressure search-type controller, the Ar search-type controller, and the Ar and Dr search-type controller, respectively. The transition of the search variables, the pressures supplied to each actuator, the tip position captured by the motion capture, and Activity of the attractor selection model were then plotted. In Figure 4, Activity increased from 3 s and became constant at 0.96. Therefore, the flexion task was almost accomplished, and the search variables converged to an attractor. However, the extension task was not completed. Activity increased from 33 s to 35 s and remained at nearly 0.8, but it decreased from 39 s and did not reach a high value. The difference in the accomplishment ratio of the task is caused by the relationship of pressure between the extensor and the flexor. When Activity increased and became constant in the flexion task, the pressure on the extensor (P e ) decreased and the pressure on the flexor responsible for driving the MP joint (P MPf ) increased. Therefore, the power for the flexion became large, and the robot flexed toward the desired position. In the extension task, pressures for the extensor (P e ) and pressure for the flexor for driving the MP joint (P MPf ) increased from 33 s to 35 s. Therefore, the power for the extension did not increase dramatically, The input voltage from the control PC is converted to the controlled pressures by the regulators, and the musculoskeletal robot is driven by the pressures supplied to each muscle of the robot. The position of the tip is measured using the marker captured by the cameras and is saved to the PC. and the robot did not accomplish the extension task. The difference between the pressure for the extensor and the flexor is important for achieving the task efficiently.

Results and Discussion
In Figure 5, Activity increased after 10 s and reached a constant value of 0.91 at 17 s. The flexion task was close to being accomplished, but compared with the pressure search-type controller, the accomplishment ratio was 5% lower and the control time for accomplishment of the task was 13 s longer. Therefore, the control was not as good as that of the pressure search-type controller. In the extension task, Activity did not attain a high value. Therefore, the extension task was not accomplished. As a result, the controller did not work well for the extension task. To efficiently accomplish the task, it was necessary that the differences between the pressures for the extensor and the flexors be easily calculated by the controller because the differences can decrease power which prevents accomplishment of tasks and can make the tip position of the robot move quickly to the desired position, but this controller cannot calculate the differences between the pressures as well as the pressure search-type controller. This controller searches for three Ar values independently. Therefore, each Ar takes on a different value, and the pressure for the flexor does not decrease satisfactorily. Thus, the power for flexion does not decrease sufficiently, and the extension task is not accomplished.
In Figure 6, Activity increased greatly and became constant at 0.95 for the flexion task. On the extension task, Activ ity increased gradually right after the change of the task and became constant at 0.97. Therefore, the flexion task and the extension task were accomplished. The Ar and Dr search- type controller is, thus, more highly adaptive for tasks than the pressure search-type controller and the Ar search-type controller. Especially, the effectiveness of the Ar and Dr search-type controller for the extension task was shown. The robot has four actuators arranged asymmetrically on the extension side and the flexion side. One actuator is installed on the extension side, and three actuators are installed on the flexion side. Therefore, the power for flexion easily becomes large compared with the power for extension. To make this robot extend sufficiently, the power for flexion must become small. In the Ar and Dr search-type controller, the difference between the power in the flexion and in the extension is determined easily because the difference between the pressures for the flexor and the extensor is calculated by searching only one Ar. Therefore, the power for flexion becomes small, and the robot can sufficiently extend. We also conducted the different task for the Ar and Dr search-type controller to show the adaptability of the controller. Figure 7 shows the results of the experiment. We set four tasks, two flexion tasks and two extension tasks, which were conducted alternately. The desired position for each flexion task was different, and each desired position for each extension task was different. Each task was changed at the following constant times: 30 s, 60 s, and 90 s after the control is started. Experimental results show that Activity became the high value for each task. Therefore, the Ar and Dr search-type controller displayed good adaptability for the tasks in the experiments. Thus, it was shown that the proposed muscle coordination control of the Ar and Dr search-type controller using the attractor selection facilitated easy control of an asymmetrically antagonistic-driven musculoskeletal robot with a polyarticular muscle.

Conclusions
This work demonstrated a muscle coordination control of an asymmetrically antagonistic-driven musculoskeletal robot using attractor selection which is a biologically inspired search method. First, muscle coordination control models of the musculoskeletal robot were built using virtual antagonistic muscle structures with a virtually symmetric arrangement of muscles, and the calculation methods of the input pressure for PAMs of the musculoskeletal robot with and without muscle coordination were shown. Next, the attractor selection was applied to both the muscle coordination control model and to another control model without the muscle coordination to compare the control performance. Finally, position control experiments were conducted, the effectiveness of the proposed muscle coordination control applied to the attractor selection was demonstrated, and it was also shown to be faster and more robust to accomplish the task by generating control commands virtually assuming a symmetrical and simpler metastructure rather than providing a control command according to the actual complex (asymmetric) antagonistic muscle structure.
Based on the virtual antagonistic muscle structure proposed in this research, we may be able to build a musculoskeletal robot that achieves a more complicated task faster by devising how to give noise [25] and adaptively updating the attractor structure [21]. In future work, the muscle coordination control method, using the attractor selection, will be applied to the multifingered robot hand formed by increasing the number of musculoskeletal robots (i.e., robot fingers). Furthermore, the effectiveness of the control method will be investigated for asymmetrically antagonistic-driven musculoskeletal robots, which have entirely different arrangements of muscles compared with our musculoskeletal robot.
The muscle coordination control method using attractor selection can be applied not only to musculoskeletal robots but also to human hands. Therefore, the control method will be applied to rehabilitation using functional electrical stimulations (FESs) [26,27] as a novel approach to controlling human hands.

Data Availability
The raw data for Figures 4-7 are available from the corresponding author upon request. Plots of the transition of the input pressures (P e , P MPf , P PIPf , and P DIPf ). In (c), the top images plot the trajectory of the tip of the robot: the first flexion task, the first extension task, the second flexion task, and the second extension task from the left image, whereas the bottom image plots the transition of the Activity. Cyan arrows and black arrows show the direction of the movement and the value of Activity, respectively.