## Abstract

This paper presents a geometric constraints driven approach to unified kinematic simulation of n-bar planar and spherical linkage mechanisms consisting of both revolute and prismatic joints. Generalized constraint equations using point, line, and plane coordinates have been proposed which unify simulation of planar and spherical linkages and are demonstrably scalable to spatial mechanisms. As opposed to some of the existing approaches, which seek to derive loop-closure equations for each type of mechanism separately, we have shown that the simulation can be made simpler and more efficient by using unified version of the geometric constraints on joints and links. This is facilitated using homogeneous coordinates and constraints on geometric primitives, such as point, line, and plane. Furthermore, the approach enables simpler programming, real-time computation, and ability to handle any type of planar and spherical mechanism. This work facilitates creation of practical and intuitive design tools for mechanism designers.

## 1 Introduction

Kinematic simulation of a mechanism requires calculation of the position and orientation of all of its constituent links during its entire range of motion. Simulation methodologies can be broadly classified into three categories: graphical, analytical, and numerical [1]. The graphical analysis method is based on dyadic decomposition, i.e., identification of four-bar loops in mechanisms [2]. Although this approach is prominently used in simulation packages like Linkages [3] and PMKS [4], its limitations are well known [5], e.g., they are unable to handle complex mechanisms like a double butterfly mechanism. Analytical methods involve solving a loop closure constraint-based system of non-linear equations [6]. Most analytical methods use the polynomial continuation method [7,8], elimination method, or Grobner bases [9] to solve the simulation problem. These methods are able to find all the possible assembly configurations of a given mechanism. However, they are not general in nature since the motion equations need to be derived for each type of mechanism separately. As a result, these approaches cannot be used to simulate n-bar planar or spherical mechanisms without manually deriving equations on a case-by-case basis.

Numerical simulation methods are iterative in nature and can handle extremely complex mechanism [10,11]. They use numerical methods such as the Newton–Raphson method to solve the system of non-linear equations for one solution only instead of all possible ones [12]. These methods accept the mechanism joints and link information as inputs. Subsequently, the algorithm repeatedly solves the finite displacement problem, i.e., the input link is iteratively moved with finite displacement, and consequently, positions of remaining links are calculated. As a result, the entire range of motion for the specified mechanism is calculated.

Hernández and Petuya have proposed a geometrical-iterative method that performs better than Newton–Rhapson method [13]. However, the approach is limited to n-bar planar mechanism with revolute joints only. Radhakrishnan and Campbell [14] have created a computational tool for planar mechanism, which carries out position analysis of planar mechanisms using a geometrically iterative algorithm. However, due to the use of dyadic decomposition, it shares the limitations of graphical methods and is limited to the planar mechanisms.

Commercial CAD software such as autodesk inventor, solidworks, adams, etc. solve differential-algebraic equations numerically to provide multibody simulation capability [1517].2 However, their use is more prominent during detail design stage rather than the conceptual design stage. Creation of feature-based assembly of planar and spherical mechanisms and initializing constraints on these systems is a nontrivial task. Changing the type of joints or the number of links for a mechanism is also more involved than carrying out the same operation on purely kinematic simulators like PMKS [4]. Additionally, their solvers model the motion problem as a set of coupled differential and algebraic equations. This type of model is more suited for dynamic simulations rather than kinematic simulations, which involves purely algebraic constraints. Also, the algebraic equations for commercial softwares are modeled using reference point representation that leads to more number of constraints when compared with other representations. Thus, use of these softwares for concept design is not ideal. SAM and GIM are two more packages that support n-bar simulation for planar linkages with both revolute and prismatic joints [18].3 Furlong et al. [19] have demonstrated a virtual reality environment for simulating spherical four-bar mechanisms, and in the academic domain, sphinx, isis, and osiris are softwares that enable the analysis and synthesis of spherical mechanisms [2022].

However, currently, there are no approaches that unify n-bar planar and spherical mechanism analysis and can be demonstrated to more complex linkage systems. The proposed approach hopes to bridge this gap and augment the capability of pure kinematic design systems like MotionGen [23]. In this paper, planar and spherical mechanisms are represented as a collection of geometric constraints spanned by points, lines, and planes. The geometric mechanism representation enables the design of unified constraint equations that are easily programmed. As a result, a simple generalized real-time framework for mechanism simulation is achieved.

Our previous work on mechanism synthesis problems includes the creation of geometric constraint equations for four-bar mechanisms with revolute or prismatic joints [2428]. In Ref. [29], we have demonstrated an algebraic fitting-based approach in the space of planar quaternions to simulate planar four-bar linkages. However, that approach does not scale for more complex planar or spherical linkages. In this paper, we show that by using homogeneous coordinates, we can derive unified geometric constraint equations for both planar and spherical linkages, which simplifies the simulation without resorting to calculations for individual types of mechanisms. The rigidity constraints imposed by the links are modeled as simple geometric constraints using points, lines and planes. Once the mechanism is specified, the solver proceeds with iteratively perturbing the input and solving the constraints for other links. To the best of authors’ knowledge, this work is the first attempt at using a point-line-plane mechanism representation and presenting unified geometric constraints for simulation.

The major intellectual contributions of this paper can be summarized as (1) presenting generalized constraint equations for planar and spherical mechanisms using point-line-plane representation and (2) enabling real-time simulation of n-bar planar or spherical mechanisms.

Rest of the paper is organized as follows. Section 2 discusses the representation and constraints required to describe the motion of a general planar and spherical mechanism. Section 3 demonstrates the algorithm required to simulate a mechanism using the iterative numerical approach. Finally, in Sec. 4, we present a few examples to demonstrate the use of proposed algorithm.

## 2 Mechanism Representation and Constraints

Selection of an apt mechanism representation and constraints is important as it has a profound effect on algorithm’s simplicity and efficiency. Conventionally, a multibody system has been specified using multiple representations, namely: relative coordinates, reference point coordinates, and natural coordinates [10]. Relative coordinates are based on parameters specifying one link relative to another; reference point coordinates are based on specifying absolute position of each link independently; while the natural coordinates are based on each link being specified by two points. Using relative coordinates enables a scalable representation while reference point coordinates tend to be more computationally efficient. Natural coordinates provide a compromise between the two approaches in terms of simplicity and efficiency. Most commercial softwares use reference point coordinate representation that usually leads to maximum number of constraint equations and subsequently high computation time.

Figure 1 shows an RRPR (R: Revolute, P: Prismatic) four-bar mechanism and its specification using different representations. For the relative coordinate representation, there are three unknown coordinates, i.e., ψ1, ψ2, L3. For the reference point coordinate representation, the mechanism has nine unknown variables, i.e., location and orientation of each link xi, yi, ψi. Similarly, for the natural coordinate representation, there are six unknown variables namely x1, y1, x2, y2, x3, y3. Since the four-bar mechanisms are a single degree-of-freedom mechanisms, each of the representation requires two, eight, and five constraint equations, respectively, to fully define the motion. In this paper, we derive unified constraint equations for all types of planar and spherical linkages consisting of both revolute and prismatic joints. We use homogeneous coordinates to write geometric constraints on points, lines, and planes. For example, our representation for the shown RRPR mechanism will require using four unknown point and line coordinates, i.e., x1, y1, a2, b2 since we can set homogenizing factors z1 = 1 and c2 = 1 without any loss in generality. Such a representation would keep the number of unknown variables smaller while also enabling construction of simpler geometric constraint equations. Computational efficiency aside, this representation also naturally maps to the geometric constraints, which for the shown RRPR mechanism is a circle-constraint on the moving pivot and a line-constraint on the fixed pivot of the RP dyad. It is well known that the time complexity of multidimensional Newton–Raphson method, which is used in this paper, is at-least O(n2) for a single iteration where n is the number of unknowns in state variable [11]. This is a direct consequence of a dense Jacobian matrix having n2 elements that need to be calculated after every iteration. Thus, more unknowns result in needing to calculate larger Jacobian matrices, which is computationally expensive.

Fig. 1
Fig. 1

Planar mechanisms can be uniquely specified using their joint and link data. A joint can be prismatic or revolute, which naturally associates with points and lines, respectively. We use homogeneous coordinates to represent both points and lines. Thus, a point P is given by homogeneous coordinates (x, y, z) whose Affine coordinates are given as (x/z, y/z), while a line L is also represented using homogeneous coordinates (a, b, c), where equation of the line passing through the point P in the projective plane is given by ax + by + cz = 0. Depending on the constraint being expressed, this line can be fixed or floating in the plane. A link can be represented by a subset of joints. The link can be binary, ternary, or n-ary depending on the number of joints it contains. An example six-bar planar mechanism is displayed in Fig. 2. Its joints are represented as points and lines while its links are defined as a group of joints as shown in Table 1.

Fig. 2
Fig. 2
Table 1

JointTypeCoordinates
J1,inputRevolute0, −1
J2Revolute1, 0.5
J3Prismatic−0.17, 0.98, −4.28
J4Revolute3.25, 1.4
J5Revolute7.72, 1.44
J6Revolute11.66, 4.17
J7Prismatic0, 1, 1.24
J8Coupler point6, −2
JointTypeCoordinates
J1,inputRevolute0, −1
J2Revolute1, 0.5
J3Prismatic−0.17, 0.98, −4.28
J4Revolute3.25, 1.4
J5Revolute7.72, 1.44
J6Revolute11.66, 4.17
J7Prismatic0, 1, 1.24
J8Coupler point6, −2
L1J1, J2
L2J2, J3, J4
L3J3, J6
L4J4, J5, J8
L5J5, J6, J7
L6,groundJ1, J7
L1J1, J2
L2J2, J3, J4
L3J3, J6
L4J4, J5, J8
L5J5, J6, J7
L6,groundJ1, J7

Similarly, spherical mechanisms can also consist of revolute and prismatic joints. A spherical prismatic joint constrains the link movement along a circular arc instead of a line. We represent a spherical revolute joint as a point P in terms of its homogeneous coordinates (x, y, z, w) with respect to the center of the unit sphere such that its Affine coordinates are (x/w, y/w, z/w). A spherical prismatic joint is defined as a plane Pl : (a, b, c, 0) passing through the center of the sphere and is given by the equation ax + by + cz = 0. The intersection of the plane and the unit sphere defines the great circle along which the constituent links are constrained to move for a spherical prismatic joint. Similar to the lines for planar mechanisms, this plane can be fixed or moving depending on the geometric constraint being expressed. An example RRPR spherical mechanism is displayed in Fig. 3. Its joints are represented as points and planes while its links are defined as a group of joints as shown in Table 2.

Fig. 3
Fig. 3
Table 2

Joint and link data for spherical RRPR linkage using Affine coordinates; the coordinates are given with respect to the fixed coordinate frame located at the center of the reference sphere

JointTypeCoordinates
J1,inputRevolute0.94, 0.24, 0.24
J2Revolute0.80, 0.27, 0.53
J3Prismatic0.68, −0.68, 0.26
J4Revolute−0.38, 0.76, 0.53
J5Coupler point0.50, −0.21, 0.84
JointTypeCoordinates
J1,inputRevolute0.94, 0.24, 0.24
J2Revolute0.80, 0.27, 0.53
J3Prismatic0.68, −0.68, 0.26
J4Revolute−0.38, 0.76, 0.53
J5Coupler point0.50, −0.21, 0.84
L1J1, J2
L2J2, J3, J5
L3J3, J4
L4,groundJ1, J4
L1J1, J2
L2J2, J3, J5
L3J3, J4
L4,groundJ1, J4

During the motion, a mechanism is subjected to a set of constraints imposed by the rigidity of its links. Thus, to simulate a mechanism, these constraint equations need to be formulated. For planar and spherical mechanisms, modeling three constraint equations are sufficient for simulation. We propose a unique constraint equation for each of the binary links RR, RP, PR, and PP. But, we will see that all of these constraints can be expressed in a single equation. Any link with more than two joints can easily be reduced to an equivalent collection of binary links. For example, a ternary link can be treated as three binary links. Thus, these constraints can successfully be used to enforce the rigidity of any link in a general mechanism. Figures 4 and 5 show different planar and spherical binary links, which are building blocks for any planar and spherical mechanism and are discussed below.

Fig. 4
Fig. 4
Fig. 5
Fig. 5
The first general constraint enforces the rigidity of a spherical binary link with two revolute joints represented by two homogeneous point coordinates of the fixed point (a1, a2, a3, a4) and floating point (b1, b2, b3, b4), where a4 and b4 are homogenizing factors. The RR link imposes the constraint that the distance between two points remains constant, i.e., dist(R1, R2) = r in Figs. 4 and 5. The constraint equation is given as
$CS,RR:2a1b1+2a2b2+2a3b3+a0b4=a4(b12+b22+b32b4)$
(1)
where a0 is given as
$a0=a4r2−a12+a22+a32a4$
Here, r is the radius of the sphere formed by the RR link with the center given by (a1, a2, a3, a4). When the z-coordinate is set to zero, the constraint equation degenerates into the one for a planar RR link. The constraint equation for a planar RR link represented by points (a1, a2, a4) and (b1, b2, b4) is thus given as
$CP,RR:2a1b1+2a2b2+a0b4=a4(b12+b22b4)$
(2)
where a0 is given as
$a0=a4r2−a12+a22a4$
The second general constraint enforces the rigidity of a binary link with one prismatic and one revolute joint represented by a homogeneous point and a plane given by (a1, a2, a3, a4) and (L1, L2, L3, L4), respectively. An RP or PR link imposes the constraint that the distance between a point and a line (planar case) or a point and a plane (spherical case) is constant, i.e., dist(R, P) = d in Figs. 4 and 5. RP and PR links are inversions of each other and are expressed by the same constraint. The general constraint equation for a spherical RP link is given as
$CS,RP:a1L1+a2L2+a3L3+a4L4=da4L12+L22+L32$
(3)
where d is the signed perpendicular distance between the revolute joint and prismatic joint. For spherical linkages, the prismatic plane always passes through the origin, i.e., L4 = 0. Thus, the spherical RP link constraint equation reduces to
$CS,RP:a1L1+a2L2+a3L3=da4L12+L22+L32$
(4)
When the z-coordinates is set to zero, the general constraint equation degenerates into a planar case. Thus, constraint equation for a planar RP or PR link represented by a point (a1, a2, a4) and a line (L1, L2, L4) is given as
$CP,RP:a1L1+a2L2+a4L4=da4L12+L22$
(5)
When the perpendicular distance d becomes zero, the equation describes a line passing through a point, i.e., constraint equation of PR or RP links. This is usually the case with the PR links where the moving joint point is constrained to be on the fixed line of the prismatic joint and in case of the RP link where the moving line is constrained to pass through the point of the fixed joint.
Finally, the third constraint enforces the rigidity of a spherical binary link with two prismatic joints represented as (L1, L2, L3, 0) and (M1, M2, M3, 0). For the PP link, the angle between two lines (planar case) or two planes (spherical case) remains constant, i.e., angle(P1, P2) = cos −1(k) in Figs. 4 and 5. The constraint equation is given as
$CS,PP:L1M1+L2M2+L3M3=kL12+L22+L32M12+M22+M32$
(6)
where k represents the cosine of angle between two prismatic joints. Similarly, for planar PP binary link represented by two lines (L1, L2, L4) and (M1, M2, M4), the constraint equation degenerates to
$CP,PP:L1M1+L2M2=kL12+L22M12+M22$
(7)
When the two prismatic joints on a binary link are defined as two parallel lines, a degree-of-freedom is added to the mechanism. This situation is impractical and will not been considered further in this paper.

It can be seen that Eqs. (2), (4)(7) are all degenerate case of the Eq. (1). In the projective plane for the planar geometric constraints, the lines and points are dual to each other; thus, their meanings can be interchanged without changing the underlying structure of the equations. In the projective three-space for the spherical constraints, the points and planes are dual to each other and thus their meanings can be interchanged. Thus, Eq. (1) is the single equation that unifies all the geometric constraints associated with all types of links for both planar and spherical mechanisms. This facilitates creation of the following metrics for computation: (1) distance between two points in space, (2) perpendicular distance between a point and a plane, and (3) angle between two planes.

For links with prismatic joints, the line or plane coordinates are homogeneous in nature, i.e., multiplying a non-zero scalar λ to prismatic coordinates (L1, L2, L3) does not change the coordinates. Thus, the magnitude of this vector can be fixed to unity without losing generality and another constraint can be written as
$CP:L12+L22+L32−1=0$
(8)
For spherical mechanisms, an additional geometric constraint is imposed on the joints due to the spherical nature of the motion. It is assumed that all the revolute joints move on the unit sphere that leads to the constraint
$CS,R:a12+a22+a32−1=0$
(9)
where (a1, a2, a3) are the coordinates of any revolute joint on a spherical mechanism. Thus, the rigidity constraints described in Eqs. (1), (2), (4)(9) are sufficient to uniquely determine the unknown coordinates of a n-bar planar or spherical mechanism. This concludes our discussion on representation and constraints for a generalized planar or spherical mechanism.

## 3 Solving Constraint Equations

In this section, we discuss the algorithmic steps required to solve the kinematic simulation problem. The general approach is to iteratively perturb the input links by a finite displacement and find the new position of the mechanism.

The simulation process involves iteratively perturbing the input link by a finite displacement. Depending on the actuating joint being revolute or prismatic, the displacement could be translation or rotational in nature. In this paper, we restrict ourselves to consider actuation at the fixed joints. The relations governing the motion of input link are derived in this subsection.

For a perturbed RR link with the actuating fixed joint (x1, y1) and moving joint (x2, y2), the new coordinates of moving revolute joint can be given as
$[X2Y21]=[T]−1[R][T][x2y21]$
(10)
where
$[Rx]=[cos(θ)−sin(θ)0sin(θ)cos(θ)0001]and[T]=[10−x101−y1001](11)$
(11)
In the above equation, (X2, Y2) represent the moving joint after perturbation and θ is the angle through which the input link is perturbed.
For a perturbed RP link with the actuating fixed joint (x1, y1) and moving joint (a, b, c), the new coordinates of moving line representing the prismatic joint can be given as
$[ABC]=([T]−1[R][T])−T[abc]$
(12)
where (A, B, C) are the moving line coordinates after perturbation, and T and Rx are the translation and rotation matrices as described in Eq. (11).
For a planar mechanism with the actuation being at prismatic joint, input link perturbation causes translation of other joints on the input link. For a perturbed PR link with the actuating joint (a, b, c) and moving joint (x, y), the new coordinates of translating revolute joint can be given as
$[XY1]=[ba2+b2d−aa2+b2d0]+[xy1]$
(13)
where (X, Y) are the moving joint coordinates after perturbation and d is the distance through which the prismatic joint is moved along the fixed line. It can be seen that in Eq. (13), the actuating line coordinate c does not affect the new position of moving joint coordinates as new position only depends on direction cosines.
For a perturbed PP link with the actuating fixed joint (a1, b1, c1) and moving joint (a2, b2, c2), the new coordinates of translating prismatic joint can be given as
$[A2B2C2]=[00a1b2−a2b1a12+b12d]+[a2b2c2]$
(14)
where (A2, B2, C2) are the moving prismatic joint coordinates after perturbation and d is the distance through which the input link has been perturbed.
Similarly, relationships determining the values of perturbed joints for spherical mechanisms can also be calculated. For spherical mechanisms with a fixed revolute actuating joint, the moving joints rotate around the axis passing through the actuation joint and the centre of sphere. The transformation matrix which rotates spherical link around an axis passing through the centre of sphere (0, 0, 0) and an arbitrary point on surface of the sphere (l, m, n) is given by
$[R](l,m,n)=[Rx]−1[Ry]−1[Rz][Ry][Rx]$
(15)
$[Rx]=[1000nm2+n2−mm2+n20mm2+n2nm2+n2]$
(16)
$[Ry]=[m2+n20−l010l0m2+n2]$
(17)
$[Rz]=[cos(θ)−sin(θ)0sin(θ)cos(θ)0001]$
(18)
where Rx, Ry, and Rz are the rotation matrix around x-, y-, and z-axis and θ is the angle by which the link is rotated around the axis.
Using Eq. (15), the new coordinates of the moving joints of a perturbed spherical RR link can be given as
$[X2Y2Z2]=[R](x1,y1,z1)[x2y2z2]$
(19)
where (x1, y1, z1) are the fixed joint coordinates, (x2, y2, z2) are the moving joint coordinates before perturbation, and (X2, Y2, Z2) are the moving joint coordinates after perturbation.
For a spherical RP link with a fixed revolute joint, the coordinates of moving prismatic joint can be given as
$[ABC]=[R](x,y,z)[abc]$
(20)
where (x, y, z) are the fixed joint coordinates, (a, b, c) are the moving joint coordinates before perturbation, (A, B, C) are the moving joint coordinates after perturbation, and as described above.
When the actuation joint is prismatic in nature, the moving joints translate on the intersection of a parallel plane and the unit sphere. This motion can also be characterized as rotation around an axis which passes through the centre of the sphere and “pole” of the prismatic joint. The poles of a great circle are defined as intersection of two circles perpendicular to the initial circle. If a spherical prismatic joint is defined as a plane (a, b, c), its pole coordinates are also given as (a, b, c). Thus, for a spherical RP link with fixed prismatic joint, the coordinates of moving revolute joint can be given as rotation, i.e.,
$[XYZ]=[R](a,b,c)[xyz]$
(21)
where (a, b, c) are the fixed prismatic joint coordinates, (x, y, z) are the moving revolute joint coordinates before perturbation and (X, Y, Z) are the moving revolute joint coordinates after perturbation.
For a spherical PP link with a fixed prismatic joint, the coordinates of a moving prismatic joint can be given as
$[A2B2C2]=[R](a1,b1,c1)[a2b2c2]$
(22)
where (a1, b1, c1) are the coordinates of fixed prismatic joint, (a2, b2, c2) are moving prismatic joint coordinates before perturbation, and (A2, B2, C2) are moving prismatic joint coordinates after perturbation.

With these expressions, we can successfully calculate the location of input link after imparting it a discrete perturbation. The next step is to find the coordinates of all the other unknown joint coordinates that are compatible with the rigidity constraints imposed on the mechanism during simulation.

### 3.2 Numerical Non-linear System of Equation Solving.

For any multi-body system, the position problem is always based on solving a system of constraint equations. This set of equations can be represented as
$Φ(q)=0$
(23)
where q is the state vector that consists of all the unknown coordinates. The well-known Newton-Raphson method can be used to solve this non-linear system of equation. It is featured by quadratic convergence in the neighborhood of the solution. Since the input link is perturbed by a small finite displacement, the previous state of mechanism serves as a good initial approximation. The number of constraint equations should be equal to or greater than the number of unknowns for this approach to work. For planar and spherical mechanisms, it is always possible to satisfy this criterion using the constraints outlined in section.
The iterative algorithm followed can be defined as
$qi+1=qi−[J−1(qi)]Φ(qi)$
(24)
where qi is the state vector at ith iteration, Φ(qi) is the vector of residuals at q = qi, and [J−1(qi)] is the inverse of Jacobian matrix evaluated at q = qi. The Jacobian matrix is of the following form
$[J(q)]=[∂ϕ1∂q1∂ϕ1∂q2⋯∂ϕ1∂qn∂ϕ2∂q1∂ϕ2∂q2⋯∂ϕ2∂qn⋯⋯⋯⋯∂ϕm∂q1∂ϕm∂q2⋯∂ϕm∂qn]$
(25)
where m is the number of constraints and n is the number of unknown coordinates. Thus to calculate the Jacobian matrix, relations describing the first-order partial derivatives of constraint equations are required.
The first-order partial derivatives for spherical RR constraint given in Eq. (1) can be given as follows:
$∂CRR∂a1=2(a1−b1)$
(26)
$∂CRR∂a2=2(a2−b2)$
(27)
$∂CRR∂a3=2(a3−b3)$
(28)
Here, the homogeneous point coordinate a0 and b4 has been assumed as unity without loss in generality. For planar RR link, only ∂CRR/∂a1 and ∂CRR/∂a2 exists.
The first-order partial derivatives for spherical RP constraint given in Eq. (4) can be given as follows:
$∂CS,RP∂a1=L1,∂CS,RP∂a2=L2,∂CS,RP∂a3=L3$
(29)
$∂CS,RP∂L1=a1−dL1L12+L22+L32$
(30)
$∂CS,RP∂L2=a2−dL2L12+L22+L32$
(31)
$∂CS,RP∂L3=a3−dL3L12+L22+L32$
(32)
Here, the homogeneous point coordinate a0 has been assumed as unity without loss in generality. For planar RR constraint given in Eq. (5), the first-order differentials are
$∂CP,RP∂a1=L1,∂CP,RP∂a2=L2$
(33)
$∂CP,RP∂L1=a1−dL1L12+L22$
(34)
$∂CP,RP∂L2=a2−dL2L12+L22$
(35)
$∂CP,RP∂L4=1$
(36)
The first-order partial derivatives for spherical PP constraint given in Eq. (6) can be given as follows:
$∂CS,PP∂L1=M1−kL1M12+M22+M32L12+L22+L32$
(37)
$∂CS,PP∂L2=M2−kL2M12+M22+M32L12+L22+L32$
(38)
$∂CS,PP∂L3=M3−kL3M12+M22+M32L12+L22+L32$
(39)
These equations degenerate to planar case when L3 = 0 and M3 = 0 and ∂CP,PP/∂L4 = 0.
The first-order partial derivatives for homogeneous Prismatic joint constraint given in Eq. (8) can be given as follows:
$∂CP∂L1=2L1,∂CP∂L2=2L2,∂CP∂L3=2L3$
(40)
The first-order partial derivatives for unit circle revolute joint constraint given in Eq. (9) can be given as follows:
$∂CS,R∂a1=2a1,∂CS,R∂a2=2a2,∂CS,R∂a3=2a3$
(41)

To automate the calculation of residual vector Φ(qi) and the Jacobian matrix [J(qi)], the constraints are handled in a sequential manner. While creating the residual vector in our implementation, first the rigidity constraints for each link are calculated and then the constraints for joints are calculated. Similarly, the Jacobian matrix is created in a column-first manner, i.e., all the partial differential equations with respect to an unknown state variable are calculated before progressing to the next variable. The outlined method is just one way of calculating Φ(qi) and [J(qi)] since their values are independent of the sequence adopted to calculate each element.

Thus, using the constraint equations and their first-order partial derivatives, it is possible to solve iteratively for the solution using Newton–Raphson method. The iterations are continued until a solution within desired accuracy is calculated.

For some input link perturbations, the Newton–Raphson method might fail to converge even after many iterations. In these instances, there does not exist a mechanism state which fulfills all the constraint equations. As a result, this input perturbation is outside the possible limits of motion of mechanism.

Thus, by iteratively perturbing the input link and solving the constraints for other joint coordinates, we are able to simulate any general planar or spherical mechanism. Numerous techniques exist that can improve the convergence and efficiency of the Newton–Raphson method. However, the basic method suffices to achieve real-time simulation. The complete algorithm has been described in

Algorithm 1.

## 4 Examples

This section presents sample examples to demonstrate the use of proposed algorithm for mechanism simulation. The simulation has been carried out in matlab on a PC running Core i5-7300 at 2.6 GHz with 8 GB RAM. The simulation is carried out within seconds for residual value of 1.0e − 8. Each closed-loop output curve is made up of 180 points while open-loop curves have less than 180. Each point corresponds to 2π/180 radian or units input perturbations.

### 4.1 Speed Comparison With a Commercial Software.

To compare the speed of commercially available CAD systems and proposed algorithm, a planar four-bar crank rocker mechanism with revolute joints is modeled and simulated. Autodesk Inventor 2020 with educational license is used as reference commercial CAD system. Simulation is performed for 180 time-steps with one-degree input link perturbation for each time-step. The simulation takes 5 s to complete Inventor while it finishes in 1.4 s when the proposed algorithm is used. Thus, even for a simple mechanism, there is a significant speed difference between the commercial solvers and proposed methodology.

A planar Stephenson-II six-bar linkage is simulated in this example. This linkage does not have a four-bar linkage and proves to be challenging to simulate using dyadic decomposition based approaches. However, our approach handles these non-dyadic mechanisms without issues.

The six-bar mechanism is displayed in Fig. 2, and its joint and link data are given in Table 1. The mechanism has J1, J7 as the fixed joints, J2 as the perturbed joint and J3, J4, J5, J6, J8 as the unknown joints defining the 11-dimensional state vector. The mechanism consists of ten rigidity constraint equations and one homogeneous coordinate equation for prismatic joint. The simulation algorithm successfully solves these constraints and plots the trajectory of the coupler point J8 as shown in Fig. 2. The run-time of this simulation was 3.79 s.

### 4.3 Planar-Modified Theo Jansen Linkage.

In this example, a planar-modified Theo Jansen linkage with one of its revolute joints replaced by a floating prismatic joints is simulated. The eight-bar mechanism is displayed in Fig. 6 and its joint and link data is given in Table 3. J1, J5 are the fixed joints, J2 is the perturbed joint and the state vector consists coordinates of J3, J4, J6, J7, J8. This results in a 11-dimensional state vector. Ten rigidity constraint equations for links and one homogeneous coordinate equation for prismatic joint are available for this mechanism. The simulation algorithm plots the trajectory of the coupler point J8 as shown in Fig. 6. Note, the length of stride for this modified mechanism is larger than that of the conventional Theo Jansen mechanism which has revolute joints only. As a result, this mechanism is a prospective candidate for walking robots. The run-time of this simulation was 3.81 s.

Fig. 6
Fig. 6
Table 3

JointCoordinates
J1,input2.77, 2.31
J22.17, 3.33
J3−0.50, 0.87, −4.80
J40.66, −1.3
J5−0.22, 1.72
J6−3.17, 0.66
J7−2.08, −2.24
J82.54, −4.64
JointCoordinates
J1,input2.77, 2.31
J22.17, 3.33
J3−0.50, 0.87, −4.80
J40.66, −1.3
J5−0.22, 1.72
J6−3.17, 0.66
J7−2.08, −2.24
J82.54, −4.64
L1J1, J2
L2J2, J3
L3J2, J4
L4J3, J5, J6
L5J5, J4
L6J6, J7
L7J4, J7, J8
L8,groundJ1, J5
L1J1, J2
L2J2, J3
L3J2, J4
L4J3, J5, J6
L5J5, J4
L6J6, J7
L7J4, J7, J8
L8,groundJ1, J5

### 4.4 Spherical RRPR Mechanism.

This example presents the simulation of a spherical RRPR mechanism which is the spherical analog of the Whitworth quick-return mechanism. The four-bar mechanism is displayed in Fig. 3, and its joint and link data are given in Table 2. The mechanism has J1, J4 as the fixed joints, J2 as the perturbed joint, and J3, J5 as the unknown joints defining the six-dimensional state vector. The mechanism consists of four rigidity constraint equations for output and coupler links which can be described using Eqs. (1) and (4). Also, one unit sphere equation for revolute joint using Eq. (9) and one homogeneous coordinate equation for prismatic joint using Eq. (8) can be written. Once the simulation is completed, the trajectory of coupler point J5 can be plotted as shown in Fig. 3. The run-time of this simulation was 1.49 s.

In this example, a spherical Watt-I six-bar linkage with prismatic input joint is simulated. Spherical Watt I type linkages have been used to design door hinges for spatial movement. The six-bar mechanism is shown in Fig. 7, and its link and joint data are given in Table 4. From the data, it is known that J1, J6 are the fixed joints, J2, J3 are the perturbed joints, and J4, J5, J7, J8 are the unknown joints representing the 12-dimensional state vector. The mechanism can be described using eight rigidity constraints for links and four unit circle constraints. Perturbing the input link along the input prismatic joint results in the motion of coupler point J8 as shown in Fig. 7. The run-time of this simulation was 3.33 s.

Fig. 7
Fig. 7
Table 4

JointCoordinates
J1,input0, 0, 1
J20.93, 0, 0.37
J30.85, −0.17, 0.51
J40.70, 0.70, 0.14
J50.73, 0.49, 0.49
J60.81, 0.41, −0.41
J70.48, −0.10, 0.87
J80.49, 0.49, 0.73
JointCoordinates
J1,input0, 0, 1
J20.93, 0, 0.37
J30.85, −0.17, 0.51
J40.70, 0.70, 0.14
J50.73, 0.49, 0.49
J60.81, 0.41, −0.41
J70.48, −0.10, 0.87
J80.49, 0.49, 0.73
L1J1, J2, J3
L2J2, J4, J5
L3J4, J6
L4J3, J7
L5J5, J7, J8
L6,groundJ1, J6
L1J1, J2, J3
L2J2, J4, J5
L3J4, J6
L4J3, J7
L5J5, J7, J8
L6,groundJ1, J6

### 4.6 Spatial 5-SS Mechanism.

In this section, we demonstrate the scalability of proposed algorithm to spatial mechanisms by simulating a 5-SS platform linkage. A 5-SS mechanism consists of five binary spherical-spherical (SS) links connected to a floating coupler link on one end and the ground link on the other [30]. Although using the Grüebler criterion [2], the mobility of this mechanism is six, five of the rotational degrees for each binary link are redundant, and therefore, 5-SS platform linkage is a one degree-of-freedom mechanism.

An example 5-SS mechanism is displayed in Fig. 8 and its joint and link data is given in Table. 5. Since the input link perturbation approaches outlined in Sec. 3.1 only cover revolute and prismatic joints, a new approach is needed for spherical input links. We attach a linear actuator between J1 and J7 to actuate the mechanism [30]. This results in the mechanism having J1, J2, J3, J4, J5 as fixed joints and J6, J7, J8, J9, J10, J11 as the unknown joints defining the 18-dimensional state vector. Five rigidity constraints for binary links, 12 independent rigidity constraints for coupler and one additional constraint for input actuator length is available for this mechanism. All these constraints are modeled using Eq. (1) and the simulation is carried out to generate a spatial trajectory of coupler point J11. The input prismatic link is perturbed by 0.01 units and run-time of this simulation was 0.28 s.

Fig. 8
Fig. 8
Table 5

JointCoordinates
J1−7.2, −5.29, 7.87
J25.72, −0.71, −8.26
J3−8.75, −4.60, 5.49
J4−10.00, −8.72, 6.09
J5−3.88, 6.61, 8.91
J6−6.61, −9.99, 3.93
J78.19, −9.05, 6.07
J87.84, −5.73, −0.62
J94.15, −0.48, −2.04
J10−0.97, −9.04, 1.92
J112.20, −7.17, 5.36
JointCoordinates
J1−7.2, −5.29, 7.87
J25.72, −0.71, −8.26
J3−8.75, −4.60, 5.49
J4−10.00, −8.72, 6.09
J5−3.88, 6.61, 8.91
J6−6.61, −9.99, 3.93
J78.19, −9.05, 6.07
J87.84, −5.73, −0.62
J94.15, −0.48, −2.04
J10−0.97, −9.04, 1.92
J112.20, −7.17, 5.36
L1J1, J6
L2J2, J7
L3J3, J8
L4J4, J9
L5J5, J10
L6J6, J7, J8, J9, J10, J11
L7,groundJ1, J2, J3, J4, J5
L1J1, J6
L2J2, J7
L3J3, J8
L4J4, J9
L5J5, J10
L6J6, J7, J8, J9, J10, J11
L7,groundJ1, J2, J3, J4, J5

## 5 Conclusion

In this paper, we have presented unified equations for motion simulation of planar and spherical n-bar mechanisms and an efficient algorithm for computation to enable real-time, interactive simulation. The approach is general and uses simple geometric primitives, such as point, line, and planes to represent the constraints inherent in mechanisms. A 5-SS mechanism is simulated to demonstrate the scalability of the proposed approach to spatial mechanisms. Once the mechanism is simulated and the path of coupler point determined, velocity and acceleration curves can easily be determined using numerical differentiation. Future research would involve finding appropriate representation and rigidity constraints for cylindrical and helical joints to further unify spatial synthesis.

## Acknowledgment

This work has been financially supported by The National Science Foundation under a research grant # CMMI-1563413 to Stony Brook University. All findings and results presented in this paper are those of the authors and do not represent those of the funding agencies.

## References

1.
Norton
,
R.
,
2011
,
Design of Machinery: An Introduction to the Synthesis and Analysis of Mechanisms and Machines
, 5th ed.,
McGraw Hill
,
New York
.
2.
Erdman
,
A. G.
, and
Sandor
,
G. N.
,
1991
,
Mechanism Design: Analysis and Synthesis
, 2nd ed., Vol.
1
,
Prentice-Hall
,
Englewood Cliffs, NJ
.
4.
Campbell
,
M.
, “
Planar Mech anism Kinematic Simulator
”, http://design.engr.oregonstate.edu/pmksintro.html
5.
Waldron
,
K.
, and
Sreenivasan
,
S.
,
1996
, “
A Study of the Solvability of the Position Problem for Multi-Circuit Mechanisms by Way of Example of the Double Butterfly Linkage
,”
ASME J. Mech. Des.
,
118
(
3
), pp.
390
395
. 10.1115/1.2826898
6.
Nielsen
,
J.
, and
Roth
,
B.
,
1999
, “
On the Kinematic Analysis of Robotic Mechanisms
,”
Int. J. Rob. Res.
,
18
(
12
), pp.
1147
1160
. 10.1177/02783649922067771
7.
Wampler
,
C. W.
,
1999
, “
Solving the Kinematics of Planar Mechanisms
,”
ASME J. Mech. Des.
,
121
(
3
), pp.
387
391
. 10.1115/1.2829473
8.
Nielsen
,
J.
, and
Roth
,
B.
,
1999
, “
Solving the Input/Output Problem for Planar Mechanisms
,”
ASME J. Mech. Des.
,
121
(
2
), pp.
206
211
. 10.1115/1.2829445
9.
Raghavan
,
M.
, and
Roth
,
B.
,
1995
, “
Solving Polynomial Systems for the Kinematic Analysis and Synthesis of Mechanisms and Robot Manipulators
,”
ASME J. Mech. Des.
,
117
(
B
), pp.
71
79
. 10.1115/1.2836473
10.
De Jalon
,
J. G.
, and
Bayo
,
E.
,
2012
,
Kinematic and Dynamic Simulation of Multibody Systems: the Real-Time Challenge
,
,
New York
.
11.
Nikravesh
,
P. E.
,
1988
,
Computer-aided Analysis of Mechanical Systems
, Vol.
186
,
Prentice-Hall
,
Englewood Cliffs, NJ
.
12.
Kreyszig
,
E.
,
2007
,
,
John Wiley & Sons
,
Hoboken, NJ
.
13.
Hernández
,
A.
, and
Petuya
,
V.
,
2004
, “
Position Analysis of Planar Mechanisms with R-Pairs Using a Geometrical–Iterative Method
,”
Mech. Mach. Theory
,
39
(
2
), pp.
133
152
. 10.1016/S0094-114X(03)00110-1
14.
,
P.
, and
Campbell
,
M. I.
,
2012
, “
An Automated Kinematic Analysis Tool for Computationally Synthesizing Planar Mechanisms
,”
ASME 2012 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
Chicago, IL
,
Aug. 12–15
,
American Society of Mechanical Engineers
, pp.
1553
1562
. https://doi.org/10.1115/DETC2012-70737
15.
Inc.
,
A.
,
2016
, “
Autodesk Inventor
,” http://www.autodesk.com/products/inventor/overview
16.
18.
Petuya
,
V.
,
Macho
,
E.
,
Altuzarra
,
O.
, and
Pinto
,
C.
,
2011
, “
Educational Software Tools for the Kinematic Analysis of Mechanisms
,”
Comp. Appl. Eng. Educ.
,
6
(
4
), pp.
261
266
.
19.
Furlong
,
T. J.
,
Vance
,
J. M.
, and
Larochelle
,
P. M.
,
1999
, “
Spherical Mechanism Synthesis in Virtual Reality
,”
ASME J. Mech. Des.
,
121
(
4
), pp.
515
520
. 10.1115/1.2829491
20.
Larochelle
,
P.
,
Dooley
,
J.
,
Murray
,
A.
, and
McCarthy
,
J. M.
,
1993
, “
Sphinx: Software for Synthesizing Spherical 4R Mechanisms
,”
Proceedings of the 1993 NSF Design and Manufacturing Systems Conference
,
Charlotte, NC
,
Jan. 6–8
, Vol.
1
, pp.
607
611
.
21.
Ruth
,
D.
, and
McCarthy
,
J.
,
1997
, “
Sphinxpc: An Implementation of Four Position Synthesis for Planar and Spherical 4r Linkages
,”
ASME Design Engineering Technical Conferences
,
Sacramento, CA
,
Sept. 14–17
, pp.
14
17
.
22.
Tse
,
D.
, and
Larochelle
,
P.
,
1999
, “
Osiris: a New Generation Spherical and Spatial Mechanism Cad Program
,”
Florida Conference on Recent Advancements in Robotics
,
Gainesville, FL
,
Apr. 28–30
.
23.
Purwar
,
A.
,
Deshpande
,
S.
, and
Ge
,
Q. J.
,
2017
, “
Motiongen: Interactive Design and Editing of Planar Four-Bar Motions Via a Unified Framework for Generating Pose- and Geometric-Constraints
,”
ASME J. Mech. Rob.
,
9
(
2
), p.
024504
. 10.1115/1.4035899
24.
Ge
,
Q. J.
,
Purwar
,
A.
,
Zhao
,
P.
, and
Deshpande
,
S.
,
2017
, “
A Task Driven Approach to Unified Synthesis of Planar Four-Bar Linkages Using Algebraic Fitting of a Pencil of G-manifolds
,”
ASME J. Comput. Inf. Sci. Eng.
,
17
(
3
), p.
031011
.
25.
Deshpande
,
S.
, and
Purwar
,
A.
,
2017
, “
A Task-Driven Approach to Optimal Synthesis of Planar Four-Bar Linkages for Extended Burmester Problem
,”
ASME J. Mech. Rob.
,
9
(
6
), p.
061005
. 10.1115/1.4037801
26.
Li
,
X.
,
Zhao
,
P.
,
Purwar
,
A.
, and
Ge
,
Q.
,
2018
, “
A Unified Approach to Exact and Approximate Motion Synthesis of Spherical Four-Bar Linkages Via Kinematic Mapping
,”
ASME J. Mech. Rob.
,
10
(
1
), p.
011003
. 10.1115/1.4038305
27.
Sharma
,
S.
,
Purwar
,
A.
, and
Ge
,
Q. J.
,
2019
, “
An Optimal Parametrization Scheme for Path Generation Using Fourier Descriptors for Four-Bar Mechanism Synthesis
,”
ASME J. Comput. Inf. Sci. Eng.
,
19
(
1
), p.
014501
. 10.1115/1.4041566
28.
Sharma
,
S.
,
Purwar
,
A.
, and
Ge
,
Q. J.
,
2019
, “
A Motion Synthesis Approach to Solving Alt-Burmester Problem by Exploiting Fourier Descriptor Relationship Between Path and Orientation Data
,”
ASME J. Mech. Rob.
,
11
(
1
), p.
011016
. 10.1115/1.4042054
29.
Li
,
X.
,
Ge
,
X.
,
Purwar
,
A.
, and
Ge
,
Q. J.
,
2015
, “
A Unified Algorithm for Analysis and Simulation of Planar Four-Bar Motions Defined with R- and P-Joints
,”
ASME J. Mech. Rob.
,
7
(
1
), p.
011014
. 10.1115/1.4029295
30.
Liao
,
Q.
, and
McCarthy
,
J. M.
,
1997
, “
On the Seven Position Synthesis of a 5-SS Platform Linkage
,”
ASME J. Mech. Des.
,
123
(
1
), pp.
74
79
. 10.1115/1.1330269