Multibody dynamics and optimal control for optimizing spinal exoskeleton design and support

In the industrial work environment, spinal exoskeletons can assist workers with heavy lifting tasks by reducing the needed muscle activity. However, the requirements for the design and control of such an exoskeleton to optimally support users with different body builds and movement styles are still open research questions. Thus, extensive testing on the human body is needed, requiring a lot of different sophisticated prototypes that subjects can wear for several hours. To facilitate this development process, we use multibody dynamics combined with optimal control to optimize the support proﬁle of an existing prototype and evaluate a new design concept (DC) that includes motors at the hip joint. A dynamic model of the prototype was developed, including its passive elements with torque generation that accounts for potential misalignment. The human-robot interaction was simulated and optimized in an all-at-once approach. The parameters that describe the characteristics of the passive elements (including beam radius, spring pretension, length of the lever arm, radius of proﬁle) and, in the case of DC, the torque proﬁles of the motors were optimized. Limits on interaction forces ensured that the exoskeleton remains comfortable to wear. Simulations without the exoskeleton allowed comparing the user’s actuation concerning joint moment and muscle activation. Our results agree well with experimental data using the prototype, making it a useful tool to optimize exoskeleton design and support and evaluate the effect of different actuation systems, mass distributions, and comfort requirements.


Introduction
Low-back pain is one of the most common musculoskeletal disorders with an estimated lifetime prevalence of 77-84% [1]. High one-time or repetitive moderate spinal compression forces can cause low-back pain [2,3] and should be avoided in the work environment. Since the calculation of spinal compression forces requires the consideration of a large number of muscles and is therefore very complex, indicators on joint moment level are preferably used. The authors of [4] showed that the moment acting at the L5/S1 lumbosacral joint (in the following referred to as lumbar joint) scales linearly with the spinal compression forces. Hence, indicators based on joint moments can be used without losing any information. The most important indicators in that regard are the cumulative low-back loads (CLBL), the area under the lumbar moment curve, and peak low-back loads (PLBL) [5,6].
Industrial work environments, in particular, often involve heavy lifting and remaining in awkward static positions for long durations [7] and carry a high risk of low-back pain [6]. To solve this problem, exoskeletons are gaining popularity in the industry to improve working conditions. Spinal exoskeletons reduce the muscle activity that the user requires to perform heavy lifting tasks, and with that, the risk of low-back pain [8,9]. A variety of spinal exoskeletons were developed such as Laevo (Laevo, Rijswijk, Netherlands), Robomate [10], wearable stop-assist device (WSAD) [11], an on-body personal lift augmentation device (PLAD) [12], muscle suit [13], backX [14], and an exoskeleton using pneumatic actuators [15].
However, the development of such a wearable robot is a challenge. It should provide forces and torques that support the user as best as possible during lifting tasks but also not hinder other movements, such as walking and stair climbing. The support depends not only on the deflection of the limbs but also on certain time events, e.g., during the lift-off of the object. However, the exoskeleton cannot exert arbitrarily high forces as this would lead to a low user acceptance rate because it is uncomfortable to wear or could even cause injuries. Extensive testing of the device is needed to meet these requirements, which takes a lot of time and requires a lot of different prototypes. Many biomechanical tests of different devices on subjects were performed in the past, e.g., [8,10,16].
Such experiments give a good overview of the current effect and wearability of the prototype. However, they cannot directly provide ideal support profiles specified by the characteristics of passive elements or motor torque profiles. Thus, it is desirable to use optimization methods that facilitate this design process by providing optimized support profiles. The authors of [17] optimize the support of a lower-limb exoskeleton by applying a human-in-theloop method. Subjects had to wear the device for a long time while an iterative optimization algorithm evaluated the performance and adjusted the support profile. It is difficult to transfer this method to spinal exoskeletons because subjects would have to lift objects repeatedly for a long time while the robot exerts high forces on them.
This work presents an alternative approach to optimizing an exoskeleton in simulation using multibody system dynamics and optimal control. A similar approach was used in [18] for optimizing a lower-limb exoskeleton for walking tasks, in [19] for optimizing a lowerlimb exoskeleton for jumping tasks, or in [20] for identifying optimal stiffness values for an ankle-foot orthosis. In [21], multibody system dynamics and optimal control was used to simulate the user-robot interaction during crutch walking with a lower-limb exoskeleton. In [22], a framework was developed allowing to create efficiently biomechanical and exoskeleton models for simulation studies. Similar to these approaches, we modeled the user and exoskeleton as multibody systems. Using optimal control, an all-at-once approach was used to optimize the exoskeleton's support profile and to compute the user's actuation still needed to reproduce recorded unassisted lifting motions. We chose unassisted lifting motions because we want to optimize the support without the user changing movements. In contrast to [18,19,21] and previous work [23,24], which used a simplified model of an exoskeleton, the prototype of an existing spinal exoskeleton [25] was modeled, simulated, and optimized. The exoskeleton model was coupled with the user in a way that allows to calculate the contact forces acting at the interfaces. In addition, a new design concept that includes motors at the hip joint was optimized and evaluated. Four different scenarios were calculated and compared: human-only simulations without the exoskeleton, simulations with the exoskeleton using the prototype configuration, optimization of the prototype configuration, and optimization of the configuration and actuation patterns of the new design concept. The current paper is an extended version of our conference paper [26]. There, we presented preliminary results on optimizing the prototype and a new design concept. We improved these results by a more thorough examination of the experimental setup and transferred the findings to the misalignment approximation. We further illustrate the method and the numerical results in more detail, which was not possible in [26] due to lack of space: Additional information is provided on the placement of markers and EMG in the experiments, on the passive elements in the prototype, and on the formulation of the optimal control problem, especially on the objective function used. The results section now includes a table with the optimized parameters and an analysis of the user actuation. The discussion section is extended based on the additional results, which are compared to experimental data using the prototype.
In Sect. 2, a brief overview of the experiments performed for the reference lifting motions and the approximation of the passive element state, taking possible misalignment between user and exoskeleton into account, is given. Then, the models of the user and the exoskeleton used for the optimal control problem are described in Sect. 3. The optimal control problem itself is formulated in Sect. 4. The results of the optimization are presented in Sect. 5. A discussion and a short conclusion in Sect. 6 and 7 conclude this paper.

Experiments for modeling and optimizing the spinal exoskeleton
Two different biomechanical experiments were the basis of this work. In the following, a short outline of the setup and procedure is given. From the first experiment, subject-specific models and the kinematics of lifting motions that are later used in the optimization were extracted. The second experiment was for gaining knowledge about the alignment between the user and the exoskeleton to achieve higher accuracy in the simulation.
An active motion capture system (Certus Optotrak, Northern Digital, Canada) was used for both experiments. Unlike passive systems that use retroreflective markers, their markers are equipped with LEDs that are actively turned on and off. This prevents the markers from being swapped by the software in the event of overlapping, significantly reduces jittering of the markers, and achieves a high resolution (down to 0.1 mm). A disadvantage is that only static movements can be recorded since the markers are not free but are connected to a system control unit by cables. Another difference to passive systems is a lower possible frame rate depending on the number of markers used. This results from the fact that not all marker positions are recorded per camera shot but only a subset. Lifting movements are static and not highly dynamic. Thus, they are not affected by these limitations and active marker-based systems have been applied successfully for decades to record them and other dynamic movements [4,[27][28][29].

Experiments for the reference lifting motions
Recorded stoop-lifts of five subjects (male; age 21-36 years; weight 60-82 kg; height 1.70-1.82 m) were used in the optimization. Stoop lift means that the movement is performed with as little knee flexion as possible, but subjects were instructed to bend their knees if it was otherwise too uncomfortable for them. A 10 kg box was placed directly in front of the subjects on a 0.3 m high pedestal. The box was equipped with handles so that the subjects could easily grab them. The subjects were instructed to grasp the box immediately and firmly with their hands to reduce possible movements with their fingers.
An active motion capture system (Certus Optotrak, Northern Digital, Canada) recorded the marker positions using two cameras at a frequency of 44 Hz. Markers grouped in clusters of three on rigid braces were placed on each segment of the subjects, excluding the hand: feet, shank, thigh, pelvis, middle trunk, upper trunk, head, upper arm, and lower arm. The box was tracked with two marker clusters. The use of marker clusters has the advantage in that they can be attached more flexibly to the segments of the subjects. This means that their position can be selected so that they interfere as little as possible with the recorded movements and are optimally visible to the cameras. Two force plates (Kistler Instrumente GmbH, Switzerland) measured the ground reaction forces of the subject and the box at 1000 Hz. In addition, uni-directional (in the vertical direction) force sensors were installed between the box and its handles, recording the forces applied by the subject. The muscle activity of the subjects was recorded with 12 EMG sensors. The considered muscles were the erector spinae muscles, the external oblique muscles, the gluteus maximus muscles, the biceps femoris muscles, and the vastus lateralis muscles, which are all strong muscles located at the torso or the thigh and are extensively used during lifting movements. For fitting the human actuation models described in Sect. 3.1.3, it is beneficial to have an estimate of the maximum muscle activation observed during the recorded lifting movements. For this purpose, maximal voluntary contraction (MVC) trials were performed during the experiments to measure the EMG signal of the muscles during full contraction. Two different MVC trials were performed to measure the strength of the lumbar extensors and flexors as well as the hip extensors: 1. The subject stood on a force plate and performed a maximum effort deadlift against a fixed barbell. 2. The subject was restrained against an instrumented backrest using heavy Velcro straps, and then we had them extend, then flex their torso against the device.
The marker and force plate data was filtered using a bidirectional second-order Butterworth filter with a cut-off frequency of 5 Hz. The EMG data was low-pass filtered with a thirdorder Chebyshev type I filter with 1 dB of passband ripple and a passband edge frequency of 440 Hz. These experiments were also presented in [24] and in [26] in shorter form.

Experiments for the alignment between user and spinal exoskeleton
Three young, healthy, male subjects (average age: 28 years, average height: 177 cm, average weight: 71.3 kg) performed various tasks (lifting, crouching, lateral bending, torso rotation) while wearing the prototype. The tasks were performed with locked and unlocked misalignment mechanisms (three-joint segment and slider, see [25]). For our regression analysis described in Sect. 3.3, we considered only the lifting tasks consisting of free lifting, lifting in a stooped position, and lifting in a squat position with unlocked configuration as this represents the default behavior of the exoskeleton. Active markers were attached to both the subject and the exoskeleton. Their positions were recorded with a three-camera active motion capture system (Certus Optotrak, Northern Digital, Canada) at a frequency of 50 Hz. Markers for the exoskeleton were attached to the slider connected to the beams at the end of the beams, the part of the pelvic module to which the beams are connected, the passive element at the hip joint, and the slider of the thigh interface. The state of the exoskeleton was additionally recorded via the built-in sensors at a frequency of 100 Hz. Similar to the previous experiment, clusters formed by three markers were placed on the shank, thigh, pelvis, and torso of the subject. Due to the exoskeleton, only one cluster could be placed on the back of the torso in contrast to two in the previous experiment. The marker data was filtered using a bidirectional second-order Butterworth filter with a cut-off frequency of 5 Hz.

Combined multibody system of the human and the spinal exoskeleton
The recorded lifting motions that were used in the optimization are fairly symmetrical. Thus, we could reduce the complexity of the system, and hence that of the optimization, by modeling the human, the exoskeleton, and the box as symmetric rigid multibody systems in the sagittal plane ( Fig. 1 Left). First, a brief overview of the human model with muscle torque generators as joint actuators is given. Then the model of the exoskeleton prototype with its passive elements and the additional actuators for DC, as well as the misalignment approximation, is described.

Modeling the human with the box
The human model consists of 11 (eight internal) degrees of freedom (DoF) (Fig. 1 Left). Both arms and legs were lumped together and the trunk was divided into three segments: pelvis, middle, and upper trunk. Every human model was adjusted to represent the properties of each subject of the used motion capture recordings. The dynamic properties (segment Fig. 2 Snapshots of the fitted motions onto the subject-specific models during the contact with the box on the pedestal. From left to right: S1, S2, S3, S4, S5 mass, center of mass, and inertia tensors) were estimated via the regression equations proposed by de Leva [30] and linearly scaled by the segment lengths of each subject measured during the experiment. All internal joints were modeled as revolute joints. The exoskeleton was connected to the human at three locations: at the thigh, slightly above the knee; the pelvis, a little below the lumbar joint (L5/S1 lumbosacral joint); and the upper trunk, slightly below the level of the shoulder joint. The location of the contact point on the hand was determined to match the contact point on the handle during contact with the box. The box was modeled to match the one used in the experiments (refer to Sect. 2).

Kinematic fitting of the reference motions
The experiments provide the location of the recorded markers in space. These must be transferred to joint positions for the models used in the optimization. For this purpose, we adopted the method of [31]. The individual segments of the human model were equipped with virtual markers. Their exact locations with respect to the associated segment were determined using a static trial in which the person remained in a given position. Due to the placement of the markers in our experiments, it was not sufficient to minimize the deviation between virtual and recorded marker position as in [31]. The markers were placed in groups of three on rigid braces to form a cluster. Since the distance between the markers on the cluster was small, small deviations between the position of virtual and recorded markers can lead to high deviation in the orientation of the virtual and recorded clusters. Therefore, we extended the method by minimizing a combination of the deviation between virtual and recorded marker positions and between virtual and recorded cluster orientations, and thereby, the joint positions q ∈ R 11 were determined: with m V i (q) and m C i ∈ R 3 the position of the virtual and motion capture markers. The matrices R V i and R C i ∈ R 3×3 describe the orientation of the body segment and the corresponding marker cluster in the global frame, respectively. N p and N r specify the set of marker and cluster indices whose position and orientation should be matched. The function f angle_axis (R) calculates the angle-axis representation of rotation matrix R. Each term is weighted by the factors a i and b i ∈ R. In Table 1 the accuracy of the kinematic fitting is given in terms of position and angle errors, and in Fig. 2, snapshots of the fitted motions to the models during contact with the box on the pedestal are shown. Please note that the same task was performed with different levels of limb deflection.

Muscle torque generators
Pairs of agonist and antagonist muscle torque generators (MTG) actuate the human model used in the optimization. One MTG represents the muscular properties of a joint in one di- rection. For models reduced to the sagittal plane, this refers to joint flexion and extension.
Since we can estimate the risk of low-back pain at the joint level, a complex system of muscle models was not needed for our computations. The benefit of MTG is that they still cover muscular properties such as the consideration of passive muscle forces, which is important for lifting motions as shown in Sect. 5. Joint damping was added to reduce vibrations in the system. The torque generated at the human joint i is calculated by The torque generated by the MTG for flexion and extension is denoted by τ FL and τ EX , respectively. The joint velocities are given byq. For each MTG, τ 0 is the maximum isometric torque and ω max is the maximum angular velocity. The damping coefficient β is scaled by a factor η.
The amount of torque a muscle can produce depends on the level of activation, the joint position, and the joint velocity. This dependency is modeled in the MTG by three curves: the active (f A ) and passive (f P ) torque-angle curve representing the active and passive forces generated by the muscles and the curve f V describing their torque-velocity relationship: with muscle activation α ∈ [0, 1] and nonlinear normalized damping term β P . The parameters s P and θ P 0 scale and shift the passive torque-angle curve. The MTG-specific angle θ and velocity ω are derived from the joint angle q and velocityq, respectively. It should be noted that only uniarticular muscles are currently supported by the MTG. Thus, the effect of biarticular muscles on movement performance was neglected during simulation. Further information regarding the MTG can be found in [23,32].

Fitting of the muscle torque generators to the experimental data
The possible strength of an MTG at a given position and velocity can be adjusted by a set of parameters. For optimization, it is vital that the muscle models can generate the torques necessary to perform the motion. This is not guaranteed when using the default properties derived from experimental data in the literature as every person differs in strength and flexibility.
A muscle-fitting routine [32] was applied to adjust the MTG so that they can reproduce the recorded motions with muscle activation not exceeding a given range [0, α max ]. Since the subjects likely did not need their full strength to lift a 10 kg-heavy box, an educated guess for the maximum activation level during the lifting motion is helpful to get a more accurate representation of the muscle group properties for each joint. For this purpose, the available EMG data of the recorded lifting motions were normalized, and its peak value was taken as the activation limit (α max ) of the MTG responsible for hip and lumbar extension during the muscle fitting process: with EMG trial the EMG data of the lifting trial, EMG quiet the mean of the EMG data during the quiet standing phase, and EMG MAX MVC the peak value observed during a set of MVC trials explained in Sect. 2. All subjects recorded were able to lift a 15 kg box without having to change their lifting style. Therefore, we assume that no muscle group was close to being fully activated when lifting the 10 kg box, and we chose 0.7 as value for α max in the case that no EMG data were available. All α max values obtained from the normalized EMG data were below 0.7, further supporting this assumption.

Modeling the spinal exoskeleton
Analogous to the human model, the exoskeleton was reduced to the sagittal plane with nine DoF (six internal) and can be organized in several modules: pelvis module, thigh bar (including the passive element), thigh interface, torso bar (set of beams), and torso interface, which are illustrated in Fig. 1. It weighs 6.7 kg, and the dynamic parameters, i.e., mass, the center of mass, and inertia tensors, were derived from CAD models of the existing prototype. The torso and thigh interface are connected by prismatic and revolute joints. The prismatic joints represent the sliders of the system and the deflection of the beams. The orientation of the interfaces and the respective human body parts are aligned via the revolute joints. Three carbon fiber beams running along the back and passive elements with a nonlinear torqueangle relationship [33] at the hip provide counter-torques. The following sections describe mathematical models that reproduce their torque-deflection relationships and were used in the optimization.
In the case of DC, the exoskeleton model was extended by a motor at each hip joint to evaluate their effect on the support and contact forces between the user and the exoskeleton. We expect they will improve the alignment by reducing the contact moment at the pelvis module. The motor is modeled as a simple torque source τ M in addition to the torques τ PH (α, p) generated by the passive element at the hip joint: with α the deflection angle and p the parameters specifying the torque-angle relationship of the passive element, which are described in the next section. These parameters were later optimized or set to the respective values used in the prototype (listed in Table 6). A detailed description of the prototype can be found in [34]. As the motors are intended to be installed on the pelvis module, its weight is increased by 3 kg. A shift of the center of mass or a change in the inertia due to the additional actuators was neglected during this study, but can be implemented easily if such information is available.

Modeling the passive element at the hip joint
Passive elements (PH) are installed at the hip joints of the exoskeleton. A detailed description of it can be found in [33]. The exoskeleton version consists of a linear spring and a profile shaped like a half-heart. A cable that is guided through a system of rollers and along the half-heart-shaped profile connects the center of the revolute joint with the linear spring. When the joint rotates, the cable is bent over the profile yielding a nonlinear torque-angle relationship. The linear spring itself can be pre-tensioned. The following five parameters p (see Eq. (5)) specify the shape of the torque profile: (S) linear spring stiffness, (P ) linear spring pretension, (B) length of the lever arm, (C) distance between rollers and cable attachment point at the joint, (R) radius of the half-heartshaped profile. Except for (S), all parameters were later optimized. The main differences between the version described in [33] and the one installed in the exoskeleton are that no motor is included to adjust the pretension (P ) automatically and that the profile has the shape of a half heart instead of a full heart, i.e., the exoskeleton only provides counter torques during hip flexion, not extension. The torque τ generated at deflection angle α is derived as follows: with Fig. 4 Cantilever beam [26] Because of geometrical reasons, the following constraint must hold to assure that the profile is fabricable: By altering the parameter values C, B, R, E, and P , different torque-angle profiles can be realized. The resulting profiles for some specific configurations are illustrated in Fig. 3.

Modeling the carbon fiber beam
The carbon fiber beam is modeled as a long, thin, cantilever beam (Fig. 4) that has a uniform, circular cross-section and is made of a linear, elastic, homogeneous, and isotropic material. At its free end, a force F acts in a direction determined by the angle α ∈ R. Since the lifting movement is limited to the sagittal plane, only sagittal deflections of the beam are taken into account. Furthermore, we assume that the Bernoulli-Euler hypothesis holds. Applying the analysis of the deflection of a cantilever beam proposed by [35], the Bernoulli-Euler relationship between bending moment and curvature for this type of beam at a point s along the beam with Cartesian coordinates (x(s), y = v(x(s))) ( Fig. 4) can be formulated as follows: where M and κ = dφ ds are the bending moment and curvature, respectively. The moment of inertia of a beam with circular cross-section specified by radius r is denoted by I = π 4 r 4 . The horizontal and vertical positions of the beam's end are given by L x and L y , and α is the angle between the direction of the force and the neutral position of the beam when no forces are applied. Between the curvature dφ ds of the deformed beam and the transverse displacement v(x(s)), the following relationship can be established: By differentiating (9) with respect to s and taking into account the relations dx ds = cos(φ) and dy ds = sin(φ), we obtain This equation together with dx ds = cos(φ) and dy ds = sin(φ) can be integrated and solved for F and L when taking the following boundary conditions into account: with L the beam length. Note that in our case, the horizontal and vertical position of the beam's end is known, whereas the generated force and the length of the beam are unknown since the connector of the torso interface can slide on a system of rollers along the beams.
Polynomial approximation of the beam deflection: The above-described method includes solving a boundary value problem that is too complex to be included in the optimal control setup described in Sect. 4. To avoid the integration of the system, we approximated the deflection of the beam by a polynomial v(x) ≈ P (x) = N i=0 a i x i of order N ∈ N. The coefficient of the polynomial can be reduced by applying the boundary conditions v(0) = 0 and dv dx (0) = 0 yielding a 0 = a 1 = 0. The remaining coefficients and the force F are computed by minimizing the distance to the deflection L y at the end and the deviation to the deflectionmoment relationship (9) at equally spaced gridpoints 0 = x 1 < x 2 < · · · < x M = L x , M ∈ N. By setting M = N − 1, the problem reduces to solving the following set of equations: with Note that the solution for N = 3 corresponds to the so-called linear model in case of The accuracy of the polynomial approximation for α = 90 • and different order N is given in Table 2. The Newton method was applied to solve the set of equations (13)- (14), and the solutions are compared to the nonlinear model solved as a boundary value problem. The solution of the linear model served as the initial guess for the Newton method. Based on these results, N = 8 was used in the optimal control setup, and the characteristics of the beam were optimized by adjusting the radius r of the beam cross-section.

Coupling of human and spinal exoskeleton
It is very difficult to design a device that remains perfectly aligned with the user throughout lifting motions which involve high deflections of the limbs. The high forces/torques that the exoskeleton applies further contribute to it. Misalignment influences the amount of support provided by the passive elements and the contact forces. Thus, it is desirable to have an accurate contact model for simulation. However, this requires suitable experimental data from the prototype, which was not available. In the experiment presented in Sect. 2.2, the movements of the arms and head, as well as the position of the box and contact forces between user and exoskeleton, were not recorded, making the data unsuitable for identifying a contact model using whole-body simulations. Thus, in this work, we assume a rigid coupling to the human at three locations (pelvis, upper trunk, and thigh) via loop closure constraints (please refer to Fig. 1), neglecting possible movement between the user and the exoskeleton. This approach was also used in our previous work [24,36,37] and is explained briefly here. Coupling the human with the exoskeleton imposes kinematic loops on the multibody system, which can be handled in the following way, starting with the equation of motion for the spanning tree derived from the closed-loop system: with M, the mass matrix, C containing the centrifugal, gravitational, and Coriolis force. The joint positions, velocities, and acceleration of the system are denoted by q,q, andq, respectively, the generalized forces by τ . The terms representing the forces exerted on the tree by loop constraints are added: with τ c constraint forces and τ a active forces, which are zero in this case since the interfaces between user and exoskeleton are not actuated. In addition, the loop-closure joints impose a set of kinematic constraints on the tree, with G the constraint Jacobian and k = −Ġ(q)q. Combining (19) and (20) leads to the following system of equation, which is solved for the joint accelerationsq and the unknown forces λ from which the contact forces are calculated: A more detailed description of the modeling of multibody systems containing kinematic loops can be found in [38]. Nonetheless, the experimental data presented in Sect. 2.2 provide an analysis of how the deflection of the human torso and thigh differs from the deflection of the corresponding interfaces of the exoskeleton. This data was used to approximate the state of the passive elements with a linear regression model. In [26], preliminary results were presented of the regression analysis. Further analysis of the data showed that there was an offset in the pelvis segment orientation that distorted the calculated deflection angles. This is now corrected.  Fig. 1 Left). In [26], the deflection of the middle and upper trunk (θ H 7 + θ H 8 ) was used as an input for θ H T in the optimization. A more detailed examination of the experimental setup led to the conclusion that the deflection of the middle trunk alone better matches the experimental data. No markers were placed on the upper trunk and can therefore not be considered during the linear regression analysis. In Table 3, the optimized linear regression parameters can be obtained, and the results show that a significant linear regression relationship exists between the human and exoskeleton deflection angles.
In addition, several studies have highlighted as an important requirement for exoskeletons that they must be comfortable to wear for hours, even though they are firmly attached to the user with straps and plates and exert push and pull forces [10,39]. We address this requirement by setting limits on the interaction forces between the user and the exoskeleton. These thresholds are based on [40], in which several subjects were exposed to constant and repetitive pressure, and in both cases rated at what level it was uncomfortable for them and at what level they felt pain. The average values for discomfort (Table 4) served as limits in the optimization. The optimal control problems for all scenarios are similar. The lifting motion was organized in three phases: The first phase begins with the user standing upright. Then the user bends down and makes contact with the box leading to the next phase. There, the user generates enough force to lift the box. The last phase starts when the box leaves the ground and ends when the user is in an upright position again while holding the box. The OCP is formulated as follows. For better readability, the time dependencies are omitted:

Exoskeleton design optimization as optimal control problem
with q,q, andq denoting the joint positions, velocities, and accelerations, respectively. The number of shooting nodes t i,0 , . . . , t i,N i of phase i is denoted by N i . The algebraic states z and the system of equations (29) define the state of the beam. The parameters p describe the design of the passive elements of the exoskeleton (parameters of PH and beam radius). The controls u are the neural excitation of the MTG. In the case of DC, the derivative of the motor torque profile is a control as well. Equation (27) is the MTG activation dynamics proposed by [41] with activation level α and (de-)activation time constant T m . The number of MTG is given by N m . The equation of motion of the constrained multibody system is given by (25) with mass matrix M, constraint Jacobian G i , and unknown force variables λ. The function C contains the centrifugal, gravitational, and Coriolis forces. The generalized forces are denoted by τ consisting of the joint torques and forces generated by the MTG and the exoskeleton.
The subjects were instructed to directly grasp and lift the box with as little finger movement as possible. This results in a certain impact when the object is grasped. Note that the impact is prescribed by the reference movement. We thus modeled the contact with the box (26) as an inelastic impulse withq + andq − the joint velocities directly before and after the impact. Inelastic refers to a vanishing of the velocities after the impact, which occurs when grabbing the box The constraints are grouped by equality (30) and inequality constraints (31). They further can be distinguished between constraints that should hold at specific moments and over the whole phase duration. At the beginning of the motion, the positions of the feet, the box, and the position and orientation of the exoskeleton interfaces with respect to the human model are specified. When making contact with the box, the hand contact points should match the contact point at the handle of the box. No forces should be applied to the box. During liftoff of the box, the normal forces between box and ground should vanish. During the phases, constraints on box-to-floor (phases 1 and 2), hand-to-box (phase 2), foot-to-floor contact forces (all phases), and interaction forces between user and exoskeleton (all phases) should be maintained. Limits on parameters, states, and controls are set as well. The OCP was discretized using direct multiple shooting, and the resulting NLP was solved with SQP and the active-set method provided by the toolbox MUSCOD-II [42]. For the rigid multibody dynamics calculations, the open-source library RBDL [43] was used.

Objective function
The objective function (24) consists of a least-squares term for tracking the recorded motion and a Lagrange term that mainly enforces the reduction of human joint moments. The motion to be tracked is given for time t i,n by the joint positions q REF i,n , and the fitting accuracy is specified by a weighting matrix W q . The translational DoF are weighted with 90 N , the rotational DoF with 30 N , N being the total number of shooting nodes. The function φ can be split into several parts: withτ a MTG andτ b MTG being the active and passive MTG-torques normalized by the maximum isometric torque, respectively. The moment acting at the pelvis contact point is denoted by m p , the torque and the torque change applied by the motor by τ M andτ M , respectively. The factors c a = 0.2 , and c u = 2×10 −7 T were applied to weigh each term of the objective function with T the duration of the motion. The matrix W τ weighs the MTG torques at the lumbar joint with 2, the remaining joints with 1. The scaling factors were selected manually in an iterative manner, yielding a high tracking accuracy and balanced force distribution while maximizing the support.
The first two terms reduce the joint moments of the user and, at the same time, shape the exoskeleton support. The separate consideration of active and passive MTG torques was selected because in a previous cost function comparison [24] this objective resulted in the best support while maintaining reasonable actuation of the human model concerning exploitation of passive MTG torques or co-contraction. The third and the last two terms are small regularization terms to eliminate redundancies. In practice, the usage ofα instead of u yielded a better convergence behavior. The fourth term results in a balanced distribution of the applied forces by the exoskeleton. In the case of H, the objective φ only consisted of the third term, in case of O, of the first three terms, in case of PO, of the first four terms, and in case of DC, of all terms. The last term was not included in [26]. It was added because sharp peaks in the motor torque profile were observed during contact with the box. The weighting was the same for all cases.

Numerical results
The exoskeletal support was optimized for several recorded lifting motions. For the prototype optimization (PO), parameters describing the behavior of the beam (beam cross-section radius) and the passive elements (see Sect and only small difference (avg. 0.14 • ) to H. This ensures that the reduction in muscle activity is due solely to the support provided by the exoskeleton and not to any variation in movement. In Table 5, the achieved support for each subject with respect to each setup is listed. Both optimized configurations, PO and DC, provide higher support than O across all subjects. However, PO increases the support by additional 2.5-4.1% regarding CLBL, 0.8-2.1% regarding the hip moment, and 0.9-3.7% regarding PLBL compared to DC. Please recall from the introduction that it is sufficient to look at the joint moment, a calculation of the spinal compression forces is not necessary to evaluate the risk of low-back pain. Compared to the previous results in [26], the optimized reduction of DC and PO remained similar, but the change in the misalignment approximation led to a significant lower support of O across all subjects, especially at the lumbar joint.
The higher support of PO is also evident when looking at the torque profiles of the beams and PH (Fig. 5). The optimized torque profiles of DC have a lower peak. With the motors, a broader torque curve is achieved, resulting in a smaller loss of support at the hip joint than at the lumbar joint compared to PO. The improvement of the misalignment approximation did not lead to a change in the magnitude of the optimized torque profiles. However, they yield a different profile during the second upright standing while holding the box compared to [26], leading to a shift of the moment at the pelvis contact point as well.
Similarly, DC results in a much higher reduction of the moment acting at the pelvis contact compared to PO (Fig. 6), staying within 4.3 Nm] in contrast to [-5.4 Nm, 6.6 Nm] and yielding a very good alignment of the exoskeleton and the user. For comparison, the contact moment of O lies in the range of 5.5 Nm] for all subjects with the highest value during bent down, which could already be reduced well in PO by adjusting the passive elements alone.
The fact that the limit on the normal contact force at the pelvis is reached for both PO and DC (Fig. 6) suggests that due to the higher weight of the pelvic module, a force/torque generation of the exoskeleton of DC at the same level as PO was not possible before the limit was reached. The higher weight is also reflected in the shear forces acting at the pelvis contact point (Fig. 6), in particular at the beginning and end of the motion when the user is standing upright and the weight of the exoskeleton mainly contributes to them. The normal and shear forces of the thigh (max. 126 N and 7 N across subjects and PO/DC) and torso (max. 118 N and 16 N) contact point stayed far below the set limits.
Analysis of the distribution of active and passive MTG moments at the lumbar and hip joints confirmed that the optimized design did not apply excessive forces. Excessively high forces could result in the human models having to flex against the device to perform the recorded motion. The two subjects with the lowest active MTG extension torque during the time of contact with the box, either at the hip or lumbar joint, are illustrated in Fig. 7 (18 Nm hip extension torque for subject 5; 27 Nm lumbar extension torque for subject 4). A reduction (with respect to H) of peak MTG activation for lumbar extension (Fig. 8)  All optimized parameters are listed in Table 6 with the configuration of the prototype on the right, which is the same for all subjects. The radius of the beam was increased significantly for all subjects. The highest value was optimized for subject 4, who was the only one to perform the stoop lift with his back strongly flexed and without knee flexion (see Fig. 2). His upper body is also significantly longer than that of the other subjects, resulting in a longer beam length and less support. For the remaining subjects, similar values for the beam radius were obtained, indicating that the lifting style and stature influence the amount of support that the exoskeleton provides at the back. There are significant differences in the  The total and the passive torque produced by the MTG at the lumbar and hip joint for subject 5 (first two pictures) and subject 4 (last two pictures) of PO optimized parameter values of PH for DC and PO. In general, they resulted in a higher pretension (P ), a lower distance to the rollers (C), a higher lever arm (B), and a lower profile radius (R). In the case of PO, higher values for B were optimized. The values for R reached the specified minimum value of 5 mm for all subjects, which produces a flatter profile. For DC, R stayed close to the original value for three subjects, and B did not change much.

Discussion
In this work, we presented a method combining multibody dynamics and optimal control to simulate, evaluate, and optimize exoskeleton design concepts during the execution of lifting motions. The optimization was performed as an all-at-once approach. The actuation of the human and the design of the exoskeleton were calculated so that they reproduce and support the recorded lifting motions best. Contact forces between the user and the exoskeleton, the box, and the ground were calculated. The method optimized the characteristics of the passive elements, which were described by a set of parameters, and motor torque profiles that can serve as guidelines for the next prototype. The modular design allows easy replacement of different exoskeleton models. It can also be used to perform optimization without reference data if the tracking term is omitted, shown in [23,44]. Note that this method does not require manufactured prototypes and biomechanical experiments with long periods of device use as in [17]. The exoskeleton was modeled based on CAD data and results of component tests of the passive elements. The misalignment approximation required experimental data, but similar findings could also be obtained by simulations using FE models.
The optimizations of the new design concept showed that better control of the contact moment at the pelvis is possible. However, this is accompanied by lower possible support. Devices with actuators are characterized by higher flexibility in the support profile but are usually significantly heavier than those that have only passive elements. It is questionable whether the improvement in wearing comfort due to the reduction in contact moment justifies the higher weight. Our calculations show that the passive elements can be adjusted so that only moderate contact moments occur during the lifting movements. A benefit of the motors is that they could provide more support during the lift-off of the box when the loading of the user peaks. The support of the passive elements is coupled with trunk and hip flexion and cannot account for that. However, this is currently restricted by the limits on the contact forces. The motors could reduce the joint moment further if higher contact forces are allowed for a short moment.
The coupling between the user and the exoskeleton is a limitation of the current approach. In reality, displacements or rotations of the interfaces occur, which affect the calculated support and contact forces. So far, these effects are only taken into account in the state calculation of the passive elements by using linear regression equations derived from experimental data. These already provide much more realistic calculations of the support but depend only on the position of the user and not on the generated torque distribution of the exoskeleton. Furthermore, they do not take into account the orientation of the upper part of the torso. This may have also influenced the calculated design and support of subject 4. Accurate modeling of human-robot interaction requires additional experimental data, which were not available for these calculations and is part of future work. These experiments are also necessary to validate the simulated force transmission. The modeling of the passive elements is based on experimental observations. For the modeling of, e.g., frictions in the joints of the exoskeleton, experiments with manufactured parts of the prototype are needed.
The distribution of the active and passive MTG torques at the lumbar and hip joint for two subjects (Fig. 7) shows another constraint that has to be considered during the design optimization. Lifting movements push you to the limits of your range of motion. In this work, we considered movements to pick up a box from a 0.3 m high pedestal. One can expect even higher passive forces when the box is placed on the ground. This emphasizes the importance of also considering the passive muscle forces during optimization so that the exoskeleton does not interfere with the performance of required work tasks or, worse, force the user to change movements, which may result in a higher risk of low-back pain.
The optimized beam radii of PO for four subjects are very similar to the beams (with a radius of 2.35 mm) used in [8]. Otherwise, the prototype was the same as in [25] and was evaluated in biomechanical experiments measuring muscle activation and estimating the joint moment and spinal compression forces during similar lifting tasks. The calculated reduction in peak MTG actuation for lumbar extension of PO is only slightly lower than the values (22 ± 5%) reported in [8]. This can be attributed to the fact that during these experiments a slight change in the performance when using the exoskeleton was observed. The trunk flexion stayed the same, which supports our approach to using recordings of unassisted lifting motions. However, a decrease in the trunk velocity occurred. This led to a reduction of the joint moment in addition to the exoskeleton support and likely contributed to an increased reduction of peak muscle activation as well. Nonetheless, the difference is so small that our approach still yields realistic results in terms of the calculated support. Our results also agree with the observation that the moment of maximum support of the exoskeleton does not coincide with the moment of peak loading of the user. As mentioned before, the motors could account for that, but this is restricted by the applied limits on the contact forces.
Our values of peak muscle activation are higher than reported in [8] (28 ± 6.3% MVC). Those values highly depend on the fitness of the user. In the experiments we used, sports students participated in contrast to luggage handlers in [8]. Luggage handlers are likely better trained to perform lifting tasks and need, therefore, less muscle activation. Thus, we can conclude that the similarity between peak muscle and peak MTG activation reduction and the consideration of passive muscle forces as passive MTG torques show that the simplified actuation via MTG is a good representation of the human muscle actuation system.

Conclusion
Multibody mechanics and optimal control are promising tools to simulate the performance of exoskeletons during task execution. For different recorded lifting motions, we were able to determine optimal values for the characteristics of the passive elements of the exoskeleton that maximize support while maintaining contact limits. The evaluation of new design concepts was exemplified by the insertion of motors at the hip joint. Similarly, the effects of different weight distributions or fits can be studied. Future work still needs to be done to validate the simulation and to improve the human-exoskeleton contact model, taking into account misalignment between the user and the exoskeleton not only in terms of the state of the passive elements but also in terms of the positioning and orientation of the interfaces.
Funding Open Access funding enabled and organized by Projekt DEAL. This work was funded by the European Commissions within the H2020 project SPEXOR (GA 687662).

Competing Interests
The authors declare no competing interests exist.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.