Higher-order constraints for higher kinematic pairs and their application to mobility and shakiness analysis of mechanisms

The mobility analysis of mechanisms rests on an adequate formulation of the constraints defining its configuration space (c-space). Whereas there is no general method for a global analysis, the higher-order mobility analysis, which locally approximates the c-space, is applicable to general mechanisms. It requires an efficient method for the evaluation of higher-order constraints, i.e. constraints on velocity, acceleration, jerk, etc. Such a method is known for linkages comprising lower pair joints only. In this paper a method for the efficient evaluation of higher-order constraints for mechanisms comprising higher pair joints is proposed. The method builds on the results for the lower pair linkages. It leads to a computationally simply recursive algorithm. This is applied to the mobility analysis that allows to determine the local finite mobility, to detect singularities, and to identify shaky mechanisms.


List of symbols n
Number of 1-DOF joints in a linkage q 2 V n Vector of joint variables Constraint mapping for a higher pair cutjoint in a kinematic loop h l Constraint mapping for a higher pair cutjoint in the fundamental cycle K l f i 0 Forward kinematic mapping of the open chain with terminal body i 0 obtained by removing the cut-joint from loop f a Forward kinematic mapping of the open chain with terminal body B a obtained by removing the cut-joint from the FC K l of which the higher pair joint J l connects to B a Y i Screw coordinate vector of joint i in the zero reference configuration q ¼ 0 S i Instantaneous screw coordinate vector of joint i V i Spatial twist vector of body i C Topological graph G, H Spanning tree, and cotree in C c Number fundamental loops of C K l Fundamental Loop (FL) of C assigned to cotree edge l 2 H

Introduction
Lower kinematic pairs (Reuleaux pairs) are mechanical generators of SE ð3Þ subgroups, and the instantaneous joint screws form a basis for the corresponding se ð3Þ subalgebras. The finite joint motion is found from the instantaneous motion via the exponential mapping on the Lie group SE ð3Þ. When studying the kinematics of linkages with lower kinematic pairs, it is thus convenient and common practice to express the kinematics in terms of joint screw coordinates and joint parameters. This is known as the product of exponentials (POE) formulation [1] for the kinematic mapping of a kinematic chain, where relative motions due to lower pairs are expressed as exponential of joint screw coordinates. The same applies to the closure constraints of kinematic loops. This gives rise to explicit algebraic relations of the derivatives of the kinematic mapping, and hence of the derivatives of the loop closure constraints for lower pair linkages. These higher-order constraints were employed to the mobility analysis of lower pair linkages [5,6,11,18,20,22].
In [17] a generally applicable computational method for mobility analysis using higher-order approximations was presented. The Lie group setting further allows for mobility estimation without resorting to higher-order approximation but rather considering intersection of motion spaces [7,8,12,21]. Higher kinematic pairs do not represent mechanical generators of motion subgroups of SE ð3Þ, and can in general not be represented by combination of such. Hence a screw-based formulation of the kinematics and loop constraints is not available. However, frequently mechanisms comprise lower as well as higher pairs. In particular, the majority of mechanisms possess kinematic loops comprising only a few higher pairs. Such mechanisms were for instance addressed in [30]. For the special case where each loop contains no more than one higher pair, a recursive formulation is presented in this paper making use of the POE formulation. To this end, the higher pair joint is removed from the loop. The motion of the two bodies connected by the higher pair is expressed by a POE. The constraints are then imposed on the relative motion of the bodies. Since the higher pair is cut open and replaced by the corresponding constraints, this is referred to as the 'cut-joint' of the loop, and this approach is known as cut-joint formulation of constraints [15,28]. This is also applicable to lower pair mechanisms. Based on the constraint formulation a computational method for the higher-order mobility analysis is introduced in the following.
The paper consists of two parts. In the first part the higher-order loop constraints are presented that are employed for mobility and singularity analysis in the second part. In Sect. 2 the geometric constraints for a single kinematic loop are introduced. The corresponding higher-order constraints are derived in Sect. 3. The extension to multi-loop mechanisms is presented in Sect. 4. This rests on the topological graph and the identification of topologically independent fundamental cycles. The application of these constraints to the mobility and singularity analysis is presented in Sect. 5. Throughout the paper, u 1 ¼ 1; 0; 0 ð Þ T , u 2 ¼ 0; 1; 0 ð Þ T and u 3 ¼ 0; 0; 1 ð Þ T denote the standard unit basis vectors.

Cut-joint versus cut-body formulation
The relative configuration of adjacent bodies is determined by the configuration of the connecting joint. The latter is determined by the joint variables (angles/dis-placements) that are referred to as relative coordinates. In terms of relative coordinates, the closure constraints for a kinematic loop can be formulated in two different ways: the cut-body formulation [9,23] and the cut-joint formulation [28]. The cut-body approach is commonly used to model the kinematics of lower pair linkages, whereas the cut-joint approach is the standard approach to derive dynamic motion equations of multibody systems (MBS) in relative coordinates. In the cut-body formulation the relative configurations of all adjacent bodies are successively combined. That is, the constraints involve all joint variables of the loop. This approach can be interpreted as cutting one body in two parts, and the reassembly of the body yields the closure condition. In the cut-joint approach, one joint (called the cut-joint) is removed and the remaining two open kinematic chains are subjected to closure constraints according to the cut-joint. The kinematic model, for a loop with n joints and bodies, thus comprises only the joint variables of the remaining n À 1 joints.
It is to be mentioned that elimination of the cut-joint may lead to artifacts. For instance consider the planar 4-bar mechanisms in Fig. 1a where joint 4 is used as cutjoint. In the shown configuration the model is in a singular configuration since the axes of joint 1, 2, and 3 are coplanar, and thus the constraint Jacobian drops rank. If joint 2 would be used as cut-joint (Fig. 1b), the Jacobian would have full rank 2. Apparently this singularity is due to the model rather than being a kinematic phenomenon of the mechanisms. The reason is that, by eliminating the cut-joint, only a n À 1 dimensional subspace of the configuration space is considered. This is, however, often ignored in MBS modeling.
Also from a computational point of view, the cutbody formulation is advantageous since the closure condition simply requires that the configuration of the virtual cut-body is the identity. In the cut-joint approach, the constraints are specific to the cut-joint. Since this is the standard approach in MBS dynamics using relative coordinates, details on the constraints for various joints can be found in the relevant literature [26,28].
The main advantage of using relative coordinates to express the configuration of lower pair joints is that they can be expressed in terms of joint screws via the exponential mapping on SE ð3Þ. This has been employed for mobility analysis in [5,6,17,20,22]. The crucial point is that the derivatives of the constraints are available explicitly in terms of joint screws. This is not possible for higher pair joints. However, if a loop comprises only one higher pair joint, this can be used as cut-joint, which yields a lower pair linkage with two open kinematic chains whose kinematics can be modelled in terms of the remaining joint screws. This is presented in the following.

Relative configurations of adjacent bodies
Consider a kinematic loop, connected to the ground (indexed with 0), comprising n lower pair joints and one higher pair joint. W.l.o.g. the lower pairs are assumed to have 1 DOF, i.e. could be revolute, prismatic, or screw joint. This is valid since any lower pair (and generally, any mechanical generator of motion subgroups) can be represented as combination of these. Then, denote with q 2 V n the vector of joint variables of the lower pairs. The higher pair is removed from the loop and serves as cut-joint of the loop, as indicated in Fig. 2. The two bodies connected by the higher pair represent the terminal bodies of the remaining two open kinematic chains, and are indexed with k 0 and l 00 , respectively. The 1-DOF lower pair joints of the respective open kinematic chain from the ground towards bodies k 0 and l 00 are indexed with 1 0 ; . . .; k 0 and 1 00 ; . . .; l 00 .
The configuration of a body is represented by the configuration of a body-fixed reference frame (RFR) relative to a spatial reference frame F 0 . This configuration is represented by a 4 Â 4 homogenous  transformation matrix C 2 SE ð3Þ (Appendix). The configurations of the bodies k 0 and l 00 , i.e. of the frame F k 0 and F l 00 , are determined by C k 0 ¼ f k 0 q ð Þ and C l 00 ¼ f l 00 q ð Þ, with the POE formulae (Appendix) Here A k 0 and A l 00 is the respective reference configuration of body k 0 and l 00 at the zero reference configuration q ¼ 0 of the mechanism, and Y i ¼ e i ; s i Â ð e i þ h i e i Þ T is the screw coordinate vector of the lower pair joint i w.r.t. to a global reference frame F 0 at q ¼ 0. Therein, e i is a unit vector along the axis of joint i, and s i is the position vector to any point on that axis, both expressed in the global reference frame [24].
The relative motion of body l 00 w.r.t. to body k 0 is thus given by C À1 k 0 C l 00 . The cut-joint constrains this relative motion. For higher pair joints this is not a screw motion, in contrast to lower pairs where it can be expressed by the exponential of a cut-joint screw.

Translation constraints for higher pairs
The translation constraints restrict the relative displacement of a point on body k 0 and a point on body l 00 . In order to express this condition, the point coordinates must be expressed in a common reference frame. The location of a point is represented by homogenous point coordinates.
Denote with p k 0 and p l 00 the position vector of a point fixed on the respective body expressed in the body-fixed reference frame F k 0 and F l 00 (Fig. 2), and with p k 0 ¼ ðp T k 0 ; 1Þ T and p l 00 ¼ ðp T l 00 ; 1Þ T their homogenous point coordinates. The homogenous coordinates expressed in the global frame F 0 are C k 0 p k 0 and C l 00 p l 00 .

Remark 1
The mapping e p : C 2 SE ð3Þ7 !Cp 2 E 3 is the evaluation mapping of SE ð3Þ that describes how the transformation group SE ð3Þ acts on the Euclidean space, i.e. how point coordinates transform under rigid body motion.
The relative displacement vector of the two points, expressed in the global frame, is thus with (1) Since the joint is fixed to the bodies, the components of the relative translation to be constrained can be identified when expressed in a body-fixed RFR. In the RFR F k 0 at body k 0 , for instance, the displacement vector is (the difference of two point coordinates is a vector). The system of m 3 (scleronomic) translation constraints can now be written in the from This accounts for general bilateral contact conditions like surface-surface, point-surface, or curve-point contact. A simple constraint for technical joints is, for instance, that the relative motion along a certain direction (which is fixed w.r.t. the body) is prohibited. Let e k 0 be a unit vector along this direction, expressed in the RFR F k 0 , then the corresponding constraint is

Orientation constraints for higher pairs
The orientation of a body w.r.t. the global frame is represented by the rotation part R 2 SO 3 ð Þ of its configuration C ¼ R; r ð Þ 2 SE ð3Þ. The relative rotation constraint of the two bodies k 0 and l 00 can be expressed as where e k 0 and e l 00 is a unit vector fixed at body k 0 and l 00 , respectively. Here R T k 0 R l 00 is the relative rotation of body l 00 w.r.t. body k 0 . The constraint (5) enforces the two unit vectors to be perpendicular. In general, the unit vectors e can depend on the configuration C. In total m 3 such rotation constraints can be introduced. An exhaustive list of constraints for various joints can be found e.g. in [25,26]. Clearly the above constraint formulations are also applicable to lower pair joints.
Example 1 (Hook Joint) A hook (or universal) joint is a 2-DOF higher kinematic pair representing a Cardanic suspension (Fig. 3). It cannot be a lower kinematic pair since there is no 2-dim motion subgroup. It could actually be modeled as combination of two (lower pair) revolute joints, but it is here considered as higher pair. The joint imposes one orientation constraint (5): the two rotation axis must remain perpendicular. By convention, the body-fixed frames are usually oriented so that e k 0 is along the 1-axis of the body-fixed RFR F k 0 on body k 0 , and e l 00 is along the 2-axis of the RFR F l 00 on body l 00 [25,26]. Then e k 0 ¼ u 1 and e l 00 ¼ u 2 , and the constraint requires the 1, 2-entry of the relative rotation matrix R T k 0 R l 00 to vanish: u T 1 R T k 0 R l 00 u 2 ¼ 0. The hook joint constrains the origin of the body-fixed frames to coincide. That is, p k 0 and p l 00 are the location vectors of the intersection point of the joint axes. The m ¼ 3 translation constraints (4) are thus simply d k0 q ð Þ ¼ 0. In total the hook joint imposes 4 constraints.
Example 2 (Pin-in-Slot Joint) A pin-in-slot joint is a 2-DOF higher kinematic pair restricting the relative motion of two bodies to a line while allowing for a relative rotation about a constant axis. Figure 4 shows a planar single-loop mechanism containing a pin-inslot and two revolute joints. The pin-in-slot joint is used as cut-joint. The remaining n ¼ 2 joint angles defined the parameter space V n ¼ T 2 . The pin-in-slot joint connects the coupler link l 00 ¼ 2 of the 4-bar mechanism with the ground body k 0 ¼ 0. The global reference frame F 0 is located at ground as shown in Fig. 4. The configuration of the ground is of course C 0 I. The configuration C 2 of link 2 is determined by f 2 q ð Þ ¼ expðY 1 q 1 Þ expðY 2 q 2 ÞA 2 with the joint screw coordinates and the reference configuration A 2 ¼ ðI;ðÀb;d;0Þ T Þ w.r.t. the world frame. The pivot point of the pin-inslot joint is the coupler point of link 2. Its position expressed in the RFR F 2 on link 2, shown in Fig. 4, is This point is restricted to move on the horizontal line fixed to the ground, and any point on the joint's sliding axis can be used as reference point on the ground. The center of the slot is used as reference point, and its coordinate vector expressed in the ground RFR (world frame) is is the relative displacement vector expressed in the RFR F 0 on link k 0 ¼ 0. The motion along the 1-axis of the RFR F 0 of link k 0 ¼ 0 is the unconstrained joint motion. The two translation constraints require the vanishing of the second and third component of d k 0 . These are with unit vectors u 2 ¼ 0; 1; 0 ð Þ T and u 3 ¼ 0; 0; 1 ð Þ T . The constraint (9) is apparently dispensable. The rotation constraints require the 3-axes of both bodies, k 0 ¼ 0 and l 00 ¼ 2 to be parallel. According to (5) the two constraints are where R 2 is the rotation matrix of body l 00 ¼ 2. Since the example is a planar mechanism these constraints are also dispensable, and their detailed expression is omitted here. Hence, the system of geometric constraints (4) consists of the single Eq. (8).
Remark 2 It should be noticed that the solution of the m 3 orientation constraints constructed with (5) is not unique. This is because the orthogonality relations do not enforce a specific direction of the unit vectors. These constraints are sufficient though as long as an admissible configuration is used for a local mobility analysis, or as initial configuration for the time integration in MBS dynamics. But when assembly configurations are determined numerically, it is common practice to impose m þ 1 4 constraints. For instance, the additional constraint for the hook joint is u T 3 R 2 u 3 ¼ 1.

Higher-order constraints for a kinematic loop
Repeated time derivatives of the constraints yield the corresponding constraints on the acceleration, jerk, jounce, etc. The above formulation of cut-joint constraints requires the time derivatives of the kinematic mappings (1) of the two lower pair chains. They become very complex. In order to handle this complexity a recursive formulation was presented in [16]. In this section they are employed to derive recursive formulations for the cut-joint constraints of higher pairs.
3.1 Time derivatives of the twist of the terminal bodies k 0 and l 00 The time derivative of the geometric constraints yields the velocity constraints for higher pairs. The orientation constraints (5) only involve the relative rotation matrix, whereas the translation constraints (4) involve the translation and orientation of the two interconnected bodies. This reflects the semidirect product structure of SE ð3Þ, meaning that derivatives of (5) involve only angular velocities, and those of (4) involve the complete relative twist. The recursive relations for the derivative of the relative twists are briefly summarized in the following.
The spatial twist of a body is defined as b V ¼ _ CC À1 2 se ð3Þ (Appendix). With the forward kinematic mapping (1), the twist of body k 0 at configuration q is given in terms of the joint rates _ q by [13] and q Therein, the vector of instantaneous screw coordinates of joint i 0 is determined by where Ad g i 0 is the frame transformation of the screw coordinate vectors from the zero reference configuration q ¼ 0 to the current configuration, with It was shown in [16] that the derivatives of S m k 0 admit the following recursive expression together with the derivatives of the screw coordinates where D ðiÞ :¼ d i dt i , and Á; Á ½ is the screw product (Lie bracket on se ð3Þ), see [24] and Appendix. The higherorder time derivatives of S 1 k 0 , and thus of the twist (11), are hence determined recursively by very simple algebraic operations. This applies analogously the other terminal body l 00 .

Translation constraints for higher pairs
The time derivative of the relative displacement vector d, i.e. the relative translational velocity, expressed in the spatial frame F 0 is where the position vectors p k 0 and p l 00 , expressed in the respective RFR F k 0 and F l 00 , are assumed to be constant. This relative velocity expressed, for instance, in the RFR on body k 0 is R T k 0 _ d, with R k 0 the rotation matrix in C k 0 (48).
Starting from the first expression in (17), higher time derivatives of d are given as They involve the derivatives of the configurations. These follow immediately from _ and analogously for C l 00 . The vector in (18) is expressed in the spatial frame. In order to formulate the higher-order constraints, they are transformed to the RFR on body k 0 as R T k 0 D i ð Þ d. The ith-order constraints are then of the form In summary, the relation (18), together with (19), (15) and (16), provide a recursive algebraic expression for the ith time derivative of the relative displacement vector d k 0 . This allows for a recursive determination of time derivatives of the constraints (4). The basic operation is the Lie bracket in (16).

Orientation constraints for higher pairs
Time derivatives of the orientation constraints (5) yield the constraints on the relative angular velocities of bodies k 0 and l 00 . This require time derivatives of the relative rotation matrix. The ith time derivative is given by derivatives of the rotation matrices of body k 0 and l 00 The relations _ R Here, x k 0 and x l 00 are the angular velocity of the spatial twists V k 0 and V l 00 . Hence, their derivatives are known from the angular part of the derivative of S 1 k 0 and S 1 l 00 , using (15). That is, the latter are determined already if the cut-joint imposes translation constraints, so that the relations (22) and (23) do not have to be evaluated.

Linkage topology
Most mechanisms comprise multiple kinematic loops. The above formulation of higher pair cut-joint constraints for a single loop can be directly adopted to multi-loop mechanisms. To this end, topologically independent loops must be identified. The kinematic mechanism topology refers to the arrangement of links and joints. This is naturally described by a linear nondirect graph referred to as the topological graph, denoted C B; J ð Þ [4,15,28]. The vertex set B represents the bodies, and the set of edges J represents the joints. In the following the basic description is recalled from [15], which provides a systematic approach. This serves to identify the kinematic loops for which the above constraints must be formulated.
Denote with G & C a spanning tree on C, and with H ¼ CnG its complement-the cotree. A fundamental cycle (FC) is a closed path containing exactly one edge of the cotree H. The number of FCs is the cyclomatic number c :¼ jJj À jBj þ 1, with |J| being the of number joints, and |B| the number of bodies. In particular, jJj ¼ n when assuming that all joints have 1 DOF. The c FCs constitute a system of topologically independent loops for which closure constraints are to be introduced.
Bodies are indexed with a; b; . . . and denoted with B a 2 B. The ground is denoted B 0 . Joints are indexed with i; j; . . . and denoted J i 2 J. Each FC contains exactly one cotree edge. The FC containing the cotree edge J l 2 H is denoted with K l . Figure 5 shows a twoloop mechanism, and its topological graph C.
In order to apply the formulation introduced in the previous section, no FC must comprise more than one higher pair. Consequently, the FCs on C must be introduced accordingly. If this is not possible, the method is not applicable. In other words, the higher pair cut-joint necessarily corresponds to the cotree edge J l of K l . This is J 3 in the example in Fig. 5. Furthermore, the two bodies connected by J l are the terminal bodies of two kinematic chains. Now, adoption of (1) requires a predecessor relation.
In G there is a unique path from any vertex (body) B a to B 0 (ground). For a given G, a root-directed tree is introduced, denoted with G 0 , such that all edges are oriented away from the ground (see Fig. 5e). Edges of G 0 are directed and represented by ordered pairs of vertices, i.e. bodies.
Let B a and B b be connected by a tree-joint. Then is the source and B a the target of the edge. This is denoted as there is a finite k, such that j ¼ i À 1 À 1. . . À 1 (k times). This is expressed as j ¼ i À k. Being a predecessor is indicated by j\i.
The tree-joint that connects B a with its predecessor is denoted with J a ð Þ. The last joint of the chain from B a within the tree that connects to the ground B 0 is denoted with J root a ð Þ. An edge of C is a tuple J i ¼ B a ; B b ð Þindicating that joint i connects bodies a and b. To indicate orientation of the joints, i.e. the meaning of relative joint motions, an oriented graph C is introduced. If J i ¼ B b ; B a À Á 2 C, then joint J i determines the motion of B a w.r.t. B b , and B a is the target and B b the source of the arc. This is represented by an arrow as in Fig. 5d.
Finally, the joint orientation must be related to the orientation of G 0 . To this end, introduce an indicator function for

Geometric constraints
As in (1), the configuration of body a is determined by the relative configurations of the preceding lower pairs in the chain towards the ground as C a ¼ f a q ð Þ, with the forward kinematic mapping e Oriented spanning tree f a q ð Þ:¼ expðr r ð ÞY r q r Þ Á . . .
where A a 2 SE ð3Þ is the reference configuration of B a . The same applies to B b , so that the relative configuration of the two bodies connected by cut-joint 2 H is known. Translation constraints Denote with p l;a and p l;b the position vector of points at B a and B b , respectively expressed in F a and F b , whose relative displacement is constrained by the cut-joint J l . The relative displacement expressed in the global frame is The system of m 3 translation constraints for the cutjoint J l ¼ B b ; B a À Á then attains the form Orientation constraints An orientation constraint for J l imposes a restriction on the relative rotation. Denote with e l;a and e l;b , specified unit vectors, fixed at B a and B b and represented in F a and F b , respectively. The constraint is introduced as where R a and R b are the rotation matrices of the two bodies. Up to three independent constraints can be introduced for the cut-joint J l .
Example 3 (Pin-in-Slot Joint, cont.) The mechanism in Fig. 5 comprises two FCs, K 3 and K 5 with cut-joints J 3 and J 5 , respectively. The remaining n ¼ 3 joint angles q 1 ; q 2 ; q 4 define the parameter space V n ¼ T 3 . The joint orientations (indicated by C in Fig. 5d) are such that r i ð Þ ¼ 1 for all i. This is in accordance with Example 2. The higher pair J 3 ¼ B 0 ; B 2 ð Þconnects B 0 and B 2 , and the corresponding constraints are formulated in terms of C 0 and C 3 (Fig. 6). The translation and orientation constraints are thus the two constraints (8) and (10), determined with C 0 ¼ I and f 2 q ð Þ in Example 2. The lower pair revolute joint Þ connects B 0 and B 3 . For sake of simplicity, the RFR F 3 on link B 3 is located at the axis of J 5 so that A 3 ¼ ðI; ða; 0; 0Þ T Þ. The configuration of B 3 is determined by f 3 q ð Þ ¼ expðY 1 q 1 Þ Á expðY 2 q 2 Þ Á expðY 4 q 4 ÞA 3 with the joint screw coordinates in (6) and The position vector to the revolute joint axis, expressed in the RFR where c 124 ¼ cosðq 1 þ q 2 þ q 4 Þ, etc. The revolute joint constraints require all three components to vanish. The third component is identically zero, and thus dispensable, for this planar mechanism. The overall system of constraints (26) consists of the Eq. (8) and the first two equations in (29).

Remark 3
The FC K 5 only contains lower pairs. The constraints for fundamental cycles only comprising lower pairs can be formulated using the POE formulation reported in [16] and the mobility analysis be pursued correspondingly [17].

Higher-order constraints
Relative twists The spatial twist of B a is again a linear combination of the instantaneous joints screws in terms of joint rates where now the tree-joints in the path from B a to B 0 are involved, so that The instantaneous joints screws are now with The relative spatial twist of B a and B b connected by cut-joint J l is thus known. Its derivatives are available adopting (15) and (16). This is straightforward, and omitted here.
Translation constraints The relations of Sect. 3.2 are directly applicable. The time derivative of the relative position vector is and in the RFR on Higher derivatives follow by adapting the results in Sect. 3.2.
Orientation constraints The relations for the orientation constraints in Sect. 3.3 can also be adopted directly and are omitted here.

Remark 4
Finally it must be remarked that the FCs are not unique. Using different systems of FCs results in different loop constraints. The FCs must be introduced such that each FC contains at most one higher pair.

Problem statement
The motivation for deriving the above recursive relations for the time derivatives is to provide a computationally efficient means for the local kinematic analysis of mechanisms comprising lower pairs and one higher pair per FC. This allows to investigate the mobility, identify singularities and motion bifurcation and reconfiguration, and to assess the shakiness of a mechanism. A local analysis aims at determining the higher-order, and eventually the finite, mobility of a mechanism in a given configuration.
Denote with h l q ð Þ ¼ 0 the system consisting of the translation constraints (26) and the orientation constraints (27) for the cut-joint J l , and with h q ð Þ ¼ 0 the overall system of all cut-joint constraints. The c-space of a multi-loop mechanism is the analytic variety A configuration q 2 V is a feasible assembly. The local geometry of V at a configuration q 2 V reveals the finite local mobility. A local analysis hence attempts to approximate the local geometry of V.
The first-order approximation of V is its tangent space T q V. This is indeed sufficient if V is a smooth manifold, i.e. q is a regular point of V. At any general (singular or regular) point q 2 V, the best local approximation is the tangent cone [27], denoted with C q V. Whereas T q V consists of tangent vectors to V, the tangent cone consists of tangent vectors to curves in V through q. It is hence the restriction of tangent vectors x 2 T q V to those that are tangent to finite motions represented by curves in the c-space. It thus provides a local model of V, and in particular its local dimension d loc q ð Þ : The concept of tangent cone has been first applied to the analysis of lower pair linkages in [11] that was continued in [18]. There the tangent cone was introduced via higher-order closure constraints, which allows for its actual determination thanks to lower pair loop constraints using the POE formulation. A summary of the computational procedure can be found in [17]. While the tangent cone provides a mathematical framework, higher-order constraints have been used for mobility analysis without resorting to this concept, e.g. [3,5]. In particular, the closed form expressions of derivatives in term of Lie brackets for linkage with lower pairs have been reported and applied in [2,6,22]. Mechanisms with higher pair joints have not yet been addressed. The above higher-order constraints now allow for a higher-order analysis of mechanisms with higher pairs while exploiting the recursive form of the derivatives of lower pair chains. It is common practice to assign the variables q ¼ 0 to the configuration to be analyzed. Then, C l 00 ¼ A l 00 and (2) and (17), and the instantaneous joint screw coordinates in (15), (16), and (19) are simply the joint screw coordinates in the reference configuration, i.e. S i 0 ¼ Y i 0 , respectively in (31).

Computational procedure
The basic idea is to identify subsets of first-order motions that satisfy the constraints up to order i. With increasing i, this yields the set of first-order motion that correspond to finite motions satisfying the geometric constraints. To this end, introduce the following mappings The constraints of order i are hence These are polynomials in _ q; . . .; q i ð Þ . Therewith the velocity constraints read H 1 ð Þ q; _ q ð Þ ¼ 0, and the vector space of feasible first-order motions is The set of vectors x 2 K 1 q satisfying the constraints up to order i is For certain order j, this sequence terminates with the tangent cone, Tangent vectors to finite curves belong to C q V and thus to all K i q . Each K i q is a cone in x, i.e. if x 2 K i q then also ax 2 K i q for a 2 R. The tangent cone is an algebraic variety of order j. That is, the c-space (an analytic variety) is locally replaced by a low-order algebraic variety.
The following remarks are in order: 1. The first-order cone K 1 q (null space of the constraint Jacobian) is not necessarily the tangent space to V at q. The tangent space is equal to the tangent cone if and only if, q 2 V is a regular point of V. The tangent space is known once the tangent cone is determined, since q , this cannot be deduced by comparing K 1 q and C q V. However, once the tangent cone is available, also the tangent space is known. Moreover, with (40), q 2 V is a c-space singularity iff C q V 6 ¼ spanC q V. This can actually be checked with the proposed computational method. For completeness it should be mentioned that (in contrast to the common belief) even in regular points q of V the dimension dim K 1 q may not be locally constant (i.e. q is a constraint singularity [19]). That is, the number of independent constraints may drop at q even if V is locally a smooth manifold at q (see the Goldberg 6R example in [17]). 2. The cone K i q consists of feasible ith-order motions. The dimension dim K i q is referred to as the ith-order DOF of the mechanism at q [17]. This is an infinitesimal DOF, iff K i then the linkage is called shaky of order i ¼ j À 1 at q 2 V. If q is a regular configuration, the mechanism is called underconstrained.
3. The inclusion in (39) is not strict. This is a critical point for the actual analysis since it is consequently not sufficient to increase the order until K i q ¼ K iÀ1 q . Moreover, the sufficient order j for a general mechanism is yet unknown.
Remark 5 A remark on the actual computation is in order. The mappings H i ð Þ q; x; y; . . . ð Þ in (38) are polynomials of order i in x, of order i À 1 in y, etc. The x satisfying the conditions in (38) are determined by eliminating y; z; . . ., etc. This is easily achieved symbolically invoking computer algebra software (e.g. Mathematica, Maple, Maxima). This is the actual advantage of the proposed method since the variety is locally replaced an algebraic variety of low degree.
Example 4 (Pin-in-Slot Joint, cont.) The planar mechanism in Fig. 5 is considered as the first example. The mechanism consists of a 4-bar linkage (bodies 0, 1, 2, 3) whose coupler point (at center of link 2) is restricted by the pin-in-slot higher pair joint so to move on a straight line. The geometry of the 4-bar linkage is set according to the parameters a ¼ 3; . Then its coupler curve is a sextic [10], i.e. it approximates a straight line up to order 5. Restricting the vertical coupler motion to the horizontal line yields an immobile but shaky linkage. The two FCs are chosen as in Fig. 5c. The configuration of the tree-topology system is described by the n ¼ 3 joint angles q 1 ; q 2 ; q 4 . The corresponding geometric constraints were derived in Sects. 2 and 3. The cut-joints only impose translation constraints. For the FC K 3 the geometric constraint is (8), with constraint mapping h 3 q ð Þ :¼ ða þ bÞs 1 þ dc 1 À bs 12 À d. The first two equations in (29) are the constraints for the FC K 5 with constraint mapping h 5 q ð Þ:¼ h 5;1 q ð Þ h 5;2 q ð Þ ¼ ða þ bÞc 1 À 2bc 12 À 2a þ ða þ bÞc 124 À ds 1 þ ds 124 dc 1 À dc 124 þ ða þ bÞs 1 À 2bs 12 þ ða þ bÞs 124

:
They are summarized in the overall system of constraints h q ð Þ ¼ 0. The joint screw coordinates in the reference configuration q 0 ¼ 0 in Fig. 5a, expressed in the world frame on B 0 , are pair so to move on a straight line. More on shakiness can be found in [17,29].
Example 5 (Spherical mechanism with 'in-line' joint) The last example is the single-loop mechanism in Fig. 7. It consists of three revolute (lower pair) joints and an higher pair 'in-line' joint. The latter constraints the relative translation to a line while not constraining the relative rotation. The three revolute joint axes intersect in one point so that the kinematic chain formed by joints 1, 2, 3 is a spherical linkage. The in-line joint 4 is used as cut-joint connecting bodies k 0 ¼ 2 and l 00 ¼ 3. The cut-joint joint constraints are imposed on the remaining joint variables q ¼ q 1 ; q 2 ; q 3 ð Þ . The RFRs F k 0 and F l 00 are at the center of the respective link as shown in Fig. 7. Shown is the reference configuration with q ¼ 0. All four links have the same length L. The world frame F 0 is located at the center of link 0, which is regarded as ground. The screw coordinates of the revolute joints represented in F 0 are The respective reference configuration of the two bodies is A k 0 ¼ ðI; À L 2 ; 0; 0 À Á T Þ and A l 00 ¼ ðI; 0; À L 2 ; 0 À Á T Þ. The configurations of the two bodies are determined by f k 0 q ð Þ ¼ expðY 1 q 1 Þ expðY 2 q 2 ÞA k 0 and f l 00 q ð Þ ¼ expðY 3 q 3 ÞA l 00 . The location vectors to the contact line defining the higher pair joint, expressed in the respective RFR, are p k 0 ¼ ð0; À L 2 ; 0Þ T ; p l 00 ¼ ðÀ L 2 ; 0; 0Þ T . The translation constraints imposed by the higher pair require the 1 and 3 component of the relative translation vector d k 0 in (3) to vanish giving rise to two geometric constraints. They are easily evaluated with (1). The two components of the constraint mapping (4) are Evaluating the second-and third-order cone shows that Furthermore, all higher-order cones are equal to K 1 q 0 , and the tangent cone is C q 0 V ¼ K 1 q 0 . This is the one-dimensional tangent vector space to V, and q 0 is a regular configuration. The local mobility of the mechanism is thus d loc q 0 ð Þ ¼ dim C q 0 V ¼ 1. This mechanism merely functions as a spherical mechanism despite the presence of a higher pair joint. Moreover, with the in-line joint the mechanism is not overconstrained. This is in contrast to the spherical 4-bar linkage (replacing the higher pair by another revolute joint whose axis intersects at the same point as the three revolute axes) that can be considered as being overconstrained since the closure constraints are redundant.

Conclusion
Any approach to determine the mobility of a holonomic mechanism requires the geometric constraints that the bodies are subjected to. Moreover, they must be formulated such that they can be evaluated efficiently. Since there is no general method to establish the global mobility, or to make global statements about the local finite mobilities, the mobility determination must generally resort to a local analysis. This relies on the higher-order constraints, i.e. the time derivatives of the geometric constraints (acceleration, jerk, jounce, etc. constraints). Approaches to the higher-order local mobility analysis for linkages comprising lower pair joints have been reported, and recently a recursive method was introduced that applies to general lower pair linkages. For such linkages the geometric loop constraints are given by the product of exponentials (POE), and the velocity constraints in terms of the instantaneous joint screws. Also the higher-order constraints are determined recursively by means of screw products of the instantaneous joint screws.
In this paper the higher-order mobility determination and the necessary constraints have been addressed for mechanisms with higher kinematic pairs. The formulation applies to general higher kinematic pairs defining bilateral geometric constraints. The formulation accounts for general higher pairs defining bilateral constraints, including cams. A recursive formulation for the time derivatives of the loop constraints is introduced for kinematic loops that contain one higher pair. To this end, the higher pair is removed and corresponding constraints are imposed. The remaining open kinematic chains contain only lower pairs. This gives rise to a recursive and computationally simple formulation of higher-order constraints. The basic operation involved is the screw product (Lie bracket), which is an extremely simple algebraic operation. For better accessibility, the formulation is first derived for a single loop and then generalized to multi-loop mechanisms.
This provides the basis for the local mobility analysis of higher pair mechanisms. A computational approach is introduced using the recursive constraint formulation. It allows for determination of the finite local mobility as well as for singularity analysis. This generally applicable approach complements the existing method for lower pair linkages. The method is applied to an immobile but shaky planar mechanism exhibiting 5th-order infinitesimal motion and to a spherical mechanism comprising an in-line higher pair joint, as examples. It should finally be mentioned that the higher-order mobility analysis in Sect. 5 is a selfcontained concept applicable as long as the higherorder constraints are available (not necessarily using the recursive formulation in Sect. 3).
Acknowledgments Open access funding provided by Johannes Kepler University Linz. The author acknowledges partial support of this work by the Austrian COMET-K2 program of the Linz Center of Mechatronics (LCM).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: rigid body motions-Lie group SE ð3Þ
The configuration of a body-fixed reference frame F 1 relative to a global frame F 0 is represented by where R 2 SO 3 ð Þ transforms the body-fixed representation of vector coordinates to their spatial representation, and r 2 R 3 is the position vector of the origin of the body-fixed frame expressed in the spatial frame. These matrices represent frame transformations. They also transform homogenous point coordinates from body-fixed RFR to global frame. let p 2 R 3 be the coordinate vector of a point expressed in the bodyfixed frame, and let p ¼ p; 1 ð Þ T be its homogeneous coordinates vector. Then the homogeneous coordinates of the point expressed in the global frame is s ¼ Cp. Indeed the corresponding vector is s ¼ Rp þ r. Whenever convenient the alternative notation C ¼ R; r ð Þ 2 SE ð3Þ is used. Details on the Lie group SE ð3Þ can be found in [24].
A general rigid body motion is a screw motion (Chasles theorem). Let e 2 R 3 be unit vector along the instantaneous screw axis, and p 2 R 3 the position vector of a point on that axis, both expressed in the world frame. The screw coordinate vector is then Y ¼ e; p Â e þ he ð Þ T . The finite motion corresponding to a screw motion about the axis e with angle u and pitch h is determined as C ¼ expðuYÞ where expðuYÞ ¼ expðueÞ ðI À expðueÞÞp þ hue 0 1 : Here expðueÞ is the Euler-Rodrigues formula for the rotation about the axis e with angle u. Y is referred to as the screw coordinate vector. The exp mapping (49) yields the finite rigid body motion corresponding to the motion of the instantaneous screw. If C 2 SE ð3Þ describes the frame transformation from frame F 1 to F 2 , and if X is the coordinate vector of a screw expressed in F 1 then the screw coordinates expressed in F 2 are determined by the Adjoint transformation Ad C X with where b r is the skew symmetric matrix associated to r. Writing screw coordinates as vector Y ¼ n; g ð Þ T , the screw product of two screws is Y 1 ; Y 2 ½ ¼ðn 1 Â n 2 ; g 1 Â n 2 À g 2 Â n 1 Þ: ð51Þ This defines the Lie bracket on se ð3Þ. To a screw coordinate vector X 2 R 6 is associated the matrix The spatial twist corresponding to a rigid body motion C t ð Þ is determined by ð 53Þ and in vector notation V ¼ x; _ r À b xr ð Þ T . Here x is the angular velocity of the body expressed in the world frame, and _ r À b xr is the translational velocity of a point on the body that is instantaneously traveling through the origin of the world frame.