Khalili Imad

Dufaud Olivier

Perrin Laurent

Nicolas Cuervo

Laurent André

Interactions turbulence-combustion lors des explosions de mélanges hybrides vapeurs/poussières -Combustion-turbulence interactions in hybrid mixtures flames

Keywords: Mélange hybride, explosions de poussières, étirement, longueur de Markstein Hybrids mixtures, dusts explosions, flame stretching, Markstein length

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

L'évaluation du risque d'explosion des nuages de poussières combustibles nécessite la connaissance préalable de leurs caractéristiques de sensibilité à l'ignition et de sévérité. Elles dépendent de nombreux paramètres dont, entre autres, la composition chimique du matériau, la granulométrie, la richesse et la turbulence du mélange combustible -comburant [START_REF] Eckhoff | Dust explosions in the process industries[END_REF]. La situation se complique lors de la prise en compte d'un nuage comportant conjointement une poussière inflammable, un gaz ou une vapeur combustible et un comburant. Ces mélanges dits hybrides sont notamment rencontrés dans les industries alimentaires (poussières de céréales et gaz de fermentation ; poussières de tourteaux et vapeur de solvant…) et les industries pharmaceutiques (excipients et solvants). Ils constituent un domaine encore insuffisamment exploré. Pourtant le risque d'explosion des mélanges hybrides est important, tant par leur probabilité d'occurrence que par leur gravité. Lors de l'évaluation des risques et en particulier lors de l'application des directives ATEX, visant à améliorer la sécurité des travailleurs en présence d'atmosphères explosives, le cas des mélanges hybrides n'est souvent pas apprécié à sa juste valeur. D'une part, lorsque la teneur en gaz ou de vapeur est faible, c'est-à-dire inférieure à la limite inférieure d'explosivité (LIE), son influence ne devrait pas être négligée. D'autre part, lorsque la concentration de poussière est inférieure à la concentration minimale d'explosion (CME), l'analyse considérant la seule présence du gaz n'est pas satisfaisante. En effet quelques études [START_REF] Bartknecht | Dynamic behavior of premixed flame fronts in laminar and turbulent flows[END_REF][START_REF] Nifuku | Ignitability assessment of shredder dusts of refrigerator and the prevention of the dust explosion[END_REF][START_REF] Dufaud | Dust/vapour explosions: Hybrid behaviours?[END_REF]) ont démontré la spécificité des explosions de mélanges hybrides solide/gaz ou vapeur par rapport aux explosions soit de gaz ou de vapeur seul, soit de solide seul. Le but de cet article est de déterminer ces caractéristiques spécifiques afin de mieux comprendre les comportements explosifs des mélanges hybrides vapeur/poussières et gaz/poussières. Durant ce travail de recherche, l'étude s'est focalisée sur des composés organiques. Dans une première étape les caractéristiques de la sensibilité à l'ignition (notamment l'Energie Minimale d'Inflammation) et de la sévérité d'explosion (pression d'explosion et vitesse de montée en pression) des mélanges hybrides ont été déterminées expérimentalement. Ensuite la mise en oeuvre d'une technique de mesure de la vitesse de flamme a permis d'examiner les mécanismes de propagation et d'accélération d'une flamme se propageant à travers des mélanges hybrides.

Dispositif expérimental et méthodes

Poudres et leurs préparations

Les tests ont été effectués sur l'amidon de blé acheté chez Sigma Aldrich, tandis que les tourteaux de tournesol ont été fournis par la société Cooperl. La distribution granulométrique des échantillons a été déterminée par granulométrie à diffraction laser (Malvern). Les échantillons ont été caractérisés par les diamètres moyens en volume d 10 , d 50 et d 90 , soit [16 ; 35 ; 76 µm] et [23, 132, 431 µm] respectivement pour l'amidon de blé et les tourteaux. Pour éviter l'effet de l'humidité, les échantillons ont été systématiquement séchés pendant deux heures à 50 ° C sous vide dynamique.

Détermination de l'Energie Minimale d'Inflammation

L'influence de l'addition de petites quantités de gaz ou de vapeur combustible sur la sensibilité de la poussière à l'inflammation a été étudiée en utilisant un tube de Hartmann (Mike 3 -Kühner AG) initialement conçu pour les poussières selon la norme EN 13821:2002. Cet appareil a été modifié afin de permettre la détermination de l'énergie minimale d'inflammation (EMI) de mélanges hybrides.

Mesure des caractéristiques de sévérité

ܵ ݑ ൌ ܵ ݏ ܣ ܣ ݂ (1)
Afin d'évaluer la valeur de A f , la surface du front de flamme peut être considérée comme celle d'un paraboloïde de révolution allongé. La vitesse S u a été corrigée en tenant compte du facteur d'expansion thermique ߙ entre les gaz chauds issus de la combustion et les gaz frais initiaux. Le facteur d'expansion est défini par le rapport entre la masse volumique des réactifs frais ρ u et celle des produits brûlés ρ b . Il a été estimé à l'aide des températures adiabatiques de flamme calculées par le logiciel CEA [START_REF] Gordon | Computer program for calculation of complex chemical equilibrium compositions and applications[END_REF].

ܵ ௨ ൌ ௌ ೞ ఈ avec ߙ ൌ ఘ ೠ ఘ ್ (2)
Plusieurs facteurs influencent la vitesse de combustion comme entre autres la composition du mélange, la pression, la température et la turbulence. L'effet de ce dernier facteur se traduit par la prise en compte d'un terme appelé « étirement ». Il s'agit de l'augmentation de la surface de flamme lorsque le front est soumis à une perturbation dans l'écoulement. Par définition, le taux d'étirement K s'écrit [START_REF] Law | Combustion physics[END_REF]:

ܭ ൌ ଵ ௗ ௗ௧ (3)
Récents Progrès en Génie des Procédés, Numéro 101 -2011ISSN 1775-335X -ISBN 2-910239-75-6, Ed. SFGP, Paris, France Axx-4 [START_REF] Markstein | Nonsteady flame propagation[END_REF]) et Clavin (1985) ont établi une relation linéaire entre la vitesse fondamentale de flamme à étirement non nul S u et le taux d'étirement K : Su = -L K + Su 0 (4) S u ° est la vitesse fondamentale à étirement nul. L'évaluation de la pente de la variation de la vitesse S u en fonction du taux d'étirement K donne la longueur de Markstein L (figure 4).

Résultats et discussion

Sensibilité à l'inflammation des mélanges hybrides

La figure 2 résume les résultats des expériences effectuées sur les mélanges hybrides tourteaux de tournesol / hexane et amidon / hexane dans l'air. Une diminution significative de l'EMI est obtenue dès la présence d'une faible concentration en hexane de l'ordre de 1% vol. Par exemple, pour l'échantillon d'amidon, l'EMI chute de 67 mJ pour la poussière pure à 4 mJ alors que seulement 1% vol. d'hexane est ajouté. Ces résultats corroborent ceux trouvés par [START_REF] Franke | Bestimmung der Minderstzudenergie von Kohlenstaub/Methan/Luft Gemisches (hybride. Gemische)[END_REF]. L'application de modèles de combustion à coeur rétrécissant suggère que l'une des hypothèses susceptible d'expliquer le comportement spécifique de ces mélanges hybrides pourrait être la modification de l'étape déterminante de la réaction de combustion, résultant d'un changement d'une réaction contrôlée par la diffusion à une réaction en phase gazeuse homogène. Nos expériences ont également confirmés qu'une explosion de mélanges hybrides est possible lorsque les concentrations en gaz/vapeurs et en poussières sont respectivement inférieures à la limite inférieure d'explosivité et à la concentration minimale explosive. Des comparaisons ont été faites avec les modèles de Bartknecht et Le Chatelier [START_REF] Bartknecht | Dynamic behavior of premixed flame fronts in laminar and turbulent flows[END_REF]. Il apparaît que ces lois ne sont pas suffisamment restrictives en termes de sécurité. Des modèles simples intégrant le diamètre critique des noyaux de flamme et la température adiabatique de flamme ont été testés avec succès afin de représenter la sensibilité à l'inflammation des ces composés.

Sévérité d'explosion des mélanges hybrides

La figure 3a illustre l'influence de la présence du méthane sur la vitesse de montée en pression. Par exemple, pour une concentration de 200 g.m -3 d'amidon la vitesse de montée en pression est de 31 bar.s -1 en l'absence de méthane, alors qu'avec une concentration de 2% en méthane elle croit jusqu'à 300 bar.s -1 . Le même phénomène est observé sur la figure 3b pour un mélange hexane -tourteaux de tournesol. Il faut aussi noter que la concentration en poussière pour laquelle la sévérité maximale est obtenue décroît nettement lors de l'ajout de quelques pourcents de solvants ou de gaz. Des tests conduits à richesse constante en mélanges hybrides ne conduisent néanmoins pas à des sévérités constantes : la concentration en combustible n'est évidemment pas le seul facteur influant. Pour des richesses inférieures à l'unité, l'addition de faibles quantités de combustible gazeux accroît donc la sévérité des explosions. Dans quelques cas particuliers, la sévérité de l'explosion peut être plus importante pour un mélange hybride Axx-5 que pour le gaz ou la vapeur seule [START_REF] Dufaud | Explosions of vapour/dust hybrid mixtures: A particular class[END_REF]. C'est le cas pour le mélange hexane -tourteaux (figure 3b -tableau 1) contrairement au mélange méthane -amidon.

Figure 3. Vitesse maximale de montée en pression des mélanges hybrides tournesol/hexane (droite -3b) et amidon/méthane (gauche -3a)

Ces indications concordent avec la tendance observée pour l'EMI (paragraphe 3.1). L'une des principales conclusions est que la présence de gaz ou de vapeur dans un nuage de poussière modifie la cinétique de combustion de poussière. Toutefois, cet effet n'est pas isolé et d'autres influences peuvent être mises en exergue. Ainsi, la présence de particules solides dans la zone de préchauffage située en amont des flammes modifie le flux thermique rayonné et influence probablement la vitesse de propagation de la flamme. Les interactions accrues entre les instabilités hydrodynamiques et la propagation des flammes dans les mélanges hybrides ne doivent également pas être négligées. Par rapport aux flammes de prémélanges gazeux, la présence de particules apporte une modification notable du champ de turbulence. Les paramètres de la turbulence dans les écoulements diphasiques gaz/solide sont alors différents de ceux existants dans les écoulements gazeux. L'étude des interactions turbulence-combustion a donc été entreprise (paragraphe 3.3). Des vitesses de flamme à étirement nul de 28 cm.s -1 sont trouvées pour les mélanges de la figure 4a et de 22 à 36 cm.s -1 pour les mélanges avec l'hexane. Ces résultats sont cohérents avec les données de la Axx-6 littérature proposant une vitesse fondamentale de combustion du méthane qui atteint au maximum 36 cm.s -1 pour une richesse unitaire et des vitesses fondamentales de l'ordre de 38 cm.s -1 pour l'hexane [START_REF] Mannan | Lee's loss prevention in the process industries[END_REF]. Ils sont également en accord avec les travaux de Proust ( 2004) et de Dahoe [START_REF] Dahoe | Determination of the laminar burning velocity and the Markstein length of powder/air flames[END_REF] donnant des vitesses laminaires de flamme à étirement nul comprises entre 15 et 30 cm.s -1 pour l'amidon seul. La longueur de Markstein est, elle aussi, spécifique à chaque mélange. Il faut d'abord souligner que les valeurs de L sont certes faibles mais positives, ce qui conduit à des nombres de Lewis (ratio de la diffusivité thermique sur la diffusivité moléculaire) supérieurs à l'unité et, par conséquent, à une atténuation du plissement de la flamme et à sa stabilisation. Les ordres de grandeurs obtenus sont également cohérents avec les valeurs de 0,1 à 0,2 mm obtenues pour les mélanges air/méthane par Dahoe [START_REF] Dahoe | Determination of the laminar burning velocity and the Markstein length of powder/air flames[END_REF]. L'apport de particules combustibles tend probablement à accroître L, dont la valeur tend vers 8,5 ou 11 mm selon les références [START_REF] Proust | Formation -Inflammation -Combustion des atmosphères explosives (ATEX) et effets associés[END_REF]) [START_REF] Dahoe | Determination of the laminar burning velocity and the Markstein length of powder/air flames[END_REF]. Ce paramètre doit être aussi comparé à l'épaisseur de flamme δ. Les essais à richesse constante reproduits en figure 4b démontrent que ce seul paramètre ne permet pas de déterminer la stabilité de la flamme et que le ratio entre les concentrations de combustibles solides et gazeux est également un paramètre influent. L'analyse plus fine et systématique de ces données permettra de déterminer les conditions de déstabilisation et donc d'extinction de ces flammes de mélanges hybrides.

La

Conclusion

Cette étude démontre le comportement spécifique des mélanges hybrides. La présence de gaz ou de vapeur dans le nuage d'une poussière combustible influence à la fois la sévérité d'explosion et la sensibilité à l'inflammation, même lorsque la concentration de vapeur ou de gaz est inférieure à sa LIE. La présence d'une faible concentration de gaz/vapeur, ici de l'ordre de 1% vol., peut conduire à un changement de l'étape déterminante de la réaction de combustion passant d'une limitation en régime de diffusion à une limitation réactionnelle ce qui implique une diminution drastique de l'EMI. Il a également été démontré que la présence simultanée des poudres et des vapeurs modifie la vitesse de la flamme qui devient beaucoup plus élevée que la vitesse de combustion de poussière pure. L'effet est principalement attribuable à une modification du processus de combustion mais aussi aux interactions turbulencecombustion ainsi qu'à la modification des modes de transfert thermique en amont du front de flamme.

2. 4

 4 Mesure et calcul de la vitesse fondamentale de flamme à étirement nul La propagation de la flamme a été étudiée à l'aide d'une enceinte d'expérimentation constituée d'un tube sous vide étanche de format parallélépipède, de section droite carrée, de volume V= 4,9 L, de longueur L=100 cm et de côté D= 7cm, soit un rapport d'allongement L/D d'environ 14 (Figure1). Le tube est fermé à l'extrémité de l'allumage et équipé d'une soupape de sécurité à l'autre extrémité. Comme deux des faces opposées du tube sont en verre (les deux autres étant en acier inoxydable), il est possible de suivre visuellement le processus de propagation de la flamme à l'aide d'une caméra vidéo à haute vitesse à 10000 images par secondes. Le système d'allumage est constitué de deux électrodes de l'appareil Mike 3. Ces électrodes ont été positionnées au milieu de l'enceinte à 11 cm au-dessus de l'extrémité inférieure du tube. La mise en suspension des poussières se fait par le système de génération de suspension du Mike 3. Le dispositif expérimental comprend aussi un système d'injection de vapeurs et de gaz permanents.

Figure 1 .

 1 Figure 1. Dispositif expérimental de visualisation de propagation de flammesLa position du front de flamme est déterminée par traitement numérique des vidéos. Elle permet de déterminer la vitesse spatiale de propagation du front de flamme. Selon[START_REF] Andrews | Determination of burning velocities: a critical review[END_REF] la relation entre la vitesse spatiale de la flamme S S et la vitesse de combustion S U peut être exprimée comme suit en fonction de la surface du front de flamme A f et de la surface de la section droite du tube A.

Figure 2 .

 2 Figure 2. Evolution des énergies minimale d'inflammation des mélanges hybrides tourteaux de tournesol/hexane et de l'amidon/hexane en fonction de la concentration de vapeur.

  vitesse fondamentale de flamme à étirement nul La figure 4 illustre l'évolution de la vitesse fondamentale de la flamme en fonction du taux d'étirement du front de flamme. Au début de la propagation de la flamme lorsque sa distance parcourue est faible, le taux d'étirement de sa surface est grand et la vitesse de propagation de la flamme est faible. On constate alors figure 4 que le taux d'étirement de la flamme diminue avec l'accroissement de la vitesse fondamentale. La vitesse fondamentale de flamme à étirement nul est obtenue en appliquant la relation (4).

Figure 4 .

 4 Figure 4. Evolution de la vitesse fondamentale de la flamme étirée en fonction de l'étirement pour des mélanges hybrides tourteaux de tournesol/hexane (droite -4b) et amidon/méthane (gauche -4a).

  Les mesures de la surpression maximale d'explosion P max , de la vitesse maximale de montée en pression (dP/dt max ) et de la concentration minimale explosive (CME) des combustibles ont été conduites dans une sphère d'explosion de 20 L (Kühner AG) normalisée (ISO 6184/1). Le tableau 1 regroupe les caractéristiques propres de sévérité de chacun des constituants des mélanges étudiés afin de définir une référence de comparaison.

		Tableau 1 : Caractéristiques de sévérité des combustibles purs
	Combustibles	CME (g.m -3 ) /	P max (bars)	dP/dt max (bar.s -1 )	K st/g (bar.m.s -1 )
		LIE (% v.)
	Tourteaux de tournesol	40	7,4	452	123
	Amidon de blé	150	8,8	397	108
	Hexane	0,3	8,3	1182	321
	Méthane	4,8	8,2	1802	489
			Axx-3


Hayder F N Al-Shuka

Burkhard Corves

Wen-Hong Zhu

The aim of this report is to derive the equations of motion for biped robot during different walking phases using two well-known formulations: Euler-Lagrange (E-L) and Newton-Euler (N-E) equations. The modeling problems of biped robots lie in their varying configurations during locomotion. They could be fully actuated during the single support phase (SSP) and overactuated during the double support phase (DSP). Therefore, first, the E-L equations of 6-link biped robot are described in some details for dynamic modeling during different walking phases with concentration on the DSP. Second, the detailed description of modified recursive Newton-Euler (N-E) formulation (which is very useful for modeling complex robotic system) is illustrated with a novel strategy for solution of the over-actuation/discontinuity problem. The derived equations of motion of the target biped for both formulations are suitable for control laws if the analyzer needs to deal with control problems. As expected, the N-E formulation is superior to the E-L concerning dealing with high degrees-of freedom (DOFs) robotic systems (larger than six DOFs).

Dynamic modeling

Humans have perfect mobility with amazing control systems; they are extremely versatile with smooth locomotion. However, comprehensive understanding of the human locomotion is entirely still not analyzed. Please see [Hay19, Hay18/1, Hay18/2, Hay18/3, Hay18/4, Hay17/1, Hay17/2, Hay16, Hay15, Hay14, Hay14/1, Hay14/2, Hay14/3, Hay14/4, Hay13/1, Hay13/2, Hay13/3, Sam08] for more details on dynamics, walking pattern generators and control of biped locomotion (biped robots, lower-extremity exoskeletons, prosthetics, etc.). To dynamically model the ZMPbased biped mechanisms, the following points should be considered:  Biped robots are kinematically varying mechanisms such that they could be fully actuated during the SSP and over-actuated during the DSP. If we assume the biped robot as fixedbase mechanism, the dynamic modeling and control strategies of fixed-base manipulators can efficiently be used.

 Dealing with unilateral contact of the foot-ground interaction as a passive joint (rigid-torigid contact) or as compliant model (penalty-based approach).

 Reducing the number of links/joints of the target biped as possible. However, they can still have more than six DOFs resulting in computational problems of advanced control systems.

 Reducing the walking phases as much as possible, e.g., most conventional ZMP-based biped robots (see Tab. 2.2 of [START_REF] Al-Shuka | walking pattern generators and adaptive control of biped robot[END_REF]) can walk with two substantial walking phases: the SSP and the DSP. Adjustments of the walking patterns are possible by modification of foot design as described in [START_REF] Sato | Trajectory planning and control for biped robot with toe and heel joint[END_REF].

 Most ZMP-based biped robots walks with flat swing /stance feet all the time; this can facilitate the analysis of biped locomotion by reducing walking phases to exactly two phases: the SSP and the DSP (see ref. [START_REF] Vanderborght | Overview of the Lucy project: Dynamic stabilization of a biped powered by pneumatic artificial muscles[END_REF]). However, heel-off/toe-off sub-phases can offer better characteristics but with careful analysis as we will see later.

In the light of the above comments, classical Euler-Lagrange (E-L) equations and recursive Newton-Euler (N-E) can be used for dynamic modeling of biped robots. For complex robotic systems, such as humanoid robots or any robot having the number of degrees of freedom (DOFs) larger than 6 DOFs, difficulties are encountered in the implementation of the control algorithms. Therefore, over 30 years, the robotics researchers have focused on the problem of computational efficiency. Many efficient O (𝑛) algorithms have been developed for inverse [Sah99] and forward dynamics [START_REF] Mohan | A recursive, numerically stable and efficient simulation algorithm for serial robots[END_REF] of robotic systems. For more literature on the efficient dynamic algorithms, refer to refs. [START_REF] Khalil | Dynamic modeling of robots using recursive Newton-Euler formulations[END_REF]. The adaptive control algorithm, however, which deals with controlling the robotic systems despite their uncertain parameters may decrease the computational efficiency of the dynamics O (𝑛) algorithms. Fu et al. [START_REF] Fu | Robotics: control, sensing, vision, and intelligence USA[END_REF] have shown that the combined identification and control algorithms can be computed in O (𝑛 3 ) time despite using the recursive (N-E) formulation.

One of the efficient tools to deal with full-dynamics-based control for complex robotic systems is the virtual decomposition control (VDC) suggested by Zhu [START_REF] Zhu | H Virtual decomposition control: towards hyper degrees of freedom Berlin[END_REF]. It is equivalent to the recursive NE formulation if the dynamic parameters of the target robotic system are known.

This work deals with the ZMP-based biped robot as a fixed-base robot with rigid foot-ground interaction. In addition, E-L equations are described in some details for dynamic modeling of the biped during different walking phases; problems of over-actuation/ discontinuity are resolved.

Then detailed description of the VDC is illustrated with a novel strategy for solution of the overactuation/discontinuity problem. The remainder of this report is organized as follows. Selection of the walking patterns suggested throughout the current work is presented in Section 1.1. Section 1.2 deals with detailed modeling of biped robot using the E-L equations and VDC. Section 2 concludes.

Selection of walking patterns

As mentioned earlier, different walking patterns can be constructed according to point view of the designer. In general, two walking patterns will be investigated throughout the current work; therefore, this report is interested with dynamic modeling of these walking patterns. The details of the DOFs for the referred walking patterns are as follows.

(i) According to Fig. 1-1 (a), the walking pattern 1 has six generalized coordinates without constraints; the biped behaves as an open-chain mechanism. Consequently, the biped has six DOFs in this walking phase with six links (neglecting the stance-foot link). Whereas, it has seven generalized coordinates with seven links during the DSP due to the rotation of the front foot, but with two constraints equations; the tips of the front and rear feet are fixed(see Fig. 1-1). Consequently, the biped has five DOFs during this constrained walking phase, DSP, with six actuators (over-actuated system). (ii) The walking pattern 2 has also six DOFs with six links during the SSP; it has the same configuration of walking pattern 2. The first sub-phase of DSP (henceforth called DSP1) has six generalized coordinates with two constraint equations; therefore, the biped has 4 DOFs with 6 actuators (over-actuated system). During the second sub-phase of DSP (henceforth called DSP2), similar configurations of that of DSP1 appear and consequently the biped has four DOFs and six actuators.

Dynamic modeling

As mentioned previously, two approaches are commonly used to obtain the differential equations of motions: E-L and N-E formulations. Depending on the purpose of the analyzer/designer, various forms of modified formulations have been conducted such as recursive E-L equations [START_REF] Hollerbach | A recursive lagrangian formulation of manipulator dynamics and comparative study of dynamics formulation complexity IEEE[END_REF], recursive N-E formulation [START_REF] Mohan | A recursive, numerically stable and efficient simulation algorithm for serial robots[END_REF], and generalized D'Alembert (G-D) formulation [START_REF] Lee | Development of the generalized d'Alembert equations of motions for mechanical manipulators[END_REF]. This work concentrates on formulating the dynamic equations that are suitable for adaptive control purposes. Throughout the current analysis, the following assumptions have been proposed.

Assumption 1-1. The stance foot, link (1), is in full contact with the ground during the SSP; therefore, its dynamics could be neglected in such case.

Assumption 1-2. The foot-ground contact is rigid-to-rigid contact. Accordingly, the tips of the feet (in case of foot rotation) are assumed passive joints.

Assumption 1-3.There are only two substantial walking phases, the SSP and the DSP, with possibly sub-phases during the DSP. The instantaneous impact event is avoided by making the swing foot contact the ground with zero velocity (disadvantage of zero end velocity at impact is the higher energy consumption due to need for braking in swing leg).

In biped systems, three important aspects should be taken into consideration:

(i) Preventing the biped legs from slippage.

(ii) Avoiding discontinuities of the ground reaction forces which can result in discontinuities of the actuator torques as detailed in Section 1.2.1.3 (iii) Concentrating on the adaptive control of the biped robot associated with less computational complexity.

The Euler-Lagrange formulation

Although the E-L equations can provide closed-form state equations suitable to advanced control strategies, their computational complexity, unless it is simplified, could be inefficient for analysis/control of complex robotic system (more than 6 DOFs) [START_REF] Zhu | H Virtual decomposition control: towards hyper degrees of freedom Berlin[END_REF]. In fact, E-L formulation could be used recursively [START_REF] Hollerbach | A recursive lagrangian formulation of manipulator dynamics and comparative study of dynamics formulation complexity IEEE[END_REF]; it has been equivalent to recursive N-E formulation in most all aspects [START_REF] Spong | Robot dynamics and control[END_REF]. In general, the computational complexity of E-L, N-E and G-D are 𝑂 (𝑛 4 )𝑜𝑟 𝑂(𝑛 3 ), 𝑂(𝑛) and 𝑂(𝑛 3 ) respectively. Below we present modeling of biped robot during the two phases: the SSP and the DSP with two different kinds of Lagrange equations.

E-L equations of the second kind (the SSP)

The E-L equations for open chain mechanism (biped robot during the SSP) can be expressed as where 𝐿 is Lagrangian function which is equal to the kinetic energy of the robotic system (𝐾) minus its potential energy (𝑃), 𝑞 𝑖 denotes the generalized coordinates of link (i), and 𝑞̇𝑖 is the derivative of the generalized coordinates.

The generalized coordinates are a set of coordinates that completely describes the location (configuration) of the dynamic systems relative to some reference configuration [START_REF] Fu | Robotics: control, sensing, vision, and intelligence USA[END_REF]. There are many choices to select these generalized coordinates; however, the joint/link displacements are proved being suitable in case of robotic systems. If the number of these generalized coordinates is equal to the degrees of freedom of the target system, then Eq. 1-1 is valid; Eq. 1-1 is called Lagrange equations of the second kind and it suitable for open-chain mechanism. Solution of Eq. 1-1 can result in the following second order differential equations. where 𝑴 ∈ ℝ 𝑛 𝑞 ×𝑛 𝑞 is the mass matrix, 𝒒, 𝒒̇ and 𝒒̈∈ ℝ 𝑛 𝑞 are the absolute angular displacement, velocity and acceleration of the robot links, 𝑪 ∈ ℝ 𝑛 𝑞 ×𝑛 𝑞 represents the Coriolis and centripetal robot matrix, 𝒈 ∈ ℝ 𝑛 𝑞 ×1 is the gravity vector, 𝑨 ∈ ℝ 𝑛 𝑞 ×𝑛 𝜏 is a mapping matrix derived by the principle of the virtual work, 𝝉 ∈ ℝ 𝑛 𝜏 ×1 is the actuating torque vector, 𝑛 𝜏 represents the number of actuators, and 𝝉 𝑓 ∈ ℝ 𝑛 𝜏 ×1 represents the dissipative torques resulted from joint friction.

In the following, some details are presented to determine the dynamic coefficient matrices of Eq. 1-2.

(i) Inertia matrix

The derivation, structure and properties of the mass matrix (𝑴) can be obtained from the total kinetic energy of the biped system. The velocity wrench, 𝒗 ̅ 𝑖 ∈ ℝ 6×1 , of link (𝑖) can be expressed in term of Jacobian matrices as follows [START_REF] Spong | Robot dynamics and control[END_REF][START_REF] Fu | Robotics: control, sensing, vision, and intelligence USA[END_REF].

𝒗 ̅ 𝑖 = [ 𝒗 𝑖 𝒘 𝑖 ] = [ 𝑱 𝑣𝑖 𝑱 𝑤𝑖 ] 𝒒̇ Eq. 1-3
with 𝒗 𝒊 ∈ ℝ 3×1 and 𝒘 𝑖 ∈ ℝ 3×1 are the translational and angular velocity components of link (𝑖), 𝑱 𝑣𝑖 ∈ ℝ 3×𝑛 𝑞 and 𝑱 𝑤𝑖 ∈ ℝ 3×𝑛 𝑞 denote the Jacobian matrices associated with translation and rotation components respectively.

The total kinetic energy of 𝑛 -DOF robotic system can be expressed as Eq. 1-6

(ii) Coriolis and centripetal terms Following the detailed derivation of [START_REF] Spong | Robot dynamics and control[END_REF][START_REF] Fu | Robotics: control, sensing, vision, and intelligence USA[END_REF], without showing the details here, the (𝑘, 𝑗) 𝑡ℎ element of the Coriolis and centripetal matrix can be defined as

𝑐 𝑘𝑗 = 1 2 ∑ [ 𝜕𝑚 𝑘𝑗 𝜕𝑞 𝑖 + 𝜕𝑚 𝑘𝑖 𝜕𝑞 𝑗 - 𝜕𝑚 𝑖𝑗 𝜕𝑞 𝑘 ] 𝑛 𝑖=1
𝑞i (𝑘, 𝑗 = 1,2, … . , 𝑛) Eq. 1-7

with 𝑚 𝑖𝑗 denotes an element of mass matrix with row index, 𝑖, and column index, 𝑗.

(iii)Gravity term

This term can be derived from the total potential energy of the biped system as follows.

𝑃 = ∑ 𝑚 𝑖 𝒈 𝑖 𝑇 𝒄 𝑖 𝑛 𝑖=1
Eq. 1-8

where 𝒈 𝑖 denotes the gravitational acceleration vector of link (𝑖); for plane system it is equal to [ 0 -9.81 0 ] 𝑇 . Every element of the gravitational term of Eq. 1-2 can be expressed as

𝑔 𝑘 = - 𝜕𝑃 𝜕𝒒 𝑘 (𝑘 = 1,2, … … . , 𝑛)
Eq. 1-9

(iv) The mapping matrix 𝐴 This mapping (coordinates transformation) matrix can be determined by the principle of the virtual work. The virtual work, 𝑊, of the generalized link torques acting on the biped system can be expressed in terms of generalized link coordinates as

𝛿𝑊 = 𝑇 1 𝛿𝑞 1 + 𝑇 2 𝛿𝑞 2 + ⋯ + 𝑇 𝑛 𝑞 𝛿𝑞 𝑛 𝑞 Eq. 1-10
with 𝑇 𝑖 denote the generalized torque of link (𝑖) associated with its generalized coordinate, 𝑞 𝑖 .

Eq. 1-10 can be re-written in a vector form as 𝛿𝑊 = 𝑻 𝑇 𝛿𝒒 Eq. 1-11

Let 𝜽 ∈ ℝ 𝑛 𝑞 be the generalized joint coordinates; thus, we can get a linear relationship between links and joint coordinates as follows.

𝛿𝒒 = 𝑨 -𝑇 𝛿𝜽 Eq. 1-12

Substituting Eq. 1-12 into Eq. 1-11, we have 𝛿𝑊 = 𝑻 𝑇 𝑨 -𝑇 𝛿𝜽 Eq. 1-13 𝑻 = 𝑨𝝉 Eq. 1-14 Remark 1-1. Alternatively, the matrix 𝑨 can be found simply according to the principle of freebody diagram as done in [START_REF] Vanderborght | Overview of the Lucy project: Dynamic stabilization of a biped powered by pneumatic artificial muscles[END_REF].

(v) Friction torques and other disturbance sources

In effect, the friction terms are complex and may be modeled approximately using the following form [START_REF] Zhu | H Virtual decomposition control: towards hyper degrees of freedom Berlin[END_REF] 𝝉 𝑓 𝑖 = 𝑐𝑜𝑢𝑙𝑜𝑚𝑏 𝑓𝑟𝑖𝑐𝑡𝑖𝑜𝑛 + 𝑣𝑖𝑠𝑐𝑜𝑢𝑠 𝑓𝑟𝑖𝑐𝑡𝑖𝑜𝑛 + 𝑆𝑡𝑟𝑖𝑏𝑒𝑐𝑘 𝑓𝑟𝑖𝑐𝑡𝑖𝑜𝑛

+ 𝑓𝑟𝑖𝑐𝑡𝑖𝑜𝑛 𝑜𝑓𝑓𝑠𝑒𝑡 𝑡𝑒𝑟𝑚 = 𝓀 𝑐𝑖 𝑠𝑖𝑔𝑛(𝜃 ̇𝑖) + 𝓀 𝑣𝑖 𝜃 ̇𝑖 + 𝓀 𝑠𝑖 𝑠𝑖𝑔𝑛(𝜃 ̇𝑖) exp(-( 𝜃 ̇𝑖 𝜂 𝑠𝑖 ⁄ )) + 𝓀 𝑜𝑖 (𝑖 = 1,2, … 𝑛)
Eq. 1-15

where 𝜃 ̇𝑖 represents the angular joint velocity of each link, 𝓀 𝑐𝑖 , 𝓀 𝑣𝑖 and 𝓀 𝑠𝑖 denote the Coulomb friction coefficient, viscous friction coefficient, and Stribeck friction effect respectively, and 𝓀 𝑜𝑖 is the friction offset term.

As we see from Eq. 1-15, friction has a local effect; the vector of friction torque is uncoupled.

In the light of the above formulae, Eq. 1-2 can be re-written as

∑ 𝑚 𝑘𝑗 (𝒒)𝑞̈𝑗 + 𝑛 𝑗=1 ∑ 𝑐 𝑘𝑗 (𝒒) 𝑞̇𝑗 + 𝑔 𝑘 (𝒒) + 𝑇 𝑓𝑘 = 𝑛 𝑗=1 𝑇 𝑘 (𝑘 = 1,2, … , 𝑛) Eq. 1-16
with 𝑇 𝑓𝑘 represents the friction torque affecting each link.

Remark 1-2.

There are several fundamental properties of the dynamic coefficient matrices, the mass matrix and the Coriolis and centripetal terms, which could be exploited in controller design of adaptive control. For more details on other properties, refer to [START_REF] Spong | Robot dynamics and control[END_REF].

Property 1-1. The mass matrix (𝑴) is symmetric and positive definite. This can be deduced from Eq. 1-4 and the property of the kinetic energy.

Property 1-2. The matrix 𝑴 ̇(𝒒) -𝟐𝑪(𝒒, 𝒒̇) is skew matrix, if the 𝑪(𝒒, 𝒒̇) matrix is described in terms of Christoffel symbols, Eq. 1-7. Proof of this property can be found in [START_REF] Spong | Robot dynamics and control[END_REF].

Property 1-3. The dynamic equations described in Eq. 1-2 are dependent linearly on certain parameters such as link masses, moment of inertia, friction coefficients etc.; consequently 𝑴(𝒒)𝒒̈+ 𝑪(𝒒, 𝒒̇)𝒒̇+ 𝒈(𝒒) + 𝑨𝝉 𝑓 = 𝒀(𝒒, 𝒒̇, 𝒒̈)𝜶 Eq. 1-17

where 𝒀(𝒒, 𝒒̇, 𝒒̈) ∈ ℝ 𝑛×𝑛 𝛼 is called the regressor matrix, a function of the known generalized coordinates and their first two derivatives, and 𝜶 ∈ ℝ 𝑛 𝛼 denotes the vector of unknown biped parameters.

Selection of 𝜶 is not unique, and it is difficult to find minimal set of these parameters [START_REF] Spong | Robot dynamics and control[END_REF]. Eq. 1-17 is very important to adaptive control.

The E-L equations of the first kind (the DSP)

As mentioned earlier, the biped mechanism constitutes a closed-chain with over-actuation during the DSP. Therefore, the Lagrange formulation of the 1 st kind, which can deals with constraints, is needed for dynamic modeling of the constrained biped. In such case, the motion equations are represented by redundant coordinates resulting in differential algebraic equations DAEs. The algebraic equations result from the constraints derived from the kinematics [START_REF] Tsai | Robot analysis: the mechanics of serial and parallel manipulators[END_REF]. The constraints can be easily incorporated into the main equations using Lagrange multipliers. The Lagrange equations of the biped robot during the DSP can be defined as where 𝜑 𝑗 denotes the constraint function of each closed loop, 𝑛 𝑐 is the number of these constraints, 𝜆 𝑗 is the Lagrange multipliers associated with each constraint. Here 𝑛 𝑞 is the number of redundant generalized coordinates and equal to the number DOFs (𝑛) of the biped systems plus the number of constraints(𝑛 𝑐 ).

Eq. 1-18 can be solved using two well-known techniques [START_REF] Pennestri | Coordinate reduction strategies in multibody dynamics: a review[END_REF]: the redundant coordinates-based techniques which are used mainly in commercial software such as MSC ADAMS, and the minimum coordinates-based techniques which could be, to some extent, suitable for control strategies and real-time applications. Many researchers have preferred the former technique due to its simplicity and ease of derivation at the expense of difficulties of numerical methods encountered in the solution [START_REF] Pennestri | Coordinate reduction strategies in multibody dynamics: a review[END_REF]. Consequently, this motivates the researchers to investigate the second technique which includes eliminating the constraint equations (Lagrange multipliers) from Eq. 1-18 to result in constraint-free differential equations [START_REF] Pennestri | Coordinate reduction strategies in multibody dynamics: a review[END_REF]. This can be implemented using one of the orthogonalization methods which are [START_REF] Pennestri | Coordinate reduction strategies in multibody dynamics: a review[END_REF]: coordinate partitioning method, zeroeigenvalue method, singular value decomposition (SVD), QR decomposition, Udwadia-Kabala formulation, PUTD method, and Schur decomposition. For more details, see [START_REF] Pennestri | Coordinate reduction strategies in multibody dynamics: a review[END_REF]. Although most researchers have written the matrices, 𝑭(𝒒), 𝑴 ̅ (𝒒), and 𝑪 ̅ (𝒒, 𝒒̇), in terms of the independent coordinates (𝒒 𝒊𝒏 ), these matrices still contain the dependent coordinates (𝒒). Therefore, we have expressed the mentioned matrices in terms of the last coordinates.

Solution of

Remark 1-5. The matrix 𝑭(𝒒) is not unique; the orthogonalization methods mentioned at the beginning of this subsection are used to get the matrix 𝑭(𝒒). Pennesri and Valentini [START_REF] Pennestri | Coordinate reduction strategies in multibody dynamics: a review[END_REF] simulated simple pendulum to compare the computational complexity of these orthogonalization methods. QR decomposition ranked best among the other methods. However, all these techniques could be computationally unsuitable to deal with the advanced adaptive control.

Remark 1-6. Eq. 1-31 has the same properties of that of Eq. 1-2 as follows [START_REF] Su | Adaptive control of robot manipulators under constrained motion[END_REF]. Since 𝑴 ̇-2𝑪 is skew-matrix according to Property 1-1, then 𝑵 is also skew-matrix.

Property 1-5. The orthogonality condition is satisfied by the matrix 𝑭(𝒒) such that Eq. 1-28 holds.

Proof. From Eq. 1-21 and substituting Eq. 1-27, we have

𝑱𝑭 𝒒̇𝑖 𝑛 = 𝟎 Eq. 1-37
Since 𝒒̇𝑖 𝑛 is linearly independent, then 𝑱𝑭 = 𝑭 𝑇 𝑱 𝑇 =0 Eq. 1-38 Property 1-6. If 𝑭(𝒒) is known, then the left hand side of Eq. 1-31 are linearly dependent on the unknown biped parameters (the same Property 1-3).

Continuous dynamic response

One of the inherent problems of legged locomotion (bipeds, quadrupeds, etc.) is the discontinuity at the transition instances due to: (i) impact events; these can be avoided by setting the foot velocity equal to zero at the instance of contact (see Chapters 5 and 6 of [START_REF] Al-Shuka | walking pattern generators and adaptive control of biped robot[END_REF]), and (ii) varying configurations of the biped from the SSP to the DSP and vice versa. As said previously in Section 1.1, the number of actuators is more than the DOFs of the biped during the constrained DSP. This means that there are infinity combinations of actuator torques to drive the biped systems as explained below.

One of the methods for determining the actuating torques and the ground reaction forces is the pseudo-inverse matrix as follows.

Eq As seen from Eq. 1-40, there is no guarantee that 𝝉 and 𝝀 have the same values at the start/end of the SSP due to this optimization solution. Therefore, the following assumption is proposed to resolve this dilemma. Accordingly, at the initial instance of DSP, 𝝀 = 𝟎, and the full ground reaction forces are supported by the rear foot, whereas, at the end of the DSP, the full support appears to be in the front foot with 𝝀 ̅ = 𝟎. On the other hand, because COG acceleration of the biped is nonlinear, the resulted ground reaction forces from Eq. 1-41 can generate nonlinear profile despite of multiplication of the latter equation with linear scaling function.

The modified recursive N-E formulation

Due to computational complexity inherent in the classical Lagrangian formulation, unless it is simplified, the researchers have resorted to the recursive N-E formulation for real time implementation. The philosophy of deriving N-E formulation is different from that of Lagrangian formulation. In the former, the translation/angular equations of motion of each link are derived sequentially using the D'Alembert principle. Due to the coupling effect between each neighbored links and appearance the translation equations of motion, the coupling force wrench appears in the derivation. Then a set of forward and backward recursive equations is used to determine the velocity and force wrenches respectively [START_REF] Spong | Robot dynamics and control[END_REF][START_REF] Fu | Robotics: control, sensing, vision, and intelligence USA[END_REF]. Various forms of N-E formulations for different robotic mechanism have been investigated; for details refer to [START_REF] Khalil | Dynamic modeling of robots using recursive Newton-Euler formulations[END_REF].

However, Fu et al. [START_REF] Fu | Robotics: control, sensing, vision, and intelligence USA[END_REF] have shown that the combined identification and control algorithms can be computed in 𝑂(𝑛 3 ) despite using recursive N-E formulation. Strictly speaking, dealing with advanced adaptive control techniques, the recursive N-E formulations could not be powerful; a modification is needed to satisfy the desired target. Zhu [START_REF] Zhu | H Virtual decomposition control: towards hyper degrees of freedom Berlin[END_REF] exploited the recursive nature of N-E equations to virtually decompose complex robotic systems into subsystems and to use the advanced adaptive techniques recursively. The derivation is exactly of that of recursive N-E formulation, but the difference is that the VDC derives the equation of motion of each link in terms of a frame attached at it first end rather than its COM.

Derivation of the dynamic equations

Now let us consider a fixed base serial-chain manipulator with revolute and prismatic joints. Thus, the links are numbered from 1 to 𝑛 𝑞 , where the base link is numbered as zero link. Fig. 123shows link (𝑖) where 𝑖 = 1,2, … . , 𝑛 𝑞 , is connected to other links via mechanical joints at its ends. This link has one driving cutting point associated with the frame {𝑇 𝑖+1 } and one driven cutting point associated with the frame {𝐵 𝑖 }. Thus, the joint (𝑖) has one driven cutting point associated with the frame {𝐵 𝑖 } and one driving cutting point associated with the frame {𝑇 𝑖 }.

Below, we will illustrate some remarks to make the derivation of dynamic equation of each subsystem (link, joint) accessible.

Eq. 1-47

Exploiting Remark 1-7 to yield

𝒇 ̅ * = 𝒇 ̅ 𝐵 𝑖 𝐵 𝑖 -𝑼 𝑇 𝑖+1 𝐵 𝑖 𝑼 𝐵 𝑖+1 𝑇 𝑖+1 𝒇 ̅ 𝐵 𝑖+1 = 𝒇 ̅ 𝐵 𝑖 -𝑼 𝐵 𝑖+1 𝐵 𝑖 𝒇 ̅ 𝐵 𝑖+1
Eq. 1-48

Remark 1-9. The velocity wrench of link (𝒊) can sequentially be determined by

𝒗 ̅ = 𝒛 𝑞i + 𝐵 𝑖 𝑼 𝐵 𝑖-1 𝐵 𝑖 𝑇 𝒗 ̅ 𝐵 𝑖-1
Eq. 1-49

with 𝒛 = [0 0 0 0 0 1] or [0 0 1 0 0 0] for revolute and prismatic joints respectively. Alternatively and simply, the velocity wrench can be calculated as

𝒗 ̅ = [ 𝑹 𝐼 𝒗 𝐵 𝑖 𝐵 𝑖 𝑹 𝐼 𝒘 𝐵 𝑖 𝐵 𝑖 ] 𝐵 𝑖
Eq. 1-50

with 𝒗 𝐵 𝑖 ∈ ℝ 3 and 𝒘 𝐵 𝑖 ∈ ℝ 3 are the absolute translational and angular velocity vector of frame 𝐵 𝑖 respect to the inertial frame {I}.

Remark 1-10. Concerning the target biped, the number of the generalized coordinates (𝒏 𝒒 ) is always equal to the number of links, e.g. both the number of links and the generalized coordinates are equal to six during the SSP. Consequently, we named the number of links as that of generalized coordinates.

(i) Dynamics of link subsystem By applying the D'Alembert principle to link (𝑖), with respect to the inertial frame about the COM of link (𝑖), we can get the following relations for the net forces 𝒇 ̅ 𝑖 * and the net moment 𝒎 ̅ 𝑖 * , as shown in Fig. 1234. where 𝒗 𝑖 ∈ ℝ 3×1 refers to the translation velocity vector of each link.

Putting Eq. 1-51 and Eq.

]

Eq. 1-61

(ii) Dynamics of revolute joint subsystem

There are two types of drive transmission systems for robotic joint systems. The first is the direct drive joints, in which the inertia of the motor is included in the corresponding links [START_REF] Zhu | H Virtual decomposition control: towards hyper degrees of freedom Berlin[END_REF][START_REF] Spong | Robot dynamics and control[END_REF], such that the dynamics of the joint is neglected. The second type of the system deals with a high gear transmission assuming that the inertial forces/torques act along the joint axis [START_REF] Zhu | H Virtual decomposition control: towards hyper degrees of freedom Berlin[END_REF]. In the latter case, the dynamic equation of the joint (𝑖) can be described as [START_REF] Zhu | H Virtual decomposition control: towards hyper degrees of freedom Berlin[END_REF] 𝐼 𝐽 𝑖 𝜃 ̈𝑖 + 𝓀 𝑐𝑖 𝑠𝑖𝑔𝑛(𝜃 ̇𝑖) + 𝓀 𝑣𝑖 𝜃 ̇𝑖 + 𝓀 𝑠𝑖 𝑠𝑖𝑔𝑛(𝜃 ̇𝑖) exp(-( 𝜃 ̇𝑖 𝜂 𝑠𝑖 ⁄ )) + 𝓀 𝑜𝑖 = 𝜏 𝑖 * (𝑖 = 1,2, … . , 𝑛 𝑗 ) Eq. 1-62

where 𝐼 𝐽 𝑖 represents the equivalent inertia of the joint (𝑖), 𝜃 ̈𝑖 denote the ith joint acceleration, 𝜏 𝑖 * represents the net torque applied to the joint (𝑖), and 𝑛 𝑗 denotes the number of joints. The net torque of the joint (𝑖) can be described as

𝜏 𝑖 * = 𝜏 𝑐𝑖 -𝒛 𝑇 𝒇 ̅ 𝐵 𝑖 (𝑖 = 1,2, … . , 𝑛 𝑗 ) Eq. 1-63
where 𝜏 𝑐𝑖 is the input control torque of the joint (𝑖) and the second term represents the output torque of the joint (𝑖) towards the link (𝑖).

Dynamics of the biped robot

Since the target biped is of a planar motion, simplifications appear in the dynamic equation of link (𝑖) derived in Eq. 1-58 to Eq. 1-61 via the following:

(i) Since the center of mass of link (𝑖) is located on the 𝑥 -𝑎𝑥𝑖𝑠 of frame {𝐵 𝑖 } with distance 𝑑 𝑖 from its its origin, 𝒓 𝐵 𝑖 𝐴 𝑖 𝐵 𝑖 = [𝑑 𝑖 , 0, 0] 𝑇 .

(ii) Since there is orientation in z-axis only, 𝒘 = [0, 0, 𝑤 𝑧 𝑖 ] 𝑇 𝐵 𝑖 .

(iii) 𝐼 𝐵 𝑖 = 𝑑𝑖𝑎𝑔(0, 0, 𝐼 𝑐 𝑖 ) .

(iv) Removal of the third to the fifth rows/columns of 𝑴 𝐵 𝑖 and 𝑪 𝐵 𝑖 .

Thus, the dynamics of planar biped robot can be expresses as in Eq. 1-58 with the following dynamic coefficient matrices

𝑴 𝐵 𝑖 = [ 𝑚 𝑖 0 0 0 𝑚 𝑖 𝑚 𝑖 𝑑 𝑖 0 𝑚 𝑖 𝑑 𝑖 𝐼 𝑐𝑖 + 𝑚 𝑖 𝑑 𝑖 2 ]
Eq. 1-64 Eq. 1-66

𝑪 𝐵 𝑖 = [ 0 -𝑚 𝑖 -

The SSP

As mentioned earlier, during this walking phase, the biped mechanism is an open chain mechanism with stance foot as fixed link; it should be in full contact with the ground. Therefore, the seven link-biped reduces to 6-link biped during its dynamic analysis. Three important points should be considered carefully when dealing with Eq. 1-56 that are:

(i) Determination of velocity wrench.

Solution of the dynamic Eq. 1-56 needs finding the velocity wrench which plays an important role in the adaptive control problem; it can be found as follows. See Fig. 1-5 for clear description of local frames. 𝑥 𝐵 2 = -0.2 , 𝑥Ḃ 2 = 0 Eq. 1-67

𝑦 𝐵 2 = 0 , 𝑦Ḃ 2 = 0 Eq. 1-68 𝒗 ̅ = [ 0 0 𝑞2 ] 𝐵 2
Eq. 1-69 Link (3) (stance thigh link):

𝑥 𝐵 3 = -0.2 + 𝑙 2 cos (𝑞 2 ) , 𝑥Ḃ 3 = -𝑙 2 sin (𝑞 2 )𝑞2
Eq. 1-70

𝑦 𝐵 3 = 𝑙 2 sin (𝑞 2 ) , 𝑦̇𝐵 3 = 𝑙 2 cos (𝑞 2 )𝑞̇2
Eq. 1-71

𝒗 ̅ = [ [ cos (𝑞 3 ) sin (𝑞 3 ) -sin (𝑞 3 ) cos (𝑞 3 ) ] [ 𝑥Ḃ 3 𝑦Ḃ 3 ] 𝑞3 ] 𝐵 3 = [ -𝑙 2 𝑞2 sin (𝑞 2 -𝑞 3 ) 𝑙 2 𝑞2 cos (𝑞 2 -𝑞 3 ) 𝑞3 ]
Eq. 1-72 Link (4) (Trunk):

𝑥 𝐵 4 = -0.2 + 𝑙 2 cos(𝑞 2 ) + 𝑙 3 cos(𝑞 3 ) , 𝑥̇𝐵 4 = -𝑙 2 sin(𝑞 2 ) 𝑞̇2 -𝑙 3 sin(𝑞 3 ) 𝑞̇3
Eq. 1-73 Eq. 1-84

𝑦 𝐵 4 = 𝑙 2 sin(𝑞 2 ) + 𝑙 3 sin(𝑞 3 ) , 𝑦̇𝐵 4 = 𝑙 2 cos(𝑞 2 ) 𝑞̇2 + 𝑙 3 cos(𝑞 3 ) 𝑞̇3

(ii) Resultant force wrench

To understand the force wrench distribution at the torso/leg interaction, see Fig. 12345. Thus, the following relations can be expressed for each link starting from the trunk.

Link (4) (trunk):

𝒇 ̅ * = 𝒇 ̅ + 𝒇 ̅ 𝐵 𝑇2 𝐵 𝑇1 𝐵 4
Eq. 1-85 with notations shown in Fig. 12345. Eq. 𝑐𝑜𝑠(𝑞 3 -𝑞 4 ) 0 -𝑙 3 sin(𝑞 3 -𝑞 4 ) 𝑙 3 cos(𝑞 3 -𝑞 4 ) 1

]

Eq. 1-93 Link (2)(stance shank):

𝒇 ̅ * 𝐵 2 = 𝒇 ̅ -𝑼 𝐵 3 𝒇 ̅ 𝐵 3 𝐵 2 𝐵 2 Eq. 1-94 with 𝑼 𝐵 3 𝐵 2 = [ 𝑐𝑜𝑠(𝑞 3 -𝑞 2 ) -𝑠𝑖𝑛(𝑞 3 -𝑞 2 ) 0 𝑠𝑖𝑛(𝑞 3 -𝑞 2 ) 𝑐𝑜𝑠(𝑞 3 -𝑞 2 ) 0 𝑙 2 sin(𝑞 3 -𝑞 2 ) 𝑙 2 cos(𝑞 3 -𝑞 2 ) 1 ]
Eq. 1-95

We note that we have 6 equations for six links(Eq. 1-85, Eq. 1-86, Eq. 1-88, Eq. 1-90, Eq. 1-90 and so on.

(iii) Actuating torques

To simplify the analysis, let us assume temporarily that the target biped has direct drive joint systems (the dynamics of the joint could be included with the corresponding links [START_REF] Zhu | H Virtual decomposition control: towards hyper degrees of freedom Berlin[END_REF][START_REF] Spong | Robot dynamics and control[END_REF]) the left hand side of Eq. 1-63 is equal to zero. Consequently, the actuating torques can be calculated from the coupling effect of the neighbored link according to Eq. 1-63.

Left shank, 𝜏 𝑐1 = 𝒛 𝑇 𝒇 ̅ 𝐵 2

Eq. 1-96

Left knee, 𝜏 𝑐2 = 𝒛 𝑇 𝒇 ̅ 𝐵 3

Eq. 1-97

Left thigh/torso interaction, 𝜏 𝑐3 = 𝒛 𝑇 𝒇 ̅ 𝐵 𝑇 1

Eq. 1-98

Right thigh/torso interaction, 𝜏 𝑐4 = 𝒛 𝑇 𝒇 ̅ 𝐵 𝑇 2

Eq. 1-99

Right knee, 𝜏 𝑐5 = 𝒛 𝑇 𝒇 ̅ 𝐵 5

Eq.

1-100

Right shank, 𝜏 𝑐6 = 𝒛 𝑇 𝒇 ̅ 𝐵 6

Eq. 1-101

The DSP1

As mentioned earlier, the biped in this walking sub-phase, DSP1, has six actuators with 4 DOFs; therefore, two redundant actuators compromise the over-actuation problem. In the following, the details of velocity and force wrenches as well as determining the redundant actuating torques are investigated.

 Velocity wrench. It has exactly the same relations described in previous subsection, Eq. 1-67 to Eq. 1-84 with replacing the word (swing) by (front), the word (stance) to (rear), and 𝑙 7𝑎 by (𝑙 7 -𝑙 7𝑎 ) for the last link.

 Force wrench. It has also the same force wrenches showed in Eq. 1-85 to Eq. 1-95.

 Actuating torques. We have three significant problems resulting from the variable configurations of the biped which include: (a) redundancy of the actuators, (b) the passive joint on the front foot (see Fig. 123456) which enforces the torque to be 𝜏 𝑐 𝑓 = 0, and (c) the discontinuity of the actuating torques. Four solutions are considered below with focus on solutions 3 and 4 which can be ranked best among the rest. Then, describing the actuating torques in terms of the design variables (𝜼).

(ii) Procedure 2-direct optimization of the torso/leg force wrench Instead of releasing internal force wrench, the actuating torques can directly be expressed in terms of 𝒇 ̅ 𝐵 𝑇 1 or 𝒇 ̅ 𝐵 𝑇 2 . Thus, we can get the same equations above but in terms of 𝒇 ̅ 𝐵 𝑇 1 as follows.

𝒇 ̅

𝐵 𝑇 1 = -(𝑩 1 𝑇 𝓦𝑩 1 ) -1 𝑩 1 𝑇 𝓦𝒃 1 -(𝑩 1 𝑇 𝓦𝑩 1 ) -1 𝓪 𝟏 ( 𝑎 2 -𝓪 𝟏 𝑇 (𝑩 1 𝑇 𝓦𝑩 1 ) -1 𝑩 1 𝑇 𝓦𝒃 1 𝓪 𝟏 𝑇 (𝑩 1 𝑇 𝓦𝑩 1 ) -1 𝓪 𝟏 )
Eq.1-111

and completing the same steps as of the procedure 1. However, the discontinuity problem has not been resolved in the above two procedures.

(iii) Procedure 3-tracking desired ground reaction forces.

Considering Assumption 1-4 and assuming the desired reaction force, see Eq. 1-41, as a constraint to yield

𝒇 ̅ = [ 𝝀 𝜏 𝑐 𝑓 ] 𝐵 7 = 𝑩 ̆𝒇 ̅ 𝐵 𝑇 1 + 𝒃 ̆ Eq.1-112
with 𝑩 ̆∈ ℝ 3×3 and 𝒃 ̆∈ ℝ 3 .

The left hand side of Eq.1-112 is known from the desired walking trajectories, so the problems of over-actuating and discontinuity are solved using the last equation without need of optimization.

Remark 1-11. If the number of constraints is equal to the design variables, no optimization of the system is necessary because the solution of equality constraints are the only candidates for the optimum design [START_REF] Arora | Introduction to optimum design[END_REF].

(iv) Procedure 4-tracking desired ground reaction forces with optimization

In this procedure, we will come back to procedure 1 representing the trunk/leg interaction force wrench in terms of the internal force wrench and 𝛼 parameter. Thus, constraint Eq. where ℴ ̅ denotes the candidate optimal solution. Because of the bounded limits of ℴ, the following procedure is proposed:

If 0 ≤ ℴ ̅ ≤ 1, then ℴ = ℴ ̅.

If ℴ ̅ ≥ 1, then ℴ = 1.

If ℴ ̅ ≤ 0, then ℴ = 0.

Eq.1-115

After finding the internal force wrench and 𝛼 parameter, it is easy to find the actuating torques in a similar way described previously (see procedure 1).

The DSP2

In this walking sub-phase, the biped robot has two redundant actuators with a slightly different configuration of that of DSP1 (see Fig. 1234567). Here, the front foot is flat on the ground with negligible dynamics while the rear foot rotates about its front tip. In similar manner to DSP1, the velocity wrench of each link can be determined as made in Subsection 1.2.2.2.1, whereas, the force wrench can be calculated successfully using procedure 1 described in the previous subsection.

 Velocity wrench. It is exactly determined in the same manner used in the SSP and DSP1; therefore, we will list only the final results of each link without details.  Actuating torques. Procedure 3 or 4 can be used to solve the problem of over-actuation and discontinuity during this walking sub-phase. After finding the internal force wrench and 𝛼 parameter as described in Eq. 1-113 to Eq.1-115, the actuating torques can be determined by Eq. 1-96 to Eq. 1-101.

Conclusions and future work

In this work, we have modeled two selected walking patterns of biped robot using Lagrangian and virtual decomposition-based N-E formulation. The problem of discontinuity is solved using linear transition ground reaction forces without impact-contact event. Lagrangian formulation, unless simplified, could require more computational complexity than that of virtual decomposition-based N-E formulation.

  Both DSP1 and DSP2 have six links as shown in Fig. 1-2 (a) and (b) respectively.

Fig. 1 - 1 :Fig. 1 - 2 :

 1112 Fig. 1-1: Walking pattern 1 with description of generalized coordinates. (a) Intermediate configuration of biped locomotion during the SSP. (b) Intermediate configuration of biped locomotion during the DSP.

  𝑴(𝒒)𝒒̈+ 𝑪(𝒒, 𝒒̇)𝒒̇+ 𝒈(𝒒) = 𝑨(𝝉 -𝝉 𝑓 )or simply,𝑴(𝒒)𝒒̈+ 𝑪(𝒒, 𝒒̇)𝒒̇+ 𝒈(𝒒) + 𝑨𝝉 𝑓 = 𝑨𝝉Eq. 1-2

  1, 2, … … . . , 𝑛 𝑞 ) Eq. 1-18

Property 1 - 4 .

 14 Let 𝑫(𝒒) = 𝑭(𝒒) 𝑻 𝑴 ̅ (𝒒), 𝑯(𝒒, 𝒒̇) = 𝑭(𝒒) 𝑻 𝑪 ̅ (𝒒, 𝒒̇), then the matrix 𝑵(𝒒, 𝒒̇) = 𝑫 ̇(𝒒, 𝒒̇) -𝟐𝑯(𝒒, 𝒒̇) is skew-matrix. Proof. Let 𝑵 = 𝑫 ̇-2𝑯 Eq. 1-35 By substituting Eq. 1-32 into Eq. 1-35 we get 𝑵 = 𝑭 ̇𝑇𝑴𝑭 + 𝑭 𝑇 𝑴 ̇𝑭 + 𝑭 𝑇 𝑴𝑭 ̇-2𝑭 𝑇 𝑴𝑭 ̇-𝑭 𝑇 𝑪𝑭 = 𝑭 𝑇 (𝑴 ̇-2𝑪)𝑭 + 𝑭 ̇𝑇𝑴𝑭 -𝑭 𝑇 𝑴𝑭 ̇= 𝑭 𝑇 (𝑴 ̇-2𝑪)𝑭 Eq. 1-36

. 1 -

 1 19 can be re-arranged to yield 𝑴(𝒒)𝒒̈+ 𝑪(𝒒, 𝒒̇)𝒒̇+ 𝒈(𝒒) + 𝑨𝝉 𝑓 = [𝑨 𝑱(𝒒) 𝑇 ] One of the possible solutions to get the actuating torques and Lagrange multipliers are [ 𝝉 𝝀 ] = [𝑨 𝑱 𝑇 ] ≠ (𝑴(𝒒)𝒒̈+ 𝑪(𝒒, 𝒒̇)𝒒̇+ 𝒈(𝒒) + 𝑨𝝉 𝑓 ) Eq. 1-40 where the notation [. ] ≠ denotes the pseudo-inverse of the referred matrix.

Assumption 1 - 4 .

 14 Because the biped robot does not have a unique solution during the DSP, a linear transition function could be proposed for the ground reaction forces [Alb12]. Thus, for the front foot 𝝀 = ( 𝑡 -𝑡 𝑠 𝑡 𝑑 -𝑡 𝑠 ) 𝑚 𝐺 (𝒄̈𝐺 + [0, 𝑔, 0] 𝑇 ) Eq. 1-41 where 𝑡, 𝑡 𝑠 and 𝑡 𝑑 are time parameter, the time of SSP, and the DSP time. Meanwhile, the ground reaction forces, 𝝀 ̅ , of the rear foot are 𝝀 ̅ = 𝑚 𝐺 (𝒄̈𝐺 + [0, 𝑔, 0] 𝑇 ) -𝝀 Eq. 1-42

Fig. 1 - 3 :

 13 Fig. 1-3: Virtual decomposition of a serial-chain manipulator Remark 1-7. [Zhu10]. The matrix of force wrench transformation, 𝑼 𝑨 ∈ ℝ 𝟔×𝟔 𝑩

Fig. 1 - 4 :

 14 Fig. 1-4: Virtual decomposition of link (𝑖) with description of the net force wrench. (a) The freebody diagram of the link (𝑖) with force wrench at its ends and gravity effect. (b) The net force wrench applied at the COG of link (𝑖) with gravity effect. (c) The net force wrench applied at one end of link (𝑖) with gravity effect.

Fig. 1 - 5 :

 15 Fig. 1-5: Biped robot during the SSP with description of assumed local frames Link (1) (stance foot link): It is assumed fixed link with negligible dynamics.

  cos(𝑞 5 -𝑞 4 ) -𝑙 5 sin(𝑞 5 -𝑞 4 ) 1 ]

Fig. 1 - 6 :

 16 Fig. 1-6: Biped robot during the DSP1 (i) Procedure 1-releasing and optimizing the internal forces

Link ( 1 Fig. 1 - 7 :

 117 Fig. 1-7: Biped robot during the DSP2

1-34 Remark 1-4.

  𝑛 𝑐 is the constraint vector, and 𝑱 ∈ ℝ 𝑛 𝑐 ×𝑛 𝑞 = To reduce the dimension size of Eq. 1-19 (to eliminate 𝝀), a relationship between the redundant generalized coordinates (𝒒 )and the independent coordinates (𝒒 𝑖𝑛 ∈ ℝ 𝑛 ) should be found. In this report, the coordinate partitioning is used for size reduction of the equation of motion[START_REF] Mitobe | Control of a biped walking robot during the double support phase[END_REF].The matrix 𝑭 ∈ ℝ 𝑛 𝑞 ×𝑛 plays an important role in eliminating 𝝀 ; the following orthogonality condition holds

	with 𝝋 ∈ ℝ 𝜕𝝋(𝒒) 𝜕𝒒	denotes the Jacobian matrix.
	Remark 1-3. The coefficient dynamic matrices (mass matrix, Coriolis and centripetal matrix etc.) of Eq. 1-19 could be determined by the same mathematical formulae defined in the open-chain 𝑱(𝒒)𝑭(𝒒) = 𝟎 Eq. 1-28
	mechanism as in Eq. 1-6 to Eq. 1-8 Differentiating Eq. 1-27 to obtain
	𝒒̈= 𝑭(𝒒) 𝒒̈𝑖 𝑛 + 𝑭 ̇(𝒒, 𝒒̇) 𝒒̇𝑖 𝑛		Eq. 1-29
	Substituting Eq. 1-27 and Eq. 1-29 into Eq. 1-19 to get
	Twice differentiating Eq. 1-20 can result in 𝑴(𝒒)(𝑭(𝒒)𝒒 ̈𝑖𝑛 + 𝑭 ̇(𝒒, 𝒒̇)𝒒 ̇𝑖𝑛 ) + 𝑪(𝒒, 𝒒̇)𝑭(𝒒)𝒒 ̇𝑖𝑛 + 𝒈(𝒒) + 𝑨𝝉 𝑓 = 𝑨𝝉 + 𝑱(𝒒) 𝑇 𝝀	Eq. 1-30
	𝑱(𝒒)𝒒̇= 𝟎				Eq. 1-21
	[	𝑱(𝒒) 𝑱 𝒓 (𝒒)	] 𝒒̇= [	𝟎 𝒒̇𝑖 𝑛	]		Eq. 1-26
	Eq. 1-18 results in Thus, it is possible to get the following important relations
	𝑴(𝒒)𝒒̈+ 𝑪(𝒒, 𝒒̇)𝒒̇+ 𝒈(𝒒) + 𝑨𝝉 𝑓 = 𝑨𝝉 + 𝑱 𝑇 𝝀 𝒒̇= [ 𝑱(𝒒) 𝑱 𝒓 (𝒒) ] -1 [ 𝟎 𝒒̇𝑖 𝑛 ] = [ 𝑭 ̆𝑭] [ 𝟎 𝒒̇𝑖 𝑛 ] = 𝑭(𝒒)𝒒̇𝑖 𝑛		Eq. 1-19 Eq. 1-27
	𝝋(𝒒) = 𝟎				Eq. 1-20

𝑱(𝒒)𝒒̈+ 𝑱 (𝒒, 𝒒̇)𝒒̇= 𝟎

Eq. 1-22

Due to the redundancy of coordinates in Eq. 1-19, it is possible to express the dependent generalized coordinates in terms of the independent ones as in Eq. 1-23. 𝒒 𝑖𝑛 = 𝝋 ̆(𝒒)𝒒 Eq. 1-23 Twice differentiating Eq. 1-23 yields 𝒒̇𝑖 𝑛 = 𝑱 𝒓 (𝒒)𝒒̇ Eq. 1-24 with 𝑱 𝒓 (𝒒) ∈ ℝ 𝑛×𝑛 𝑞 = 𝜕𝝋 ̆(𝒒) 𝜕𝒒 𝒒̈𝑖 𝑛 = 𝑱 𝒓 (𝒒)𝒒̈+ 𝑱̇𝒓(𝒒)𝒒̇ Eq. 1-25 Blocking together Eq. 1-21 and Eq. 1-24 to get Alternatively, Eq. 1-30 can be re-written as 𝑴 ̅ (𝒒)𝒒̈𝑖 𝑛 + 𝑪 ̅ (𝒒, 𝒒̇)𝒒̇𝑖 𝑛 + 𝒈(𝒒) + 𝑨𝝉 𝑓 = 𝑨𝝉 + 𝑱(𝒒) 𝑇 𝝀 Eq. 1-31 with 𝑴 ̅ (𝒒) = 𝑴(𝒒)𝑭(𝒒), 𝑪 ̅ (𝒒, 𝒒̇) = 𝑴(𝒒)𝑭 ̇(𝒒, 𝒒̇) + 𝑪(𝒒, 𝒒̇)𝑭(𝒒) Eq. 1-32 Using Eq. 1-27 and Eq. 1-30 can yield 𝝀 = 𝑭 ̆(𝒒) 𝑇 (𝑴 ̅ (𝒒)𝒒̈𝑖 𝑛 + 𝑪 ̅ (𝒒, 𝒒̇)𝒒̇𝑖 𝑛 + 𝒈(𝒒) + 𝑨𝝉 𝑓 -𝑨𝝉) Eq. 1-33 Exploiting Eq. 1-28 and pre-multiplying Eq. 1-31 by 𝑭(𝒒) 𝑇 to obtain 𝑭(𝒒) 𝑇 𝑴 ̅ (𝒒) 𝒒̈𝑖 𝑛 + 𝑭(𝒒) 𝑇 𝑪 ̅ (𝒒, 𝒒̇)𝒒̇𝑖 𝑛 + 𝑭(𝒒) 𝑇 𝒈(𝒒) + 𝑭(𝒒) 𝑇 𝑨𝝉 𝑓 = 𝑭(𝒒) 𝑇 𝑨𝝉 Eq.

  In similar manner, the velocity wrench can be represented in terms of the frame {𝐵 𝑖 } as where 𝑰 𝐵 𝑖 represents the configuration-free inertia tensor expressed in frame {𝐵 𝑖 }

	𝒗 ̅ 𝑖 = [ 𝒘 𝑖 𝒗 𝑖	] = [	𝑹 𝐴 𝑖 𝟎 3×3 𝐼		𝟎 3×3 𝑹 𝐴 𝑖 𝐼	] 𝑼 𝐵 𝑖	𝐴 𝑖 𝑇	𝐵 𝑖	𝒗 ̅ 𝑖	Eq. 1-56
	Differentiating Eq. 1-56 results in
	𝒗 ̅ ̇𝑖 = [	𝒗̇𝑖 𝒘̇𝑖	] = [	(𝒘 𝑖 ×) 𝑹 𝐴 𝑖 𝐼 𝟎 3×3	𝟎 3×3 (𝒘 𝑖 ×) 𝑹 𝐴 𝑖 𝐼	] 𝑼 𝐴 𝑖 𝑇 𝐵 𝑖	𝐵 𝑖	𝒗 ̅ 𝑖	+ [	𝑹 𝐴 𝑖 𝟎 3×3 𝐼	𝟎 3×3 𝑹 𝐴 𝑖 𝐼	] 𝑼 𝐴 𝑖 𝑇 𝐵 𝑖	𝑑 𝑑𝑡	( 𝒗 ̅ 𝑖 𝐵 𝑖 )	Eq. 1-57
	Substituting Eq. 1-55 and Eq. 1-57 into Eq. 1-54 results in
	𝑴 𝐵 𝑖 𝒗 ̅ ̇𝑖 + 𝑪 𝐵 𝑖 𝐵 𝑖	( 𝒘 𝐵 𝑖	𝑖 ) 𝒗 ̅ 𝑖 𝐵 𝑖	+ 𝒈 𝐵 𝑖 = 𝒇 ̅ * 𝐵 𝑖	Eq. 1-58
	with
	𝑴 𝐵 𝑖 = [ 𝑚 𝑖 ( 𝒓 𝑚 𝑖 𝑰 3 𝐵 𝑖 𝐵 𝑖 𝐴 𝑖 ×) 𝑰 𝐵 𝑖 -𝑚 𝑖 ( 𝒓 -𝑚 𝑖 ( 𝒓 𝐵 𝑖 𝐵 𝑖 𝐴 𝑖 ×) 𝐵 𝑖 𝐵 𝑖 𝐴 𝑖 ) 2 ]	Eq. 1-59
	𝑪 𝐵 𝑖 ( 𝒘 𝐵 𝑖	𝑖 ) = [	𝑐 11 𝑐 12 𝑐 21 𝑐 22	]	Eq. 1-60
	[ with 𝑚 𝑖 𝑰 3 𝟎 3×3 𝟎 3×3 𝓘 𝒊 (𝑡) 𝑐 11 = 𝑚 𝑖 ( 𝒘 𝑖 ×) ] [ 𝒘̇𝑖 𝒗̇𝑖 𝐵 𝑖 , 𝑐 12 = -𝑚 𝑖 ( 𝒘 𝑖 ×) ] + [ 𝑚 𝑖 𝒈 (𝒘 𝑖 ×)𝓘 𝒊 (𝑡)𝒘 𝑖 𝐵 𝑖 ( 𝒓 ] = [ 𝐵 𝑖 𝐵 𝑖 𝐴 𝑖 ×), 𝑐 21 = 𝑚 𝑖 ( 𝒓 * 𝒇 𝑖 𝒎 𝑖 * ] 𝐵 𝑖	𝐵 𝑖 𝐴 𝑖 ×)( 𝒘 𝑖 ×) 𝐵 𝑖	Eq. 1-53
	or [ 𝑐 22 = ( 𝒘 𝑖 ×) 𝑚 𝑖 𝑰 𝟑 𝟎 3×3 𝟎 3×3 𝓘 𝒊 (𝑡) 𝐵 𝑖 𝑰 𝐵 𝑖 + 𝑰 𝐵 𝑖 ( 𝒘 𝑖 ×) ] 𝒗 ̅ ̇𝑖 + [ 𝐵 𝑖 𝑚 𝑖 𝒈 (𝒘 𝑖 ×) 𝓘 𝒊 (𝑡)𝒘 𝑖 -𝑚 𝑖 ( 𝒓 𝐵 𝑖 ] = 𝒇 ̅ 𝑖 * with 𝑰 3 is 3 × 3 identity matrix. 𝐵 𝑖 𝑚 𝑖 𝑹 𝐼 𝒈 𝒈 𝐵 𝑖 = [ 𝑚 𝑖 ( 𝒓 𝐵 𝑖 𝐵 𝑖 𝐴 𝑖 ×) 𝑹 𝐼 𝒈 𝐵 𝑖	𝐵 𝑖 𝐴 𝑖 ×)( 𝒘 𝑖 ×)( 𝒓 𝐵 𝑖 𝐵 𝑖	𝐵 𝑖 𝐴 𝑖 ×)	Eq. 1-54
	Exploiting Remark 1-8, the net force wrench on the right hand side of Eq. 1-54 can be expressed
	(transformed) in terms of the frame {𝐵 𝑖 } as follows.
	𝐵 𝑖	𝒇 ̅ * = 𝑼 𝐵 𝑖	𝐴 𝑖	𝐴 𝑖	𝒇 ̅ 𝑖 *	= 𝑼 𝐵 𝑖	𝐴 𝑖 [ 𝟎 3×3 𝑹 𝐼 𝐴 𝑖	𝟎 3×3 𝑹 𝐼 𝐴 𝑖	] [	* 𝒇 𝑖 𝒎 𝑖 * ]	Eq. 1-55

1-52 into block matrix to deal with velocity and force wrenches

  ). Because the swing foot does not have force wrench at the frame {𝐵 7 }, so 𝒇 ̅ 𝐵 7 = 𝟎. Thus, 𝒇 ̅ 𝐵 6 can recursively be calculated from Eq.

	1-92,
	Eq. 1-94), with 7 unknowns ( 𝒇 ̅ 𝐵 𝑇 1 , 𝒇 ̅ 𝐵 𝑇 2 , 𝒇 ̅ 𝐵 5 , 𝒇 ̅ 𝐵 6 , 𝒇 ̅ 𝐵 7 , 𝒇 ̅ 𝐵 3 , 𝒇 ̅ 𝐵 2

  with 𝑩 2 ∈ ℝ 3×3 , 𝒃 2 ∈ ℝ 3 and 𝒃 3 ∈ ℝ 3 . Re-arranging Eq. 1-113 and using the pseudo-inverse definition

	[	𝜼 ℴ ̅	] = [𝑩 2 𝒃 2 ] ≠ ( 𝒇 ̅ 𝐵 7 -𝒃 3 )	Eq.1-114
							1-112 can be
	expressed as follows
	𝐵 7	𝒇 ̅ = [	𝝀 𝜏 𝑐 𝑓	]	= 𝑩 2 𝜼 + 𝒃 2 ℴ + 𝒃 3	Eq. 1-113

𝝉 = 𝑩 1 𝜼 + 𝒃 1

Eq.

1-104

with the constraint of passive joint

Eq.

1-105

By defining the objective function

Eq.

1-106

where 𝓦 ∈ ℝ 6×6 is a symmetric weighting matrix; it is assumed as an identity matrix in our solution.

Substituting Eq. 1-104 into Eq. 1-106 and incorporating the constraint of Eq. 4-105 to get

Eq.

1-107

Differentiating Eq. 1-107 with respect to 𝜼 and setting it to zero

1-108

Substituting Eq. 1-108 into Eq. 1-105 to yield

Eq.

1-109

Substituting Eq. 1-109 into Eq. 1-108 to get the internal force wrench

Eq.

1-110

Thus, the force wrench at the torso can be determined from Eq. 1-85, and sequentially finding the rest force wrenches and the required torques via the following Eq. 1-86, Eq. 1-88, Eq. 1-90,Eq. 1-92,Eq. 1-94, Eq. 1-96 to Eq. 1-101. The disadvantages of this procedure is that ℴ is a free parameter; it has not been considered as design variable, and there is also no guarantee to satisfy continuous dynamic response related to actuating torques.