# Combined Feedback Feedforward Control of a 3-Link Musculoskeletal System Based on the Iterative Training Method

Authors: Amin Valizadeh; Ali Akbar Akbari
Journal: BioMed Research International
 (2021)

---

## Abstract

The investigation and study of the limbs, especially the human arm, have inspired a wide range of humanoid robots, such as movement and muscle redundancy, as a human motor system. One of the main issues related to musculoskeletal systems is the joint redundancy that causes no unique answer for each angle in return for an arm’s end effector’s arbitrary trajectory. As a result, there are many architectures like the torques applied to the joints. In this study, an iterative learning controller was applied to control the 3-link musculoskeletal system’s motion with 6 muscles. In this controller, the robot’s task space was assumed as the feedforward of the controller and muscle space as the controller feedback. In both task and muscle spaces, some noises cause the system to be unstable, so a forgetting factor was used to a convergence task space output in the neighborhood of the desired trajectories. The results show that the controller performance has improved gradually by iterating the learning steps, and the error rate has decreased so that the trajectory passed by the end effector has practically matched the desired trajectory after 1000 iterations.


---

## Body

## 1. Introduction

The reaching movement is accounted for a huge part of hand movements. In all these activities, a swift and complex process occurs in the brain, and after processing, the generated control signals are transmitted to body motors, namely, muscles. This complex process in the brain comprises some levels. First, the desired trajectory is determined for reaching an object, and in the second step, the coordinates of the specified trajectory estimated by vision are converted into the body coordinates; in the last step, control commands are sent to the muscular system to go through the desired trajectory. Investigating the body’s musculoskeletal system’s control mechanism can lead us to develop a robust control technique that can be applied to rehabilitation robotics. The design process and application of the actuators in such robots are similar to the simulation of the human body’s neural control system. Many controllers have been introduced and employed to control such systems and produce motions similar to the human movement, which are of different design methods and performances based on their design space (robot task space, joint space, and muscle space) [1, 2]. Each of these spaces has its features and complexity, and as we move from the task space toward the muscle space, it will be difficult to design the controller because of the increasing space order. It should be noted that the controller design in the muscle space should be carried out carefully so that the forces out of the body are ignored, and the model approaches reality. The joint redundancy causes each angle to have no unique solution in return for an arbitrary trajectory of the arm’s end effector [3]. Another problem caused by redundancy is the lack of a unique solution for model forces to generate unique torque [4]. Many optimization techniques have been proposed to overcome this problem in classification [5], biology and robotics [6–9].On the other hand, in everyday life, we can easily carry out the most complex movements with the highest possible accuracy in the presence of the same redundancies with the least possible. One solution to overcome this complexity is that the central nervous system activates a small group of muscles—called muscle synergy [10]—which allows the control of body movements with less computational cost by reducing the required independent degrees of freedom. In recent years, Suetani and Morimoto [11] presented an innovative hypothesis under the title of “Virtual Spring-Damper Hypothesis,” where there is no need to optimize the redundancy criterion to overcome the redundancy problem. However, the previous problem is solved by applying this hypothesis, but due to the redundancy of muscles, we will deal with other problems that require the application of the muscle nonlinear models. In 2013, aiming to find the synergies of reaching and balancing movements on the musculoskeletal system of the hand, Tahara et al. [12] conducted a research study to investigate muscular integrity force data and the data of body receptors (proprioception and vision). The canonical correlation analysis (CCA) method, which follows the natural behavior of the body, was used to obtain the relation between the data related to muscles and body receptors. In a real system, the time delay and noise should be considered in the body’s actuator and sensor systems. Hence, this study is aimed at examining the effects of time delay and noise on determining the synergies of the hand’s musculoskeletal system. The results of this study not only can be used to understand the biological data of the motor control system but also can be applied as an artificial controller for a high-DOF robot. In this study, the motion equations for the 3-link musculoskeletal system of the human arm and the iterative learning controller are presented in Section 2. The results obtained from the simulation of the ILC (iterative learning control) with the neuro-fuzzy controller are discussed in Section 4. In Section 6, the remarks concluded from this study are described.


## 2. DOF Human Musculoskeletal Arm Model

The 3-DOF human musculoskeletal arm model used in this study consisting of three solid links and six monoarticular muscles is shown in Figure1. Since this arm moves on the horizontal plane, the effect of the gravity force can be ignored. As shown in Figure 1, this model consists of six muscles that can only apply tensile forces so that each joint moves by some of these related muscles. Muscles are assumed to be without weight and designed based on the Hill model, which are directly connected to links as
(1)fm=Pα¯−PAα¯C+C0l̇,α¯=α¯1,α¯2,⋯,α¯6ϵR6,P=diagp1,p2,⋯,p6ϵR6×6,Aα¯=diagα¯1,α¯2,⋯,α¯6ϵR6×6,C=diagc1,c2,⋯,c6ϵR6×6,C0=diagc01,c02,⋯,c06ϵR6×6.Figure 1
Schematic view of the 3-DOF musculoskeletal model for the hand.From Ref. [13], parameter fm is the muscles’ contractile force, which is the nonlinear function of the muscle’s contractile velocity and the control input produced in the central nervous system.ConsideringL1, L2, and L3 to be the first, second, and third links, respectively, as well as their relation angle with respect to the x-axis, first link, and second link as θ1, θ2, and θ3, respectively, the arm’s end effector position to the joint angles is expressed by the following equation:
(2)X=Jθ=L1cosθ1+L2cosθ1+θ2+L3cosθ1+θ2+θ3L1sinθ1+L2sinθ1+θ2+L3sinθ1+θ2+θ3̇∈R2.Defining length vectors of the muscles asl=l1l2l3l4l5l6T results in
(3)l=Wθ=r1π−θ1−acosr1s1+s12−r12r2π2+θ1−acosr2s2+s22−r22r3π−θ2−acosr3s3+s32−r32r4π2+θ2−acosr4s4+s42−r42r5π2+θ3−acosr5s5+s52−r52r6π2+θ3−acosr6s6+s62−r62∈R6,wherer1−6 and s1−6 represent the torque levers, as shown in Figure 1. The following equations are obtained by taking the time derivative of equations (2) and (3) with respect to time:
(4)Ẋ=Jθ̇,(5)l̇=WTθ̇.Ẋ=ẋẏT is the end effector velocity of the arm, θ̇=θ̇θ̇2θ̇3T is the angular velocity of joints, and l̇=l̇1l̇2l̇3l̇4l̇5l̇6T represents the stretch rate of muscles. Also, J∈R2×3 is the Jacobian matrix that shows the relation between the linear velocities of the arm’s end effector and angular velocities while WT∈R6×3 is the Jacobian matrix, which relates the contractile rate of muscles to the angular velocity of the joints as
(6)J=J11J12J13J21J22J23,J11=−L1sinθ1−L2sinθ1+θ2−L3sinθ1+θ2+θ3,J12=−L2sinθ1+θ2−L3sinθ1+θ2+θ3,J13=−L2sinθ1+θ2−L3sinθ1+θ2+θ3,J21=L1cosθ1+L2cosθ1+θ2+L3cosθ1+θ2+θ3,J22=L2cosθ1+θ2+L3cosθ1+θ2+θ3,J23=L3cosθ1+θ2+θ3,W=−r1r2000000−r3r4000000−r5r6.By assumingJ as a full-rank matrix, the inverse of equations (2) and (4) is obtained as follows:
(7)θ=Gx−1x∈R2,(8)θ̇=J−1ẋ∈R2.Gx−1x represents a vector with nonlinear functions, which shows the inverse kinematics from the task space to joint space. Also, J−1 shows the inverse kinematics from the task space velocity to the joint’s angular velocity. By substituting equation (8) into equations (3) and (5), we can state that
(9)l=GιGx−1x∈R6,(10)l̇=WTJ−1ẋ∈R6.Equation (9) demonstrates the inverse kinematics from the task space to the muscle space, which is applied to the controller’s feedforward behavior.By applying the principle of virtual work, the work done by muscle torque is defined as follows:(11)T=Wfm∈R2,where fm=f1f2f3f4f5f6T is the vector representing the tensile forces of muscles and T=T1T2T3T is the joint torque vector.By assuming thatW∈R3×6 is a row full-rank matrix during movement, the inverse of equation (11) is expressed as follows:
(12)fm=W+T∈R6,(13)W+=WTWWT−1∈R6×3.Besides, the static relation betweenT and the output vector of forces applied to the arm’s endpoint in the space F∈R2 is expressed as follows:
(14)T=JTF∈R3.By substituting equation (14) into equation (13), it is concluded that
(15)fm=W+JTF∈.Equation (15) demonstrates the static inverse relation between fm and F.


## 3. Iterative Learning Control

An ILC strategy of the PI type has been introduced in Reference [14] to trace an arbitrary time-dependent trajectory using the robotic arm model. The errors related to the position and velocity in a test are stored to be tuned for the next test by an input correction strategy. The data stored in the first step are multiplied by a factor and added to the input in the next test. Implementing a simple task space feedback control for a 2-DOF arm is considered by Tahara et al. to address the muscle space redundancy problem on the contractile output force [15]. They also studied multiple space variables to enhance the robustness of the 2-DOF arm exposed to sensory noises. Despite the nonlinear equations of the human arm’s motion, the suggested method sufficiently improves the system’s robustness regarding the traditional ILC methods [16]. Therefore, the proposed method is considered in our study. As discussed in the previous section, to compensate for the iterative learning controller’s input, there are three representatives of the state space, namely, muscle space, joint space, and task space. Therefore, any space that can better compensate for the control input is of great importance in achieving the desired performance. Furthermore, it should be noted that many noises cause damage to sensory information, and its huge impact on the movement of the musculoskeletal system is inevitable. Therefore, the system’s robustness to deal with the noise varies depending on the space in which the system is modeled. A new control strategy based on iterative learning, which uses the sum of state-space variables, is employed to improve the robustness of the system against noise. In the present paper, a case study is performed by considering the task space and muscle space as the spaces for feedback and feedforward behaviors, respectively. The control input of the muscles in the ith test is defined as follows:
(16)ui=−Wi+JiTKpΔxi−KυΔẋi+υi,where index i represents the test number, Kp=diagkp1,kp2∈R2×2>0 and Kυ=diagkυ1,kυ2∈R2×2>0 are the feedback coefficients of position and velocity in the task space, and υi is the feedforward parameter obtained from the iterative learning process. The error of position and velocity is defined as Δxi=xi−xd∈R2 and Δẋi=ẋi−ẋd∈R2, respectively; xd and ẋd also represent the end effector’s position and velocity, respectively. The feedforward parameter, υi∈R6, is not designed in the task space similar to feedback behavior, but it is modeled in the muscle space and updated as follows:
(17)υi=0,i=1,1−βυi−1−ΦΔιi−1+ΨΔι̇i−1,i>1,where Φ=diagϕ1,ϕ2,⋯,ϕ6∈R6×6>0 and Ψ=diagψ1,ψ2,⋯,ψ6∈R6×6>0 are the coefficient matrices of position and velocity, respectively; besides, the position error is defined as Διi=ιi−ιd and the velocity errors in the muscle space are expressed as Δι̇i=ι̇i−ι̇d. ιd∈R6 and ι̇d∈R6 are the length of muscles and their contraction rate relative to the position and velocity of the end effector, respectively. These parameters are obtained by calculating the inverse dynamic as
(18)Διi=GιGx−1xi−GιGx−1xd,Δι̇i=WiTJi−1Δẋi.In this study, the Gaussian noise is used as a noise which is applied to sensory information. An error in the initial conditions of two consecutive tests and dynamic oscillations due to different types of noises causes the general system to be unstable using the iterative learning controller. Therefore, to overcome these noises, Suetani and Morimoto [11] introduced a forgetting factor to update the iterative learning controller. Using this forgetting factor ensures that the final converged trajectory after good learning is in the desired trajectory neighborhood. In equation (17), β is the forgetting factor that has to satisfy the condition of 0<β<1. It is assumed that the muscle’s length and end effector position and velocity signals include Gaussian noise individually. Due to Refs. [15, 17], the magnitude of the noise existing in the end effector’s position and velocity is 4% of real data, and the magnitude of noise existing in the length of the muscle and its contraction rate is 50% of the real data. This is because the data related to the end effector’s position and velocity are obtained through observation, which is relatively accurate. However, the data related to the muscle’s length and contraction rate are received through the muscular bulk which has large electrical noise leading to inaccurate results [18].


## 4. Results

The simulation results are presented in this section. Tables1 and 2 demonstrate the numerical values associated with the 3-link model and the values related to the muscles’ physical properties, respectively. Also, the coefficients of the controller are listed in Table 3.Table 1
Numerical parameters of the model.
Length (m)Mass (kg)Inertial moment (kg·m2)CoM position (m)1st link0.311.930.01410.1652nd link0.271.320.01200.1353nd link0.150.350.00100.075Table 2
Geometric parameters of the muscles.
MuscleValue (m)l1r1=0.055s1=0.080l2r2=0.055s2=0.080l3r3=0.030s3=0.120l4r4=0.030s4=0.120l5r5=0.035s5=0.220l6r6=0.040s6=0.250Table 3
The parameters associated with the controller.
ParameterValueFeedback gainKp=8080Feedback gainKv=5050Learning gainΦ1=⋯=Φ6=250Learning gainΨ1=⋯=Ψ6=140Forgetting factorβ=0.3The controller is aimed at tracing a semicircular trajectory. Therefore, we consider the following trajectory:(19)x=0.2+0.1cost,y=0.55+0.1sint.The simulation’s total time is assumed to beT=πs, and the hand is initially located at point −0.1,0.55. Therefore, during the aforementioned period, the robot is expected to cover the semi-semicircular trajectory fully. For evaluating the robustness of the presented model against uncertainties, the simulation parameters have been changed by 5%. To compare the controller’s performance with similar counterparts, the model control results are compared to the neuro-fuzzy control method presented in our other paper [19]. In the cited article, the similar given trajectory was precisely followed by the muscle optimization, and the results concluded appropriate compliances with the hand’s natural motion. The model simulation was performed in MATLAB version 2021a running on an Intel Core i7 (2.8 GHz and 16 Gb RAM). For all simulations, the variable-step MATLAB ODEs solver ode45 with relative solver tolerance 1×10−4 was implemented, which took 4.817 s for the ILC controller compared to 10.045 s for the neuro-fuzzy controller.Figure2 depicts different trajectories that the controller has taken over 1000 iterations to reach the desired trajectory. As can be seen, as the controller’s performance is improved, the error is reduced gradually. Therefore, in iteration no. 1000, the trajectory is adjusted to the desired trajectory. Such a process is similar to learning and muscle memory that can perfectly go through a trajectory with practice and repetition. The trajectories that both controllers have gone through at the same time are shown in Figure 3. The results show that the proposed controller has better performance. In other words, if we exceed 1000 iterations in training the controller, we will observe a further improvement in the controller results. However, it should be noted that the simulation time increases with increasing the number of iterations.Figure 2
The trajectories passed by the model per 1000 iterations.Figure 3
Comparison of trajectories covered by the neuro-fuzzy adaptive controller and ILC.Figure4 displays the displacement of different joints during the movement scenario. The displacement of joints is similar to another. The adaptive controller performance is based on the optimization of the cost function and the iterative controller performance on learning; hence, Figure 4 designates that the proposed controller performance is acceptable compared to the adaptive controller’s performance using the neuro-fuzzy adaptive controller.Figure 4
Comparison of joint displacement between the neuro-fuzzy adaptive controller and ILC.Finally, the magnitude of forces applied to each muscle during the desired trajectory is illustrated in Figure5. The neuro-fuzzy controller uses muscle optimization; therefore, its force diagrams are much more ideal. On the other hand, the ILC controller has a more smooth bell-shaped profile similar to the agonist-antagonist paired muscles involved in the natural movements of the human body. However, both controllers have similar patterns. In addition, muscle forces in both methods are in the adequate range for the human body, where the ILC method has almost fewer force values.Figure 5
Comparison of the generated forces using the neuro-fuzzy adaptive controller and ILC.


## 5. Discussion

The controller has two pieces, the first of which is feedback input comprising task space variables, while the other part is the feedforward input, which is made up of muscle space parameters gained through the iterative learning algorithm. Although all controller gains are tuned in all iterations, the simulation results demonstrated that the hand endpoint in the first iteration is significantly different from the 1000th iteration. Despite the nonlinear equations of muscles, the controller could pass through the desired trajectory after the 1000th iteration. In addition, the path tracking error has considerably been mitigated by increasing the repetition number. The use of different variable spaces in conjunction with the learning algorithm was the primary reason for the sufficient accuracy of path tracking in the proposed controller.Furthermore, the travel time of the simulation was reduced to half using the ILC controller, compared to the neuro-fuzzy one following the 1000th repetition of the desired path. Force values for the given path were also in the sufficient force ranges of the human hand muscles. In most muscles, the endpoint passed through the given trajectory with a much lower force than the neuro-fuzzy controller. These results characterize the efficiency of this controller for musculoskeletal modeling in the human body. As a future study, we intend to conduct movement trials in actual and uncontrolled conditions using EMG signals plus effective technologies such as user-friendly contactless path recognition to increase the method’s productivity [5].


## 6. Conclusion

The controller’s performance was improved by iterating learning, and subsequently, the related error was reduced so that the final trajectory that has gone through simulation is practically adjusted to the desired trajectory. Such a process is similar to learning and muscle memory that can lead to perfectly going through a trajectory with practice and repetition. The quantitative comparison between the iterative learning controller and neuro-fuzzy controller results suggested that the proposed controller has a better performance. In other words, if we exceed 1000 iterations in training the controller, we will observe a further improvement in the controller results. However, it should be noted that the time required for solving the problem increases by increasing the number of iterations. By comparing the forces generated in the muscles for both controllers, it was observed that the maximum value of these forces for the current controller was less than that of the adaptive controller, although the average of generated force is higher for the current controller. Considering that the muscle forces’ optimization is one of the design indicators in adaptive controllers, it was not considered in the proposed method. Here, it was important that the controller can successfully guide the model on the desired trajectory in the presence of system uncertainties, and the forces applied to the muscles are in the desired range.



---

*Source: 8701869-2021-11-08.xml*