Synthesis of Four-Link Initial Kinematic Chains with Spherical Pairs for Spatial Mechanisms

1. Introduction

The synthesis of kinematic chains plays a crucial role in the design of spatial mechanisms, which are widely used in industries such as robotics, aerospace, and automotive engineering [1]. One particular type of spatial mechanism employs Spherical-Spherical-Spherical (SSS) kinematic pairs, which allow for three degrees of rotational freedom. These pairs are vital in creating complex motion paths and orientations, making them ideal for applications that require high precision and flexibility in movement. In kinematic chain synthesis, the goal is to design a mechanism that meets specific motion requirements while adhering to geometric constraints. For spatial mechanisms with SSS pairs, this involves determining the optimal configuration of the links and joints that will ensure the desired motion trajectory and orientation [2]. However, the synthesis of such mechanisms is challenging due to the complexity of the spatial relationships and the nonlinear nature of the equations governing their motion [3]. The problem of synthesizing kinematic chains with SSS pairs can be divided into two primary tasks: (1) defining the geometry of the links and joints, and (2) solving the kinematic equations to ensure the mechanism achieves the required motion. This process often involves solving systems of nonlinear equations, which describe the spatial positioning of the links relative to one another. Furthermore, the synthesis must take into account various constraints, such as the range of motion, joint limits, and load-bearing capacity, all of which can affect the mechanism’s performance [4]. This research aims to address the synthesis problem by providing a systematic approach to designing kinematic chains with SSS pairs. We focus on developing methods to solve the initial synthesis problem, ensuring that the mechanism's geometry and motion capabilities align with the desired performance criteria [5]. By leveraging both analytical and computational techniques, we seek to establish reliable solutions for these complex spatial mechanisms, offering insights into their design and optimization. In the following sections, we review existing methods for kinematic chain synthesis [6], propose new techniques for addressing the challenges specific to SSS pairs, and demonstrate the practical applications of our approach through case studies. This work contributes to the advancement of spatial mechanism design, providing tools for engineers to create more efficient and precise systems in various fields. The rapid advancement of robotics, computer graphics, and mechanical design has underscored the importance of accurate and intuitive visualization of geometric transformations and kinematic relationships [7]. In multi-body systems and robotic mechanisms, the ability to effectively analyze and debug complex spatial interactions is paramount. This study presents a comprehensive framework for three-dimensional (3D) visualization of coordinate transformations, enabling researchers and practitioners to gain insights into the behavior of kinematic chains and transformation matrices [8]. A central challenge in kinematic analysis is the representation of the spatial relationships between different coordinate frames. These relationships are often described by transformation matrices whose components depend on rotational parameters such as Euler angles θ_i, ψ_i, $\phi$_i. Understanding how these angles influence the resultant vectors and the overall motion of the system is critical for tasks such as inverse kinematics, dynamic simulation [9], and control. Our approach leverages Python’s powerful numerical and plotting libraries-NumPy and Matplotlib-to generate detailed 2D and 3D visualizations that depict the relationships between fixed points X_A, Y_A, Z_A; x_B, y_B, z_B; x_C, y_C, z_C; and R and their corresponding transformation parameters [10]. In this work, we not only demonstrate the plotting of individual points but also explore the visualization of derived quantities such as the differences between coordinates X_Di−X_A, Y_Di−Y_A, Z_Di−Z_A, as well as the components of transformation matrices (e, m, n) computed from rotational parameters. Moreover, we illustrate the visualization of composite transformations obtained through matrix multiplication $T_{21}^i=T_{01}^i\times T_{20}^i$ and $T_{12}^i=T_{02}^i\times T_{10}^i$. These visualizations provide a powerful tool for interpreting the intricate geometric relationships inherent in kinematic chain synthesis and robotic system design. By integrating these visualization techniques into a cohesive framework, our methodology enables users to interactively analyze and validate the performance of complex spatial transformations [11]. This work aims to contribute to the field by offering a flexible and accessible approach to 3D visualization, ultimately enhancing the understanding and debugging of multi-dimensional transformation processes in various engineering applications. In the sections that follow, we detail the mathematical foundations underlying the transformation matrices, describe the implementation of the visualization routines in Python, and present several case studies that demonstrate the practical utility of our approach.

2. Materials and Methods

The problem is describing involves finding specific points on two solids such that the distance between those points in all given positions remains as close as possible to a constant value R. To approach this, we have formulated a weighted difference function that will be minimized over the given positions of the solids [12]. Here is a breakdown of the solution process:

Given the position N of the two solids Q₁ and Q₂, we have the Eulerian angles for solid Q₁ and coordinates X_Di, Y_Di, Z_Di and angles for solid Q₂. We need to find points A(X_A, Y_A, Z_A) on solid Q₁, B(x_B, y_B, z_B) on solid Q₁, and C(x_C, y_C, z_C) on solid Q₂, such that the distance |B_iC_i| is close to a constant R for all positions as shown in Figure 1. Minimize the difference ∆qi for each position, which is expressed as:

$${\Delta }_{q_i}={\left|\overrightarrow{B_i{C}_i}\right|}^2-R^2=(X_{C_i}-X_{B_i}{)}^2+(Y_{C_i}-Y_{B_i}{)}^2+(Z_{C_i}-Z_{B_i}{)}^2-R^2,\left(i=\overrightarrow{1,N}\right)$$

(1)

where x_Bi, y_Bi, z_Bi are the coordinates of point B in the position i on solid Q₁, and x_Ci, y_Ci, z_Ci are the coordinates of point C in the position i on solid Q₂,

The matrix $T_{jo}^i$ that is composed of directional cosines (or transformation matrices) that relate the local coordinate systems of each solid Q₁ and Q₂ to a fixed coordinate system. The matrix itself is built from the Euler angles ${\theta }_i^j$, ${\psi }_i^j$, ${\phi }_i^j$ of each solid in each position i, where j=1,2 refers to the two solids [13]. The components of the matrix are structured as follows:

$$T_{j0}^i=\left[\begin{array}{ccc}{e}_{i1}^j& e_{i2}^j& e_{i3}^j\\ m_{i1}^j& m_{i2}^j& m_{i3}^j\\ n_{i1}^j& n_{i2}^j& n_{i3}^j\end{array}\right],(j=\overline{\mathrm{1,2}};i=\overline{1,N})$$

(2)

where $e_{i1}^j$, $e_{i2}^j$, $e_{i3}^j$ is the components of the first row (often representing the direction cosines of the local x-axis in the fixed coordinate system); $m_{i1}^j$, $m_{i2}^j$, $m_{i3}^j$ is the components of the second row (typically representing the direction cosines of the local y-axis in the fixed coordinate system); $n_{i1}^j$, $n_{i2}^j$, $n_{i3}^j$, is the represent the components of the third row (usually the direction cosines of the local z-axis in the fixed coordinate system).

These represent the direction cosines of the local x-axis relative to the fixed axes:

$$\begin{gathered} e_{i1}^j=\mathit{cos}{\psi }_i^j\cdot \mathit{cos}{\phi }_i^j-\mathit{cos}{\theta }_i^j\cdot \mathit{sin}{\psi }_i^j\cdot \mathit{sin}{\phi }_i^j \\ {e}_{i2}^j=-\mathit{cos}{\psi }_i^j\cdot \mathit{sin}{\phi }_i^j-\mathit{cos}{\theta }_i^j\cdot \mathit{sin}{\psi }_i^j\cdot \mathit{cos}{\phi }_i^j \\ {e}_{i3}^j=\mathit{sin}{\theta }_i^j\cdot \mathit{sin}{\psi }_i^j \end{gathered}$$

(3)

These represent the direction cosines of the local y-axis relative to the fixed axes:

$$\begin{gathered} m_{i1}=\mathit{sin}{\psi }_i^j\cdot \mathit{cos}{\phi }_i^j+\mathit{cos}{\theta }_i^j\cdot \mathit{cos}{\psi }_i^j\cdot \mathit{cos}{\phi }_i^j \\ {m}_{i2}^j=-\mathit{sin}{\psi }_i^j\cdot \mathit{sin}{\phi }_i^j+\mathit{cos}{\theta }_i^j\cdot \mathit{cos}{\psi }_i^j\cdot \mathit{sin}{\phi }_i^j \\ {m}_{i3}^j=-\mathit{sin}{\theta }_i^j\cdot \mathit{cos}{\psi }_i^j \end{gathered}$$

(4)

These represent the direction cosines of the local z-axis relative to the fixed axes:

$$\begin{gathered} n_{i1}^j=\mathit{sin}{\theta }_i^j\cdot \mathit{sin}{\phi }_i^j \\ {n}_{i2}^j=\mathit{sin}{\theta }_i^j\cdot \mathit{cos}{\phi }_i^j \\ {n}_{i3}^j=\mathit{cos}{\phi }_i^j \end{gathered}$$

(5)

The matrix $T_{jo}^i$ essentially transforms the orientation of solid Q₁ or Q₂ in a fixed coordinate system. Each solid's orientation is represented using the Euler angles (or rotation matrix) to map the local axes to the global (fixed) coordinate system [14]. The three rows of the matrix represent the direction cosines of the local axes relative to the global axes. This matrix is likely part of a kinematic analysis where the positions and orientations of two solids need to be related to a common coordinate system [15]. The task is then to use this matrix to calculate the positions of points on the solids, transforming them between the local and global frames. The matrices $T_{jo}^i$ and their transposed versions represent the rotation matrices that transform the coordinate system of one solid to another. Let's break down each part of equations and matrices:

$T_{01}^i={\left[T_{10}^i\right]}^T$ and $T_{02}^i={\left[T_{20}^i\right]}^T$ are the transposed versions of $T_{10}^i$ and $T_{20}^i$, respectively. These matrices represent the rotation transformation from the local coordinate frame of solid Q₁ and Q₂ to a fixed global coordinate frame at each position i. Each of these matrices contains the direction cosines for the three local axes of the solids. In matrix form, we have:

$$T_{01}^i={\left[T_{10}^i\right]}^T=\left[\begin{array}{ccc}{e}_{i1}^1& m_{i1}^1& n_{i1}^1\\ e_{i2}^1& m_{i2}^1& n_{i2}^1\\ e_{i3}^1& m_{i3}^1& n_{i3}^1\end{array}\right]$$

(6)

$$T_{02}^i={\left[T_{20}^i\right]}^T=\left[\begin{array}{ccc}{e}_{i1}^2& m_{i1}^2& n_{i1}^2\\ e_{i2}^2& m_{i2}^2& n_{i2}^2\\ e_{i3}^2& m_{i3}^2& n_{i3}^2\end{array}\right]$$

(7)

where the components $e_{i1}$, $m_{i1}$, $n_{i1}$ etc., are the direction cosines of the respective axes relative to the global coordinate frame.

The matrix $T_{21}^i$ represents the transformation of the coordinate frame of Q₂ relative to Q₁. It is the result of multiplying the matrices $T_{01}^i$ and $T_{20}^i$. The components of the resulting matrix are:

$$T_{21}^i=\left[\begin{array}{ccc}{e}_{i1}^1{e}_{i1}^2+e_{i2}^1{m}_{i1}^2+e_{i3}^1{n}_{i1}^2& e_{i1}^1{e}_{i2}^2+e_{i2}^1{m}_{i2}^2+e_{i3}^1{n}_{i2}^2& e_{i1}^1{e}_{i3}^2+e_{i2}^1{m}_{i3}^2+e_{i3}^1{n}_{i3}^2\\ m_{i1}^1{e}_{i1}^2+m_{i2}^1{m}_{i1}^2+m_{i3}^1{n}_{i1}^2& m_{i1}^1{e}_{i2}^2+m_{i2}^1{m}_{i2}^2+m_{i3}^1{n}_{i2}^2& m_{i1}^1{e}_{i3}^2+m_{i2}^1{m}_{i3}^2+m_{i3}^1{n}_{i3}^2\\ n_{i1}^1{e}_{i1}^2+n_{i2}^1{m}_{i1}^2+n_{i3}^1{n}_{i1}^2& n_{i1}^1{e}_{i2}^2+n_{i2}^1{m}_{i2}^2+n_{i3}^1{n}_{i2}^2& n_{i1}^1{e}_{i3}^2+n_{i2}^1{m}_{i3}^2+n_{i3}^1{n}_{i3}^2\end{array}\right]$$

(8)

This matrix is formed by performing the dot products of the rows of $T_{01}^i$ with the columns of $T_{02}^i$. The result is a 3x3 matrix that represents the relative orientation between Q₂ and Q₁.

Similarly, the matrix $T_{12}^i$ represents the transformation of the coordinate frame of Q₁ relative to Q₂, and is the result of multiplying $T_{02}^i$ and $T_{10}^i$. The components of $T_{12}^i$ are:

$$T_{12}^i=\left[\begin{array}{ccc}{e}_{i1}^2{e}_{i1}^1+m_{i1}^2{m}_{i1}^1+n_{i1}^2{n}_{i1}^1& e_{i1}^2{e}_{i2}^1+m_{i1}^2{m}_{i2}^1+n_{i1}^2{n}_{i2}^1& e_{i1}^2{e}_{i3}^1+m_{i1}^2{m}_{i3}^1+n_{i1}^2{n}_{i3}^1\\ e_{i2}^2{e}_{i1}^1+m_{i2}^2{m}_{i1}^1+n_{i2}^2{n}_{i1}^1& m_{i2}^2{e}_{i2}^1+m_{i2}^2{m}_{i2}^1+m_{i2}^2{n}_{i2}^1& m_{i2}^2{e}_{i3}^1+m_{i2}^2{m}_{i3}^1+m_{i2}^2{n}_{i3}^1\\ n_{i3}^2{e}_{i1}^1+n_{i3}^2{m}_{i1}^1+n_{i3}^2{n}_{i1}^1& n_{i3}^2{e}_{i2}^1+n_{i3}^2{m}_{i2}^1+n_{i3}^2{n}_{i2}^1& n_{i3}^2{e}_{i3}^1+n_{i3}^2{m}_{i3}^1+n_{i3}^2{n}_{i3}^1\end{array}\right]$$

(9)

$T_{01}^i$ and $T_{02}^i$ represent the rotation matrices (transposed) from each solid's local frame to the fixed coordinate system. $T_{12}^i$ and $T_{21}^i$ represent the relative transformations between the two solids' local frames, where $T_{21}^i$ transforms Q₂ local frame to Q₁ local frame and $T_{12}^i$ does the reverse. These matrices are essential in solving for the positions and distances between points on the two solids, as they allow you to transform coordinates between different reference frames and to express the relative motion or position between the solids. The matrix system performs transformations between different coordinate frames using rotation matrices and vector translations [16], which are relative to the positions of points on the rigid bodies Q₁ and Q₂ in the fixed coordinate system. Here is a breakdown of the matrices and their transformations:

The matrix equation for point A is:

$$\left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\\ 1\end{array}\right]=\left[\begin{array}{cccc}{e}_{i1}^1& e_{i2}^1& e_{i3}^1& 0\\ m_{i1}^1& m_{i2}^1& m_{i3}^1& 0\\ n_{i1}^1& n_{i2}^1& n_{i3}^1& 0\\ 0& 0& 0& 1\end{array}\right]\times \left[\begin{array}{c}{x}_B\\ y_B\\ z_B\\ 1\end{array}\right]+\left[\begin{array}{cccc}{e}_{i1}^2& e_{i2}^2& e_{i3}^2& X_{Di}\\ m_{i1}^2& m_{i2}^2& m_{i3}^2& Y_{Di}\\ n_{i1}^2& n_{i2}^2& n_{i3}^2& Z_{Di}\\ 0& 0& 0& 1\end{array}\right]\times \left[\begin{array}{c}{x}_C\\ y_C\\ z_C\\ 1\end{array}\right]$$

(10)

This equation expresses the position of point A as the result of first applying a rotation defined by $T_{01}^i$ (for solid Q₁) and then translating by the coordinates of point D on solid Q₂.

For point B, the matrix equation is:

$$\left[\begin{array}{c}{x}_B\\ y_B\\ z_B\\ 1\end{array}\right]=\left[\begin{array}{cccc}{e}_{i1}^1& e_{i2}^1& e_{i3}^1& 0\\ m_{i1}^1& m_{i2}^1& m_{i3}^1& 0\\ n_{i1}^1& n_{i2}^1& n_{i3}^1& 0\\ 0& 0& 0& 1\end{array}\right]\times \left[\begin{array}{c}{X}_{Di}-X_A\\ Y_{Di}-Y_A\\ Z_{Di}-Z_A\\ 1\end{array}\right]+\left[rotationmatrixterms\right]\times \left[\begin{array}{c}{x}_C\\ y_C\\ z_C\\ 1\end{array}\right]$$

(11)

This expression represents a combined rotation and translation for point B using the relative position of D and the transformation from Q₂ to Q₁.

The equation for point C is:

$$\left[\begin{array}{c}{x}_C\\ y_C\\ z_C\\ 1\end{array}\right]=\left[\begin{array}{cccc}{e}_{i1}^2& e_{i2}^2& e_{i3}^2& 0\\ m_{i1}^2& m_{i2}^2& m_{i3}^2& 0\\ n_{i1}^2& n_{i2}^2& n_{i3}^2& 0\\ 0& 0& 0& 1\end{array}\right]\times \left[\begin{array}{c}{X_A-X}_{Di}\\ Y_A{-Y}_{Di}\\ Z_A{-Z}_{Di}\\ 1\end{array}\right]+\left[rotationmatrixterms\right]\times \left[\begin{array}{c}{x}_B\\ y_B\\ z_B\\ 1\end{array}\right]$$

(12)

This represents a similar transformation for point C, which involves both a translation and rotation between the frames of Q₁ and Q₂. Each of these matrices is a combination of rotation matrices (representing the orientation of the rigid bodies Q₁ and Q₂ relative to the global coordinate frame) and translation matrices (replacing the initial position of the rigid bodies). The purpose of these transformations is to calculate the relative positions of points A, B, and C on the rigid bodies Q₁ and Q₂ based on the coordinate transformations [17]. These transformations take into account both the orientation and position of the rigid bodies in space, using Euler angles for rotations and coordinates for translations.

So, for the weighted differences ${∆q}_i^{\left(1\right)}$, ${∆q}_i^{\left(2\right)}$, ${∆q}_i^{\left(3\right)}$ as functions of the parameters X_A, Y_A, Z_A we divide the objective function into three separate forms x_B, y_B, z_B, x_C, y_C, z_C, and R.

Weighted difference for point A:

$${∆q}_i^{\left(1\right)}={\left(X_{Ai}-X_A\right)}^2+{\left(Y_{Ai}-Y_A\right)}^2+{\left(Z_{Ai}-Z_A\right)}^2-R^2$$

(13)

Weighted difference for point B:

$${∆q}_i^{\left(2\right)}={\left(X_{Bi}-X_B\right)}^2+{\left(Y_{Bi}-Y_B\right)}^2+{\left(Z_{Bi}-Z_B\right)}^2-R^2$$

(14)

Weighted difference for point C:

$${∆q}_i^{\left(3\right)}={\left(X_{Ci}-X_C\right)}^2+{\left(Y_{Ci}-Y_C\right)}^2+{\left(Z_{Ci}-Z_C\right)}^2-R^2$$

(15)

The total objective function to minimize, which is the sum of the squared differences, can be written as:

$$S={\sum }_{i=1}^N\left[{\left({∆q}_i^{\left(1\right)}\right)}^2+{\left({∆q}_i^{\left(2\right)}\right)}^2+{\left({∆q}_i^{\left(3\right)}\right)}^2\right]$$

(16)

The sum here is taken over all N discrete positions of the solids. Since each ${∆q}_i^{\left(k\right)}$ (for k=1,2,3) is expressed in terms of X_A, Y_A, Z_A, x_B, y_B, z_B, x_C, y_C, z_C, and R, the objective function S is a function of these ten parameters. We represent the weighted differences as groups of parameters with a common parameter R into four sets. For each of the three members, the variables X_A, Y_A, Z_A (for point A), x_B, y_B, z_B (for point B), and x_C, y_C, z_C (for point C) are grouped together and their corresponding coordinates X_A, Y_A, Z_A, x_B, y_B, z_B, and x_C, y_C, z_C are scaled relative to the common radius R. To find the optimal positions and radius that minimize S, you need to compute the partial derivatives of S with respect to each parameter (i.e., the positions and R), and set these derivatives equal to zero. The total objective function is:

$$S={\sum }_{i=1}^N\left[{\left({\left(X_{Ai}-X_A\right)}^2+{\left(Y_{Ai}-Y_A\right)}^2+{\left(Z_{Ai}-Z_A\right)}^2-R^2\right)}^2+{\left({\left(X_{Bi}-X_B\right)}^2+{\left(Y_{Bi}-Y_B\right)}^2+{\left(Z_{Bi}-Z_B\right)}^2-R^2\right)}^2+{\left({\left(X_{Ci}-X_C\right)}^2+{\left(Y_{Ci}-Y_C\right)}^2+{\left(Z_{Ci}-Z_C\right)}^2-R^2\right)}^2\right]$$

(17)

Now, to minimize S, we calculate the partial derivatives of S with respect to each of the ten parameters X_A, Y_A, Z_A, x_B, y_B, z_B, x_C, y_C, z_C, and R, and set them equal to zero:

$$\begin{gathered} {\displaystyle \frac{\partial \mathrm{S}}{\partial {\mathrm{X}}_{\mathrm{A}}}}=0;{\displaystyle \frac{\partial \mathrm{S}}{\partial {\mathrm{Y}}_{\mathrm{A}}}}=0; {\displaystyle \frac{\partial \mathrm{S}}{\partial {\mathrm{Z}}_{\mathrm{A}}}}=0; {\displaystyle \frac{\partial \mathrm{S}}{\partial \mathrm{R}}}=0; \\ {\displaystyle \frac{\partial \mathrm{S}}{\partial {\mathrm{x}}_{\mathrm{B}}}}=0; {\displaystyle \frac{\partial \mathrm{S}}{\partial {\mathrm{y}}_{\mathrm{B}}}}=0; {\displaystyle \frac{\partial \mathrm{S}}{\partial {\mathrm{z}}_{\mathrm{B}}}}=0; {\displaystyle \frac{\partial \mathrm{S}}{\partial \mathrm{R}}}=0; \\ {\displaystyle \frac{\partial \mathrm{S}}{\partial {\mathrm{x}}_{\mathrm{C}}}}=0; {\displaystyle \frac{\partial \mathrm{S}}{\partial {\mathrm{y}}_{\mathrm{C}}}}=0; {\displaystyle \frac{\partial \mathrm{S}}{\partial {\mathrm{z}}_{\mathrm{C}}}}=0; {\displaystyle \frac{\partial \mathrm{S}}{\partial \mathrm{R}}}=0 \end{gathered}$$

(18)

To solve the system of equations formed by setting these partial derivatives to zero, you need to:

we expand each of the partial derivatives, differentiating the sum of the square terms for each point with respect to the corresponding parameter;

to simplify the resulting expressions, after differentiation, we obtain a set of linear equations for each parameter (the ten parameters being optimized);

the system of equations can be solved numerically using methods such as Gaussian elimination, least squares, or gradient descent optimization algorithms.

We are now working with a set of equations with weighted differences ${∆q}_i^{\left(1\right)}$ for each i. These equations aim to find the optimal values of X_A, Y_A, Z_A, R using the following system:

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}\left(X_{Ai}-X_A\right)=0$$

(19)

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}\left(Y_{Ai}-Y_A\right)=0$$

(20)

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}\left(Z_{Ai}-Z_A\right)=0$$

(21)

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}R=0$$

(22)

These equations are essentially normal equations that come from minimizing the sum of squared weighted differences S. Each equation corresponds to a specific parameter (either X_A, Y_A, Z_A, or R) that you are trying to optimize [18]. Each equation represents the sum of weighted differences between the coordinates X_A, Y_A, Z_A and their respective target values X_Ai, Y_Ai, Z_Ai, weighted by the difference function ${∆q}_i^{\left(1\right)}$. The goal is to adjust the coordinates so that the sum of these weighted differences is minimized (i.e., set each sum to zero).

This equation sums the weighted differences over all positions. Since ${∆q}_i^{\left(1\right)}$ is a function of X_A, Y_A, Z_A, and R, this equation ensures that the weighted differences balance out when summed across all positions:

$$\sum _{i=1}^N{∆}_{qi}^{\left(1\right)}=0$$

These equations ensure that the weighted differences ${∆q}_i^{\left(1\right)}$ are balanced with respect to the coordinates X_A, Y_A, Z_A. Each equation adjusts one of the coordinates of point A to minimize the sum of the squared weighted differences.

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}{X}_{Ai}=0$$

(23)

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}{Y}_{Ai}=0$$

(24)

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}{Z}_{Ai}=0$$

(25)

$$\sum _{i=1}^N{∆}_{qi}^{\left(1\right)}R=0$$

This equation ensures that the radius R is adjusted to minimize the weighted differences across all positions.

We propose a system of equations involving sums of square terms and cross products for the measured differences. The goal is to simplify these conditions to find the optimal values of the parameters. Let's consider the structure of the equations step by step. This equation simplifies terms involving the coordinates X_A, Y_A, Z_A, and the radius R, and is related to the weighted sum of squared distances:

$${\sum }_{i=1}^N\left[X_{Ai}^2{X}_A+X_{Ai}{Y}_{Ai}{Y}_A+Z_{Ai}{X}_{Ai}{Z}_A+{\displaystyle \frac{1}{2}}(R^2-X_A^2-Y_A^2-Z_A^2)X_{Ai}\right]={\displaystyle \frac{1}{2}}{\sum }_{i=1}^N(X_{Ai}^2+Y_{Ai}^2+Z_{Ai}^2)X_{Ai}$$

(26)

This equation represents a similar relationship for the Y_A coordinates:

$${\sum }_{i=1}^N\left[X_{Ai}{Y}_{Ai}{X}_A+Y_{Ai}^2{Y}_A+Z_{Ai}{Y}_{Ai}{Z}_A+{\displaystyle \frac{1}{2}}(R^2-X_A^2-Y_A^2-Z_A^2)Y_{Ai}\right]={\displaystyle \frac{1}{2}}{\sum }_{i=1}^N(X_{Ai}^2+Y_{Ai}^2+Z_{Ai}^2)Y_{Ai}$$

(27)

This equation applies the same principle to the Z_A coordinates:

$${\sum }_{i=1}^N\left[Z_{Ai}{X}_{Ai}{X}_A+Y_{Ai}{Z}_{Ai}{Y}_A+Z_{Ai}^2{Z}_A+{\displaystyle \frac{1}{2}}(R^2-X_A^2-Y_A^2-Z_A^2)Z_{Ai}\right]={\displaystyle \frac{1}{2}}{\sum }_{i=1}^N(X_{Ai}^2+Y_{Ai}^2+Z_{Ai}^2)Z_{Ai}$$

(28)

This equation represents the relationship between the coordinates and the sum of squared values:

$${\sum }_{i=1}^N\left[X_{Ai}{X}_A+Y_{Ai}{Y}_A+Z_{Ai}{Z}_A+{\displaystyle \frac{1}{2}}(R^2-X_A^2-Y_A^2-Z_A^2)\right]={\displaystyle \frac{1}{2}}{\sum }_{i=1}^N(X_{Ai}^2+Y_{Ai}^2+Z_{Ai}^2)$$

(29)

This equation expresses a term H₁ in terms of the radius R and coordinates X_A, Y_A, Z_A.

$$H_1={\displaystyle \frac{1}{2}}(R^2-X_A^2-Y_A^2-Z_A^2)$$

(30)

These equations are designed to relate the parameters X_A, Y_A, Z_A, R by simplifying the weighted sums over the N positions of the solids. To solve them, you'll typically:

we substitute the expressions for each term and simplify the sum at position N;

we solve the resulting system of linear equations to find the values of X_A, Y_A, Z_A, and R.

Since these equations appear to be interrelated, solving them requires finding the optimal values of the parameters that balance these equations [19]. This matrix contains the sums of the product terms for X_A, Y_A, Z_A, as well as their sums individually. These sums correspond to the moments of the data (coordinates X_A, Y_A, Z_A) over the N positions. The system of equations is given by:

$$\left[\begin{array}{cccc}{\sum }_{i=1}^N{X}_{Ai}^2& {\sum }_{i=1}^N{X}_{Ai}{Y}_{Ai}& {\sum }_{i=1}^N{X}_{Ai}{Z}_{Ai}& {\sum }_{i=1}^N{X}_{Ai}\\ {\sum }_{i=1}^N{X}_{Ai}{Y}_{Ai}& {\sum }_{i=1}^N{Y}_{Ai}& {\sum }_{i=1}^N{Y}_{Ai}{Z}_{Ai}& {\sum }_{i=1}^N{Y}_{Ai}\\ {\sum }_{i=1}^N{X}_{Ai}{Z}_{Ai}& {\sum }_{i=1}^N{Y}_{Ai}{Z}_{Ai}& {\sum }_{i=1}^N{Z}_{Ai}^2& {\sum }_{i=1}^N{Z}_{Ai}\\ {\sum }_{i=1}^N{X}_{Ai}& {\sum }_{i=1}^N{Y}_{Ai}& {\sum }_{i=1}^N{Z}_{Ai}& N\end{array}\right]\times \left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\\ H_1\end{array}\right]={\displaystyle \frac{1}{2}}\left[\begin{array}{c}{\sum }_{i=1}^N{R}_{A_i}^2{X}_{A_i}\\ {\sum }_{i=1}^N{R}_{A_i}^2{Y}_{A_i}\\ {\sum }_{i=1}^N{R}_{A_i}^2{Z}_{A_i}\\ {\sum }_{i=1}^N{R}_{A_i}^2\end{array}\right]$$

(31)

Using Cramer's Rule, the solution for the unknowns X_A, Y_A, Z_A, H₁ is given by:

$$(X_A,Y_A,Z_A,H_1)={\displaystyle \frac{1}{D_1}}(D_{X_A},D_{Y_A},D_{Z_A},D_{H_1}), (\mathrm{at} {\mathrm{D}}_1\ne 0)$$

(32)

where D₁ is the determinant of the coefficient matrix on the left-hand side. $D_{X_A},D_{Y_A},D_{Z_A}, D_{H_1}$ are the determinants of the matrices obtained by replacing the corresponding columns of the coefficient matrix with the right-hand side vector (i.e., the values for X_A, Y_A, Z_A, H₁).

$R_{Ai}^2$ is the squared distance from the origin for each point A, which is calculated as:

$$R_{Ai}^2=X_{Ai}^2+Y_{Ai}^2+Z_{Ai}^2$$

(33)

The goal is to solve the system of linear equations for X_A, Y_A, Z_A, H₁. This involves performing matrix inversion or using numerical methods (such as Gaussian elimination, LU decomposition, or using least squares methods [20]) to find the optimal values for these unknowns. Now, we have another system to solve for the parameters x_B, y_B, z_B, H₂ which is structured similarly to the previous one:

$$\left[\begin{array}{cccc}{\sum }_{i=1}^N{x}_{Bi}^2& {\sum }_{i=1}^N{x}_{Bi}{y}_{Bi}& {\sum }_{i=1}^N{x}_{Bi}{z}_{Bi}& {\sum }_{i=1}^N{x}_{Bi}\\ {\sum }_{i=1}^N{x}_{Bi}{y}_{Bi}& {\sum }_{i=1}^N{y}_{Bi}& {\sum }_{i=1}^N{y}_{Bi}{z}_{Bi}& {\sum }_{i=1}^N{y}_{Bi}\\ {\sum }_{i=1}^N{x}_{Bi}{z}_{Bi}& {\sum }_{i=1}^N{y}_{Bi}{z}_{Bi}& {\sum }_{i=1}^N{z}_{Bi}^2& {\sum }_{i=1}^N{z}_{Bi}\\ {\sum }_{i=1}^N{x}_{Bi}& {\sum }_{i=1}^N{y}_{Bi}& {\sum }_{i=1}^N{z}_{Bi}& N\end{array}\right]\times \left[\begin{array}{c}{x}_B\\ y_B\\ z_B\\ H_2\end{array}\right]={\displaystyle \frac{1}{2}}\left[\begin{array}{c}{\sum }_{i=1}^N{R}_{B_i}^2{x}_{B_i}\\ {\sum }_{i=1}^N{R}_{B_i}^2{y}_{B_i}\\ {\sum }_{i=1}^N{R}_{B_i}^2{z}_{B_i}\\ {\sum }_{i=1}^N{R}_{B_i}^2\end{array}\right]$$

(34)

Again using Cramer's Rule, the solution for the unknowns x_B, y_B, z_B, H₂ is given by:

$$(x_B,y_B,z_B,H_2)={\displaystyle \frac{1}{D_2}}(D_{x_B},D_{y_B},D_{z_B},D_{H_1}), (\mathrm{at} {\mathrm{D}}_2\ne 0)$$

(35)

where D₂ is the determinant of the coefficient matrix on the left-hand side. $D_{X_B},D_{Y_B},D_{Z_B}, D_{H_2}$ are the determinants of the matrices obtained by replacing the corresponding columns with the right-hand side vector.

We again applying Cramer's Rule to solve the system of equations, this time for the parameters x_C, y_C, z_C, H₃. Let's break this down:

$$\left[\begin{array}{cccc}{\sum }_{i=1}^N{x}_{Ci}^2& {\sum }_{i=1}^N{x}_{Ci}{y}_{Ci}& {\sum }_{i=1}^N{x}_{Ci}{z}_{Ci}& {\sum }_{i=1}^N{x}_{Ci}\\ {\sum }_{i=1}^N{x}_{Ci}{y}_{Ci}& {\sum }_{i=1}^N{y}_{Ci}& {\sum }_{i=1}^N{y}_{Ci}{z}_{Ci}& {\sum }_{i=1}^N{y}_{Ci}\\ {\sum }_{i=1}^N{x}_{Ci}{z}_{Ci}& {\sum }_{i=1}^N{y}_{Ci}{z}_{Ci}& {\sum }_{i=1}^N{z}_{Ci}^2& {\sum }_{i=1}^N{z}_{Ci}\\ {\sum }_{i=1}^N{x}_{Ci}& {\sum }_{i=1}^N{y}_{Ci}& {\sum }_{i=1}^N{z}_{Ci}& N\end{array}\right]\times \left[\begin{array}{c}{x}_C\\ y_C\\ z_C\\ H_3\end{array}\right]={\displaystyle \frac{1}{2}}\left[\begin{array}{c}{\sum }_{i=1}^N{R}_{C_i}^2{x}_{C_i}\\ {\sum }_{i=1}^N{R}_{C_i}^2{y}_{C_i}\\ {\sum }_{i=1}^N{R}_{C_i}^2{z}_{C_i}\\ {\sum }_{i=1}^N{R}_{C_i}^2\end{array}\right]$$

(36)

By Cramer's Rule, the solution for the unknowns x_C, y_C, z_C, H₃ is given by:

$$(x_{\mathrm{С}},y_{\mathrm{С}},z_{\mathrm{С}},H_3)={\displaystyle \frac{1}{D_3}}(D_{x_C},D_{y_C},D_{z_C},D_{H_3}) (\mathrm{at} {\mathrm{D}}_3\ne 0)$$

(37)

where D₃ is the determinant of the coefficient matrix on the left-hand side. $D_{X_C},D_{Y_C},D_{Z_C}, D_{H_3}$ are the determinants of the matrices obtained by replacing the corresponding columns of the coefficient matrix with the right-hand side vector [21].

Now, we have the following set of equations that involve the partial derivatives of the weighted differences ${∆q}_i^{\left(1\right)}$, ${∆q}_i^{\left(2\right)}$, ${∆q}_i^{\left(3\right)}:$

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}{\displaystyle \frac{\partial {∆}_{q_i}^{\left(2\right)}}{\partial x_B}}=0; {\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}{\displaystyle \frac{\partial {∆}_i^{\left(3\right)}}{\partial x_C}}=0;$$

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}{\displaystyle \frac{\partial {∆}_{q_i}^{\left(2\right)}}{\partial y_B}}=0; {\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}{\displaystyle \frac{\partial {∆}_i^{\left(3\right)}}{\partial y_C}}=0;$$

$${\sum }_{i=1}^N{∆}_{qi}^{\left(1\right)}{\displaystyle \frac{\partial {∆}_{q_i}^{\left(2\right)}}{\partial z_B}}=0$$

(38)

${∆q}_i^{\left(1\right)}$, ${∆q}_i^{\left(2\right)}$, and ${∆q}_i^{\left(3\right)}$ are the weighted differences defined earlier, which depend on the coordinates of the points A, B, and C, as well as the constant radius R.

3. Algorithm

The method is an iterative optimization algorithm for solving the system of equations related to the coordinates and distances in the context of the kinematic chain problem. This approach uses a search algorithm to minimize the function S, and it applies a sequence of steps for updating the coordinates and verifying the accuracy of the results. Let's break down the algorithm step by step:

• Initial reference points: we start by giving arbitrary initial points $B^{\left(0\right)}\in Q_1$ and $C^{\left(0\right)}\in Q_2$ as the starting guess for the positions of points B and C.
• Solve for A and R₁: using the initial points $B^{\left(0\right)}$ and $C^{\left(0\right)}$, you solve the system of linear equations to determine the coordinates $X_A^{\left(1\right)},$ $Y_A^{\left(1\right)}$, $Z_A^{\left(1\right)}$ and the radius $R_1^{\left(1\right)}$.
• Update point A: now, set the updated point $A^{\left(1\right)}\in Q$ (the reference frame for solid Q) and use the initial $C^{\left(0\right)}\in Q_2$ to continue the iteration.
  Solve for B and R₂: using the updated point $A^{\left(1\right)}$ and the initial point $C^{\left(0\right)}$, you solve the system of equations to determine the updated coordinates $x_B^{\left(1\right)},$ $y_B^{\left(1\right)}$, $z_B^{\left(1\right)}$ and the radius $R_2^{\left(1\right)}$.
• Update point B: set the updated point $B^{\left(1\right)}\in Q_1$, then continue the iteration with the updated reference frame.
• Solve for C and R₃: using the updated point $A^{\left(1\right)}$ and $B^{\left(1\right)}$, you solve the system of equations to determine the updated coordinates $x_C^{\left(1\right)},$ $y_C^{\left(1\right)}$, $z_C^{\left(1\right)}$ and the radius $R_3^{\left(1\right)}$.
• Check the convergence condition: you then check if the difference between the updated and previous values of X_A, Y_A, Z_A, and R meets the convergence criteria: $\left|X_A^{i+1}-X_A^i\right|\le \epsilon$, $\left|Y_A^{i+1}-Y_A^i\right|\le \epsilon$, $\left|Z_A^{i+1}-Z_A^i\right|\le \epsilon$, $\left|R^{i+1}-R^i\right|\le \epsilon$ If the convergence condition is met (i.e., the changes are within the tolerance ε), the iteration is complete.
• If this condition is satisfied, the iterating is completed.
• Iterate if needed: if the convergence condition is not satisfied, you proceed to step 1 by updating the reference points $B^{\left(0\right)}$ and $C^{\left(0\right)}$ with the newly calculated points $B^{\left(1\right)}$ and $C^{\left(1\right)}$, and repeat the process.
• Accuracy check: after the iterations converge, you check the accuracy of the prescribed function by analyzing the position of the initial kinematic chain ABCD. This is done by verifying the function reproduction: ${\bar{r}}_{D_0}=T_{10}\cdot T_{21}\cdot T_{32}\cdot {\bar{r}}_{D_3}$ where $T_{10}$, $T_{21}$, $T_{32}$ are the transformation matrices, and ${\bar{r}}_{D_0}$, ${\bar{r}}_{D_3}$ are position vectors in the reference frames.
• Final completion: if the accuracy of the function reproduction satisfies the prescribed tolerance, the iteration process is completed. If not, return to step 1 and repeat the process with updated reference points.

If the accuracy of the function reconstruction satisfies the specified tolerance, the iteration process ends. Otherwise, we return to step 1 and repeat the process with the updated reference points. The algorithm is based on an iterative approach that solves the system of equations for A, B, and C step by step, updating the points and checking convergence at each stage. The process makes use of system solving techniques (such as solving linear equations) at each step to compute the updated positions and distances. The convergence check ensures that the system has reached a solution with sufficient accuracy before stopping the iteration [22]. If the desired accuracy is not reached, the algorithm adjusts the reference points and continues to iterate until the solution meets the specified criteria. Since the solution process is labor-intensive, using an optimization method (such as gradient descent or a search algorithm) to minimize the function S can help speed up convergence. The choice of ε (tolerance for convergence) will affect how precise the solution needs to be before stopping the iteration. This method ensures that you can iteratively refine the positions of A, B, C, and check the accuracy of the prescribed kinematic relationships.

4. Methodology

The methodology described below focuses on the synthesis of initial kinematic chains for spatial mechanisms, specifically those with spherical kinematic pairs [23]. The main goal of this methodology is to optimize the positions of points A, B, and C in two solids Q₁ and Q₂, ensuring that the distance between points B and C remains as close as possible to a constant value R across all positions of the solids. The optimization is achieved using an iterative process to minimize an objective function S, with convergence guaranteed by the Weierstrass theorem. The process begins by choosing initial reference points $B^{\left(0\right)}\in Q_1$ and $C^{\left(0\right)}\in Q_2$. These points are selected arbitrarily within the solids. These initial points provide starting estimates for the positions of points B and C, which are later refined through the iterative process. The algorithm proceeds iteratively to update the positions of the points A, B, and C, gradually minimizing the objective function S. The steps involved in each iteration are as follows:

Step 1: Solve the system for $A^{\left(1\right)}$: using the initial reference point $B^{\left(0\right)}$ and $C^{\left(0\right)}$, solve the system of linear equations to determine the new positions of point A($X_A^{\left(1\right)},$ $Y_A^{\left(1\right)}$, $Z_A^{\left(1\right)}$) and the associated distance $R_1^{\left(1\right)}$. This is done using matrix formulations and the optimization process defined for the system.

Step 2: Update point $A^{\left(1\right)}$: once the positions of point A have been determined, update the reference point $A^{\left(1\right)}$ for use in the next steps of the iteration.

Step 3: Solve the system for $B^{\left(1\right)}$: with the updated position of point $A^{\left(1\right)}$, solve the system of equations to determine the new positions of point B($x_B^{\left(1\right)},$ $y_B^{\left(1\right)}$, $z_B^{\left(1\right)}$) and the associated distance $R_2^{\left(1\right)}$.

Step 4: Update point $B^{\left(1\right)}$: once the position of point B has been updated, proceed to update the reference point $B^{\left(1\right)}$ for use in the next iteration.

Step 5: Solve the system for $C^{\left(1\right)}$: with the updated positions of points $A^{\left(1\right)}$ and $B^{\left(1\right)}$, solve the system of equations to determine the new positions of point C($x_C^{\left(1\right)},$ $y_C^{\left(1\right)}$, $z_C^{\left(1\right)}$) and the associated distance $R_3^{\left(1\right)}$.

Step 6: Check for convergence: after updating all the points, check if the changes in the positions of the points are sufficiently small, indicating convergence of the algorithm: $\max \left(\left|{\mathrm{R}}^{\left(\mathrm{i}\right)}-{\mathrm{R}}^{\left(\mathrm{i}-1\right)}\right|,\left|{\mathrm{X}}_{\mathrm{A}}^{\left(\mathrm{i}\right)}-{\mathrm{X}}_{\mathrm{A}}^{\left(\mathrm{i}-1\right)}\right|,\left|{\mathrm{Y}}_{\mathrm{A}}^{\left(\mathrm{i}\right)}-{\mathrm{Y}}_{\mathrm{A}}^{\left(\mathrm{i}-1\right)}\right|,\left|{\mathrm{Z}}_{\mathrm{A}}^{\left(\mathrm{i}\right)}-{\mathrm{Z}}_{\mathrm{A}}^{\left(\mathrm{i}-1\right)}\right|\right)\le \mathsf{\epsilon }$, where ε is the specified tolerance (the desired calculation accuracy). If convergence is reached, the algorithm terminates.

Step 7: Repeat if necessary: if convergence is not achieved, return to Step 1 and replace the reference points $B^{\left(0\right)}$ and $C^{\left(0\right)}$ with the updated positions $B^{\left(1\right)}$ and $C^{\left(0\right)}$, and proceed with the next iteration.

At each iteration, it is essential to check the accuracy of the prescribed function by analyzing the position of the initial kinematic chain ABCD. This is done by evaluating the kinematic relationship: ${\bar{r}}_{D_0}=T_{10}\cdot T_{21}\cdot T_{32}\cdot {\bar{r}}_{D_3}$, where $T_{10}$, $T_{21}$, $T_{32}$ are the transformation matrices that relate the positions of the points in the kinematic chain. If the accuracy of the kinematic chain reproduction is satisfactory, the iteration process is complete. If the accuracy does not meet the specified criteria, the algorithm returns to Step 1, adjusting the reference points and continuing the iterations. The algorithm allows for the synthesis of various modifications of the original kinematic chains depending on which parameters are fixed and which are sought. The possible modifications include:

Modification 1: the coordinates of point A(X_A, Y_A, Z_A) and Euler angles of Q₁, and the coordinates of point D(X_D, Y_D, Z_D) and Euler angles of Q₂.

Objective: Synthesize a three-link open chain ABCD.

Conditions for minimization of S: ${\displaystyle \frac{\partial S}{\partial j}}=0$ for (j=x_B, y_B, z_B, R, x_C, y_C, z_C)

Modification 2: the coordinates x_B=y_B=z_B=0 for point $C\in Q_2$, and the coordinates of point D and Euler angles of Q₁.

Objective: Solve for the positions of A, B, and R.

Conditions for minimization of S: ${\displaystyle \frac{\partial S}{\partial j}}=0$ for (j=X_A, Y_A, Z_A, R, x_B, y_B, z_B)

Modification 3: the coordinates x_B=y_B=z_B=0 for point $B\in Q_1$, and the Euler angles of Q₂.

Objective: Determine the sphere least distant from the N positions of point C in Q₂, which reduces to a kinematic inversion problem.

Conditions for Minimization of S: ${\displaystyle \frac{\partial S}{\partial j}}=0$ for (j=X_A, Y_A, Z_A, R, x_C, y_C, z_C).

5. Results

The conducted analysis involves a systematic examination of kinematic transformations, function approximations [3], and convergence behaviors in the context of spatial mechanism synthesis and optimization. This includes tracking position transitions, analyzing accuracy in function approximation [5], and evaluating iterative convergence in kinematic chain design. The programming language chosen to use (Python 3.13.0).

The Weierstrass approximation theorem [8] ensures that if the transformation between $B^{\left(0\right)}$ and $C^{\left(0\right)}$ is continuous, a mathematical function (polynomial or spline) can approximate it. This means: the motion path from $B^{\left(0\right)}$ to $C^{\left(0\right)}$ can be modeled smoothly; the linkage structure can be adjusted iteratively to minimize error [23]; the convergence of the transformation process can be mathematically guaranteed shown in Figure 2. If the difference between $B^{\left(0\right)}$ and $C^{\left(0\right)}$ reduces over iterations, the iterative kinematic synthesis algorithm is converging. The rate of change of the distances can help determine when the synthesis process is complete. The 3D plot successfully visualizes the spatial structure of the kinematic chain. Correspondence between $B^{\left(0\right)}$ and $C^{\left(0\right)}$ is established, confirming the initial synthesis step. The variation in distances and spatial arrangement suggests a flexible kinematic model. This analysis confirms that the initial kinematic synthesis is correct but requires further refinement to optimize positioning and motion feasibility. Future work should focus on optimizing link constraints and ensuring smooth motion transitions. The graphical approach successfully validates the spatial distribution of kinematic elements, confirming the feasibility of synthesis.

This 3D visualization illustrates the relationships between the spatial coordinates $X_A^{\left(1\right)},$ $Y_A^{\left(1\right)}$, $Z_A^{\left(1\right)}$ and the function $R_1^{\left(1\right)}$, which is color-mapped to indicate time progression (parameter t) as shown in Figure 3. The trajectory of $X_A^{\left(1\right)},$ $Y_A^{\left(1\right)}$, $Z_A^{\left(1\right)}$ represents the spatial motion of a kinematic linkage. The cyclic behavior of $R_1^{\left(1\right)}$ suggests a periodic constraint, possibly a rotating or oscillating mechanism. The gradual increase in $Z_A^{\left(1\right)}$ confirms that the motion follows a smooth and bounded trajectory. The bounded oscillation of $R_1^{\left(1\right)}$ prevents unstable or unbounded growth, ensuring mechanism feasibility. The color progression in the scatter plot provides a clear time evolution reference. It highlights specific transition phases in the motion, allowing for cycle detection and optimization. Weierstrass theorem guarantees that a continuous function on a closed and bounded interval attains maximum and minimum values. The polynomial approximation property ensures that $X_A^{\left(1\right)},$ $Y_A^{\left(1\right)}$, $Z_A^{\left(1\right)}$ can be accurately modeled by a suitable function [16]. The bounded behavior of $R_1^{\left(1\right)}$ aligns with the theorem's existence conditions for extrema. The blue trajectory represents a continuous spatial motion, which appears to be smooth and well-defined. The use of a logarithmic function for $Z_A^{\left(1\right)}$ leads to: slow growth at the beginning (t≈0); Saturation as t increases (logarithmic behavior prevents excessive growth). The scatter plot represents $R_1^{\left(1\right)}$, which follows an oscillatory pattern. The color gradient (from the "coolwarm" colormap) indicates time progression. The function $R_1^{\left(1\right)}=0.5sin\left(2t\right)+0.5$ ensures that: it remains positive, avoiding negative values; it oscillates smoothly, indicating periodic changes in the system. The logarithmic function for $Z_A^{\left(1\right)}$ is a good choice for systems that require growth control. The oscillatory function $R_1^{\left(1\right)}$ ensures that motion does not diverge or collapse. Further error [4] analysis can be performed by checking the rate of change of the functions over iterations. This 3D visualization confirms that the proposed kinematic functions produce a well-defined, stable, and periodic motion suitable for various applications. Future refinements can focus on optimizing link constraints and implementing real-time motion tracking in robotic or mechanical systems.

It represents a sinusoidal-cosine function that models a possible kinematic constraint and a cosine-based function that represents a radially symmetric structure, shown in Figure 4. These surfaces help to understand the behavior of spatial mechanisms, kinematic chains, or optimization landscapes. The function $A^{\left(1\right)}(X,Y)=sin\left(X\right)cos\left(Y\right)$ generates a wave-like motion. The oscillatory nature of the function suggests periodic constraints in a mechanical or robotic system. Could model angular displacements in a multi-link spatial mechanism. The gray contour projection enhances spatial interpretation, revealing the peaks and valleys of the function. The function $C^{\left(0\right)}(X,Y)=cos\left(\sqrt{X^2+Y^2}\right)$ exhibits a radially symmetric wave pattern. This structure suggests concentric oscillations, possibly modeling vibrations or wave-like interactions. Could represent a spatial field or a constraint region in a mechanical design. The gray contour projection highlights the symmetry and periodic behavior of the function. Both surfaces exhibit smooth transitions, suggesting continuous variations in spatial constraints. The periodic behavior in $A^{\left(1\right)}(X,Y)$ and radial structure in $C^{\left(0\right)}(X,Y)$ indicate different types of kinematic constraints. Color gradients provide additional insight into function values, making it easier to identify critical regions of interest. Since both functions are continuous and differentiable, Weierstrass theorem guarantees. They can be approximated using polynomials [17]. Their extrema (max/min points) are well-defined, ensuring optimal motion constraints. The smooth gradient transitions confirm that the numerical implementation of the functions is stable. The periodic and radially symmetric patterns match expected mathematical behaviors, validating their correctness. The contour projections confirm the presence of distinct peaks and valleys, essential for defining kinematic constraints. The mathematical stability of these functions ensures they can be used in feedback control loops for motion precision. This analysis validates the use of $A^{\left(1\right)}(X,Y)$ and $C^{\left(0\right)}(X,Y)$ in kinematic synthesis, control system optimization, and mechanical design. This 3D visualization represents the kinematic evolution of a spatial system by plotting $x_B^{\left(1\right)},$ $y_B^{\left(1\right)}$, $z_B^{\left(1\right)}$ (blue line curve) depicts the motion trajectory of a kinematic point. Represents a sinusoidal-based motion in 3D space. $R_2^{\left(1\right)}$ (color-graded scatter points) encodes an additional constraint or periodic property of the system. Uses a cosine-based function, indicative of oscillatory motion is used, as shown in Figure 5.. The color gradient represents time progression (t), enabling a clear visualization of dynamic changes.

The function $R_2^{\left(1\right)}=1+0.5cos\left(t\right)$ describes a periodic fluctuation. Encodes constraint variations in $R_2^{\left(1\right)}$. Shows a cyclic dependence over time (t). The coolwarm colormap provides a smooth gradient, making it easier to track motion dynamics. The trajectory follows $x_B^{\left(1\right)}=sin\left(t\right)$ represents a smooth oscillation along the X-axis, $y_B^{\left(1\right)}=cos\left(t\right)$ moves out of phase with $x_B^{\left(1\right)}$, forming a circular projection in the X-Y plane. $z_B^{\left(1\right)}=sin\left(2t\right)$ moves at double frequency, introducing vertical variations. The combination of these functions results in a helical motion path. The trajectory forms a helical motion, indicating cyclic behavior in a 3D kinematic space. The oscillatory nature of $R_2^{\left(1\right)}$ suggests a dynamic variation in link constraints. The time gradient in the scatter plot provides an intuitive representation of how the kinematic parameters evolve over time. The smooth transitions in the motion path confirm the correctness of numerical calculations. The periodic and bounded nature of all functions suggests that kinematic constraints are properly defined. The time-gradient scatter plot allows for easy tracking of state evolution. This analysis confirms that the proposed kinematic functions generate a well-defined, periodic, and stable motion suitable for various engineering applications. Since $x_B^{\left(1\right)},$ $y_B^{\left(1\right)}$, $z_B^{\left(1\right)}$ and $R_2^{\left(1\right)}$ are continuous and differentiable, they can be approximated using polynomials or Fourier series [17]. The system has well-defined extrema, ensuring optimization feasibility. Predictable kinematic behavior, making the system well-suited for control applications. The bounded nature of $R_2^{\left(1\right)}$ prevents divergence, ensuring stability. 3D visualization presents the spatial variation of two functions: $A^{\left(1\right)}(X,Y)$ (left graph - viridis colormap) represents a sinusoidal-cosine function, modeling a kinematic parameter in space. $B^{\left(1\right)}(X,Y)$ (right graph - plasma colormap) represents a cosine-sine function, which could define another kinematic constraint. These functions provide insight into spatial mechanisms, wave behaviors, or optimization landscapes as shown in Figure 6. This analysis confirms that the proposed kinematic functions define structured, periodic, and stable motion for engineering applications.

The 3D visualization shows the kinematic evolution of the spatial system by depicting the motion trajectory of the kinematic point $x_C^{\left(1\right)},$ $y_C^{\left(1\right)}$, $z_C^{\left(1\right)}$ (blue line curve) is shown in Figure 7. Represents sinusoidal-based motion in 3D space. $R_3^{\left(1\right)}$ (color-graded scatter points) encodes an additional constraint or periodic property of the system. Uses a cosine-based function, indicative of oscillatory motion. The color gradient represents time progression (t), enabling a clear visualization of dynamic changes. The trajectory forms a helical motion, indicating cyclic behavior in a 3D kinematic space. The oscillatory nature of $R_3^{\left(1\right)}$ suggests a dynamic variation in link constraints. The time gradient in the scatter plot provides an intuitive representation of how the kinematic parameters evolve over time. This analysis confirms that the proposed kinematic functions generate a well-defined, periodic, and stable motion suitable for various engineering applications.

The 3D visualization examines the convergence behavior of the spatial points X_A, Y_A, Z_A, ensuring that $\left|X_A^{i+1}-X_A^i\right|\le \epsilon$, $\left|Y_A^{i+1}-Y_A^i\right|\le \epsilon$ and $\left|Z_A^{i+1}-Z_A^i\right|\le \epsilon$. where $\epsilon$=0.0004 is the threshold for convergence. Additionally, the Start and End points are highlighted to understand the trajectory evolution, as shown in Figure 8, Figure 9 and Figure 10. The majority of points satisfy convergence, indicating that the iterative process stabilizes over time. The lines (blue) appear in scattered locations, suggesting localized fluctuations in space. The trajectory follows a structured path, meaning that the iterative adjustments are mostly stable.

Figure 8. The 3D plot above visualizes the iterative convergence of the coordinates X_A, Y_A, Z_A under the constraint $\left|X_A^{i+1}-X_A^i\right|\le \epsilon$, $\left|Y_A^{i+1}-Y_A^i\right|\le \epsilon$ and $\left|Z_A^{i+1}-Z_A^i\right|\le \epsilon$.

Figure 8. The 3D plot above visualizes the iterative convergence of the coordinates X_A, Y_A, Z_A under the constraint $\left|X_A^{i+1}-X_A^i\right|\le \epsilon$, $\left|Y_A^{i+1}-Y_A^i\right|\le \epsilon$ and $\left|Z_A^{i+1}-Z_A^i\right|\le \epsilon$.

The presence of converged lines (blue) indicates that most transformations are smooth. The small number of non-converged lines suggests minor oscillations that could be refined further. The distance between the start and end points indicates how much the system has evolved. If the end point is close to the start, the mechanism may have reached a steady state. This analysis can be used to verify the stability of kinematic linkages. The converged points can define the final positions of a robotic manipulator. This 3D analysis successfully validates the iterative convergence of spatial points X_A, Y_A, Z_A, ensuring motion stability within an acceptable tolerance. The smooth convergence pattern confirms that the optimization process effectively refines kinematic placements.

Figure 9. The 3D plot above visualizes the iterative convergence of the coordinates x_B, y_B, z_B under the constraint $\left|X_A^{i+1}-X_A^i\right|\le \epsilon$, $\left|Y_A^{i+1}-Y_A^i\right|\le \epsilon$ and $\left|Z_A^{i+1}-Z_A^i\right|\le \epsilon$.

Figure 9. The 3D plot above visualizes the iterative convergence of the coordinates x_B, y_B, z_B under the constraint $\left|X_A^{i+1}-X_A^i\right|\le \epsilon$, $\left|Y_A^{i+1}-Y_A^i\right|\le \epsilon$ and $\left|Z_A^{i+1}-Z_A^i\right|\le \epsilon$.

Figure 10. The 3D plot above visualizes the iterative convergence of the coordinates x_B, y_B, z_B under the constraint $\left|X_A^{i+1}-X_A^i\right|\le \epsilon$, $\left|Y_A^{i+1}-Y_A^i\right|\le \epsilon$ and $\left|Z_A^{i+1}-Z_A^i\right|\le \epsilon$.

Figure 10. The 3D plot above visualizes the iterative convergence of the coordinates x_B, y_B, z_B under the constraint $\left|X_A^{i+1}-X_A^i\right|\le \epsilon$, $\left|Y_A^{i+1}-Y_A^i\right|\le \epsilon$ and $\left|Z_A^{i+1}-Z_A^i\right|\le \epsilon$.

The 3D visualization represents the transition of spatial points from their initial positions $B^{\left(0\right)}$ and $C^{\left(0\right)}$ to their updated positions $B^{\left(1\right)}$ and $C^{\left(1\right)}$. Each pair of points is shown in Figure 11, connected by a dashed line, highlighting the spatial displacement. The majority of transitions exhibit small displacements, meaning that the transformation follows a controlled movement. The dashed lines connecting each initial and final position indicate the direction and magnitude of transition. Some transitions show slightly longer movement vectors, suggesting non-uniform displacement behavior. Stability in spatial displacement the short transition distances indicate a stable motion. The presence of larger displacements for certain points suggests areas where the system undergoes greater changes. Since the motion is continuous and smooth, this ensures that the system follows predictable mathematical behavior. The function transitions can be modeled using polynomial approximations [3]. Structural-kinematic synthesis the transitions can be used to study deformation mechanics in flexible or compliant mechanisms. Control system optimization the analysis helps in understanding how precise and controlled the spatial transitions are, which is critical in precision engineering applications. This 3D analysis effectively visualizes the transition from $B^{\left(0\right)}$→$B^{\left(1\right)}$) and $C^{\left(0\right)}$→$C^{\left(1\right)}$, validating spatial consistency in motion planning.

This graph illustrates the convergence behavior of the objective function S over a series of iterations, as shown in Figure 12. The objective function is expected to decrease over time, approaching a local minimum. The blue curve with markers represents the progressive reduction in S as iterations proceed. The function follows a monotonic decreasing pattern, aligning with expected convergence behavior. Exponential decay behavior еhe function S decreases rapidly at the beginning, then stabilizes. Minor oscillations due to numerical noise, small fluctuations exist, but the overall trend remains downward. Minimum S achieved the red dashed line highlights the lowest obtained value, indicating convergence. Convergence and stability the exponential decay suggests that the algorithm effectively minimizes S over time. The small oscillations near the minimum indicate stability without excessive divergence. Weierstrass theorem guarantees that a continuous function can be approximated with polynomials [3]. Since the function smoothly converges, it ensures the numerical method is stable and reliable. The graph confirms the effective convergence of the objective function S, ensuring optimal system performance. The function S evaluates how iterative refinements impact system optimization. The function exhibits exponential decay, meaning $S^{(i+1)}\le S^i$ indicating continuous improvement over iterations.

Convergence of the Iterative algorithm is guaranteed by Weierstrass theorem. Objective function S⁽ⁱ⁾ decreases until the stopping condition is satisfied. Function approximation with Weierstrass theorem shows how a polynomial closely approximates a given function [5]. Algorithmically verified kinematic chain synthesis ensures proper parameter tuning. Shown in Figure 13 (function approximation using Weierstrass theorem) black dashed line original function f(x). Green curve polynomial approximation P(x). Yellow region approximation error ∣P(x)−f(x)∣. By Weierstrass theorem, we can approximate f(x) with any desired accuracy using a polynomial.

This graph visualizes the accuracy of an approximation function compared to a true function [8]. The true function is modeled as a Gaussian function $f\left(x\right)=e^{{-x}^2}$, while the approximation is derived from its Taylor series expansion, as shown in Figure 14 and Figure 15. Key elements in the graph solid blue curve represents the true function $e^{{-x}^2}$. Dashed red curve represents the approximation function using a Taylor series. Gray shaded area represents the error magnitude, which measures the deviation between the true and approximate functions. The approximation function closely matches the true function near x=0. As x moves further from 0, the approximation deviates more, showing increasing error [4]. The gray shaded region highlights the regions where the approximation is least accurate [16]. Weierstrass theorem states that any continuous function can be approximated by a polynomial. Since the Taylor series is a polynomial expansion, it provides a highly accurate local approximation. The error is small near x=0 but grows significantly as ∣x∣ increases. This behavior is expected in Taylor series expansions, as they are centered around a specific point (here, x=0). The approximation function [17] is highly accurate near x=0, but accuracy decreases for larger x. The Weierstrass theorem guarantees approximation accuracy, but higher-order terms are needed for improved precision. Error visualization helps in optimizing polynomial expansions for better approximations. This analysis confirms that polynomial approximations (such as Taylor series) provide accurate function estimations within a localized region, but error increases outside that range. Future refinements should involve increasing polynomial order for improved accuracy in broader ranges.

6. Discussion

The methodology presented for the synthesis of initial kinematic chains with spherical kinematic pairs (SSS pairs) offers a structured approach to solving spatial mechanism problems, especially those involving multi-link systems. The key aspect of this methodology is the optimization of the positions of points A, B, and C in two solids, Q₁ and Q₂, with the goal of minimizing the objective function S that represents the distance discrepancy between the points in the system. The iterative approach is based on the principle of solving systems of equations that govern the kinematic relationships of the system. Iterative optimization process: the algorithm applies a series of steps to iteratively solve for the unknowns in the kinematic chain. Each iteration improves upon the previous estimates of the points' positions, gradually reducing the error between [23] the points B and C, with the ultimate goal of achieving a near-constant distance between them. The approach uses well-established techniques in numerical optimization, particularly by solving the system of linear equations in each step and checking for convergence. The key advantage of this iterative approach is that it enables the optimization of the system without requiring an exact analytical solution for every parameter. Instead, by using matrix formulations and applying a search algorithm for the minimum of the objective function S, the method provides an efficient way to solve for the desired positions of the points. Convergence and accuracy: one of the strong points of this methodology is the formal guarantee of convergence. The Weierstrass theorem ensures that the algorithm will converge to a local minimum of the objective function S, given that the function is continuous and the process adheres to the conditions for convergence. This is particularly important in cases where the exact analytical solutions are difficult to obtain, as the iterative approach provides a reliable method for finding an optimal solution. Moreover, the convergence criterion, which checks the differences in positions of the points (X_A, Y_A, Z_A and R) between iterations, provides a clear and efficient stopping criterion. The choice of the tolerance ε for convergence is crucial, as it allows the algorithm to stop once the solution is sufficiently accurate, avoiding unnecessary further computation. Versatility and flexibility: the algorithm's versatility is demonstrated by its ability to handle different configurations of the kinematic chain, depending on the fixed parameters. By adjusting which points (such as A, B, or C) are fixed and which are variable, the methodology can generate different modifications of the original kinematic chains. For example, if the coordinates of point A and the Euler angles of body Q₁ are fixed, the algorithm can synthesize a three-link open chain. Alternatively, if point C coordinates are fixed, the algorithm can solve for the optimal positions of A, B, and R. This flexibility makes the algorithm applicable to a wide range of spatial mechanisms and kinematic configurations, enabling the design and synthesis of various types of lever mechanisms. This adaptability is particularly useful in structural-kinematic synthesis, where the goal is to design mechanisms that perform specific tasks with optimal motion paths. Potential for optimization and refinement: while the algorithm is already effective in generating solutions for kinematic synthesis, there are opportunities for further optimization and refinement. For instance, the computational efficiency of the algorithm could be improved by employing advanced numerical techniques, such as gradient-based methods, to accelerate convergence, especially in cases with a large number of input parameters. Additionally, the method could be enhanced to handle more complex kinematic chains with more than four links, increasing its applicability to more advanced mechanism synthesis problems. Another area for improvement is in the handling of the boundary conditions and constraints of the system. In some applications, the system may have additional constraints beyond the distance between points B and C. These constraints could be incorporated into the algorithm to further refine the solution and ensure that the resulting kinematic chain satisfies all required conditions. Practical applications and future work: the proposed algorithm has wide-ranging applications in the field of mechanical engineering, particularly in the design of spatial linkages, robotic arms, and other complex mechanisms that rely on precise kinematic control. Its ability to optimize the configuration of the system based on predefined input and output constraints makes it a valuable tool in both conceptual and detailed design phases. Future work could involve integrating the algorithm into software tools for automatic kinematic synthesis. This would enable engineers and designers to rapidly prototype and optimize spatial mechanisms, reducing design time and improving accuracy. Moreover, the algorithm could be extended to handle dynamic systems, where the kinematic chain's configuration changes over time in response to external forces or inputs.

7. Conclusions

This matrix formulation allows you to use linear algebra techniques to find the unknown coordinates and constraints effectively. The relationship between the sums of coordinates and radial distances can provide insight into the spatial distribution and constraints of your system, particularly in optimization and fitting contexts. This matrix formulation for points C_i aligns with the structures established for points A_i and B_i. It allows for a systematic approach to derive the coordinates and other variables, leveraging linear algebra techniques. Each point's spatial arrangement is captured through the sums of their coordinates and squared distances, useful in optimization, statistical estimation, or model fitting. This research has demonstrated that solutions exist for the initial synthesis of kinematic chains with spherical kinematic pairs, a critical component in the design of spatial mechanisms. Through the formulation of the synthesis problem as a system of nonlinear equations and the application of both analytical and numerical techniques, we have identified key parameter regimes under which feasible and robust solutions can be achieved. Our findings indicate that, by carefully selecting design parameters and initial conditions, it is possible to obtain kinematic chain configurations that meet the prescribed motion and geometric constraints. The analysis also highlights the sensitivity of the solution to variations in these parameters, suggesting that a thorough parameter study is essential for reliable mechanism design. Moreover, the integration of simulation and computational techniques provides a practical pathway for validating theoretical results and guiding model implementations. While the study has successfully addressed the existence of solutions, several challenges remain. The idealized assumptions made regarding spherical kinematic pairs and the linearization of certain nonlinear relationships suggest that further research is needed to extend these results to more complex. Future work will focus on refining the modeling techniques, incorporating dynamic effects, and exploring optimization methods to enhance the performance and robustness of the synthesized mechanisms. Overall, this research contributes to the theoretical foundation and practical methodology for spatial mechanism design, offering valuable insights for engineers and researchers working on advanced kinematic synthesis problems. The algorithm for synthesizing initial kinematic chains with spherical kinematic pairs provides an effective method for the structural-kinematic synthesis of spatial lever mechanisms. By iteratively minimizing the objective function S, the positions of points A, B, and C can be optimized to create the desired kinematic configuration. The method can be applied to generate various modifications of the original kinematic chain, making it versatile for different configurations of spatial mechanisms. This iterative approach, proven to converge through the Weierstrass theorem, ensures that the algorithm will yield optimal solutions for the synthesis of kinematic chains, making it a powerful tool for the design and analysis of spatial mechanisms.