In this work, a new method is introduced for solving large polynomial systems for the kinematic synthesis of linkages. The method is designed for solving systems with degrees beyond 100,000, which often are found to possess quantities of finite roots that are orders of magnitude smaller. Current root-finding methods for large polynomial systems discover both finite and infinite roots, although only finite roots have meaning for engineering purposes. Our method demonstrates how all infinite roots can be precluded in order to obtain substantial computational savings. Infinite roots are avoided by generating random linkage dimensions to construct startpoints and start systems for homotopy continuation paths. The method is benchmarked with a four-bar path synthesis problem.

## Introduction

The kinematic synthesis of linkages often requires finding the finite isolated roots of large polynomial systems. The size of a polynomial system is measured by its number of roots, of which its degree provides a maximum. It is interesting that the kinematic design of relatively simple mechanical systems involves formulating polynomial systems with degrees beyond 100,000. For example, to find all of the four-bar linkages which can pass a coupler point through nine arbitrary points in the plane involves solving a polynomial system with degree 286,720 [1]. Furthermore, more complex mechanical systems seem to result in exponential growth of the degree of their design equations [2,3].

This growth arises from additional design parameters that increase the number of dimensions in a nonlinear design space. Finding complete solution sets to the polynomial systems which comprise kinematic design equations is a major challenge which tends to demand large-scale computations. Many methods have been devised, which avoid computing complete solution sets [4], however, the advantage of obtaining complete sets is to provide a full survey of design options that might not be found through methods based on optimization or differential evolution. The state-of-the-art for solving large polynomial systems is polynomial homotopy continuation. Homotopy algorithms compute all the roots to general polynomial systems in time proportional to the degree of the system. These methods are quite powerful but still have practical limits based on available computing resources.

Although kinematic polynomial systems tend to have large degrees, they also tend to have sparse monomial structures that indicate the number of finite roots they possess is much smaller than their degree. Returning to the four-bar example, it has been shown in Ref. [1] that 8652 of 286,720 roots are finite and represent linkage designs. The polynomial degree totals the number of finite roots, infinite roots, and their multiplicities. For engineering purposes, only finite roots are sought. Observations of the monomial structure of a system can often be made in order to intelligently multihomogenize a system and remove many roots at infinity leading to a new *multihomogeneous root count*. For smaller systems, this strategy can eliminate all the roots at infinity, but for larger systems it is not always clear how the multihomogenization should proceed. Eliminating roots at infinity is important when applying homotopy algorithms because each root requires a path tracking computation. Counting the percentage of infinite roots for the four-bar example, it is estimated that 97% of the computational effort is dedicated to discovering roots at infinity, which are filtered out in the end.

In this work, we introduce a new method built on homotopy continuation termed *Finite Root Generation* (FRG) that obtains all or most of the finite roots of kinematic polynomial systems while avoiding all computations of roots at infinity. Traditionally, homotopy proceeds with a single start system and multiple startpoints, whereas FRG uses multiple start systems each with one startpoint, which are constructed from random linkage dimensions so that each startpoint is guaranteed to track to a finite root. Since startpoints are randomly generated, the possibility exists for tracking to a nonsingular finite root multiple times. This duplication rate is modeled by the *coupon collector problem* from probability theory and can be used to estimate the number of finite roots without actually computing all of them.

We apply FRG to a numerically general target system that is suitable for constructing parameter homotopies such that an FRG solution set only needs to be computed once for a particular family of polynomial systems. Parameter homotopies provide a means for efficiently solving all subsequent systems belonging to that family. Applying FRG to the benchmark four-bar path generator problem found over 95% of the roots in 26,658 homotopy paths and 99.988% of the roots in 87,549 paths. This is compared to 286,720 for the best known multihomogeneous homotopy [1] and 152,224 for the best known regeneration homotopy [5]. The benchmark slightly underperformed expectations calculated from modeling roots as equally probable coupons.

## Literature Review

Continuation methods were pioneered by Roth and Freudenstein [6], who developed the *bootstrap method* for solving the nonlinear path synthesis equations of a geared five-bar. They were motivated to provide an alternative to the Newton–Raphson method which required a good initial approximation in order to converge. In their work, the authors indeed constructed continuation startpoints from arbitrary starting mechanisms. Since then, the field of numerical algebraic geometry has flourished resulting in sophisticated homotopy continuation algorithms [5,7–10] that apply to general polynomial systems where startpoints are constructed by combinatoric procedures that are advantageously blind to the physical systems represented by the equations. In this work, we return to the schema of constructing startpoints by arbitrary mechanisms, but instead of finding one or two mechanism solutions, we aim to find complete root sets of numerically general synthesis systems. These general root sets are suitable for constructing parameter homotopies to efficiently solve later synthesis systems of the same family [11], a concept that appeared 26 years after Roth and Freudenstein [6].

Other solution techniques used to solve polynomial equations include resultant elimination methods, Gröbner bases, and interval analysis. Resultant elimination methods include modifications of Sylvester's dialytic elimination method [12] and involve procedures to transform root finding into a generalized eigenvalue problem. Su and McCarthy [13] showed how to use this method to solve a synthesis problem with 64 roots. Masouleh et al. [14] used Gröbner bases combined with resultants to compute 220 roots for the forward kinematics of a spatial parallel manipulator. Gröbner bases transform multivariate root finding into root finding for a sequence of univariate polynomials. Interval analysis takes advantage of interval arithmetic to find roots within a range of the solution space. Lee et al. [15] provided an example of interval analysis applied to the synthesis of an RRR spatial serial chain.

There have been modifications of homotopy methods to eliminate unwanted solutions. Tari et al. [16] appended additional equations to synthesis systems to eliminate degenerate solutions, but their method does not improve computation times. Tsai and Lu [17] devised a procedure for nine-point path synthesis of a four-bar where startpoints were constructed from five-point solutions without attempting to obtain the complete solution set. Perhaps today's most powerful method is called regeneration. Regeneration can greatly reduce the number of homotopy paths to track by solving a system “equation by equation” and discarding singular and infinite roots along the way [5]. However, regeneration still requires many roots at infinity to be tracked. Plecnik and McCarthy [2] used regeneration to approximate a solution set of 1.5 × 10^{6} roots.

In this work, we introduce a method which constructs startpoints and start systems from randomly generated mechanisms in order to avoid tracking homotopy paths to infinite roots. We benchmark FRG on the four-bar path synthesis problem and find that the rate at which finite roots are collected is predicted by the coupon collector problem of probability theory, providing a means for estimating the size of the root set being discovered. Finally, we use the parameter homotopy method to apply the obtained root set to an example linkage design problem.

## Description of Finite Root Generation

The FRG method generates startpoints and start systems that track to the finite roots of a target system using polynomial homotopy continuation. In order to ensure this is the case, a start system must possess a monomial structure that is no more or no less general than the target system, and a startpoint must be a finite root to that start system. These features are ensured for a single startpoint/start system pair by constructing the startpoint from a randomly generated linkage and constructing the start system from the motion of that linkage. A parameter homotopy then tracks that startpoint to a single finite root of the target system.

In order to find another finite root to the target system, a new random linkage can be generated so that another startpoint/start system pair is constructed and tracked to (hopefully) another finite root of the target system. If there are *N* finite roots of the target system to be discovered, and there is equal probability of happening upon any root, then the probability of the second iteration happening upon a different finite root from the first is (*N* − 1)/*N*, which for large *N* are good odds. If the root of the second iteration duplicates the first, then it is discarded, or if it is unique, then it is stored along with the root from the first iteration. The goal is to obtain all or most of the finite roots of the target system. Iterations (also called trials) continue in this manner as illustrated in Fig. 1 until a sufficient number of roots have been collected. As more and more roots are collected, the possibility of finding another unique root decreases, indicating a situation of diminishing returns. By tracking the rate of success of finding new roots between iterations, the size of the total root set can be estimated.

### Root Collection.

The decreasing probability of finding new roots is modeled by the coupon collector problem from probability theory [18]. That is, given a set of *N* unique coupons, how many trials *T _{N}* of picking random single coupons are expected in order to obtain all the coupons. Coupons are sampled with replacement and we assume all the coupons have equal probability of being picked. The problem statement is slightly modified to ask how many trials

*T*are expected to obtain

_{n}*n*coupons where

*n*≤

*N*. In our case,

*T*is the number of homotopy paths and

_{n}*n*is the number of finite roots we wish to obtain.

*n*− 1 roots have been collected, the probability

*p*of collecting the

_{n}*n*th root is

*X*denote the expected number of trials to take place in the interval between finding

_{n}*n*− 1 roots and

*n*roots. Because

*X*is a random variable with a geometric distribution, it is computed as

_{n}*T*is computed as

_{n}Equation (3) provides an estimation of the expected number of trials to collect *n* roots when the total number *N* is known.

### Estimation Procedure.

*N*is not known but can be estimated by the ratio of success

*α*of finding new roots defined as

*n*is the number of roots collected, and

*T*is the total number of trials it took to collect those roots. As the probability of finding new roots decreases,

_{n}*α*tends to decrease. The aim is to derive an expression for estimating the percentage of roots obtained by the success ratio

*α*. First, replace

*T*in Eq. (8) with its approximation from Eq. (7). Then, substitute $n\u0302=n/N$ as the percentage of finite roots found. The result is

_{n}*N*has been eliminated from the calculation. Inversion of Eq. (9) requires the principle branch of the Lambert function

*W*and takes the form

Equation (10) provides an estimation of the percentage of roots obtained $n\u0302$ from the current success ratio *α* during the FRG solution process. Note that this equation only provides an estimate of the number of finite roots based on the success ratio, unlike Bézout numbers which provide conclusive statements on the upper bound of finite roots. FRG estimates might indicate smaller root sets, but without certainty.

## Benchmark

In order to benchmark the FRG method, it was applied to the nine-point path synthesis problem for four-bar linkages. The complete solution to this problem involves finding 8652 finite roots but the smallest known formulation of this system counts 286,720 roots (including infinite roots). Multihomogeneous homotopy has been applied to solve this problem by tracking 286,720 homotopy paths to discover the 8652 finite roots [1] (see Fig. 2). As an improvement, regeneration homotopy finds all the finite roots while only tracking 152,224 paths [5]. From Eq. (3), it is expected for FRG to find 95% of the roots in 25,922 paths and 100% in 83,430 paths (see Figs. 3 and 4). In other words, the expectation is for 31% of the computational effort to find 95% of the roots and 70% of the effort to find the final 5% of the roots. In this section, we formulate the synthesis equations and test the expectations.

### Formulation of Four-Bar Path Synthesis.

The objective of path generation is to find a four-bar linkage that can move a trace point attached to its coupler link through *m* points *P _{j}*,

*j*= 0,…,

*m*− 1. As shown in Fig. 5, the locations of the four-bar's ground pivots are represented by

*A*=

*A*+

_{x}*A*and

_{y}i*B*=

*B*+

_{x}*B*. The locations of the moving pivots in position

_{y}i*j*= 0 are represented by

*C*=

*C*+

_{x}*C*and

_{y}i*D*=

*D*+

_{x}*D*. The rotations of links

_{y}i*AC*,

*BD*, and

*CDP*from position 0 to

*j*are measured by $\varphi j$,

*ψ*, and

_{j}*θ*, respectively. Using the exponential rotation operators

_{j}*j*

*Q*are eliminated by solving for them in Eq. (12) and substituting into Eq. (14). Similarly,

_{j}*S*is eliminated by solving Eq. (13) and substituting into Eq. (15). These substitutions obtain

_{j}Equation (19) contains the eight unknowns {*A*, *B*, *C*, *D*, $A\xaf,\u2009B\xaf,\u2009C\xaf,\u2009D\xaf$} and represents a square system of eight equations for the case *m* = 9.

Our formulation is similar to those presented in Refs. [1] and [5], which provide additional information on the maximum number of roots in order to determine the number of homotopy paths to be tracked. Advantageously, FRG paths are tracked without the knowledge of root ceilings, similar to the regeneration method. For the sake of comparison, sorting the variables of $T$ into groups ${A,A\xaf},\u2009{B,B\xaf},\u2009{C,C\xaf}$, and ${D,D\xaf}$ finds a four-homogeneous Bézout number of 645,120, which was the same calculated in Ref. [5]. In Ref. [1], intermediate variables were introduced to further reduce the multihomogeneous Bézout number to 286,720.

### Methods.

*i*, 1 +

*i*), assigning these numbers to

*P*and $P\xafj,\u2009j=0,\u2026,8$,

_{j}and then substituting these parameters into $T$. To obtain a numerically general system, conjugate relationships were broken.

*i*, 0.5 + 0.5

*i*). Four of these numbers define a startpoint

and the fifth is assigned to *P*_{0}. Here, conjugate relationships remain intact so that these five numbers define a four-bar mechanism. Eight configurations are selected from this mechanism with disregard to whether they belong to the same circuit. These eight configurations define points $(Pj,P\xafj),\u2009j=1,\u2026,8$, which along with $(P0,P\xaf0)$ are substituted into $T$ to obtain a start system that corresponds to *s _{p}*. Random distributions were uniform.

Path tracking was completed using the algorithm within the software bertini [20]. This path tracking algorithm utilizes random numbers to form projective patches and defines homotopies based on a random parameter *γ* [5] (unrelated to the constant of Eq. (6)). Occasionally, bertini would experience a path tracking failure, e.g., reaching a minimum step size. In these cases, an FRG iteration would be considered invalid. Table 1 shows a few startpoints and start system parameters used to find roots of the target system defined by the parameters of Eq. (21).

## Results

The FRG method was applied to find all 8652 roots of a numerically general version of $T$. The algorithm was executed for 100,000 trials. By the 87,549th trial, the algorithm found 8651 roots, but did not find the final root beyond that. The expectation was to find all the roots on the 84,430th trial within a standard deviation of 11,092 trials. Curtailing results before the brunt of diminishing returns, the algorithm found 95% of the roots by 26,658 trials, underperforming the expected value by 736 trials. Similar to expectations, FRG found over 95% of the roots in 30% of the iterations required to find 99.988% of the roots. This is depicted in Fig. 6, which shows the number of roots obtained as trials progressed. Roots accumulated quickly in the beginning leading to a diminishing rate of returns such that at least 70% of the calculations were dedicated for obtaining the final 5% of the roots. Table 2 gives the number of trials performed versus expected to acquire various percentages of the complete solution set.

As trials progressed, an estimate of the percentage of roots obtained was calculated from the success ratio using Eq. (10). This is a useful feature for systems in which the total number of roots is unknown. These estimates are compared to the actual percentages for the known benchmark in Table 2, and the running difference between the two is graphed in Fig. 7. Estimation error was greatest at the beginning of FRG, however, as the number of sample trials that comprise the calculation of the success ratio increased, the estimation became more accurate. Predictions of the total number of roots are given in Table 2. The success ratio as a function of the percentage of roots obtained for the benchmark is shown in Fig. 8 compared to the theoretical function from Eq. (9). This function illustrates the sharp decrease in success near the end of FRG, tending toward zero for systems with more roots.

The error estimates of this paper are based on the assumption that roots appear with equal probability. Conducting Pearson's chi-squared test on tabulated root appearances revealed that this was not the case. Nonetheless, Figs. 6–8 and Table 2 demonstrate reasonable accuracy of the estimates. The path failure rate was 5% for the benchmark.

## Example

The motivation behind FRG is to obtain solution sets to large polynomial systems for kinematic design. In this paper, the method is introduced and only applied to four-bar path synthesis for benchmarking. In order to bridge FRG with tangible kinematic design, this section presents an example. The objective is to design a linkage that traces a D-shaped path potentially for a walking mechanism.

Note that the straight portion of the “D” was curved (Figs. 9 and 10) in the hope that this would result in four bars with smaller link lengths. Substituting Eq. (23) into $T$ forms a synthesis system which was solved by a parameter homotopy constructed from the 8651 roots obtained in Sec. 5. The results include 554 roots that represent physical link geometry of which 16 were found to be free of branch and circuit defects. Note that a finite root is necessary but not sufficient for a solution to represent physical geometry. In other words, the number of finite roots provides a bound for useful engineering solutions.

The number of linkage solutions found should be divisible by six because of known symmetries in the synthesis system. This was not the case due to a path failure rate of 7% during parameter homotopy. Solutions should exist in symmetric pairs with values $A\u2194B$ and $C\u2194D$ swapped, as well as in sets of three cognates [21], creating a six-way symmetry. With this in mind, the benchmark problem may be reduced to collecting *N* = 1442 sets of roots.

Examples of defect-free and defective linkages are given in Figs. 9 and 10. The linkage shown in Fig. 9(a) has the potential to be used as a walking mechanism although its large link lengths are not ideal. The linkages shown in Fig. 10 are similar to Chebyshev cognate linkages, which are known to trace a similar “D”. A third Chebyshev-like linkage was found as well, which is the horizontal mirror of Fig. 10(a). The linkage shown in Fig. 10(b) is branch-defective meaning it cannot be actuated from a base joint. Linkage designs were analyzed for branch and circuit defects by computing mechanism singularities and ranges of motion. This analysis did not consider order defects. The lack of strong design candidates from the example's results motivates a search into more complex linkage topologies, which is out of the scope of this paper.

## Conclusion

The FRG technique was introduced for solving large polynomial systems for the kinematic design of linkages. The technique provides a means of generating homotopy continuation startpoints that track to the finite roots of a numerically general target system, precluding infinite roots which often consume the large majority of computations. The rate at which FRG accumulates unique roots is modeled as a coupon collector problem which also provides means for estimating the entire root set. Applying FRG to a benchmark four-bar design problem found all but one root in 42% less paths than the best known regeneration method. Characteristic of FRG is a diminishing rate of returns for finding the final roots. FRG found 95% of roots in 82% less paths than the best known regeneration method. The results of the algorithm were applied to design a few four bars.

## Acknowledgment

The authors gratefully acknowledge the support of the National Science Foundation Award Nos. CMMI-1549667 and CMMI-1636302. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.