An Improved Fourier-Based Method for Path Generation of Planar Four-Bar Linkages without Prescribed Timing

Abstract

Four-bar linkages are critical fundamental elements of many mechanical systems, and their design synthesis is often mathematically complicated with iterative numerical solutions. Analytical methods based on Fourier coefficients can circumvent these difficulties but have issues with time parameters assignment for path generation without prescribed time in previous studies. In this paper, an improved Fourier-based point-to-point combination method is presented, which generates more points by Fourier approximation and assigns the time parameters to the given points while allowing discarding solutions with order defects. This method relies on the Fourier coefficients, improving the accuracy of synthesis solutions, and simplifying the computational procedure. Time parameters are assigned directly to the given points, which avoids the complex calculations to find intersection points in the given path, eliminates combinations that would lead to solutions with order defects, and simplifies the assessment process of synthesis results. The parameters obtained by the point-to-point combination method can be used as the parameters of the input dyad, skipping the first set of design equations for faster calculation. Several examples are presented to demonstrate the advantages of the proposed synthesis method, which is easy-understanding, computationally efficient, and yields more accurate solutions than available synthesis methods.

1. Introduction

Dimensional synthesis, which determines the dimensions of the mechanism to satisfy specific motion characteristics, plays an important role in engineering design. It can be classified into three categories: function generation, path generation, and motion generation, depending on the goal of the mechanism. Planar four-bar linkage is one of the simplest mechanisms, and its synthesis is thus a vital kinematics problem to address.

Path generation ensures that a point on the coupler link can move along a desired path. Solving path synthesis problems requires minimizing the differences between the generated and desired path [1]. Three types of approaches are usually employed for this problem: graphical, analytical, and optimization methods. Graphical methods [2] are simple and intuitive but lack accuracy. Analytical methods [3,4,5] are accurate and easy to understand, but they are only applicable to paths with a few numbers of precision points, and their complexity increases significantly with the number of precision points. Optimization methods [6,7,8,9] overcome the limitation of analytical methods with accurate solutions but suffer from slow optimization, difficulty in determining suitable initial guesses, and possibly resulting in non-convergence of the iterative solutions. Researchers are continuously making efforts to improve these methods.

Numerical atlas methods, as developed by graphical methods with advances in computer graphics, compile and collect a large number of coupler curves and generate linkages by comparing the similarity of the curves [10]. In building up a numerical atlas database, different methods are used to extract the feature information of the coupler curves [10], including the coupler-angle function curve [11], Fourier descriptors [12,13,14,15,16,17], wavelet transforms [10,18], and cumulative angular function or curvature [19,20]. Analytical methods are mathematically complicated for the equations to be solved. For path synthesis of planar four-bar linkages based on loop closure equations, the number of precision points is nine, which is limited to matching a number of design variables [21]. The equations can be solved in closed form without iteration for up to 5 precision points and using iterative methods for 6 to 9 precision points [22]. Another analytical method, the coupler curve equation method that traces a continuous path using an algebraic coupler curve equation proposed by Blechschmidt and Uicker [23] and extended by Ananthasuresh and Kota [24], was once considered as over-determined without analytical solutions [25]. However, Bai and Angeles [26] demonstrated the existence of exact solutions for the equations. This research was followed up by Bai [27] and Wu et al. [21]. Li et al. [28,29] introduced Fourier coefficients to the analytical method for a novel approach for function, path, and motion generation of four-bar linkages. This method obtains approximate solutions based on loop closure equations. Optimization methods are used by most of the literature for solving path generation problems [30]. Research relating to optimization methods can be classified into two types [31]: the former is research that focuses on the description of the optimization problem [1,31,32,33,34,35,36], and the latter concentrates on modifications to optimization techniques [9,37,38]. Commonly used optimization methods are traditional gradient-based [35,39] and metaheuristic methods [31,37,40,41,42,43]. Among all these methods, the Fourier-based analytical method overcomes the limitation of precision point numbers and reduces computation time by avoiding building databases and iterative calculations. The method is easy to understand and ready for programming. However, there remains scope for further improvement of this method.

In the works of Yu [11] and Li et al. [29], a two-DOF mechanism is used to establish the relationships between the given points and the input angles (i.e., the time parameters). However, this method requires complex calculations to identify intersection points along the given path with the floating link of the auxiliary mechanism. Furthermore, it may result in solutions with order defects. Assessing the accuracy of the generated results can also be challenging, as it requires finding the closest coordinates to the given points in the generated path. The objective of this paper is to improve the time parameter assignment process. Firstly, Fourier approximation is used instead of spline interpolation to better capture the curve of the given points, improving the accuracy of the synthesis, and simplifying the computational procedure at the programming level. Secondly, time parameters are assigned directly to the given points to avoid complex calculations. Decision conditions are proposed to find valid sets of point-to-point combinations, ensuring discarding combinations that lead to solutions with order defects and simplifying the accuracy evaluation process without seeking the corresponding generated path point for the given points. Additionally, the parameters obtained by the point-to-point combination method are used as the parameters of the input dyad to skip the first set of design equations for faster calculation.

This paper is organized as follows. Section 2 provides background concepts and the loop closure equations, from which the design equations are derived for the planar four-bar linkages. Section 3 presents the numerical procedure to obtain the Fourier coefficients for three types of paths, i.e., paths with uniform spacing time, non-uniform spacing time, and without prescribed time. An improved method to obtain the Fourier coefficients for path generation without prescribed time is proposed. Section 4 presents the synthesis procedure and criteria to assess the accuracy of the synthesis results. Three numerical examples are presented in Section 5 to demonstrate the advantages of the proposed methodology. The conclusions of this work are summarized in Section 6.

2. Planar Four-Bar Linkage Kinematic Formulation

2.1. Fourier Representation for Path Generation

A planar four-bar linkage for path generation is illustrated in Figure 1, where $P$ is a coupler point fixed in the coupler link BC. The coupler link dimensions PB and PC are denoted as $f$ and $g$, and the angle between PB and BC is $\alpha$. The length of each link is denoted as $a$, $b$, $c$, and $d$ for the input link AB, the coupler link BC, the output link CD, and the fixed link AD, respectively. $r_A$, $\mu ,$ and $\beta$ describe the position of the fixed link AD relative to the coordinate origin. The input rotation angle is $\phi$ with ${\phi }_0$ as the initial input angle. $\theta$ and $\psi$ are the angles links. Path generation is defined to ensure that the coupler point $P$ can move along a desired curve. Assume that the input link AB rotates at a constant speed $\omega$. The curve traced by $P$ can be represented in rectangular notation as:

$$r_p\left(t\right)=x\left(t\right)+iy\left(t\right)$$

(1)

The exponential form of the Fourier series of $r_p\left(t\right)$ is given by:

$$r_p\left(t\right)={\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{c}_n{e}^{in\omega t}={\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{c}_n{e}^{in\phi }$$

(2)

where the Fourier coefficients $c_n$ are:

$$c_n=\frac{1}{T}{\int }_0^T\left(x\left(t\right)+iy\left(t\right)\right)e^{-in\phi }dt\left(n=0,\pm 1,\pm 2,\pm 3\dots \right)$$

(3)

where $c_0$ is the coefficient of the fundamental harmonic, $c_n$ and $c_{-n}$ are the Fourier coefficients of the $n$-th harmonic.

As shown in Figure 1, The coupler curve ${\mathit{r}}_{\mathit{p}}$ can be expressed using phasor notation as:

$${\mathit{r}}_{\mathit{p}}=r_A{e}^{i\mu }+a{e}^{i\left({\phi }_0+\phi \right)}+f{e}^{i\left(\alpha +\beta \right)}{e}^{i\theta }$$

(4)

where $e^{i\theta }$ relates the relative motion between the input angle $\phi$ and coupler angle $\theta$ with a period $T=2\pi /\omega$, the exponential form of the Fourier series of which is given by:

$$e^{i\theta }={\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{c}_n^{\prime }{e}^{in\omega t}={\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{c}_n^{\prime }{e}^{in\phi }$$

(5)

Substituting Equations (2) and (5) into Equation (4) and expanding it according to $e^{in\phi }$ yields:

$$\begin{gathered} {\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{c}_n{e}^{in\omega t}=r_A{e}^{i\mu }+c_0^{\prime }f{e}^{i\left(\alpha +\beta \right)}+(a{e}^{i{\phi }_0}+c_1^{\prime }f{e}^{i\left(\alpha +\beta \right)})e^{i\phi }+ \\ {\sum }_{n\ne \mathrm{0,1}}{c}_n^{\prime }{f{e}^{i\left(\alpha +\beta \right)}e}^{in\phi } \end{gathered}$$

(6)

This gives the relationships between $c_n$ and $c_n^{\prime }$ as follows:

$$\begin{array}{c}{c}_0^{\prime }=\frac{c_0-r_A{e}^{i\mu }}{f{e}^{i\left(\alpha +\beta \right)}}=\frac{c_0-r_A{e}^{i\mu }}{f}{e}^{-i\left(\alpha +\beta \right)}\\ c_1^{\prime }=\frac{c_1-a{e}^{i{\phi }_0}}{f{e}^{i\left(\alpha +\beta \right)}}=\frac{c_1-a{e}^{i{\phi }_0}}{f}{e}^{-i\left(\alpha +\beta \right)}\\ c_n^{\prime }\left(n\ne \mathrm{0,1}\right)=\frac{c_n}{f{e}^{i\left(\alpha +\beta \right)}}=\frac{c_n}{f}{e}^{-i\left(\alpha +\beta \right)}\end{array}$$

(7)

2.2. Design Equations for the Input Dyad of Path Generation

A planar four-bar linkage for path generation is made of two dyads, of which the input dyad is illustrated in Figure 2. The input dyad involves five design variables $\left(r_A, \mu , a, {\phi }_0, f\right)$. Design equations to solve these variables can be obtained based on the loop closure equation of this dyad.

As shown in Figure 2, the vector loop closure equation can be written as:

$${\mathit{r}}_{\mathit{A}}+\mathit{a}+\mathit{f}={\mathit{r}}_{\mathit{p}}$$

(8)

Using phasor notation and rearranging, Equation (8) can be presented as:

$${\mathit{r}}_{\mathit{p}}-r_A{e}^{i\mu }-a{e}^{i\left({\phi }_0+\phi \right)}=f{e}^{i\left(\alpha +\beta +\theta \right)}$$

(9)

The variable $\alpha$, $\beta ,$ and $\theta$ can be eliminated by multiplying Equation (9) by its conjugate equation to yield:

$${\mathit{r}}_{\mathit{p}}\overline{{\mathit{r}}_{\mathit{p}}}-{\mathit{r}}_{\mathit{p}}\overline{x}-{\mathit{r}}_{\mathit{p}}\overline{y}{e}^{-i\phi }-\overline{{\mathit{r}}_{\mathit{p}}}x+x\overline{x}+x\overline{y}{e}^{-i\phi }-\overline{{\mathit{r}}_{\mathit{p}}}y{e}^{i\phi }+\overline{x}y{e}^{i\phi }+y\overline{y}-f^2=0$$

(10)

where $x=r_A{e}^{i\mu }$, $y=a{e}^{i{\phi }_0}$. By introducing the Fourier series with harmonics ${\mathit{r}}_{\mathit{p}}={\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{c}_n{e}^{in\phi }$ and $\overline{{\mathit{r}}_{\mathit{p}}}={\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}\overline{c_n}{e}^{-in\phi }$, the loop closure equation in the form of complex numbers can be represented in the following form:

$${\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{H}_n{e}^{in\phi }=0$$

(11)

where:

$$\begin{array}{c}{H}_{-1}=\overline{x}{c}_{-1}+\overline{y}{c}_0+x\overline{c_1}+y\overline{c_2}-x\overline{y}-k_{-1}\\ H_0=\overline{x}{c}_0+\overline{y}{c}_1+x\overline{c_0}+y\overline{c_1}+f^2-x\overline{x}-y\overline{y}-k_0\\ H_1=\overline{x}{c}_1+\overline{y}{c}_2+x\overline{c_{-1}}+y\overline{c_0}-\overline{x}y-k_1\\ H_n=\overline{x}{c}_n+\overline{y}{c}_{n+1}+x\overline{c_{-n}}+y\overline{c_{-n+1}}-k_n \left(n=0, \pm 1\right)\end{array}$$

and:

$$k_n=\sum _{m=0}^{+\mathrm{\infty }}{c}_m\overline{c_{m-n}}$$

Applying the orthogonality condition of the harmonics, i.e., Equation (11) is satisfied if $H_n=0$, and considering $n=0, \pm 1, \pm 2$ harmonics, the following five design equations are obtained to solve for $(r_A, \mu , a, {\phi }_0, f)$:

$$\begin{array}{c}{H}_{-2}=\overline{x}{c}_{-2}+\overline{y}{c}_{-1}+x\overline{c_2}+y\overline{c_3}-k_{-2}=0\\ H_{-1}=\overline{x}{c}_{-1}+\overline{y}{c}_0+x\overline{c_1}+y\overline{c_2}-x\overline{y}-k_{-1}=0\\ H_0=\overline{x}{c}_0+\overline{y}{c}_1+x\overline{c_0}+y\overline{c_1}+f^2-x\overline{x}-y\overline{y}-k_0=0\\ H_1=\overline{x}{c}_1+\overline{y}{c}_2+x\overline{c_{-1}}+y\overline{c_0}-\overline{x}y-k_1=0\\ H_2=\overline{x}{c}_2+\overline{y}{c}_3+x\overline{c_{-2}}+y\overline{c_{-1}}-k_2=0\end{array}$$

(12)

2.3. Design Equations for Variables $\left(b,c,d,\alpha ,\beta \right)$

As shown in Figure 1, the vector loop closure equation of the four-bar linkage in the Cartesian coordinate system $x\prime o\prime y\prime$ can be expressed as:

$$\mathit{a}+\mathit{b}=\mathit{d}+\mathit{c}$$

(13)

Using phasor notation and rearranging, Equation (13) can be represented as:

$$e^{i\left({\phi }_0+\phi -\beta \right)}+b{e}^{i\theta }-d=c{e}^{i\psi }$$

(14)

The variable $\psi$ can be eliminated by multiplying Equation (14) by its conjugate equation to yield:

$$h_{-3}{e}^{-i\theta }+h_{-2}{e}^{-i\theta }{e}^{i\phi }+h_{-1}{e}^{-i\phi }+h_0+h_1{e}^{i\phi }+h_2{e}^{i\theta }{e}^{-i\phi }+h_3{e}^{i\theta }=0$$

(15)

where $h_{-3}=-bd$, $h_{-2}=ab{e}^{i\left({\phi }_0-\beta \right)}$, $h_{-1}=-ad{e}^{-i\left({\phi }_0-\beta \right)}$, $h_0=a^2+b^2+d^2-c^2$, $h_1=-ad{e}^{i\left({\phi }_0-\beta \right)}$, $h_2=ab{e}^{-i\left({\phi }_0-\beta \right)}$, and $h_3=-bd$. By introducing the Fourier series with harmonics $e^{i\theta }={\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{c}_n^{\prime }{e}^{in\omega t}$ and $e^{-i\theta }\&={\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}\overline{c_n^{\prime }}{e}^{-in\omega t}$, the loop closure equation in the form of complex numbers can be represented in the following form:

$${\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{H}_n^{\prime }{e}^{in\phi }=0$$

(16)

where:

$$\begin{array}{c}{H}_n^{\prime }=h_{-3}\overline{c_{-k}^{\prime }}+h_{-2}\overline{c_{-k+1}^{\prime }}+h_2{c}_{k+1}^{\prime }+h_3{c}_k^{\prime } \left(n\ne 0,\pm 1\right)\\ H_n^{\prime }=h_{-3}\overline{c_{-k}^{\prime }}+h_{-2}\overline{c_{-k+1}^{\prime }}+h_2{c}_{k+1}^{\prime }+h_3{c}_k^{\prime }+h_n \left(n=0,\pm 1\right)\end{array}$$

Applying the orthogonality condition of the harmonics, i.e., Equation (16) is satisfied if $H_n=0$, and considering $n=0, \pm 1, \pm 2$ harmonics, five design equations to solve $(b, c, d, \alpha , \beta )$ can be obtained by introducing Equation (7) as follows:

$$\begin{array}{c}-uv\overline{c_2}+uy\overline{c_3}+\overline{u}\overline{y}{c}_{-1}-\overline{u}\overline{v}{c}_{-2}=0\\ -uv\left(\overline{c_1}-\overline{y}\right)+uy\overline{c_2}+\overline{u}\overline{y}\left(c_0-x\right)-\overline{u}\overline{v}{c}_{-1}-v\overline{y}f=0\\ -uv(\overline{c_0}-\overline{x})+uy(\overline{c_1}-\overline{y})+\overline{u}\overline{y}(c_1-y)-\overline{u}\overline{v}(c_0-x)+\left(y\overline{y}+u\overline{u}+v\overline{v}-c^2\right)f=0\\ -uv\overline{c_{-1}}+u\left(\overline{c_0}-\overline{x}\right)+\overline{u}\overline{y}{c}_2-\overline{u}\overline{v}\left(c_1-y\right)-\overline{v}yf=0\\ -uv\overline{c_{-2}}+uy\overline{c_{-1}}+\overline{u}\overline{y}{c}_3-\overline{u}\overline{v}{c}_2=0\end{array}$$

(17)

where $u=b{e}^{i\alpha }$, $v=d{e}^{i\beta }$.

3. Numerical Procedure to Obtain the Fourier Coefficients

A numerical procedure must be applied to evaluate the Fourier coefficients, $c_n$, in Equations (12) and (17). The desired path for path generation is usually given as a sequence of points, with or without prescribed timing. In cases where the timing is not prescribed, time parameters should be assigned to the given points along the path to obtain the Fourier coefficients.

3.1. For Path Generation with Prescribed Timing

There are two types of paths with prescribed time: those with uniform spacing timing and those with non-uniform spacing timing. The Fourier coefficients can be directly calculated using Equation (3) for paths with uniform spacing timing. However, for paths with non-uniform spacing timing, a method based on the complex least-square method is adopted [44]. Assuming the number of the given paths is $M$ and the order of the Fourier coefficients to be obtained is $N$. The method is as follows:

$$\mathrm{C}={\left(\overline{X}X\right)}^T{\overline{X}}^{T}Y$$

(18)

where $\mathrm{C}={\left[{\mathrm{c}}_{\mathrm{o}}, c_1, c_{-1}, \dots , c_N, c_{-N}\right]}^T$ is a $(2N+1)\times 1$ column matrix of the Fourier coefficients of the given path, $Y={\left[x_1+i{y}_1, x_2+i{y}_2, \dots , x_M+i{y}_M\right]}^T$ is a $M\times 1$ column matrix of the given path that represented in rectangular notation, and $X$ is a $(2N+1)\times M$ matrix denoted as follows:

$$X=\left[\begin{array}{cccccc}1& e^{i{\phi }_1}& e^{-i{\phi }_1}& \dots & e^{i{N\phi }_1}& e^{-i{N\phi }_1}\\ 1& e^{i{\phi }_2}& e^{-i{\phi }_2}& \dots & e^{i{N\phi }_2}& e^{-iN{\phi }_2}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& e^{i{\phi }_M}& e^{-i{\phi }_M}& \dots & e^{iN{\phi }_M}& e^{-iN{\phi }_M}\end{array}\right]$$

The Fourier coefficients can be obtained by Equation (18) for any given path of $M$ points with non-uniform spacing timing.

3.2. For Path Generation without Prescribed Timing

As was found by McGarva and Mullineux [12], different time parametrization can result in different Fourier coefficients, it is crucial to select appropriate time parameters for path generation synthesis. A two-DOF auxiliary mechanism, proposed by Yu. et al. [11] to transform the given path to a coupler-angle curve, is used and improved to assign time parameters to the given points.

An auxiliary mechanism, ABP, is illustrated in Figure 3. The input link AB acts as a crank, and the input angle is denoted as $\phi$. The floating link BP moves accordingly when link AB rotates, with point P following the given path, and point P completes a closed trajectory as the input link rotates a whole circuit around point A. Crank AB is colinear with the floating link BP twice, extended or overlapping, during the rotating. The furthest and the closest distance between point A and the given path are denoted as $l_{max}$ and $l_{min}$ for $l_{A{P}_i}$ and $l_{A{P}_j}$, satisfying the following conditions:

$$\begin{array}{c}{l}_{BP}+l_{AB}=l_{max}\\ l_{BP}-l_{AP}=l_{min}\end{array}$$

(19)

To find the coordinates of point A, which is the intersection of the two normal lines passing through points $P_i$ and $P_j$ to the path, a point-to-point combination method is utilized. The following are the steps:

• Generate points along the given path.
• Choose any two points ($P_i$ and $P_j$) on the given path and note the intersection of their normal lines as point E.
• Calculate the distance between point E and the remaining points noted as $P_k\left(k\ne i,j\right)$ on the given path.
• Check if the calculated distance satisfies the condition $l_{E{P}_j}<l_{E{P}_k}<l_{E{P}_i}$. If it does, then select point E as the crank center A. If not, re-select the combination of points $P_i$ and $P_j$.

However, it is essential to note that the location of point A can only be determined with sufficient accuracy as there are enough points on the given path. Spline interpolation is used by Yu [11] and Peng and Sodhi [45] to refine the path with additional points for problems with a limited number of points that form a given path by inserting points between neighboring points. In this paper, Fourier approximation is used to refine the path to more points. The time parameters are not considered as the generated points are not used to calculate the mechanism parameters in subsequence steps. This method relies on the Fourier coefficients instead of spline interpolation, resulting in improved accuracy of synthesis solutions and simplified computational procedure at the programming level. A comparison between these two point-generating methods is presented later in the paper by example 1 in Section 5. Once the position of point A has been determined, $l_{max}$ and $l_{min}$ are obtained, from which the $l_{AB}$ and $l_{BP}$ variables of the auxiliary mechanism are calculated as follows: $l_{AB}=\frac{\left(l_{max}-l_{min}\right)}{2}$, $l_{BP}=\frac{\left(l_{max}+l_{min}\right)}{2}$.

Once the parameters of the auxiliary mechanism are obtained, the relationship between the points on the given path and time parameters can be established. The method used by Yu [11] and Li et al. [29] involves generating new points while rotating the crank AB counterclockwise at an interval of a certain constant angle to make point P pass through the given path in one direction. However, this approach presents two issues. Firstly, it is difficult to find the point of intersection between the floating link BP and the given path. Secondly, assessing the accuracy of the generated results is challenging since it requires finding the closest coordinates to the given path in the generated path. To address these issues, a method is proposed that directly assigns the time parameters to the given points by calculating the corresponding value of $\phi$ for each point.

Assuming that the number of points for the given paths is $M$ and the given points are noted as $P_{di}(i=1, 2,\dots , M)$. As illustrated in Figure 4, the intersection points of the circle with radius $l_{AB}$ at point A, and the circle with radius $l_{BP}$ at point $P_{di}$ are denoted as $B_{di1}$ and $B_{di2}$. One of these points is the proper $B_{di}$ that corresponds to $P_{di}$, and the angle of the vector $\overrightarrow{A{B}_{di}}$ represents the corresponding time parameter. To ensure the appropriate selection of $B_{di}$, the following conditions must be met: if $P_{di}$ lies between $P_j$ and $P_i$, then the link $A{B}_{di}$ should be positioned in the counterclockwise interval of $A{P}_j$ and $A{P}_i$, and vice versa. By satisfying these conditions, M number of $B_{di}$ can be obtained. If the series of $B_{di}$ obtained do not rotate counterclockwise, discard this set of $P_i, P_j$ point combination. Otherwise, the time parameters are successfully assigned to the given points, and the Fourier coefficients can be obtained using Equation (18).

Figure 5 shows the flowchart for calculating the Fourier coefficients for path generation without prescribed timing presented in this section. According to this process, multiple sets of Fourier coefficients can be obtained along with the parameters $(x_A, y_A, l_{AB}, l_{BP})$ of the auxiliary mechanism in correspondence. This paper presents an improved method upon the work of Yu [11] and Li [29]. Rather than using point interpolation, the proposed method employs Fourier to obtain more accurate points along the given path. Moreover, time parameters are assigned directly to the given points, which avoids the need for complex calculations to find points in the given path and eliminates combinations that would lead to solutions with order defect. This method simplifies the whole synthesis procedure at the programming level by avoiding complex computation procedures in point generation and time parameter assignment. Additionally, it offers more accurate solutions.

4. Synthesis Procedure for Path Generation of Planar Four-Bar Linkages

The problem of path generation synthesis for four-bar linkages can be simplified into a problem of solving two polynomial equations, as explained in Section 2. For paths with prescribed timing, the synthesis procedure is straightforward. It involves obtaining Fourier coefficients, solving two sets of design equations (Equations (12) and (17)), and verifying results in kinematic analysis.

However, for paths without prescribed timing, the Fourier coefficients are obtained following the procedure in Figure 5, together with the parameters of the corresponding auxiliary mechanism. This auxiliary mechanism is part of the desired four-bar linkage [11], which means that the parameters of the input dyad are already obtained. Therefore, the calculation of the design equations Equation (12) can be skipped for faster calculation speed. The flowchart of the synthesis procedure is shown in Figure 6. The synthesis procedure is performed in MATLAB R2020a in Microsoft Windows 10 on an Intel Core i7-7500U CPU at 2.70 GHz.

The comparison between the proposed synthesis procedure and the one that does not consider auxiliary mechanism parameters and adopts the Fourier coefficients obtained by the auxiliary mechanism to follow the exact procedure as paths with prescribed timing is presented in Section 5.

To quantify the accuracy of the synthesis results, TE (tracking error) is used as follows:

$$\mathrm{T}\mathrm{E}=\frac{1}{M}{\sum }_{i=1}^M\sqrt{{\left(x_i-x_{di}\right)}^2+{\left(y_i-y_{di}\right)}^2}$$

(20)

where $(x_{di}, y_{di})$ are the given points and $(x_{di}, y_{di})$ are the generated points of the synthesized linkage corresponding with the prescribed time or the assigned time. The tracking error is the mean of the Euclidean distances for all given and generated points.

5. Applications and Discussions

5.1. Example 1: An Ellipse of 10 Points without Prescribed Timing

In this example, the desired path is an ellipse with 10 given points without prescribed timing. It was first presented and solved by Acharyya et al. [46] using evolutionary algorithms. This problem has been studied by many researchers using various methods since then. The coordinates of the given points are listed in

Table 1. Coordinates of the given points for example 1.

$\mathbf{Point} \mathbf{Number}, \mathit{i}$	1	2	3	5	6	7	8	9	10
$x_{di}$	20	17.66	11.736	5	0.60307	0.60307	5	11.736	17.66
$y_{di}$	10	15.142	17.878	16.928	12.736	7.2638	3.0718	2.1215	4.8577

.

The example is solved with and without considering time parameters to demonstrate the importance of time parameter assignment. The given path is refined to 100 points using Fourier approximation and spline interpolation to determine the auxiliary mechanism. When disregarding time parameter, a valid solution and its cognate are obtained and are listed in

Table 2. Synthesis results for example 1.

Design Variables	Solution 1 1	Solution 1 Cognate 1	Solution 1 2	Solution 2 2	Solution 1 3	Solution 2 3
$a$	8.9999	9.9998	8.0283	8.0283	9.2163	9.2163
$b$	9.0008	60,385.8909	104.9061	107.5520	296.7365	315.8426
$c$	54,347.8187	10.0008	88.0198	96.7131	278.4022	359.2052
$d$	54,347.8182	60,385.8903	167.1844	178.8989	489.2287	592.9863
$r_A$	14.1419	6028.0809	17.5653	17.5653	60.3629	60.3629
$f$	1.0000	6038.0722	22.5292	22.5292	57.5619	57.5619
$\mu$	0.7854	3.1399	−0.7829	−0.7829	−0.4673	−0.4673
$\alpha$	3.1415	0.0000	2.8611	−2.3222	−1.7631	−0.7260
$\beta$	0.0000	0.0000	−0.7099	−2.7514	−1.5479	2.6307
TE	0.5690	0.5690	0.0038	0.0021	0.0005	0.0003
Conf.	crossed	open	crossed	open	crossed	open

¹ Solutions obtained disregarding time parameters. ² Solutions obtained considering time parameters using spline interpolation. ³ Solutions obtained considering time parameters using Fourier approximation.

. The generated path does not fit the given points well. By considering non-uniform time parameters, 100 generated points result in 4950 point-to-point combinations for finding potential auxiliary mechanisms. In this example, we only consider reaching the farthest point first, while the auxiliary mechanism reaching the closest point first will be discussed in example 2. Using spline interpolation, 597 sets of point-to-point combinations are obtained, of which 273 sets are valid for generating given points counterclockwise. Among these sets, 111 sets have a point A outside the path, while 162 have it inside. Similarly, using Fourier approximation, 629 sets of combinations are obtained, of which 283 sets are valid, with 112 having a point A outside of the path and 171 sets having it inside. For comparison, the combination that can best perform a Fourier approximation to the given points is made from the valid sets of point-to-point combinations of each method for path synthesis. The synthesis results and corresponding time parameters using the proposed method are presented in Table 2 and

Table 3. Time parameters for example 1.

$\mathbf{Point} \mathbf{Number}, \mathit{i}$	1	2	3	5	6	7	8	9	10
${\mathsf{\phi }}_{\mathrm{i}}$ 1	−0.0168	0.6766	1.3661	2.0585	2.7542	3.4555	4.1600	4.8648	5.5677
${\mathsf{\phi }}_{\mathrm{i}}$ 2	−0.1107	0.5873	1.2870	1.9876	2.6880	3.3872	4.0847	4.7807	5.4763

¹ Time parameters obtained using spline interpolation. ² Time parameters obtained using Fourier approximation.

, respectively. Despite both methods yielding two solutions of high accuracy, it should be noted that using Fourier approximation results in a tracking error that is almost ten times smaller than that of spline interpolation.

The synthesis results with time parameters in Table 2 are calculated using the proposed method based on considering the auxiliary mechanism as the exact input dyad for synthesis solutions. The solutions are achieved by solving only one set of design equations (Equation (17)). In comparison, the method that requires solving two sets of design equations (Equations (12) and (17)) without considering the auxiliary mechanism parameters is also presented. Six solutions are obtained, including three solutions and their cognates. The results are listed in

Table 4. Synthesis results for example 1 without considering the auxiliary mechanism.

Design Variables	Solution 1	Solution 2	Solution 1 Cognate	Solution 3	Solution 2 Cognate	Solution 3 Cognate
$a$	9.7177	9.7177	9.2163	9.2163	8.0379	8.0379
$b$	293.5478	304.9371	296.7365	315.8426	305.7470	313.2784
$c$	312.8796	369.6413	278.4022	359.2052	252.2273	275.4601
$d$	515.8436	586.6166	489.2287	592.9863	485.2171	517.1691
$r_A$	41.2041	41.2041	60.3629	60.3629	77.1553	77.1553
$f$	54.0054	54.0054	57.5619	57.5619	65.4647	65.4647
$\mu$	−2.8664	−2.8664	−0.4673	−0.4673	1.3290	1.3290
${\phi }_0$	0.1816	0.1816	0.0000	0.0000	0.1392	0.1392
$\alpha$	1.1969	0.1993	−1.7631	−0.7260	−2.8999	2.2764
$\beta$	−1.3663	0.7777	−1.5479	2.6307	0.7353	2.7699
TE	0.0005	0.0007	0.0005	0.0003	0.0007	0.0003
Conf.	open	crossed	crossed	open	open	crossed

, where the two solutions (solution 1 cognate and solution 3) are the same as those obtained in Table 2. The cognate mechanisms share the same tracking error, i.e., only one more solution can be obtained using this method, but almost 10 times the computational time is used, as shown in

Table 5. Computational time for example 1.

Method	Computational Time (s) 1
Method considering the auxiliary mechanism	1.9500
Method without considering the auxiliary mechanism	17.1092

¹ The computational time is the average of 100 calculations for each method.

. Reselecting a set of point-to-point combinations using the proposed method for more synthesis solutions is the better option with high efficiency.

Great efforts have been made to solve path generation problems. Acharyya and Mandal [46] first solved the problem using evolutionary algorithms, Lin and Hsiao [6] presented a CMDE (combined-mutation differential evolution), Eqra et al. [47] used adaptive PSO (particle swarm optimization), Li et al. [29] presented a novel analytical method, and Torres-Moreno et al. developed an open-source tool based on NSDVs (normalized shape-descriptor vectors).

Table 6. Comparison of synthesis results for example 1.

Design Variables	Proposed Method	Acharyya. 2009 [46]	Lin. 2017 [6]	Eqra. 2018 [47]	Li. 2020 [29]	Moreno. 2022 [1]
$a$	9.2163	8.6834	8.0457	8.5320	8.8355	8.0700
$b$	315.8426	34.3186	50.8190	31.4848	109.0051	50.5900
$c$	359.2052	79.9962	42.2080	33.2131	118.0160	42.0100
$d$	592.9863	54.3609	80.0000	54.7218	199.3538	79.5700
$r_A$	60.3629	15.5770	8.5286	8.7906	23.4929	8.2257
$f$	57.5619	1.4653	10.8809	6.0145	20.4649	10.3833
$\mu$	−0.4673	0.7909	−0.0890	0.5897	−0.2535	−0.0791
$\alpha$	−0.7260	1.5707	−2.9294	1.5004	−1.0982	−2.9937
$\beta$	2.6307	2.1297	3.8892	0.0965	2.8291	3.9300
TE	0.0003	2.3432	0.0051	0.0145	0.0037	0.4490
Conf.	open	open	open	crossed	open	open

compares the solutions generated by these approaches, and the generated paths are shown in Figure 7. Recent studies have addressed this problem properly. However, the proposed method possesses the most minor tracking error. The kinematic diagram of the proposed method for example 1 is shown in Figure 8.

5.2. Example 2: A Teardrop Shape of 25 Points without Prescribed Timing

This example comes in the shape of a teardrop with a cusp and is presented by McGarva [13] as a path for a particular piece of packaging machine. The coordinates of the given points are listed in

Table 7. Coordinates of the given points for example 2.

$\mathbf{Point} \mathbf{Number}, \mathit{i}$	1	2	3	4	5	6	7	8	9	10	11	12	13
$x_{di}$	7.03	6.95	6.77	6.4	5.91	5.43	4.93	4.67	4.38	4.04	3.76	3.76	3.76
$y_{di}$	5.99	5.45	5.03	4.6	4.03	3.56	2.94	2.6	2.2	1.67	1.22	1.97	2.78
	14	15	16	17	18	19	20	21	22	23	24	25
$x_{di}$	3.76	3.76	3.76	3.76	3.8	4.07	4.53	5.07	5.45	5.89	6.41	6.92
$y_{di}$	3.56	4.34	4.91	5.47	5.98	6.4	6.75	6.85	6.84	6.83	6.8	6.58

.

For this example, no valid point-to-point combination is obtained considering the coupler point reaching the farthest point first in the path. However, by considering the coupler point reaching the closest point first, the parameters assigned to the given points are listed in

Table 8. Time parameters for example 2.

$\mathbf{Point} \mathbf{Number}, \mathit{i}$	1	2	3	4	5	6	7	8	9
${\phi }_i$	3.3675	3.6823	3.9022	4.1316	4.3962	4.6129	4.8593	4.9919	5.1496
	10	11	12	13	14	15	16	17	18
${\phi }_i$	5.3695	5.6202	0.1458	0.4443	0.7182	0.9979	1.2165	1.4516	1.6822
	19	20	21	22	23	24	25
${\phi }_i$	1.9178	2.1792	2.3998	2.5443	2.7071	2.8931	3.1029

. The synthesis result is listed in

Table 9. Synthesis results for example 2.

Design Variables	Proposed Method	Zhang. 2019 [37]	Li. 2020 [29]	Hernández. 2021 [35]
$a$	1.7869	1.92715	2.01440	2.44870
$b$	4.2405	4.65983	5.54970	5.07450
$c$	4.0146	7.27223	6.71420	6.61420
$d$	6.2684	10.00000	10.10860	8.18640
$r_A$	10.5236	3.82046	4.54090	2.35603
$f$	9.2977	8.18655	9.74180	8.53541
$\mu$	−0.2698	2.28210	2.64080	2.95265
$\alpha$	0.1621	0.19171	0.10880	−0.35591
$\beta$	2.6661	5.66686	−0.45260	−0.06010
TE	0.0371	0.05201	0.04951	0.05568
Conf.	crossed	open	open	open

, and the kinematic diagram is shown in Figure 9.

This example has been studied by many researchers over the years. Table 9 compares the solutions generated by different approaches. Zhang et al. [29] presented an HLIDE (Lagrange interpolation differential evolution) algorithm for path synthesis of four-bar mechanisms. Hernández et al. [35] reformulated the error function for the gradient method for better results. As shown in Table 9, each method demonstrates competitive accuracy, but the proposed method presents better performance. The comparison of the paths generated by different approaches is shown in Figure 10 intuitively.

5.3. Example 3: A Complex Triple Loop Path of 90 Points without Prescribed Timing

This example is a complex triple loop path with two crunodes, presented by Hadizdeh Kafash and Sahand [48] to test the performance of their method. The given points are generated by an existing linkage with the following parameters: $a=3.1, b=5, c=8.6, d=10.4, r_a=0, \mu =0, r_f=0, \alpha =1, \beta =0$. The crank rotates counterclockwise at an interval of 4° for generating the points. The generated points are given as a path synthesis problem without prescribed timing. The solutions obtained with and without considering time parameters, Hadizdeh Kafash and Sahand’s solution [48], and the parameters of the original linkage are listed in

Table 10. Synthesis results for example 3.

Design Variables	Proposed Method	Disregarding Time Parameters	Hadizdeh Kafash 2017 [48]	Original Linkage
$a$	3.1609	3.1076	3.1092	3.1
$b$	5.0157	5.0012	5.0394	5
$c$	9.3046	8.6050	8.7281	8.6
$d$	11.0445	10.4199	10.5479	10.4
$r_A$	0.1542	0.0309	0.0488	0
$f$	5.9161	5.9765	6.0387	6
$\mu$	2.9681	2.4404	−1.4744	0
${\phi }_0$	0.0000	−0.0030	0.0000	0
$\alpha$	0.9741	1.0026	0.9993	1
$\beta$	−0.0256	−0.0019	0.0000	0
TE	0.0125	0.0096	0.0194	0
Conf.	open	open	open	open

. The synthesis results are intuitively compared in Figure 11.

These approaches produce solutions with decent accuracy. The method that does not consider time parameters results in a minimal tracking error because other methods treat the generated points as a path without prescribed timing, although they are generated using uniform spacing timing. The proposed method is very effective in solving this problem, with a smaller tracking error compared to Hadizdeh Kafash and Sahand’s [48] method. Figure 12 shows the kinematic diagram of the proposed method, which has location and design parameters that are very similar to those of the original linkage. This example is a challenging problem as it is a complex path with two crunodes, thus demonstrating the efficacy of the proposed method.

6. Conclusions

In this paper, an analytical synthesis method based on the Fourier series for path generation is improved. Two sets of design equations are generated based on loop closure equations for the input dyad and the linkage itself, and in solving these equations, the solutions are obtained. A more practical method to assign time parameters is introduced to address issues in previous papers for generating paths without prescribed timing. Three numerical examples are presented to assess the synthesis approach. The approach is easy to understand and computationally efficient, producing more accurate solutions compared with other known methods. The contributions of this paper are as follows:

• Fourier approximation is used instead of spline interpolation to better capture the curve of the given points, although it does not contain the given points. The synthesis solutions had improved accuracy compared to those using spline interpolation. In addition, as the foundation of this analytical synthesis is the Fourier series, using Fourier approximation, consistent with the foundation, simplifies the computational procedure at the programming level. It can be concluded that it is computationally efficient and accurate to generate points readily for the point-to-point combination method using Fourier approximation.
• Time parameters are assigned directly to the given points, avoiding complex calculations to find intersection points along the given path with the floating link of the auxiliary mechanism. The conditions proposed to find valid sets of point-to-point combinations whilst discarding combinations that lead to solutions with order defects before calculating design parameters. In addition, the accuracy evaluation process is simplified by the fact that the time parameters for the given points are known, thus eliminating the need to find the corresponding generated path point for the given points.
• The parameters obtained by the point-to-point combination method are used as the parameters of the input dyad, skipping the first set of design equations presented in Section 2.2 for faster calculation. Examples demonstrate that this method is highly effective in the synthesis process with competitive accuracy. The computation time is significantly reduced by one order of magnitude compared to that of solving two sets of design equations.