## 2 Materials and Methods
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:  
\[J=\sum_{i}\int_{0}^{T}e_{i}(t)^{2}dt,\quad e_{i}(t)=X_{i,ref}(t)-X_{i}(t) \tag{1}\]  
Figure 1: Processes used in OpenSim and MATLAB. (**a**) Process for obtaining tracking objective data. (**b**) Optimization process using MATLAB and OpenSim. where \(e_{i}(t)\) is the error between the target variables \(X_{i,ref}(t)\) and design variables \(X_{i}(t)\) for the \(i\)th 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 \(\underline{q}\), generalized speeds \(\underline{q}\), muscle activations \(\underline{a}\), fiber lengths \(\underline{l}\), and control variables: muscle excitations \(\underline{u}\), external actuator torques \(\underline{\tau}\) or forces \(\underline{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 (\(\alpha,\beta,\gamma\) 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 \(\underline{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 (\(\underline{\tau}\)\(\underline{min}\) and \(\underline{f}\)\(\underline{min}\)) and maximum (\(\underline{\tau}\)\(\underline{max}\) and \(\underline{f}\)\(\underline{max}\)) values were set considering the location and type of the actuator. \[(1-\alpha)\underline{q}_{ref}(t) \leq\underline{q}(t)\leq(1+\alpha)\underline{q}_{ref}(t) \tag{2}\] \[(1-\beta)\underline{q}_{ref}(t) \leq\underline{q}(t)\leq(1+\beta)\underline{q}_{ref}(t)\] (3) \[(1-\gamma)\underline{a}_{ref}(t) \leq\underline{a}(t)\leq(1+\gamma)\underline{a}_{ref}(t)\] (4) \[0.25\underline{l}_{opt} \leq\underline{l}(t)\leq 1.25\underline{l}_{opt}\] (5) \[0 \leq\underline{u}(t)\leq 1\] (6) \[\underline{\tau}_{min} \leq\underline{\tau}(t)\leq\underline{\tau}_{max}\] (7) \[\underline{f}_{min} \leq\underline{f}(t)\leq\underline{f}_{max} \tag{8}\]  
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].