Research Article

INTERNATIONAL JOURNAL OF ADVANCED ROBOTIC SYSTEMS

Robot motion synthesis using ground reaction forces pattern: Analysis of walking posture

1 2 Ramamoorthy Luxman and Teresa Zielinska

Abstract

Gait generation in its realization stage in majority of bipeds involves inverse dynamics model with the need of pre-determined joint trajectories and the zero moment point method for final postural stability adjustments. This requires generation of such joint trajectories, where the zero moment point criterion is fulfilled and the human-like posture is kept. Such task is complex and still misses the universal solution. Preserving natural (human-like) and stable body posture under acting disturbances is an issue too. High computational costs, problem with disturbances rejection and controller complexity are the major disadvantages associated with such approaches. This article presents a method for motion control with forward dynamics model and ground reaction forces as the reference values. The ground reaction forces measured during real human walking are utilized here. A MATLAB/Simulink framework is used for testing the proposed synthesis of bipedal locomotion. The ground reaction forces recorded during the walk with different footwear are applied as the references. The obtained motion patterns and postures are compared with those of the human. Obtained results were satisfactory and justified the proposed approach.

International Journal of Advanced Robotic Systems July-August 2017: 1-12 © The Author(s) 2017 DOI: 10.1177/1729881417720873 journals.sagepub.com/home/arx

(DSAGE

Keywords

Bipedal locomotion, biped dynamics, human ground reaction forces, joint torques, motion synthesis

Date received: 31 August 2016; accepted: 20 June 2017

Topic: Climbing and Walking Robots Topic Editor: Manuel Armada

Introduction

Bipedal locomotion (human gait) is an eminent topic in research. Creating human-like walking or running in robots is still in its infancy and has not been satisfactorily realized. Many theories have been proposed to establish standardized models for biped dynamics, gait generation as well as gait control. Bipedal locomotion is one type among many locomotion methods observed among living organisms. Birds (when on the ground), kangaroos, ostriches, chimpanzees and humans are some of the vertebrates that exhibit bipedal locomotion. Ostriches and kangaroos can reach a speed as high as 70 km/h and have very robust and highly efficient gaits. Ostriches are facilitated with highly compliant tendons in their legs, which enable them to adopt an effective bouncing gait in the absence of aerial phase

(ground running).1 This concept is employed for developing bionic boots with springs and elastic properties similar to that of an ostrich's leg. With these bionic boots, it is possible to attain higher running speed. An et al.2 demonstrate the analytical modelling of elbow using the muscle length-tension relationship. The model determines the muscle force distribution across the elbow joint for various

1 Politechnika Warszawska, Warszawa, Poland MEiL, Politechnika Warszawska, Warszawa, Poland

Corresponding author:

Ramamoorthy Luxman, Politechnika Warszawska, Nowowiejska 24,

Warszawa 00-665, Poland.

Email: l.ramamoorthy0l@gmail.com

Creative Commons CC BY: This article is distributed under the terms of the Creative Commons Attribution 4.0 License (http://www.creativecommons.org/licenses/by/4.0/) which permits any use, reproduction and distribution of the work without further permission provided the original work is attributed as specified on the SAGE and Open Access pages (https://us.sagepub.com/en-us/nam/ open-access-at-sage).

configurations. Andrews et al.3 present Euler's and Lagrange's equations for three-dimensional rigid links human motion model. A technique applying inverse dynamic model for estimating the joint forces and force moments with kinematic and dynamic constraints consideration has been proposed in the study by Koopman et al.4 Ren et al.5 present a three-dimensional whole body multi-segmental model for inverse dynamics analysis over a complete gait cycle. The relationship between leg-end forces and gait properties is not yet sufficiently studied. This issue is important for gait synthesis of real bipedal robots. Furthermore, the reaction forces influence the robot posture which is needed for the stable gait. In recent years, greater focus is paid to the dynamic constraints on a walking robot motion.6 Especially when the robot is moving over a complex terrain, the fulfilment of the contact force constraints7 must be focused. In the study by Yu etal.,8 the human ground reaction forces are used for obtaining the reference trajectory of zero moment point (ZMP) which is next applied for synthesis of humanoid gait pattern. In humanoids besides designing the legs displacement, there is a need to adjust the movement of other parts of the body for obtaining human-like and stable walking posture. It is obvious that the time courses of contact forces result from motion dynamics of the whole body parts. Many diverse approaches have been applied for the mapping of human motion data in humanoids.910 Some of these approaches are synthesis of motion by dynamics simulation, optimization method, model predictive control, breaking the motion into sequence of motion primitives and using backstepping integral control, and so on.

Having the above in mind, we decided to apply the ground reaction force (GRF) trajectories and forward dynamics model for gait generation. In the recent years, greater focus is paid to the dynamic constraints existing during robot motion 11 for motion generation. This approach is reverse to that commonly used where first the joint (or body parts) trajectories are established and next the reaction forces and the ZMP are monitored. This approach allows to obtain the human-like body posture easily. Following the observation that the contact forces are inherently related with the whole body posture maintained during locomotion,12 we propose the method for motion synthesis using forward dynamics model and GRF trajectories. The validation of the method is made by testing it on how the body posture is properly adjusted to maintain human-like walking. For that purpose, we have considered the gaits with different footwear (what obviously affects the posture). This work also provides novel insights on how joint trajectories, needed torques and trunk position of a biped have to be adjusted depending on the robot foot in order to achieve human-like walking.

The method was validated using two-dimensional model of the robot body.

Assumptions

The modelling of the robot motion is done by projecting the whole biped system onto a single two-dimensional X-Yplane. It is assumed that mass is uniformly distributed across the body segments. The shape of the links (segments) is approximated to a perfect cylinder, and the centre of mass (CoM) and moments of inertia are evaluated accordingly. Joints are assumed to be frictionless. Legend is assumed to be point feet. Using the current state of the system containing joint angles, the instantaneous moment is calculated using the GRF trajectories and knowing the lengths of the links. A database containing recorded human joint trajectories, CoM trajectories (for results validation) and GRF for different footwear types (sneakers, high heel shoes and barefoot) are considered. Recorded GRFs are used as the input feed forces in our simulation models. Suitable Propotional Integral Derivative (PID) controllers are designed for joint torques generation during leg transfer phase, where the GRF is equal to zero. The results obtained from the simulation for different scenarios are compared (particularly the trunk inclination). An easy to use MATLAB graphical user interface (GUI) is developed with inputs for delivering the GRF and with adjustable control gains.

Firstly, kinematic model of the biped system with definition of reference frames and Degree Of Freedoms (DOFs) is established. Denavit-Hartenberg (DH) parameters are tabu-larized, considering each limb as a separate single open loop kinematic chain. Direct kinematics problem (DKP) and inverse kinematics problem (IKP) solutions are then derived with these DH parameters. Subsequently, dynamic modelling of the biped system is carried out. Dynamical modelling is performed using Lagrangian mechanics. The model is actuated by moments about the hip and knee axes. Incorporation of the model in Simulink was carried out using GRFs in real human walking as the input. The data were collected using clinically validated systems for the three-dimensional gait assessment of children and adults - VICON MX T-Series system. The system uses passive optical motion capturing technology, which is the most accurate and common method for professional recording of human data in hospitals and research centres. Force plates, one of the elements of this system, measure the leg-end forces. The advanced software of this system automatically processes the data, eliminating the noise and smoothening the data. The data used in this project were recorded for a healthy female of 28 years old with no motion disorders. The subject weighed 53.4 kg at the time of experiment and is 1.66 m tall. Before gathering the experimental data, the barefoot gait of the tested person was analyzed and diagnosed as normal, without any deviations by the vision package of the used VICON system. Simulink model is developed with appropriate definition of ground surface properties. A GUI with an easy-to-use interface is also developed, accompanied by animation of the simulated biped locomotion for

^CoM ^CoM

Z 0,Z 1,Z2,Z 3 Normal to paper

n0, n1

Hip joint '■. 'i

Knee Joint '2

<XL1,yii>l

(XL2,^L2)

CR yn i)

left & right thighs

(XR2, yR2)

left & right shanks

Figure 1. Kinematics diagram. (a) Five links biped system model. (b) For analysis of system dynamics.

(XoM - fsina)

- ^cosa)

IKP is solved geometrically. The solution for joint angles for a given foot position is obtained as follows

Y = a tan2

X 2 IV2 _/2_ /2

where C = Xfoot ^ 1 /2

P = atan2(ifoot,Xfoot)- atan2(W2, Wi)

where W1 and respectively.

are (/1 + /2 cosy) and (/2 sinY),

visualization and control/ configuration tools for simulation. Mathematical model of the five link biped system is built according to Hyotyniemi and Haavisto13 with modifications according to our objective.

Kinematics

The biped system consists of two open loop kinematic chains (left and right limbs) - Figure 1(a). DKP transformation matrices are obtained with notations a, / and y as joint rotation angles (subscripts r and / are put to represent right and left leg, respectively)

cosP — sinP 0 0

sinP cosP 0 0

0 0 10

0 0 0 1

' cos(01 + 02) sin(p + y) 0 0

cos(p + y) sin(p + y)

— sin(p + y) 0 /1 cosp

cos(p + y) 0 /1 sinp

— sin(p + y) cos(p + y) 0 0

/2 cos(P + y) /2 sin(P + y) 0 1

Dynamics

The link dimensions and masses are taken from anthropo-metric data as tabulated in Table 1.

The seven independent, complete and holonomic coordinates constitute state vector q7- (refer Figure 1(b))

q = [x0; y0 > a pr; y/ , Yr ] (7)

Generalized forces include all the non-conservative forces acting on the biped system. Each co-ordinate is associated with the corresponding generalized force.

Generalized forces, Q, are

Q = [FX0, Fy0, Fa, FP/, FPr , FY/, FYr]

- /1 cosp'

- /1 sinp

The Trunk CoM position in the system state vector can hence be transformed to hip, knee and foot position using the above derived transformation matrix

CoM of left leg (xLg, yLg), thigh and shank are expressed as follows. For right leg, equations are the same with corresponding joints rotation angles

xL1 = x0 — r0 sina — r1 sin(a — /L) (9)

Jli = yo — ro sina — ri cos(a — /¿) (10)

xL2 = x0 — r0 sina — li sin(a — /L) — r2 sin(a — + yl)

JL2 = X0 — r0 cosa — /1 cos(a — /¿) — r2 cos(a — + yz)

xLg = x0 — r0 sina — /1 sin(a — /L) — /2 sin(a — + yl)

= x0 — r0 cosa — /1 cos(a — /L) — /2 cos(a — + yl)

Ankle joint

Table 1. Biped system configuration (total body mass, m = 55.33 kg) (magnitudes of the subject human body in the real walking experiments).

Torso Thigh Shank

Link length (m) lo = tÍT = 0.788 l, = £ = 0.406 li = & = 0.473

CoM distance (m) ro = 2 = 0.394 r i = 2 = 0.203 r2 = = 0.236

Link masses (kg) mo - = 0.68m = 37.63 mi = 0. im = 5.53 m2 = = 0.06m = 3.32

Link inertia (kgm2) Jo = m0l2 = 1.9483 J, = m,|2 = 0.076 J2 = m2l2 = 0.0619

CoM: centre of mass.

The total kinetic energy of the system T is

T = 1 (m0(x20 + 3>2)+mi(iL1 + fLl + + yRl

+ m2(xL2 + 3_L2 + xR2 + y2Rl)

Generalized forces are

Fx2 = FL + FR

displacement. The detailed dynamic equations are given in Appendix 1.

(15) Estimation of joint torques

Having obtained the GRF [FJ, FR, FJ, FR]T (Left (L) and right (R) leg end forces, along x and y axes), the joint

(16) torques [MLhip, MRip, MLknee, MR"ee]

Fyo = -(mo + 2m i + 2m2)g + FL + FR (17)

/SyLi . S/L2 . Sy—- Sy—- SyR2 . Fq = — —- mi +—- m2 mi +—- mi + —- m2 g

Sa Sa Sa Sa Sa

+ ^ fL + ^ FR + SîLl Fl + SXRg Fr

F, = - ^ m, + ^ m^ g + g. FL + H FL + M?

SyL2 , SyL^y , SxL

fyl = - qL m2g+qL FL+iYLr F + MLr (22)

(equations for Fo (right leg) are similar to (19) and (18)) R

where joint torques vector (M) and GRF vector (F) are

m = [MLhip, M—ip, mL^, mR^]

F = [FL, FL, FR, FR]

(2i) (22)

d / ST \ ST = Q

dt lSq J Sq,

The set of seven equations of second differentials can be written in matrix form as

I(q)q = b(q, q, M, F)

where vector b is a (7 x 1) function of joint torques and GRFs at the leg end. It represents time derivative of partial of L with respect to velocity and partial of L with respect to

are estimated simultaneously. Mathematical equations for estimating joint torques are obtained as follows:

Support and swing phases of each leg are detected simply by checking the instantaneous value of the GRF. If GRF (Fx, Fy) of a leg are zero, then it is apparent that the leg is in swing phase. A non-zero GRF indicates that the leg is in contact with the ground. Hence, estimation of joint torques in each case is carried out separately.

In order to estimate internal joint torques, net joint torques are calculated from the GRF and the torques due to the external gravitational and frictional forces that are subtracted from that. The internal torques are not known.

Net joint torque for hip and knee joints are calculated as given below

T hep = F grf (z2 sin( p - arctan ^ Ff - « - P + y)

Placing these appropriately, we obtain the Lagrangian of the system as

P /Fgr

+ l, sin ( — — arctan I —

i 1 2 Vf gr

in /Fgrf Fgrf 12 sin--arctan

V 2 V Fxgrf

a — ß

a — ß + Y

Net torque derived as above is the resultant of internal torques and the torque caused by external forces including gravitational pull and friction force. Hence internal torques, M1 and M2, are calculated by equating these forces

Tnet -(m-g sin(a + ß))

M2 = Tkntee -(m2g sin(« + $ + )) (28)

where is static friction coefficient. A detailed derivation is given in Appendix 2.

(Support phase) Calculated moments [ML, 0, 0]

[0, 0, MR, MR ]

E(s) PID

-' C(s)

(Swing phase)

Dynamic simulator

Figure 2. Scheme of the system control.

In swing phase, when the legs are off the ground, the GRFs are equal to zero. Hence, it is not possible to apply the above derived expressions for calculating the torques.

In order to calculate the needed torque in transfer phase, a PID controller c(s) is defined. The current state of the system is given directly as the feedback and reference state corresponding to the current time instant that is given as input to the PID controller, which gives the estimated

Figure 3. Considered leg-end forces; rows are for different footwear (barefoot, sneakers and high heels) and columns are for left foot (left column) and right foot (right column).

Figure 4. Stick diagrams illustrating animation of synthesized motion.

Figure 5. Position of trunk CoM depends on footwear: vertical (Y) coordinate - left column, trunk inclination (angle measured towards Y axis) - right column. CoM: centre of mass.

swing phase torques as output. Proportional, integral and Figure 2 illustrates the joint torques obtained using the derivative gains are tuned using Ziegler-Nichols method in PID controller, which are furnishing the biped dynamic the simulation. transfer function block.

Figure 6. Joint angles - angles plotted as referenced in Figure 1(b). Output of the PID controller, M(t) is

M(t)=Kp£(t)+K,-

E(t) dt + Kd

dE(t) dt

The proportional, integral and derivative gains Kp, K,-and Kd are chosen suitably based on the control results. The chosen values are Kp = between (2.5 and 3), K,- = between (0.75 and 1.2) and Kd = between (0.2 and 0.6).

Results

The simulations are done using MATLAB (Simulink). Experimental data sets containing the GRFs measured during one complete gait cycle of a real human walking are used as force input to this simulation model. The plots

below represent the results obtained in simulations. The presented work demonstrates a robust human-like gait generation using GRFs in contrast to the widely used methods based on inverse dynamic modelling. Joint torques evaluation from the GRFs is straightforward and thus obtaining the joint movements using the forward dynamics is computationally inexpensive. Considered GRFs are illustrated in Figure 3

Discussions on results

Torso positioning

In the plots showing the trunk behaviour (Figure 4) obtained during simulations, it was detected that the

maximal trunk inclination occurs by the end of single support phase and the minimal inclination is during the double support phase close to the beginning of the next single support phase. The trunk behaviour slightly alters depending on the supporting leg (left or right). It is apparent that this phenomenon varies from person to person. It is observed that it is more symmetrical for walking with sneakers than for walking barefoot.

There is a notable difference in the range of trunk inclinations exhibited in the three walking scenarios (Figure 5). Walking barefoot

f Right foot take off left foot landing : 2 .7° to 2 .5° range = <

[ Left foot take off right foot landing : 1 .2° to 3 .7° Walking with sneakers

f Right foot take off left foot landing : 3 .2° to 5 .9°

range =

[ Left foot take off right foot landing : 2 .8° to 6 .2° Walking with high-heel footwear

Right foot take off left foot landing : 2 :2° to 5 :2° Left foot take off right foot landing : —1 .8° to 5 .7°

This concludes that the torso is subjected to maximal range of inclinations (sway) when walking on high heels (3°L/7.5°R), and minimal when walking in barefoot (1.8°L/2.5°R), the range of inclination in case of sneakers is found to be (2.9°L/3.2°R).

In case of walking barefoot and with sneakers, the trunk is subjected only to forward sway. Whereas, walking with a high-heeled footwear, the trunk is leaned backward (negative angle). The highest forward inclination is observed in case of sneakers and the least in case of barefoot.

We can observe local fluctuations in torso trajectories. In case of walking with high heels, this fluctuation is prominent. This indicates that the role of torso is vital for maintaining the postural dynamic stability, particularly in the case when the gait stability is obstructed, which in our case was the walking with high heels.

Knee rotation angle

From Figure 6, it is observed that the range of knee rotation angle is —22° to 52°. In the case of high heels walking, the range is slightly smaller than walking barefoot and with sneakers (—22° to 42°). It is also seen that the rotation is relatively smooth in the case of barefoot and sneakers, while in the case of high heel walking, it shows high variability.

Stride length

From the results, it is observed that the stride length of the walking gait is relatively less in case of high heels as compared to the cases of sneakers and barefoot. However, the difference is not very significant. Among sneakers and barefoot, barefoot walking exhibits slightly lesser stride length.

Also, the less oscillatory trajectories of torso inclination and of knee joint in the cases of barefoot and sneakers as compared to high heel shows better postural stability. Those results are confirming that in the cases of GRF-based motion control, the reference forces must be carefully defined. Moreover, the force trajectories influence the trunk posture and the smoothness of joint trajectories.

Conclusions

Obtained results are comparable to the real human walking. Figure 6(b) shows the human data recorded for walking with sneakers and presents the simulation and the real data. The real and simulated joint trajectories are very similar. The case is the same for the other considered examples. This proves the correctness of the presented formal descriptions and confirms the usability of direct dynamic modelling for force-based gait generation in real bipeds. The method delivers not only the needed torques, human-like joint trajectories but also the human-like trunk position adjustments. Applied by us, data illustrate how the footwear influences the motion properties. Such information can be useful for the designers of the robotic feet. The footwear affects the gait, dynamic postural equilibrium and trunk motion. The permanent use of some type of the shoes results in temporary or permanent body adjustments, depending on the time and frequency of the shoes used.14 According to the author's knowledge, very few works are currently available that have studied the relation between GRF and dynamics of the body. These are mostly the research works with bio-medical context and are rather general. Alternative methods for swing phase torque estimation, postural balance considering human walking GRF, investigation of gait with GRF analysis are the themes of future research both for robotics and for bio-mechanics investigating the reasons of gait disorders.

Authors' note

We had full access to all of the data in this study and we take complete responsibility for the integrity of the data and the accuracy of the data analysis

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

1. Rubenson J, Heliams DB, Lloyd DG, et al. Gait selection in the ostrich: mechanical and metabolic characteristics of walking and running with and without an aerial phase. Proc Biol

Sci 2004; 271(1543): 1091-1099. D01:10.1098/rspb.2004. 2702.2004.

2. An KN, Kaufman KR and Chao EYS. Physiological considerations of muscle force through the elbow joint. Biomech J 1989; 22(11-12): 1249-1256. Pergamon Press.

3. Andrews JG. Euler's and Lagrange's equations for linked rigid body models of three dimensional human motion. In: Allard P, Stokes IAF and Blanchi JP (eds) Three dimensional analysis of human movement. Champaign, IL: Human Kinetics, pp. 145-176.

4. Koopman B, Grootenboer HJ and de Jongh HJ. An inverse dynamics model for the analysis, reconstruction and prediction of bipedal walking. JBiomech 1995; 28(11): 1369-1376.

5. Ren L, Jones RK and Howard D. Whole body inverse dynamics over a complete gait cycle based only on measured kinematics. JBiomech 2008; 41(12): 2750-2759. D0I:10. 1016/j.jbiomech.2008.06.001.

6. Bin Hamman G, Wensing PM, Dariush B, et al. Kinodynami-cally consistent motion retargeting for humanoids. /nt J Human Robot 2015; 12(04): 27. D0I:10.1142/S0219843615500176.

7. Dai H, Valenzuela A and Tedrake R. Whole-body motion planning with centroidal dynamics and full kinematics. In: 2014 /EEE-RAS International conference on humanoid robots, Madrid, Spain, 18-20 November 2014, pp. 295-302. IEEE. DOI: 10.1109/HUMAN0IDS.2014.7041375.

8. Yu Z, Chen X, Huang Q, et al. Humanoid Walking Pattern Generation Based on the Ground Reaction Force Features of Human Walking. In: HMmanoid walking pattern generation based on the ground reaction force features of human walking proceeding of the /EEE International Conference on information and Automation Shenyang, Shenyang, China, 16-18 June 2012, pp. 573-578. IEEE. DOI: 10.1109/ ICInfA.2012.6246919.

9. Yamane K and Hodgins J. Control-aware mapping of human motion data with stepping for humanoid robots. In: 2010 /EEE/RSJ International Conference on Intelligent Robots and Systems, 18-22 October 2010, Taipei, Taiwan, pp. 726-733. IEEE. DOI: 10.1109/IR0S.2010.5652781.

10. Forsyth DA, Arikan O, Ikemoto L, et al. Computational studies of human motion: part 1, tracking and motion synthesis. Found Trend Banan Comput Graph Vis 2005; 1(2-3): 77-254. D0I:10.1561/0600000005.

11. Zielinska T and Chmielniak A. Biologically inspired motion synthesis method of two-legged robot with compliant feet. Robotica 2011; 29: 1049-1057. Bananamark® Cambridge University Press 2011. D0I:10.1017/S0263574711000300.

12. Augusto F Barbieri and Vitrio R. Locomotion andposture in older adults: The role of aging and movement disorders. Springer, 2017.

13. Hyotyniemi H and Haavisto O. Simulation tool of a biped walking robot model. PhD thesis, Helsinki University of Technology, Helsinki, Finland, 2004. https://www.research gate.net/publication/237512585_SIMULATI0N_T00L_ 0F_A_BIPED_WALKING_R0B0T_M0DEL.

14. Gueguen N, Morio C, Lakec MJ, et al. The influence of footwear on foot motion during walking and running. JBiomech 2009; 42(13): 2081-2088.

Appendix 1

Dynamic model of the biped system considered in our simulator is formulated as

I(q) = b(q, q, M, F)

Inertia matrix I(q) and vector b are derived to be as follows I(q) :

111 = m + 2(m1 + «2)

112 = 2

113 = (—2m1ro — 2m2ro) cosa

+ (—11m2 — m1r1) cos(a — fa) — 11 m2 cos(a — )

114 = (l1m2 + m1r1) cos(a — fa) + m2r2 cos(a — fa + yi)

115 = (l1m2 + m1r1) cos(a — fa)+ m2r2 cos(a — fa + Yr)

116 = —m2r2 cos(a — fa + Yi)

117 = —m2r2 cos(a — fa + Yr)

I21 = 2

123 = (2^1^ + 2^2^) sina + ('1^2 + «1^) sin(a — fa)

+ l1m2 sin(a — fa)

124 = (—/1^2 — m1n) sin (a — fa)— m2r2 sin(a — fa + Yl)

125 = (—'1^2 — M1n) sin (a — fa)— m2r2 sin(a — fa + Yr)

126 = «2 r2 sin (a — fa + Y' )

127 = M2T2 sin(a — fa + Yr)

131 = (—2m1r2 — 2m2r2) cosa +(—l1m2 — m1r1) cos(a — fa)

— l1m2 cos(a — fa) — m1r1 cos(a — fa)

— «2^2 cos(a — fa + y') — «2^2 cos(a — fa + yr)

132 = (2^1^ + 2^2^) sina + (l1m2m1 n) sin(a — fa)

+ Z1m2 sin(a — fa) + m1r1 cos(a — fa) + m2r2 sin(a — fa + Y') + m2r2 sin(a — fa + Yr)

133 = 2('12TO2 + + «2r22 + «1^1 2 + «2r22)

+ 2(l1m2r2 + «1^1) cosfa + 2('1rn2r2 + «1^1) cosfa. + 2m2 r2r2cos(fa — yl) + 2'^^cosy' + 2m2r2r2 cos(fa — yr) + 2'1 m2r2 cosyr

134 = — '12m2 — m1^12 — «2r22 + (—'1^2^ — «1^1) cosfa

— «2^ 2^ 2 cos(fa — y') — 2l1m2r2cosy'

135 = — '1 2«2 — «1^12 — M2^22 + (—^«2^ — «1^1) cosfa

— «2r 2r 2 cos(fa — yr) — 2l1«2r2 cosyr

136 = «2^2 + ^ cos(fa — Yl) + '1 cosYl

137 = «2r2(r2 + r2 cos(fa — Yr) + '1 cosYr

I41 = C1OT2 + «1^) cos(a — fa) + «2r 2 cos(a — fa + yi )

142 = (-1i«2 - miri) sin (a - ß; )- m2r2 sin(a - ß; + y; )

143 = — li2m2 — miri2 — m2r22 + (—lim2ro — mirori) cosß;

- m2ror2cos(^; — 2lim2r2COSy;

144 = li2m2 + miri2 + m2r22 + 2lim2r2 cosy;

145 = 0

146 = m2r2(-r2 - li cosyi)

147 = 0

151 = (lim2 + miri) cos(a — ßr)+ m2r2 cos(a — ßr + Yr)

152 = (—lim2 — miri) sin (a — ßr) — m2r2 sin(a — ßr + yr)

153 = — li2m2 — miri2 — m2r22 + (—lim2ro — mirori) cosßr

— m2ror2 cos(ßr — yr) — 2lim2r2 cosyr

154 = o

155 = li2m2 + miri2 + m2r22 + 2lim2r2 cosYr b

156 = o

157 = m2r2(-r2 - li cosYr)

161 = -m2r2 cos(a - ßl + Yl)

162 = m2r2 sin (a - ß; + Yl )

163 = m2r2(r2 + ro cos(ßl - Yl) + li cosYl

164 = m2r2(-r2 - li cosYl)

165 = o

166 = m2r22

167 = o

- Yr2m2r2 sin(a - ßr + Yr) — ß 12m2r2 sin(a - ßl + y;)

I7i = -m2r2 cos(a - ßr + Yr) -

r!(n — r)!

172 = m2r2 sin(a - ßr + Yr)

173 = m2r2(r2 + ro cos(ßr - Yr) + li cosYr

174 = o

175 = m2r2(-r2 - li cosYr)

176 = o

177 = m2r22

b(q, q, F, M) :

bi = — 2<â2miro sina + FRR — ßr2miri sin(a — ßr)

— a2lim2 sin(a — ßl)+ FLx

— Y12m2r2 sin(a — ßl + Yl) — a2miro sin(a — ßr)

— ßr2lim2 sin(a — ßr) — i2lim2 sin(a — ßr)

— ß12miri sin(a — ßl)ßr2m2r2 sin(a — ßr + yr)

— ai2miri sin(a — ßl)— ß12lim2 sin(a — ßl)

— a2m2r2 sin(a — ßr + yr)

+ ci2m2r2 sin(a — ßl + y;) + 2ißrlim2 sin(a — ßr) + 2ßrYrr2m2 sin(a - ßr + yr)

+ 2c(rßlr2m2 sin(a — ßl + y;) + 2iß^m^ sin(a — ßl) + 2aßrr2m2 sin(a — ßr + yr) + 2ißll1m2 sin(a — ß;)

— 2am2ro sina + 2ißrlim2 sin(a — ßr) + 2ß;Y;m2r2 sin(a - ß; + y;)

— 2a;ylm2r2 sin(a - ß; + y;) 2 = -ßl2m2r2 cos(a - ß; + y;)+ FR - 2gm2

— Yr 2m2r2 cos (a - ßr + yr )- ß 2m2r2 cos(a - ßr + yr )

— a2m2r2 cos(a - ßr + yr)- Y;2m2r2 cos(a - ß; + y;)

— a2rimi cos(a — ßr)— ßr2rimi cos(a — ßr)

— a2(2mi + 2m2)ro cosa — ci2lim2 cos(a — ßr)

— ßr2lim2 cos(a — ßr)+ FlL + ci2m2r2 cos(a — ß; + y;)

— 2gm1 + 2cißrl1m2 cos(a — ßr)

— 2ßr-yrr2m2 cos(a - ßr + yr)

— 2a Yrm2 r2 cos(a — ßr + yr ) + ßlY;m2r2 cos(a - ß; + y;)

+ 2ciß;m2r2 cos(a — ß; + y;) + 2ciß^m^ cos(a — ßr)

— (ci ß ; (—2/im2 — 2miri) + i2(lim2 + miri) + 'ßi2(lim2 + miri)) cos(a - ß;)- gmo

= Yr2lim2r2 sinyr + FRRl2 sin(a - ßr + yr)

— Fjl1 cos(a - ß;)- FRl1 cos(a - ßr)

— FJ;2 cos(a - ßr + y;) + F£ro sina + FRro sina + F^Ll1 sin (a - ß; ) + F^ sin (a - ßr + y;) 2gm2r2 sin(a - ßr + yr)

— Y r 2m2ror2 sin(ßr - Yr )- ßr 2m2ror2 sin(ßr - Yr)

— gm2r2 sin (a - ß; + y; )+ yyl 2lim2 r2 siny;

— Y;2m2ror2 sin(ßl - y;)- ßr2mirori sin(ßr)

+ P/2m2r0r2 sin(P/ — Y/) — Pr2/1m2r0 sin(Pr)

— P12m2r0r2 sin(P/ — Y/) — P/2m1r0r1 sinP/

— g/1m2 sin (a — Pr) + FRR/1 sin(a — Pr)

— gm1r1 sin(a — Pr)— g/1m2 sin(a — P/)

— gmm sin (a — P;) — P ; 2/1 m2r0 sinP/

— (FJr0 + FRRr0) cosa — 2aY/M2r0r2 sin(P/ — Y/)

+ 2/3/Yy/m2r0r2 sin(P/ — y/)+ 2aP/M2r0r2 sin(P/ — y/) + 2aPrM1r0r1 sinPr — 2gm2r0 sina + 2<â;yr/1m2r2 sinYr + 2aPrm2r0r2 sin(Pr — Yr)

— 2aYr«2r0r2 sin(Pr — Yr) + 2aY^^^ sinY/

— 2P5/Y//1 m2r2 sinY/ — 2gm1r0 sina

+ 2aPr/1m2r0 sinPr + 2aP//1m2r0 sinP/ + 2aP/m1r0r1 sinP/ — 2P5r;yr/1m2r2 sinYr + 2/àr'yrm2r0r2 sin(Pr — Yr)— FRR/2 cos(a — Pr + Yr) b4 = MLP + FJ/1 cos(a — P/)+ FL/2 cos(a — P/ + y/ )

— sin(a — P/)+ g/1m2 sin(a — P/) + gm1r1 sin(a — P/) — a2/1m2r0 sinP/

— a2m1r0r1 sinP/ — a2m2r0r2 sin(P/ — Y/) + 2aY//1m2r2 sinY/ + 2PP/Y//1M2r2 sinY/

— Y / 2/1«2 r 2 sinY/ — FLL/2 sin (a — P/ + y/ ) + gm2r2 sin(a — P/ + y/)

b5 = Mr^'' + FR/1 cos(a — Pr)

+ FRx/2 cos(a — Pr + Yr )— FRR/1 sin(a — Pr ) + g/1m2 sin (a — Pr)+ gm1r1 sin(a — Pr )

— a2/1m2r0 sinPr — a2m1 r0r1 sinPr

— a2m2r0r2 sin(Pr — Yr) + 2<â-7r/1m2r2 sinYr

+ 2P5rYr/1m2r2 sin-r — Yr2/1m2r2 sin-r

— FRR/2 sin (a — Pr + Yr )+ gm2r2 sin(a — Pr + -r ) b6 = M/"ee — FL/2cos(a — P/ + y/)

+ <i2m2r0r2 sin(P/ — y/)— a2/1 m2r2 sin-/ + 2aP//1m2r2 sin(-/)— P//1m2r2 sin-/ + FLL/2 sin(a — P/ + y/)— gm2r2 sin(a — P/ + y/)

Figure 2A. Geometrical interpretation - for net torque calculation.

b7 = M/"ee — FR/2 cos (a — Pr + Yr)

+ a2m2r0r2 sin(Pr — Yr) — a2/1m2r2 sin-/ + 2aPr/1M2r2 sin(Yr)— Pr/1M2r2 sin-r + FRR/2 sin (a — Pr + Yr )— gm2r2 sin(a — Pr + -r )

Appendix 2

Support phase

In order to estimate internal joint torques, net joint torques are calculated from the GRF and the torques due to external gravitational and frictional forces are subtracted from it. Forces are identified as known forces and unknown forces. Known forces include force due to body weight, force induced due to friction with the ground and GRFs. Unknown forces include the forces associated with the internal torques about the hip and knee joints. All these forces are carefully resolved into their x and y components for formulating the torque equations.

Net joint torques in Figures 2A and 2B represent the geometrical interpretation of the net torque calculation. The net joint torques about hip and knees are calculated as follows

Torque, t = F x R

/ Thigh\

cW J /

(ß fay Shank H a+ß + Y

\Thigh

1_____\ psFgt-L -J

m j g m 2 g

(a) (b)

_ net — Fgrf T hip — F

_net _ Fgr/

T knee = F

r1 and r2 are calculated geometrically. In the diagram, 6>2 and 6>3 are known (from state vector). They are

02 = a + ß

03 = y

01 = arctan I ^g-f

VFf-/y

AABG-AEG7, hence

0i = 06

Finally

=>■ 06 = arctan —

04 = 2 - 06 - 02

05 = 04 + 03

) 05 = 2 - 0i - 02 + 03

-2 = /2 Sin05

-2 sin04

sin04 -2

+ /1 + /1

Figure 2B. Torque around (a) hip and (b) knee.

It is to be noted that at knee joint, the force is just being transferred between leg end and hip. There is no additional force produced in the knee joint. Hence, the net moment about hip and knee are

) -1 = -2 + (/1 sin04)

Therefore

/ • /p t

-2 = /2 sin--arctan

2 2 v2 ^f/

a — ß + y

( (p /Ff-A

-1 = -2 + /1 sin--arctan — a — ß

1 2 v1 v2 ^F/ ßy

From (31) and (32)

Fg-/ ^/2 sin( p - arctan ^ F/j - a - ß + y^

+ /1 sin I 2 - arctan

(!)- a - ß)

net rknee

Fg-//2 sinf p - arctan f F/ j - a - ß + y^

Interna/ torques. The net torques derived above are the result of internal torques and the torque driven by external forces including gravitational pull and friction. Hence, internal torques, M1, M2, are calculated by equating these forces.

About hip: Equating the summation of all torques about hip to the net hip torque

T nep = M1 + (m1g sin (a + fa) ) M1 = tJ? -(m, x g x sin(a + fa) (34)

About knee: equating the torques about the knee

: M2 + (m2g sin (a + ß)) + (MiFff )

) M2 = tknL -(m2g sin(a + ^) + (MsFff)) (35)

Equations (34) and (35) are the two required equations for calculating the internal torques during the support phase.