2.1. Human-exoskeleton model

A 2D human skeletal model is considered in this work [Reference Zaman, Arefeen, Quarnstrom, Barman, Yang and Xiang31], as shown in Fig. 1(a). The human skeletal model is symmetric along the sagittal plane and has $n=10$ DOFs, including three global DOFs. In Fig. 1(b), it is illustrated that two 1-DOF powered exoskeletons are attached at the knee joints using straps for the experiment. Exoskeleton’s mass and inertia are added to the upper and lower legs mathematically. Therefore, the coupled human-exoskeleton mass and inertia of the upper and lower legs are calculated for the simulation. The human skeletal model is constructed using the robotic formulation of the DH method [Reference Denavit and Hartenberg30]. Each DOF represents relative rotation/translation of two body segments connected by a revolute/prismatic joint. For translational and rotational joints, the DOF is given in the local $z$ -direction. It is noted that the global rotation joint ( $z_{3}$ ), spine joint ( $z_{4}$ ), and hip joint ( $z_{7}$ ) coincide at the same location. The positive directions for all the local rotation joints ( $z_{3}\sim z_{10}$ ) are clockwise in the global Y-Z plane. The two physical branches are the spine-arm branch and leg branch. Two arms are represented in the spine-arm branch by a single branch, since only 2D symmetric lifting is studied. In the leg branch, two legs are combined as a single branch, including thigh, tibia, and foot. In this study, the subject’s body segment length data are generated from the Xsens motion capture system (MVN Analyze Pro software) based on the measured height and foot length. The DH parameters for the human model are described in Table I, where L₁ to L₇ are the human joint link lengths.

Figure 1. (a) The 2D human skeletal model with knee exoskeleton and (b) Knee exoskeleton.

Table I. DH table for 2D human model.

2.2. Equations of motion

In this paper, the kinematics and dynamics of the human model are studied using recursive kinematics and Lagrangian dynamics, and details refer to refs. [Reference Xiang and Arefeen32, Reference Xiang, Arora, Rahmatalla and Abdel-Malek33]. The process is comprised of two parts: forward kinematics and backward dynamics. Forward kinematics transmits motion from the root to the end effectors, while backward dynamics transfers forces from the end effectors to the root. The human dynamics equation can be written as [Reference Xiang and Arefeen32, Reference Xiang, Arora, Rahmatalla and Abdel-Malek33]:

(1) \begin{equation}\tau _{{h_{i}}}=\mathrm{tr}\left(\frac{\partial \mathbf{A}_{\mathrm{i}}}{\partial q_{i}}\mathbf{D}_{i}\right)-\mathbf{g}^{\mathrm{T}}\frac{\partial \mathbf{A}_{i}}{\partial q_{i}}\mathbf{E}_{i}-\mathbf{f}_{k}^{\mathrm{T}}\frac{\partial \mathbf{A}_{i}}{\partial q_{i}}\mathbf{F}_{i}-\mathbf{G}_{i}^{\mathrm{T}}\mathbf{A}_{i-1}\mathbf{z}_{0}\end{equation}

where $\tau _{{h_{i}}}$ is the human torque at ith joint. On the right-hand side of Eq. (1), the first term is inertia and Coriolis torque, the second term is the torque due to gravity, the third term is the torque due to external forces, and the fourth term is the torque due to external moments:

(2) \begin{equation}\mathbf{D}_{i}=\mathbf{I}_{i}\mathbf{C}_{i}^{\mathrm{T}}+\mathbf{T}_{i+1}\mathbf{D}_{i+1}\end{equation}

(3) \begin{equation}\mathbf{E}_{i}=m_{i}\mathbf{r}_{i}+\mathbf{T}_{i+1}\mathbf{E}_{i+1}\end{equation}

(4) \begin{equation}\mathbf{F}_{i}=\mathbf{r}_{k}{\unicode[Arial]{x03B4}} _{ik}+\mathbf{T}_{i+1}\mathbf{F}_{i+1}\end{equation}

(5) \begin{equation}\mathbf{G}_{i}=\mathbf{h}_{k}{\unicode[Arial]{x03B4}} _{ik}+\mathbf{G}_{i+1}\end{equation}

where $tr(\cdot )$ is the trace of a matrix. $\mathbf{A}_{i}$ , $\mathbf{C}_{i}$ are global position and acceleration transformation matrices, $\mathbf{I}_{i}$ is the inertia matrix for link i, $\mathbf{D}_{i}$ is the recursive inertia and Coriolis matrix, $\mathbf{E}_{i}$ is the recursive vector for gravity torque calculation, $\mathbf{F}_{i}$ is the recursive vector for external force–torque calculation, $\mathbf{G}_{i}$ is the recursive vector for external moment torque calculation, $\mathbf{g}$ is the gravity vector, $m_{i}$ is the mass of link i, $\mathbf{r}_{i}$ is the COM of link i in the $i$ th local frame, $\mathbf{f}_{k}=[0 \quad f_{ky} \quad f_{kz} \quad 0]^{\mathrm{T}}$ is the external force applied on link k, $\mathbf{r}_{k}$ is the position of the external force in the $k$ th local frame, $\mathbf{h}_{k}=[h_{x} \quad 0 \quad 0\quad 0]^{\mathrm{T}}$ is the external moment applied on link k, $\mathbf{z}_{0}=[0 \quad 0 \quad 1 \quad 0 ]^{\mathrm{T}}$ is for a revolute joint, $\mathbf{z}_{0}=[0 \quad 0 \quad 0 \quad 0]^{\mathrm{T}}$ is for a prismatic joint, and finally, ${\unicode[Arial]{x03B4}} _{ik}$ is Kronecker delta, and the starting conditions are $\mathbf{D}_{n+1}=[\mathbf{0}]$ and $\mathbf{E}_{n+1}=\mathbf{F}_{n+1}=\mathbf{G}_{n+1}=[\mathbf{0}]$ . The sensitivities to state variables are detailed in [Reference Xiang and Arefeen32, Reference Xiang, Arora, Rahmatalla and Abdel-Malek33].

2.2.1. Exoskeleton dynamics

The electromechanical dynamics of DC motors of the exoskeletons are modeled in this study. The dynamics equations can be expressed as follows [Reference Nguyen, LaPre, Price, Umberger and Sup34]:

(6) \begin{equation}L\frac{dI}{dt}=V-K\frac{d\theta }{dt}-RI\end{equation}

(7) \begin{equation}T_{\textit{motor}}=KI\end{equation}

(8) \begin{equation}T_{l}=T_{\textit{motor}}-J_{m}\frac{d^{2}\theta }{dt^{2}}-b\frac{d\theta }{dt}\end{equation}

where V, I, L, and R are the voltage input, current, inductance, and resistance, respectively. The mechanical terms $J_{m}$ , b, K, and $\theta$ are the rotor moment of inertia, coefficient of viscous friction of the motor, motor torque constant, and rotor angle, respectively. $T_{\textit{motor}}$ is the motor output torque, and $T_{l}$ is the load torque of the exoskeleton. The gearbox ratios ( $GB_{r}$ ) are chosen so that the devices can provide the required torque output. Here, the exoskeleton includes the motor and the gearbox. So, the output torque ( $\tau _{e})$ of the exoskeleton can be expressed as follows:

(9) \begin{equation}\tau _{e}=GB_{r}\times T_{l}\end{equation}

Furthermore, we assume that the exoskeleton’s motion aligns with human joint motion. As a result, the derivative of the rotor angle $\dot{\theta }$ and $\ddot{\theta }$ can be expressed in terms of human joint angle derivative $\dot{q}$ and $\ddot{q}$ with gear ratio, respectively:

(10) \begin{equation}\dot{\theta }=GB_{r}\times \dot{q}\end{equation}

(11) \begin{equation}\ddot{\theta }=GB_{r}\times \ddot{q}\end{equation}

where $\dot{q}$ and $\ddot{q}$ are human angular velocity and acceleration for a joint, respectively.

2.2.2. Human-exoskeleton EOM

The coupled human-exoskeleton EOM and sensitivity analysis are formed using a recursive Lagrangian dynamics formulation. The overall dynamics can be expressed as follows:

(12) \begin{equation}\tau _{{h_{i}}}+\tau _{{e_{i}}}=\mathrm{tr}\left(\frac{\partial \mathbf{A}_{\mathrm{i}}}{\partial q_{i}}\mathbf{D}_{i}\right)-\mathbf{g}^{\mathrm{T}}\frac{\partial \mathbf{A}_{i}}{\partial q_{i}}\mathbf{E}_{i}-\mathbf{f}_{k}^{\mathrm{T}}\frac{\partial \mathbf{A}_{i}}{\partial q_{i}}\mathbf{F}_{i}-\mathbf{G}_{i}^{\mathrm{T}}\mathbf{A}_{i-1}\mathbf{z}_{0}\end{equation}

where $\tau _{{e_{i}}}$ is the exoskeleton output torque for $i$ th joint. The ground reaction forces (GRFs) are calculated from Eq. (1) using a two-step active-and-passive algorithm in refs. [Reference Xiang and Arefeen32, Reference Xiang, Arora, Rahmatalla and Abdel-Malek33].

3.1. Design variables

The motor current $I(t)$ profiles and human joint profiles $q(t)$ are discretized using cubic B-splines [Reference Xiang, Arora, Rahmatalla, Marler, Bhatt and Abdel-Malek35]. The design variables ( $\mathbf{x}$ ) are human joint angle control points $\mathbf{P}_{\textit{human}}$ and exoskeleton current control points $\mathbf{P}_{\textit{current}}$ . As a result, the design variables for the human-exoskeleton motion prediction are $\mathbf{x}=\left[\begin{array}{c@{\quad}c} \mathbf{P}_{\textit{human}}^{\mathrm{T}} & \mathbf{P}_{\textit{current}}^{\mathrm{T}} \end{array}\right]^{\mathrm{T}}$ .

3.2. Objective functions

The objective function is to minimize the summation of the normalized human joint torque squares [Reference Arefeen and Xiang36, Reference Arefeen and Xiang37]:

(13) \begin{equation}\min _{\mathbf{x}} J_{1}\left(\mathbf{x}\right)=\sum _{i=3}^{n}\int _{0}^{T}\left[\frac{\tau _{hi}(\mathbf{x})}{\left(\tau _{i}^{U}-\tau _{i}^{L}\right)}\right]^{2}dt\end{equation}

where $T$ is the specified total time for the lifting task, and $\tau _{i}^{U}$ and $\tau _{i}^{L}$ are upper and lower torque limits for the $i$ th joint, respectively.

3.3. Constraints

Time-dependent constraints include (1) human joint angle limits, (2) human joint torque limits, (3) human feet contacting position, (4) box forward position, (5) collision avoidance, (6) human dynamic stability condition, and (7) exoskeleton torque limits. Time-independent constraints include (8) initial and final hand positions, (9) static conditions at the beginning and end of the motion, and (10) initial, mid-time, and final joint angles of the spine, shoulder, elbow, hip, knee, and ankle. Time-dependent constraints are calculated sequentially in the optimization process at every time discretization point. In contrast, the optimization calculates the time-independent constraints at a specific time.

(1) Human joint angle limits:

(14) \begin{equation}\mathbf{q}^{L}\leq \mathbf{q}\left(\mathbf{x},t\right)\leq \mathbf{q}^{U}\end{equation}

where $\mathbf{q}^{L}$ and $\mathbf{q}^{U}$ are the lower and upper joint limits for human [Reference Xiang, Zaman, Rakshit and Yang38–Reference Zaman, Xiang, Cruz and Yang41], respectively.

(2) Human joint torque limits:

(15) \begin{equation}\boldsymbol{\tau }^{L}\leq \boldsymbol{\tau }\left(\mathbf{x},t\right)\leq \boldsymbol{\tau }^{U}\end{equation}

where $\boldsymbol{\tau }^{L}$ and $\boldsymbol{\tau }^{U}$ are the lower and upper joint torque limits for human [Reference Xiang, Zaman, Rakshit and Yang38–Reference Zaman, Xiang, Cruz and Yang41], respectively.

(3) Human feet contacting position:

(16) \begin{equation}p_{feet}\left(\mathbf{x},t\right)=p_{feet}^{s}\end{equation}

where $p_{feet}^{s}$ is the specified feet contact position on level ground.

(4) Box forward position:

(17) \begin{equation}Z_{\textit{wrist}}(\mathbf{x},t)-Z_{\textit{pelvis}}(\mathbf{x},t)\geq 0\end{equation}

where $Z_{\textit{wrist}}$ and $Z_{\textit{pelvis}}$ are the global Z coordinates of the human wrist and pelvis points [Reference Xiang and Arefeen32].

(5) Collision avoidance:

(18) \begin{equation}d_{\textit{human}}\left(\mathbf{x},t\right)\geq r_{\textit{human}}+\frac{dep}{2}\end{equation}

where $d_{\textit{human}}$ is the calculated distance between the hand and the circle center on the body segment representing the body thickness, $dep$ is the box depth, and $r_{\textit{human}}$ is the radius of the circle filled on human limbs as shown in Fig. 2(a) [Reference Luo, Androwis, Adamovich, Nunez, Su and Zhou27]. There are seven circles on the human limb: five are on the leg (ankle, shank, knee, thigh, and hip) and two on the spine as shown in Fig. 2(b).

Figure 2. Collision avoidance constraint.

(6) Dynamic stability condition:

(19) \begin{equation}p_{\textit{human}\_ ZMP}\left(\mathbf{x},t\right)\in \mathrm{FSR}\end{equation}

where zero moment point (ZMP) position is inside the foot support region (FSR) for human [Reference Xiang and Arefeen32].

(7) Exoskeleton torque limits:

(20) \begin{equation}\boldsymbol{\tau }_{exo}^{L}\leq \boldsymbol{\tau }_{exo}\left(\mathbf{x},t\right)\leq \boldsymbol{\tau }_{exo}^{U}\end{equation}

where $\boldsymbol{\tau }_{exo}^{L}$ is lower torque limit, and $\boldsymbol{\tau }_{exo}^{U}$ is the upper limit for the exoskeleton.

(8) Initial and final hand positions:

(21) \begin{equation}p_{\textit{human}\_ hand}\left(\mathbf{x},t\right)=p_{\textit{human}\_ hand}^{s}\left(t\right);\, t=0, T\end{equation}

where, $p_{\textit{human}\_ hand}^{s}$ is the specified hand position at initial and final times.

(9) Initial and final static conditions:

(22) \begin{equation}\dot{\mathbf{q}}_{\textit{human}}\left(\mathbf{x},t\right)=\mathbf{0};\, t=0,T\end{equation}

(10) Initial, mid-time, and final joint angles for the spine, shoulder, elbow, hip, knee, and ankle:

(23) \begin{equation}\left| q_{i\_ \textit{human}}\left(\mathbf{x},t\right)-q_{i\_ \textit{human}}^{E}\left(t\right)\right| \leq \varepsilon ;\, t=0,\frac{T}{2},T\end{equation}

where ${\unicode[Arial]{x03B5}} =0.2$ rad, and $q_{i\_ \textit{human}}^{E}$ is the experimental joint angle for the spine, shoulder, elbow, hip, knee, and ankle joints.

4.1. Control strategy

We propose a new exoskeleton control strategy which uses the offline solved optimal simulation results under the controlled environment in real time to assist human lifting. This control algorithm has two phases as shown in Fig. 3 and consists of two offline optimizations: lifting optimization and assistive torque interpolation optimization.

Figure 3. The diagram of the proposed lifting assistance control algorithm. Here, $\theta _{{k_{0}}}$ and $\theta _{{k_{f}}}$ are initial and final knee angles, respectively. $\tau _{l}$ and $\tau _{u}$ are the lower and upper limits of the exoskeleton torque. $\tau _{exo}^{*}(\theta _{e})$ is the optimal exoskeleton torque in joint angle domain, and $\tau _{exo}(\theta _{e})$ is the real exoskeleton assistive torque.

In Section 2.1, we modeled a 2D human model with a knee exoskeleton. For the optimization in the first phase, an inverse dynamics optimization formulation is used to find the optimal exoskeleton torque, $\tau _{exo}^{*}(t)$ and optimal knee angle $\theta _{k}^{*}(t)$ for the lifting task, as discussed in Section 3. However, for exoskeleton control, the exoskeleton reads the knee encoder angle ( $\theta _{e}$ ), and outputs the exoskeleton torque, which is a function of the knee joint encoder angle $\tau _{exo}(\theta _{e}).$ Therefore, we need to transfer the discretized optimal exoskeleton torque $\tau _{exo}^{*}(t)$ obtained from the simulation in the time domain to a continuous function of the encoder angle in joint angle domain. In this study, we use a B-spline interpolation to represent $\tau _{exo}(\theta _{e})$ as in Eq. (24):

(24) \begin{equation}\tau _{exo}\left(\mathbf{s},\mathbf{P}_{{\tau _{exo}}},\theta _{e}\right)=\sum _{i=1}^{m}\mathrm{N}_{i}\left(\mathbf{s},\theta _{e}\right){\left(\mathbf{P}_{{\tau _{exo}}}\right)}_{i} \theta _{{e_{0}}}\leq \theta _{e}\leq \theta _{{e_{f}}}\end{equation}

where s = $\left[\theta _{{e_{0}}},\ldots .\,.\,, \theta _{{e_{f}}}\right]^{\mathrm{T}}$ is the encoder angle knot vector, $\mathbf{P}_{{\tau _{exo}}} = \left[{(P_{{\tau _{exo}}}})_{1},\ldots \ldots, {(P_{{\tau _{exo}}}})_{m}\right]^{\mathrm{T}}$ is the exoskeleton torque control point vector, m is the number of control points, $\theta _{e}$ is the current encoder angle, and $\mathrm{N}_{i}(\mathbf{s},\theta _{e})$ is the basis function. Here, $\theta _{{e_{0}}}$ and $\theta _{{e_{f}}}$ are initial (squat position) and final (standing position) exoskeleton encoder angles, respectively. These angles of encoder were obtained from the experiment.

The second optimization is to find the exoskeleton torque control points $\mathbf{P}_{{\tau _{exo}}}$ , to minimize the summation of the error squares of the optimal exoskeleton torque (simulation) and exoskeleton torque from the encoder angle (Eq. (24)):

(25) \begin{equation}\min _{\textbf{P}_{{\tau _{{exo}}}}} J_{2}\left(\textbf{P}_{{\tau _{exo}}}\right)=\sum _{\theta _{e}=\theta _{{e_{0}}}}^{\theta _{{e_{f}}}}\left[{\tau _{exo}}\left({\textbf{s},\textbf{P}_{{\tau _{exo}}}},{\theta _{e}}\right)-\tau _{exo}^{*}\left(t\left(\theta _{e}\right)\right)\right]^{2}\end{equation}

(26) \begin{equation}s.t.\colon \tau _{exo}\left(\textbf{P}_{{\tau _{exo}}},\theta _{e}\right)-\tau _{exo}^{*}\left(t\left(\theta _{e}\right)\right)=0, \theta _{e}=\theta _{{e_{0}}},\theta _{{e_{f}}}\end{equation}

This problem is solved using an SQP method in SNOPT [Reference Shepherd and Rouse23]. The second optimization is able to find the optimal torque control points for the exoskeleton torque $\tau _{exo}(\theta _{e})$ in terms of the encoder angle for lifting tasks.

In phase 2, based on the knee encoder angle during lifting, the optimal assistive torque is calculated from the B-spline interpolation function in real time using Eq. (24) in joint angle domain. This gives human knee joint an assistive torque. It is noted that the encoder angle must be within the given initial and final boundary conditions to provide optimal assistance. If it is outside the knee joint boundaries, a randomly selected small and safe constant assistive torque $\tau _{exo}=2.1$ Nm is provided.

A comparison between exoskeleton torque obtained from B-spline interpolation ( $\tau _{exo}^{*}$ ) and optimal exoskeleton torque profiles ( $\tau _{exo}$ ) are illustrated in Fig. 4. Figure 5 presents the optimal real-time exoskeleton torques with respect to knee joint angles with the maximum assistive torque as 6 Nm and 16 Nm, respectively.

Figure 4. Comparison between B-spline interpolation and optimal exoskeleton torques as a function of human knee angles.

Figure 5. The optimal exoskeleton torque as a function of human knee angles during lifting. The knee joint angle starts from a squat position.

Figure 6. Lifting experiment setup without the exoskeleton.

4.2. Experimental procedure

The Institutional Review Board (IRB) approved the optimal control of the powered knee exoskeleton lifting experiment was conducted at Biodynamics Optimization and Motion Control Lab at Oklahoma State University. The test subject was a healthy 22-year-old male with no previous injuries, and a written consent form from the subject was taken before the experiment. The subject’s height and weight were 1.90 m and 81.36 Kg, respectively.

For this study, the subject had two lab visits. For the first visit, the subject’s anthropometric measurements were taken. Xsens motion capture system was used to obtain 3D kinematic data at 60 Hz, and the whole body sensors protocol was considered, as shown in Fig. 6. The Xsens system was first calibrated. Delsys Trigno EMG sensors were attached to the subject’s left leg and spine. For this experiment, the EMG activities of vastus medialis, vastus lateralis, and biceps femoris from the left leg and latissimus dorsi from the upper body were measured at 2000 Hz. The subject’s maximum voluntary contractions (MVCs) for the four muscles were measured for muscle activation normalization. In this study, three trials were conducted for each muscle at every MVC testing position [Reference Hamzaid, Smith and Davis42, Reference Wagner, Forister and Huot43]. A 30-s rest period was observed between two trials. The participant was instructed to exert his maximum force during each MVC testing session. For lifting task, the subject was instructed to lift an 11-Kg box without the exoskeleton with two feet standing on two Bertec force plates (Bertec, Columbus, Ohio, USA). However, the participant was allowed to choose his own lifting strategy. The lifting motion was recorded in Xsens’s MVN Analyze Pro software (MVN Awinda, Xsens Technologies BV, Enschede, the Netherlands). The GRFs were recorded using OptiTrack Motive 3.0 software (Natural Point, Corvallis, OR, USA) at 1000 Hz. In addition, EMGworks Acquisition (Delsys Inc., Natick, MA, USA) was used to collect and record the raw data of the EMG. The lifting task was repeated three times with a 3-min break between them. Next, the processed data from the MVN Analyze and Motive 3.0 were used to validate the simulation results against the experimental results without the exoskeleton. The recorded EMG data series were transferred to EMGworks Analysis (Delsys Inc., Natick, MA, USA) to obtain EMG activity in millivolts. Then EMG data series were further processed in MATLAB (MathWorks Inc., Natick, MA, USA), where data were bandpass-filtered with the cutoffs at 10–450 Hz, rectified, and then lowpass-filtered with a cutoff frequency of 6 Hz. Finally, the filtered EMG data series were normalized by the corresponding MVC EMG values. For processed GRFs, the average of two force plates was taken at each trial. Finally, we compared the average of the three experimental trials against the predicted data.

For the second visit, a light-weight, highly-back-drivable wearable knee exoskeleton developed by Picasso Intelligence, LLC, was used in this study. The unilateral mass of the exoskeleton was 2.25 Kg, and it could provide the maximum 25 Nm output torque. The exoskeleton’s smart motor included a motor, gear, driver, encoder, and controller. It used CAN bus protocol and was controlled by a teensy-based PCB board with a CAN bus chip. The proposed control strategy from Section 4.1 was implemented, and the code was uploaded to the teensy board of the exoskeleton control board. Finally, the subject was asked to wear Xsens whole body sensors, EMG sensors, and the exoskeleton to perform the lifting task again, as shown in Fig. 7. Similar to the first visit, the lifting motion was recorded in Xsens’s MVN Analyze Pro software, EMG raw data series were recorded in EMGworks Acquisition, and GRFs were recorded using OptiTrack Motive 3.0 software. A bluetooth-connected customized MATLAB GUI interface was used to collect the exoskeleton knee joint angle and knee joint torque. The lifting task was repeated three times, each with a 3-min break between them. The EMG data series were processed in EMGworks Analysis and MATLAB. Next, the processed data from the MVN Analyze, Motive 3.0, and Matlab GUI were used to validate the simulation results against the experimental results with the exoskeleton.

Figure 7. Lifting experiment setup with the exoskeleton.