Higher-order derivatives of rigid body dynamics with application to the dynamic balance of spatial linkages

Dynamic balance eliminates the ﬂuctuating reaction forces and moments induced by high- speed robots that would otherwise cause undesired base vibrations, noise and accuracy loss. Many balancing procedures, such as the addition of counter-rotating inertia wheels, increase the complexity and motor torques. There exist, however, a small set of closed-loop linkages that can be balanced by a speciﬁc design of the links’ mass distribution, poten- tially leading to simpler and cost-effective solutions. Yet, the intricacy of the balance conditions hinder the extension of this set of linkages. Namely, these conditions contain com- plex closed-form kinematic models to express them in minimal coordinates. This paper presents an alternative approach by satisfying all higher-order derivatives of the balance conditions, thus avoiding ﬁnite closed-form kinematic models while providing a full solution for arbitrary linkages. The resulting dynamic balance conditions are linear in the iner- tia parameters such that a null space operation, either numeric or symbolic, yield the full design space. The concept of inertia transfer provides a graphical interpretation to retain intuition. A novel dynamically balanced 3- RSR spatially moving mechanism is presented together with known examples to illustrate the method. an eigenvalue to be solved. The non-linear nature of the eigenvalue a closed-form for arbitrary linkages. study is done on an individual basis. for simple linkages the partitioning and multipole-rod representation interpretation and feasibility study eigenvalues some the linkages design


Introduction
Fluctuating reaction forces and moments generated by fast moving robots cause unwanted base vibrations and accuracy loss at the end-effector [1] . These shaking forces and moments may be reduced or even eliminated by a specific design of the robot's structure and inertia parameters [2] . Such mechanisms, that emit neither shaking forces nor shaking moments, are termed dynamically balanced, or force-balanced when only the shaking forces are zero. We distinguish three major approaches to design mechanisms with this feature. Firstly, one may add supplementary counter-mechanisms to a given mechanism, such as counter-rotating wheels [3,4] or idler loops [5][6][7] . Secondly, various synthesis methods combine and recombine elementary dynamically balanced modules such as four-bar linkages [8,9] or pantograph-like structures [10][11][12] into force balanced or dynamically balanced mechanisms with more degrees of freedom (DOFs). Thirdly, such an elementary module itself is obtained for the analysis of its dynamic balance conditions. By inspecting the equations that describe its motion and dynamics, a range of inertia parameters, i.e. masses, centres of mass (COMs), and moments of inertia (MOIs) may be found that balance this specific mechanism [8,13,14] .
For the viability of dynamic balance, it is essential to find simple and low-weight mechanisms that still fulfil the desired kinematic task. In this view, the addition of counter-mechanisms seems undesirable since it will increases the mass, complexity and the required motor torque. The synthesis approaches, on the other hand, have proven to be versatile for force balance, yet, incomplete for full dynamic balance [15] . In order to expand the scope of dynamic balance and to enable new synthesis methods, the focus of this paper lies in the improvement of the third approach, i.e. the generalization and the automation of the analysis approach.
The necessary and sufficient dynamic balance conditions are a constant linear and angular momentum, as their derivatives, the shaking forces and moments, then will be zero [16] . In practice, when the system is initially at rest, a zero linear and angular momentum suffices. A set of inertia parameters that satisfy these conditions is said to be a dynamically balanced solution, whereas the full description of all solutions is termed the design space of dynamically balanced inertia parameters, or design space for sake of brevity. It should be noted that open-chain linkages cannot be dynamically balanced without additional counter-mechanisms as they require zero or negative moments of inertia [4] . Closed-chain linkages, on the other hand, can in some cases by dynamically balanced by a suitable choice of inertia parameters. However, obtaining the complete design space for these linkages is not trivial as the dynamic balance conditions are to be expressed in minimal coordinates [17] . This involves kinematic loop-closure models, which may be intricate, even for relatively simple linkages [8] , or unavailable in closed-form for more complex mechanisms [18] . The lack of symbolic transparency and closed-form description renders the process of solving the balance conditions an arduous task.
This complexity of the balance conditions is partly overcome by the Linear Independent Vector Method [19] and derived methods [20][21][22] and the Inertia Flow Method [23] . These methods eliminate only a subset of dependent coordinates, leading to simpler balancing conditions while still yielding the complete design space in general. However, in special kinematic cases, such as parallelograms, these incomplete kinematic models lead to spurious force or moment balance conditions and therefore to an incomplete description of the design space [17] . Such special kinematic cases are of particular interest for dynamic balance as they permit more solutions than the general case. For instance, Ricard and Gosselin [8] showed that the kite-type and antiparallogram-type of the planar 4 R four-bar linkage may be fully dynamically balanced by a specific mass distribution. This in contrast to the general four-bar linkage that does require additional counter-rotating measures.
To prove the work of [8] formally, Moore et al [24] . factorized the balance conditions and loop-closure equations by means of toric geometry and Laurent polynomials. Subsequently they showed through a similar algebraic approach that the spherical four-bar linkage cannot be dynamically balanced without additional structures [25] . Currently, these algebraic methods still require a tailored inspection per mechanism and are yet to be extended to multi-DOF mechanisms. An alternative method to deal with the kinematic complexity of the loop-closure equations was adopted in [26] . There, screw theory was applied to find instantaneous dynamic balance, yielding a single pose in which mechanism accelerations will not induce shaking forces and moments. Since outside this pose the balance quality is not guaranteed, this method yields and solves only the necessary but not sufficient conditions for dynamic balance.
To summarize; in literature several systematic analysis methods were presented that solve the dynamic balance conditions for given linkages. Yet, no method yields the complete dynamically balanced design space of arbitrary linkages without a tailored manipulation of the loop-closure equations. Such a method is desired to advance our understanding of dynamic balance and to find new, simple and lightweight balance solutions.
In this paper, this long-standing problem is tackled by extending the instantaneous approach of [26] over the complete workspace by including and solving a sufficient number of higher-order derivatives of the dynamic balance conditions. The higher-order kinematic and dynamic models are readily available through recursive application of the implicit function theorem [27] , thus avoiding the use of closed-form kinematic models. This method leads to the necessary and sufficient dynamic balance conditions, and an automatic and complete characterization of all dynamically balanced designs of any given nonsingular mechanism consisting of lower kinematic pairs. To that end, this paper for the first time presents an algorithm to compute the derivatives of the bodies' mass matrices and momentum equations in open-and closed-loop linkages up to arbitrary order. This paper starts with a synopsis of the method to guide the reader through the following sections ( Section 2 ). Thereafter the higher-order derivatives of kinematics is outlined ( Section 3 ), followed by a recapitulation of the rigid body dynamics in the screw theory framework ( Section 4 ). This leads to a recursive algorithm that yields the higher-order derivatives of the linear and angular momentum equations (the dynamic balancing conditions) of open and closed-chain linkages ( Section 5 ). The resulting higher-order momentum equations are then recast into the parameter-linear form [28,29] to provide dynamic balance conditions that are linear in the inertia parameters and solvable by null space algorithms ( Section 6 ). General null space algorithms, i.e. singular value decomposition or Gaussian elimination, yield a complete description of all dynamically balanced mass distributions. This description however, is strongly mixed in the inertia parameters, causing a loss of interpretation and design intuition. Therefore, an alternative, meaningful description of the design space of open-chain and closed-chain linkages is presented ( Section 7 ). This description is derived from the concept of inertia transfer and the multipole representation of the inertia parameters , as used in the parameter identification of robots [30] . This interpretation, here termed the multipole-rod representation, is shown to aid the feasibility study of dynamic balanced linkages. It should be noted that open chains receive quite some attention in this work. Although they cannot be dynamically balanced themselves, they provide insight into the solution space of closed-loop linkages. More specifically, it will be shown that a large portion of the design space of a closed-loop linkage is build up from the open-chain equivalents into which the linkage may be decomposed. Case studies of a 6-DOF serial robot, a 4 R planar four-bar linkage, and a 3-RSR mechanism illustrate the higher-order dynamic balance method ( Section 8 ). This results in a novel 3-RSR mechanism design that is dynamically balanced for the 2-DOF that lie on three planes of mirror symmetry. Refer to Table A.1 for a list of symbols used in this paper.

Synopsis of the higher-order dynamic balance method
Ground-based open-chain linkages are dynamically balanced if the momentum h is zero for all n b joint coordinates q and all joint velocities ˙ q . Note that this h is a combination of the linear and angular momentum and thus a 6-dimensional vector. Since the momentum must be zero for all joint velocities and since these joint velocities are linear in the momentum we obtain the following balancing condition } denote the collection of the basis vectors of h with respect to ˙ q . z denotes the collection of the inertia parameters of the linkage, i.e. masses, centres of mass, and moments of inertia. The aim of this paper now lies in the derivation of the complete set of inertia parameters z that guarantee dynamic balance for any given open-or closed-chain linkage.
For closed-loop linkages the joint coordinates q are no longer independent since a set of loop closure constraint equations hold for all motion This leads to dependencies in q and ˙ q , and consequently, in a reduced set of balancing conditions. Conventionally, this dependency would be incorporated into Eq. 1 by selecting a set of minimal coordinates u and solving Eq. 2 for the dependent coordinates q = c( u ) . However, this approach is not always applicable as there is in general no closed form solution to the loop-closure equation, i.e. c is not always known explicitly. Furthermore, if the solution is found nevertheless, it is typically a set of involved equations that are hard to use in the balancing procedure.
In this paper we take a different, Taylor-based approach. It relies on three features. Firstly we leverage the fact that, although c might not be available in closed form, its higher-order derivatives D k u ( c ) are available in the reference configuration u 0 through a recursive application of the implicit function theorem ( [27] and recapitulated in Section 3 ). The resulting higher-order derivatives enable a Taylor expansion of the dynamic balance conditions such that with a slight abuse of notation these read As this must hold for all motion u , all Taylor coefficients (the higher partial derivatives D k u h ) are required to be zero for dynamic balance.
Since h and f are analytic functions in a non-singular configuration, we obtain the following necessary and sufficient conditions Note that this does not require an explicit solution to the loop-closure equations ( Section 5 ). Furthermore since the balance conditions are analytic, only a finite (but unknown) number k max of partial derivatives is sufficient to ensure dynamic balance, enabling an algorithmic treatment of this problem. The second feature is that these Taylor coefficients ( Eq. 4 ) are linear in inertia parameters z ( Section 4 ). The balance conditions can therefore, with the help of the regression matrix X k , be written in the form (5) and consecutively solved by null space algorithms ( Section 6 ), leading to a full description of the dynamically balanced mass distributions z ∈ ker ( X k ) , provided that a sufficient number of derivatives is used.
Thirdly, in this paper we present a systematic partitioning and interpretation of the resulting design space in order to retain some insight in of the design dependencies and the feasibility of the solution ( Section 7 ). We illustrate the method on known and new examples ( Section 8 ).

Kinematics
In this section the groundwork of the method is laid by describing the notation and kinematics, and by recapitulating the higher-order derivatives of kinematics.

Kinematics of open-chain and closed-chain linkages using lie group and screw theory
Screw theory is used throughout this work as it gives a concise representation of the higher-order derivatives of kinematics and dynamics. This screw theory framework interprets motion of a body as a combination of an angular velocity ω around an axis in space, passing through point r t , and a velocity along that axis, termed the pitch λ t The twist t is a function of the angular velocity ω and the velocity v of the body's particles that instantaneously pass through the origin of the reference frame. Two special cases exist: 1) a pure rotation, i.e. the angular velocity is orthogonal to the velocity, resulting in a zero pitch ( λ t = 0 ), and 2) a pure translation, when the angular velocity is zero and the pitch is infinite ( λ t = ∞ ).
The twist of the n b bodies in a open chain is linearly dependent on the joint velocities ˙ q of the joints in the chain. The twist basis vector s j associated to each joint j is termed unit twists or instantaneous screw axis (ISA). In the current context s j is always taken in the reference configuration q 0 . The joints in a single chain are numbered 1 to n b , from the base to the end-effector. Therefore, the Jacobian J i of body i is formed by the ISAs lower in the chain These ISA are pure geometric quantities that solely dependent on the joint location r s , the orientation of the joint, encoded by unit vector n , and the pitch of the joint λ s . In this case we treat the three single DOF lower kinematic pairs: revolute R , helical H , or prismatic P . For an R -joint the pitch is zero λ s = 0 , while for a P -joint the pitch is infinite λ s = ∞ , resulting in a limit case s ∞ and is therefore treated separately. Multi DOF joints, such as ball-socket joints, are treated as instantaneously identical to a set of serially connected single DOF joints.
In closed-chain linkages, the body twists are restricted by a set of loop closure conditions f . The resulting twists may be found by regarding each loop as a connection of multiple open chains. A single loop for example is opened by cutting an arbitrary body, resulting in two open chains of which the last 'virtual' bodies have the same twists. These loop-closure conditions constrain the twists of the bodies, as encoded by matrix K . By selecting an independent set of n d input coordinates u , this system is solved and all dependent joint velocities determined The linkage Jacobian J = [ J 1 · · · J n b ] is the collection of all n b body Jacobians. The n b × n d C -matrix denotes the first-order solution to the loop-closure equations.
To express finite motion, a reference frame ψ i is associated to each body i . A homogeneous transformation matrix H i , consisting of a rotation matrix R i and a translation vector o i , express a point a from a body-fixed frame into the inertial frame of reference a 0 = H i a i . In this convention the a i vector consist of four values; 3 Cartesian coordinates and a 1. The superscripts denote the frames of expression of the vector. These transformation matrices relate again to the ISAs in a open chain through a product of matrix exponentials, leading to the general forward kinematics of open-chain linkages [31] in which exp (q j s j × ) denotes the matrix exponential of the 4 × 4 matrix of the ISA in the reference (initial) configuration ( q = q 0 = 0 ) , and n × the 3 × 3 skew symmetric matrix of n .
The ISA are expressed in another coordinate frame by the adjoint transformation matrix Ad H . The ISA is expressed from the body fixed reference frame s i i in the inertial frame of reference s 0 i according to The time derivative of the transformation matrix relates to the body twist through the matrix form of the adjoint twist transformation, here termed adjoint twist matrix ad( t i ) d dt

Higher-order derivatives of kinematics
For parallel mechanism a closed-form solution to the kinematic loop-closure equations does not exist in general. Yet, a higher-order approximation of the motion is available by treating the closed loop as a connection of several open chains. For such a connection, the higher-order derivatives of the loop-closure equations are found and solved yielding a Taylor approximation of finite motion [27] . In that approach, the higher-order partial derivatives of the body twists are found from the adjoint twist matrices corresponding to the ISA that are lower in the open-chain equivalent linkage [32] . Since each ISA is constant when expressed in a local body-fixed frame, all these derivatives follow from a repetitive application of Eq. 11 to Eq. 10 , such that D α In here D α denotes the higher-order partial derivatives with respect to the elements of q . Vector α = (α 1 , . . . , α n ) comprises the order of the derivatives corresponding to q , running from the base to the end-effector. Hence we assume an ordered sequence, i.e. α i corresponds to the joint q i . The k = α 1 + . . . + α n = | α| is the total order, see Appendix A . The joints higher in the chain have no contribution to the motion of the lower joints, such that this derivative ( Eq. 12 ) is set to zero, i.e. if α j = 0 for j ≥ i . By this, all the higher-order partial derivatives of the body Jacobians D α This procedure is used for the solution of the higher-order closed-loop constraints [27] by recasting it into the matrix derivative framework of Vetter for bookkeeping [33] and Appendix A . In this notation the collection of all first-order partial derivatives of matrix A = [ a 1 · · · a m ] are sorted according to 1 (13) With this, the derivatives of the loop-closure solution D u ( C ) are found through application of the chain rule and product rule ( Appendix A ) to Eq. 8 . The collection of second-order loop-closure constraints read In here A B denotes the Kronecker product of two matrices ( Appendix A ). From this equation D u ( C ) is determined. A recursive application leads to the k -th order constraints from which D k u ( C ) may be solved through the algorithm presented in [27] . The Kronecker power is denoted by a k superscript. The exact composition of the C k collection matrix is found through repetitive application of the chain and product rules, but is omitted here due to space limitation.

Rigid body dynamics
The rigid body dynamics of spatially moving objects and mechanisms is concisely written with the use of screw theory [34,35] . This section briefly introduces the use of screw and Lie group theory in rigid body dynamics, followed by the presentation of the multipole-rod representation of the inertia parameters as used in the interpretation of the dynamically balanced solution later on.

Momentum wrench and mass matrix
The momentum of a body is the product of the body's spatial mass matrix M and the twist t associated to it. The momentum is a co-screw or a wrench-like entity and therefore termed momentum wrench hereafter The mass matrix of a body is formed by the integral over the body volume This gives rise to the classical description with a mass m , a centre of mass c and inertia matrix E with respect to the inertial frame of reference. The inertia matrix E contains 3 inertia moments and 3 products of inertia, respectively on its diagonal e d = [ e 1 e 2 e 3 ] and its off diagonal e o = [ e 4 e 5 e 6 ] . The matrix I 3 denotes a 3 × 3 identity matrix. Due to the frame invariance of kinetic energy K = 1 / 2 t M t , the mass matrix transforms with an adjoint transformation matrix on the right and its transposed on the left. By choosing a frame that is located at the centre of mass and aligned with the principal axis The multipole representation [30] with parameters that are linear in the mass matrix; a monopole m at r , a dipole δ in the direction a , and a quadripole η in the direction of b . One monopole, three dipoles and six quadripoles are sufficient to describe arbitrary bodies. (c) The multipolerod representation reduces the number of graphical elements by interpreting the quadripole as an infinitely long, slender rod, termed 'pure-inertia rod' and depicted as a striped bar. The monopole is termed 'point mass', whereas the dipole is treated as a 'displacement' of the point mass with negative pure-inertia rod in the same direction.
of inertia, any mass matrix can be diagonalized. The corresponding transformation matrix from this principal axes frame to the current frame is Ad H p . This gives rise to three principal MOIs g 1 , g 2 , and g 3 In this body-fixed frame, the mass matrix is constant, i.e. ˙ m and ˙ g i = 0 , due to the rigid body assumption. Since the mass matrix is formed by a collection of positive mass particles, the mass matrix itself is symmetric positive definite, leading to 7 inequality conditions on the mass and the principal MOIs

Momentum wrench basis
Similar to the twist basis, we define a linkage's momentum basis that spans all possible momentum wrenches at a given pose. The basis vectors, termed the instantaneous momentum wrenches (IMW) and denoted with ˆ h i , are the momentum wrenches generated by unit actuation of each joint. The total momentum wrench of a open chain is therefore given by In here M = [ M 1 · · · M n ] denotes the collection of all mass matrices in the chain. For dynamic balance all the IMWs must be zero for arbitrary motion. For closed-chain linkages the momentum wrench basis is computed by applying the first order loop-closure solution C

Multipole-rod interpretation of the mass matrix
In the current dynamic balancing procedure we will use the fact that the balancing conditions are linear in the elements of the mass matrix such that they can be solved through a set of linear operations. The conventional mass matrix parametrisation, consisting of masses m , COMs c and principal MOIs g , is not suitable for the interpretation of the resulting design space, since it is not linear in the elements of the mass matrix. Therefore we will use a slight adaptation of the multipole concept of Ros et al [30] ., termed the multipole-rod representation ( Fig. 1 ). This interpretation relies on the fact that a mass matrix can be decomposed into three primitive elements; 1) a single point mass at r , denoted with a subscript m , 2) a displacement of the point mass in the direction of a unit vector a combined with a pure-inertia rod of opposite magnitude, denoted with a subscript δ, and 3) a pure-inertia rod in the direction of a unit vector b , denoted with a subscript η. These pure-inertia rods are interpreted as infinitely long slender rods in the direction of their unit vector. Their mass is assumed zero such that only the rotational velocity component in a perpendicular direction generates angular momentum. A rotation around their longitudinal axis generates no angular momentum. The sole difference with the method of [30] is the graphical representation. This reduces the larger number of point masses (poles), which otherwise might congest the figures. Now, any mass matrix can be represented by choice of 10 of these primitive elements, one point mass, three displacements, and six pure-inertia rods, as long as the unit vectors a i and b i are unique The mass, the moment of mass, and the moment of inertia of these elements are given by For the planar case, this representation requires one point mass, two displacements and one pure-inertia rod, of which the elements reduce to For feasibility of each body, they must consist of at least one positive point mass, and three non-coplanar positive pureinertia rods ( Eq. 19 ), since two pure-inertia rods represent an infinitely flat object. A negative pure-inertia rod requires at least 3 arbitrarily oriented positive pure-inertia rods (or two positive coplanar pure-inertia rods) of sufficient magnitude to represent a feasible body. A closed-form feasibility description of an arbitrary collection of these elements can be found through eigendecomposition of the resulting mass matrix, but lies outside the scope of this paper.

Higher-order derivatives of the momentum equations and of the dynamic balance conditions
The previously presented higher-order analysis of the kinematics is extended to rigid body dynamics in this section. The aim is to find and solve the necessary and sufficient dynamic balance conditions of arbitrary linkages without invoking the closed-form solution to the loop-closure equations. For dynamic balancing purposes this study is confined to the change of rigid body momentum. Other effects such as gravity, elasticity, or external forces are not taken into account. The dynamic balance conditions are obtained from the partial derivatives of the mass matrices and momentum equations of open-chain linkages, which are extended thereafter to closed-chain linkages by including the higher-order derivatives of the loop-closure solution. It should be noted that although open-chain linkages cannot be dynamically balanced without additional countermechanisms, their description is important for dynamic balance since closed-loop linkages can be regarded as connected open chains.

Derivatives of the mass matrix in a open chain
The mass matrix of a body i in a open chain depends on the pose of the joints that are lower in the kinematic chain according to Eq. 9 and Eq. 18 . Therefore, its partial derivative with respect to a joint j , lower in the chain ( j ≤ i ), is found by applying Eq. 11 to Eq. 18 Here we have used the fact that the mass matrix is constant in the body-fixed frame. For all partial derivatives with respect to joints higher in the chain ( j > i ) this derivative is zero. A second (non-zero) partial derivative is either with respect to a joint higher ( j ≤ l ≤ i ) or joint lower ( l ≤ j ≤ i ) in the chain. In the first case ( j ≤ l ≤ i ) the partial derivative becomes Here the Jacobi identity only the derivative of the mass matrix has to be taken into account as a higher joint does not influence a lower ISA . This results in a similar equation as Eq. 26 , with the sole difference that the indices j and l are swapped.
This also follows from the symmetry of partial derivatives. This nested structure, i.e. the pre-and postmultiplication of adjoint twist matrices, is preserved for the higher orders, leading to a recursive formula for all partial derivatives of the mass matrix in here j is the lowest joint to which a partial derivative is taken, i.e. α l = 0 for all l < j . In case α l = 0 for any l > i , this equation is set to zero.

Derivatives of the momentum wrench in a open chain
Now that the derivatives of the mass matrix up to arbitrary order are available, we consider the partial derivatives of the momentum wrench with the aim of obtaining all higher-order dynamic balance conditions. Consider the momentum wrench generated by the j th body due to unit actuation of joint i , which is lower in the chain. Two types of non-zero partial derivatives appear. Either joint l -with respect to which the derivative is taken -is below the joint i , or between the joint i and the j th body. In the first case ( l ≤ i ≤ j ), the partial derivative of both the mass matrix and the ISA have to be taken into account, partially canceling out In the second case ( i < l ≤ j ), the partial derivative of the ISA vanishes ∂ / ∂q l ( s i ) = 0 . Therefore, the partial derivative of the momentum wrench becomes The higher-order partial derivatives are found similarly by making a split between the partial derivatives related to joints lower than the momentum generating ISA, and the ones related to the joints between the ISA and the body. Therefore. a second multi-index is introduced for which holds β l = α l for all i < l ≤ j and β l = 0 for all l ≤ i . The partial derivatives of the momentum wrench are found from Eq. 27 according to Again this equation is zero if α l = 0 for any l > j . These partial derivatives may be summed to obtain the derivatives of the total momentum of the linkage. Notice that in this equation the momentum derivatives are formulated as a sequence of matrix operations, which are linear in the mass matrix.

Derivatives of the dynamic balance conditions of a open-chain linkages
The dynamic balance conditions dictate that the momentum wrench of a linkage is zero for all motion. Therefore also all higher-order derivatives of the momentum wrench must be zero. With a large enough number of derivatives k m ax these are not only the necessary but also the sufficient dynamic balance conditions for nonsingular linkages. In fact, here it will be shown that for open-chain linkages only derivatives up to the second order are needed ( k max ≤ 2). When these are satisfied, all the higher-order dynamic balance conditions satisfied, and the linkage is dynamically balanced for finite motion.
For zeroth-order dynamic balance, the condition ( Eq. 20 ) imposed on each IMW is The aggregated mass matrix ˜ M i is the sum of the mass matrices belonging to bodies higher in the chain than s i . Consider now the following momentum derivatives of ˆ h j and ˆ h l , involving any triplet s l , s j , and s i of zero or finite pitch ISA, which are arranged in ascending order ( l ≤ j ≤ i ) Notice that these dynamic balancing conditions impose constraints on the same aggregated mass matrix ˜ M i since q i is higher in the chain than q l and q j such that ∂ / ∂q i ( M j ) = 0 for j ≤ i . As the first-order balancing conditions ( Eq. 32 .b) ensure that A recursive application shows that this extends to the higher orders, such that all balance conditions are of the form Moreover, the zeroth-order balance conditions ( Eq. 31 ) satisfies all higher-order force balancing conditions since is a function of the linear momentum and the mass is assumed to be constant Therefore, only the following first-and second-order moment balance conditions remain: In the general case, when n j ࢲ n l , this imposes 9 independent constraints on the derivative of the inertia matrix, requiring , thus directly satisfying all higher-order partial derivatives ( Eq. 35 ). This shows that derivatives of a higher order than k max = 2 impose no new dynamic balance conditions for open-chain linkages. When, in the special case, all noninfinite pitch ISA lower in the chain are parallel, i.e. n j n i for all j < i , the moment balance conditions ( Eq. 37 ) vanish or become equivalent. Then, only three higher-order constraints are imposed on the aggregated inertia matrix ˜ E i . Prismatic joints (infinite pitch ISA) lower in the chain impose no higher-order moment balance conditions as their angular velocities n j are zero.
To summarize: for open-chain linkages the zero-order force and moment balance conditions ( Eq. 31 ) and the first-and second-order moment balance conditions ( Eq. 37 ) are necessary and sufficient, leading to a k max = 2 .

Derivatives of the dynamic balance conditions of closed-chain linkages
The dynamic balance conditions of closed-chain linkages dictate a zero momentum wrench ( Eq. 21 ) for all independent velocities ˙ u . Therefore the zeroth-order balancing conditions read Also all higher-order partial derivatives with respect to u should be zero for dynamic balance. These conditions are found by repetitive application of the chain rule, the product rule and derivatives of the Kronecker product. Similar to Eq. 15 , the first-order dynamic balancing conditions become This generalizes to higher-orders by a repetitive application of the chain and product rules From the analyticity of the momentum equations it may be deduced that there is finite k max which renders these conditions not only necessary but also sufficient for the dynamic balance for closed chains in nonsingular poses. Refer to Section 9 for a discussion on the necessity and sufficiency of these conditions.
It should be noted that these higher-order dynamic balance conditions are linear in the mass matrices and can be obtained through a series of matrix multiplications and linear operations. This method is therefore able to treat symbolic or numerical input.

Dynamic balance solution using the parameter-linear form
Now, to solve these higher-order dynamic balance conditions, we recast Eq. 40 in the parameter-linear form [28,29] as used in the parameter identification. This enables null space procedures to extract the dynamically balanced mass distributions.

Parameter-linear form
Since the m, m c and E ( Eq. 17 ) are linear in the momentum equation, the following parameter-linear form holds in which the z -vector is formed by the inertia parameters of the body. The twist dependent 'regression' matrix is given by Notice that the ordering of the inertia parameter slightly differs from [29] . The parameter-linear form of the momentum basis of a open-chain linkages is directly computed from Eq. 20 in here h and z denote the concatenation of all IWM and all inertia parameters in the chain, respectively. To obtain the parameter-linear form of closed-chain linkages, the vectorization of matrix products ( Appendix A ) is applied to Eq. 38 where I 6 is a 6 × 6 identity matrix.

Higher-order dynamic balance conditions in the parameter-linear form
The parameter-linear form also applies to higher-order derivatives of the balance conditions as they are formed through a sequence of matrix operations that are linear the inertia parameters. The higher-order open chain regression matrices W k can be found accordingly, i.e. by the application of Eq. 41 to Eq. 27 and Eq. 30 , resulting in the following condition vec D k For closed chains the parameter-linear form is found by applying the matrix vectorization to Eq. 40 , such that vec D k in which W k = W 1 · · · W k . Now we have arrived at the parameter-linear form of the higher-order derivatives of the momentum equations of open-and closed-chain linkages. It should be observed that all these steps solely rely on matrix operations suitable for algorithmic treatment.

Solving the dynamic balance condition
Dynamic balance requires inertia parameters z that are on the intersection of the null spaces of all the X i matrices z ∈ ker ( X i ) , z ∈ ker ( X k max ) , z = N y (47) in which X k max = X 1 · · · X k max is the collection of all regression matrices up to order k max . It should be emphasised that there is a finite k max , which makes the approach practically feasible. The columns of the N matrix form a basis that span this null space and therewith describe the full design space of the dynamically balanced inertia parameters. This N matrix is termed the design space matrix and may be found through numeric or symbolic null space algorithms such as Gauss-Jordan elimination or singular value decomposition. The corresponding design parameters are collected in y . With this the complete set of dynamically balanced inertia parameters of any given nonsingular linkage may be found.

Partitioning and interpretation of the dynamic balance solution
The application of null space algorithms to the dynamic balance problem ( Eq. 47 ) may result in a design space description that is strongly mixed in the inertia parameters, compromising structure and design intuition. To aid the designer, a partitioning of the design space with respect to the joint topology is presented alongside a multipole-rod representation ( Fig. 1 ) of these partitions. We shall show that 6 types of inertia transfer matrices completely describe the design space of open-chain linkages. These inertia transfer matrices contain all inertia parameters that may be transferred between two hinged bodies, i.e. subtracted from one body and added to the other, without changing the momentum generated by the linkage. This partitioning will lead to a general description of the design space of open-chain linkages that, more importantly, also covers a large part of the design space of closed-loop linkages. Closed-loop linkages may be seen as a connection of multiple open chains. A balancing solution that is valid for open-chain linkages is therefore also valid for closed-chain linkages. Although the open-chain design space itself is always unfeasible, in combination with a closed-chain design space it allows for more feasible solutions as shown later in the examples.

Partitioning the design space of open-chain linkages
The dynamic balancing conditions of open-chain linkages ( Eq. 31 and Eq. 37 ) are formulated in terms of aggregated mass matrices ˜ M i . Before presenting the general solution it may already be observed that solution to these equations will also be in terms the aggregated mass matrices. From these aggregated solutions each individual mass matrix can be found accordingly Therefore, the complete design space matrix N of an open chain ( Eq. 47 ) may be partitioned as a band diagonal matrix in here the submatrix N i describes all inertia parameters that can be exchanged between the two bodies connected by joint i without changing the dynamic behavior of the chain. These N i submatrices are therefore termed inertia transfer matrices. In Section 7.3 it is shown that there exist actually 6 types of inertia transfer matrices depending on the type of joint and parallelism with the joint axes lower in the chain. Table 1 The dimensions of the 6 inertia transfer matrices. Each joint i in a chain extends the design space depending on the type of joint; revolute ( R ), helical ( H ), or prismatic ( P ) and the alignment with all non-prismatic joints j < i lower in the chain; a) skew or b) parallel. * With a prismatic joint, the prismatic joint direction applies n i = n i, ∞ .

Joint type
It should be noted that a similar concept is used in the context of parameter identification to describe the set of unidentifiable inertia parameters [30,36] . Broadly speaking, inertia parameters are said to be unidentifiable if they do not contribute to the kinetic energy of the linkage. The dynamically balanced design space of open-chain linkages, as found here, is formed by unidentifiable inertia parameters as zero momentum in this case also implies zero kinetic energy. The inverse is not true in general. This also shows that the inertia parameters in this design space do not affect the required motor effort of the linkage. (51)

Interpretation of the design space via the concept of inertia transfer
In Section 5.3 , it was shown that dynamic balance imposes two conditions on the aggregated mass matrices of openchain linkages: Firstly, each aggregated mass matrix ˜ M i should be chosen such that its IMW vanishes ( ˜ M i s i ≡ 0 ). Secondly, the actuation of the corresponding joint q i should not change the angular momentum generated by aany joint lower in the From the first condition three cases arise; an ISA of zero, finite or infinite pitch, while for the second condition two cases exist; either all axes up to n i are parallel ( n j n i for all j < i ) or at least one is skew ( n j ࢲ n i for j < i ). This gives rise to 6 types of design space for 1-DOF lower kinematic pairs, and, consequently, 6 types of inertia transfer matrices N i ( Eq. 49 ). These are discussed now. For higher-DOF joints and joints in planar linkages a similar representation exist as shown subsequently.

Six inertia transfer matrices
Here, the multipole-rod representation of these six inertia transfer matrices are given ( Fig. 2 ). In this notation the point mass, displacement and pure-inertia rod elements of the multipole-rod representation ( Eq. 23 ) are respectively denoted by z m ( r ), z δ ( r , n ), and z η ( n ). The dimensions of these inertia transfer matrices are in Table 1 . Starting from a revolute joint, whose axis has no particular alignment, the six cases are discussed and interpreted. N 0, ࢲ The inertia transfer matrix associated to a revolute joint ( λ = 0 ) -whose joint axis is skew ( ࢲ ) with respect to one or more preceding revolute or helical joints -comprise of three inertia parameters. These three parameters can be freely exchanged (added to one and subtracted from the other) between the two bodies hinged by this joint without Fig. 2. The interpretation of the six sets of inertia parameters that can be exchanged between the two (grey) bodies attached to joint i (subtracted from one and added to the other) without changing the dynamic behaviour of the chain as a whole. These six cases arise from the three types of 1-DOF lower pairs, and parallelism with all preceding revolute and helical joints. The orientation of the preceding prismatic joints have no influence. It should be noted that for clarity sake the effect of the displacement δ on the MOIs is not shown, as it can be compensated by, or absorbed in η 1 . Since the pure-inertia rods have no application point, they are displayed at an arbitrary location. affecting the momentum generated by the chain as a whole. These parameters are: 1) a point mass z m on the joint axis r s , 2) a displacement of this point z δ in the direction of the joint axis n , 3) a pure-inertia rod η in the direction of the joint axis n . The corresponding inertia transfer matrix therefore reads The reason for these three inertia transfers is that the actuation of a joint with a point mass m anywhere on its axis r s does not induce any linear or angular momentum, nor does it change the IWM of lower joints ( Eq. 37 ), since and equal and opposite point mass is attached to the connecting body. This yields the design freedoms z m and z δ . The third design freedom, a pure-inertia rod η, generates no angular momentum as it is aligned with the joint axis. This alignment also makes sure that the rotation of this pure-inertia rod by the joint will not cause a change in the inertia matrix felt by the lower joints. Any other exchange of mass or inertia between the two bodies connected by this joint will either change the momentum generated by this joint or by the joints lower in the chain. N 0, When the revolute joint ( λ = 0 ) is parallel with respect to all preceding revolute and helical joints, two additional parameters are obtained, in comparison to N 0, ࢲ . These parameters are two perpendicular pairs of pure-inertia rods. All these all four rods are on a single plane perpendicular to n . These pure-inertia rod are of opposite magnitude in a pair wise manner ( Fig. 2 ). These four additional pure-inertia rods allow for a modification of the inertia tensor without changing the dynamics of the chain. The first of the pure-inertia rod pairs η 2 is in the direction of b 2 , which is perpendicular to n . The angular momentum induced by b 2 is cancelled by an equal and negative pure-inertia rod in a direction perpendicular to both n and b 2 . This also holds for a second pair η 3 with corresponding b 3 . This additional pure-inertia rods arise since their common plane which is perpendicular to n is not changing by actuation of the joints lower in the chain. The inertia transfer matrix is therefore parametrized according to (53) N f, ࢲ For a helical ( λ = finite ), non-parallel joint any point mass will generate a linear momentum through its pitching motion such that its inertia transfer matrix only contains a pure-inertia rod in the direction of the joint axis. The displacement ( z δ ) would cause a pose dependent inertia matrix and non-constant IWMs associated to the lower joints. The sole inertia change is therefore N f, When a helical joint is parallel to all preceding revolute and parallel joints it has a similar inertia transfer space as N 0, ( Eq. 53 ) with the sole difference that the mass should therefore equate to zero as the pitching motion would generate a linear momentum. The displacement ( z δ ) on the other hand does not induce linear momentum and therefore remains The inertia transfer of a prismatic joint ( λ = ∞ ) whose joint axis is not aligned with all preceding revolute or helical joint axes has a size 6, since its E can be selected freely. As these moments and products of inertia will not induce angular momentum (and are constant) they can be selected as desired. Here this choice is parameterized by 6 pure-inertia rods (56) N ∞ , When the prismatic joint ( λ = ∞ ) is aligned with preceding zero and finite pitch joints it gains a displacement z δ in the direction of the joint axis, leading to an inertia transfer matrix with size 7 (57)

Multi-DOF joints
This approach also holds for multi-DOF joints that can locally be modelled as a serial connection of 1-DOF joints, e.g., cylindrical, planar, universal or spherical joints. These multi-DOF joints can the transmit inertia parameters that are common in the lower kinematic pair analogue. For example, a cylindrical joint can be modelled as prismatic and a revolute joint in series such that its inertia transfer matrix is the intersection of the image of N 0 and N ∞ , which is the inertia transfer associated to a helical joint N f . Depending on the parallelism with the joints lower in the chain, either type is to be selected. A planar joint is serial connection of two prismatic joints such that its inertia transfer is N ∞ . Universal and spherical joints locally behave as a serial connection of multiple non-parallel, intersecting revolute joints. The associated inertia transfer is therefore a point mass z m ( r ) on the intersection point r of these axes. We assume here that the joint itself does not contain intermediate bodies with mass or inertia.

Joints in planar linkages
In the planar case only zero and infinite pitch ISA exist. Therefore, three types of inertia transfer matrices appear, 1) the ISA is revolute and therefore automatically parallel to all other revolute joints, 2) the ISA and all lower ISA are prismatic joints or 3) the ISA is prismatic but at least one of the lower joints is a revolute joint In the first case, the point mass should be on the revolute joint and the MOI around that point should be zero. In the second case, the body solely translates, therefore the mass should be zero and the first and second moments of mass are free, as parameterized by a displacement and a pure-inertia rod. In the third case, when the ISA under inspection is prismatic and one or more lower joints are revolute, this displacement will causes a changing MOI associated to the rotation of the lower joints. This the displacement should therefore be zero.
With this description of the inertia transfer matrices N i of common joints, the dynamically balanced design space of any open-chain linkage may be obtained. Also for closed-chain linkages, the open-chain equivalent design space matrix N O is completely determined, generalizing [23] to spatial linkages. It should be noted that open linkages cannot be dynamically balanced without addition of counter-mechanisms. This can also be established from the inertia transfer matrices as none of them permit both a positive mass and positive MOIs. Therefore, a specific N C is required to render a feasible dynamically balanced design space. The existence of this additional design space is found on a case-by-case basis in the next section.

Case studies
The higher-order dynamic balance approach is illustrated here with case studies of a serial 6-DOF robot, a planar 4 R four-bar linkage, and a 3-RSR mechanism. In all cases an interpretation of the closed-chain design space bases N C will be given, if present.

Serial 6-DOF robot
Here we study a 6-DOF robot to show how the open-chain design space matrix is found. We consider a HPRS -robot ( Fig. 3 ) of which the first joint axis remains parallel to the third joint axis n 1 = n 3 . The second, prismatic joint has no particular alignment. The fourth, spherical joint is modelled as three intersecting twists s 4 , s 5 , and s 6 with no intermediate bodies. Therefore this linkages consists of four bodies in total. The zeroth-order regression matrix and the associated design space matrix become This zeroth-order regression matrix has rank ( W 0 ) = 23 . Through application of the algorithm presented in Section 6 , it is found that the regression matrix terminates at k max = 1 with rank ( W 1 ) = 24 leaving n b × 10 − rank ( W 1 ) = 16 design parameters free.
This can be explained by further investigation of W 1 . It is found that only one partial derivative ( ∂ / ∂q 2 ( ˆ h 1 ) = ∂ / ∂q 2 ( ˜ M 2 ) s 1 ) leads to additional dynamic balancing conditions in comparison to W 0 . All other derivatives are 'covered' by the zeroth-order conditions. Therefore, the only relevant rows of the first-order regression matrix are It should be noted that the first 10 columns are zero, due to ∂ / ∂q 2 ( M 1 ) = 0 . Since s 2 is not aligned with s 1 , the actuation of q 2 will cause a displacement of the aggregated COM of bodies 2, 3, and 4, leading to a change of the aggregated mass matrix felt by s 1 . As the other joints are either aligned ( n 1 = n 3 ) or spherical ( S , Section 7.3.2 ), no additional constraints are placed. The corresponding inertia transfer matrices in Eq. 59 are therefore The first joint is a helical joint fixed to the base such that it inertia transfer matrix has rank ( N f , ) = 4 ( Eq. 55 ). The second joint is a prismatic joint whose axis is not parallel to joint 1 such that is inertia transfer matrix has rank ( N ∞ , ∦ ) = 6 ( Eq. 56 ).
In fact, if only the zeroth-order balancing condition would be taken into account N 2 would gain an additional 'displacement' and it would be N ∞ , . The revolute joint, joint 3, is parallel and stays parallel to the preceding finite pitch ISA (joint 1) such that it has an inertia transfer matrix of rank ( N 0 , ) = 5 ( Eq. 53 ). The spherical joint is seen as the intersection of three perpendicular revolute joints such that it solely transmits a point mass at the intersection point z m ( r 4 ). This leads to a total design freedom of rank ( N ) = 16 , which again is in agreement with the order of the regression matrix W k max that is found with the presented algorithm. This shows that the complete design space is found through a combination of the inertia transfer matrices that are readily available by inspecting the joint types and alignment of these joints in the chain.

Planar four-bar linkage
The planar 4 R four-bar linkage is a classical example in the dynamic balance literature. Although studied for decades [3,37,38] , a solution that does not rely on additional counter-mechanisms was not found until the year 20 0 0 by Ricard and Gosselin [8] . Thereafter, this linkage has been used as dynamic balancing module to obtain 3-DOF and 6-DOF dynamically balanced mechanisms [9] . Here, this linkage is studied to show that the current method replicates the results from literature.
The linkage consist of four revolute joints at r 1 up to r 4 . The first and last joints are attached to the base ( Fig. 4 ). The moving links have lengths l 1 , l 2 , and l 3 and the base length is l 4 . We analyse this linkage by cutting the third link between joint 3 and 4. Two kinematic branches I and II are obtained, consisting of joint 1, 2, and 3, and joint 4, respectively. Therefore, the open-chain equivalent design space matrix becomes Since all joints are planar revolute joints, the inertia transfer matrices consist of a point mass at the joint location ( Eq. 58 ), i.e. N i = N 0 = z m ( r i ) . This shows that the inertia parameter of link 3, at which the cut has been made, will lie on the span of N 3 and N 4 . The loop opening can also be performed at another link resulting in another basis for the same null space. Therefore, rank ( N O ) = 4 .
Through application of the higher-order derivatives algorithm, the X -matrix is obtained. Inspection of the rank X -matrix reveals that in the general case, this matrix has a maximum rank of rank ( X k max ) = 8 and requires k max = 3 derivatives to reach this limit. Since there are 12 inertia parameters, this linkage has a design freedom rank ( N ) = 4 . This shows that with the loop-opening method the complete design space of the general four-bar linkage has been found. This open-chain equivalent design space has no feasible solutions, since the addition of a positive point mass to one link requires a negative point mass on the connecting link. A feasible link with a positive MOI requires at least two positive point masses, which is not supported by this general solution.
Fortunately, there exist two special kinematic conditions in which the dimension of the design space increases and the linkage permits feasible dynamically balanced solutions [8,14] . These are the anti-parallelogram ( l 1 = l 3 , l 2 = l 4 ) and the kite-type four-bar linkage ( l 1 = l 2 , l 3 = l 4 ), as shown in Fig. 4 . In both cases, the base link has to be wider than two equal moving links. Here we study only the anti-parallelogram case, for which the rank of the X matrix drops to rank ( X ) = 7 and the design space is extended with a basis This may be interpreted as three pure-inertia rods in which the inertia connected to the coupler is opposite to the other two [39] . This solution arises since the links of the linkage follow a linear relation on the angular velocities ω 1 − ω 2 + ω 3 = 0 [40] . The solution is parameterized by y 1 up to y 5 . The first four parameters deal with the open-chain equivalent N O , while the variable y 5 is associated to N C . The resulting body masses and local MOIs become From the mass positivity condition ( Eq. 19 ) it follows that the design parameters should be selected decreasingly through the chain The positivity of the MOI places upper and lower bounds on the choice of y 5 (66) Fig. 5. The kinematic definition of the 3-RSR . During motion the end-effector stays tangential to an expanding sphere touching the platform centre and the base centre [41] . This mechanism is mirror symmetric with respect to a plane s passing through the spherical joints r 2 , r 5 , and r 8 ..
A similar result emerges for the kite-type four-bar linkage, with the sole difference that the crank and the coupler have an equal positive moment of inertia whereas the rocker has a negative moment of inertia. With this example we have shown that dynamic balancing results from literature could be replicated using the proposed higher-order dynamic balancing method. Furthermore we have shown an example in which the loop opening solution aids in interpreting the whole dynamic balance solution and the feasibility conditions. The special kinematic cases, as found by Moore and Gosselin [14] , may also be derived from rank reduction of X -matrix, but this lies outside the scope of this paper.

3-RSR spatial parallel mechanism
To show that this procedure also holds for multi-DOF linkages, we investigate a symmetric 3-RSR mechanism [41] . In this section we will show that the design space of the 3-RSR does not contain feasible designs when considering the whole of its workspace. Fortunately, it permits a dynamically balanced solution for 2-DOF movements over three planes of symmetry.

Kinematics of the 3-RSR mechanism
The 3-RSR under investigation ( Fig. 5 ) consist of 3 arms whose upper and lower links are of equal length l u = l l . The platform and the base have identical dimensions l b = l p . The first arm consist of two revolute joints 1 and 3, respectively hinged at the base at r 1 and at the platform at r 3 . Their joint axes n 1 and n 3 are co-linear in the reference configuration. The spherical joint, which connect link 1 and 2, is located at r 2 . Arms 2 and 3 are similar yet rotated with φ and −φ around the vertical axis. The platform is link 7. The three base joints are chosen to be the independent coordinates u 1 , u 2 , and u 3 . Throughout motion the mechanism is mirror symmetric with respect to a plane s passing through the three S -joints [41,42] . Due to this symmetry, the platform of this robot remains tangential ( r p ) to a variable sized sphere whose south pole is fixed at the origin ( r b ) of the robot.

General dynamic balance solution
Similar to previous example, the open-chain equivalent design space is found by cutting the three kinematic loops at the platform leaving three RSR chains, respectively consisting of joints 1-3, 4-6 and 7-9. The corresponding design space matrix is . 6. The 3-RSR mechanism can be dynamically balanced for motion over three additional planes of symmetry. For movement over the first plane, indicated by 1 , one may add a pure-inertia rods to the first arm and a negative pure-inertia rod to the platform η 7,1 . Also on arms 2 and 3 extra pureinertia rods may be placed in combination with a negative pure-inertia rod on the platform η 7,2 . Although not shown here, arm 3 has to be identical to arm 2. Later the balance for the other planes 2 and 3 is solved.
The base revolute joints have fixed axes while the platform joints have moving axes, giving rise to their respective inertia transfer matrices. Since a spherical joint ( S ) is instantaneously equivalent to three perpendicular intersecting revolute joints, its inertia transfer matrix N i is a point mass at the joint location r i , for i = 2, 5 and 8. This yields a open-chain equivalent design space with a rank of 3 × (5 + 1 + 3) = 27 . The X -matrix, which is found through the presented algorithm, terminates at k max = 2 and has a rank ( X k max ) = 43 . With 7 moving links and 70 inertia parameters to determine, it is shown that the complete design space is given by the open-chain equivalent design space ( 70 = 43 + 27 ) and no closed-chain design space exists. Therefore, additional counter-mechanisms are required for the dynamic balance of this 3-RSR mechanism.

Dynamic balance for 2-DOF motion on planes of symmetry
Fortunately, the mechanism can be dynamically balanced when only a part of the workspace is considered. Here we will investigate the case when the platform centre r p moves on another plane 1 that normal to joint axis n 1 ( Fig. 6 ). This occurs when the two other base joints move in unison u 2 = u 3 . During this 2-DOF motion the mechanism maintains kinematic mirror symmetry with respect to this plane. Later on, we will show that this permits a symmetric mass distribution over the 3 arms and the platform such that the mechanism is also dynamically balanced when moving over the other two planes of symmetry, i.e 2 when u 1 = u 3 and 3 when u 1 = u 2 .
When the mechanism's end-effector is constrained to move on plane 1 , the rank of the regression matrix drops by 8 to rank ( X k max ) = 35 , yielding a design space of dimension 35. A part of this enlarged design space can be explained by the fact that the S -joint of arm 1 now behaves as a revolute joint ( n 1 n 2 n 3 ) enlarging the open-chain design space by 6. Now the first arm's inertia transfer matrices (Eq. ) become Furthermore, an additional closed-chain design space of rank 2 appears, similar to that of the kite-type planar four-bar linkage ( Section 8.2 ). When the first joint ( u 1 ) is actuated and the others ( ˙ u 2 = ˙ u 3 = 0 ) are fixed, the first arm and the platform moves as a kite-type four-bar linkage due to the mirror symmetry with respect to the plane s through the Sjoints ( Fig. 5 ). Due to this symmetry for all poses, the angular velocities of the bodies have a particular linear relation ω 1 + ω 2 = ω 7 . The first column of the closed-chain design space matrix is therefore For the second 'symmetric' DOF u 2 = u 3 a second mirror symmetry is of importance; the mirror symmetry with respect to the plane 1 perpendicular to n 1 and passing through the centre of the base ( Fig. 6 ). This is encoded in an angular velocity relation ω 1 + ω 2 + ω 7 = ω 3 + ω 4 + ω 5 + ω 6 . Since arms 2 and 3 move in mirror symmetry with this plane, the second column is found to be So arms 2 and 3 move in opposite directions with respect to this additional symmetric plane 1 ( Fig. 6 ), cancelling the angular momentum components in this plane. In the perpendicular direction ( n 1 ) the angular momentum is counter-acted by the negative pure-inertia rod of the platform.
The principal MOIs follow from the summation of the pure-inertia rods. Here the first principal axis is n 1 , the second is a 1 , and the third is a 1 × n 1 . The principal MOIs of the first link are For the second link the principal axes are oriented similarly; in the directions of n 3 , a 2 , and n 3 × a 2 , respectively The MOIs of the links of the other 2 arms are equal g 1 = g 3 = g 5 and g 2 = g 4 = g 5 . The MOIs of the platform are formed by the closed-chain design space basis y 34 , the added point masses at the shoulder y 11 , the inertia transfer y 13 over the connecting joints (i.e. in the direction of n i for i = 3, 6, and 9), and the modification of by inertia transfer over the first chain y 14 . The first principal MOI is around a 7 , the second around n 3 and in the third out-of-plane direction n 3 × a 7 g 7 , 1 = 3 2 l 2 p y 11 + y 13 − 1 2 y 8 − y 14 , g 7 , 2 = 3 2 l 2 p y 11 + y 13 − y 8 − 9 4 y 34 , g 7 , 3 = g 7 , 1 + g 7 , 2 + 2 y 14 .
Now the platform is chosen axisymmetric g 7 , 1 = g 7 , 2 by setting This axisymmetric solution extends the dynamically balanced motion over the two other planes of mirror symmetry 2 and 3 , which are perpendicular to the base joints n 4 and n 7 , respectively. It should be noted that the resulting mechanism is force balanced for its complete workspace.
By invoking the feasibility conditions ( Eq. 19 ) we obtain the following set of feasibility conditions for the 3-RSR mechanism moving in a single plane 1 . It should be noted that several shadowing inequality conditions have been removed 0 < y 11 < y 6 < y 1 , y 13 < y 8 < y 3 , 2 l 2 u y 6 y 11 m 2 < y 34 , 0 < y 34 + y 5 + l 2 u y 1 1 − It should be noted that similar feasibility conditions, yet more complex, can be obtained for the full 35-dimensional design. Furthermore, the satisfaction of the feasibility conditions does not guarantee a practical implementation of the resulting balance solution. The solution, for example, may require very slender links or centres of mass at a long distance from the joint, which may be hard to realize in practice.

Numerical study
Based on these findings a geometry and design parameters y are selected for simulation ( Table 2 , and 3 ). The corresponding mass distribution is found in Table 4 . In here the COM location of the arms are defined as c i = r i − c i a i . Whereas c i denote the distance to the connecting joint ( Fig. 7 ). The COM of platform is at the centre of the platform c 7 = r p .
Simulations of this mechanism in a multibody software package confirm the dynamically balanced mass distribution by showing shaking forces and moments that are effectively zero for movement over the three symmetric planes ( Fig. 8 ). The remaining shaking force and moments are attributed to round off errors. Table 3 The design parameters of the dynamically balanced 3-RSR mechanism as used in the numerical example.

Discussion
The proposed method achieves dynamic balance by eliminating a finite but sufficient number of higher-order partial derivatives of the momentum wrench. This results in the necessary and sufficient balancing conditions since the momentum is an analytic function outside of singularities. Therefore, if the higher-order derivatives of such an analytic function are zero, the function itself is identically equal to zero and no shaking forces and moments are exerted on the base.
In this paper it is shown ( Section 5.3 ) that for open-chain linkages up to two derivatives ( k max = 2 ) are required to obtain the necessary and sufficient conditions for dynamic balance. For closed-chain linkages such a general upper-bound on the number of derivatives is not found. Therefore, each case requires an inspection of the X -matrix, as done in the examples. Although a constant rank of the X -matrix for higher orders it is a strong indicator, verification of the balancing solution is still formally required, for example by numerical simulation. Furthermore, different linkages, but also different geometries lead to different k max values. For example, for the general planar four-bar k max = 3 while for the (anti-)parallelogram geometry k max = 2 .
By selection of a sub-matrix of the X -matrix, a larger design space may appear and the conditions for exact dynamic balance are relaxed in favour of other requirements. For example, force balance only may be obtained by selecting the relevant rows. Alternatively, approximate balancing may be obtained by satisfying only derivatives up to a certain order, yielding dynamic balance in a neighbourhood of the reference pose. Path balance may be obtained by only regarding paths or planes through the reference pose, as shown by the 3-RSR example. This method therefore generalizes the result of [26] .
The feasibility conditions ( Eq. 19 ) dictate whether or not a desired mass distribution can be constructed. To check the feasibility, an eigenvalue problem is to be solved. The non-linear nature of the eigenvalue problem however, prevents a closed-form treatment for arbitrary linkages. Currently this study is done on an individual basis. Yet, for simple linkages the partitioning and multipole-rod representation allowed for an interpretation and feasibility study without eigenvalues decomposition. Furthermore, even though some of the linkages studied here have no feasible design space, the complete characterization of the design space is expected to assist in the selection of the most favourable application dependent solution, i.e. with minimal motor torque and/or minimal number of counter-mechanisms.

Conclusion
This paper presented a higher-order dynamic balancing method that yields and solves the necessary and sufficient dynamic balance conditions of open-or closed-chain, planar or spatial linkages in nonsingular configurations. This method relies on a recursive algorithm to obtain the higher-order derivatives of the rigid body momentum equations up to an arbitrary order. In the parameter-linear form these higher-order derivatives furnish the dynamic balance conditions that arise from a rank deficiency of the regression matrix. This enables generic null space procedures to yield the complete set of dynamically balanced inertia parameters, i.e. the complete dynamically balanced design space of a given linkage. It was shown that for open-chain linkages, derivatives up to the second order yield the necessary and sufficient balance conditions. For such linkages the partial derivatives of higher orders will pose no new dynamic balance conditions. For closed-chain linkages such an upper bound was found per case. Partitioning and interpretation of the design space of open-chain linkages was presented based on the concept of inertia transfer. Six inertia transfer matrices where found that, depending on the alignment of the joints and the types of joints in the chain, allow for an exchange of inertia parameters without affecting the momentum generated by the linkage as a whole. For general closed-loop linkages this open-chain equivalent design space fully describe the dynamically balanced design space, resulting in non-feasible designs due to negative masses or moments of inertia. Fortunately, in special cases a the loop-closure permits a larger design space enabling feasible and constructable balance solutions. Also, a larger design space is found when only a part of the workspace is to be dynamically balanced. For the 3-RSR mechanism, for instance, three dynamically balanced planes of symmetric 2-DOF motion were obtained without resorting to counter-mechanisms.
With this method a systematic, closed-form dynamic balance algorithm is presented that provides a uniform framework to study linkages that require no or a minimal number of counter-mechanisms, paving the way for the dynamic balancing of high-speed robots.

Declaration of Competing Interest
None.

Acknowledgements
The second author acknowledges that this work has been partially supported by the Austrian COMET-K2 program of the Linz Center of Mechatronics (LCM).

Appendix. Multivariate matrix derivatives using Kronecker product
The higher-order partial derivatives of matrices are managed in this paper with the use of the Kronecker product [33] . In this appendix the Kronecker product notation and some elementary properties are listed along with the application to the bookkeeping of the higher-order partial derivatives. The differences between the repeated higher-order partial derivatives D α x ( A ) and the collection of higher-order partial derivatives D k x ( A ) are highlighted.
The Kronecker product is defined as the collection of the element-wise multiplication of all elements in the respective matrices. Consider the following set of matrices: A ∈ R n ×m , B ∈ R p×q A B = ⎡ ⎣ a 11 B · · · a 1 m B . . . . . . . . . a n 1 B · · · a nm B ⎤ ⎦ . (A.1) The vectorization of a triple matrix product is written by means of the Kronecker product vec ( ABC ) = ( C A ) vec ( B ) .

(A.2)
A set of nested higher-order partial derivatives is denoted with a vector in the superscript D α Consider for example the following repeated partial derivatives 3) The α-vector is an ordered multi-index corresponding to x . The total order is given by i α i = | α| . These repeated higherorder partial derivatives form the elements (rows and submatrices) of the collection of higher-order partial derivatives, i.e. D k x ( A ) is formed by all nested partial derivatives with respect to α for which | α| = k .
The collection of partial derivatives of a given vector a with respect to x ∈ R r are organized according to

(A.4)
This is extended to the partial derivatives of a matrix A = a 1 . . . a m , which are organized according (A.5) The following entities may be found for the derivatives of matrices A ( x ) ∈ R n ×m , B ( x ) ∈ R m ×q and x ∈ R r Product rule For a matrix product the following derivative holds Kronecker product The derivative of the Kronecker product is written using the permutation matrix P q,r [33] D A recursive application of these equations allow for the extension of these derivatives up to arbitrary order.