Simulation Evaluation for Methods Used to Determine Muscular Internal Force Based on Joint Stiffness Using Muscular Internal Force Feedforward Controller for Musculoskeletal System

This study proposes two novel methods for determining the muscular internal force (MIF) based on joint stiffness, using an MIF feedforward controller for the musculoskeletal system. The controller was developed in a previous study, where we found that it could be applied to achieve any desired end-point position without the use of sensors, by providing the MIF as a feedforward input to individual muscles. However, achieving motion with good response and low stiffness using the system, posed a challenge. Furthermore, the controller was subject to an ill-posed problem, where the input could not be uniquely determined. We propose two methods to improve the control performance of this controller. The first method involves determining a MIF that can independently control the response and stiffness at a desired position, and the second method involves the definition of an arbitrary vector that describes the stiffnesses at the initial and desired positions to uniquely determine the MIF balance at each position. The numerical simulation results reported in this study demonstrate the effectiveness of both proposed methods.


INTRODUCTION
Humans are able to easily achieve specific complex movements by manipulating the structure of the human body. The human body structure can be described in terms of a musculoskeletal structure in which muscles, corresponding to actuators, cover bones and joints. Joint angle and stiffness can be adjusted by contracting the muscles, thus allowing the individual to perform complex motions.
A number of robots that can imitate this musculoskeletal structure have recently been investigated. Marques et al. (Marques et al., 2010) developed a human coexistence robot that imitates the musculoskeletal structure, while a research group including Inaba (Mizuuchi et al., 2006;Asano et al., 2016) developed another musculoskeletal structure-imitating humanoid robot. Verrelst et al. (Verrelst et al., 2005) developed a biped robot actuated with antagonistic pneumatic artificial muscles. This robot can control joint stiffness and achieve stable walking. Niiyama et al. (Niiyama et al., 2012) developed a double legged robot with a musculoskeletal structure that emulates an athlete's physique and style of running. Their robot had the capability of running stably for 4 m. Niiyama et al. (Niiyama et al., 2007) developed a jumping robot with a musculoskeletal structure that could perform big leaps. The results of these studies indicate that the musculoskeletal structure functions as a feedback system when undergoing instantaneous movements; thus, an understanding of the musculoskeletal structure can be used to develop robots that can perform previously difficult-to-achieve motions.
In a previous study, we proposed a muscular internal force (MIF) feedforward controller for a musculoskeletal system (Kino et al., , 2013. This controller could be used to attain any desired position without the use of sensors, by providing a MIF as a feedforward input to the associated muscle. The controller input would create a potential field in the musculoskeletal system; allowing the system to converge at the desired position when this potential field achieved a stable equilibrium point at the desired position. This end-point convergence is influenced by muscle arrangement . The MIF must be increased to obtain a desirable control response, through increase in system stiffness. Subsequently, reducing the MIF to lower the stiffness causes the responsiveness to deteriorate. As a result of the correlation between response and stiffness, it is challenging to simultaneously achieve high-response motion and low system stiffness. Furthermore, the controller is affected by an ill-posed problem where, its input cannot be uniquely determined due to the muscle redundancy of the musculoskeletal system.
To tackle this muscle redundancy, Buchler et al. (Buchler et al., 2018) proposed a variance-regularized control that can achieve high accuracy trajectory tracking in the musculoskeletal system. In addition, Jantsch et al. (Jantsch et al., 2011) proposed a new controller based on computed torque control and a method for obtaining the muscle Jacobian by using machine-learning techniques for the musculoskeletal system. Their method indicates that trajectory control and force control can be achieved. Katayama and Kawato (Katayama and Kawato, 1990) proposed a parallelhierarchical neural network model based on a feedback-errorlearning scheme, in the field of exercise physiology. This model requires repeated trials; however, it is capable of achieving rapid movements that have otherwise been difficult to achieve using a feedback controller. A tendon driven robot equivalent to the musculoskeletal system has been explored (Viau et al., 2015). However, position control methods, like the feedforward controller method from our previous study (Kino et al., , 2013, have not been fully explored. Further, various studies have been conducted on methods for contributive employment of redundancies (Klein and Huang, 1983). For instance, it has been observed that redundancy can be used to avoid specific robot attitudes (Mayorga and Wong, 1988) and obstacles (Baillieul, 1986;Chirikjian and Burdick, 1990) during motion. In addition, redundancies can be used to optimize torques in the robotic structure (Suh and Hollerbach, 1987), avoid restrictive joint movement (Shimizu et al., 2008), and apply stiffness (Svinin et al., 2002) and impedance control (Albu-Schaffer et al., 2003). Yoshikawa (Yoshikawa, 1985) suggested that the operability of a manipulator can be improved using redundancy; while, Hanafusa et al. (Hanafusa et al., 1981) and Nakamura et al. (Nakamura and Hanafusa, 1987) proposed control methods for achieving sub-tasks, based on the prioritization of robot actions. In parallel wire-driven systems (Sui and Zhao, 2004;Behzadipour and Khajepour, 2006), system actuation redundancy can be used to generate internal wire forces. The stiffness of the system can be changed by controlling these internal forces, thus indicating that robot performance and versatility can be improved by manipulating the robot's redundancy. Hence, redundancy has several potential advantages. Meanwhile, sensors are mostly utilized in methods that employ redundancy. The musculoskeletal system continues to encounter muscle redundancy. The musculoskeletal system can achieve various kinds of motion by utilizing muscle redundancy without the addition of a new actuator.
The primary objective of this study is, to resolve the issues in the MIF feedforward controller proposed in previous studies. We began by eliminating the foremost issue of correlation between the response and stiffness. Further, the MIF that separately sets the response and stiffness was determined, and the ill-posed problem was solved. In this study, a new generic method for determining MIF using joint stiffness was developed, as an approach to improving the control performance of our previously developed MIF feedforward controller. The proposed method involves the application of an approach derived from Nakamura (Nakamura and Hanafusa, 1987), which includes prioritized application of subtasks to increase the stiffness of joints at desired positions by way of redundancy. This approach was used to develop two specific MIF determination methods: one involving the independent setting of response and stiffness, and the other involving the unique determination of MIF. The first method was used to identify the MIF that can independently determine the response and stiffness at a desired position, while the second was used to determine an arbitrary vector that reflects the stiffness at the initial and desired positions, and uniquely defines the balance of MIFs at each position.
The remainder of this paper is organized as follows. Section 2 describes the kinematics of the musculoskeletal system. Section 3 describes the MIF feedforward system used to control the musculoskeletal system. In Section 4, describes the development of a joint stiffness matrix of the musculoskeletal system. The two methods and their simulation results are presented in Sections 5 and Section 6, respectively.

MUSCULOSKELETAL SYSTEM
The musculoskeletal system considered in this study is shown in Figure 1. The system comprised of two links and six muscles. The joint had one degree of freedom and could rotate when muscle tension was applied through the link. A joint with one degree of freedom requires at least two muscles to drive its link in the clockwise and counter-clockwise directions because muscles can only produce shrinkage forces; as a result of which, this system has an actuation redundancy. Additionally, we assumed that the movement of the system was unaffected by the force of gravity because the system was set in the x-y plane.
The muscle consisted of a combination of wire and motor instead of a pneumatic actuator . The muscles were connected to an anchored base, and to each link. The anchoring was set at an explicit offset that influenced the convergence of the MIF feedforward controller (Kino et al., , 2013. In our analysis of the muscular arrangement, we employed a set of parameters that converged at a desired position. Although, as mentioned above, the influence of gravity is also not considered; however, the muscular internal force feedforward controller is able to achieve system movement under the influence of gravity from a result of an other study (Matsutani et al., 2019b).

Relation Between Muscle and Joint Spaces
The relation between the muscle length vector q θ q 1 , . . . , q 6 T and the joint angle vector θ θ 1 , θ 2 T within the musculoskeletal system is given by where L i (i 1, 2) is the link length, a j (j 1, . . . , 4), b j , d j , h j , u j , s j are the muscular arrangement parameters, as shown in Figure 1, C i and S i are defined in terms of i as C i cos θ i and S i sin θ i , respectively, and C 12 and S 12 are defined as C 12 cos (θ 1 + θ 2 ) and S 12 sin (θ 1 + θ 2 ), respectively. By differentiating Eq. 1 with respect to time, the relation between the muscle contractile vector _ q θ ∈ R 6 and the joint angular velocity vector _ θ ∈ R 2 can be obtained: where W θ ∈ R 2×6 is the Jacobi matrix showing the relation between the respective vectors, which is given by The relation between the joint torque vector τ θ τ 1 θ , τ 2 θ T and the muscle tension vector α θ ∈ R 6 can be obtained based on the principle of virtual work as The inverse relationship to that given in Eq. 4 is expressed as where I 6 ∈ R 6×6 is a unit matrix, k e ∈ R 6 is an arbitrary vector, v θ ∈ R 6 is the MIF vector, and W θ The MIF vector v(θ) corresponds to the internal force applied within the muscle. This vector belongs to the null space, that is, it is a muscle tension component that cannot generate joint torque and satisfies the following equation As a result, the musculoskeletal system was subjected to an illposed problem because the MIF vector corresponding to a given joint torque could not be uniquely determined.

Relation Between Task and Joint Spaces
The relation between the end-point position vector x(θ) [x,y] T and the joint angle vector θ of the musculoskeletal system is given by By differentiating Eq. 9 with respect to time, the relation between the end-point velocity vector _ x θ ∈ R 2 and the joint angular velocity vector _ θ can be obtained as Where J θ ∈ R 2×2 is the Jacobi matrix showing the relation between the respective vectors.

MUSCULAR INTERNAL FORCE FEEDFORWARD CONTROLLER
The MIF feedforward controller, developed in our previous study (Kino et al., , 2013) was used to positionally control the musculoskeletal system. In this method, the controller inputs a balancing MIF at a desired joint angle of θ d into the muscles of the musculoskeletal system. This control input has been expressed using Eqs 5, 6 as Under this method, control was obtained using only kinematics; the input of any physical parameters or sensor information of the musculoskeletal system was not required. The Jacobi matrix W (θ d ) had a constant value at all times, across all joint angular displacements because the matrix was obtained using the desired joint angle θ d . Additionally, the control input was constant when the arbitrary vector k e was constant.
When the shrinkage direction of the muscle was set at a positive, the muscle tension α(θ) also had a positive value because the muscle can produce only tension or shrinkage. Therefore, the elements v dl (l 1, . . . , 6) of the Under this control method, the arbitrary vector k e was set freely within a range satisfying Eq. 12.
Joint torque cannot be generated at a joint angle of θ θ d to the muscle, by applying the control input detailed in Eq. 11: The control input is a muscle tension component that does not generate joint torque because the MIF vector v (θ d ) belongs to the null space of the Jacobi matrix W (θ d ). By contrast, at joint angles of θ ≠ θ d the control input v (θ d ) cannot be considered to be a vector belonging to the null space of W(θ). It is possible for a control input v (θ d ) to generate a joint torque because v (θ d ) can be decomposed into two components: a muscle tension component that does not generate joint torque and the vector v(θ) that belongs to the null space: Using the abovementioned characteristics, the MIF controller generated joint movement.
The convergence of the motion end-point was related to the potential P generated by the MIF v (θ d ) (Kino et al., , 2013, which is defined as A musculoskeletal system movement converged on a desired joint angle θ d when the angle dependent on muscular arrangement, corresponded to a stable equilibrium point in P (Kino et al., , 2013. Moreover, the muscular arrangements that were used to attain convergence were identified based on the results of our previous studies (Kino et al., , 2013

JOINT STIFFNESS MATRIX
The joint stiffness matrix K j θ ∈ R 2×2 of the musculoskeletal system (Nakiri and Kino, 2009) is given by Here, because the joint torque τ(θ) was generated by the MIF controller, the joint stiffness matrix was derived from Eqs 4, 11 as follows: For convenience, the Jacobi matrix W(θ) from Eq. 3 is rewritten as The joint stiffness matrix in Eq. 17 can be re-expressed in terms of its matrix elements as Equation 19 can then be substituted into Eq. 16 to obtain Also, from Eq. 18, The joint stiffness matrix K j (θ) is a diagonal matrix for which the elements can be obtained using from which the following is derived as In this study, Eq. 23 was used to model two methods of MIF determination. Using the first method, an MIF that eliminates the correlation between the response and stiffness at the desired position was determined. Using the second method, an arbitrary vector describing the stiffness at the initial and desired positions to uniquely determine the balance of MIFs at a given position, was defined.

METHOD FOR DETERMINING MUSCULAR INTERNAL FORCE USING JOINT STIFFNESS AT THE DESIRED POSITION
In this section, a novel method for determining the MIF, which can eliminate the correlation between the response and stiffness at a desired position, is proposed. Under this method (Method 1), the MIF at the desired position was determined using the joint stiffness.
At a joint angle of θ θ d , the matrix in Eq. 23 becomes Where The arbitrary vector k e is obtained as the inverse of Eq. 24: where the matrix Z θ d + ∈ R 6×3 is a pseudo inverse matrix of the matrix Z (θ d ) defined by where β ∈ R 6 is an arbitrary vector that was set freely within the range satisfying Eq. 12. A block diagram of the application of the MIF controller using Method 1 is shown in Figure 2.

Simulation Results
The simulated positioning control of the musculoskeletal system shown in Figure 1 was conducted to assess the ability of Method 1. This simulation shows that the MIF can eliminate the correlation between the response and stiffness at the desired position using Method 1. In these simulations, an end-point stiffness ellipsoid (Tsuji et al., 1995) was used to obtain a visual representation of the stiffness. Under this approach, the desired stiffness was represented as the end-point stiffness K e x θ ∈ R 2×2 , and a joint stiffness K j (θ) corresponding to the end-point stiffness was calculated. The relation between the end-point and joint stiffnesses is given by The MIF can be determined by substituting the elements of the desired joint stiffness, K 11 (θ d ), K 12 (θ d ), and K 22 (θ d ), calculated from the end-point stiffnesses K e (x (θ d )), respectively, into Eq. 26.
The results of positioning control using each input determined by three desired end-point stiffness matrixes K e1 (x (θ d )), K e2 (x (θ d )), and K e3 (x (θ d )) were compared in this simulation. The simulation employed a conventional method (CM) (Kino et al., 2013), which minimized the norm of the MIF as a comparison target. The arbitrary vector using CM is given by The result of the simulation performed by setting c 70 is expressed as CM1 and that for c 115 is expressed as CM2. Table 1 lists the physical parameters and muscular arrangements of the musculoskeletal system, respectively, and Table 2 lists the initial and desired parameters. The arbitrary vector β was set as follows:  The results of the simulations are shown in Figure 3 ∼ 4. Figure 3A shows the transient responses of the end-point position; Figure 3B compares the control inputs; Figure 3C ∼ Figure 3F show the end-point stiffness ellipsoids and loci of the end-point positions in the task space, respectively; Figure 4A compares the end-point stiffness ellipsoids; Figure 4B compares the end-point stiffness ellipsoids of CM2 and K e2 (x (θ d )); Figure 4C compares the loci of the end-point position in the   task space of CM2 and K e2 (x (θ d )). From Figure 3A, Figure 3C ∼ Figure 3F, it is seen that the end-point positions converge at the desired position. The end-point stiffness ellipsoids corresponding to 0.1-s intervals on the end-point trajectories are shown in Figure 3C ∼ Figure 3F. In Figure 3D and Figure 3F, the major axis of the ellipsoid is shortened and the end-point trajectory bulges outward. In Figure 3E, the major axis of the ellipsoid is lengthened and the end-point trajectory suppresses the outward bulge. Although Figure 3F shows an end-point trajectory that is equivalent to the trajectory in Figure 3D, the end-point stiffness ellipsoids at the desired positions differ and the control inputs v d3 , v d4 are significantly reduced, as shown in Figure 3B. Those results show that Method 1 can achieve low stiffness without sacrificing the response. A comparison of the stiffness ellipsoids K eR1 , K eR2 , K eR3 obtained using the target end-point stiffnesses K e1 , K e2 , K e3 , respectively, is shown in Figure 4A. Figure 4A also shows the stiffness ellipsoids of CM1. In addition, joint stiffness vectors K r which were used in determining the MIF, and joint  stiffness vectors K r which were realized by control, are shown in Table 2. These results depict that the size of the ellipsoid changed, and the desired end-point stiffness used for determining the arbitrary vector k e was achieved. However, the ellipsoids showed almost no change in the direction of the major axis. As the direction of the major axis is dependent on the posture of the musculoskeletal system, it is difficult to have a significant impact through the controller input. However, when the arbitrary vector c was set by constructing an end-point stiffness ellipsoid that had the same length as the major axis at the K e2 ellipsoid, CM2 incurred an overshoot, whereas high stiffness was realized, as shown in Figure 4B and Figure 4C. Method 1 can realize high stiffness without causing an overshoot and can eliminate the correlation between the response and stiffness at the desired position.

Evaluation of End-point Stiffness Ellipsoid
The characteristics of the ellipsoid were evaluated by changing the direction of the major axis using trial and error, to uncover the relation between the end-point stiffness K e (x (θ d )) and axis direction, because, as discussed in the preceding section, it is difficult to significantly change the direction of the major axis of the end-point stiffness ellipsoid using control input. The stiffness matrix at the time of inputting the MIF, determined using a conventional method, was used as reference. The elements of the standardized stiffness matrix were changed to rotate the longest axis of the ellipsoid in the clockwise direction. This change was made using trial and error because the relationship between the elements and rotational direction remains unidentified. The MIF calculated based on the stiffness matrix indicates that it satisfied Eq. 12 and was less than 500. If the conditions are satisfied, same elements are changed; otherwise, other elements are modified to rotate the axis in a counter-clockwise direction. Furthermore, the axis length was not evaluated because the lengths of the major and minor axes of the ellipsoid were altered simultaneously while maintaining a constant length ratio. The response of the end-point position was also not considered. The results of varying the direction of the major axis of the end-point stiffness ellipsoid are shown in Figure 5 and listed in Table 3. It can be observed through the results that the parameters were changed within a range at par with the norm of MIF being less than 500. Figure 5  These results indicate that the direction of the end-point stiffness ellipsoid generated by the application of Method 1 is dependent on the posture of the musculoskeletal system. In this case, the relation between the stiffness matrix parameter K e (x (θ d )) and the major axis direction cannot be clarified due to position difference.

Evaluation of Modeling Error
Next, the influence of the musculoskeletal system, including modeling error, was evaluated. The modeling error that occurred during the application of Method 1 to an actual musculoskeletal system could not be ignored. However, the error was justified, as it was difficult to model the actual system accurately. Therefore, the influence of muscular arrangements, including modeling error, was evaluated. In this evaluation, it was assumed that errors in muscular arrangements occurred because of an assembly error in the musculoskeletal system. The proposed method was evaluated as a true value using the parameters listed in Table 1. The evaluation function of the end-point position and joint stiffness vector is defined as where x f is an end-point position of a simulation finish time and K rf is the joint stiffness vector at x f .  One cause for the error in muscular arrangements is the assembly error of the musculoskeletal system. In this case, it was considered that large errors do not occur in the measurement values because muscular arrangements can be easily measured. Therefore, randomly selected values from a range of error rates of 5% were used as actual muscular arrangements. The MIF was determined using five muscular arrangements and three stiffness matrices, and the evaluation function was also calculated. Table 4 lists the average values of the five evaluation functions. The evaluation function of the joint stiffness vector was larger than that of the end-point. This is because the joint stiffness vector was calculated using the end-point and muscular arrangement. The evaluation function of the joint stiffness vector became larger because it was difficult to estimate the accurate joint stiffness vector when modeling errors occurred. In particular, the evaluation function of the joint stiffness vector using K e2 (x (θ d )) was larger than all the other values. This function synergistically became larger because the evaluation function of the end-point was larger than the other values. We propose a robust position control method for the arrangement error at this point (Matsutani et al., 2019a). A decrease in the error of the end-point position as a result of combining the proposed method with the discussed control method is expected. It was assumed that the decrease in the error of the joint stiffness vector was caused by the impact of the end-point position.

Evaluation of Arbitrary Vector β
We then investigated the influence of the arbitrary vector β on the control result. In this evaluation, the MIF was determined using the end-point stiffnesses K e1 (x (θ d )) listed in Table 2 and the arbitrary vectors β derived from Eq. 33 ∼ (Eq. 35). The simulation results obtained by applying the resulting MIFs were then compared.
Table 5 lists the arbitrary vectors k e derived from application of Method 1, using the outputs of Eq. 33 ∼ (35); the simulation results are shown in Figure 6. Figure 6A shows the transient responses of the end-point position; Figure 6B and Table 5 compare the control inputs; Figure 6C and Figure 6D show the end-point stiffness ellipsoids and loci of end-point positions in the task space, respectively. The results obtained using the arbitrary vector β 1 were identical to those shown in Figure 3D and were therefore, not repeated. It can be observed from Figure 3D and 6 that, while the stiffness ellipsoids at the desired position were the same, the transient responses of the end-point position differed. While analyzing the changes in end-point stiffness during movement, it was observed that the ellipsoids varied significantly near the initial position, but converged completely near the desired position. This occurs because the arbitrary vector β affected only the second term of Equation 26, a term that has no effect on the endpoint stiffness matrix K e1 (x (θ d )). This is so, because it becomes a vector that belongs to the null space of the matrix when the end-point position is the desired position. This indicates that Method 1 can be applied to change the transient responses of the end-point position, to achieve the desired end-point stiffness at the same position.
Changing the arbitrary vector β from its value in Figure 6B confirms that the arbitrary vector has no effect on the MIFs v d3 ∼ v d6 . Similarly, changing the end-point stiffness K e (x (θ d )) from its value inFigure 3B confirms that this factor has no significant effect on the MIFs v d1 , v d2 . These results suggest that muscles 1, 2 contribute to the maintenance of the transient response of the end-point position while muscles 3 ∼ 6 contribute to the maintenance of end-point stiffness. However, determination of the specific roles of the muscles will require further investigation and is left as a subject for future analysis.

METHOD OF DETERMINING MUSCULAR INTERNAL FORCE USING STIFFNESS AT THE INITIAL AND DESIRED POSITIONS
In this section, a novel method for uniquely determining the MIF based on the stiffnesses at the initial and desired positions (Method 2) has been introduced. By applying Method 2, the stiffnesses at the respective positions can be controlled while performing repetitive movements. An example of this is the pick-and-place operation, in which the stiffness is set to low when the robot picks up an object, and set to high when the robot places the object.    164.9, 135.2, 83.0, 18.4, 55.6, 20.9] T β 2 [5.5, 4.5, 110.6, −12.0, 60.9, 0.7] T β 3 [549.8, 450.6, 16.3, 91.5, 42.8, 69.5 154.6, 145.5, 69.0, 33.6, 36.1, 47.3] T β 2 [9.5, 0.5, 69.0, 33.5, 36.1, 47.3 The arbitrary vector k e can be obtained as the inverse of Eq. 38 as follows: The MIFs v (θ 0 ) and v (θ d ) calculated using this arbitrary vector correspond to the musculoskeletal system stiffening to values of K j (θ 0 ) and K j (θ d ) at positions θ 0 and θ d , respectively. The application of Method 2 using the MIF controller is shown as a block diagram in Figure 7.

Simulation Results
Method 2 was used to apply individual control inputs to the musculoskeletal system as it carried out repetitive movements. By changing the values of the arbitrary vector k e , different target endpoint stiffnesses were applied at various initial positions A, to serve as control inputs to the muscle for moving the system from A to a desired position B. Arbitrary values of k e were also applied as desired end-point stiffnesses at B, using the method. The simulated positioning control of the musculoskeletal system shown in Figure 1 was conducted to evaluate the abilities of Method 2. A simulation was performed to infer that the MIF determined using Method 2, can realize the stiffness at two positions. Table 6 lists the arbitrary vector determined using Method 2 and the set of initial and desired positions, A and B, used in the simulations, respectively.
In each simulation, the musculoskeletal system underwent controlled movement from the initial position A to the desired position B and then back from B to A. The results validated the effectiveness of Method 2 in terms of controlling the musculoskeletal system based on control of stiffness at the end-point and other positions. The physical and muscular parameters of the musculoskeletal system used in this simulation are listed in Table 1.

Position Control for Moving From Position A to B
Further, the position control process for motion between positions A and B was then simulated. The results of the simulation are shown in Figure 8. Figure 8A shows the endpoint stiffness ellipsoids and loci of the end-point position in the task space; Figure 8B shows the transient responses of the endpoint position; Figure 8C compares the respective control inputs. It can be observed from Figure 8 that the end-point position converged at the desired position x B , and that end-point stiffness K eB was achieved at the desired position x B . As shown in Equation 36, a balanced MIF of v (θ A ) at the initial position was required to achieve an end-point stiffness of K eA at this position. The end-point stiffness at the initial position x A in Figure 8A differed from the end-point stiffness K eA because the position control applied a balanced MIF v (θ B ) at the desired position x B . The control input v dl determined by Method 2 ( Figure 8C) satisfied Eq. 12.

Position Control for Moving From Position B to A
We then simulated position control for the movement from position B to A. The results of this simulation are shown in Figure 9. Figure 9A shows the end-point stiffness ellipsoids and loci of the end-point positions in the task space. Figure 9B shows the transient responses of the end-point position. While, Figure 9C compares the control inputs. It can be observed from Figure 9 that the end-point position converged at the desired position x A and that an end-point stiffness K eA was achieved at the desired position. The control input determined by applying Method 2 ( Figure 9C) satisfied Eq. 12. The results in this section and the preceding one demonstrate that Method 2 can be applied to uniquely determine the MIF. It is difficult to significantly alter the direction of the major axis of the end-point stiffness ellipsoid using the control input because the end-point  stiffness depends on the posture of the musculoskeletal system. However, this method can still be used to control the end-point stiffnesses at the two positions effectively. In this study, the control input can be uniquely determined using Method 2 because the control target is a two-link system with six muscles. If Method 2 were used in a system with a greater number of redundancies, it would not be possible to determine the control input uniquely. Meanwhile, it was considered to increase the selectable range of the direction of the end-point stiffness ellipsoid. In addition, the control input can be uniquely determined by setting subtasks according to the redundancies in this case.

CONCLUSION
In this study, two new methods for determining MIF using joint stiffness were developed, to solve the issues in a previously developed musculoskeletal system MIF feedforward controller. Using the first method (Method 1), a MIF that eliminates the  Frontiers in Robotics and AI | www.frontiersin.org September 2021 | Volume 8 | Article 699792 correlation between the response and stiffness at a desired position was determined. Using the second method (Method 2), an arbitrary vector that reflects the stiffnesses at the initial and desired positions was determined, and applied to uniquely determine the balance of MIFs at each position. Although it was found to be difficult to significantly alter the stiffness using the control input, as the end-point stiffness depends on the posture of the musculoskeletal system, we believe that the proposed system can be adapted to expand the region of feasible stiffness control ; and this will be a subject of our future analysis. The proposed methods can realize sensorless position and stiffness control. In addition, by combining the existing control method with the proposed method, runaway of movement at sensor crashes can be prevented. Further, we plan to experimentally verify this in a future study. Furthermore, the MIF is expected to be combined with feedback control in a future work. This paper describes the evaluation of the musculoskeletal system, including the modeling error. Feedforward control is incapable of handling unanticipated contact. A combination of the MIF and feedback control complement each other to improve accuracy. For example, a stable region of feedback control that tends to become unstable can be assured by combining feedforward control based on the potential no needing an inverse dynamics model (Matsutani et al., 2018). In future works, this combination needs to be discussed in detail.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.