Abstract:
An exercise apparatus is disclosed. The exercise apparatus comprises: a three-dimensional coordinate system having an x-axis, a y-axis, and a z-axis; lower limb segments; a moving arm; a multi-axis force sensor; at least one footplate attached to said multi-axis force sensor for supporting a human foot or shoe; a compact real-time ankle motion measurement device, said motion measurement device attached between the leg and foot and capable of measuring the ankle joint motion in at least one degree of freedom; and a central processor, said central processor capable of calculating lower-limb joint moments in real time. The central processor calculates lower-limb joint moments with a modified inverse dynamic method, which does not involve computation of location of center of pressure.

Description:
BACKGROUND OF INVENTION 
     More than 21% of adult Americans (46 million) have some form of diagnosed arthritis, and the number is anticipated to rise to 67 million by 2030 (Hootman and Helmick, 2006). The most common type of arthritis, osteoarthritis (OA), affected more than 27 million Americans (Lawrence et al., 2008). Nearly one-half of all adults and two-thirds of those who are obese can now be expected to develop symptomatic knee OA at some point in their lives (Murphy et al., 2008). Knee OA is a major contributor to functional impairment and reduced independence in older adults (Peat et al., 2001). 
     There is a growing recognition of the role of local neuromechanical factors in influencing physical function and knee OA disease progression. One such factor is varus-valgus laxity and dynamic instability. However, few exercise treatment strategies to target this particular biomechanical and neuromuscular impairment exist. 
     Although it is well recognized that knee OA is associated with excessive knee adduction moment (KAM) in the frontal plane, there has been a lack of convenient and effective ways of training people with knee OA in better dealing with the frontal plane knee adduction moment loading. 
     In general, one or more components of lower-limb joint moments may need to be better controlled through feedback training for the purposes of post-injury rehabilitation and injury prevention. There is a general need for real-time feedback training of selected joint moment(s). 
     The multi-axis rehabilitation system described here addresses such a need and provides us a useful tool to better control a certain joint moment through real-time biofeedback of the targeted joint moment(s). 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  shows schematics of real-time knee moment computation and feedback training  FIG. 1(   a ) shows a multi-DOF goniometer used as the compact lower-limb motion measurement unit.  FIG. 1(   b ) shows an inertial measurement unit used as the compact lower-limb motion measurement unit. 
         FIG. 2  ( a ) is a sagittal plane schematic diagram of a right-side lower limb with the proposed robotic multi-axis elliptical trainer (ET).  FIG. 2(   b ) is a frontal plane free body diagram (FBD) of a foot (triangle), and a 6-axis force/torque (F/T) sensor, located underneath the foot; the top vertex of the triangle, P a , represents ankle JC.  FIG. 2(   c ) is a FBD of the shank (thick vertical line) with the ankle JC (lower end of the rod) and the knee JC (upper end of the rod). 
         FIG. 3  demonstrates 3D knee moment computation methods, among which  3 ( a ) shows well-known typical methods (Winter, 2005; Zatsiorsky, 2002), and  3 ( b ) shows the proposed method, which does not require COP computation thanks to the features of the proposed robotic multi-axis ET. 
         FIG. 4  is a schematic diagram of the robotic multi-axis ET with a lower limb. 
         FIG. 5  shows a six-DOF goniometer used to measure ankle dorsiflexion/planar-flexion and inversion/eversion angles and shank internal/external rotation angle accurately in real-time. 
         FIG. 6(   a ) shows knee adduction moment (KAM) estimated during normal stepping;  FIG. 6(   b ) shows knee adduction moment (KAM) estimated during adducted stepping. In the figures, KAM estimated with the 6-DOF goniometer (blue line) and with the Optotrak (red dashed line) matched each other closely. Peak KAM locations of the two methods in each cycle are almost identical. 
         FIG. 7  shows variables needed to be measured on-line for calculation of moments of lower limb joints 
     
    
    
     COMPONENT LIST 
     
         
         
           
               1 : Multi-axis force sensor 
               2 . Footplate 
               3 G: Multi-DOF goniometer 
               3 I: Inertial measurement unit (IMU) 
               4 : feedback device 
               5 : foot strap 
               6 : pivoting 
               7 : sliding 
               8 : motor 
               9 : position sensor 
               10 : potentiometer 
               11 : oblique ramp 
               12 : driving wheel 
               13 : long beam 
           
         
       
    
     DETAILED DESCRIPTION OF THE INVENTION 
     An apparatus and method of 3-dimensional (3-D) knee and ankle moments and forces estimation method are described, which can be used for real-time feedback training on controlling one or more of the 3-D knee and ankle joint moments and forces without requiring the computation of center of pressure (COP)—the location of the vertical ground reaction force (GRF) vector from a force platform ( FIG. 1 ). 
     The 3-D knee moments under real-time feedback control can be one or more of the moment components including knee flexion-extension, abduction-adduction, and internal-external rotation moments; and the 3-D knee forces can be one or more of the force components including lateral, anterior, and proximal translation forces. The 3-D ankle moments under real-time feedback control can be one or more of the moment components including dorsiflexion-plantar flexion, inversion-eversion and internal-external rotation moments; and the 3-D ankle forces can be one of more force components including lateral, anterior, and proximal translation forces. 
     For example, knee adduction moment (KAM), a widely used surrogate measure of medial knee compression (Foroughi et al., 2009; Komistek et al., 1999; Miyazaki et al., 2002; Pollo et al., 2002; Reeves and Bowling, 2011; Shelburne et al., 2008; Simic et al., 2011), is significantly correlated with the severity of knee osteoarthritis (OA) (Foroughi et al., 2009; Reeves and Bowling, 2011; Simic et al., 2011), a disease affecting more than 27 million Americans (Lawrence et al., 2008); and is the cause of 60˜80% of the compressive load at the medial tibiofemoral knee compartment (Foroughi et al., 2009). Especially, the largest peak of the KAM, which occurs during the load acceptance phase, is the strong predictor of medial compartment OA presence (Simic et al., 2011), radiographic disease severity (Sharma et al., 1998; Simic et al., 2011), rate of progression (Miyazaki et al., 2002), and the presence of OA symptoms (Thorp et al., 2007). 
     Thus, KAM estimated in real-time may enable the patients with knee OA to use it as a biofeedback in order to control the KAM during locomotion. Considering the wide usage of the ET, the real-time estimation of KAM during ET exercise may be beneficial to millions of patients with knee OA. Until now, there are no reports on real-time biofeedback of the KAM during the training but there are some reports on real-time feedback of dynamic knee frontal alignment (Barrios et al., 2010; Reeves and Bowling, 2011; Simic et al., 2011). 
     Since the emergence of the biomechanical and gait study, there has been extensive efforts on the estimation of 3D moments (as well as forces) with the advance of measurement technologies (Baker, 2006; Elftman, 1939; Koontz et al., 2006; Shiavi et al., 1987; Whittle, 1996; Winter, 2005; Zhang et al., 2003; Zhou and Hu, 2008), because the 3D moments of lower limb joints are important and provide valuable information in many disciplines including biomechanics, orthopedics, physical rehabilitation, and sports science. Examples are, to name just a few, assessing patients with a walking disorder (Whittle, 1996) based on the normal gait patterns (Elftman, 1939; Riley et al., 2001; Winter, 1984, 1980) and investigating association between KAM and biomechanical variables of patients with knee osteoarthritis (OA) (Foroughi et al., 2009; Hurwitz et al., 2002). 
     However, most current estimation methods need off-line post processing in order to obtain the moments and forces of lower limb joints. Moreover, there may be other potential difficulties in using readily available methods including a long duration of preparation of attaching markers on the parts of the lower limbs of interest for visual/non-visual motion tracking, occlusion at some points (phases) of gait with visual tracking system (e.g., Optotrak™), drift in the computation of velocity and position from accelerometer data, and distortion of sensed information of magnetic sensors due to the magnetic fields of other devices (e.g., electromagnetic motors and/or ferrite structures) (Zhou and Hu, 2008). 
     Besides, one of the differences of the proposed robotic multi-axis ET exercise and the ground walking (or the running) is that there is no relative motion between the foot and the footplate, because the bootstrap prevents motion of the foot relative to the footplate. Thus, the foot and the footplate always move together. In other words, the distance between the center of mass (COM) of the foot (P Ft ) and the center of the 6-axis Force/Torque (F/T) sensor (P F/T ), which is firmly mounted underneath the footplate, is time-invariant during the exercise. Moreover, due to the movement constraint imposed on the foot and the footplate, non-zero pure moment (also called couple), which is usually ignorable in normal walking and running in a gait lab setting (Cohen et al., 1980), is no longer negligible and may be applied to the footplate. Thus, the typical well established inverse dynamics methods (Winter, 2005; Zatsiorsky, 2002), which require the location of COP, is not applicable to the proposed robotic multi-axis ET exercise. Specifically, with a regular force plate in a gait lab setting, which gives same measured forces/moments as the 6-axis F/T sensor mounted on the proposed robotic multi-axis ET, the COP is computed under the assumption that there is no pure moment exerted on the force plate in the horizontal plane (i.e., zero T ax  and T ay ) (Cohen et al., 1980; Winter, 2005; Zatsiorsky, 2002). Therefore, it might not be applicable to use existing methods developed for the analysis of the over-ground walking and the running, because of the constrained foot movement. 
     Hence, in order to compute the 3D knee moments as well as other moments including ankle joint moments, a modified inverse dynamic method is described below, which does not involve the computation of location of COP and is tailored for the robotic multi-axis ET exercise. 
     1. Computation of COP 
     For the computation of the 3D knee moments and other lower limb moments, COP is generally needed in order to use well-established inverse dynamics methods based on Newton-Euler equation (Cohen et al., 1980; Winter, 2005; Zatsiorsky, 2002). A typical 3D knee moments computation method (Cohen et al., 1980; Winter, 2005; Zatsiorsky, 2002) may be summarized as follows ( FIG. 3(   a )): from the anthropometry data given in (Winter, 2005; Zatsiorsky, 2002) or some other reliable sources, kinematic and dynamic parameters including coordinates of COM, mass and inertia with respect to the COM of each segment (in this case, foot and shank) are obtained; COP is obtained using forces and moments measured with the force plate; translational and rotational accelerations are obtained from a motion tracking device (e.g. Optotrak); one can then compute moment (M a ε           ) and force (F a ε         ) vectors acting on the ankle using the inverse dynamics with GRF and COP; consequently, from M a  and F a  thus computed, knee moment vector (M k ε         ) can be obtained.
       FIG. 2  ( a ) is the sagittal plane schematic diagram of a right-side lower limb with the proposed robotic multi-axis elliptical trainer (ET). Positive x o  direction is the lateral direction; Positive y o  direction the anterior direction. A circle centered at P 4  represents circular disk of the ET. A thick line from P 2  to P 3  represents the rod where the footplate is mounted on. World frame denotes the frame which is a fixed Lab Ankle frame denotes a frame attached to foot and origin of frame is ankle joint center (JC), P a ; knee frame denotes a frame attached to shank and origin of the frame is knee JC, P k . P Sk  and P Ft  denote positions of center of mass (COM) of the shank and the foot, respectively.  FIG. 2(   b ) is the frontal plane free body diagram (FBD) of a foot (triangle), and a 6-axis force/torque (F/T) sensor, located underneath the foot; the top vertex of the triangle, P a , represents ankle JC. Positive y a  direction is the anterior direction; positive x a  direction the lateral direction. Since the footplate, which is located between the foot and the 6-axis F/T sensor, is highly stiff, it is omitted for clarity.  FIG. 2  ( c ) is FBD of the shank (thick vertical line) with the ankle JC (lower end of the rod) and the knee JC (upper end of the rod). From the forces thus measured with the 6-axis F/T sensor with the 3-dimensional (3D) translational and rotational accelerations of COM of the foot, forces and moments acting on the ankle JC are obtained. Consequently, using forces and moments acting on the ankle JC with the 3D translational and rotational accelerations of COM of the shank, 3D knee moments can be obtained. d 1  and h 1  denote the x a  and z a  direction coordinate differences between the F/T sensor (P F/T ), which is the point of action of measured forces/torques, and the COM of the foot (center of white circle, P Ft ), respectively; d 2  and h 2  denotes the x a  and z a  direction coordinate differences between the ankle JC(P a ) and the COM of the foot (P Ft ), respectively;  a a Ft     —     x  and  a a Ft     —     z  the x a  and z a  direction translational accelerations of the COM of the foot, respectively;  a {umlaut over (θ)} Ft     —     y  the y a  (frontal) direction rotational acceleration of COM of the foot;  a F F/T     —     x ,  a F F/T     —     z ,  k F F/T     —     x , and  k F  F/T     —     z  the x a , z a , x k  and z k  direction forces action on the bottom of the foot, respectively;  a M F/T     —     y  and  a M ay  the y a  (frontal) the direction moments measured with the 6-axis F/T sensor and computed at the ankle JC, respectively.  k M ky  and  k M ay  denote the y k  (frontal) direction moments at the knee JC and the ankle JC, respectively;  k a Sk     —     x  and  k a Sk     —     z  the x k  and z k  direction translational accelerations of COM of the shank, respectively;  a F ax ,  a F az ,  k F ax , and  k F az  denote the x a , z a , x k , and z k  direction forces acting on the ankle JC, respectively;  k F kx , and  k F kz  denote the x k  and z k  direction forces acting on the knee JC, respectively. Specifically, COP is required to compute the net moment acting on the foot. With the measured forces and moments, which correspond to F F/T     —     x , F F/T     —     z  and M F/T     —y    in  FIG. 2 , and COP, net resultant foot y direction moment (M y     —     net ), which is equated to the same direction foot inertial torque, can be expressed as follows (Winter, 2005) (Here, for simplicity, it is assumed that the footplate is in horizontal plane (i.e. β=0) and all equations are written with respect to world frame without loss of generality. Thus, left-superscripts, which represent frame of references, are omitted.):
 
 M   y     —     net   =F   F/T     —     z ( d   1   +x   cop )+ F   F/T     —     x   h   1   +M   ay   −F   az   d   2   +F   ax   h   2   (1)
 
where
 
 h   1   =z   F/T   −z   Ft   (2)
 
 h   2   =z   a   −z   Ft   (3)
 
 d   1   =x   F/T   −x   Ft   (4)
 
 d   2   =x   a   −x   Ft   (5)
 
P Ft =[x Ft  y Ft  z Ft ] T  and P a =[x a  y a  z a ] T  denote coordinates of the COM of the foot, and that of ankle, respectively; x cop  denotes the x direction coordinate difference between F/T sensor center and the COP, and can be computed as below
 
 x   cop   =x   cp   −x   F/T   (6)
 
where P F/T =[x F/T  y F/T  z F/T ] T  and P cp =[x cp  y cp  z cp ] T  denote coordinates of center of F/T sensor and that of COP.
 
     Thus, computing location of COP is an important part of the knee moment computation. In (Winter, 2005), of course, a COP computing method with regular force plate when there is no pure moment (also called couple) in horizontal direction (i.e., x a  and y a  direction in  FIG. 2(   b )) can be found. Again, outputs of regular support force plate are the same as those of the 6-axis F/T sensor mounted on the proposed robotic multi-axis ET. 
     In general, the y direction moment measured with the 6-axis F/T sensor M F/T     —     y  can be represented with forces F F/T     —     x  and F F/T     —     y  and pure moment T y  as follows (e.g., see AMTI or Bertec technical note or User manual):
 
 M   F/T     —     y   =F   F/T     —     x   z   off   −F   F/T     —     z   x   cop   +T   y   (7)
 
where z off  denotes the z direction distance between center of the 6-axis F/T sensor and the top surface of the sensor (plus footplate thickness). Note that (7) is a generalization of the relation of the COP with the forces and moments acting on a force plate in that, compared with the COP computation method given in Winter (2005), non-zero pure moment T y  is included. It is assumed in (Winter, 2005) that T y  is zero (i.e. T y =0), thus from (7), the x direction component of COP, x cop , was obtained as follows:
 
 x   cop =( F   F/T     —     x   z   off   −M   F/T     —     y )/ F   F/T     —     z   (8)
 
One disadvantage of using (8) is that if F F/T     —     z  is small (&lt;2% of body weight), x cop  is sensitive to the change of F F/T     —     z  and thus can be incorrect due to the measurement noise accompanied with F F/T     —     z  (Winter, 2005). Moreover, if T y  is not zero—which is the case of the proposed robotic multi-axis ET training—COP cannot be computed using (8). It should be noted that the aim of this specific research is not to find COP but to compute 3D knee and ankle moments.
 
     In order to remove the requirement of the COP location from net moment calculation in (1), advantages of the proposed robotic multi-axis ET are utilized. As was mentioned, the distance between the COM of foot (P Ft ) and the 6-axis F/T sensor (P F/T ) is fixed; so is the distance between the ankle joint (P a ) and the COM of foot (P Ft ). In other words, once a subject&#39;s anthropometric data is obtained, d 1 , d 2 , h 1  and h 2  in  FIG. 2(   b ) are known and time-invariant during the exercise. Thus, net y direction moment acting on the COM of the foot, M y     —     net  given in (1), can be rewritten as follows:
 
 M   y     —     net   =F   F/T     —     x   h   1   −F   F/T     —     z   d   1   +M   F/T     —     y   +F   ax   h   2   −F   az   d   2   +M   ay   (9)
 
Obviously, (9) does not require the COP location. Moreover, all the distances (d 1 , d 2 , h 1  and h 2 ) are time-invariant (i.e., constants) during the exercise. Thus, these distances can be computed off-line, whereas COP requires a real-time online computation for the real-time computation of the 3D knee and ankle moments (and forces). Accordingly, the problem of potential incorrectness of x cop  computation in (8) due to high sensitivity under the small F F/T     —     z  (Winter, 2005) is removed.
 
     Based on (9), the 3D knee moments and forces (and ankle moments and forces) computation can be simplified as shown in  FIG. 3(   b ). Compared with the typical methods (Cohen et al., 1980; Winter, 2005; Zatsiorsky, 2002), the COP computation is no longer needed. 
     In short, due to the design of the robotic multi-axis ET, 3D knee and ankle moments and forces can be computed without computing COP location. In the following, the details of the proposed 3D knee and ankle moments computation method are provided by exemplifying the computation of the knee abduction/adduction moment. In case, if the foot is not strapped to footplate, then the usual COP calculation can still be utilized for the computation of knee and ankle forces and moments. 
     2. Derivation of Dynamic Equations for Real-Time Computation of 3D Moments of Lower Limb Joints 
     In this section, dynamic equations for real-time computation of 3D moments and forces of lower limb joints are provided with a simple 6-DOF goniometer based kinematic measurement technique. 
     1) Computation of Kinematic Variables 
     For the computation of 3D knee and ankle moments, kinematic information, which include positions, velocities, and accelerations of the joints and those of COM of segments, and joint angles, angular velocities and accelerations, are required. Specifically, one can compute the 3D knee moments with footplate angle (β), ankle dorsiflexion/planar-flexion (θ 1 ), inversion/eversion (θ 2 ), and shank internal/external rotation (θ 3 ) angles, and shank length (L Sk ) between knee joint and ankle joint, length between knee joint and COM of shank (L Sk     —     com ), and footplate (or ankle) position, velocity, and acceleration. 
     A. Kinematic Variables from the Robotic Multi-Axis Elliptical Trainer 
       FIG. 4  is a schematic diagram of the robotic multi-axis ET with a lower limb. The distance from P 3  to P 4  is l 1 ; that of road from P 2  to P 3  is l 2 ; and that of from P F/T  to P 3  l 2FT ; that of from P pot , the location of potentiometer, to P 2  l ob . l 1 , l 2 , and l 2FT  are known constants in advance once the ET dimension is decide whereas l ob  is measured with the potentiometer mounted at fixed point P pot . Footplate angle (β) and footplate position can be computed using P 2  position ( FIG. 4 ) measured along the oblique ramp direction with a precision potentiometer (l ob ). Considering the serial linkages from P 2  to P 3  with length l 2  and from P 3  to P 4  with length l 1 , because the motion of P 2  is constrained by the oblique ramp, once the position of P 2  is measured then the two angles (φ 1 , φ 2 ) can be obtained using a simple inverse kinematics. With the l ob  measured and y and z direction coordinates of P pot  with respect to the world frame (i.e.,  o y pot  and  o z pot ), one can express they and z coordinates of P 2  with respect to the world frame ( o y 2  and  o z 2 ) as follows: 
                   {             y   2               o     =         -     l   ob       ⁢     cos   ⁡     (   α   )         +     y   pot               o                     z   2               o     =         -     l   ob       ⁢     sin   ⁡     (   α   )         +     z   pot               o                       (   10   )               
where α denotes the angle between the ground and the oblique ramp in the front part of the robotic multi-axis ET. Since α,  o y pot , and  o z pot  are constants, once l ob  is measured with the potentiometer located at P pot ,  o y 2  and  o z 2  can be obtained.
 
       o y 2  and  o z 2  can be also expressed using the two angles (φ 1 , φ 2 ) and lengths (l 1  and l 2 ) 
                   {             y   2               o     =         l   1     ⁢     cos   ⁡     (     ϕ   1     )         +       l   2     ⁢     cos   ⁡     (       ϕ   1     +     ϕ   2       )                         z   2               o     =         l   1     ⁢     sin   ⁡     (     ϕ   1     )         +       l   2     ⁢     sin   ⁡     (       ϕ   1     +     ϕ   2       )                           (   11   )               
From (11), φ 1 , φ 2  can be obtained as follows (Craig, 1989):
 
                   {               ϕ   1     =         tan     -   1       (       z   2               o       y   2               o       )     -       tan     -   1       (         l   2     ⁢     sin   ⁡     (     ϕ   2     )             l   1     +       l   2     ⁢     cos   ⁡     (     ϕ   2     )             )                     ϕ   2     =       tan     -   1       ⁡     (       sin   ⁡     (     ϕ   2     )         cos   ⁡     (     ϕ   2     )         )               ⁢     
     ⁢   where     ⁢     
             (   12   )               {             cos   ⁡     (     ϕ   2     )       =         y   2   2           o     +     z   2   2           o     -     (       l   1   2     +     l   2   2       )         2   ⁢     l   1     ⁢     l   2                       sin   ⁡     (     ϕ   2     )       =     {           -       1   -       cos   2     ⁡     (     ϕ   2     )                         if   ⁢           ⁢   footplate   ⁢           ⁢   is   ⁢           ⁢   in   ⁢           ⁢   upper               position   ⁢           ⁢     (       180   ⁢   °     &lt;     ϕ   2     ≤     360   ⁢   °       )                         1   -       cos   2     ⁡     (     ϕ   2     )                       if   ⁢           ⁢   footplate   ⁢           ⁢   is   ⁢           ⁢   in   ⁢           ⁢   lower               position   ⁢           ⁢     (       0   ⁢   °     &lt;     ϕ   2     ≤     180   ⁢   °       )                                     (   13   )               
Thus, φ 1  and φ 2  can be obtained with  o y 2  and  o z 2  directly computed from l ob  measured with a potentiometer located at P 2 . Although there are two inverse kinematics solutions for a given P 2  position, however, a correct one can be easily selected with initial condition such as starting position.
 
With φ 1 , φ 2  thus measured, footplate angle (β) can be easily derived as below
 
β=φ 1 +φ 2   (14)
 
     From simple forward kinematics, P F/T  the 6-axis F/T sensory and z coordinates with respect to world frame can be computed with φ 1  and φ 2  as below 
                   {             y     F   /   T                 o     =         l   1     ⁢     cos   ⁡     (     ϕ   1     )         +       l     2   ⁢   FT       ⁢     cos   ⁡     (       ϕ   1     +     ϕ   2       )                         z     F   /   T                 o     =         l   1     ⁢     sin   ⁡     (     ϕ   1     )         +       l     2   ⁢   FT       ⁢     sin   ⁡     (       ϕ   1     +     ϕ   2       )                           (   15   )               
If the footplate and the 6-axis F/T sensor below the plate is rotating or sliding, angle of rotation or sliding distance can be incorporated into the computation of location of the F/T sensor (P F/T ) in (15).
 
     B. Kinematic Variables from the 6-DOF Goniometer 
       FIG. 5  shows a six-DOF goniometer used to measure ankle dorsiflexion/planar-flexion, inversion/eversion, and internal/external rotation angles accurately in real-time. The proximal end is firmly strapped to the bony frontal-medial surface of shank and the distal end is attached to the footplate, which has no motion relative to foot. The ankle joint center (JC), midpoint between the medial and lateral malleoli, was aligned vertically with the F/T sensor center. By using the 6-DOF goniometer, the ankle dorsiflexion/planar flexion angle (θ 1 ), inversion/eversion angle (θ 2 ), and internal/external rotation angle (θ 3 ) can be measured with proper calibration. Rotation angles are defined by following the sequence of x-y-z Euler angle. In order to measure the ankle dorsiflexion/planar-flexion angle (θ 1 ), inversion/eversion angle (θ 2 ), and internal/external rotation angle (θ 3 ) in real-time, a low-cost 6 degree-of-freedom (DOF) goniometer ( FIG. 5 ), which can be conveniently connected to any analog to digital (A/D) data acquisition devices together with other analog signals including the 6-axis forces/torques and locomotion variables, can be utilized without time-consuming preparation (Shiavi et al., 1987; Zhang et al., 2003). The only preparation required is to calibrate the initial angles with a simple calibration plate, on which all potentiometer angles of the 6-DOF goniometer are at known positions (Shiavi et al., 1987; Zhang et al., 2003). Note that the zero angles of 6-DOF goniometer, calibrated with the calibration plate, are not the same as zero anatomical angles of dorsiflexion/planar-flexion, inversion/eversion and internal/external rotation. 
     In order to calculate the ankle dorsiflexion/planar-flexion, inversion/eversion, and internal/external rotation angle using the 6-DOF goniometer, before start training on ET, subject is asked to stand straight by trying to align knee joint with laser pointer lights coming from lateral (and medial) side and from front. The laser pointer at the lateral side of footplate is aligned with center of F/T sensor (i.e., y coordinate of the laser pointer with respect to ankle frame is the same as that of 6-axis F/T sensor) and can be tilted vertically up and down in x a -z a  plane with respect to ankle frame. The laser pointer at the front size of footplate is located in front of 2 nd  toe (i.e., x coordinate of center of ankle with respect to ankle frame and that of the laser pointer is identical) and also can be tilted vertically up and down in y a -z a  plane with respect to ankle frame. By this way, one can find a posture that can be regarded as ‘zero degree posture’ (i.e., 0° dorsiflexion/planar-flexion, 0° inversion/eversion, 0° internal/external rotation). 
     Once the knee joint alignment is confirmed by checking the laser lights, the 6 angles of goniometers and the rotation matrix from the last (6 th ) frame (attached to the bony surface of shank) to the base frame of the 6-axis goniometer (attached to the footplate) at that time ( lg   bg R i ) are saved. Thus, one can compute rotation matrix from knee frame to ankle frame as follows (Shiavi et al., 1987; Zhang et al., 2003):
 
 k   a   R=   lg   bg   R ( lg   bg   R   i ) −1   (16)
 
where  bg   lg R denotes rotation matrix from base frame to last frame of the 6-axis goniometer. From ankle dorsiflexion/planar-flexion (θ 1 ) and inversion/eversion (θ 2 ) angles and shank internal/external rotation angle (θ 3 ),  k   a R can be also obtained as follows:
 
                                 k   a     ⁢   R     =       ⁢     [           r   11           k   a           r   12           k   a           r   13           k   a               r   21           k   a           r   22           k   a           r   23           k   a               r   31           k   a           r   32           k   a           r   33           k   a           ]                 =       ⁢     [           C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3               -   C     ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3             S   ⁢           ⁢     θ   2                   S   ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +     C   ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   3                   -   S     ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +     C   ⁢           ⁢     θ   1     ⁢   C   ⁢           ⁢     θ   3                 -   S     ⁢           ⁢     θ   1     ⁢   C   ⁢           ⁢     θ   2                     -   C     ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +     S   ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   3                 C   ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +     S   ⁢           ⁢     θ   1     ⁢   C   ⁢           ⁢     θ   3               C   ⁢           ⁢     θ   1     ⁢   C   ⁢           ⁢     θ   2             ]                   (   17   )               
where  k   a r ij  denotes i,j th  element of  k   a R obtained from (16). Note that, although the 6-DOF goniometer has a one translational DOF and a potentiometer for measuring the translation, this information is not required for computing  k   a R. Clearly, from (16) and (17), one can obtain the three angles (θ 1 , θ 2 , θ 3 ) as follows:
 
                       θ   1     =       tan     -   1       ⁡     (       -     r   23           k   a         r   33           k   a       )         ,     
     ⁢       θ   2     =       tan     -   1       ⁡     (       r   13           k   a           -   S     ⁢           ⁢     θ   1     ⁢     r   23           k   a       +     C   ⁢           ⁢     θ   1     ⁢     r   33           k   a           )         ,     
     ⁢       θ   3     =       tan     -   1       ⁡     (       -     r   12           k   a         r   11           k   a       )                 (   18   )               
From (18), one can select a proper solution because the dorsiflexion/planar-flexion, inversion/eversion and internal/external rotation angles must be between ±90° during normal exercise.
 
     By using a digitizer (e.g., MicroScribe or probe of Optotrak) or referring anthropometric data (Winter, 2005; Zatsiorsky, 2002), the length of shank, L Sk , from ankle to knee joint can be measured. Moreover, with L Sk  thus measured, the length between knee JC and COM of shank, L Sk     —     com  can be obtained by referring anthropometric data. Note that the two lengths (L Sk , L Sk     —     com ) are constant values for each subject, thus only one time off-line measurement is required for each subject unless the subject is growing during the training period. Depending on the joint and forces/moments of interest, not all of the angles are required for the computation as shown in  FIG. 7 . 
     2) Derivation of Equations of Motion for Computation of 3D Moments of Lower Limb Joints 
     For the derivation of knee moment, available variables, which can be known in advance or measured/estimated, are listed below. 
     1) Ankle position with respect to the world frame
 
 o   P   a =[ o   x   a   o   y   a   o   z   a ] T   (19)
 
with respect to the ankle frame
 
 a   P   a =[0 0 0] T   (20)
 
and with respect to the knee frame
 
 k   P   a =[0 0− L   Sk ] T   (21)
 
2) Knee position with respect to the knee frame
 
 k   P   k =[ k   x   k   k   y   k   k   z   k ] T =[0 0 0] T .  (22)
 
3) Shank COM position with respect to the knee frame
 
 k   P   Sk =[ k   x   Sk   k   y   Sk   k   z   Sk ] T =[0 0− L   Sk     —     com ] T   (23)
 
4) COM of foot position with respect to the ankle frame
 
 a   P   Ft =[ a   x   Ft   a   y   Ft   a   z   Ft ] T   (24)
 
which is a constant vector (i.e.,  a x Ft ,  a y Ft , and  a z Ft  are constants).
 
5) Center of 6-axis F/T sensor position with respect to the ankle frame
 
 a   P   F/T =[ a   x   F/T   a   y   F/T   a   z   F/T ] T   (25)
 
which is a constant vector (i.e.,  a x F/T ,  a y F/T  and  a z F/T  are constants).
 
6) 6-axis measurement of forces and torques exerted on foot by robotic multi-axis ET with respect to the ankle frame
 
 a   F   F/T =[ a   F   F/T     —     x   a   F   F/T     —     y   a   F   F/T     —     z ] T   (26)
 
and
 
 a   M   F/T =[ a   M   F/T     —     x   a   M   F/T     —     y   a   M   F/T     —     z ] T   (27)
 
7) Angular acceleration of foot with respect to the world frame and the ankle frame
 
 a {umlaut over (θ)} Ft = o {umlaut over (θ)} Ft =[{umlaut over (β)}0 0] T .  (28)
 
8) Angular velocity and acceleration of shank with respect to the knee frame
 
                       θ   .     Sk               k     =     [               (         θ   .     1     +     β   .       )     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +         θ   .     2     ⁢   S   ⁢           ⁢     θ   3                       -     (         θ   .     1     +     β   .       )       ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +         θ   .     2     ⁢   C   ⁢           ⁢     θ   3                       (         θ   .     1     +   β     )     ⁢   S   ⁢           ⁢     θ   2       +       θ   .     3             ]             (   29   )                   θ   ¨     Sk               k     =     [               (         θ   ¨     1     +     β   ¨       )     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +         θ   ¨     2     ⁢   S   ⁢           ⁢     θ   3       -       (         θ   .     1     +     β   .       )     ⁢     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3         )       +         θ   .     2     ⁢       θ   .     3     ⁢   C   ⁢           ⁢     θ   3                       -     (         θ   ¨     1     +     β   ¨       )       ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +         θ   ¨     2     ⁢   C   ⁢           ⁢     θ   3       +       (         θ   .     1     +     β   .       )     ⁢     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       -         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3         )       -         θ   .     2     ⁢       θ   .     3     ⁢   S   ⁢           ⁢     θ   3                       (         θ   ¨     1     +     β   ¨       )     ⁢   S   ⁢           ⁢     θ   2       +       θ   ¨     3     +       (         θ   .     1     +     β   .       )     ⁢       θ   .     2     ⁢   C   ⁢           ⁢     θ   2               ]             (   30   )               
9) Rotation matrix from the ankle frame to the world frame
 
                           a   o     ⁢   R     =     [         1       0       0           0         C   ⁢           ⁢   β             -   S     ⁢           ⁢   β             0         S   ⁢           ⁢   β           C   ⁢           ⁢   β           ]             (   31   )               
10) Rotation matrix from the world frame to the ankle frame
 
                           o   a     ⁢   R     =     [         1       0       0           0         C   ⁢           ⁢   β           S   ⁢           ⁢   β             0           -   S     ⁢           ⁢   β           C   ⁢           ⁢   β           ]             (   32   )               
11) Rotation matrix from the ankle frame to the knee frame
 
                           a   k     ⁢   R     =     [           C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3               S   ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +     C   ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   3                   -   C     ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +     S   ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   3                     -   C     ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3                 -   S     ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +     C   ⁢           ⁢     θ   1     ⁢   C   ⁢           ⁢     θ   3                 C   ⁢           ⁢     θ   1     ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +     S   ⁢           ⁢     θ   1     ⁢   C   ⁢           ⁢     θ   3                   S   ⁢           ⁢     θ   2               -   S     ⁢           ⁢     θ   1     ⁢   C   ⁢           ⁢     θ   2             C   ⁢           ⁢     θ   1     ⁢   C   ⁢           ⁢     θ   2             ]             (   33   )               
12) Rotation matrix from the knee frame to the world frame
 
                           k   o     ⁢   R     =     [           C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3               -   C     ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3             S   ⁢           ⁢     θ   2                     S   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +       C   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   3                   -     S   ⁡     (     β   +     θ   1       )         ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +       C   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   3                 -     S   ⁡     (     β   +     θ   1       )         ⁢   C   ⁢           ⁢     θ   2                     -     C   ⁡     (     β   +     θ   1       )         ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +       S   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   3                   C   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +       S   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   3                 C   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   2             ]             (   34   )               
13) Gravitational acceleration vector with respect to the world frame
 
 o   g=[ 0 0 −g]   T   (35)
 
with respect to the ankle frame
 
 a   g=   o   a   R   o   g=[ 0 −gSβ−gCβ]   T   (36)
 
with respect to the knee frame
 
                           k     ⁢   g     =         R           a       a   k     ⁢   g     =       [           -     g   [         -     C   ⁡     (     β   +     θ   1       )         ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +       S   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   3                       -     g   ⁡     [         C   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +       S   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   3         ]                     -     gC   ⁡     (     β   +     θ   1       )         ⁢   C   ⁢           ⁢     θ   2             ]     .               (   37   )               
From (19), (23) and (34), COM of shank position with respect to the world frame can be computed as follows:
 
                     P   Sk               o     =       [             x   a               o     -       L   Sk_com     ⁢   S   ⁢           ⁢     θ   2                     y   a               o     +       L   Sk_com     ⁢     S   ⁡     (       θ   1     +   β     )       ⁢   C   ⁢           ⁢     θ   2                     z   a               o     -       L   Sk_com     ⁢     C   ⁡     (       θ   1     +   β     )       ⁢   C   ⁢           ⁢     θ   2               ]     .             (   38   )               
By taking 2 nd  order time derivative of both side (38) and multiplying  o   k R, the transpose of  k   o R given in (34), on both side, acceleration of COM of shank,  k a Sk , can be obtained as follows:
 
                     a   Sk               k     =     [             C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3   o     ⁢       x   ¨     a       +         [         S   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +       C   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   3         ]     o     ⁢       y   ¨     a                     +     [         -     C   ⁡     (     β   +     θ   1       )         ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +       S   ⁡     (     β   +     θ   1       )       ⁢   S     ⁢           -     θ   3       ]       ⁢       z   ¨     a               o                 +       L   Sk_com     [       C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢       θ   3     ⁡     (         θ   ¨     1     +     β   ¨       )         -     C   ⁢           ⁢     θ   3     ⁢       θ   ¨     2                             +   0.5     ⁢           ⁢     S   ⁡     (     2   ⁢           ⁢     θ   2       )       ⁢   C   ⁢           ⁢         θ   3     ⁡     (         θ   .     1     +     β   .       )       2       -     2   ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢       θ   3     ⁡     (         θ   .     1     +     β   .       )       ⁢       θ   .     2         ]                                   -   C     ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3     ⁢       x   ¨     a               o       +       [         C   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   3       -       S   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3         ]     ⁢       y   ¨     a               o                     +     [         S   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   3       +       C   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3         ]       ⁢       z   ¨     a               o                   +     L   Sk_com       ⁢     {         +   C     ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢       θ   3     ⁡     (         θ   ¨     1     +     β   ¨       )         +     S   ⁢           ⁢     θ   3     ⁢       θ   ¨     2                             +   0.5     ⁢     S   ⁡     (     2   ⁢           ⁢     θ   2       )       ⁢   S   ⁢           ⁢         θ   3     ⁡     (         θ   .     1     +     β   .       )       2       -     2   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢       θ   3     ⁡     (         θ   .     1     +     β   .       )       ⁢       θ   .     2         }                                     x   ¨     a               o     ⁢   S   ⁢           ⁢     θ   2       -         y   ¨     a               o     ⁢     S   ⁡     (     β   -     θ   1       )       ⁢   C   ⁢           ⁢     θ   2       +         z   ¨     a               o     ⁢     C   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   2                   +       L   Sk_com     ⁡     [         θ   .     2   2     +         (         θ   .     1     +     β   .       )     2     ⁢           ⁢     C   2     ⁢     θ   2         ]               ]             (   39   )               
Similarly, from (19), (24), and (31), COM of foot position with respect to the world frame can be obtained as follows:
 
                     P   Ft               o     =         P   a               o     +       R       a       a   o     ⁢     P   Ft         =     [             x   Ft               a     +     x   a               o                     y   Ft               a     ⁢   C   ⁢           ⁢   β     -       z   Ft               a     ⁢   S   ⁢           ⁢   β     +     y   a               o                     y   Ft               a     ⁢   S   ⁢           ⁢   β     +       z   Ft               a     ⁢   C   ⁢           ⁢   β     +     z   a               o             ]               (   40   )               
Thus, from (40), one can get acceleration of COM of foot with respect to the ankle frame.
 
                     a   Ft               a     =     [             x   ¨     a               o                     y   ¨     a               o     ⁢   C   ⁢           ⁢   β     +         z   ¨     a               o     ⁢   S   ⁢           ⁢   β     -         β   .     2     ⁢           ⁢     y   Ft               a       -         β   ¨     a     ⁢     z   Ft                       -       y   ¨     a               o       ⁢   S   ⁢           ⁢   β     +         z   ¨     a               o     ⁢   C   ⁢           ⁢   β     +         β   ¨     a     ⁢     y   Ft       -         β   .     2     ⁢           ⁢     z   Ft               a               ]             (   41   )               
Because  a P Ft =[ a x Ft   a y Ft   a z Ft ] T  is a constant vector,  a {dot over (x)} Ft ,  a {dot over (y)} Ft ,  a ż Ft  and  a {umlaut over (x)} Ft ,  a ÿ Ft ,  a {umlaut over (z)} Ft  are not shown in (41). Note that the aim of this derivation is to compute the knee and ankle moments and forces.
 
     From  FIG. 2(   c ), 6 DOF equation of motion of shank with respect to knee frame can be written as follows:
 
 m   Sk   k   a   Sk =− k   F   a + k   F   k   +m   Sk   k   g   (42)
 
and
 
 k   I   Sk   k {umlaut over (θ)} Sk + k {dot over (θ)} Sk × k   I   Sk   k {dot over (θ)} Sk =( k   P   a − k   P   Sh )×(− k   F   a )+( k   P   k − k   P   Sh )× k   F   k − k   M   a + k   M   k   (43)
 
(43) is written with respect to the knee frame, however, in (43),  k I Sk  denotes the inertia matrix of shank about its COM, and the moment arms ( k P a − k P Sh ), and ( k P k − k P Sh ) are the vectors from COM of shank to ankle, and to knee, respectively. Thus, (43) can be regarded as the equation of rotational motion of shank with respect to the COM of shank frame—a frame which has the same orientation of the knee frame but its origin is located at COM of shank.
 
     Similarly, from  FIG. 2(   b ), 6 DOF equation of motion of foot with respect to the ankle frame can be obtained.
 
 m   Ft   a   a   Ft = a   F   F/T + a   F   a   +m   Ft   a   g   (44)
 
 a   I   Ft   a {umlaut over (θ)} Ft =( a   P   F/T − a   P   Ft )× a   F   F/T +( a   P   a − a   P   Ft )× a   F   a + a   M   F/T + a   M   a   (45)
 
     Similar to (43), (45) is written with respect to ankle frame, however, in (45),  a I Ft  denotes the inertia matrix of foot about its COM, and the moment arms ( a P F/T − a P Ft ), and ( a P a − a P Ft ) are the vectors from COM of foot to center of 6-axis F/T sensor, and to ankle, respectively. Thus, (45) can be regarded as the equation of rotational motion with respect to the COM of foot frame—a frame which has the same orientation of ankle frame but its origin is located at COM of foot. Further, because  a {dot over (θ)} Ft ==[{dot over (β)}0 0] T ,  a {dot over (θ)} Ft × a I Ft   a {dot over (θ)} Ft  is always a zero vector (we assumed that foot inertia matrix is a diagonal matrix). Thus it is omitted from left-hand-side of (45). Solving (43) for  k M k , we have
 
 k   M   k = k   I   Sk   k {umlaut over (θ)} Sk + k {dot over (θ)} Sk × k   I   Sk   k {dot over (θ)} Sk −( k   P   a − k   P   Sh )×(− k   F   a )−( k   P   k − k   P   Sh )× k   F   k + k   M   a   (46)
 
Substituting (20), (22), (23), and (29) into (46) and rearranging it yields
 
                     [           M   kx               k               M   ky               k               M   kz               k           ]     =       [             I   Sk_x               k     ⁢       θ   ¨     Sk_x               k                   I   Sk_y               k     ⁢       θ   ¨     Sk_y               k                   I   Sk_z               k     ⁢       θ   ¨     Sk_z               k             ]     +     
     ⁢           ⁢     [             (       I   Sk_z               k     -     I   Sk_y               k       )     ⁢     {         (         θ   .     1     +     β   .       )     ⁡     [         -     (         θ   .     1     +     β   .       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       -         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3         )       ]       +         θ   .     2     ⁢       θ   .     3     ⁢   C   ⁢           ⁢     θ   3         }                         (       I   Sk_x               k     -     I   Sk_z               k       )       ⁢     {         (         θ   .     1     +     β   .       )     ⁡     [         (         θ   .     1     +     β   .       )     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3         )       ]       +         θ   .     2     ⁢       θ   .     3     ⁢   S   ⁢           ⁢     θ   3         }                   (       I   Sk_y               k     -     I   Sk_x               k       )     ⁡     [           (         θ   .     1     +     β   .       )     ⁡     [         -     (         θ   .     1     +     β   .       )       ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3     ⁢   C   ⁢           ⁢     θ   3       +         θ   .     2     ⁢     C   ⁡     (     2   ⁢           ⁢     θ   3       )           ]       ⁢   C   ⁢           ⁢     θ   2       +         θ   .     2   2     ⁢   S   ⁢           ⁢     θ   3     ⁢   C   ⁢           ⁢     θ   3         ]             ]     +     
     ⁢           ⁢     [               (       L   Sk     -     L   Sk_com       )     k     ⁢     F   ay                   -       (       L   Sk     -     L   Sk_com       )     k       ⁢     F   ax               0         ]     +     [             L   Sk_com     ⁢           ⁢     F   ky               k                   -     L   Sk_com       ⁢           ⁢     F   kx               k               0         ]     +     [           M   ax               k               M   ay               k               M   az               k           ]               (   47   )               
From (47), one can have the following expression of  k M k 
 
                   {             M   kx               k     =         I   Sk_x               k     ⁢       θ   ¨     Sk_x               k       +     M   ax               k                     +     (       I   Sk_z               k     -     I   Sk_y               k       )       ⁢     {         (         θ   .     1     +     β   .       )     ⁡     [         -     (         θ   .     1     +     β   .       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       -         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3         )       ]       +         θ   .     2     ⁢       θ   .     3     ⁢   C   ⁢           ⁢     θ   3         }                     +       (       L   Sk     -     L   Sk_com       )     k       ⁢     F   ay       +       L   Sk_com     ⁢           ⁢     F   ky               k                       M   ky               k     =         I   Sk_y               k     ⁢       θ   ¨     Sk_y               k       +     M   ay               k         ⁢                         +     (       I   Sk_x               k     -     I   Sk_z               k       )       ⁢     {         (         θ   .     1     +     β   .       )     ⁡     [         (         θ   .     1     +     β   .       )     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3         )       ]       +         θ   .     2     ⁢       θ   .     3     ⁢   S   ⁢           ⁢     θ   3         }                     -       (       L   Sk     -     L   Sk_com       )     k       ⁢     F   ax       -       L   Sk_com     ⁢           ⁢     F   kx               k                     M   kz               k     =         I   Sk_z               k     ⁢           ⁢       θ   ¨     Sk_z               k       +     M   az               k                   +       (       I   Sk_y               k     -     I   Sk_x               k       )     ⁡     [           (         θ   .     1     +     β   .       )     ⁡     [         -     (         θ   .     1     +     β   .       )       ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3     ⁢   C   ⁢           ⁢     θ   3       +         θ   .     2     ⁢     C   ⁡     (     2   ⁢           ⁢     θ   3       )           ]       ⁢   C   ⁢           ⁢     θ   2       +         θ   .     2   2     ⁢   S   ⁢           ⁢     θ   3     ⁢   C   ⁢           ⁢     θ   3         ]                       (   48   )               
Substituting (30) into (48) yields
 
 k   M   kx = k   I   Sk     —     x [({umlaut over (θ)} 1 +{umlaut over (β)}) Cθ   2   Cθ   3 +{umlaut over (θ)} 2   Sθ   3 −({dot over (θ)} 1 +{dot over (β)})({dot over (θ)} 2   Sθ   2   Cθ   3 +{dot over (θ)} 3   Cθ   2   Sθ   3 )+{dot over (θ)} 2 {dot over (θ)} 3   Sθ   3 ]
 
+( k   I   Sk     —     z − k   I   Sk     —     y ){({dot over (θ)} 1 +{dot over (β)})[−({dot over (θ)} 1 +{dot over (β)}) Sθ   2   Cθ   2   Sθ   3 +({dot over (θ)} 2   Sθ   2   Cθ   3 −{dot over (θ)} 3   Cθ   2   Sθ   3 )]+{dot over (θ)} 2 {dot over (θ)} 3   Cθ   3 }
 
+( L   Sk   −L   Sk     —     com ) k   F   ay   +L   Sk     —     com   k   F   ky + k   M   ax  
 
 k   M   ky = k   I   Sk     —     y [−({umlaut over (θ)} 1 +{umlaut over (β)}) Cθ   2   Sθ   3 +{umlaut over (θ)} 2   Cθ   3 +({dot over (θ)} 1 +{dot over (β)})({dot over (θ)} 2   Sθ   2   Sθ   3 −{dot over (θ)} 3   Cθ   2   Cθ   3 )−{dot over (θ)} 2 {dot over (θ)} 3   Sθ   3 ]
 
+( k   I   Sk     —     x − k   I   Sk     —     z ){({dot over (θ)} 1 +{dot over (β)})[({dot over (θ)} 1 +{dot over (β)}) Sθ   2   Cθ   2   Cθ   3 +({dot over (θ)} 2   Sθ   2   Sθ   3 +{dot over (θ)} 3   Cθ   2   Cθ   3 )]+{dot over (θ)} 2 {dot over (θ)} 3   Sθ   3 }
 
−( L   Sk   −L   Sk     —     com ) k   F   ax   −L   Sk     —     com   k   F   kx + k   M   ay  
 
 k   M   kz = k   I   Sk     —     z [({umlaut over (θ)} 1 +{umlaut over (β)}) Sθ   2 +{umlaut over (θ)} 3 +({dot over (θ)} 1 +{dot over (β)}){dot over (θ)} 2   Cθ   2 ]+ k   M   az  
 
+( k   I   Sk     —     y − k   I   Sk     —     x )[({dot over (θ)} 1 +{dot over (β)})[−({dot over (θ)} 1 +{dot over (β)}) Cθ   2   Sθ   3   Cθ   3 +{dot over (θ)} 2   C (2θ 3 )] Cθ   2 +{dot over (θ)} 2   2   Sθ   3   Cθ   3 ]  (49)
 
Substituting (37) into (42) and solving for  k F k , we have
 
                     F   k               k     =     [               m   Sk     ⁢     {       a   Sk_x               k     +     g   ⁡     [         -     C   ⁡     (     β   -     θ   1       )         ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +       S   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   3         ]         }       +     F   ax               k                     m   Sk     ⁢     {       a   Sk_y               k     +     g   ⁡     [         C   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +       S   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   3         ]         }       +     F   ay               k                     m   Sk     ⁢     {       a   Sk_z               k     +       gC   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   2         }       +     F   az               k             ]             (   50   )               
Substituting (50) into (49), one can rewrite (49) as follows:
 
                   {             M   kx               k     =       I   Sk_x               k     ⁡     [         (         θ   ¨     1     +     β   ¨       )     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +         θ   ¨     2     ⁢   S   ⁢           ⁢     θ   3       -       (         θ   .     1     +     β   .       )     ⁢     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3         )       +         θ   .     2     ⁢       θ   .     3     ⁢   S   ⁢           ⁢     θ   3         ]                     +     (       I   Sk_z               k     -     I   Sk_y               k       )       ⁢     {         (         θ   .     1     +     β   .       )     ⁡     [         -     (         θ   .     1     +     β   .       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       -         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3         )       ]       +         θ   .     2     ⁢       θ   .     3     ⁢   C   ⁢           ⁢     θ   3         }                     +     L   Sk_com       ⁢     n   Sk     ⁢     {       a   Sk_y               k     +     g   ⁡     [         C   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +       S   ⁡     (     β   +     θ   1       )       ⁢   C   ⁢           ⁢     θ   3         ]         }       +       L   Sk     ⁢           ⁢     F   ay               k       +     M   ax               k                                   M   ky               k     =       I   Sk_y               k     ⁡     [         -     (         θ   ¨     1     +     β   ¨       )       ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +         θ   ¨     2     ⁢   C   ⁢           ⁢     θ   3       +       (         θ   .     1     +     β   .       )     ⁢     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       -         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3         )       -         θ   .     2     ⁢       θ   .     3     ⁢   S   ⁢           ⁢     θ   3         ]                     +     (       I   Sk_x               k     -     I   Sk_z               k       )       ⁢     {         (         θ   .     1     +     β   .       )     ⁡     [         (         θ   .     1     +     β   .       )     ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +     (           θ   .     2     ⁢   S   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3       +         θ   .     3     ⁢   C   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3         )       ]       +         θ   .     2     ⁢       θ   .     3     ⁢   S   ⁢           ⁢     θ   3         }                     -     L   Sk_com       ⁢     m   Sk     ⁢     {       a   Sk_x               k     +     g   ⁡     [         -     C   ⁡     (     β   +     θ   1       )         ⁢   S   ⁢           ⁢     θ   2     ⁢   C   ⁢           ⁢     θ   3       +       S   ⁡     (     β   +     θ   1       )       ⁢   S   ⁢           ⁢     θ   3         ]         }       -       L   Sk     ⁢           ⁢     F   ax               k       +     M   ay               k                                   M   kz               k     =       I   Sk_z               k     ⁡     [         (         θ   ¨     1     +     β   ¨       )     ⁢   S   ⁢           ⁢     θ   2       +       θ   ¨     3     +       (         θ   .     1     +     β   .       )     ⁢       θ   .     2     ⁢   C   ⁢           ⁢     θ   2         ]                   +       (       I   Sk_y               k     -     I   Sk_x               k       )     ⁡     [           (         θ   .     1     +     β   .       )     ⁡     [         -     (         θ   .     1     +     β   .       )       ⁢   C   ⁢           ⁢     θ   2     ⁢   S   ⁢           ⁢     θ   3     ⁢   C   ⁢           ⁢     θ   3       +         θ   .     2     ⁢     C   ⁡     (     2   ⁢           ⁢     θ   3       )           ]       ⁢   C   ⁢           ⁢     θ   2       +         θ   .     2   2     ⁢   S   ⁢           ⁢     θ   3     ⁢   C   ⁢           ⁢     θ   3         ]                   +     M   az               k                     (   51   )               
Solving (44) for  a F a  yields
 
 a   F   a   =m   Ft   a   a   Ft − a   F   F/T   −m   Ft   a   g   (52)
 
Substituting (36) and (41) into (52), we have
 
                     [           F   ax               a               F   ay               a               F   az               a           ]     =     [               m   Ft     ⁢       x   ¨     a               o       -     F     F   ⁢     /     ⁢   T_x                 a                     m   Ft     ⁡     (           y   ¨     a               o     ⁢   C   ⁢           ⁢   β     +         z   ¨     a               o     ⁢   S   ⁢           ⁢   β     -         β   .     2     ⁢           ⁢     y   Ft               a       -         β   ¨     a     ⁢     z   Ft         )       -     F     F   ⁢     /     ⁢   T_y                 a     +       m   Ft     ⁢   gS   ⁢           ⁢   β                     m   Ft     ⁡     (         -       y   ¨     a               o       ⁢   S   ⁢           ⁢   β     +         z   ¨     a               o     ⁢   C   ⁢           ⁢   β     +         β   ¨     a     ⁢     y   Ft       -         β   .     2     ⁢           ⁢     z   Ft               a         )       -     F     F   ⁢     /     ⁢   T_z                 a     +       m   Ft     ⁢   gC   ⁢           ⁢   β             ]             (   53   )               
Solving (45) for  a M a  yields
 
 a   M   a = a   I   Ft   a {umlaut over (θ)} Ft −( a   P   F/T − a   P   Ft )× a   F   F/T −( a   P   a − a   P   Ft )× a   F   a − a   M   F/T   (54)
 
Substituting (20), (24)-(28), and (53) into (54) and rearranging it yield
 
 a   M   ax = a   F   F/T     —     y   a   z   F/T − a   F   F/T     —     z   a   y   F/T − a   M   F/T     —     x  
 
+ m   Ft [−( a   y   Ft   Sβ+   a   z   Ft   C β) o   ÿ   a +( a   y   Ft   Cβ−   a   z   Ft   S β)( o   {umlaut over (z)}   a   +g )]
 
+[ a   I   ft     —     x   +m   Ft ( a   y   Ft   2 + a   z   Ft   2 )]{umlaut over (β)}  (55)+
 
 a   M   ay =− a   F   F/T     —     x   a   z   F/T + a   F   F/T     —     z   a   x   F/T − a   M   F/T     —     y  
 
+ m   Ft { a   z   Ft   o   {umlaut over (x)}   a + a   x   Ft   Sβ   o   ÿ   a − a   x   Ft   C β( o   {umlaut over (z)}   a   +g )− a   x   Ft   a   y   Ft {umlaut over (β)}+ a   x   Ft   a   z   Ft {dot over (β)} 2 }  (56)
 
and
 
 a   M   az = a   F   F/T     —     x   a   y   F/T − a   F   F/T     —     y   a   x   F/T − a   M   F/T     —     z  
 
+ m   Ft { a   x   Ft   Cβ   o   ÿ   a + a   x   Ft   S β( o   {umlaut over (z)}   a   +g )− a   x   Ft   a   z   Ft {umlaut over (β)}− a   x   Ft   a   y   Ft {dot over (β)} 2 }  (57)
 
By multiplying rotation matrix  k   a R, given in (33), on both side of (53) and  a M a , elements of which are given in (55), (56), and (57), one can get  k F a  and  k M a  as follows:
 
 k   F   ax =− a   F   F/T     —     x   Cθ   2   Cθ   3 − a   F   F/T     —     y ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )+ a   F   F/T     —     z ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )
 
+ m   Ft   {Cθ   2   Cθ   3   o   {umlaut over (x)}   a   +[S (β+θ 1 ) Sθ   2   Cθ   3   +C (β+θ 1 ) Sθ   3 ] o   ÿ   a  
 
+[−C(β+θ 1 ) Sθ   2   Cθ   3   +S (β+θ 1 ) Sθ   3 ]( o   {umlaut over (z)}   a   +g )
 
+[− Sθ   2   Cθ   3 ( a   y   Ft   Cθ   1 + a   z   Ft   Sθ   1 )+( a   y   Ft   Sθ   1 − a   z   Ft   Cθ   1 ) Sθ   3 ]{umlaut over (β)}
 
+ [Sθ   2   Cθ   3 (− a   y   Ft   Sθ   1 + a   z   Ft   Cθ   1 )− Sθ   3 ( a   y   Ft   Cθ   1 + a   z   Ft   Sθ   1 )]{dot over (β)} 2 }  (58)
 
 k   F   ay = a   F   F/T     —     x   Cθ   2   Sθ   3 + a   F   F/T     —     y ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )− a   F   F/T     —     z ( Cθ   1   Sθ   2   Sƒ   3   +Sƒ   1   Cθ   3 )
 
 m   Ft {[−( a   y   FT   Cθ   1 + a   z   Ft   Sθ   1 ) Cθ   3 +( a   y   Ft   Sθ   1 − a   z   Ft   Cθ   1 ) Sθ   2   Sθ   3 ]{dot over (β)} 2  
 
+[( a   z   Ft   Sθ   1 + a   y   Ft   Cθ   1 ) Sθ   2   Sθ   3 +( a   y   Ft   Sθ   1 − a   z   Ft   Cθ   1 ) Cθ   3 ]{umlaut over (β)}
 
− o   {umlaut over (x)}   a   Cθ   2   Sθ   3   +[−S (β+θ 1 ) Sθ   2   Sθ   3   +C (β+θ 1 ) Cθ   3 ] o   ÿ   a  
 
+[ C (β+θ 1 ) Sθ   2   Sθ   3   +S (β+θ 1 ) Cθ   3 ]( o   {umlaut over (z)}   a   +g )}  (59)
 
 k   F   az = a   F   F/T     —     x   Sθ   2 + a   F   F/T     —     y   Sθ   1   Cθ   2 − a   F   F/T     —     z   Cθ   1   Cθ   2  
 
+ m   Ft   {+Sθ   2   o   {umlaut over (x)}   a   −S (β+θ 1 ) Cθ   2   o   ÿ   a   +C (β+θ 1 ) Cθ   2 ( o   {umlaut over (z)}   a   +g )+
 
( a   y   Ft   Cθ   1 + a   z   Ft   Sθ   1 ) Cθ   2 {umlaut over (β)}+( a   y   Ft   Sθ   1 − a   z   Ft   Cθ   1 ) Cθ   2 {dot over (β)} 2 }  (60)
 
 k   M   ax =( a   I   Ft     —     x {umlaut over (β)}− a   M   F/T     —     x ) Cθ   2   Cθ   3  
 
− a   M   F/T     —     y ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )− a   M   F/T     —     z ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )
 
− a   F   F/T     —     x [ a   y   F/T ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )+ a   z   F/T ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )]
 
+ a   F   F/T     —y   [ a   x   F/T ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )+ a   z   F/T   Cθ   2   Cθ   3 ]
 
+ a   F   F/T     —     z [ a   x   F/T ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )− a   y   F/T   Cθ   2   Cθ   3 ]
 
+ m   Ft   o   {umlaut over (x)}   a [ a   y   Ft ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )+ a   z   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )]
 
+ m   Ft   o   ÿ   a { a   x   Ft   [−C (β+θ 1 ) Sθ   2   Cθ   3   +S (β+θ 1 ) Sθ   3 ]−( a   y   Ft   Sβ+   a   z   Ft   C β) Cθ   2   Cθ   3 }
 
− m   Ft ( o   {umlaut over (z)}   a   +g ){ a   x   Ft   [S (β+θ 1 ) Sθ   2   Cθ   3   +C (β+θ 1 ) Sθ   3 ]−( a   y   Ft   Cβ−   a   z   Ft   S β) Cθ   2   Cθ   3 }
 
− m   Ft {umlaut over (β)}{ a   x   Ft [ a   y   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )+ a   z   Ft (− Cθ   1   ,Sθ   2   Cθ   3   +Sθ   1   Sθ   3 )]
 
−( a   y   Ft   2 + a   z   Ft   2 ) Cθ   2   Cθ   3 }
 
+ m   Ft {dot over (β)} 2a   x   Ft [ a   y   Ft ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )+ a   z   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )]  (61)
 
 k   M   ay =−( a   I   Ft     —     x {umlaut over (β)}− a   M   F/T     —     x ) Cθ   2   Sθ   3 − a   M   F/T     —     y (− Sθ   1   Sθ   2   Sθ   3   +Cθ   1   Cθ   3 )
 
− a   M   F/T     —     z ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )
 
+ a   F   F/T     —     x [ a   y   F/T ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )− a   z   F/T (− Sθ   1   Sθ   2   Sθ   3   +Cθ   1   Cθ   3 )]
 
− F   F/T     —     y [ a   x   F/T ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )+ a   z   F/T   Cθ   2   Sθ   3 ]
 
+ a   F   F/T     —     z [ a   x   F/T (− Sθ   1   Sθ   2   Sθ   3   +Cθ   1   Cθ   3 )+ a   y   F/T   Cθ   1   Sθ   3 ]
 
− m   Ft   o   {umlaut over (x)}   a [ a   y   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )− a   z   Ft (− Sθ   1   Sθ   2   Sθ   3   +Cθ   1   Cθ   3 )]
 
+ m   Ft   o   ÿ   a { a   x   Ft   [C (β+θ 1 ) Sθ   2   Sθ   3   +S (β+θ 1 ) Cθ   3 ]+( a   y   Ft   Sβ+   a   z   Ft   C β) Cθ   2   Sθ   3 }
 
− m   Ft ( o   {umlaut over (z)}   a   +g ){ a   x   Ft   [−S (β+θ 1 ) Sθ   2   Sθ   3   +C (β+θ 1 ) Cθ   3 ]+( a   y   Ft   Cβ−   a   z   Ft   S β) Cθ   2   Sθ   3 }
 
− m   Ft {umlaut over (β)}{ a   x   Ft [ a   y   Ft (− Sθ   1   Sθ   2   Sθ   3   +Cθ   1   Cθ   3 )+ a   z   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )]
 
+( a   y   Ft   2 + a   z   Ft   2 ) Cθ   2   Sθ   3 }
 
− m   Ft {dot over (β)} 2a   x   Ft [ a   y   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )− a   z   Ft (− Sθ   1   Sθ   2   Sθ   3   +Cθ   1   Cθ   3 )]  (62)
 
 k   M   az =( a   I   Ft     —     x {umlaut over (β)}− a   M   F/T     —     x ) Sθ   2 +( a   M   F/T     —     y   Sθ   1 − a   M   F/T     —     z   Cθ   1 ) Cθ   2  
 
+ a   F   F/T     —     x ( a   y   F/T   Cθ   1 + a   z   F/T   Sθ   1 ) Cθ   2 − a   F   F/T     —     y ( a   x   F/T   Cθ   1   Cθ   2 − a   z   F/T   Sθ   2 )
 
− a   F   F/T     —     z ( a   x   F/T   Sθ   1   Cθ   2 + a   y   F/T   Sθ   2 )
 
− m   Ft   o   {umlaut over (x)}   a ( a   y   Ft   Cθ   1 + a   z   Ft   Sθ   1 ) Cθ   2   −m   Ft   o   ÿ   a [ a   x   Ft   C (β+θ 1 ) Cθ   2 +( a   y   Ft   Sβ+   a   z   Ft   C β) Sθ   2 ]
 
+ m   Ft ( o   {umlaut over (z)}   a   +g )[ a   x   Ft   S (β+θ 1 ) Cθ   2 +( a   y   Ft   Cβ−   a   z   Ft   S β) Sθ   2 ]
 
+ m   Ft {umlaut over (β)}[( a   y   Ft   2 + a   z   Ft   2 ) Sθ   2 + a   x   Ft ( a   y   Ft   Sθ   1   Cθ   2 − a   z   Ft   Cθ   1   Cθ   2 )]
 
 −m   Ft {dot over (β)} 2a   x   Ft ( a   y   Ft   Cθ   1 + a   z   Ft   Sθ   1 ) Cθ   2   (63)
 
     Substituting (39), (58), (59) and (60) into (50), one can rewrite (50) as follows:
 
 k   F   kx =− a   F   F/T     —     x   Cθ   2   Cθ   3 − a   F   F/T     —     y ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )+ a   F   F/T     —     z ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )
 
+( m   Sk   +m   Ft ) Cθ   2   Cθ   3   o   {umlaut over (x)}   a +( m   Sk   +m   Ft )[ S (β+θ 1 ) Sθ   2   Cθ   3   +C (β+θ 1 ) Sθ   3 ] o   ÿ   a  
 
−( m   Sk   +m   Ft )[ C (β+θ 1 ) Sθ   2   Cθ   3   −S (β+θ 1 ) Sθ   3 ]( o   {umlaut over (z)}   a   +g )
 
− m   Ft {umlaut over (β)}[ a   y   Ft ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )+ a   z   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )]
 
+ m   Sk   L   Sk     —     com   {Cθ   2   Sθ   3 ({umlaut over (θ)} 1 +{umlaut over (β)})− Cθ   3 {umlaut over (θ)} 2   −Sθ   2 [2 Sθ   3 ({dot over (θ)} 1 +{dot over (β)}){dot over (θ)} 2   +Cθ   2   Cθ   3 ({dot over (θ)} 1 +{dot over (β)}) 2 ]}
 
− m   Ft {dot over (β)} 2 [ a   y   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )− a   z   Ft ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )]  (64)
 
 k   F   ky = a   F   F/T     —     x   Cθ   2   Sθ   3 + a   F   F/T     —     y ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )− a   F   F/T     —     z ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )
 
−( m   Sk   +m   Ft ){ Cθ   2   Sθ   3   o   {umlaut over (x)}   a   +[S (β+θ 1 ) Sθ   2   Sθ   3   −C (β+θ 1 ) Cθ   3 ] o   ÿ   a }
 
+( m   Sk   +m   Ft )[ C (β+θ 1 ) Sθ   2   Sθ   3   +S (β+θ 1 ) Cθ   3 ]( o   {umlaut over (z)}   a   +g )
 
+ m   Ft {umlaut over (β)}[ a   y   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )+ a   z   Ft ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )]
 
+ m   Sk   L   Sk     —     com   {Cθ   2   Cθ   3 ({umlaut over (θ)} 1 +{umlaut over (β)})+ Sθ   3 {umlaut over (θ)} 2   +Sθ   2 [−2 Cθ   3 ({dot over (θ)} 1 +{dot over (β)}){dot over (θ)} 2   +Cθ   2   Sθ   3 ({dot over (θ)} 1 +{dot over (β)}) 2 ]}
 
+ m   Ft {dot over (β)} 2 [ a   y   Ft ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )− a   z   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )]  (65)
 
 k   F   kz =− a   F   F/T     —     x   Sθ   2 + a   F   F/T     —     y   Sθ   1   Cθ   2 − a   F   F/T     —     z   Cθ   1   Cθ   2  
 
+( m   Sk   +m   Ft )[ Sθ   2   o   {umlaut over (x)}   a   −S (β+θ 1 ) Cθ   2   o   ÿ   a ]+( m   Sk   +m   Ft ) C (β+θ 1 ) Cθ   2 ( o   {umlaut over (z)}   a   +g )
 
+ m   Ft {umlaut over (β)}( a   z   Ft   Sθ   1 + a   y   Ft   Cθ   1 ) Cθ   2   +m   Ft {dot over (β)} 2 ( a   y   Ft   Sθ   1 − a   z   Ft   Cθ   1 ) Cθ   2  
 
 +m   Sk   L   Sk     —     com [{dot over (θ)} 2   2   +C   2 θ 2 ({dot over (θ)} 1 +{dot over (β)}) 2 ]  (66)
 
     Substituting (39), (58), (59), (60), (61), (62), and (63) into (51) yields
 
 k   M   kx =− a   M   F/T     —     x   Cθ   2   Cθ   3 − a   M   F/T     —     y ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )+ a   M   F/T     —     z ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )
 
− a   F   F/T     —     x [ a   y   F/T ( Cθ   1   Sθ   2   Sθ   3   −Sθ   1   Sθ   3 )+ a   z   F/T ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )− L   Sk   Cθ   2   Sθ   3 ]
 
+ a   F   F/T     —     y [ a   x   F/T ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )+ a   z   F/T   Cθ   2   Cθ   3   +L   Sk ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )]
 
+ a   F   F/T     —     z [ a   x   F/T ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )− a   y   F/T   Cθ   2   Cθ   3   −L   Sk ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )]
 
+{−( m   Ft   L   Sk   +m   Sk   L   Sk     —     com ) Cθ   2   Sθ   3  
 
+ m   Ft [ a   y   Ft ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )+ a   z   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )]} o   {umlaut over (x)}   a  
 
+           ( m   Ft   L   Sk   +m   Sk   L   Sk     —     com )[− S (β+θ 1 ) Sθ   2   Sθ   3   +C (β+θ 1 ) Cθ   3 ]
 
+ m   Ft { a   x   Ft   [−C (β+θ 1 ) Sθ   2   Cθ   3   +S (β+θ 1 ) Sθ   3 ]−( a   y   Ft   Sβ+   a   z   Ft   C β) Cθ   2   Cθ   3 }           o   ÿ   a  
 
+         ( m   Ft   L   Sk   +m   Sk   L   Sk     —     com )[ C (β+θ 1 ) Sθ   2   Sθ   3   +S (β+θ 1 ) Cθ   3 ]
 
− m   Ft { a   x   Ft   [S (β+θ 1 ) Sθ   2   Cθ   3   +C (β+θ 1 ) Sθ   3 ]−( a   y   Ft   Cβ−   a   z   Ft   S β) Cθ   2   Cθ   3 }         ( o   {umlaut over (z)}   a   +g )
 
+( k   I   Sk     —     x   +m   Sk   L   Sk     —     com ) Cθ   2   Cθ   3 ({umlaut over (θ)} 1 +{umlaut over (β)})+( k   I   Sk     —     x   +m   Sk   L   Sk     —     com ) Sθ   3 {umlaut over (θ)} 2  
 
+         [ a   I   Ft     —     x   +m   Ft ( a   y   Ft   2 + a   z   Ft   2 )] Cθ   2   Cθ   3  
 
+ m   Ft   {L   Sk [ a   y   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )+ a   z   Ft ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )]
 
− a   x   Ft [ a   y   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cƒ   1   Sƒ   3 )+ a   z   Ft (− Cθ   1   Sθ   2   Cθ   3   +Sθ   1   Sθ   3 )]}         {umlaut over (β)}
 
+( k   I   Sk     —     y − k   I   Sk     —     z   +m   Sk   L   Sk     —     com ) Cθ   2   Sθ   2   Cθ   3 {dot over (θ)} 1   2 +         ( k   I   Sk     —     y − k   I   Sk     —     z   +m   Sk   L   Sk     —     com ) Cθ   2   Sθ   2   Sθ   3  
 
+ m   Ft   {L   Sk [ a   y   Ft ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )− a   z   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )]
 
+ a   x   Ft [ a   y   Ft ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )+ a   z   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )]}         {dot over (β)} 2  
 
−( k   I   Sk     —     x + k   I   Sk     —     y − k   I   Sk     —     z +2 m   Sk   L   2   Sk     —     com ) Sθ   2   Cθ   3 ({dot over (θ)} 1 +{dot over (β)}){dot over (θ)} 2  
 
+( k   I   Sk     —     x − k   I   Sk     —     y + k   I   Sk     —     z )( Sθ   2   Cθ   3 {dot over (θ)} 2   −Cθ   2   Sθ   3 ({dot over (θ)} 1 +{dot over (β)})){dot over (θ)} 3  
 
+( k   I   Sk     —     y − k   I   Sk     —     z   +m   Sk   L   Sk     —     com   2 ) S (2θ 2 ) Sθ   3 {dot over (θ)} 1 {dot over (β)}  (67)
 
 k   M   ky = a   M   F/T     —     x   Cθ   2   Sθ   3 + a   M   F/T     —     y ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )− a   M   F/T     —     z ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )
 
+ a   F   F/T     —     x [ a   y   F/T ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )+ a   z   F/T ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )+ L   Sk   Cθ   2   Cθ   3 ]
 
− F   F/T     —     y [ a   x   F/T ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )+ a   z   F/T   Cθ   2   Sθ   3   −L   Sk ( Sθ   1   Sθ   2 θ 3   +Cθ   1   Sθ   3 )]
 
+ a   F   F/T     —     z [ a   x   F/T (− Sθ   1   Sθ   2   Sθ   3   +Cθ   1   Cθ   3 )+ a   y   F/T   Cθ   2   Sθ   3   −L   Sk ( Cθ   1   Sθ   2 θ 3   −Sθ   1   Sθ   3 )]
 
−{( m   Sk   L   Sk     —     com   +m   Ft   L   Sk ) Cθ   2   Cθ   3  
 
+ m   Ft [ a   y   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )− a   z   Ft (− Sθ   1   Sθ   2   Sθ   3   +Cθ   1   Cθ   3 )]} o   {umlaut over (x)}   a  
 
+         −( L   Sk     —     com   m   Sk   +m   Ft   L   Sk )[ S (β+θ 1 ) Sθ   2   Cθ   3   +C (β+θ 1 ) Sθ   3 ]
 
+ m   Ft { a   x   Ft   [C (β+θ 1 ) Sθ   2   Sθ   3   +S (β+θ 1 ) Cθ   3 ]+( a   y   Ft   Sβ+   a   z   Ft   C β) Cθ   2   Sθ   3 }           o   ÿ   a  
 
−         ( m   Sk   L   Sk     —     com   +m   Ft   L   Sk )[− C (β+θ 1 ) Sθ   2   Cθ   3   +S (β+θ 1 ) Sθ   3 ]
 
+ m   Ft { a   x   Ft   [−S (β+θ 1 ) Sθ   2   Sθ   3   +C (β+θ 1 ) Cθ   3 ]+( a   y   Ft   Cβ−   a   z   Ft   S β) Cθ   2   Sθ   3 }         ( o   {umlaut over (z)}   a   +g )
 
         −[ a   I   Ft     —     x   +m   Ft ( a   y   Ft   2 + a   z   Ft   2 )] Cθ   2   Sθ   3  
 
+ m   Ft   {L   Sk [ a   y   Ft ( Cθ   1   Sθ   2   Cθ   3   −Sθ   1   Sθ   3 )+ a   z   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )]
 
+ a   x   Ft [ a   y   Ft ( Sθ   1   Sθ   2   Sθ   3   −Cθ   1   Cθ   3 )− a   z   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )]}         {umlaut over (β)}
 
+( k   I   Sk     —     y   +m   Sk   L   Sk     —     com )[− Cθ   2   Sθ   3 ({umlaut over (θ)} 1 +{umlaut over (β)})+ Cθ   3 {umlaut over (θ)} 2 ]
 
+( k   I   Sk     —     x − k   I   Sk     —     z   +m   Sk   L   Sk     —     com   2 ) Cθ   2   Sθ   2   Cθ   3 [({dot over (θ)} 1   2 +{dot over (β)} 2 )+2{dot over (θ)} 1 {dot over (β)}]
 
+ m   Ft   {L   Sk [ a   y   Ft ( Sθ   1   Sθ   2   Cθ   3   +Cθ   1   Sθ   3 )+ a   z   Ft (− Cθ   1   Sθ   2   Cθ   3   +Sθ   1   Sθ   3 )]
 
− a   x   Ft [ a   y   Ft ( Cθ   1   Sθ   2   Sθ   3   +Sθ   1   Cθ   3 )− a   z   Ft (− Sθ   1   Sθ   2   Sθ   3   +Cθ   1   Cθ   3 )]}{dot over (β)} 2  
 
+( k   i   Sk     —     x − k   I   Sk     —     z +2 m   Sk   L   Sk     —     com   2 ) Sθ   2   Sθ   3 {dot over (θ)} 2 ({dot over (θ)} 1 +{dot over (β)})
 
+( k   I   Sk     —     x − k   I   Sk     —     y − k   I   Sk     —     z ){dot over (θ)} 3   [Sθ   3 {dot over (θ)} 2   +Cθ   2   Cθ   3 ({dot over (θ)} 1 +{dot over (β)})]
 
− k   I   Sk     —     y   Sθ   2   Sθ   3 {dot over (θ)} 2 {dot over (β)}  (68)
 
 k   M   kz =− a   M   F/T     —     x   Sθ   2 +( a   M   F/T     —     y   Sθ   1 − a   M   F/T     —     z   Cθ   1 ) Cθ   2  
 
+ a   F   F/T     —     x ( a   y   F/T   Cθ   1 + a   z   F/T   Sθ   1 ) Cθ   2  
 
− a   F   F/T     —     y ( a   x   F/T   Cθ   1   Cθ   2 − a   z   F/T   Sθ   2 )
 
− a   F   F/T     —     z ( a   x   F/T   Sθ   1   Cθ   2 + a   y   F/T   Sθ   2 )
 
− m   Ft   o   {umlaut over (x)}   a ( a   y   Ft   Cθ   1 + a   z   Ft   Sθ   1 ) Cθ   2  
 
− m   Ft   o   ÿ   a [ a   x   Ft   C (β+θ 1 ) Cθ   2 +( a   y   Ft   Sβ+   a   z   Ft   C β) Sθ   2 ]
 
+ m   Ft ( o   {umlaut over (z)}   a   +g )[ a   x   Ft   S (β+θ 1 ) Cθ   2 +( a   y   Ft   Cβ−   a   z   Ft   S β) Sθ   2 ]
 
+ k   I   Sk     —     z   [Sθ   2 ({umlaut over (θ)} 1 +{umlaut over (β)})+{umlaut over (θ)} 3 ]
 
+{ a   I   Ft     —     x   Sθ   2   +m   Ft [ a   x   Ft ( a   y   Ft   Sθ   1 − a   z   Ft   Cθ   1 ) Cθ   2 +( a   y   Ft   2 + a   z   Ft   2 ) Sθ   2 ]}{umlaut over (β)}
 
 +[−m   Ft   a   x   Ft ( a   y   Ft   Sθ   1 + a   z   Ft   Cθ   1 )+( k   I   Sk     —     y − k   I   Sk     —     x ) Cθ   2   Sθ   3   Cθ   3   ]Cθ   2 {dot over (β)} 2  
 
+ Cθ   2 [ k   I   Sk     —     z +( k   I   Sk     —     y − k   I   Sk     —     x ) C (2θ 3 )]({dot over (θ)} 1 +{dot over (β)}){dot over (θ)} 2  
 
+( k   I   Sk     —     x − k   I   Sk     —     y ) Sθ   3   Cθ   3 [{dot over (θ)} 1   2   C   2 θ 2 −{dot over (θ)} 2   2 ]
 
+2( k   I   Sk     —     x − k   I   Sk     —     y ) C   2 θ 2   S (2θ 3 ){dot over (θ)} 1 {dot over (β)}  (69)

     Thus, all moments and forces at knee and ankle are obtained. With close observation of (53), (55), (56), (57), (64), (65), (66), (67), (68) and (69), it is clear that the knee and ankle moments and forces can be computed in real-time with online measurements of four angles (i.e., footplate angle (β), dorsiflexion/planar-flexion angle θ 1 , inversion/eversion angle θ 2 , shank internal/external rotation angle θ 3 ), which can be measured with the 6-DOF goniometer in  FIG. 5  and  FIG. 6 , and three forces and three moments acting on foot (i.e.,  a F F/T     —     x ,  a F F/T     —     y ,  a F F/T     —     z ,  a M F/T     —     x ,  a M F/T     —     y ,  a M F/T     —     z ), which can be measured with the 6-axis F/T sensor. Note that, two of the moments applied from footplate to foot ( a M F/T     —     x  and  a M F/T     —     y ) are not required for calculation of moments of lower limb joints, if foot is not strapped to footplate (Winter, 2005). 
     Of course, mass of shank and foot (m Sk , m Ft ), inertias of shank ( k I Sk     —     x ,  k I Sk     —     y ,  k I Sk     —     z ) and foot ( a I Ft     —     x ), and segment lengths (L Sk , L Sk     —     com ) are also required. However, these are constant values for each subject and can be measured/estimated offline. Masses and inertias can be estimated by referring anthropometric data (Winter, 2005; Zatsiorsky, 2002), and segment lengths can be either measured using a simple digitizer (e.g., MicroScribe or a digitizing probe of Optotrak) or estimated by referring anthropometry (Winter, 2005; Zatsiorsky, 2002). 
     Variables needed to be measured from the goniometer and force/torque sensor below the foot in real-time for calculation of each joint moment and for are listed in  FIG. 7  as a summary. Clearly, one can see depends on the joint and forces/moments of interest, number of required variables vary. 
     3. Feasibility Testing Experiment 
     1) A Sample Experimental Setup 
     A modified ET (e.g., Reebok Spacesaver RL) was instrumented with 6-axis F/T sensors (e.g., JR3 Inc., Woodland, Calif.) on both sides underneath the footplate. The ankle joint center (JC), midpoint between the medial and lateral malleoli, was aligned with the F/T sensor center. One end of the 6-DOF goniometer was attached to the footplate and the other end was strapped to the bony frontal-medial surface of shank. To corroborate the estimated Knee Adduction Moment (KAM) using the 6-DOF goniometer, an Optotrak 3020 system with 3 sets of markers (4 markers/set attached to a rigid shell) attached on thigh, shank, and foot was also used. The F/T and goniometer data were sampled at 1 KHz and the Optotrak data at 50 Hz. The foot was strapped to the footplate with no relative motions to the footplate. As an initial calibration, the subject stood still and initial knee JC coordinates with respect to the ankle frame was determined. From anthropometry data (Winter, 2005; Zatsiorsky, 2002), lengths—between knee JC and COM of shank and between knee JC and ankle JC—and masses and inertias with respect to COM of foot and to that of shank were obtained. With all the measured/estimated kinematic and kinetic variables, KAM was computed with respect to the knee local coordinate system and displayed as external KAM following the convention in (Foroughi et al., 2009). 
     2) Results 
     The KAM during stepping movement was estimated on a healthy volunteer (male, age: 32; height: 181 cm; mass: 97 kg). Strong KAM was generated during the stepping movement and the moment varied systematically with the stepping cycle ( FIG. 6  ( a )). When the subject adducted the knee during stepping, the KAM increased markedly ( FIG. 6(   b )). Furthermore, KAM estimated using the 6-DOF goniometer matched closely with that using the Optotrak and especially the peak KAM locations in each cycle with the 6-DOF goniometer and that with Optotrak are almost identical, thereby confirming the effectiveness of the proposed method.