A Unified Real-Time Motion Generation Algorithm for Approximate Position Analysis of Planar N-Bar Mechanisms

The kinematic simulation of multi-body systems is crucial for the rapid design and evaluation of mechanisms in computer-aided design (CAD) systems. Faster simulation algorithms can significantly reduce the time required for mechanism design and enable the creation of high-quality datasets for machine learning research, which has garnered interest in recent years [1–9]. A recent position paper by Purwar and Chakraborty [10] emphasizes the importance of generating high-quality datasets for designing robot mechanisms that rely on robust and fast kinematic simulation. Planar linkage mechanisms are widely utilized in products due to their practicality, simplicity, and performance characteristics. While several commercially available CAD systems have general-purpose multi-body dynamics capabilities, the development of stand-alone kinematic simulation programs with comprehensive capabilities has been limited, with most software applications originating from academic research groups. Examples of such programs include lincages [11,12], kinsyn III [13], Kihonge et al. [14], Spades [15], Sphinx [16], Sphinxpc [17], Osiris [18], and Synthetica [19]. Although these projects are no longer active, a few notable commercially available systems include SAM [20], MechGen [21], Linkages [22], Ch Mechanism Toolkit [23], MechDesigner [24], and Universal Mechanism [25]. Geometer’s Sketchpad [26] and Geogebra [27], which are primarily geometry tools, have also been used for simulating mechanisms. Several other kinematic simulation tools, such as Linkage Mechanism Simulator [28], PMKS+ [29,30], GIM [31], Simionescu [32,33], and Schmidt and Lax [34], provide varying levels of capabilities and have filled a critical need for the mechanism simulation. This paper presents a novel simulation algorithm implemented in a web-based kinematic design and simulation application called motiongen,² which provides real-time simulation of N-bar planar linkage mechanisms with an arbitrary number of revolute and prismatic joints, rotary and linear actuators, and topological structures. Real-time simulation requires prompt feedback and data updates with low latency. While there are no rigid rules for the time to update, it is generally expected that this updating will happen quickly and typically at a rate that mimics real-world time, such as seconds or milliseconds per simulation time-step.

Kinematic simulation of planar linkages requires computing the unknown positions of moving joints, which leads to solving a set of geometric constraints—a system of equations constructed from the mechanism. Generally, there are two main issues in reducing the simulation time. The first is to create a fast computation method, or solver, to solve a specific set of geometric constraints. One of the simplest analytical methods for planar linkages is the dyadic decomposition [35,36], which solves one planar joint position, or two variables (x and y), at a time. While this works for simpler mechanisms, researchers have developed other methods for handling more complex structures. Dhingra et al. [37] employed symbolic computation and followed the computation of a Gröbner basis with a solution based on a Sylvester-type determinant. On the other hand, Wampler [38,39] proposed a combination of the Dixon determinant procedure from Nielsen and Roth [40] with a complex plane formulation. He presented a method for formulating kinematic equations in complex planes using isotropic coordinates, which results in loop equations in a simplified form that can be solved using the Sylvester-type formula. Hernández and Petuya [41] have introduced a novel geometric iterative (GI) technique for solving the position problem with only revolute joints, which does not require an initial guess, as in the case of local optimization methods. However, the Newton–Raphson (NR) method is still a popular choice for solving systems of nonlinear equations [42] due to its numerical search for a particular solution, with quadratic convergence in the neighborhood of the solution.

The second issue is the identification and joint position analysis of subassemblies within a mechanism using a divide-and-conquer approach. This technique involves identifying subsets within the geometric constraints of the system, which can be used to solve some of the unknown variables, or unknown joints. The subset can then be substituted back into the system to solve the remaining unknown variables. This approach enables the solution of a set of kinematic constraints to be broken down into a sequence of solutions of the kinematic constraint subsets, a process that typically takes less time than solving the entire set of constraints together. This process of identifying subsets is also referred to as the decomposition of the system. Researchers have proposed various methods for implementing the divide-and-conquer approach in system analysis. Ait-Aoudia et al. [43,44] demonstrated a decomposition method that utilizes bipartite graphs under-laid by systems of geometric constraint equations. Bouma et al. [45] and Fudos and Hoffmann [46] proposed a graph-theory-based method to reduce the number of variables that need to be solved in each step. Additionally, Fudos and Hoffmann developed a phase called the reduction of the system that simplifies the system of equations by merging fixed links together into a single link.

Most simulation software apply one or more of these aforementioned techniques. For example, linkages [22] uses the dyadic decomposition extensively, while pmks/pmks+ [29,30] and gim [31] use the GI technique. autocad based simulation [47] and mekin2d [32] use the modular approach to directly identify the subassemblies with a saved Assur group library.

Compared to these simulation packages, our method’s originality lies in reducing all geometric constraints to a unified form and converting all elements of a mechanism into revolute only joints and links, thereby allowing us to require only two types of solvers for all degrees-of-freedom (DOF) planar N-bar mechanisms. We also developed a method to analyze the complexity of a certain mechanism. Specifically,

1. The prismatic joints and actuators are converted to a combinations of revolute joints and links. This novel reinterpretation means that an additional solver for the prismatic joint is not necessary. Furthermore, this also means a simpler optimization function to program.

2. Our algorithm requires only two solvers. One provides the fastest computation speed possible using dyadic decomposition and the other is the optimization-based method for solving complex subassemblies.

3. This algorithm also allows arbitrary placement of actuators and arbitrary combinations of revolute joints and prismatic joints.

4. A mobility analysis method is presented to analyze the system of equations, which includes the reduction phase (static link merge) and the decomposition phase (dynamic and static link merge).

A flowchart of the complete simulation process can be seen in Fig. 1. The details of various steps in this figure are described in the following sections. Section 2 presents a matrix representation of the kinematic structures. Following that, Sec. 3 presents the reinterpretation trick and a derivation of the accuracy error incurred during the reinterpretation. Next, the general formulation of the system of equations is presented in Sec. 4. Then, Sec. 5 provides an overview of the two types of solvers. After that, the analysis of the kinematic system, an estimate of the time cost, and a measure of the structural complexity of mechanisms are presented in Sec. 6. Finally, we present a few typical example mechanisms in Sec. 7, which seek to illustrate the overall method.

In summary, the main contributions of this work are in (1) proposing a reinterpretation trick to reduce all geometric constraints to a single type, (2) topological modification of N-bar mechanisms to only include revolute joints, (3) formulation of a general-purpose solver which can automatically select the fastest method for computing joint positions, and (4) implementation in a cross-platform browser-based software to satisfy the needs of users with varied expertise.

The mechanisms in this paper are described by their topological configurations. To fully represent a planar linkage mechanism, four basic components are considered:

1. Joints: Joints in this work are considered as points of interest characterized by their positional variables, denoted as x_p and y_p whether or not there are two links connected at that joint. Each joint is represented by j_p, where p is the joint identifier. Coupler or tracer points, which are points only on one link, are also considered joints since they can potentially function as joints. They are not distinguished from joints shared by multiple links.

2. Links: A link is a rigid body consisting of one or more joints. The relative positions of the joints in a link remain constant throughout the mechanism’s operation. A link can be binary, containing two joints, or n-ary, containing n specified number of joints. Additionally, a binary link with joints i and j can be denoted as l_i,j.

3. Slots: Slots are assumed to be straight and finite constraints placed on binary links to represent linear constraints. A slot is defined by its two endpoint joints, which determine the slot axis.

4. Actuators: Actuators introduce changes to the mechanism by varying either the angle or the length. There are two types of actuators: rotary and linear. An actuator is defined by its type and three joints; see Fig. 2 for an illustration.

Throughout this paper, subscripts are used to represent matrix elements, with indexing starting from 0. For example, K_i,j refers to the element in the ith row and jth column of the matrix K. If a single number represents a link, it indicates the index of that link. For instance, l_n represents the (n + 1)th link.

It should be noted that the joints are order sensitive. For example, for the rotary actuator α_i,j,k shown in the first row of I in Eq. (8), the angle constraint is between vector $v→j,i$ and $v→j,k$⁠. Furthermore, we use the cross-product value to represent this constraint so that the angle is not flipped. For the linear actuator α_p,q,r in the last row of I, the direction is defined by vector $v→p,r$⁠, and the output distance is defined by the length of $v→p,q$⁠. Figure 2 follows this description.

Last, it should be pointed out that each actuator brings one dynamic link (a similar concept is the active pair in Ref. [49]). For the rotary actuator α_i,j,k, when the link lengths of l_i,j and l_j,k and the input angle are given, l_i,k is the dynamic link and can be computed according to the law of cosine. For the linear actuator α_p,q,r, l_p,q is the dynamic link, and the length is the input itself.

The idea of reinterpreting prismatic joints as revolute joints arises from the fact that PR and RP links can be considered as sufficiently large RR links with one of their joints located either proximate or distal from the workspace of the mechanism. In the case of a PR link with a fixed line, the RR link can be chosen to be large enough within a desired accuracy with a distal fixed pivot while for the RP link with a fixed pivot, the link length of the approximated RR link should be large enough, but this time with the fixed pivot proximate. By using homogeneous coordinates, a planar joint position is represented in a projective plane as (x, y, w) with its position in Euclidean space given by (x/w, y/w). By setting the homogenizing factor w to a small value, the joint is moved far away from the origin (0, 0) and effectively converted into a prismatic joint. In the kinematic simulation, we perform the inverse operation by converting the prismatic joint back into revolute joints and links. This allows solvers that can only handle revolute joints to solve systems that include prismatic joints.

For ease of implementation and intuitive sketching, we use joints in a slot to represent prismatic joints. Figure 4(a) shows an example mechanism with a slot, which restricts the slot joint to a fixed-line segment. Slots can be fixed or can be moving so as to represent RP or PR links.

Here, the slot joint j₂ performs a linear motion along the slot. A straight line can be reinterpreted as the arc of a circle whose radius is large enough such that the curvature is near zero. Consequently, the straight line is represented as the arc of a circle, and the intersections between the straight line and the arc occur at the slot’s edges (end joints). Furthermore, the center of the circle lies on the perpendicular bisector of the straight line.

The resultant reinterpretation can be seen in Fig. 4(b). First, the slot is removed, and a new far joint (j₅, or j_f) is added to the slotted link. Then, two constraining links (l_3,6 and l_2,6) and a constraining joint (j₆, or j_c) are added to constrain the motion along the straight line. These two binary links together can also be seen as an RRR robot arm so that when the arm is stretched, the slot joint j₂ goes near j₄, and when two binary links overlap, j₂ goes near j₃. During simulation, the slot is removed, and the joint previously in the slot now moves along an arc. The farther the joint j₅, the flatter the arc is. Therefore, the difference between the arc and the straight line can be controlled by changing the distance of j₅ from the two end joints j₃ and j₄.

With these two link lengths and the positions of j₂ and j₃ in the initial state, the initial position of j₆ can be computed through arc intersection. Of course, there are usually two solutions, but both work for the simulation, and choosing either is sufficient.

We utilize two types of solvers for kinematic analysis: (1) the analytical (arc intersection) and (2) the numerical-based (Newton’s method). These methods are employed to handle mechanisms with varying complexity and optimize computational efficiency.

The analytical method, specifically the arc intersection method, is highly time-efficient. However, it is not capable of solving mechanisms with complex loop structures. On the other hand, the numerical-based method, i.e., Newton’s method, can handle mechanisms with arbitrary complexity, but requires more computational time compared to the arc intersection method, especially when the number of joints increases.

Here, F(x_t) = f(x_t)^Tf(x_t), which is the sum of squared error (SSE). This value is close to zero when vector x is a good numerical approximation, which also serves as the termination of the iteration. Term $∇F(xt)$ is the gradient of F(x_t) and $[HF(xt)]$ is the Hessian matrix of F(x_t).

These two methods are both computationally expensive because they require matrix multiplications. Furthermore, the optimization methods take several iterations to converge. So, for the same number of joints to solve, they are much slower compared to the arc intersection solver. Through our experiments, we found that LM tended to converge faster in most scenarios.

According to the description above, to ensure a specific set of links having 0DOF, we must start the examination from the smallest subset of this link set. In the worst-case scenario, all subsets except the empty set and one-member sets need to be examined. This leads to a possibly unacceptable time cost, which is one of the reasons why some modern simulation software uses the Assur Group library. However, the Assur Group is infinite, and the library saves the most common kinematic assembly. As for our algorithm, we employed the aforementioned criterion to maintain generality.

In this section, we present a few examples to illustrate the working of our algorithm. The first example includes all steps of the solving procedure for a mechanism with a combination of binary, ternary, quaternary links, and a moving slot connected by revolute and prismatic joints. Then, this section will go through several cases as subsections to discuss each aspect of this simulation algorithm. Furthermore, in Sec. 7.6, we present the time complexity of each of the example mechanisms.

Figure 5 shows the example mechanism to be analyzed in this section. We note that this is purely an academic example chosen to illustrate several steps of position analysis presented earlier.

The length of the triad, l₈, l₉ and l₁₀ are computed according to Sec. 3 in the paper.

Static Link Merge: In the reduction phase, some of the links in the mechanism can be merged, resulting in less number of links. Furthermore, through this merge operation, two new binary links are found, which are l_0,3 and l_6,9. The subscripts of the new binary links represent the endpoint joint location of the links. Figure 7 shows the complete process.

Following the conversion, the next step is to find the least number of joints to solve per step, as shown in Fig. 8. The decomposition starts from the known joints j₀, j₄, and j₉, as shown in Fig. 8(a).

It is worth noting that a dynamic link does not replace its corresponding actuator. Instead, it is a necessary condition to satisfy the angle constraint. Although not shown in this example, a dynamic link can be useful sometimes to help decompose the system, i.e., giving more possible solving paths.

The next steps from Figs. 8(d)–8(h) use the same logic: if two joints are solved and one unknown joint connects to these two joints, this unknown joint can be solved through arc intersection. Table 1 shows the details.

Solving the positions of the joints is in the positional analysis stage, where the actual input θ_(0,4,5) is given. Iteratively solving all these steps with different input values gives the path of all joints as shown in Fig. 9.

With the reinterpretation of prismatic joints, our algorithm can simulate all combinations of RR, RP, PR, and PP links, whether fixed or moving; see Fig. 10. Notably, a true slider in a slot is represented by two joints. This is how a prismatic joint is drawn in motiongen. Comparatively, the three-joint relationship for a slot joint and two end joints is equivalent to an RP link.

Our algorithm supports moving slots and their combinations as well. Furthermore, a slot can be embedded in other slots. An example is shown in Fig. 11. For this mechanism, j₈ in the right figure shows that when the reinterpreted far joints (j₉, j₁₀, and j₁₁) are avoidably placed too close to the two end joints, the resulting path is curved and incurs large approximation error.

In this example, we present a complex kinematic structure, which requires that all the moving joints are solved simultaneously. Figure 12 shows such a mechanism, where there are 12 “legs” that connect with the ground and 12 serial links that each connect two legs one after another. In this case, all the distance constraints (links) and all the angle constraints (the actuator between l₆ and link l₂₅) need to be accounted for in the system of equations in order to solve for positions of movable joints. Additionally, removing any moving link will result in an under-constrained system. Therefore, this system cannot be further decomposed. Furthermore, this example also shows that an actuator can be placed on a moving joint and that an Assur group can be an infinite group because there can be n number of legs in a certain subassembly.

Pantographs are harmonic graphs that use electric motors and links to move a pen to create artistic drawings. Figure 13 is an example of such a machine. This mechanism is a rotating two-DOF five-bar mechanism with revolute joints only, where the path of j₅ is a long curve. This case is a simple system (only the arc intersection method is needed). However, it has many states (i.e., a big s in Eq. (30)). Therefore, the simulation time is comparatively long.

Table 2 presents the time complexity statistics for the mechanisms demonstrated in Secs. 7.2– 7.5. Again, for simplicity, we use C_a and C_n to represent the time constants for the arc intersection solver and Newton’s solver (regardless of which one), respectively.

From this table, it is clear that when a mechanism can be solved with the arc intersection solver only, the time to compute is short, usually less than 10 ms. The outlier is the Pantograph machine, which has a comparatively much larger number of states. If an optimization-based solver is used, the computation time is always more than 10 ms. Furthermore, its path is not long because the millipede structure is mechanically inefficient. However, it takes more than 200 ms for this small number of states.

We also tested the speed difference between the arc intersection solver and Newton’s solvers for the same question (system of two/three equations). For the pantograph machine in Fig. 13, the LM takes around 1750 ms to finish, and the Newton–Raphson solver takes around 1450 ms to complete, compared to 60 ms for the arc intersection solver.

In this paper, we have introduced a simulation algorithm that reinterprets prismatic joints in planar linkage mechanisms as a combination of links and revolute joints. We provided representations of the modified mechanisms and showed that there are only two methods needed to perform positional analysis of any planar mechanism. One of the methods, dyadic decomposition is computationally efficient, but can perform simulation only on relatively simpler mechanisms, while the optimization-based method works in all of the other cases. Both of the methods leverage a unified form of geometric and actuation constraints to find the unknown positions of moving joints, while the error incurred in approximating prismatic elements with revolute joints and links is user-controllable. Additionally, this paper presented expressions for estimating simulation time and structural complexity based on mobility analysis. Several examples demonstrated the algorithm’s efficiency and effectiveness in simulating mechanisms with various joint types and actuators. Moreover, this work highlights the conversion potential of Assur groups from mixed revolute joints and prismatic joints to exclusively revolute joints. Future extensions of this work include simulation of one and multi-degrees-of-freedom spatial and spherical mechanisms.

This work has been financially supported by the National Science Foundation under research grant STTR phase II #2126882 to co-author and Co-PI Purwar who also holds stocks in Mechanismic Inc. The research findings included in this publication may or may not necessarily relate to the interests of Mechanismic Inc. The terms of this arrangement have been reviewed and approved by Stony Brook University in accordance with its policy on objectivity in research.

All findings and results presented in this paper are those of the authors and do not represent those of the funding agencies.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.