Simulation of a Lower Extremity Assistive Device for Resistance Training in a Microgravity Environment

: Lower extremity assistive devices (LEADs) have been developed in various ﬁelds, such as rehabilitation, military, and industry, in the form of exoskeleton robots or treadmills, and most of them are aimed at supporting muscle strength. However, unlike the aforementioned ﬁelds, the objective of LEADs developed in the space ﬁeld is to provide resistance training to prevent muscle atrophy, which is a problem that arises in astronauts during long-duration space ﬂights. Because the purpose of a LEAD is di ﬀ erent from those of systems that are intended for use under Earth gravity (1 g) condition, other factors should be considered for the system design. In this study, the appropriate locations and types of actuators for reproducing the kinematics and muscle-related state variables observed in 1 g normal walking in a microgravity environment were proposed, and the corresponding control inputs obtained using a dynamic optimization simulation method. In detail, two actuation types were proposed, considering the characteristics of a microgravity environment in which both the magnitude of the gravitational acceleration and the ground reaction force were decreased. Moreover, by using the proposed actuating system, the control inputs required to track kinematics data and muscle activity were obtained. A human lower-limb model, with six degrees of freedom, i.e., an 18-muscle model with the pelvis ﬁxed, was used with ideal actuators to apply torques or forces to joints or soles. Dynamic optimization was performed to solve these problems using direct collocation with OpenSim and MATLAB. Using the two proposed types of actuation, the results agreed with the kinematics and muscle activity of 1 g normal walking, and the total joint torques by the muscles also exhibited similar curves to that of the net joint torques under 1 g normal walking. The results of this study suggested an actuation method and its control input that can be used in the design of a LEAD for resistance training in microgravity.


Introduction
Lower extremity assistive devices(LEADs) have been developed in various forms such as a type in which an actuator is attached to each joint position [1][2][3][4] or a treadmill-based exoskeleton [5,6]. Such assistive devices, which have been developed for the purpose of muscle strength enhancement or power assistance, are intended for use under Earth gravity (1 g) conditions with a gravitational acceleration of approximately 9.81 m/s 2 (1 g) in the vertical direction. During the movement of a human body, this gravitational acceleration prevents muscle atrophy by continuously generating a resistance of the movement and makes it possible to maintain contact with the ground during walking, thereby producing a ground reaction force (GRF). However, the resistance force is barely generated when moving in a microgravity environment. The microgravity means 10 −6 times the gravitational acceleration on earth (1 g) and is used synonymously with zero gravity in general. A lack of resistance force in microgravity results in problems such as decrease in bone density and muscle atrophy during 2 of 19 long-duration space flights [7]. When walking in microgravity conditions, the kinematic pattern changes [8][9][10], and even if kinematically the same gait is performed as 1 g walking, a difference in net joint torques is caused by the lack of a GRF. To address these issues, several devices have been developed to assist muscle-strengthening exercises by generating resistance to movement [11,12]. Astronauts staying at the International Space Station (ISS) use these devices to perform more than two hours of exercise a day. The types of equipment for exercise include a cycle ergometer for pedaling, treadmill for walking or jogging, and a weight-bearing resistance device (RED) [11]. These devices are so designed that forces or torques are applied to the shoulders, feet, or joints of the human body to facilitate a user to perform actions such as squatting and walking under the resistance force. To develop an assistive device for resistance training in a microgravity environment, it is necessary to select the types and locations of appropriate actuators and determine the control inputs to implement the desired task.
However, some issues such as the microgravity environment, time required to produce a prototype, and measuring muscle activation make it difficult to determine the feasibility of the designed actuation method and obtain corresponding control inputs by experimental methods. These problems can be dealt with more efficiently by using simulations. Some related studies have been performed on the prediction of the muscular force generated using given kinematics data and external force conditions [13], and on optimal motion prediction to achieve a specific motor task [14][15][16][17]. In the design of a resistance training device under microgravity conditions, musculoskeletal simulations can be used to examine the actuator locations and acquire the corresponding control input, which enables the tracking of the desired motion data and muscle forces.
The simulation of the human musculoskeletal system is categorized into data tracking and predictive simulation based on its purpose. The main purpose of a data tracking simulation is to predict the muscle force using noninvasive methods based on experimentally obtained motion data and the GRF data. The goal of the predictive simulation is to obtain optimal movement to perform a specific motor task without experimental data. Both methods exploit optimization methods for different purposes. In data tracking, optimization is used to determine the muscle forces in the redundant nature of a human body, where multiple muscles can be involved in the actuation of a joint. In the case of predictive simulation, optimization is used to predict the behavior of minimizing or maximizing the objective function composed of specific motor tasks.
Optimization methods consist of dynamic optimization and static optimization. The dynamic optimization is a method that finds the optimal solution for the entire motion cycle by considering muscle contraction dynamics. Static optimization solves the optimization problem for each divided moment of the entire motion cycle without considering muscle contraction dynamics [18]. Data tracking simulation aims to predict muscle force, for which static optimization is used with inverse dynamics [19]. In addition, the computed muscle control (CMC) method that uses both static optimization and forward dynamics is used to obtain muscle excitation [20]. In predictive simulation (or optimal control), optimization and forward dynamics are used to obtain muscle excitations and corresponding actions to maximize or minimize the objective function, and integration is required over time [19]. As a method for solving the optimal control problem numerically, direct shooting, which uses the control variables as the design variables [21,22], and direct collocation, which optimizes the state and control variables [14,23,24], are widely used [25]. With a good initial guess, it is known that direct collocation is superior to the shooting method in terms of computational performance, which is approximately 100 to 250 times better when reaching the same solution [14]. Direct collocation can be used not only for predictive simulation but also for solving data tracking problems [13]. The human musculoskeletal system, which consists of muscles and bones, can be represented by a system that is driven by the governing equations of the muscle excitation-contraction dynamics, musculotendon dynamics, musculoskeletal dynamics, and skeletal dynamics [26]. It is possible to model the musculoskeletal system and perform simulations considering these governing equations by using software such as OpenSim [27], Anybody, and SIMM. OpenSim, which is free software, provides an interface with other software via the OpenSim application programming interface (API). It allows the human models provided by OpenSim to be accessed from external programs, which adopt other programming languages such as MATLAB and Python. This feature can be exploited, not only in data tracking simulation that relies on experimental data, but also for dynamic optimization, which enables predictive simulation that incorporates muscle contraction dynamics. Related simulation studies that feature a combined use of MATLAB and OpenSim include predictive simulation using a dynamic optimization method [14,28] and data tracking simulation of a 3-D walking motion [13].
Studies on the control inputs required when external actuator forces or torques are applied to a human body include CMC-based simulation, which yields the optimal locations of the external actuators and their assistive torques to minimize the metabolic cost of a human body. This is done while tracking the kinematics data of running under loading conditions [29]. A similar study was also performed using a unidirectional actuator under the load gait condition [30].
In the present study, the objective of the simulation is to obtain the control input of an external actuator to allow the musculoskeletal model to track the state variables generated during a 1 g normal gait in a microgravity environment. The reason for selecting walking motion as the target motion is that it is not only one of the motions included in the exercise routine during long-term space flight, but is also advantageous to verify the results as it is relatively easy to access experimental data compared to other movements. In order to achieve our research goals, we propose two actuating methods to compensate for the changes in the governing equation of the human musculoskeletal system due to the reduction in gravitational force. Second, for the proposed actuating methods, the control input required to track the state variables generated in a 1 g normal gait is obtained. We used OpenSim (version 3.3) and MATLAB with OPTI Toolbox, which provides the IPOPT solver for solving large-scale optimization problems [31]. The direct collocation method was used for obtaining the optimization solution.

Materials and Methods
The model of 'Gait10dof18musc' [32], which is provided with the OpenSim distribution, was modified and used for the simulation. The marker set data, GRF, and torque data during walking (about 1.2 m/s walking speed), which are also provided with the model, were used as the tracking objective data [33]. The data from one heel strike (HS) to the next HS, when considered at the right leg, were taken as one gait cycle. The original model has 10 degrees of freedom (pelvic translation in x and y directions (2), pelvic tilt (1), flexion/extension in bilateral hip, knee, and ankle (6), and flexion/extension of lumbar (1)) without foot-ground contact model. To simulate resistance training in microgravity, the kinematic degrees of freedom of the model were modified to a total of six degrees of freedom (flexion/extension in bilateral hip, knee, and ankle) with the pelvis fixed to the ground. The upper body located above the pelvis was also removed so that only the lower body remained. The gravitational acceleration was set to zero. For other damping effects on the joints of the human body, we added dampers by referring to the model of the relevant study [13]. Actuators representing the assistive device were added to the model using the "coordinate actuator", "point actuator" and "torque actuator" functions. The functions apply forces or torques ideally to the joints of the human body model. It enables to exclude undesirable effects that occur in actual devices such as the power loss caused by joint axis misalignment between human and device.
The purpose of the simulation is to allow the human model to track the state variables (generalized coordinate, generalized speed, muscle activation, and muscle fiber length) measured and predicted during 1 g normal walking with the assistance of an external actuator under microgravity conditions. Therefore, it is necessary to obtain the state variables from the marker set and GRF data in order to generate target tracking data. To do this, we performed inverse kinematics (IK), residual reduction algorithm (RRA), and CMC consecutively [27], and obtained the state variables ( Figure 1a). The CMC uses static optimization to determine human muscle force at the given kinematic data. In slow motions such as walking, static optimization yields nearly the same results as using dynamic optimization, and computation is much faster [34]. Therefore, we used CMC to obtain target tracking data efficiently. dynamically match the modified human body model for simulation in the microgravity situation. Therefore, using this data directly as an initial guess can adversely affect the convergence. A method that uses a solution of the system of defect functions, which were set as the constraints of the optimization problem, as an initial guess, can be exploited to improve convergence [13]. For this, we solved the defect functions by using the MATLAB 'fsolve' function that solves a system of non-linear equations and used the obtained solution as an initial guess for the optimization problem. The OpenSim API was used to link the human musculoskeletal model and corresponding governing equations (muscle contraction dynamics, musculoskeletal dynamics, and musculotendon dynamics [35]) provided by OpenSim with MATLAB, which has advantages in optimization and matrix calculation. The optimization design variables, states and controls, in the MATLAB optimization solver, can be made to satisfy the governing equations of the OpenSim model by exploiting the defect constraints ( Figure 1b). In the process of solving the optimization problem using the direct collocation method, all the state and control variables are discretized by mesh interval, resulting in a large number of optimization design variables. To solve this large-scale problem, we used the IPOPT solver, which is more efficient than the MATLAB optimization solver, fmincon [28].
The gait period was divided into several nodes with equal intervals. Starting with a mesh of 25 nodes, 51, 75, and 101 meshes were used for the simulation. The grid refinement method, which uses the state and control variables obtained from the simulation using the previous node density, was used as an initial guess, except for the mesh of 25 nodes [28]. The objective function was set as the sum of the integral of the squared error: Initial guess is also required to solve the optimization problem. Because the state and control variables described in (Figure 1a) are datasets obtained from 1 g normal walking data, they do not dynamically match the modified human body model for simulation in the microgravity situation. Therefore, using this data directly as an initial guess can adversely affect the convergence. A method that uses a solution of the system of defect functions, which were set as the constraints of the optimization problem, as an initial guess, can be exploited to improve convergence [13]. For this, we solved the defect functions by using the MATLAB 'fsolve' function that solves a system of non-linear equations and used the obtained solution as an initial guess for the optimization problem.
The OpenSim API was used to link the human musculoskeletal model and corresponding governing equations (muscle contraction dynamics, musculoskeletal dynamics, and musculotendon dynamics [35]) provided by OpenSim with MATLAB, which has advantages in optimization and matrix calculation. The optimization design variables, states and controls, in the MATLAB optimization solver, can be made to satisfy the governing equations of the OpenSim model by exploiting the defect constraints ( Figure 1b). In the process of solving the optimization problem using the direct collocation method, all the state and control variables are discretized by mesh interval, resulting in a large number of optimization design variables. To solve this large-scale problem, we used the IPOPT solver, which is more efficient than the MATLAB optimization solver, fmincon [28].
The gait period was divided into several nodes with equal intervals. Starting with a mesh of 25 nodes, 51, 75, and 101 meshes were used for the simulation. The grid refinement method, which uses the state and control variables obtained from the simulation using the previous node density, was used as an initial guess, except for the mesh of 25 nodes [28]. The objective function was set as the sum of the integral of the squared error: where e i (t) is the error between the target variables X i,re f (t) and design variables X i (t) for the ith state variable. In the objective function, the target state variables of 6 generalized coordinates, 6 generalized speeds, and 18 muscle activations were included. The design variables consisted of state variables: generalized coordinates q, generalized speeds . q, muscle activations a, fiber lengths l, and control variables: muscle excitations u, external actuator torques τ or forces f . The bound constraints of these optimization variables are constructed as expressed in (Equations (2)-(8)) so that the state variables of generalized coordinate (Equation (2)), generalized speed (Equation (3)), and muscle activation (Equation (4)) do not deviate beyond a certain range from the tracking target data value. Even if the assistive (or resistance) external torques or forces are applied to the body, it is difficult to accurately track the state variables produced in a 1 g environment due to the gravity and external force changes. Therefore, we set the range of the design variables as the upper and lower bounds. The bound values were set based on the original state variable values during one gait cycle and were empirically provided (α, β, γ expressed in Equations 2-4). The muscle activation was set between 0 and 1 (Equation (4)), and the fiber length was limited to the range of 0.25 to 1.25 times the optimal fiber length l opt as expressed in (Equation (5)) [13]. For the control variables (Equations (6)-(8)), the muscle excitation was set to be between 0 to 1 (Equation (6)), and for the external torque or force inputs (Equations (7) and (8)), the appropriate minimum (τ min and f min ) and maximum (τ max and f max ) values were set considering the location and type of the actuator.
To ensure that the optimization variable complied with the governing equations of the OpenSim human body model, the defect constraints were defined as (Equation (9)). In this part, the dynamics of the system was transformed into an algebraic equation between adjacent nodes. In the process of approximating the state variable differential by a finite difference, we used the midpoint Euler method, which uses the state variable value at the midpoint of two adjacent nodes [15]. In (Equation (9)), x, v, and t represents state, control variables, and time, respectively, and f is the time derivatives of the state variables. The subscript i represents the ith node and its maximum value is N − 1: The state variable x in (Equation (9)) means q, . q, a, l mentioned in (Equations (2)-(5)), and v denotes u, τ, f used in (Equations (6)-(8)).
Since human body walking is periodic, the actuator input should also be input periodically. To satisfy this, additional periodicity constraints were imposed so that the actuator input matches at the first and last node of the walk (Equations (10) and (11)): where the τ i,1 means ith torque actuator input at 1st node, and the τ i,N means input at Nth node. The f can be interpreted in the same manner. The relationship between the forces acting on the human body due to the muscle, GRF, and external actuator, and the motion caused thereby is expressed in (Equation (12)) [26,36]: where M q is the system mass matrix, C q, .
q is Coriolis and centripetal term, G q is gravity term, q is external force term, and q, q is the vector of generalized coordinates, velocities, and accelerations.
In the microgravity environment, the gravitational force is reduced, resulting in a decrease in body weight (BW). This causes a reduction in the GRF, which is an external force generated in the stance phase. Thus, when compared to a 1 g gait, kinematically different gait patterns and joint moments are generated. The GRF and BW correspond to the E q, . q and G q terms in (Equation (12)). If external actuators are used to compensate for the effects caused by the change in the E q, . q and G q terms in (Equation (12)), then the same governing equation as that in the case of 1 g normal walking can be established, and the tracking can be achieved. For the model modified for microgravity simulation in this study, (Equation (12)) can be simplified as (Equation (13)), and the goal of the simulation in the view of the musculoskeletal system can be summarized as finding L to satisfy the desired q, . q, and F MT .
where BL represents a vector of six externally applied loads (force/moment on each foot, or six exoskeleton torques), L is the magnitudes of the actuator inputs, and B is an invertible (6 × 6) matrix. Muscle activation dynamics (Equation (14)) and musculotendon dynamics (Equation (15)) were used to obtain F MT in (Equation (13)) [37,38]: . .
A simulation was performed for Type 1, in which the actuators are attached to the position similar to that when the GRF acts on the human body during walking ( Figure 2a). In the simulation for Type 2, the rotary actuator was attached to the generalized coordinates of the model (Figure 2b).
Appl. Sci. 2020, 10, x FOR PEER REVIEW  6 of 19 where the , means th torque actuator input at 1st node, and the , means input at th node. The can be interpreted in the same manner. The relationship between the forces acting on the human body due to the muscle, GRF, and external actuator, and the motion caused thereby is expressed in (Equation (12)) [26,36]: where ( ) is the system mass matrix, ( , ) is Coriolis and centripetal term, ( ) is gravity term, ( ) is moment arm matrix, is a vector of muscle forces, ( , ) is external force term, and , , is the vector of generalized coordinates, velocities, and accelerations.
In the microgravity environment, the gravitational force is reduced, resulting in a decrease in body weight (BW). This causes a reduction in the GRF, which is an external force generated in the stance phase. Thus, when compared to a 1 g gait, kinematically different gait patterns and joint moments are generated. The GRF and BW correspond to the ( , ) and ( ) terms in (Equation (12)). If external actuators are used to compensate for the effects caused by the change in the ( , ) and ( ) terms in (Equation (12)), then the same governing equation as that in the case of 1 g normal walking can be established, and the tracking can be achieved. For the model modified for microgravity simulation in this study, (Equation (12)) can be simplified as (Equation (13)), and the goal of the simulation in the view of the musculoskeletal system can be summarized as finding to satisfy the desired , , and .
where represents a vector of six externally applied loads (force/moment on each foot, or six exoskeleton torques), is the magnitudes of the actuator inputs, and is an invertible (6 × 6) matrix.
Muscle activation dynamics (Equation (14)) and musculotendon dynamics (Equation (15)) were used to obtain in (Equation (13)) [37,38]: = , , A simulation was performed for Type 1, in which the actuators are attached to the position similar to that when the GRF acts on the human body during walking ( Figure 2a). In the simulation for Type 2, the rotary actuator was attached to the generalized coordinates of the model (Figure 2b). Type 1: Actuating method in which the x and y -directional forces and torque actuator act on the sagittal plane at the midpoint of the right and left soles of the model with fixed pelvis. The conceptual underlying motivation of this system is that the difference from the 1 g environment can be compensated by an external torque corresponding to the GRF, which presents the largest change in the governing equation when the gravity condition changes. Even though the contact point between a foot and ground gradually moves from the heel to the toe in the stance phase during walking, most of the assistive devices apply forces to the entire sole. To reflect this, two forces and one torque (F x , F y , τ z ) actuators acting on the midpoint of the sole were considered (Figure 2a). The force vectors F x and F y are expressed in global coordinates.
Type 2: It is an actuation method in which the assistive torques are applied by one degree of freedom rotary actuators attached to the hip, knee, and ankle (HKA) joints. This type of assistance is widely adopted in lower-limb exoskeleton robots, which have been developed so far and is easy to use independently; therefore, it has excellent extensibility for other applications. A total of six actuators were constructed (Figure 2b, τ h , τ k , τ a for left and right).
OpenSim and MATLAB were used for the simulation, and IPOPT was used as the optimization solver. An i5-6500 3.2 GHz CPU and 16 GB for RAM were used.

Results
The state variables were closely tracked by the model, and the control inputs of the external actuators were obtained simultaneously. Simulations were performed on both legs, and the results on both sides were nearly identical, except for the timing. Therefore, for readability, only the results for the right leg representing HS to the next HS are shown in this section. For results on the opposite side, see Supplementary Materials.

Actuation Type 1
In the case of Type 1, the state tracking results closely tracked the target state variables regardless of the mesh density used in the simulation. The state variables that were included in the objective function, which are the generalized coordinates (Figure 3a [39,40] were added to ( Figure 4) for comparison with the target muscle activations obtained from CMC. The EMG patterns were scaled to match the CMC results. All of the muscle activations had similar patterns with EMG, except for the bifemsh and iliopsoas muscles. The fiber length and muscle excitation, which were not included in the objective function in (Equation (1)), showed practically identical results to the target values ( Figure 5).
The assistive torques and forces were obtained as shown in ( Figure 6). For comparison, GRF, which occurs during 1 g normal gait, are also shown (red dashed lines in Figure 6). Due to the periodicity constraints assigned to the actuators, the actuator inputs showed the same values at both ends ( Figure 6). F x showed peak force of 400 N around 55% and showed about 2 times magnitude difference from GRFx data. Except for that point, the entire curves are located in the range −100 N to 100 N (Figure 6a). The maximum force of F y was about 1 BW and 1.3 BW at 13% and 50% of the walking cycle, which is similar to the peak force of GRFy (Figure 6b). The slopes in the 0-10% and 50-60% of the gait cycle were in good agreement with the GRF. The magnitude of the GRFy was positive in all sections, whereas F y showed negative direction assistance at 90%. The assistive torque in the z-direction was significantly different from the ground reaction moment in the same direction ( Figure 6c). The smoother control input curves were obtained as the mesh density increased. We can see that when using the smallest number of 25 nodes, the results are not smooth ( Figure 6). Therefore, it is expected that the reliability of the results would be decreased significantly for node mesh smaller than 25.       . The difference between the mesh densities is small. The shaded areas are surface EMG patterns obtained from natural speed walking in healthy adults [39,40]. It is scaled to match the simulation data for comparison. See Supplementary Materials for results on the left limb ( Figure S2). obtained from natural speed walking in healthy adults [39,40]. It is scaled to match the simulation data for comparison. See Supplementary Materials for results on the left limb ( Figure S2).   Figure S4).

Actuation Type 2
In the case of Type 2, which provides torque assistance to all the lower limb joints, the target tracking results also showed a good match, and there was no significant difference according to mesh density. The generalized coordinate (Figure 7a  obtained from natural speed walking in healthy adults [39,40]. It is scaled to match the simulation data for comparison. See Supplementary Materials for results on the left limb ( Figure S2).   Figure S4).

Actuation Type 2
In the case of Type 2, which provides torque assistance to all the lower limb joints, the target tracking results also showed a good match, and there was no significant difference according to mesh density. The generalized coordinate (Figure 7a-c), generalized speed (Figure 7d-f), and muscle activation (Figure 8a-i) showed good agreement with tracking values. The experimental EMG patterns during a natural walking speed of adults [39,40] were added to ( Figure 8) for comparison with the target muscle activations obtained from CMC. The EMG patterns were scaled to match the CMC results. The fiber length (Figure 9a-i) and muscle excitation (Figure 9j-r) which were not Figure 6. Assistive forces and torques of the external actuators (Type 1). The actuator inputs are shown for the right limb ((a-c), solid lines). For comparison, the GRF data are also shown for the 1 g walking (red dashed lines). Due to the periodicity constraints assigned to the actuators, the actuator inputs showed the same values at both ends. The smoother control input curves were obtained as the mesh density increased. See Supplementary Materials for results on the left limb ( Figure S4).

Actuation Type 2
In the case of Type 2, which provides torque assistance to all the lower limb joints, the target tracking results also showed a good match, and there was no significant difference according to mesh density. The generalized coordinate (Figure 7a included in the objective function (Equation (1)) also showed practically the same results as the tracking values.   included in the objective function (Equation (1)) also showed practically the same results as the tracking values.   . The difference between the node densities is small. The shaded areas are surface EMG patterns obtained from natural speed walking in healthy adults from the literature [39,40]. It is scaled to match the simulation data for comparison. See Supplementary Materials for results on the left limb ( Figure S6). obtained from natural speed walking in healthy adults from the literature [39,40]. It is scaled to match the simulation data for comparison. See Supplementary Materials for results on the left limb ( Figure  S6).  Figure S7). Figure 10 shows the torque input of the Type 2 actuators. For comparison, the differences between the results of inverse dynamics in the two environments, i.e., in the absence of GRF data at 0 g (microgravity) and in the presence of GRF data at 1 g (earth), were also represented together. (ID difference in Figure 10). The knee and ankle joint actuators had similar results to 'ID difference' (Figure 10b,c). The hip actuator (red dashed line in Figure 10a) showed fluctuating results, and the deviation of maximum and minimum values was larger than that of 1 g GRF data. The maximum magnitude was about 50 Nm at 10% and 60% of the gait cycle, and the minimum was about -50 Nm at 85%.   Figure S7). Figure 10 shows the torque input of the Type 2 actuators. For comparison, the differences between the results of inverse dynamics in the two environments, i.e., in the absence of GRF data at 0 g (microgravity) and in the presence of GRF data at 1 g (earth), were also represented together. (ID difference in Figure 10). The knee and ankle joint actuators had similar results to 'ID difference' (Figure 10b,c). The hip actuator (red dashed line in Figure 10a) showed fluctuating results, and the deviation of maximum and minimum values was larger than that of 1 g GRF data. The maximum magnitude was about 50 Nm at 10% and 60% of the gait cycle, and the minimum was about -50 Nm at 85%. the simulation data for comparison. See Supplementary Materials for results on the left limb ( Figure  S6).  Figure S7). Figure 10 shows the torque input of the Type 2 actuators. For comparison, the differences between the results of inverse dynamics in the two environments, i.e., in the absence of GRF data at 0 g (microgravity) and in the presence of GRF data at 1 g (earth), were also represented together. (ID difference in Figure 10). The knee and ankle joint actuators had similar results to 'ID difference' (Figure 10b,c). The hip actuator (red dashed line in Figure 10a) showed fluctuating results, and the deviation of maximum and minimum values was larger than that of 1 g GRF data. The maximum magnitude was about 50 Nm at 10% and 60% of the gait cycle, and the minimum was about -50 Nm at 85%.   (Table 1). In both types of simulation, the design variables comprised of 48 state variables (6 + 6 + 18 + 18) and 24 control variables (18 + 6). As the mesh density increases, the value of the objective function becomes smaller. In both Type 1 and Type 2, the number of iterations tended to decrease with increasing mesh density. Type 1 required 1.9-4.0 h of CPU time and Type 2 required 0.9-3.0 h.

Discussion
We obtained the actuator inputs needed to generate the load of walking motion during resistance training in a microgravity environment and the tracking results of state variables using the direct collocation optimization method. The state variables showed good tracking results, and the net joint moments calculated from the obtained results were similar to those of inverse dynamics at 1 g normal walking. Fiber length and muscle excitation were not included in the objective function but results equivalent to the target data were obtained. In addition, periodic inputs were obtained due to the periodicity constraint included in the actuator input.
In the case of Type 1 assistance (Figure 2a), the actuation system was designed with the underlying intention to compensate for the external torque term, which is a decreasing term in the governing equation in reduced gravity condition. Therefore, the result was expected to be similar to that from GRF data for 1 g normal walking. As expected, the simulation results showed similar results to 1 g GRF data. The negative values observed in the F y input corresponding to the swing phase (gait cycle during 80-95% in Figure 6b) may be to generate a force by the weight of the leg at 1 g walking under a microgravity condition. The torque in the z-direction (red dashed line in Figure 6c) is the torque produced by the GRF on the feet and is calculated using COP x × GRF y . Type 1 has two linear actuators and one torque actuator acting on the center of the soles. However, during normal walking, the center of pressure (CoP) of a sole moves towards the toes from the heel, and the resulting effect is compensated by the external torque input at a fixed position in the center of the sole. In Figure 6c, the results of the simulations and the calculated torque results have been similarly obtained, so it is possible to confirm that the effect produced by the movement of the COP, the point at which the GRF is applied, has been compensated by the τ z actuator. Figure 10 shows the torque input of the external actuators required to follow the state variables in 1 g normal walking using Type 2. Because Type 2 has actuators at the hip, knee, and ankle joints, it can be presumed that the external torque input will appear as the difference between the net joint torques obtained by performing inverse dynamics at 0 g and 1 g. Therefore, for comparison, the "ID difference" (red dashed lines) is shown as the difference in the net joint torques obtained by the inverse dynamics using the same kinematics data for the 1 g normal gait (1 g with GRF) and microgravity (0 g without the GRF). The knee and ankle assist (Figure 10b,c) show similar results as expected, but the hip shows relatively different tendency (Figure 10a). The potential reason for this discrepancy may be due to differences between the original and modified model. The original model was originally capable of tilt, x and y motions in pelvis and lumbar flexion/extension. Thus, in inverse dynamics, there were forces or moments in the pelvis and lumbar acting in those directions. Resistance training in microgravity requires the body to be fixed, so the pelvis is fixed in the modified model to reflect it. As a result, the forces and moments mentioned above can no longer act on the corresponding body segment. However, the goal of the simulation was to track the kinematics data and muscle related state variables obtained from the original model at 1 g normal walking in the modified model. To achieve this, the forces and moments that existed in the original model must be compensated by the external actuators in the modified model, which is why this difference appears to have occurred. Figure 11 shows the net joint torques by the muscles calculated from the state variables, as obtained from the Type 1 simulation, and the joint torques acquired by inverse dynamics using the 1 g gait data.
The movement of the ankle joint (solid lines in Figure 11c) shows a tendency opposite to that of the Type 1 torque actuator (solid lines in Figure 6c). Therefore, it can be regarded that the torque actuator is mainly used for generating a resistance torque against the ankle torque.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 19 expected, but the hip shows relatively different tendency (Figure 10a). The potential reason for this discrepancy may be due to differences between the original and modified model. The original model was originally capable of tilt, and motions in pelvis and lumbar flexion/extension. Thus, in inverse dynamics, there were forces or moments in the pelvis and lumbar acting in those directions. Resistance training in microgravity requires the body to be fixed, so the pelvis is fixed in the modified model to reflect it. As a result, the forces and moments mentioned above can no longer act on the corresponding body segment. However, the goal of the simulation was to track the kinematics data and muscle related state variables obtained from the original model at 1 g normal walking in the modified model. To achieve this, the forces and moments that existed in the original model must be compensated by the external actuators in the modified model, which is why this difference appears to have occurred. Figure 11 shows the net joint torques by the muscles calculated from the state variables, as obtained from the Type 1 simulation, and the joint torques acquired by inverse dynamics using the 1 g gait data. The movement of the ankle joint (solid lines in Figure 11c) shows a tendency opposite to that of the Type 1 torque actuator (solid lines in Figure 6c). Therefore, it can be regarded that the torque actuator is mainly used for generating a resistance torque against the ankle torque.    expected, but the hip shows relatively different tendency (Figure 10a). The potential reason for this discrepancy may be due to differences between the original and modified model. The original model was originally capable of tilt, and motions in pelvis and lumbar flexion/extension. Thus, in inverse dynamics, there were forces or moments in the pelvis and lumbar acting in those directions. Resistance training in microgravity requires the body to be fixed, so the pelvis is fixed in the modified model to reflect it. As a result, the forces and moments mentioned above can no longer act on the corresponding body segment. However, the goal of the simulation was to track the kinematics data and muscle related state variables obtained from the original model at 1 g normal walking in the modified model. To achieve this, the forces and moments that existed in the original model must be compensated by the external actuators in the modified model, which is why this difference appears to have occurred. Figure 11 shows the net joint torques by the muscles calculated from the state variables, as obtained from the Type 1 simulation, and the joint torques acquired by inverse dynamics using the 1 g gait data. The movement of the ankle joint (solid lines in Figure 11c) shows a tendency opposite to that of the Type 1 torque actuator (solid lines in Figure 6c). Therefore, it can be regarded that the torque actuator is mainly used for generating a resistance torque against the ankle torque.   The direct collocation method divides state and control variables by the number of node meshes. In Table 1, it can be seen that the total number of variables ranges from 1800 to 7272 depending on the number of node meshes used. In addition, the number of iterations can be seen to decrease gradually. This may be because the grid refinement method allowed the simulation to start closer to the optimal solution as the node mesh increased.
The most computationally expensive process during simulation is the numerical computation of constraint Jacobian. Therefore, as the mesh density increases, more CPU time was required (CPU time (s), Table 1). The minimum objective function value means the sum of the integral of the squared error between the target tracking data and the optimal solution (Equation (1)). This value tends to get smaller as a large number of node meshes are used. For the 101 nodes in TYPE 2 (Table 1), this value increased slightly, which will be reduced by setting small convergence criteria.
In this study, the actuator position and corresponding control input, which can be used in the design of a LEAD for resistance training in microgravity, were obtained by an optimization method using OpenSim and MATLAB. We calculated all the state variables from the 1 g normal gait data and obtained the position and control input of the external actuator to track the state variables. We did not perform the muscle force predictions from the experimental data similar to the previously performed data tracking methods. Comparing the existing data tracking method [13] with the method used in this study, there is a similarity in that it tracks data obtained during movement. However, there is a difference in that the objective is not to conduct muscle force prediction based on the human walking strategy such as metabolic cost minimization. One of the main goals is to obtain the actuator input for tracking. Accordingly, it can be used to check whether the desired target data can be tracked with the actuator configuration presented in the design of LEAD, and the corresponding control input can be obtained so that the capacity of the actuator can be selected or utilized for other controls.
The Type 1 and Type 2 actuation method used in this study was chosen based on the actuation method of the strength assist system that has been developed so far. Systems such as Type 1 are called end-effector types used in the field of rehabilitation [41,42], and the form of Type 2 is actually similar to the X1 [11] that was developed by NASA for long-term flight resistance training for astronauts. The form of Type 1 has the advantage of easy setup since there are relatively few fasteners, and Type 2 is typically smaller in volume than Type 1 and suitable for measuring the user motion.
One of the most important factors in a robot that transmits force between the human body and the robot is user safety and comfort. In order to apply the simulation results obtained from this study to actual resistive training devices, appropriate actuators and control methods should be used, taking these factors into account. For this purpose, the Series Elastic Actuator (SEA) can be used. The SEA has a structure in which spring elements are attached to the transmission part of the actuator and calculates force or torque from spring displacement [43,44]. The SEA, a type of compliant actuator, has been used for torque control of wearable exoskeleton robot [43,45] because it improves safety due to its low impedance characteristics and also has excellent force control performance [46]. Hard-stops and motion limiting controllers are also required to prevent hyperextension or hyperflexion of the joints. In the control method, torque and force profiles should not be applied directly to the human joints as a periodic function over time. Because if the user does not produce torques of exactly the same magnitude with opposite direction, it can cause severe injury. Therefore, in practice, desired torque curves should be created based on the angle of the joint and then applied to the user by considering walking phase [46,47]. It will require further experimentation and verification for actual implementation.
Although both systems were able to achieve the simulation's objective of muscle activity and kinematics tracking (Figures 3, 4, 7 and 8), the advantages and disadvantages of each system are prominent from different perspectives. In addition to the muscle atrophy focused on this study, the Type 1 system is advantageous from the point of view of preventing bone loss. Because human bones are generally generated in the direction of compressive load, Type 1, which has a structure in which y-directional force acts on the sole, will cause more bone-loading in the longitudinal direction. However, Type 2 may have more advantages in performing partial joint movements or various motions.
We have obtained a solution to the problem by performing the simulation for the design of a LEAD in a microgravity environment; however, there are also some limitations of this study. First, the model only has degrees of freedom in the sagittal plane, even though a three-dimensional model would more accurately reflect the human body. Although the motion in the sagittal plane is more dominant than in the other directions, the influence of the muscles that simultaneously affect two or more degrees of freedom may not be correctly evaluated. Second, the ideal actuator is used without regard to the effects of the mass or inertia, power transmission loss, and axis misalignment between the human body and assistive device. Therefore, problems arising in the actual world such as power loss caused by the misalignment of the rotation axis between the actuator and human joint could not be considered. Third, the interaction between the human and device is not considered. If the human body is assisted by an external device, a different walking strategy can be adopted even under the same joint moment because it changes the walking pattern to adapt to the device or to use the assistive force better [15,17]. Experiments or additional predictive simulations should be performed to deal with these areas. Fourth, we use target tracking data from the CMC results. The use of the result data of CMC may not have a significant effect on the simulation results in that the previous study suggests that dynamic optimization and static optimization have the same results in low-speed gait motion [34]. However, it is better to use the state variables obtained by predictive simulation for rapid or novel movements.
In this study, the direct collocation method used to obtain the dynamic optimization solution transforms the optimization problem into nonlinear programming (NLP) problem via the parametrization of the state variables and control variables. This results in a sizable constraint Jacobian matrix, which requires a significant computation time; however, it can be shortened by using the sparsity of the matrix [25]. In addition, because the governing equations of the system are transformed into algebraic equations and integration is not necessary, the calculation speed is higher than that of the shooting method, which uses the method of discretizing only the control variables and integrating the governing equations. The direct collocation method can be used for both data tracking and predictive simulation, and studies using this method such as maximum height jumps [14] and 3D data tracking [13] have been performed.
The purpose of the existing data tracking problem study was to predict the muscle force from the given kinematics data and external force, GRF, and a cost function, such as the muscle activation squared or metabolic cost, which is widely used as the physiological term of the human gait, is used [13,29,30]. References [29,30] predicted the muscle force while tracking the gait data without an external assist device using CMC. In addition, when an external actuator was used for assistance, the muscle force change and corresponding control input were obtained to analyze the metabolic cost reduction and the relation between the actuator position for assistance.
CMC is a tool that can efficiently track a target dataset by using static optimization, forward dynamics, and proportional-derivative (PD) control [20]; however, in this study, it is difficult to use in the simulation. The reason is that the CMC tool calculates the muscle force to track it with kinematics data and an external force (GRF), whereas the simulation performed in this study requires tracking both kinematics data and muscle force. In order to obtain the Type 1 and Type 2 actuator inputs presented in this study using the method of CMC in OpenSim, the objective function of the static optimization needs to be changed, which is initially set to the square of the muscle activity. However, it cannot be easily changed in MATLAB. In this study, because the optimization was performed outside of OpenSim by performing dynamic optimization by the DC method, there was no restriction on the setting of the objective function; therefore, it could be constructed to minimize the tracking error of the state variables. Because muscle force prediction was not the primary objective, the physiological term was not used in the objective function as in previous studies.
In this study, we used the six degrees of freedom lower-limb model with six actuators. Thus, the input L could have been solved much faster algebraically by multiplying the inverse of B in (Equation (13)). Although we focused a lot on methodological challenges, we were able to obtain muscle-related variables with less mathematical effort required. The direct collocation approach we used in this paper can also be easily applied to predictive simulation by altering the objective function and constraints. In fact, one of the challenges of using MATLAB and OpenSim in conjunction with the large-scale optimization solver IPOPT is to calculate the numerical derivatives of the objective function and constraints if they could not be expressed analytically. We provide the code used in this simulation with this paper, which is expected to be useful for researchers to conduct related research.

Conclusions
Simulation for a LEAD in microgravity was performed using MATLAB and OpenSim. The net joint torques by the muscles calculated from the results of Type 1 and Type 2 (Figures 11 and 12) showed similar results to the net joint torques of a 1 g normal gait. From this, it can be concluded that tracking of the state variables corresponding to the 1 g walking level can be achieved by using the actuation method and control input proposed in the study. It was broadly presumed that the control input for the Type 1 and Type 2 actuators would have a similar tendency with the GRF data in 1 g for the former and the difference in ID results between 0 g and 1 g for the latter. As a result of simulation, we confirmed that control inputs similar to our expectation was obtained.
As mentioned in the Introduction section, previous studies have performed data tracking to predict muscle forces in given motion, or predictive simulation to predict motion and corresponding muscle forces to minimize or maximize a given motor task. However, unlike those previous studies, the problem we dealt with in this study was to simultaneously track the joint angles and muscle forces that occur in certain movements (1 g of normal walking). Therefore, we constructed the objective function to include the angle, angular velocity, and muscle activity terms. This made it possible to track kinematics data and muscle forces and get corresponding control inputs.
To the best of our knowledge, no previous studies have been conducted on simulations of systems such as Type 1 and Type 2 for resistance training in space. The results of this study provide a control input that can be used as a reference for the design of resistance training devices in microgravity environments.
In this study, we used a simplified model with one degree of freedom in each joint of the lower extremities to reduce the simulation time. Therefore, it could be approached by the algebraic method or static optimization method. However, the dynamic optimization method used in this study may contribute more to future works. For example, it is a case of extending the degree of freedom of a model to predict the effect of the degree of freedom occurring in directions other than the sagittal plane. In addition, it also includes a simulation to predict changes in motion due to adaptation when the human body is wearing the device. By using the dynamic simulation code that we distribute with this study, readers will be able to perform simulations for the aforementioned cases.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3417/10/3/1160/s1, Figure S1: Tracking results of the generalized coordinates and generalized speeds (Type 1), Figure S2: Tracking results of the muscle activation (Type 1), Figure S3: Results of the fiber length and muscle excitation, Figure S4: Assistive forces and torques of the external actuators (Type 1), Figure S5: Tracking results of the generalized coordinates and generalized speeds (Type 2), Figure S6: Tracking results of the muscle activation (Type 2), Figure S7: Results of the fiber length and muscle excitation (Type 2), Figure S8: Assistive torques of the external actuators (Type 2), Figure S9: Net joint torques by the muscles (Type 1), Figure S10: Net joint torques by the muscles (Type 2), Raw data S1 Simulation results of Type 1 and Type 2, Code S1 Code for regenerating the optimization results.