This paper studies the problem of spherical four-bar motion synthesis from the viewpoint of acquiring circular geometric constraints from a set of prescribed spherical poses. The proposed approach extends our planar four-bar linkage synthesis work to spherical case. Using the image space representation of spherical poses, a quadratic equation with ten linear homogeneous coefficients, which corresponds to a constraint manifold in the image space, can be obtained to represent a spherical RR dyad. Therefore, our approach to synthesizing a spherical four-bar linkage decomposes into two steps. First, find a pencil of general manifolds that best fit the task image points in the least-squares sense, which can be solved using singular value decomposition (SVD), and the singular vectors associated with the smallest singular values are used to form the null-space solution of the pencil of general manifolds; second, additional constraint equations on the resulting solution space are imposed to identify the general manifolds that are qualified to become the constraint manifolds, which can represent the spherical circular constraints and thus their corresponding spherical dyads. After the inverse computation that converts the coefficients of the constraint manifolds to the design parameters of spherical RR dyad, spherical four-bar linkages that best navigate through the set of task poses can be constructed by the obtained dyads. The result is a fast and efficient algorithm that extracts the geometric constraints associated with a spherical motion task, and leads naturally to a unified treatment for both exact and approximate spherical motion synthesis.

## Introduction

A spherical four-bar mechanism is a closed chain formed by two spherical dyads, which are connected by revolute joints and rotate about the same fixed point. This paper deals with the problem of synthesizing a spherical four-bar mechanism from the perspective of kinematic extraction of geometric constraints from the task motion. Kinematic motion synthesis is a well-known concept in computational kinematics for planar and spatial mechanisms, and it plays an important role in the fields of mechanism design, robotics, dynamic control, motion planning, and automatic assembly; see Refs. [1–8] for established methods and applications in motion synthesis of mechanisms. While the exact motion synthesis problem of four or five prescribed poses has been well studied, most of the current solutions to the synthesis of mechanisms for approximate motion (an arbitrary number of poses), however, require complex and time consuming algorithms. In our previous work, we proposed a simple and efficient algebraic method for planar four-bar linkages synthesis [9–11], which uncovers the geometric constraint hidden in the given motion via a linear, two-step method. During the synthesis process, kinematic mapping was employed to express the given motion and define the constraint equations. The concept of kinematic mapping of SE(2) was proposed by Blaschke [12] and Grünwald [13] almost a century ago, and was explored by Ravani and Roth [14] for motion approximation. In the kinematic mapping approach to kinematic synthesis, both planar and spherical poses in Cartesian space can be mapped into points in a three-dimensional projective space (called *image space* of planar or spherical kinematics), while workspace constraints of a mechanism map into algebraic manifolds in the same space. In this way, a single degree-of-freedom (DOF) motion of a planar or spherical mechanism is represented by the intersection curve of two algebraic manifolds. The problem of motion approximation is transformed into ab algebraic curve fitting problem in the image space, where various methods in approximation theory may be applied. This includes the definition of the approximation error (called fitting error) in the image space, formulation of a least-squares problem, and application of appropriate numerical methods to find values of the design variables for minimization of the error. Pursuant to Ravani and Roth's kinematic mapping approach for mechanism synthesis, further research has been done by Bodduluri and McCarthy [15], Bodduluri [16], Ge and Larochelle [17], Husty et al. [6], Wu et al. [18], and Purwar et al. [19], both of which have demonstrated a visual, computer graphics approach for multi-degrees-of-freedom mechanism design that exploits the connection between constraint manifold geometry and its apparent effect on the parameters of a mechanism to interactively perform kinematic synthesis. Hayes and Rusu [20] have presented preliminary results for combining type and dimensional synthesis of planar mechanisms for multipose rigid body guidance.

In this paper, we seek to extend our approach for planar case to spherical four-bar linkage synthesis. Existing works on this topic includes [17,21–24], in which Ruth and McCarthy [23] described the implementation of spherical four-orientation synthesis in the software sphinxpc, which encodes a new formulation of classical Burmester theory based on the equation of a spherical triangle, which yields a convenient parameterized equation for the central-axis cone. Purwar et al. [21] brought together the kinematics of spherical robot arms and freeform rational motions to study the problem of synthesizing constrained rational motions for Cartesian motion planning, and realized the synthesis of spherical 2R and 3R robot arms. Also based on kinematic mapping, Husty et al. [22] proposed an approach to the five spherical position synthesis by converting the design problem into a polynomial of degree six. Most of these works either focus only on the exact motion synthesis or involve a great amount of computation for approximate synthesis.

In our approach, we first formulate the data fitting problem into a linear fitting system through kinematic mapping, with the goal to find the quadratic general manifolds in the image space that best fit the given image points in the least-squares sense. Using singular value decomposition (SVD), the singular vectors associated with the smallest singular values are chosen to form the null-space solution of general manifolds. Second, four additional constraint equations are then imposed to identify the constraint manifolds qualified for representing a spherical dyad from the previously obtained general manifolds. After the inverse computation that converts the constraint manifolds to the design parameters of the corresponding dyads, spherical four-bar linkages that best guide through the set of task poses can be obtained with the resulting dyads as building blocks. Both exact and approximate syntheses are addressed in a unified framework, and the main computational cost (the null-space computation of the linear system) remains nearly constant irrespective of the number of given poses.

The organization of the paper is as follows: Section 2 reviews the concept of kinematic mapping and image space in so far as necessary for the development of this paper. Section 3 relates the spherical circular constraint to the spherical dyad motion and represents the spherical circular constraint in the image space as the constraint manifold. Section 4 deals with the spherical dyad synthesis for both exact and approximate motion realization with our approach, in which a pencil of general manifolds are first found that best fit the task image points in least-squares sense, and the constraint manifold qualified for spherical circular constraint are identified by imposing additional constraint equations. In Sec. 5, an example of classic spherical four-bar linkage coupler motion is first presented to show the validity of our approach, followed by two examples which address an approximate motion and an exact motion, respectively, to illustrate the efficacy of our unified approach.

## Parameterizing a Spherical Pose

**s**= (

*s*,

_{x}*s*,

_{y}*s*) denote a unit vector in the direction of the rotation axis, and

_{z}*θ*denote the rotation angle, we can represent the spatial rotation with Euler–Rodrigues parameters

*q*

_{1},

*q*

_{2},

*q*

_{3},

*q*

_{4}).

**q**is called a unit quaternion when $q12+q22+q32+q42=1$. See Bottema and Roth [1] and McCarthy [2] for more details on quaternions and the ensuing discussion.

**q**= (

*q*

_{1},

*q*

_{2},

*q*

_{3},

*q*

_{4}) is said to define a point in a projective three-dimensional space called the

*image space*of spherical displacement. A spherical displacement is represented by a point in image space in this manner, and therefore a one or two degree-of-freedom spherical motion is represented by a curve or surface in image space.

## Constraining a Spherical Pose With Spherical Dyad

As illustrated in Ref. [7], the analysis and synthesis of spherical linkages are examined on the unit sphere. A spherical RR dyad, as shown in Fig. 1, consists of a floating link connected to a crank by a revolute joint, which in turn is connected to ground by a revolute joint. To define the spherical pose of the floating link on the unit sphere, we introduce a fixed frame *F* and a moving frame *M* attached to the floating link. In the figure, the moving frame *M* is drawn on the surface of the unit sphere, though it is understood that the origins of both the moving and fixed frames are located at the center of the sphere. Let the unit vector directed along the fixed joint axis be denoted by *A*, the unit vector along the moving joint axis be *B*, and angular length of the crank be *α*.

*x*,

*y*,

*z*) denote the moving frame coordinates of a point on the moving link and satisfy the relation

*x*

^{2}+

*y*

^{2}+

*z*

^{2}= 1. Then, the coordinates of the point in fixed frame

*F*are given by

where [*R*] is given by Eq. (3) and represents the rotation of *M* relative to *F*. Note that the fixed frame coordinates (*X*, *Y*, *Z*) satisfy the relation *X*^{2} + *Y*^{2} + *Z*^{2} = 1 as well, since [*R*] is an orthonormal matrix.

*X*,

*Y*,

*Z*) satisfy the following equations:

*R*] always keeps the coordinates (

*X*,

*Y*,

*Z*) on the unit sphere. Therefore, what needs to be imposed is only the plane equation given by Eq. (5). In other words, the trajectory of (

*X*,

*Y*,

*Z*) will be a spherical circle as long as they satisfy Eq. (5). The center of the spherical circle lies exactly on the fixed joint axis A of the spherical RR dyad, which is given by the unit vector representation in the fixed frame

*F*

*d*≠ 0 and sgn(

*d*) denote the sign function of

*d*. In case of

*d*= 0, (

*X*,

_{A}*Y*,

_{A}*Z*) is either given as $((a/a2+b2+c2),(b/a2+b2+c2),(c/a2+b2+c2))$ or $(\u2212(a/a2+b2+c2),\u2212(b/a2+b2+c2),\u2212(c/a2+b2+c2))$. The radius of the circle constraint is given by

_{A}*image space*with ten homogeneous terms

*a*,

*b*,

*c*,

*d*,

*x*,

*y*,

*z*} of the spherical RR dyad are included in ten homogeneous coefficients

*p*($i=1,2,\u202610$), with (

_{i}*a*,

*b*,

*c*,

*d*) as the homogeneous coordinates of the plane constraint and (

*x*,

*y*,

*z*) as the coordinates of a point on the moving link relative to the moving frame

*M*. The numbers of independent design parameters and independent coefficients

*p*are five and nine, respectively. Thus, four relationships should exist among those ten homogeneous coefficients. Aside from the stand-alone

_{i}*p*

_{10}, the rest of coefficients

*p*($i=1,2,\u20269$) satisfy the following conditions:

_{i}*p*

_{1}to

*p*

_{9}. In this paper, we take the following four relations to represent Eq. (11):

*p*

_{2}

*p*

_{5}

*p*

_{8}= 0, these four relations become necessary, but not a sufficient, conditions for Eq. (11).

To sum up, Eq. (9) defines a quadric surface algebraically in the image space with ten homogeneous coefficients *q _{i}* ($i=1,2,\u202610$). In this paper, we call this quadric a generalized constraint manifold, or

*G*-

*manifold*in short. For this generalized manifold to become the candidate for the constraint manifold or

*C*-

*manifold*of a spherical RR dyad, one must impose the four additional conditions in Eq. (12) on the coefficients

*p*.

_{i}**p**= {

*p*

_{1},

*p*

_{2},…,

*p*

_{10}}, which satisfy the relations in Eq. (12) (If

*p*

_{2}

*p*

_{5}

*p*

_{8}= 0, a further test is needed to verify the relations given by Eq. (11)). Then, they are qualified to be the coefficients of C-manifold as in Eq. (10), and hereby utilized to determine the design parameters of the corresponding spherical RR dyad. The inverse calculation from the coefficients

*p*to the design parameters {

_{i}*a*,

*b*,

*c*,

*d*,

*x*,

*y*,

*z*} could be obtained by

Thus, a spherical circle constraint can now be represented in terms of the C-manifold defined by Eqs. (9) and (10). To visually examine the C-manifold, we project its quadric equation into the three-dimensional space parameterized by (*q*_{1}/*q*_{4}, *q*_{2}/*q*_{4}, *q*_{3}/*q*_{4}). It can be shown that the projected C-manifold is a hyperboloid with one sheet, with one specific example shown in Fig. 2.

## Spherical Dyad Synthesis for N-Pose Realization

In Sec. 3, we have cast the representation of the spherical circular constraint associated with the spherical RR dyad in the *image space* as the constraint manifold, or C-manifold. Based on this formulation, we propose a simple and fast approach, which can be decomposed into two steps: first is finding a pencil of G-manifolds that best fit the given poses in least-squares sense, and second is identifying the C-manifolds from the pencil of G-manifolds whose coefficients are qualified to be the combinations of dimensional parameters of a spherical circular constraint as indicated by Eq. (10).

### Least-Squares Fitting of a Pencil of G-Manifolds.

*N image space*points that correspond to

*N*spherical poses in Cartesian space. This problem can be formulated as a linear algebraic fitting problem [

*A*]

**p**= 0 obtained by substituting the values of

*N*task image points for [

*A*] using Eq. (9), and representing the homogeneous coefficients

*p*(

_{i}*i*= 1…10) as the column vector

**p**. Therefore, the coefficient matrix [

*A*] of this linear problem is given by

*i*th row of matrix [

*A*] is defined as

**p**= [

*p*

_{1},

*p*

_{2},…,

*p*

_{10}] is obtained from the linear system (14) to yield the solution of a pencil of G-manifolds. The null-space problem is linear and can be readily solved using various algebraic methods such as QR decomposition and Gaussian elimination. In this paper, we apply singular value decomposition method to solve this linear fitting problem and measure the fitting error in the least-squares sense. In linear algebra, the SVD of the

*N*× 10 matrix [

*A*] is a factorization of the form

where [*U*] is an *N* × *N* orthonormal matrix, whose columns are called the *left singular vectors* of [*A*], and are also the eigenvectors of [*A*][*A*]^{T}; [*S*] is an *N* × 10 rectangular diagonal matrix with nonzero singular values of [*A*] on the diagonal, which are also the square roots of nonzero eigenvalues of [*A*][*A*]^{T} (or equivalently [*A*]^{T}[*A*]); [*V*] is an 10 × 10 orthonormal matrix, whose columns are called the *right-singular vectors* of [*A*], and are also the eigenvectors of [*A*]^{T}[*A*].

The application of the SVD algorithm is that it provides an explicit representation of the null space of the matrix [*A*]. The right-singular vectors corresponding to vanishing singular values of [*A*], or zero eigenvalues of [*A*]^{T}[*A*], span the null space of [*A*]. For the known five-pose synthesis problem that a spherical RR dyad can go through at most five arbitrary spherical poses exactly (see Ref. [7]), the size of [*A*] is 5 × 10, indicating that [*S*] has five nonzero singular values on the diagonal and thus the number of vanishing or zero singular values is five. Consequently, the right-singular vectors corresponding to those zero singular values form a five-dimensional exact null space. If the number of poses increases, there could be no zero singular values. In this case, we take the singular vectors associated with the least singular values to form the approximate null space, since the squared value of the fitting error of a singular vector is equal to its corresponding singular value (see our previous work [25] for proof).

*k*singular vectors

**v**

*of the least singular values are used to form the basis for the approximate null space, with their corresponding singular values denoted by*

_{k}*σ*then, any solution vector in this null space is given by

_{k}### Identification of C-Manifolds From a Pencil of G-Manifolds.

The least-squares solutions of [*A*] **p** = 0 might not be qualified to represent the C-manifolds of a spherical RR dyad, because they may not satisfy the four relations in Eq. (12). We now consider the problem of identifying those special quadrics within the resulting pencil of quadrics that correspond to the C-manifold.

**v**

_{1}–

**v**

_{5}of the least singular values

*α*($i=1,2,\u20265$) denote five real homogeneous coefficients. Substituting Eq. (19) into Eq. (12), we obtain four homogeneous quadratic equations

_{i}where each coefficient *k _{ij}* is a constant determined by components of the five singular vectors

**v**

_{1}–

**v**

_{5}. Note that

*F*are homogeneous in

_{i}*α*, which means that if {

_{k}*α*

_{1},

*α*

_{2},

*α*

_{3},

*α*

_{4},

*α*

_{5}} is a solution of these equations, then {

*kα*

_{1},

*kα*

_{2},

*kα*

_{3},

*kα*

_{4},

*kα*

_{5}} is a solution as well. Thus, we can set one of

*α*equal to one, say

_{k}*α*

_{1}= 1, such that these four homogeneous equations become nonhomogeneous. A number of existing numerical algorithms could handle a group of these quadratic equations. One example is the NSolve command in mathematica software, which can solve our four equations instantly.

*p*

_{2}

*p*

_{5}

*p*

_{8}of the identified solutions is equal to zero, they may not be qualified to become C-manifolds. As discussed in Sec. 3, the four relations in Eq. (12) are sufficient and necessary to capture the nine relations in Eq. (11), except for the case that

*p*

_{2}

*p*

_{5}

*p*

_{8}= 0 and (12) are just necessary but not sufficient conditions for Eq. (11). To readily exclude those solutions, which satisfy Eq. (12) but not Eq. (11), we have the following structural error:

**p**

^{T}

**p**= 1 since

**p**is homogeneous.

### Realization of Exact Motion.

Burmester theory [26] states that a planar dyad can go through at most five arbitrary planar poses exactly. The spherical version of Burmester theory, as the planar counterpart, states that up to five arbitrary spherical poses can be realized exactly by a spherical dyad (see Chiang [27] for details).

In our formulation, the number of vanishing, or zero eigenvalues of the matrix [*A*] is 10 − 5 = 5 with five specified poses. Therefore, [*A*] has a five-dimensional exact null space, which is defined by its five orthonormal right-singular vectors (also the eigenvectors of [*A*]^{T}[*A*]), associated with the five zero singular values. Thus, a G-manifold solution **p** can be expressed by a linear combination of those singular vectors as shown in Eq. (19). By observing Eq. (18), we find that the solution's fitting error *e _{f}* is equal to zero regardless of the values of

*α*. Furthermore, for obtaining the solutions of G-manifold

_{k}**p**that are qualified to become C-manifolds of spherical RR dyad, the coefficients

*α*has to satisfy Eq. (20). As mentioned in Sec. 4.2, we can use NSolve function in Mathematica to numerically solve for

_{k}*α*that eventually are used to construct the C-manifolds, yielding zero fitting error for the five arbitrary spherical poses.

_{k}### Realization of Approximate Motion.

Toward six or more arbitrarily given spherical poses, in general, one can only find approximate solutions, i.e., spherical dyads that could only go through the given poses approximately, which is classified as approximate motion synthesis.

In this problem, the number of zero eigenvalues of the matrix [*A*] is 10 − *N* with *N* > 5 specified poses. Therefore, [*A*] has a five-dimensional approximate null space with the basis vectors given by the five orthonormal singular vectors associated with the five near-zero singular values, and a G-manifold solution **p** can be expressed by their linear combination. Again, for obtaining the solutions of G-manifold **p** that are qualified to become C-manifolds of spherical RR dyad, we impose the relations in Eq. (20) to determine the coefficients *α _{k}*. By checking Eq. (18) against the resulting

*α*, the algebraic fitting error is found to be larger than but meanwhile close to zero, as the chosen basis vectors correspond to the smallest (near-zero) singular values. As a consequence, the resulting dyads will navigate through the

_{k}*N*specified poses as closely as possible.

It should be pointed out that since the size of [*A*]^{T}[*A*] remains 10 × 10 no matter how many poses are prescribed; thus, the main computational cost (i.e., the computing of eigenvectors of [*A*]^{T}[*A*] remains constant even if the realization of more task poses are required. In short, this approach leads to a unified algorithm for both exact (*N* < 5) and approximate (*N* > 5) spherical motion synthesis.

### Circuit/Assembly Mode of Resulting Linkages.

In practical design applications, aside from the realization of task positions, the circuit/assembly mode should also be considered for spherical four-bar linkages, i.e., it should be investigated that if all the task positions belong to the same circuit/assembly mode.

As discussed in Sec. 3, a spherical dyad constraint is converted to a quadric in the form of Eq. (9), which could be represented by a hyperboloid-shaped C-manifold in the image space. The coupler motion of a spherical four-bar linkage, which is constrained by two spherical dyads, could be represented by the intersection curve of two such C-manifolds. Different circuits/assembly modes are illustrated in the image space as separate intersection curves (see Fig. 3). Since the task positions are converted to a set of points in the image space, it is clear to identify the circuits/assembly modes of the coupler motion by investigating the intersection curves, and hereby determine if all task positions lie within the same circuit of the spherical four-bar linkage. If only one intersection curve exists, as in the right figure of Fig. 3, it means that there exists only one circuit/assembly mode, i.e., it is a non-Grashof linkage; if two intersection curves exist as in the left figure, it implies that the resulting spherical four-bar linkage contains two circuits, and further test is needed to tell whether task image points locate on the same intersection. Circuit defect happens if the task image points are located separately on two intersection curves.

To avoid the circuit defect of the resulting spherical four-bar linkage, during the above design process in Secs. 4.1–4.5, one needs to make sure that all the task positions lie on the same intersection curve once a four-bar is formed by the resulting spherical RR dyads. Considering that the C-manifold of a spherical RR dyad is hyperboloid-shaped, same as the C-manifold of a planar RR dyad, the approach presented in our previous work [28] can be used to carry out the test.

## Examples

In this section, we present three examples for spherical four-bar motion synthesis using the geometric constraint acquisition method. In the first example, the task motion is a sequence of poses generated from a spherical four-bar mechanism with known dimensions to demonstrate the validity of the proposed method. In the last two examples, an approximate motion with twelve poses and an exact motion with five poses are given to illustrate the effectiveness of our method to both approximate and exact motion synthesis.

### Example 1: Synthesis of a Known Spherical Four-Bar Motion.

First presented is an example of generation of a known spherical four-bar coupler motion whose linkage dimensions are given as follows: the two fixed pivots are located at (−1, 0, 0) and (0, −1, 0) on the sphere, respectively, and the angular link lengths of the input, coupler and output link are 30 deg, 60 deg, and 75 deg, respectively. We sample 12 poses from this motion as listed in Table 1, which are also plotted in Fig. 4.

Rotation axis | Angle (deg) | (q_{1}, q_{2}, q_{3}, q_{4}) | |
---|---|---|---|

1 | (0.273, 0.484, 0.831) | 128.36 | (0.246, 0.436, 0.749, 0.436) |

2 | (0.211, 0.454,0.866) | 122.54 | (0.185, 0.398, 0.759, 0.481) |

3 | (0.160, 0.426, 0.891) | 114.76 | (0.135, 0.358, 0.750, 0.539) |

4 | (0.129, 0.400, 0.907) | 105.82 | (0.103, 0.319, 0.724, 0.603) |

5 | (0.122, 0.383, 0.916) | 96.44 | (0.091, 0.286, 0.683, 0.666) |

6 | (0.140, 0.381, 0.914) | 87.15 | (0.097, 0.263, 0.630, 0.725) |

7 | (0.185, 0.403, 0.896) | 78.47 | (0.117, 0.255, 0.567, 0.775) |

8 | (0.249, 0.460, 0.852) | 71.14 | (0.145, 0.268, 0.500, 0.813) |

9 | (0.322, 0.556, 0.767) | 66.32 | (0.176, 0.304, 0.420, 0.837) |

10 | (0.378, 0.671, 0.638) | 65.36 | (0.204, 0.362, 0.345, 0.842) |

11 | (0.400, 0.769, 0.499) | 68.87 | (0.226, 0.435, 0.282, 0.825) |

12 | (0.400, 0.825, 0.400) | 75.72 | (0.246, 0.506, 0.246, 0.790) |

Rotation axis | Angle (deg) | (q_{1}, q_{2}, q_{3}, q_{4}) | |
---|---|---|---|

1 | (0.273, 0.484, 0.831) | 128.36 | (0.246, 0.436, 0.749, 0.436) |

2 | (0.211, 0.454,0.866) | 122.54 | (0.185, 0.398, 0.759, 0.481) |

3 | (0.160, 0.426, 0.891) | 114.76 | (0.135, 0.358, 0.750, 0.539) |

4 | (0.129, 0.400, 0.907) | 105.82 | (0.103, 0.319, 0.724, 0.603) |

5 | (0.122, 0.383, 0.916) | 96.44 | (0.091, 0.286, 0.683, 0.666) |

6 | (0.140, 0.381, 0.914) | 87.15 | (0.097, 0.263, 0.630, 0.725) |

7 | (0.185, 0.403, 0.896) | 78.47 | (0.117, 0.255, 0.567, 0.775) |

8 | (0.249, 0.460, 0.852) | 71.14 | (0.145, 0.268, 0.500, 0.813) |

9 | (0.322, 0.556, 0.767) | 66.32 | (0.176, 0.304, 0.420, 0.837) |

10 | (0.378, 0.671, 0.638) | 65.36 | (0.204, 0.362, 0.345, 0.842) |

11 | (0.400, 0.769, 0.499) | 68.87 | (0.226, 0.435, 0.282, 0.825) |

12 | (0.400, 0.825, 0.400) | 75.72 | (0.246, 0.506, 0.246, 0.790) |

*A*] constructed by Eq. (14) is 12 × 10. The first step is to find the singular values and their corresponding singular vectors of [

*A*] via SVD. The singular values are listed below in ascending order of magnitude

Note that there are two singular values that are almost zero, resulting from the fact that the given data are perfectly four-bar motion data as their image points lie exactly on the intersection curve of two hyperboloids that correspond to the two dyads of the known four-bar linkage.

Next, we take five singular vectors **v**_{1}–**v**_{5} associated with the five smallest singular values, as shown listed in Table 2, to form the null-space basis of [*A*] and define the G-manifold vector as in Eq. (19). By further requiring that the coefficients *α _{k}* ($i=1,2,\u20265$) satisfy the four equations in Eq. (20) and excluding the solutions that would cause big structural error

*e*, we obtain the valid solutions of

_{s}*α*qualified for the coefficients of the C-manifold, and their corresponding fitting errors

_{k}*e*as well as structural errors

_{f}*e*are listed in Table 3. From the column of

_{s}*e*, we can see that all the solutions satisfy Eq. (12) perfectly, and thus qualify to be the coefficients of C-manifold; also according to the fitting errors

_{s}*e*, the C-manifolds defined by the first two groups of solutions fit the given data perfectly, while the last two groups result in some error. Therefore, we adopt the first two groups as our final choices and compute their C-manifold vectors by Eq. (19). The two resulting C-manifolds are plotted in Fig. 5 where the 12 image points in the figure lie on the intersection of the two manifolds and represent 12 task poses in Table 1.

_{f}p_{1} | p_{2} | p_{3} | p_{4} | p_{5} | p_{6} | p_{7} | p_{8} | p_{9} | p_{10} | |
---|---|---|---|---|---|---|---|---|---|---|

v_{1} | 0 | −0.3532 | 0.6117 | 0 | 0.1135 | 0.1966 | 0 | 0 | 0 | −0.6705 |

v_{2} | 0 | −0.1486 | 0.2574 | 0 | −0.4774 | −0.8269 | 0 | 0 | 0 | −0.0102 |

v_{3} | −0.1588 | 0.4380 | −0.2040 | 0.5451 | 0.0734 | −0.1789 | −0.1230 | 0.1713 | 0.3920 | −0.4568 |

v_{4} | −0.3286 | 0.2743 | 0.3329 | −0.2232 | 0.1685 | −0.0452 | 0.2349 | −0.5358 | 0.5029 | 0.1745 |

v_{5} | −0.5259 | −0.3223 | −0.2177 | −0.1826 | 0.5505 | −0.3273 | 0.1820 | 0.2983 | −0.0748 | −0.0317 |

p_{1} | p_{2} | p_{3} | p_{4} | p_{5} | p_{6} | p_{7} | p_{8} | p_{9} | p_{10} | |
---|---|---|---|---|---|---|---|---|---|---|

v_{1} | 0 | −0.3532 | 0.6117 | 0 | 0.1135 | 0.1966 | 0 | 0 | 0 | −0.6705 |

v_{2} | 0 | −0.1486 | 0.2574 | 0 | −0.4774 | −0.8269 | 0 | 0 | 0 | −0.0102 |

v_{3} | −0.1588 | 0.4380 | −0.2040 | 0.5451 | 0.0734 | −0.1789 | −0.1230 | 0.1713 | 0.3920 | −0.4568 |

v_{4} | −0.3286 | 0.2743 | 0.3329 | −0.2232 | 0.1685 | −0.0452 | 0.2349 | −0.5358 | 0.5029 | 0.1745 |

v_{5} | −0.5259 | −0.3223 | −0.2177 | −0.1826 | 0.5505 | −0.3273 | 0.1820 | 0.2983 | −0.0748 | −0.0317 |

α_{1} | α_{2} | α_{3} | α_{4} | α_{5} | e_{s} | e_{f} | |
---|---|---|---|---|---|---|---|

Solution 1 | 1 | 0.2378 | −9.08 × 10^{−9} | −3.06 × 10^{−11} | −1.24 × 10^{−12} | 3.13 × 10^{−9} | 2.50 × 10^{−16} |

Solution 2 | 1 | −2.3769 | −1.01 × 10^{−8} | 5.66 × 10^{−10} | −1.26 × 10^{−10} | 1.21 × 10^{−8} | 2.40 × 10^{−16} |

Solution 3 | 1 | 10.4409 | 6.4721 | −0.3098 | −0.3305 | 3.63 × 10^{−9} | 0.0006 |

Solution 4 | 1 | 0.0782 | 0.0162 | −0.6113 | −0.3588 | 6.79 × 10^{−9} | 0.0064 |

α_{1} | α_{2} | α_{3} | α_{4} | α_{5} | e_{s} | e_{f} | |
---|---|---|---|---|---|---|---|

Solution 1 | 1 | 0.2378 | −9.08 × 10^{−9} | −3.06 × 10^{−11} | −1.24 × 10^{−12} | 3.13 × 10^{−9} | 2.50 × 10^{−16} |

Solution 2 | 1 | −2.3769 | −1.01 × 10^{−8} | 5.66 × 10^{−10} | −1.26 × 10^{−10} | 1.21 × 10^{−8} | 2.40 × 10^{−16} |

Solution 3 | 1 | 10.4409 | 6.4721 | −0.3098 | −0.3305 | 3.63 × 10^{−9} | 0.0006 |

Solution 4 | 1 | 0.0782 | 0.0162 | −0.6113 | −0.3588 | 6.79 × 10^{−9} | 0.0064 |

Finally, by applying the inverse computation Eq. (13) on the two resulting C-manifolds, the dimensions of the two resulting spherical circular constraints are also listed in Table 4. It is easily verified that the two dyads combined together constrain the motion of the spherical four-bar linkage as given.

C-manifold | a:b:c:d | x | y | z | Angular link length (deg) |
---|---|---|---|---|---|

p_{1} | 1:0.0000:0.0000:0.8660 | 0.0000 | 0.5000 | −0.8660 | 30.0029 |

p_{2} | 0.0000:1:0.0000:0.2588 | 0.0000 | −0.5000 | −0.8660 | 75.0011 |

C-manifold | a:b:c:d | x | y | z | Angular link length (deg) |
---|---|---|---|---|---|

p_{1} | 1:0.0000:0.0000:0.8660 | 0.0000 | 0.5000 | −0.8660 | 30.0029 |

p_{2} | 0.0000:1:0.0000:0.2588 | 0.0000 | −0.5000 | −0.8660 | 75.0011 |

### Example 2: Synthesis of an Approximate Spherical Motion (12 Poses).

In this section, an approximate example is demonstrated, whose prescribed poses are obtained by truncating the rotation angle in Table 1 to the whole number, as shown in Table 5.

Rotation axis | Angle (deg) | (q_{1}, q_{2}, q_{3}, q_{4}) | |
---|---|---|---|

1 | (0.273, 0.484, 0.831) | 128 | (0.245, 0.435, 0.747, 0.438) |

2 | (0.211, 0.454, 0.866) | 123 | (0.185, 0.397, 0.757, 0.485) |

3 | (0.160, 0.426, 0.891) | 115 | (0.134, 0.357, 0.747, 0.545) |

4 | (0.129, 0.400, 0.907) | 106 | (0.103, 0.320, 0.725, 0.602) |

5 | (0.122, 0.383, 0.916) | 96 | (0.090, 0.285, 0.680, 0.669) |

6 | (0.140, 0.381, 0.914) | 87 | (0.098, 0.265, 0.635, 0.719) |

7 | (0.185, 0.403, 0.896) | 78 | (0.116, 0.254, 0.564, 0.777) |

8 | (0.249, 0.460, 0.852) | 71 | (0.147, 0.271, 0.501, 0.809) |

9 | (0.322, 0.556, 0.767) | 66 | (0.175, 0.303, 0.418, 0.839) |

10 | (0.378, 0.671, 0.638) | 65 | (0.206, 0.366, 0.348, 0.839) |

11 | (0.400, 0.769, 0.499) | 69 | (0.224, 0.430, 0.279, 0.829) |

12 | (0.400, 0.825, 0.400) | 76 | (0.246, 0.508, 0.246, 0.788) |

Rotation axis | Angle (deg) | (q_{1}, q_{2}, q_{3}, q_{4}) | |
---|---|---|---|

1 | (0.273, 0.484, 0.831) | 128 | (0.245, 0.435, 0.747, 0.438) |

2 | (0.211, 0.454, 0.866) | 123 | (0.185, 0.397, 0.757, 0.485) |

3 | (0.160, 0.426, 0.891) | 115 | (0.134, 0.357, 0.747, 0.545) |

4 | (0.129, 0.400, 0.907) | 106 | (0.103, 0.320, 0.725, 0.602) |

5 | (0.122, 0.383, 0.916) | 96 | (0.090, 0.285, 0.680, 0.669) |

6 | (0.140, 0.381, 0.914) | 87 | (0.098, 0.265, 0.635, 0.719) |

7 | (0.185, 0.403, 0.896) | 78 | (0.116, 0.254, 0.564, 0.777) |

8 | (0.249, 0.460, 0.852) | 71 | (0.147, 0.271, 0.501, 0.809) |

9 | (0.322, 0.556, 0.767) | 66 | (0.175, 0.303, 0.418, 0.839) |

10 | (0.378, 0.671, 0.638) | 65 | (0.206, 0.366, 0.348, 0.839) |

11 | (0.400, 0.769, 0.499) | 69 | (0.224, 0.430, 0.279, 0.829) |

12 | (0.400, 0.825, 0.400) | 76 | (0.246, 0.508, 0.246, 0.788) |

*A*]:

The singular vectors that are associated with the five smallest singular values are listed in Table 6, forming the approximate null-space basis of [*A*]. The valid solutions of *α _{k}* in Table 7 that qualify to become the coefficients of C-manifold. In terms of the fitting error

*e*, all four solutions can be used to construct C-manifolds that fit the given data near perfectly. By applying the inverse computation given by Eq. (13) on the four resulting C-manifold vectors, the dimensions of the four resulting spherical circular constraints are listed in Table 8, and their corresponding C-manifolds are plotted in Fig. 6.

_{f}p_{1} | p_{2} | p_{3} | p_{4} | p_{5} | p_{6} | p_{7} | p_{8} | p_{9} | p_{10} | |
---|---|---|---|---|---|---|---|---|---|---|

v_{1} | −0.2746 | −0.0940 | 0.1309 | −0.3094 | −0.1319 | −0.0108 | 0.0317 | 0.0535 | −0.1062 | 0.8776 |

v_{2} | −0.3530 | −0.0014 | −0.0475 | 0.2533 | 0.7921 | 0.0959 | −0.2180 | −0.2065 | −0.2708 | 0.0937 |

v_{3} | −0.4655 | 0.0549 | −0.2255 | 0.6909 | −0.4450 | 0.1240 | 0.1659 | −0.0782 | −0.0152 | 0.0690 |

v_{4} | −0.0050 | 0.3005 | −0.4250 | −0.0202 | 0.0129 | −0.5885 | 0.0246 | 0.4217 | −0.4513 | 0.0004 |

v_{5} | −0.2053 | 0.1573 | 0.4045 | −0.2247 | −0.2718 | 0.2279 | −0.0435 | −0.1330 | −0.6957 | −0.2994 |

p_{1} | p_{2} | p_{3} | p_{4} | p_{5} | p_{6} | p_{7} | p_{8} | p_{9} | p_{10} | |
---|---|---|---|---|---|---|---|---|---|---|

v_{1} | −0.2746 | −0.0940 | 0.1309 | −0.3094 | −0.1319 | −0.0108 | 0.0317 | 0.0535 | −0.1062 | 0.8776 |

v_{2} | −0.3530 | −0.0014 | −0.0475 | 0.2533 | 0.7921 | 0.0959 | −0.2180 | −0.2065 | −0.2708 | 0.0937 |

v_{3} | −0.4655 | 0.0549 | −0.2255 | 0.6909 | −0.4450 | 0.1240 | 0.1659 | −0.0782 | −0.0152 | 0.0690 |

v_{4} | −0.0050 | 0.3005 | −0.4250 | −0.0202 | 0.0129 | −0.5885 | 0.0246 | 0.4217 | −0.4513 | 0.0004 |

v_{5} | −0.2053 | 0.1573 | 0.4045 | −0.2247 | −0.2718 | 0.2279 | −0.0435 | −0.1330 | −0.6957 | −0.2994 |

α_{1} | α_{2} | α_{3} | α_{4} | α_{5} | e_{s} | e_{f} | |
---|---|---|---|---|---|---|---|

Solution 1 | 1 | 2.7477 | −19.3633 | 50.9312 | −0.3269 | 9.86 × 10^{−13} | 1.45 × 10^{−6} |

Solution 2 | 1 | 4.9167 | −1.9715 | −2.7507 | −1.2318 | 2.01 × 10^{−15} | 1.99 × 10^{−4} |

Solution 3 | 1 | 1.8588 | −4.6311 | 3.1037 | 2.8302 | 1.28 × 10^{−15} | 5.64 × 10^{−4} |

Solution 4 | 1 | 1.2748 | −3.2178 | −2.1117 | −0.9790 | 2.86 × 10^{−16} | 8.45 × 10^{−4} |

α_{1} | α_{2} | α_{3} | α_{4} | α_{5} | e_{s} | e_{f} | |
---|---|---|---|---|---|---|---|

Solution 1 | 1 | 2.7477 | −19.3633 | 50.9312 | −0.3269 | 9.86 × 10^{−13} | 1.45 × 10^{−6} |

Solution 2 | 1 | 4.9167 | −1.9715 | −2.7507 | −1.2318 | 2.01 × 10^{−15} | 1.99 × 10^{−4} |

Solution 3 | 1 | 1.8588 | −4.6311 | 3.1037 | 2.8302 | 1.28 × 10^{−15} | 5.64 × 10^{−4} |

Solution 4 | 1 | 1.2748 | −3.2178 | −2.1117 | −0.9790 | 2.86 × 10^{−16} | 8.45 × 10^{−4} |

C-manifold | a:b:c:d | x | y | z | Angular link length (deg) |
---|---|---|---|---|---|

p_{1} | 0.0068:1:0.0052:0.2766 | −0.0135 | −0.5009 | −0.8654 | 73.9432 |

p_{2} | 1:0.5360:−0.5835:1.0359 | −0.5159 | 0.5621 | −0.6464 | 35.7145 |

p_{3} | 1:−4.5301:1.2694:1.7383 | 0.3434 | −0.4068 | −0.8465 | 68.8128 |

p_{4} | 1:−0.0143:0.0475:0.8749 | 0.0223 | 0.4754 | −0.8795 | 29.0937 |

C-manifold | a:b:c:d | x | y | z | Angular link length (deg) |
---|---|---|---|---|---|

p_{1} | 0.0068:1:0.0052:0.2766 | −0.0135 | −0.5009 | −0.8654 | 73.9432 |

p_{2} | 1:0.5360:−0.5835:1.0359 | −0.5159 | 0.5621 | −0.6464 | 35.7145 |

p_{3} | 1:−4.5301:1.2694:1.7383 | 0.3434 | −0.4068 | −0.8465 | 68.8128 |

p_{4} | 1:−0.0143:0.0475:0.8749 | 0.0223 | 0.4754 | −0.8795 | 29.0937 |

### Example 3: Application of Exact Spherical Motion Synthesis (Five Poses).

In this section, a practical application is used to illustrate that our approach can also be applied to exact motion synthesis of spherical four-bar mechanism, especially for the five-pose case. This application seeks to design a driving mechanism to enable the head of a cooling fan to move along the surface of a sphere in a three-dimensional pattern. Larochelle et al. [29] designed and patented a *infinity fan* driven by a spherical four-bar mechanism (see Fig. 8), which navigates the fan head through the five task poses as listed in Table 9.

Rotation axis | Angle (deg) | (q_{1}, q_{2}, q_{3}, q_{4}) | |
---|---|---|---|

1 | (−0.470, 0.032, −0.882) | 48.39 | (−0.193, 0.013, −0.362, 0.912) |

2 | (−0.194, 0.964, 0.180) | 82.23 | (−0.128, 0.634, 0.118, 0.753) |

3 | (−0.158, 0.841, −0.518) | 98.27 | (−0.119, 0.636, −0.392, 0.654) |

4 | (0.899, 0.436, −0.034) | 22.28 | (0.174, 0.084, −0.007, 0.981) |

5 | (−0.225, 0.938, −0.263) | 96.27 | (−0.167, 0.699, −0.196, 0.667) |

Rotation axis | Angle (deg) | (q_{1}, q_{2}, q_{3}, q_{4}) | |
---|---|---|---|

1 | (−0.470, 0.032, −0.882) | 48.39 | (−0.193, 0.013, −0.362, 0.912) |

2 | (−0.194, 0.964, 0.180) | 82.23 | (−0.128, 0.634, 0.118, 0.753) |

3 | (−0.158, 0.841, −0.518) | 98.27 | (−0.119, 0.636, −0.392, 0.654) |

4 | (0.899, 0.436, −0.034) | 22.28 | (0.174, 0.084, −0.007, 0.981) |

5 | (−0.225, 0.938, −0.263) | 96.27 | (−0.167, 0.699, −0.196, 0.667) |

*A*] of the size 5 × 10. Computing the singular values of [

*A*] yields

It could be observed that the first five singular values are equal to zero. Their associated right-singular vectors in Table 10 constitute the exact null-space basis of [*A*] and define the G-manifold vector by Eq. (19). Again, requiring additionally that *α _{k}* satisfy the four equations in Eq. (20) and excluding the solutions that would cause big structural error

*e*, we obtain the valid solutions of

_{s}*α*that are qualified for the coefficients of the C-manifold, and their corresponding fitting errors

_{k}*e*as well as structural errors

_{f}*e*are listed in Table 11. Based on the fitting error

_{s}*e*, all six solutions can be used to construct the C-manifolds that fit the given data perfectly. Therefore, we can apply the inverse computation Eq. (13) on the six resulting C-manifolds and get the dimensions of the six corresponding spherical RR dyads.

_{f}p_{1} | p_{2} | p_{3} | p_{4} | p_{5} | p_{6} | p_{7} | p_{8} | p_{9} | p_{10} | |
---|---|---|---|---|---|---|---|---|---|---|

v_{1} | 0.2234 | 0.0530 | 0.7581 | 0.2218 | −0.1996 | −0.0143 | 0.5229 | 0.0277 | −0.0960 | 0 |

v_{2} | 0.1339 | 0.6653 | −0.1295 | 0.5484 | 0.2319 | 0.0360 | −0.1414 | −0.1120 | −0.3665 | 0 |

v_{3} | 0.2947 | −0.3799 | −0.1256 | 0.3168 | −0.1259 | −0.0824 | −0.0298 | −0.7868 | 0.1008 | 0 |

v_{4} | 0.5844 | 0.1188 | 0.1790 | −0.2843 | −0.1944 | 0.0092 | −0.4809 | 0.0648 | −0.0157 | −0.5079 |

v_{5} | 0.2550 | −0.0341 | −0.1814 | −0.0869 | −0.5671 | 0.0034 | −0.1031 | 0.1491 | −0.3989 | 0.6163 |

p_{1} | p_{2} | p_{3} | p_{4} | p_{5} | p_{6} | p_{7} | p_{8} | p_{9} | p_{10} | |
---|---|---|---|---|---|---|---|---|---|---|

v_{1} | 0.2234 | 0.0530 | 0.7581 | 0.2218 | −0.1996 | −0.0143 | 0.5229 | 0.0277 | −0.0960 | 0 |

v_{2} | 0.1339 | 0.6653 | −0.1295 | 0.5484 | 0.2319 | 0.0360 | −0.1414 | −0.1120 | −0.3665 | 0 |

v_{3} | 0.2947 | −0.3799 | −0.1256 | 0.3168 | −0.1259 | −0.0824 | −0.0298 | −0.7868 | 0.1008 | 0 |

v_{4} | 0.5844 | 0.1188 | 0.1790 | −0.2843 | −0.1944 | 0.0092 | −0.4809 | 0.0648 | −0.0157 | −0.5079 |

v_{5} | 0.2550 | −0.0341 | −0.1814 | −0.0869 | −0.5671 | 0.0034 | −0.1031 | 0.1491 | −0.3989 | 0.6163 |

α_{1} | α_{2} | α_{3} | α_{4} | α_{5} | e_{f} | e_{s} | |
---|---|---|---|---|---|---|---|

Solution 1 | 1 | −5.7367 | −1.4498 | −1.5185 | 1.0993 | 2.19 × 10^{−8} | 8.11 × 10^{−17} |

Solution 2 | 1 | −1.5642 | 0.0819 | −0.0991 | −1.4850 | 2.12 × 10^{−8} | 1.31 × 10^{−17} |

Solution 3 | 1 | −4.8092 | −1.3928 | −1.8058 | 0.6413 | 2.17 × 10^{−8} | 4.30 × 10^{−17} |

Solution 4 | 1 | −0.4316 | −0.8223 | 0.7149 | −0.5728 | 2.22 × 10^{−8} | 9.02 × 10^{−17} |

Solution 5 | 1 | −0.5317 | 1.2974 | −0.4174 | 0.3263 | 2.21 × 10^{−8} | 6.32 × 10^{−17} |

Solution 6 | 1 | −1.2468 | −0.6894 | −1.2776 | 0.2009 | 2.14 × 10^{−8} | 1.02 × 10^{−16} |

α_{1} | α_{2} | α_{3} | α_{4} | α_{5} | e_{f} | e_{s} | |
---|---|---|---|---|---|---|---|

Solution 1 | 1 | −5.7367 | −1.4498 | −1.5185 | 1.0993 | 2.19 × 10^{−8} | 8.11 × 10^{−17} |

Solution 2 | 1 | −1.5642 | 0.0819 | −0.0991 | −1.4850 | 2.12 × 10^{−8} | 1.31 × 10^{−17} |

Solution 3 | 1 | −4.8092 | −1.3928 | −1.8058 | 0.6413 | 2.17 × 10^{−8} | 4.30 × 10^{−17} |

Solution 4 | 1 | −0.4316 | −0.8223 | 0.7149 | −0.5728 | 2.22 × 10^{−8} | 9.02 × 10^{−17} |

Solution 5 | 1 | −0.5317 | 1.2974 | −0.4174 | 0.3263 | 2.21 × 10^{−8} | 6.32 × 10^{−17} |

Solution 6 | 1 | −1.2468 | −0.6894 | −1.2776 | 0.2009 | 2.14 × 10^{−8} | 1.02 × 10^{−16} |

According to the results, we find that the third and fourth dyad solutions are very close to the dimensions of the two dyads that Larochelle et al. used to construct their infinity fan. Therefore, instead of listing the dimensions of our six dyads individually, we employ the third and fourth dyad to build a four-bar linkage, whose dimensions are compared against those of Larochelle's linkage given in Ref. [29]. Table 12 presents the comparison which shows excellent match.

Ours | Larochelle's | |
---|---|---|

Input link length | 47.3497 deg | 47.3502 deg |

Output link length | 53.8243 deg | 53.8062 deg |

Coupler link length | 88.3134 deg | 88.2839 deg |

Fixed link length | 87.8610 deg | 87.8943 deg |

Input fixed pivot (in frame F) | (0.7332, 0.0022, −0.6800) | (0.7327, 0.0000 −0.6805) |

Input moving pivot (in frame M) | (0.9257, 0.3757, 0.0445) | (0.9263, 0.3742, 0.0437) |

Output fixed pivot (in frame F) | (0.7033, −0.1042, 0.7032) | (0.7032, −0.1043, 0.7032) |

Output moving pivot (in frame M) | (0.2729, −0.6753, 0.6852) | (0.2727, −0.6750, 0.6856) |

Ours | Larochelle's | |
---|---|---|

Input link length | 47.3497 deg | 47.3502 deg |

Output link length | 53.8243 deg | 53.8062 deg |

Coupler link length | 88.3134 deg | 88.2839 deg |

Fixed link length | 87.8610 deg | 87.8943 deg |

Input fixed pivot (in frame F) | (0.7332, 0.0022, −0.6800) | (0.7327, 0.0000 −0.6805) |

Input moving pivot (in frame M) | (0.9257, 0.3757, 0.0445) | (0.9263, 0.3742, 0.0437) |

Output fixed pivot (in frame F) | (0.7033, −0.1042, 0.7032) | (0.7032, −0.1043, 0.7032) |

Output moving pivot (in frame M) | (0.2729, −0.6753, 0.6852) | (0.2727, −0.6750, 0.6856) |

## Conclusions

In this paper, we seek to extend our previous work for planar four-bar synthesis case to spherical four-bar linkage synthesis. Our approach can be distinctly divided into two steps. By formulating the data fitting process in a linear form via kinematic mapping, the first step is finding a pencil of general quadratic manifolds, called G-manifolds, in the image space that best fit the task image points in the least-squares sense, which is done by using singular value decomposition and solving for the singular vectors. The singular vectors associated with the smallest singular values are then linearly combined to define the null-space solution of a pencil of G-manifolds. Second, four additional constraints on the components of G-manifold vectors are then imposed to identify the C-manifolds from the pencil of G-manifolds that are qualified for representing a spherical circular constraint. After the inverse computation converting the C-manifold vectors to the design parameters of the spherical RR dyad, spherical four-bar linkages that best guide through the set of given poses can be constructed with the obtained spherical RR dyads as building blocks. The resulting approach for spherical four-bar linkage synthesis is simple and efficient, which leads to a unified algorithm for both exact and approximate synthesis of spherical dyads.

## Acknowledgment

All findings and results presented in this paper are those of the authors and do not represent those of the funding agencies. The authors would also like to thank Kartik Thakkar for the contribution of Mathematica codes implementation work in Example 2.

## Funding Data

National Natural Science Foundation of China (Grant Nos. 51405128 and 51775155).

Division of Civil, Mechanical and Manufacturing Innovation, U.S. National Science Foundation (Grant No. CMMI-1563413).

Fundamental Research Funds for the Central Universities of China (Grant No. 2682017CX037).