Recursive inverse dynamics sensitivity analysis of open-tree-type multibody systems

We present a first-order recursive approach to sensitivity analysis based on the application of the direct differentiation method to the inverse Lagrangian dynamics of rigid multibody systems. Our method is simple and efficient and is characterized by the following features. Firstly, it describes the kinematics of multibody systems using branch connectivity graphs and joint-branch connectivity matrices. For most mechanical systems with an open-tree kinematic structure, this method turns out to be more efficient compared to other kinematic descriptions employing joint or link connectivity graphs. Secondly, a recursive sensitivity analysis is presented for a dynamic system with an open-tree kinematic structure and inverse dynamic equations described in terms of the Lagrangian formalism. Thirdly, known approaches to recursive inverse dynamic and sensitivity analyses are modified to include dynamic systems with external forces and torques acting simultaneously at all joints. Finally, the proposed method for sensitivity analysis is easy to implement and computationally efficient. It can be utilized to evaluate the derivatives of the dynamic equations of multibody systems in gradient-based optimization algorithms. It also allows less experienced users to perform sensitivity analyses using the power of high-level programming languages such as MATLAB. To illustrate the method, simulation results for a human body model are discussed. The shortcomings of the method and possible directions for future work are outlined.


Introduction
Recent advancements in the analysis of nonlinear dynamics of mechanical systems has been significantly facilitated by the progress in numerical methods, optimal control, nonlinear control, optimal design, system identification, and related fields. Combination of new efficient techniques with a rapidly growing computational power enables researchers to tackle problems deemed unattainable a decade ago. The problem of optimal control of a nonlinear dynamic sys-tem such as human body is one such example. Our main motivation in this paper is the implementation of the dynamic and sensitivity analysis of the human body.
Equations of motion for multibody systems can be rough-ly divided into two large classes in terms of their mathematical appearance or form: closed-form and recursive equations [24] . 1 In the closed form, also called a bulk form, dynamic equations of motion are written as a single vector differential equation containing all generalized coordinates. In the recursive form, all terms in the dynamic equations for generalized coordinates of one joint are presented as functions of generalized coordinates of neighboring joint. The main advantage of the closed form is that is it easier to handle and analyze with analytical tools, primarily for the design of controllers in the state space. However, in this case equations of motion become complicated and highly nonlinear, especially as the degrees of freedom (DOF) of the system increase. The attractive advantages of the recursive formulation are its computational efficiency and the ease of coding. The main disadvantage is the "distortion" of the structure of the dynamic model to the form that excludes the possibility of gaining the insight into the dynamics of the system and exploiting it efficiently [24]. However, some recursive formulations do not have this shortcoming and still allow to gain insight into dynamics [33,48].
Multibody dynamics can be also classified as forward dynamics (kinematic variables are not known) and inverse dynamics (force/torque variables are not known) [12,26,51]. However, in applied problems neither positions and velocities nor input forces and torques are known in advance; these variables should be determined during simulation. The problems where both kinematic and dynamic variables are not known fall into the third category; they are formulated as optimization problems and solved using optimization tools. In various engineering fields this class of problems is called trajectory optimization, motion synthesis, motion control, or predictive dynamics [1]. Efficient solution of optimization problems with gradientbased methods depends largely on the gradients of the objective function and constraints. Calculation of gra-dient matrices and vectors is called sensitivity analysis. In essence, sensitivity analysis is the inverse of the optimization problem. In the optimization problem, the task is to minimize the objective function. This is equivalent to the following problem: given the value for the gradient of the objective function, find the value of the independent variable which satisfies it. In the sensitivity analysis the problem is the opposite: given the value for the independent variable, evaluate the gradient of the objective function. Sensitivity analysis can be performed with respect to states (generalized coordinates, velocities, and accelerations) [1,57]. Sensitivity analysis with respect to physical system parameters (mass, inertia terms, lengths, etc.) or any other design parameters (such as initial conditions) is termed design sensitivity [27, 53,55]. Sensitivity analysis is important for problems involving design optimization, parameter estimation, and model correlation. It is also crucial for implicit numerical integration of differential-algebraic systems and kinematic workspace analysis of multibody systems [49].
Sensitivity analysis can be performed using analytical, semi-analytical, and numerical methods [10]. Methods for numerical sensitivity analysis can be further classified as finite difference and automatic differentiation methods. Analytical methods include direct differentiation and adjoint variable method [5]. We draw the reader's attention to the fact that the above classification is not standardized and accurate, because there are hybrid methods that combine the aforementioned techniques. For example, there might be direct and/or adjoint variable method combined with numerical and/or automatic differentiation [9,10,37]. The above classification is provided only for information purposes. The finite difference method is commonly employed for evaluating sensitivity matrices. It is simple and easy to implement, but its accuracy depends on the magnitude of the perturbation and is affected by truncation and rounding errors. Therefore, finite difference method may lead to ill-conditioned problems and unreliable results; its computational costs are prohibitively high for large DOF systems [19]. The automatic differentiation method is based on the application of the symbolic chain rule of differentiation; its accuracy and speed are superior to the finite difference method. Open-source software based on automatic differentiation method and written in C++ (DRAKE [52], MBSlib [56], and RobCoGen [45]) can be used to perform sensitivity analysis of dynamic equations with respect to generalized coordinates, velocities, and parameters. The advantages of the automatic differentiation method include possibilities to simulate forward, inverse, and hybrid dynamics, parameter estimations, and trajectory optimization. The main limitation is the need to couple automatic differentiation software with that employed by the user; this might require additional learning. For instance, it is known that the wrong application of the chain rule of differentiation can lead to erroneous results [16]. Furthermore, computation of implicit derivatives often requires costly time integration. The direct differentiation method is based on a straightforward application of differentiation rules in the equations of motion [17,49]. Its advantages include easiness of implementation, accuracy, and higher numerical stability [10]; the method is especially efficient when the number of sensitivity variables is small. Examples of the application of the direct differentiation method for sensitivity analysis of multibody mechanical systems can be found in [6,15,47]. In the adjoint variable method, the objective function is modified with dynamic equations and constraints and additional adjoint conditions that simplify the computation of the Jacobian of the objective function are obtained [28,44]. This method is accurate and numerically efficient, but can become complex and lengthy since it requires backward integration in time [7,46]; it is efficient when the number of objective functions is small. Relevant examples of the application of the adjoint variable method for sensitivity analysis of multibody mechanical systems can be found in [4,7,16,29,44]. One can also find in the literature approaches based on the combination of existing methods. For instance, a semianalytical sensitivity analysis replaces symbolic derivatives with respect to design variables in analytic expressions by finite differences combining the simplicity of the finite difference method with the accuracy of the adjoint variable and direct differentiation methods in [46].
Most research on the sensitivity analysis uses a closed-form dynamic formulation where all equations are written using a single vector-matrix form [6,7,10,[15][16][17][18]29,44,46,47,49]. To the best of our knowledge, only a few papers on the recursive sensitivity analysis of mechanical systems with open-tree topology structure are available, although it is often computationally more efficient compared to "bulk" formulations for large multibody systems [5]. 2 Closed-form expressions for the sensitivity of multibody system dynamics equations using the spatial notation provided in [25] are abstract and their correct implementation requires substantial user effort. A useful framework for generating motion sequences in musculoskeletal systems was designed by combining the computational tools of MATLAB with the modeling capabilities of OpenSim [39]. However, the sensitivity analysis is not presented there; it was performed numerically using forward finite differences. As expected, the authors noted that the fmincon solver was very slow. Sensitivity analysis for simple open-chain multibody structures using inverse recursive Newton-Euler dynamic formulation and spatial notation (as in [21]) was performed in [20]. Using the aggregate force and momentum expressions (changing the order in which terms are computed), linear time analytical first derivatives of the objective function were obtained. However, this method provides true values not for all joint momenta and torques but only for the base joint. We would like to emphasize that the direct differentiation sensitivity analysis using the inverse recursive Lagrangian formulation [1,57], forward recursive Kane's notation [3], and transfer matrix method [54] are valid for simple open-chain structures.
We use Lagrange's formulation of the system dynamics for improving the sensitivity analysis method suggested in [1,57] for open-tree topological structures. Our work complements the recursive sensitivity analysis techniques based on the Kane's method [8,31,43] and on the Newton-Euler method with the spatial notation [33,50]. Contrary to the sensitivity analysis derived in [8,31,43] for forward recursive dynamics, we perform it for inverse recursive dynamics. In comparison with the forward dynamics formulation, the inverse dynamics does not require time integration and increases the efficiency of the sensitivity analysis. The main differences between our proposed method to describe system topology and the method used in [33] are the following. Firstly, we propose a simple method to describe system topology using the branch connectivity graph and the joint-branch connectivity matrix. In contrast, the square adjacency matrix and the block-weighted adjacency matrix are used in [33]. The adjacency matrix carries the same information as the branch connectivity graph and the joint-branch connectivity matrix but requires more space and so is less efficient. The block-weighted adjacency matrix has a larger size than the adjacency matrix. Secondly, in this work, a homogeneous transformation matrix is used rather than the spatial operator algebra notation (6-D vector formulation) as in [33,50]. Our work extends the possibilities of recursive dynamics and sensitivity analysis presented in [1,57] by allowing the inclusion of external forces and torques acting simultaneously at all joints and/or links. This is necessary, for example, for the modeling of two reaction forces acting on both feet during the double support phase of walking. We expect that the proposed sensitivity analysis will reduce the time to evaluate sensitivity information compared to traditional numerical methods. This in turn might reduce the time required by the conventional fmincon solver for computations required in optimization problems for nonlinear dynamic systems with many DOF, such as the human body. Since using MATLAB for testing and debugging is easier than, for instance, employing the interior point optimizer IPOPT [39], our contribution should positively impact research on multibody dynamics.
In summary, we present a first-order recursive sensitivity analysis based on the application of the direct differentiation method to the inverse Lagrangian dynamics of rigid multibody systems which is not affected by the perturbations of design variables. Major contributions in this paper are: -a new method for describing the topology of mechanical systems with an open-tree structure is proposed; -known recursive dynamic and sensitivity analyses are modified for the use with dynamic systems having an open-tree structure where external forces and torques act simultaneously on all joints; -the proposed algorithm can be easily implemented in MATLAB thus allowing the use of high-level programming capabilities for the human motion synthesis.
The paper is organized as follows. Section 2 introduces kinematics of the multibody mechanical system. Dynamic equations of motion are defined in Sect. 3, followed by the sensitivity analysis in Sect. 4. Simulation model is presented in Sect. 5 and simulation results are discussed in Sect. 6.

Joint connectivity graph
A mechanical system can be represented as a finite number of connected links and joints. A graph, on the other hand, is a collection of nodes and edges [32]. To correctly represent the kinematics of a mechanical system, the connectivity of the joints should be described. To this end, we utilize the joint connectivity graph, an undirected graph similar to the one used in [22,23] where a node represents a link and an arc represents a joint. Since every type of joint (prismatic, revolute, spherical, universal joint, etc.) has to be defined separately, this approach is less efficient in practice. Similarly, each joint type requires its own kinematic description in [33]. Our modification of the connectivity graph employs nodes (vertices) for representing joints, while the arcs (edges) represent connections between neighboring joints. If a mechanical system has an open-tree structure (that is, it does not have closed-loop chains), then its joint connectivity graph is a joint topological tree. On a joint connectivity graph, all joints are numbered so that a given joint i ∈ N has a lower number than any of its children, and a higher number than any of its parents, see Fig. 1a. In other words, if a number λ j (i) denotes a parent and ν j (i) denotes a child of a joint i, then This indexing is opposite to (or inverse of) the canonical tree notation used in [33,34]. It follows from (1) that the root node (joint) has the number 1, while in [22,23] the root node (link) starts with number 0. For a joint topological tree, a joint has only a single parent but can have multiple children. Let the number of children of the node i be denoted by κ j (i). We call a node that has at most one child a simple node and the one that has more than one child a complex node (also called junction node in [32] ). A simple node that does not have children nodes will be called an external node. For example, nodes 3 and 5 in Fig. 1a are complex, nodes 6, 7, 9 and 11 are external, and all remaining nodes are simple. For a root node (joint) λ j (i) = 0, and for an external node ν j (i) = 0. A joint numbering system is not unique, and a given joint topological tree can have multiple correct ways of doing this. Ideally, the knowledge of λ j (i) for each i furnishes complete information about the topological structure of the sys- Fig. 1 Topological tree of a mechanical system with 11 joints and 6 branches: a joint connectivity graph with node numbers representing joints, b branches of the joint connectivity graph, c branch connectivity graph with node numbers representing branches tem which can be expressed, for example, in the form of the so-called joint parent matrix Υ ∈ R n j ×2 , where n j is the total number of joints. The first column Υ (i, 1) would specify the current joint number i, and the second column would specify the number of its parent Compared to other descriptions of the system topology in terms of adjacency and incidence matrices 3 [32,34,35,42], the joint parent matrix is already compact. However, the description of the information in the joint parent matrix form can be further compressed and made more efficient.

Branch connectivity graph and joint-branch connectivity matrix
In this section, we suggest a more efficient (compact) way to describe the topological structure of a mechanical system. Analyzing joint topological trees, we conclude that most nodes (joints) are simple and their description within the joint connectivity graph is rather straightforward, whereas difficulties usually arise in the description of complex nodes. For this purpose, a branch of the joint connectivity graph is defined as a set of its nodes that have simple linear connections, Fig. 1b. One might notice that any branch starts either with a root node or with a child of a complex node and ends either with an external node or with a complex node. Definitions of nodes and branches are analogous to those used in electric circuit theory. A branch con-nectivity graph is defined similarly to the joint connectivity graph, but the difference is that each node represents a single branch, Fig. 1c, that is, a node in the former graph represents a single graph branch of the latter one. If a mechanical system has an opentree structure, its branch connectivity graph is a branch topological tree. Similarly to the process of joint numbering, nodes of a branch connectivity graph (open-tree branches of a mechanical system) are numbered so that a given tree branch i ∈ N has a lower number than any of its child branches and a higher number than any of its parent branches. For a branch k with a parent λ b (k) and a child ν b (k), the identity similar to (1) holds, For a root node (branch) one has λ b (k) = 0 and for an external node ν b (k) = 0. Simply speaking, a branch connectivity graph is a "lighter" version of a joint connectivity graph where all simple nodes are removed. Therefore, for any node in a branch connectivity graph, one has ν b (k) = 1, which also implies that the branch connectivity graph cannot be further simplified by a similar procedure. In other words, a branch connectivity graph is the simplest representation of the joint connectivity graph which preserves all information about complex nodes. Unless stated otherwise, in what follows the index i refers to the joint number and the index k denotes the branch number. Let κ b (k) denote the number of children of a node in a branch connectivity tree. From the definition of a branch, where i is the last node of the branch k in the joint connectivity graph, that is, the number of children of the branch k is equal to the number of children of its last joint. In general, i = k, because the number of joints is not equal to the number of branches. Since our focus is on a branch connectivity graph and not on simple joint nodes, for the sake of simplicity, we slightly abuse the notation denoting by κ(i) the number of children of a node. Complete information about the branch topological tree of a given system is provided if for each branch either λ b (k) or ν b (k) is specified. This can be achieved by specifying a branch parent matrix Φ ∈ R n b ×2 , where n b is the total number of branches in the branch topological tree for the given system. The first column of the matrix shows that the current branch number Φ(k, 1) = k, and the second column specifies its parent branch number Φ(k, 2) = λ b (k). The branch numbering is not unique and can be performed in multiple ways.
When the branch connectivity graph is specified, one cannot construct the kinematic tree of the dynamic system yet because the knowledge about the joints that compose each branch is missing. This information can be provided in a joint-branch connectivity matrix Ψ ∈ R n b ×2 . Under the assumption that the joints in a given branch are numbered consecutively, the first column in Ψ (k, 1) specifies the joint number where the branch k starts, whereas the second column Ψ (k, 2) specifies the joint number where it terminates. There are two alternative ways to describe system's kinematics structure, using either the joint connectivity graph (this requires a single matrix Υ ) or the branch connectivity graph in conjunction with the joint-branch connectivity matrix (this requires two matrices Φ and Ψ ). If a mechanical system possesses many joints and a few branches, using the latter approach is simpler and more efficient whereas for a system with many branches and joints the former description is preferable. For example, for a system with one hundred DOF and only four branches, the branch topological tree has four nodes, and the joint-branch connectivity matrix is a 4 × 2 matrix. Recall that a node represents a joint in a joint connectivity graph and a branch in a branch connectivity graph. Therefore, for a system whose kinematic structure is represented by a full binary tree, the joint connectivity graph requires less information. Our method to describe the system's topology is also different (simpler) than in [33]. However, this simplicity is only related to the description of the system's topology, and it does not affect the sensitivity analysis. In other words, the sensitivity analysis does not depend much on the method to describe the topology of the system.

Denavit-Hartenberg convention
The Denavit-Hartenberg (D-H) convention is used to formulate the kinematics of the model. It specifies four parameters for each DOF, θ i , d i , a i , and α i , where the variables θ i and d i are associated with the cases where the joint is revolute or prismatic, respectively. We note that there are multiple variations of the D-H convention. In this paper, the distal version of the D-H convention is utilized. This means that the i-th local frame is attached to the distal end of the link i [24], and a joint i is located at the proximal end of the same link. Alternatively, the same local frame can be attached to the proximal end of the link [12]. Even though all D-H conventions describe the kinematics of arbitrary mechanical systems equally well, care should be taken when working with the D-H convention because of the aforementioned notational differences.

Dynamics
Assume that we have a mechanical n q -DOF system, then the number of joints n q = n j . For a fully actuated system, the number of input (driving) torques is equal to the number of DOF. If the number of input torques is less than the number of joints, then one has an underactuated system; otherwise the system is over-actuated. A human body is an over-actuated system, which is a beneficial feature in terms of redundancy, because a single muscle can actuate multiple joints and a single joint can be actuated by multiple muscles. However, this redundancy is difficult to model, especially when our goal is to generate a dynamically feasible human body motion. To achieve this, it suffices to consider a fully actuated human body model which determines our main focus on fully actuated multibody systems. The detailed derivation of the dynamic equations of motion is outside the scope of this paper. Our work is based on the Lagrangian formulation, so a reader interested in more details about the derivation (such as Lagrangian function), can consult [1]. Dynamic equations of motion for a mechanical n q -DOF system are often represented in the form where q = [q 1 q 2 . . . q n q ] T ∈ R n q andq = [q 1q2 . . .q n q ] T ∈ R n q are the vectors of generalized coordinates and generalized velocities, respectively, M(q) ∈ R n q ×n q is the inertia matrix, K (q,q) ∈ R n q is the vector of Coriolis and normal inertial forces, W (q) ∈ R n q is the vector of gravitational forces, and Q = [Q 1 Q 2 . . . Q n q ] T ∈ R n q is the vector of generalized forces. The superscript (·) T denotes the transpose operator. These equations can be referred to as "bulk" form, because all DOF are included in the equation, and are therefore computed simultaneously [24]. However, the direct computation of this set of equations is burdensome, especially for high-DOF systems with many joints and links. Instead, there are efficient "recursive" techniques of calculation where each DOF is computed recursively from the previous DOF [11,30,41]. We use the recursive Lagrangian formulation to derive the dynamic equations of motion based on the extension of ideas reported in [1,57,58] to the case of multiple external loads and topological tree structures.

Recursive kinematics
For each DOF, the corresponding homogeneous transformation matrix A i ∈ R 4×4 is defined as where i = 1, . . . , n q [1]. The matrix transforms homogeneous coordinates from the frame i to its parent frame i − 1. The first, second and third-order derivatives of A i with respect to q i are defined accordingly [24] where the matrix Z ∈ R 4×4 has all zero entries except for the following: Z (1, 2) = −1 and Z (2, 1) = 1 for a revolute joint, and Z (3, 4) = 1 for a prismatic joint. The global homogeneous transformation matrix is introduced using the local transformation matrices as . The first and second-order time derivatives of the global transformation matrix are . These matrices can be calculated in recursive form starting from the first, root joint, until the last, end-effector joint, in the kinematic chain: whereq i = dq i /dt. For a system which has only simple nodes in the joint topological tree or for a single branch within the joint topological tree, the subscript p = i −1. For a node that has a complex node parent, if the joint connectivity graph is used, the subscript p = λ j (i) = Υ (i, 2) and if the branch connectivity graph is used, p = Ψ (λ b (k), 2) = Ψ (Φ(k, 2), 2), where k is the index of the branch to which the joint i belongs. For the first joint (i = 1), one sets p = 0. In summary, The initial values of the global transformation matrix T i and its time derivatives are T 0 = I , S 0 = 0, and R 0 = 0, where I ∈ R 4×4 and 0 ∈ R 4×4 are, correspondingly, the identity and zero matrices. With the help of homogeneous matrices, the position r g , velocity v g , and acceleration a g of any point of the mechanical system with respect to the global frame can be specified as where r i is the position of the point in the local coordinate system i. The set of equations in (5) defines what is called the recursive forward kinematics.

Recursive dynamics
After the forward kinematics matrices are computed, dynamic matrices can be obtained. Firstly, we define the mass m i ∈ R and inertia matrix J i ∈ R 4×4 of the link associated with a joint i and expressed in the local frame i. The inertia matrix is defined as where x, y, z, are the moments and products of inertia of link i in its local frame of reference, while c x i , c y i , and c z i denote the location of the center of mass of link i in the coordinate system i. Secondly, the inertia H i ∈ R 4×4 and the external force F i ∈ R 4×n q matrices, the gravity E i ∈ R 4 and the external torque G i ∈ R 4 vectors are defined. For simplicity, we refer to these matrices as dynamic transformation matrices, because they transform dynamic quantities from a given joint to the next one. These matrices are also defined recursively, but this time from the last joint (end-effector) to the first joint (root) If a system does not have complex nodes in its joint topological tree or within a branch κ(i) = 1, the subscript p = i + 1. For a complex node in the joint topological tree p = ν j (i), and in the branch topological tree p = Ψ (ν b (k), 1) = Ψ (Φ( j (k), 1), 1), where, as above, k is the index of the branch to which the joint i belongs. The index notation j (k) is used for the children of the branch k found through the branch parent matrix by using the identity k = Φ( j, 2). For the last joint in the branch without child joints, κ(i) = 0; we can also assume that the initial values of the dynamic matrices and vectors are D n q +1 = 0 and E n q +1 = F n q +1 = G n q +1 = 0. The same initial values are used for any external node in the joint connectivity graph. In summary, The generalized force Q i acting on the link i (i-th component of the vector Q from (3)) can be found by a similar recursive process starting from the last link to the first link, for i = 1, . . . , n q , tr[·] is the trace operator, and tr[·] i:n q denotes the sum of diagonal matrix elements from i until n q , g = [g x g y g z 0] T is the gravity vector expressed in the global frame of reference, 5 the vector w 0 is defined as w 0 = [0 0 1 0] T for a revolute joint and as w 0 = 0 ∈ R 4 for a prismatic joint. The subscript p in (9) is defined as in (5): The external force f i and the torque h i are both expressed in the global coordinate system. The process of calculating the torque and 4 If there are multiple contact and/or distributed forces and torques acting on the link, then we can always find their net resultants. 5 in our case g = [0 0 − 9.81 0] T . the terms defined in (8) is called the recursive backward dynamics.

Sensitivity analysis
In gradient-based optimization algorithms, convergence to the local minimum is highly dependent on the gradient information. It is known that the computation of the derivatives of the objective function and constraints in some cases consumes over 90% of the CPU time required for each iteration of the algorithm [2]. The computation time depends on the selection of the set of generalized coordinates, constraints, constraints enforcement schemes, the set of sensitivity parameters, differentiation methods, and other factors. Thus, for solving an NLP problem efficiently, it is desirable to provide gradients of the object function and constrains with respect to optimization variables. Furthermore, the gradient and Jacobian information enable a coherent and easy extension of algorithms to more complex models making algorithms well scalable. To this end, a sensitivity analysis is performed, i.e., partial derivatives of the desired expressions with respect to the generalized coordinates, velocities, and accelerations are determined. There are two reasons for the utilization of the direct differentiation method for sensitivity analysis in this paper. The first is that the adjoint variable method requires integration of the adjoint sensitivity equations. However, we use the inverse dynamics formulation to avoid the integration of dynamic equations of motion. Thus, the adjoint variable method would neutralize our effort to avoid the time integration. Secondly, as already mentioned, our work is extension of the work in [1], and so we continue to use the direct differentiation methods as they did.

Kinematic sensitivity analysis
The sensitivity of homogeneous transformation matrices with respect to the generalized coordinates, velocities and accelerations is written as where i = 1, . . . , n q , m = 1, . . . , n q . Since T i depends only on the generalized coordinates, its partial derivatives with respect to generalized velocities and accelerations are zeros. Similarly, partial derivatives of S i with respect to generalized accelerations are zeros. The homogeneous transformation matrices can be found from the following relations. For m < i, For the case (m = i): (10) and (11) if, in addition to the conditions imposed at the beginning of each equation, i = Ψ (k, 1) ∧ λ b (k) = 0. These additional conditions are used to determine whether a parent branch exists. The index k refers to the branch number associated with the joint number i, e.g., it is understood that the joint i is located within the branch k. If i = Ψ (k, 1)∧λ b (k) = 0, then p = 0 in (11) for m = i, but for m < i (12) applies. If i = Ψ (k, 1), then p = i − 1. Alternatively, in the joint connectivity graph notation, p = λ j (i) = Υ (i, 2).

Dynamic sensitivity analysis
Similarly to the kinematics sensitivity analysis, sensitivity of the dynamic transformation matrices with respect to generalized coordinates, velocities and accelerations is written as where i = 1, . . . , n q , m = 1, . . . , n q . These matrices can be found from the following relations. For m ≤ i, For m ≥ i + 1, For the case m ≥ i where p = Ψ (ν b (k), 1) = Ψ (Φ( j (k), 1), 1) provided that, in addition to the conditions imposed at the beginning of each equation, i = Ψ (k, 2) ∧ ν b (k) = 0. These conditions verify the existence of children of the branch k and the sum includes all its child branches κ(i). (13) for m ≤ i, but for m ≥ i + 1 (15) applies. In the case i = Ψ (k, 2), one has p = i + 1 and κ(i) = 1. The coefficient β = 1 if m = p, which means that m = Ψ (ν b (k), 1) for the case i = Ψ (k, 2) and m = i + 1 for the case i = Ψ (k, 2), otherwise, β = 0. In the joint connectivity graph notation, the index p = ν j (i).
Finally, the torque sensitivity can be found as follows.
where, as in (9) The corresponding pseudo-algorithms are summarized in Algorithms 1 -3, which require three loops, the kinematic loop from the root to the external nodes, the dynamic loop from external nodes to the root, and, finally, the torque loop back again from the root to external nodes. It is possible to include computations in Algorithm 3 in Algorithm 2 so that only two loops are required. For the sake of clarity, however, we present the sensitivity analysis in three loops.

Human model
A rigid human body model consisting of rigid links and revolute joints is used to test the proposed sensitivity method. A complete human body has about 350 joints [40]. For our purposes, the human model has 43 DOF and 20 links. The number of joints exceeds the number of links due to the presence of virtual links which have zero length and mass. For example, a spherical joint is represented as a combination of three consecutive rotational joints. The 43-DOF human model is shown schematically in Fig. 2 where the green dot indicates the origin of the global coordinate system. The red dot indicates the pelvis point used as the reference point for the human body model. The global coordinate system has axes x 0 -y 0 -z 0 , and the local coordinate system i, associated with the joint i and attached to the distal end of the link i, has the unit vectors x i -y i -z i . The mobility of the human model with respect to the global frame of reference is achieved by the global virtual DOF of the joints connected to the pelvis point. Global virtual DOF consists of three translational and three rotational DOF. All other human links and joints are referenced through the pelvis point. The joint and branch connectivity graphs for the human model are shown in Fig. 3. There are 43 joints and 7 branches. The corresponding joint parent matrix Υ ∈ R 43×2 , the branch parent matrix Φ ∈ R 7×2 , and the joint-branch connectivity matrix Ψ ∈ R 7×2 are defined by (18)- (20). Even for this modestly-sized system, the branch connectivity graph allows a more compact and efficient description of the system topology in comparison with the joint connectivity graph. The values of D-H parameters for the human model are collected in Table 1

Implementation
We performed the sensitivity analysis using the MAT-LAB software. Due to the inverse dynamics formulation utilized in this work, the reference trajectories      a simple time differentiation of the generalized coordinates. These reference trajectories were discretized with a sampling time of 10 ms, whereas the duration of the simulation was set to 3 s, resulting in a total of 301 sampling instances. At each sampling instance, both numerical and suggested analytical sensitivity values were evaluated. Numerical sensitivity was estimated from (9) using the jacobianest function of MAT-LAB, which is a fully adaptive and robust numerical differentiation tool [14]. The analytical sensitivity was obtained from (16)- (17 ). Note that for a system with 43 DOF, the total number of different torque sensitivity values ∂ Q i ∂q m , i, m = 1, . . . , 43, is equal to 43 2 = 1849. Similarly for the torque sensitivity values with respect to velocities and accelerations, but some values would Note that the index j in this table corresponds to the indexes in columns four-five in Table 1 be zero. For example, for torque and joint angles from different branches, such as i = 42 and m = 18 in Fig. 3, the torque sensitivity is zero because these branches do not influence each other. The root mean squared (RMS) error was also calculated. The simulations were performed on a ThinkPad notebook with an Intel Core i5-10310U CPU processor and 16 GB RAM. The sensitivity analysis results are provided in the next section.

Results
The results for randomly selected torque sensitivity values ∂ Q 20 /∂q 5 , ∂ Q 30 /∂q 25 , and ∂ Q 42 /∂q 22 are shown in Fig. 4 with the obtained sensitivity trajectories in the left column and the difference between the numerical and analytical methods in the right column. In other words, the error plots on the right show a simple difference between analytical and numerical sensitivities. The results obtained by using the analytical method agree well with the numerical results. The example where the numerical and analytical results differ is shown in Fig. 5. Observe that there are spikes in the numerical torque sensitivity values for ∂ Q 18 /∂q 14 at t = 0.5 and t = 1 s, but similar spikes are not noticeable for ∂ Q 18 /∂q 14 and ∂ Q 18 /∂q 14 . This indicates that the observed spike is due to the inaccuracy in the numerical estimation of the sensitivity ∂ Q 18 /∂q 14 . Otherwise, similar spikes should have been also observed in the plots of velocity and acceleration sensitivities. To further analyze the numerical behavior, we plotted the trajectories of the generalized coordinates q 5 and q 14 , as well as computed the torques Q 20 and Q 18 in Fig. 6.
One can see that the numerical errors observed in the sensitivity ∂ Q 18 /∂q 14 occur exactly when q 14 = 0 and Q 18 = 0. Thus, numerical errors occur exactly when both the numerator and the denominator vanish. From the trajectories of q 5 and Q 20 , we observe that, when the coordinate q 5 = 0, the torque Q 20 = 0. Similarly, the trajectories ofq 25 and Q 30 demonstrate that at times when the former (angular velocity) is zero, the latter (torque) differs from zero, (see Fig. 6). As expected, the corresponding numerical sensitivity ∂ Q 30 /∂q 25 , Fig. 4, is computed correctly. Thus, the numerical algorithm copes well with the sensitivity calculation when only the denominator vanishes or is close to zero, but it might give erroneous results when both the numerator and the denominator are close to zero. Therefore, numerical evaluation of the sensitivity values close to a singular-ity should be exercised with care. The cumulative time for obtaining all of the analytical sensitivity expressions (measured using tic and toc commands at each sampling time and summed up for all 301 sampling times) was 16.6 ms, whereas the cumulative time for computing the numerical sensitivity was 2.62·10 5 s, amounting roughly to 3 days and 48 min. Thus, the numerical sensitivity analysis takes seven orders of magnitude more time than its analytical counterpart. We expect that for more complex, nonlinear systems with a larger number of DOF the magnitude of errors and the computational time for the numerical sensitivity analysis would be even larger. These results suggest that the numerical sensitivity analysis is error-prone, time-consuming, and computationally expensive. As a result, the numerical sensitivity analysis should be used only as the very last option when other methods are not available for some reason. Our proposed recursive sensitivity analysis is accurate and efficient.

Conclusions
The first-order recursive sensitivity analysis based on the application of the direct differentiation method to the inverse Lagrangian dynamics of rigid multibody systems was presented. We suggested a new description of the system topology based on the branch connectivity graph. The inverse dynamic algorithm based on the Lagrangian formulation was extended to systems with external forces and torques acting at each joint. The advantages of the proposed sensitivity analysis are the following. The method does not suffer from numerical issues associated with the perturbation magnitude in design variables. It has low computational costs and is easy to implement in high-level programming languages such as MATLAB which makes programming and debugging easy. The method is valid for open-chain and open-tree type topological systems. Finally, it is stable, accurate, and numerically efficient. As a result, the proposed algorithm becomes a useful tool in the sensitivity analysis of mechanical multibody systems. Future work could include the incorporation of the divide-and-conquer algorithm into the sensitivity analysis to speed up computation time through parallelization. An extension of our sensitivity analysis for closed- Fig. 4 Results of the sensitivity analysis for ∂ Q20 ∂q 5 , ∂ Q30 ∂q 25 , and ∂ Q42 ∂q22 Fig. 5 Results of the sensitivity analysis for ∂ Q18 ∂q14 , ∂ Q18 ∂q14 , and ∂ Q18 ∂q14 Fig. 6 Time variation of the generalized coordinates q 5 , q 14 , velocityq 25 , and torques Q 18 , Q 20 , Q 30 loop mechanical systems and its application to the human body motion synthesis problem are pending.
Author contributions All authors contributed to writing the manuscript. Simulations were performed by AZ. All authors read and approved the final manuscript.
Funding Open access funding provided by University of Agder. This work has been carried out within the scope of CareWell project funded by the Research Council of Norway under the Grant number 300638.
Data and code availability The MATLAB codes used to generate the results presented in the manuscript are available upon request from the corresponding author.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose. The authors declare that there is no conflict of interests regarding the publication of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.