Dynamics and Control of Humanoid Robots: A Geometrical Approach

his paper reviews modern geometrical dynamics and control of humanoid robots. This general Lagrangian and Hamiltonian formalism starts with a proper definition of humanoid's configuration manifold, which is a set of all robot's active joint angles. Based on the `covariant force law', the general humanoid's dynamics and control are developed. Autonomous Lagrangian dynamics is formulated on the associated `humanoid velocity phase space', while autonomous Hamiltonian dynamics is formulated on the associated `humanoid momentum phase space'. Neural-like hierarchical humanoid control naturally follows this geometrical prescription. This purely rotational and autonomous dynamics and control is then generalized into the framework of modern non-autonomous biomechanics, defining the Hamiltonian fitness function. The paper concludes with several simulation examples. Keywords: Humanoid robots, Lagrangian and Hamiltonian formalisms, neural-like humanoid control, time-dependent biodynamics


Introduction
Humanoid robots, being the future of robotic science, are becoming more and more human-like in all aspects of their functioning. Both human biodynamics and humanoid robotics are governed by Newtonian dynamical laws and reflex-like nonlinear controls [1,10,11,8,2].
Although motion of humanoid robots increasingly resembles human motion, we still need to emphasize that human joints are (and will probably always remain) significantly more flexible than humanoid robot joints. Each joint of a humanoid robot consists of a pair of coupled segments with only Eulerian rotational degrees of freedom. Each human synovial joint, on the other hand, not only exhibits gross rotational movement (roll, pitch and yaw) but is also capable of exhibiting some hidden and restricted translations along (X, Y, Z) axes. For example, in the knee joint, patella (knee cap) moves for about 7-10 cm from maximal extension to maximal flexion. It is well-known that translational amplitudes in the shoulder joint are even greater. In other words, within the realm of rigid body mechanics, a segment of a human arm or leg is not properly represented as a rigid body fixed at a certain point, but rather as a rigid body hanging on rope-like ligaments. More generally, the whole skeleton mechanically represents a system of flexibly coupled rigid bodies, technically an anthropomorphic topological product of SE(3)-groups. This implies more complex kinematics, dynamics and control than in the case of humanoid robots [3].
This paper reviews modern geometrical approaches to humanoid robot's dynamics and control. It is largely based on authors' own research in closely related fields of human biodynamics, biomechanics and humanoid robotics. This general approach starts with a proper definition of humanoid's configuration manifold M , which is a set of all active degrees-of-freedom (DOF). Based on the covariant force law, the general humanoid's dynamics with large number of DOF is developed. The tangent bundle of the manifold M (called the velocity phase space) is the stage for autonomous Lagrangian formulation of humanoid's dynamics, while the cotangent bundle of the manifold M (called the momentum phase space) is the stage for autonomous Hamiltonian formulation of humanoid's dynamics. This purely rotational and autonomous robot dynamics is then generalized along the two main lines of modern non-autonomous biomechanics: (i) more flexible joints, and (ii) time-dependent energy function (with energy sources and sinks).
In contrast to our previously published papers, the present article provides full technical details of both autonomous and non-autonomous (time-dependent) biodynamics and robotics, including the new neuro-muscular fitness dynamics. This thorough theoretical background would provide an interested reader with superb capability to develop their own non-autonomous humanoid simulator.

Configuration Manifold and the Covariant Force Law
Representation of an ideal humanoid-robot motion (with human-like spine, see Figure 1) is rigorously defined in terms of rotational constrained SO(3)-groups of motion [5,6,7,9] in all main robot joints. Therefore, the configuration manifold M rob for humanoid dynamics is defined as a topological product of all included SO(3) groups, M rob = i SO(3) i . Consequently, the natural stage for autonomous Lagrangian dynamics of robot motion is the tangent bundle T M rob . 1 System's Lagrangian (energy function) is a natural energy function on the tangent bundle [8]. Similarly, the natural stage for autonomous Hamiltonian dynamics of robot motion is the cotangent bundle T * M rob . 2 The Hamiltonian is a natural energy function on the tangent bundle [10,1,11].
More precisely, the three-axial SO(3)−group of humanoid-robot joint rotations depends on three parameters, Euler joint angles q i = (ϕ, ψ, θ), defining the rotations about the Cartesian coordinate triedar (x, y, z) placed at the joint pivot point. Each of the Euler angles are defined in the constrained range (−π, π), so the joint group space is a constrained sphere of radius π [2,9].
be the group of rotations in R 3 . It is a Lie group and dim(G) = 3. Let us isolate its one-parameter joint subgroups, i.e., consider the three operators of the finite joint rotations R ϕ , R ψ , R θ ∈ SO(3), given by corresponding respectively to rotations about x−axis by an angle ϕ, about y−axis by an angle ψ, and about z−axis by an angle θ.
However, the order of these matrix products matters: different order products give different results, 1 Recall that in multibody mechanics, to each n−dimensional (nD) configuration manifold M there is associated its 2nD velocity phase-space manifold, denoted by T M and called the tangent bundle of M . The original smooth manifold M is called the base of T M . There is an onto map π : T M → M , called the projection. Above each point x ∈ M there is a tangent space TxM = π −1 (x) to M at x, which is called a fibre. The fibre TxM ⊂ T M is the subset of T M , such that the total tangent bundle, T M = m∈M TxM , is a disjoint union of tangent spaces TxM to M for all points x ∈ M . From dynamical perspective, the most important quantity in the tangent bundle concept is the smooth map v : M → T M , which is an inverse to the projection π, i.e, π • v = Id M , π(v(x)) = x. It is called the velocity vector-field. Its graph (x, v(x)) represents the cross-section of the tangent bundle T M . This explains the dynamical term velocity phase-space, given to the tangent bundle T M of the manifold M . The tangent bundle is where tangent vectors live, and is itself a smooth manifold. Vector-fields are cross-sections of the tangent bundle. 2 Recall that in multibody mechanics, a dual notion to the tangent space TmM to a smooth configuration manifold M at a point m is its cotangent space T * m M at the same point m. Similarly to the tangent bundle, for a smooth manifold M of dimension n, its cotangent bundle T * M is the disjoint union of all its cotangent spaces T * m M at all points m ∈ M , i.e., T * M = as the matrix product is noncommutative product. 4 This is the reason why Hamilton's quaternions 5 are today commonly used to parameterize the SO(3)−group, especially in the field of 3D computer graphics.
The autonomous humanoid dynamics (both Lagrangian and Hamiltonian), is based on the postulate of conservation of total mechanical energy. It can be derived from the covariant force law [2,9], which in 'plain English' states: Force 1-form = Mass distribution × Acceleration vector-field, 4 The one-parameter rotations Rϕ, R ψ , R θ define curves in SO(3) starting from Their derivatives in ϕ = 0, ψ = 0 and θ = 0 belong to the associated tangent Lie algebra so (3). That is, the corresponding infinitesimal generators of joint rotations -joint angular velocities vϕ, v ψ , v θ ∈ so(3) -are respectively given by Moreover, the elements are linearly independent and so Then we have the following identities: . The exponential map exp : so(3) → SO(3) is given by Rodrigues relation The the dual, cotangent Lie algebra so(3) * , includes the three joint angular momenta pϕ, p ψ , p θ ∈ so(3) * , derived from the joint velocities v by multiplying them with corresponding moments of inertia. 5 Recall that the set of Hamilton's quaternions H represents an extension of the set of complex numbers C. We can compute a rotation about the unit vector, u by an angle θ. The quaternion q that computes this rotation is and formally reads (using Einstein's summation convention over repeated indices): Here, the force 1-form F i = F i (t, q, p) = F ′ i (t, q,q), (i = 1, ..., n) denotes any type of actuator torques; m ij is the material (mass-inertia) metric tensor, which gives the total mass distribution of the robot (including all segmental masses and their individual inertia tensors); a j is the total acceleration vector-field, including all segmental vector-fields, defined as the absolute (Bianchi) derivativev i of all segmental angular velocities v i =ẋ i , (i = 1, ..., n), where n is the total number of active DOF with local coordinates (x i ).
More formally, this central Law of robotics represents the covariant force functor F * constructed over robot's configuration manifold M rob = M and defined by the following commutative diagram: There is a unique smooth map from the right-hand branch to the left-hand branch of the diagram (2): It is called the Legendre transformation, or fiber derivative (for details see, e.g. [9]). The fundamental covariant force functor F * :

Lagrangian vs. Hamiltonian Approach to Humanoid Robotics
The humanoid's configuration manifold M rob = M is coordinated by local joint angular coordinates x i (t), i = 1, ..., n = total number of active DOF. The corresponding joint angular velocitiesẋ i (t) live in the velocity phase space T M (the tangent bundle of the configuration manifold M ), 7 which has the Riemannian geometry with the local metric form: is the material metric tensor defined by humanoid's mass-inertia matrix (composed of individual segmental masses m µ ) and dx i are differentials of the local joint coordinates x i on M . Besides giving the local distances between the points on the manifold M, the Riemannian metric form g defines the system's kinetic energy: 7 On the velocity phase-space manifold T M exists: T M is an orientable manifold, admitting the standard volume given by giving the Lagrangian equations of the conservative skeleton motion with kinetic-minus-potential energy Lagrangian L = T − V , with the corresponding geodesic form [9] d dt where subscripts denote partial derivatives, while Γ i jk are the Christoffel symbols of the affine Levi-Civita connection of the humanoid manifold M , given by The general form of autonomous Lagrangian humanoid robotics on the corresponding Riemannian tangent bundles T M rob and T M hum of the configuration manifolds M rob and M hum (precisely derived in [8]) can be formulated in a unified form as: where F i are all possible torque 1-forms, including robot's actuators, joint dissipations and external disturbances.
On the other hand, we develop the autonomous Hamiltonian robotics on humanoid's configuration manifold M rob = M in three steps, following the standard symplectic geometry prescription (see [2,9]): Step A Find a symplectic momentum phase-space (P, ω).
Recall that a symplectic structure on a smooth manifold M is a nondegenerate closed 8 2−form ω on M , i.e., for each x ∈ M , ω(x) is nondegenerate, and dω = 0.
Let T * x M be a cotangent space to M at m. The cotangent bundle T * M represents a union ∪ m∈M T * x M , together with the standard topology on T * M and a natural smooth manifold structure, the dimension of which is twice the dimension of M . A 1−form θ on M represents a section θ : M → T * M of the cotangent bundle T * M . P = T * M is our momentum phase-space. On P there is a nondegenerate symplectic 2−form ω is defined in local joint coordinates x i , p i ∈ U , U open in P , as ω = dx i ∧ dp i . In that case the From this condition one can see that the closed form (the kernel of the exterior derivative operator d) is conserved quantity. Therefore, closed p−forms possess certain invariant properties, physically corresponding to the conservation laws.
Also, a p−form β that is an exterior derivative of some is called exact (the image of the exterior derivative operator d). By Poincaré lemma, exact forms prove to be closed automatically, dβ = d(dα) = 0.
Since d 2 = 0, every exact form is closed. The converse is only partially true, by Poincaré lemma: every closed form is locally exact.
Technically, this means that given a closed p−form α ∈ Ω p (U ), defined on an open set U of a smooth manifold M any point m ∈ U has a neighborhood on which there exists a (p − 1)−form β ∈ Ω p−1 (U ) such that dβ = α| U . In particular, there is a Poincaré lemma for contractible manifolds: Any closed form on a smoothly contractible manifold is exact. coordinates x i , p i ∈ U are called canonical. In a usual procedure the canonical 1−form θ is first defined as θ = p i dx i , and then the canonical 2-form ω is defined as ω = −dθ.
A symplectic phase-space manifold is a pair (P, ω).
Step B Find a Hamiltonian vector-field X H on (P, ω). Let (P, ω) be a symplectic manifold. A vector-field X : P → T P is called Hamiltonian if there is a smooth function F : P → R such that X⌋ω = dF (X⌋ω ≡ i X ω denotes the interior product or contraction of the vector-field X and the 2-form ω). X is locally Hamiltonian if X⌋ω is closed.
Let the smooth real-valued Hamiltonian function H : P → R, representing the total humanoid energy H(x, p) = T (p) + V (x) (T and V denote kinetic and potential energy of the system, respectively), be given in local canonical coordinates x i , p i ∈ U , U open in P . The Hamiltonian vector-field X H , condition by X H ⌋ω = dH, is actually defined via symplectic matrix J, in a local chart U , as where I denotes the n × n identity matrix and ∇ is the gradient operator.
Step C Find a Hamiltonian phase-flow φ t of X H . Let (P, ω) be a symplectic phase-space manifold and X H = J∇H a Hamiltonian vector-field corresponding to a smooth real-valued Hamiltonian function H : P → R, on it. If a unique oneparameter group of diffeomorphisms φ t : An integral curve is said to be maximal if it is not a restriction of an integral curve defined on a larger interval of R. It follows from the standard theorem on the existence and uniqueness of the solution of a system of ODEs with smooth r.h.s, that if the manifold (P, ω) is Hausdorff, then for any point x = (x i , p i ) ∈ U , U open in P , there exists a maximal integral curve of X H = J∇H, passing for t = 0, through point x. In case X H is complete, i.e., X H is C p and (P, ω) is compact, the maximal integral curve of X H is the Hamiltonian phase-flow φ t : 10 of ω upon X H ). 9 Given a map f : X − → X ′ between the two manifolds, the pullback on X of a form α on X ′ by f is denoted by f * α. The pullback satisfies the relations for any two forms α, β ∈ Ω p (X). 10 The Lie derivative Luα of p−form α along a vector-field u is defined by Cartan's 'magic' formula (see [9]): Symplectic phase-flow φ t consists of canonical transformations on (P, ω), i.e., diffeomorphisms in canonical coordinates x i , p i ∈ U , U open on all (P, ω) which leave ω invariant. In this case the Liouville theorem is valid: φ t preserves the phase volume on (P, ω). Also, the system's total energy Recall that the Riemannian metrics g = <, > on the configuration manifold M is a positivedefinite quadratic form g : T M → R, in local coordinates x i ∈ U , U open in M . Given the metrics g ij , the system's Hamiltonian function represents a momentum p-dependent quadratic form where T * M is an orientable manifold, admitting the standard volume form For Hamiltonian vector-field, X H on M , there is a base integral curve γ 0 (t) = x i (t), p i (t) iff γ 0 (t) is a geodesic, given by the one-form force equatioṅ The l.h.sṗ i of the covariant momentum equation (8) represents the intrinsic or Bianchi covariant derivative of the momentum with respect to time t. Basic relationṗ i = 0 defines the parallel transport on T N , the simplest form of humanoid's dynamics. In that case Hamiltonian vector-field X H is called the geodesic spray and its phase-flow is called the geodesic flow.
For Earthly dynamics in the gravitational potential field V : M → R, the Hamiltonian H : T * M → R (7) extends into potential form with Hamiltonian vector-field X H = J∇H still defined by canonical equations (6).
A general form of a driven, non-conservative Hamiltonian equations reads: where F i = F i (t, x, p) represent any kind of joint-driving covariant torques, including active neuromuscular-like controls, as functions of time, angles and momenta, as well as passive dissipative and elastic joint torques. In the covariant momentum formulation (8), the non-conservative Hamiltonian equations (9) becomė The general form of autonomous Hamiltonian robotics is given by dissipative, driven Hamiltonian equations on T * M :ẋ including contravariant equation (10) -the velocity vector-field, and covariant equation (11) -the force 1-form (field), together with initial joint angles and momenta (12). Here R = R(x, p) denotes the Raileigh nonlinear (biquadratic) dissipation function, and F i = F i (t, x, p) are covariant driving torques of robot's actuators. The velocity vector-field (10) and the force 1−form (11) together define the generalized Hamiltonian vector-field X H ; the Hamiltonian energy function H = H(x, p) is its generating function. As a Lie group, the humanoid's configuration manifold M = j SO(3) j is Hausdorff. 11 Therefore, for x = (x i , p i ) ∈ U p , where U p is an open coordinate chart in T * M , there exists a unique one-parameter group of diffeomorphisms φ t : T * M → T * M , that is the autonomous Hamiltonian phase-flow: given by (10)(11)(12) The general form of Hamiltonian humanoid robotics on the symplectic cotangent bundle T * M rob of the configuration manifold M rob (as derived in [11,3,4]) is based on the affine Hamiltonian function H a : T * M → R, in local canonical coordinates on T * M given by where H 0 (x, p) is the physical Hamiltonian (kinetic + potential energy) dependent on joint coordinates x i and canonical momenta p i , H j = H j (x, p), (j = 1, . . . , m ≤ n are the coupling Hamiltonians corresponding to the system's active joints and u i = u i (t, x, p) are (reflex) feedback-controls.
Using (14) we come to the affine Hamiltonian control HBE-system, in deterministic forṁ p)), as well as in the fuzzy-stochastic form In (15)

Generalization to Human Biodynamics
If we neglect anatomy and physiology of human sensors and effectors, that is, from purely mechanical perspective, there are two main dynamical differences between robots and humans: (i) human joints are more flexible than robot joints (effectively many more degrees-of-freedom), and (ii) human dynamics is usually non-autonomous, or time-dependent. We will explain both differences in some detail in the following subsections.

Realistic Configuration Manifold of Human Motion
Every rotation in all synovial human joints is followed by the corresponding micro-translation, which occurs after the rotational amplitude is reached [3]. So, representation of human motion is rigorously defined in terms of Euclidean SE(3)-groups of full rigid-body motion [7,2,9] in all main human joints (see Figure 2). Therefore, the configuration manifold M hum for human dynamics is defined as a topological product of all included constrained SE(3) groups, M hum = i SE(3) i . Consequently, the natural stage for autonomous Lagrangian dynamics of human motion is the tangent bundle T M hum [8], and for the corresponding autonomous Hamiltonian dynamics is the cotangent bundle T * M hum [10,1,11].

Time-Dependent Biodynamics
Recall that in ordinary autonomous mechanics we have a configuration manifold M (which denotes both M rob and M hum ), coordinated by (x i ), and the corresponding velocity phase-space manifold is its tangent bundle T M , coordinated by (x i ,ẋ i ). However, in modern geometrical settings of non-autonomous mechanics, the configuration manifold of time-dependent mechanics is a fibre bundle π : M → R, called the configuration bundle, coordinated by (t, x i ), where t ∈ R is a Cartesian coordinate on the time axis R with the transition functions t ′ = t+const. The corresponding velocity phase-space is the 1-jet space J 1 (R, M ), which admits the adapted coordinates (t, x i , x i t ) = (t, x i ,ẋ i ). Every dynamical equation ξ defines a connection on the affine jet bundle J 1 (R, M ) → M , and vice versa [9].
Given the configuration fibre bundle M → R over the time axis R, we say that the 1−jet manifold J 1 (R, M ) is defined as the set of equivalence classes j 1 t s of sections s i : R → M of the bundle M → R, which are identified by their values s i (t), as well as by the values of their partial derivatives ∂ t s i = ∂ t s i (t) at time points t ∈ R. The 1-jet manifold J 1 (R, M ) is coordinated by (t, x i ,ẋ i ), so the 1-jets are local coordinate maps Similarly, the 2−jet manifold J 2 (R, M ) is the set of equivalence classes j 2 t s of sections s i : R → M of the configuration bundle π : M → R, which are identified by their values s i (t), as well as the values of their first and second partial derivatives, ∂ t s i = ∂ t s i (t) and ∂ tt s i = ∂ tt s i (t), respectively, at time points t ∈ R. The 2-jet manifold J 2 (R, M ) is coordinated by (t, x i ,ẋ i ,ẍ i ), so the 2-jets are local coordinate maps , any dynamical equation ξ on the configuration bundle M → R, which generalizes Lagrangian equation (4), is equivalent to the geodesic equation with respect to some affine connection Γ on the tangent bundle which fulfills the conditions A holonomic connection ξ is represented by the horizontal vector-field on J 1 (R, M ), A dynamical equation ξ is said to be conservative if there exists a trivialization M ∼ = R×M such that the vector-field ξ (19) on J 1 (R, M ) ∼ = R × T M is projectable onto T M . Then this projection Conversely, every second-order dynamical equation Ξ (20) on a manifold M can be seen as a conservative dynamical equation on the trivial fibre bundle R × M → R.

Nonautonomous Dissipative Hamiltonian Dynamics
We can now formulate the time-dependent biomechanics [13,14,15] in which the biomechanical phase space is the Legendre manifold 12 Π, endowed with the holonomic coordinates (t, y i , p i ) with the transition functions p ′ i = ∂y j ∂y ′ i p j . Π admits the canonical form Λ given by We say that a connection on the bundle Π → X is locally Hamiltonian if the exterior form γ⌋Λ is closed and Hamiltonian if the form γ⌋Λ is exact [16]. A connection γ is locally Hamiltonian iff it obeys the conditions: Note that every connection Γ = dt ⊗ (∂ t + Γ i ∂ i ) on the bundle Y − → X gives rise to the Hamiltonian connection Γ on Π − → X, given by

The corresponding Hamiltonian form H Γ is given by
Let H be a dissipative Hamiltonian form on Π, which reads: We call H and H in the decomposition (21) the Hamiltonian and the Hamiltonian function respectively. Let γ be a Hamiltonian connection on Π − → X associated with the Hamiltonian form (21). It satisfies the relations [16] γ⌋Λ = dp i ∧ dy i + γ i dy i ∧ dt − γ i dp i ∧ dt = dH, From equations (22) we see that, in the case of biomechanics, one and only one Hamiltonian connection is associated with a given Hamiltonian form. Every connection γ on Π − → X yields the system of first-order differential equations: They are called the evolution equations. If γ is a Hamiltonian connection associated with the Hamiltonian form H (21), the evolution equations (23) become the dissipative time-dependent Hamiltonian equations:ẏ In addition, given any scalar function f on Π, we have the dissipative Hamiltonian evolution equation relative to the Hamiltonian H. On solutions s of the Hamiltonian equations (24), the evolution equation (25) is equal to the total time derivative of the function f :

Neuro-Muscular Fitness Dynamics
The dissipative Hamiltonian system (24) Physiologically, the active muscular drives F i = F i (t, y, p, f ) consist of [2]): 1. Synovial joint mechanics, giving the first stabilizing effect to the conservative skeleton dynamics, is described by the (y,ẏ)-form of the Rayleigh-Van der Pol's dissipation function where α i and β i denote dissipation parameters. Its partial derivatives give rise to the viscousdamping torques and forces in the joints which are linear inẏ i and quadratic in y i .

Muscular mechanics,
giving the driving torques and forces F musc i = F musc i (t, y,ẏ) with (i = 1, . . . , n) for human biomechanics, describes the internal excitation and contraction dynamics of equivalent muscular actuators [2].
(a) The excitation dynamics can be described by an impulse force-time relation where F 0 i denote the maximal isometric muscular torques and forces, while τ i denote the associated time characteristics of particular muscular actuators. This relation represents a solution of the Wilkie's muscular active-state element equation [17] µ + Γ µ = Γ S A, µ(0) = 0, 0 < S < 1, where µ = µ(t) represents the active state of the muscle, Γ denotes the element gain, A corresponds to the maximum tension the element can develop, and S = S(r) is the 'desired' active state as a function of the motor unit stimulus rate r. This is the basis for biomechanical force controller.
(b) The contraction dynamics has classically been described by Hill's hyperbolic force-velocity relation [18] where a i and b i denote Hill's parameters, corresponding to the energy dissipated during the contraction and the phosphagenic energy conversion rate, respectively, while δ ij is the Kronecker's δ−tensor.
In this way, human biomechanics describes the excitation/contraction dynamics for the ith equivalent muscle-joint actuator, using the simple impulse-hyperbolic product relation

Spinal Control Level
The force HBE servo-controller is formulated as affine control Hamiltonian-systems (15-16) (with possible extensions along the lines of the previous section), which resemble an autogenetic motor servo [20], acting on the spinal-reflex level of the human locomotion control. A voluntary contraction force F of human skeletal muscle is reflexly excited (positive feedback +F −1 ) by the responses of its spindle receptors to stretch and is reflexly inhibited (negative feedback −F −1 ) by the responses of its Golgi tendon organs to contraction. Stretch and unloading reflexes are mediated by combined actions of several autogenetic neural pathways, forming the so-called 'motor servo.' The term 'autogenetic' means that the stimulus excites receptors located in the same muscle that is the target of the reflex response. The most important of these muscle receptors are the primary and secondary endings in the muscle-spindles, which are sensitive to length change -positive length feedback +F −1 , and the Golgi tendon organs, which are sensitive to contractile force -negative force feedback −F −1 .
The gain G of the length feedback +F −1 can be expressed as the positional stiffness (the ratio G ≈ S = dF/dx of the force-F change to the length-x change) of the muscle system. The greater the stiffness S, the less the muscle will be disturbed by a change in load. The autogenetic circuits +F −1 and −F −1 appear to function as servoregulatory loops that convey continuously graded amounts of excitation and inhibition to the large (alpha) skeletomotor neurons. Small (gamma) fusimotor neurons innervate the contractile poles of muscle spindles and function to modulate spindle-receptor discharge.

Cerebellum-Like Velocity and Jerk Control
Nonlinear velocity and jerk (time derivative of acceleration) servo-controllers [2], developed using the Lie-derivative formalism [1], resemble self-stabilizing and adaptive tracking action of the cerebellum [21]. By introducing the vector-fields f and g, given respectively by we obtain the affine controller in the standard nonlinear MIMO-system form (see [9]) Finally, using the Lie derivative formalism [9] 13 and applying the constant relative degree r to all HB joints, the control law for asymptotic tracking of the reference outputs o j R = o j R (t) could be formulated as (generalized from [22]) where c s−1 are the coefficients of the linear differential equation of order r for the error function e(t) = x j (t) − o j R (t) e (r) + c r−1 e (r−1) + · · · + c 1 e (1) + c 0 e = 0.
The control law (27) can be implemented symbolically in M athematica T M in the following three steps: 13 Let F (M ) denote the set of all smooth (i.e., C ∞ ) real valued functions f : M → R on a smooth manifold M , V (M )the set of all smooth vector-fields on M , and V * (M ) -the set of all differential one-forms on M . Also, let the vector-field ζ ∈ V (M ) be given with its local flow φ t : M → M such that at a point x ∈ M , d dt | t=0 φ t x = ζ(x), and φ * t representing the pull-back by φ t . The Lie derivative differential operator L ζ is defined: (i) on a function f ∈ F (M ) as (ii) on a vector-field η ∈ V (M ) as -the Lie bracket, and (iii) on a one-form α ∈ V * (M ) as In general, for any smooth tensor field T on M , the Lie derivative L ζ T geometrically represents a directional derivative of T along the flow φ t .
The affine nonlinear MIMO control system (26) with the Lie-derivative control law (27) resembles the self-stabilizing and synergistic output tracking action of the human cerebellum [24]. To make it adaptive (and thus more realistic), instead of the 'rigid' controller (27), we can use the adaptive Lie-derivative controller, as explained in the seminal paper on geometrical nonlinear control [23].

Cortical-Like Fuzzy-Topological Control
For the purpose of our cortical control, the dominant, rotational part of the human configuration manifold M N , could be first, reduced to an N -torus, and second, transformed to an N -cube ('hyperjoystick'), using the following topological techniques (see [9]). 15 Let S 1 denote the constrained unit circle in the complex plane, which is an Abelian Lie group. Firstly, we propose two reduction homeomorphisms, using the Cartesian product of the constrained SO(2)−groups: (2) and SO(2) ≈ S 1 .
Next, let I N be the unit cube [0, 1] N in R N and '∼' an equivalence relation on R N obtained by 'gluing' together the opposite sides of I N , preserving their orientation. Therefore, M N can be represented as the quotient space of R N by the space of the integral lattice points in R N , that is an oriented and constrained N -dimensional torus T N : Its Euler-Poincaré characteristic is (by the De Rham theorem) both for the configuration manifold T N and its momentum phase-space T * T N given by (see [9]) where b p are the Betti numbers defined as Conversely by 'ungluing' the configuration space we obtain the primary unit cube. Let '∼ * ' denote an equivalent decomposition or 'ungluing' relation. According to Tychonoff's producttopology theorem [9], for every such quotient space there exists a 'selector' such that their quotient models are homeomorphic, that is, T N / ∼ * ≈ A N / ∼ * . Therefore I N q represents a 'selector' for the configuration torus T N and can be used as an N -directional 'q-command-space' for the feedback control (FC). Any subset of degrees of freedom on the configuration torus T N representing the joints included in HB has its simple, rectangular image in the rectifiedq-command space -selector I N q , and any joint angle q i has its rectified imageq i .
In the case of an end-effector,q i reduces to the position vector in external-Cartesian coordinates z r (r = 1, . . . , 3). If orientation of the end-effector can be neglected, this gives a topological solution to the standard inverse kinematics problem.
Analogously, all momentap i have their images as rectified momentap i in thep-command space -selector I N p . Therefore, the total momentum phase-space manifold T * T N obtains its 'cortical image' as the (q, p)-command space, a trivial 2N -dimensional bundle I N q × I N p . Now, the simplest way to perform the feedback FC on the cortical (q, p)-command space I N q × I N p , and also to mimic the cortical-like behavior, is to use the 2N -dimensional fuzzy-logic controller, in much the same way as in the popular 'inverted pendulum' examples (see [19]).
We propose the fuzzy feedback-control map Ξ that maps all the rectified joint angles and momenta into the feedback-control one-forms so that their corresponding universes of discourse,Q i = (q i max −q i min ),P i = (p max i −p min i ) and U i = (u max i − u min i ), respectively, are mapped as Ξ : The 2N -dimensional map Ξ (29,30) represents a fuzzy inference system, defined by [24]: 1. Fuzzification of the crisp rectified and discretized angles, momenta and controls using Gaussianbell membership functions where χ ∈ D is the common symbol forq i ,p i and u i (q, p) and D is the common symbol for Q i ,P i and i ; the mean values m k of the nine partitions of each universe of discourse D are defined as m k = λ k D + χ min , with partition coefficients λ k uniformly spanning the range of D, with weights ̟ kl ∈ {0.6, 0.6, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 1.0}.

2.
Mamdani inference is used on each FAM-matrix µ kl for all human joints: (i) µ(q i ) and µ(p i ) are combined inside the fuzzy IF-THEN rules using AND (Intersection, or Minimum) operator, (ii) the output sets from different IF-THEN rules are then combined using OR (Union, or Maximum) operator, to get the final output, fuzzy-covariant torques, 3. Defuzzification of the fuzzy controls µ[u i (q, p)] with the 'center of gravity' method to update the crisp feedback-control one-forms u i = u i (t, q, p). Figure 4: The HBE simulating a jump-kick: calculating joint angles and muscular torques. Now, it is easy to make this top-level controller adaptive, simply by weighting both the above fuzzy-rules and membership functions, by the use of any standard competitive neural-network (see, e.g., [19]). Operationally, the construction of the cortical (q, p)-command space I N q × I N p and the 2N -dimensional feedback map Ξ (29,30), mimic the regulation of the motor conditioned reflexes by the motor cortex [21].
It has been implicitly assumed that close resemblance of hierarchical control structures presented in this section with the corresponding human neuro-physiological control mechanisms would assure the necessary overall stability of biodynamics. However, in future work, these control structures need to be properly analyzed, starting with Lyapunov stability criteria.

Simulation Examples
In this section we give several simulation examples of the sophisticated virtual humanoid called Human Biodynamics Engine (HBE), designed at Defence Science & Technology Organisation, Australia. The HBE is a sophisticated human neuro-musculo-skeletal dynamics simulator, based on generalized Lagrangian and Hamiltonian mechanics and Lie-derivative nonlinear control. It includes 270 active degrees of freedom (DOF), while fingers are not separately modelled: 135 rotational DOF are considered active, and 135 translational DOF are considered passive. The HBE incorporates both forward and inverse dynamics, as well as two neural-like control levels. Active rotational joint dynamics is driven by 270 nonlinear muscular actuators, each with its own excitation-contraction dynamics (following traditional Hill-Hatze biomechanical models). Passive translational joint dynamics models visco-elastic properties of inter-vertebral discs, joint tendons and muscular ligaments as a nonlinear spring-damper system. The lower neural control level resembles spinal-reflex positive and negative force feedbacks, resembling stretch and Golgi reflexes, respectively. The higher neural control level mimics cerebellum postural stabilization and velocity target-tracking control. The HBE's core is the full spine simulator, considering human spine as a chain of 26 flexiblycoupled rigid bodies (formally, the product of 26 SE(3)-groups). The HBE includes over 3000 body parameters, all derived from individual user data, using standard biomechanical tables. The HBE incorporates a new theory of soft neuro-musculo-skeletal injuries, based on the concept of the local rotational and translational jolts, which are the time rates of change of the total forces and torques localized in each joint at a particular time instant. Figure 6: The HBE simulating an effect of an aircraft pilot-seat ejection to human spine compression: before the seat ejection (left) and after ejection (right).
The first version of the HBE simulator had the full human-like skeleton, driven by the generalized Hamiltonian dynamics (including muscular force-velocity and force-time curves) and two levels of reflex-like motor control (simulated using the Lie derivative formalism) [1,2]. It had 135 purely rotational DOF, strictly following Figure 1. It was created for prediction and prevention of musculo-skeletal injuries occurring in the joints, mostly spinal (intervertebral). Its performance looked kinematically realistic, while it was not possible to validate the driving torques. It included a small library of target movements which were followed by the HBE's Lie-derivative controllers with efficiency of about 90% (see Figures 3 and 4).
The HBE also includes a generic crash simulator, based on the simplified road-vehicle impact simulator (see Figure 5). While implementing the generic crash simulator, it became clear that purely rotational joint dynamics would not be sufficient for the realistic prediction of musculo-skeletal injuries. In particular, to simulate the action of a Russian aircraft ejection-seat currently used by the American space shuttle, we needed to implement micro translations in the intervertebral joints (see Figure 7: The HBE calculating translational forces distributed along the spinal joints during the seat ejection. Figures 6 and 7). This is because the seat provides full body restraint and hence the ejection rockets firing with 15 g per .15 s cause pure compression of the spine (without any bending).
Finally, the HBE includes the defence-specific land-mine crash simulator. It is calibrated on a hypothetical double-impact under the armor-protected military vehicle, including: 1. A land-mine blast of 350g with a duration of 5ms; 2. A 1s pause when the hypothetical vehicle is in the air; and 3. The vehicle hard landing with an acceleration of 100g and a duration 1s.
The HBE calculates full rotational and translational dynamics caused by the land-mine doubleimpact in extreme force/time scales (including both linear and angular displacements, velocities, forces and jolts in all human joints (see Figure 8)). The variations of the applied g-forces and durations of the two impacts can be simulated, to see the differences in their effects on the hypothetical passenger's body.
In this way a full rotational + translational biodynamics simulator has been created with 270 DOF in total (not representing separate fingers). The 'HBE-simulator' has been kinematically validated (see Figure 9) against the standard biomechanical gait-analysis system 'Vicon' [25].