Iterative Synthesis of Initial Kinematic Chains with Rotational Pairs Based on Stepwise Minimization of Geometric Deviations

1. Introduction

This paper explores the challenges of kinematic synthesis for four-link initial kinematic chains, which are fundamental in the design of planar (flat) lever mechanisms. Such synthesis is especially important in developing high-class mechanisms [1,2,3], where specific motion laws are imposed on both input and output links [4]. The core of the synthesis problem lies in defining a plane-parallel motion of two rigid bodies (or their planar cross-sections) with respect to a third, fixed body [5,6]. The approach is grounded in the theory of approximate kinematic geometry, which allows for practical realization of mechanisms based on discrete positional data. The synthesis of four-link initial kinematic chains with rotational pairs is a fundamental problem in the design of planar lever mechanisms, especially high-class mechanisms with prescribed laws of motion for input and output links. This problem arises frequently in the design of machinery where precise coordination of moving components is required, such as robotic manipulators, mechanical linkages [7,8], and automated systems. The challenge lies in finding appropriate geometric configurations that satisfy desired motion criteria while ensuring kinematic compatibility and smooth operation [9,10]. Traditional methods for mechanism synthesis often rely on graphical or analytical approaches, which can be cumbersome and limited in their application. As mechanisms become more complex, particularly those involving rotational and translational pairs [11], the design process requires more robust mathematical tools. Recent advancements in computational techniques have enabled the development of more efficient and accurate methods for solving synthesis problems [12,13,14], particularly when formulating them as optimization problems. In this work, we focus on the synthesis of four-link initial kinematic chains formed by connecting three moving planes: $Q$ (fixed plane), $Q_1$ (rotating plane), and $Q_2$ (plane-parallel moving plane). The objective is to determine the optimal positions of points $A\left(X_A,Y_A\right)\in Q$, $B\left(x_B,y_B\right)\in Q_1$, and $C\left({x\prime }_C,{y\prime }_C\right)\in Q_2$ such that the distance between B and C remains approximately constant across various configurations. This problem formulation is particularly important when designing open kinematic chains of the form ABCD with three degrees of freedom, which serve as the foundation for more complex mechanisms [1]. The design and synthesis of planar mechanisms with specified motion characteristics is a fundamental problem in mechanical engineering and robotics. In particular, the synthesis of four-link kinematic chains plays a crucial role in developing precise and efficient linkage systems used in a wide range of industrial and technological applications, such as manipulators, tracing devices, and automation equipment. One of the central tasks in kinematic synthesis is to determine the optimal geometric configuration of links and joints that can reproduce a prescribed set of positions and orientations of moving bodies. This paper focuses on the synthesis of initial four-link kinematic chains with rotational pairs [15,16,17], which are formed by connecting two moving planes through hinges to a fixed base, creating an open-loop mechanism with three degrees of freedom [18]. Traditional methods of mechanism synthesis often involve solving highly nonlinear systems of algebraic equations derived from loop-closure conditions or geometric constraints. These methods, while exact in theory, become computationally intensive and unwieldy when applied to systems with multiple positional constraints. Moreover, the multiplicity of solutions (including local maxima, saddle points, and non-physical configurations) requires careful selection and verification of valid designs. To address these challenges, we propose an efficient and theoretically grounded approach based on approximate kinematic geometry [1,20]. The method relies on minimizing a weighted sum-of-squares error function that quantifies the deviation of the actual mechanism from the desired motion behavior. By expressing the geometric relationships in matrix form and breaking down the full system into a sequence of simpler linear subsystems, the synthesis problem becomes tractable and can be solved iteratively. A key result of this work is a theorem that guarantees the monotonic decrease of the objective function across successive iterations, ensuring convergence to a local minimum. This allows for practical implementation of the algorithm without resorting to the full solution of high-degree polynomial systems [21]. Furthermore, the algorithm’s flexibility permits it to be applied to a broad class of mechanisms, as long as no strict parameter constraints are imposed. The remainder of the paper details the mathematical formulation of the problem, derivation of the objective function, linearization strategy, the minimization theorem, and the iterative synthesis algorithm. Visualization of the iterative process is also presented to demonstrate convergence behavior and support practical implementation. The problems of mechanism synthesis have mainly developed in two directions: algebraic methods of approximation synthesis of lever mechanisms and methods of synthesis of lever mechanisms based on approximation kinematic geometry [22,23]. The development of these two directions of mechanism synthesis has had a decisive influence on the formation of methods of approximation synthesis of lever mechanisms [24]. Kinematic synthesis is one of the most important stages in the process of designing a mechanism, since it is at this stage that the basic kinematic properties necessary for the mechanism to perform its assigned functions are formed. The advantages of optimization methods for mechanism synthesis are especially evident in cases where "classical" methods of kinematic synthesis based on kinematic geometry or various approximation methods are inapplicable or ineffective. The success of searching for an optimal mechanism largely depends on the choice of the initial approximation determined by classical methods [11], while a mechanism designed by classical methods often requires optimization taking into account additional synthesis conditions. Based on the above, it can be concluded that the development of methods for synthesizing flat lever mechanisms was aimed at creating a theory of synthesis that would make it possible to synthesize any mechanism schemes based on a single principle. In particular, many studies at the heart of solving problems of structural and kinematic synthesis of mechanisms are essentially the idea of forming a mechanism from certain initial kinematic chains, as from component parts. The purpose of this work is to develop a unified set of numerical and analytical methods for the structural and kinematic synthesis of planar lever mechanisms, based on the prescribed laws of motion of the input and output links. The methods are designed to be applicable across all stages of the synthesis process and are grounded in the systematic use of initial kinematic chains (IKC) and their structural modifications [25]. This approach enables the automated generation of mechanism topologies and dimensional parameters that satisfy both geometric and functional motion constraints, thereby facilitating the efficient and accurate design of high-class planar mechanisms through a computationally supported framework. Reliability of theoretical results is ensured by correct use of theoretical provisions and methods of theoretical and applied mechanics and applied mathematics. The efficiency of the method and algorithms is confirmed by the results of numerical modeling on a computer. Numerical results of the synthesis of specific mechanisms confirm the operability and efficiency of the developed method of structural and kinematic synthesis of flat lever mechanisms according to the specified laws of motion of the input and output links of the mechanism.

2. Materials and Methods

Let there be given N discrete positions of two planes, denoted as Q₁ and Q₂, which represent sections of two rigid bodies. Their positions relative to a fixed reference plane Q are specified as follows: plane Q₁ is defined by coordinates (X_A, Y_A) and a rotation angle ${\phi }_i$ about point A, where $i=1,\dots ,N$; plane Q₂ is defined by coordinates ($X_{Di}, Y_{Di}$) and a rotation angle ${\psi }_i$ about point D, for the same $i=1,\dots ,N$:

$$Q_1\left(X_A,Y_A,{\phi }_i\right), Q_2\left(X_{Di},Y_{Di},{\psi }_i\right), i=1,\dots ,N$$

(1)

To formalized the motions: A rectangular coordinate system OXY is rigidly attached to the fixed plane Q. The moving planes Q₁ and Q₂ each carry their own coordinate systems: $A_{xy}$, fixed in Q₁; $D_i{x}_i^{\prime }{y}_i^{\prime }$, fixed in Q₂ for each position $i$. In this setup: plane Q₁ undergoes pure rotational motion about the fixed-point A; plane Q₂ executes a general plane-parallel motion [1], governed by the trajectory of point D and rotation $\psi$. The aim is to synthesize a kinematic chain-typically a four-bar linkage-capable of guiding these two bodies through the given set of positions with high precision, ideally minimizing deviations using methods of approximate synthesis. Where X_A, Y_A - is the fixed-point A on the base plane Q; $x_B, y_B$ - is the coordinates of hinge B on Q₁; $x_C, y_C$ - is the coordinates of hinge C on Q₂; $X_{Di}, Y_{Di}$ - is the path (trajectory) of point D on Q₂; ${\phi }_i$ - is the rotation angle of Q₁ in position $i$; ${\psi }_i$ - is the rotation angle of Q₂ in position $i$; a, c - is the coefficients of straight line equation (e.g. for translational guidance); R - length of connecting link. Research manuscripts reporting large datasets that are deposited in a publicly available database should specify where the data have been deposited and provide the relevant accession numbers. If the accession numbers have not yet been obtained at the time of submission, please state that they will be provided during review. They must be provided prior to publication. To achieve coordinated motion of the two moving planes Q₁ and Q₂, geometric connections must be established between them. These connections can be implemented in two principal ways, leading to different types of four-link kinematic chains: 1. Kinematic chain with a two-hinged link. In this case, a rigid link with rotational (revolute) pairs is used to directly connect the two moving planes: The resulting mechanism is a closed four-link kinematic chain ABCD, consisting entirely of rotational joints and having three degrees of freedom [11]; The configuration is illustrated schematically in Figure 1c; The task is to determine the positions of three key points: $A\in Q$ - a fixed hinge on the stationary plane; $B\in Q_1$ - a hinge point on the first moving plane; $C\in Q_2$ a hinge point on the second moving plane. These points must be chosen such that the mechanism ensures the specified laws of motion (i.e., prescribed positions and orientations of Q₁ and Q₂) across all N configurations. 2. Kinematic chains with a link group (open chains). Alternatively, coordination can be achieved using a link group that combines rotational and translational pairs: The guide link is attached to one of the moving planes, and the stone link (slider or connecting link) is pivotally attached to the other. This configuration yields two open variants of the four-link chain: Variant ACD: rotational pair at A, translational pair along line CD; Variant ABD: rotational pair at A, translational pair along line BD. Both variants retain three degrees of freedom [2], but they differ in how motion is transferred between Q₁ and Q₂. In this approach, synthesis requires: Determining the coefficients (slope and intercept) of the guiding line equation, along which translational motion occurs. Specifying the position of the hinge: $B\in Q_1$ in the ABD configuration; $C\in Q_2$ in the ACD configuration. Identifying the fixed hinge point $A\in Q$ on the stationary plane. This Figure 2 presents the geometric configurations of initial four-link kinematic chains with rotational pairs formed by connecting two moving planes Q₁ and Q₂ through a fixed base plane $Q$. The subfigures illustrate the different problem formulations and corresponding parameter groupings used in the synthesis algorithm: Full configuration of the four-bar chain ABCD, showing the fixed base point $A\in Q$, and moving points $B\in Q_1$, $C\in Q_2$, and the constant inter-hinge distance R. Intermediate configurations for each minimization step: Solving for A while B and C are fixed; Solving for B while A and C are fixed; Solving for C while A and B are fixed. Reduced configurations corresponding to boundary cases (e.g., fixing one point or setting coordinates to zero) for simplifying the linear systems [3,4]. Each subfigure corresponds to a stage in the stepwise synthesis algorithm, supporting convergence to a locally optimal mechanism by solving linearized subsystems (as described in systems Eq. (17), Eq. (21), and Eq. (24)).

This Figure 2 presents geometric configurations for synthesizing kinematic chains of type ACD and AD, where the connection between the moving planes Q₁ and Q₂ is implemented using one rotational and one translational pair. The aim is to constrain the motion of one point (typically on $Q_2$) to follow an approximately linear path relative to another plane. Initial configurations with known base point $A\in Q$ and a point $C\in Q_2$ moving along a path defined implicitly by the guiding line. Simplified geometric construction of the motion constraint imposed on point $C$. Line equation $ax+by+c=0$ is imposed on point $C$ trajectory relative to plane Q₁. Plane-parallel motion and angular orientation of plane $Q_2$ shown relative to base $Q$.

Reduced cases where specific coordinates are zeroed $x_C=y_C=0$ for simplification in solving the associated linear systems. These configurations correspond to problem formulations using expressions for ${∆}_{qi}^{\left(3\right)}$, and systems Eq. (24) and Eq. (23), leading to identification of optimal line parameters a, c, and base point coordinates. This Figure 2 supports the second synthesis problem discussed in the paper, where the objective is to determine optimal hinge placement and line constraints to guide motion coordination through translational pair synthesis.

This Figure 3 depicts the geometric arrangements for the synthesis of ABD and AD kinematic chains, where the moving planes Q₁ and Q₂ are connected through one fixed base plane $Q$ using a combination of rotational and translational pairs. In this formulation, the trajectory of point $B\in Q_1$ is constrained relative to plane Q₂ via an approximate straight-line path [5,6,7,8], while point $A\in Q$ serves as the rotation base. Initial geometric setup where the straight-line guiding point $B$ is implicitly defined through transformations relative to Q₂. Local coordinate system and motion vector orientation of point $D_i\in Q_2$.

Explicit formulation of the constraint line $ax+by+c=0$, through which point $B$ trajectory is minimized. Graphical construction of synthesis problem to determine the optimal hinge coordinates and line parameters on Q₁ and Q₂. Reduced boundary cases for solving simplified systems by setting $x_B=y_B=0$, enabling direct solution of line parameters a, c for motion planning in the translational domain. These configurations support the third synthesis problem of the study [9,10,11], where the objective is to determine optimal placement of the base point $A\in Q$, line coefficients a, c, and hinge location $B\in Q_1$ such that the movement is constrained according to translational behavior. The Figure 3c corresponds to solving systems such as Eq. (24) and Eq. (23) with appropriate linearization.

3. Synthesis of Initial Kinematic Chain with Rotation Pairs (ABCD Chain)

To determine the linear dimensions of the closed four-link initial kinematic chain ABCD is shown in Figure 1c, we formulate the synthesis task as a geometric optimization problem. The aim is to ensure that the distance between two moving points - one on each of the mobile planes $Q_1$ and $Q_2$ - remains approximately constant across multiple specified positions [5,12].

Problem 1. Synthesis of ABCD initial kinematic chain with rotation pairs.

$A\left(X_A, Y_A\right)\in Q$ – on the stationary base plane;

$B\left(x_B,y_B\right)\in Q_1$ - fixed in the local coordinate system of the first moving plane;

$C\left({x\prime }_C,{y\prime }_C\right)\in Q_2$ - fixed in the local coordinate system of the second moving plane.

Such that the distance between $B$ and $C$ remains close to a constant $R$ for all given finite relative positions $i=\mathrm{1,2},\dots ,N$ of the planes $Q_1$ and $Q_2$.

For each position $i$, we are given:

plane $Q_1$: defined by rotation angle ${\phi }_i$ around point $A$;

plane $Q_2$: defined by rotation angle ${\psi }_i$ and displacement $X_{Di}$, $Y_{Di}$ of point $D$.

The coordinates of points $B$ and $C$ at each position $i$, relative to the global system $OXY$, are obtained using homogeneous transformation matrices.

Point $B$ in $Q_1$ (rotation about $A$):

$$\left[\begin{array}{c}{x}_{Bi}\\ y_{Bi}\\ 1\end{array}\right]=\left[\begin{array}{ccc}cos{\phi }_i& -sin{\phi }_i& X_A\\ sin{\phi }_i& cos{\phi }_i& Y_A\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_B\\ y_B\\ 1\end{array}\right]$$

(2)

Point $C$ in $Q_2$ (undergoing planar motion):

$$\left[\begin{array}{c}{x}_{Ci}\\ y_{Ci}\\ 1\end{array}\right]=\left[\begin{array}{ccc}cos{\psi }_i& -sin{\psi }_i& X_{Di}\\ sin{\psi }_i& cos{\psi }_i& Y_{Di}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x\prime }_C\\ {y\prime }_C\\ 1\end{array}\right]$$

(3)

For each position $i$, the constraint is:

$${\left(x_{Ci}-x_{Bi}\right)}^2+{\left(y_{Ci}-y_{Bi}\right)}^2\approx R^2$$

(4)

The deviation from this constraint across all $N$ positions is minimized using an objective function, typically the sum of squared deviations:

$$S=\sum _{i=1}^N{\left[{\left(x_{Ci}-x_{Bi}\right)}^2+{\left(y_{Ci}-y_{Bi}\right)}^2-R^2\right]}^2$$

Unknown parameters to determine (7 total):

$X_A, Y_A$ – are the coordinates of the base point $A\in Q$;

$x_B,y_B$ – are the coordinates of the hinge point $B\in Q_1$ in its local frame;

${x\prime }_C,{y\prime }_C$ – are the coordinates of the hinge point $C{\in Q}_2$ in its local frame;

$R$ – is the approximate link length between points $B$ and $C$;

Problem 2. Synthesis of ACD initial kinematic chain with rotational and translational pairs.

This problem concerns the synthesis of an open initial kinematic chain [13,14] of the form ACD is shown in Figure 2c, in which one link endpoint is constrained to move along a straight line (a prismatic pair), and the other endpoint is pivotally connected to the second moving body.

A base hinge point $A\left(X_A, Y_A\right)\in Q$;

A straight line $ax+by+c=0$ - fixed in the coordinate system of plane $Q_1$;

A point $C\left({x\prime }_C,{y\prime }_C\right)\in Q_2$ – defined in the local frame of $Q_2$.

Such that the trajectory of point $C$ (from plane $Q_2$) approaches the desired straight line (defined in $Q_1$) as closely as possible across N given motion positions.

This task defines the synthesis of an open initial kinematic chain ACD with mixed kinematic pairs:

A is a rotational pair between the fixed base and plane $Q_1$;

$C$ is a point constrained to slide along a line (i.e., a prismatic pair) defined in the frame of $Q_1$.

This is appropriate for mechanisms such as guide-link couplers or line-follower linkages, where one endpoint traces a line and the opposite pivot is free to rotate.

For each position $i$, point $C\in Q_2$ (in its local coordinates ${x\prime }_C$, ${y\prime }_C$) is transformed into the global frame using a homogeneous transformation:

$$\left[\begin{array}{c}{x}_{Ci}\\ y_{Ci}\\ 1\end{array}\right]=\left[\begin{array}{ccc}cos{\psi }_i& -sin{\psi }_i& X_{Di}\\ sin{\psi }_i& cos{\psi }_i& Y_{Di}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x\prime }_C\\ {y\prime }_C\\ 1\end{array}\right]$$

We define the desired guide line on plane $Q_1$ in its local coordinates by $ax+by+c=0$. To express this line in the global frame $OXY$, we apply the inverse transformation of plane $Q_1$ motion for each position $i$. The condition that transformed point $Ci$ lies approximately on the transformed line yields:

$$a\left[\left(x_{Ci}-X_A\right)cos{\phi }_i+\left(y_{Ci}-Y_A\right)sin{\phi }_i\right]+b\left[-\left(x_{Ci}-X_A\right)sin{\phi }_i+\left(y_{Ci}-Y_A\right)cos{\phi }_i\right]+c\approx 0$$

This expression represents the signed distance from point $Ci$ to the moving guide line (transformed to the global system), and we seek to minimize the sum of squared deviations from this line across all positions:

$$S=\sum _{i=1}^N{\left[a·{\xi }_i+b·{\eta }_i-c\right]}^2$$

$${\xi }_i=\left(x_{Ci}-X_A\right)cos{\phi }_i+\left(y_{Ci}-Y_A\right)sin{\phi }_i; {\eta }_i=-\left(x_{Ci}-X_A\right)sin{\phi }_i+\left(y_{Ci}-Y_A\right)cos{\phi }_i$$

(5)

Unknown parameters (7 total):

$X_A, Y_A$ – is the position of hinge $A$ on base plane $Q$;

$a, b, c$ – are the coefficients of the line equation (defined in $Q_1$ frame);

${x\prime }_C,{y\prime }_C$ – are the local coordinates of point $C$ on plane $Q_2$.

Problem 3. Synthesis of ABD initial kinematic chain with translational pairs.

This problem addresses the design of an open initial kinematic chain ABD is shown in Figure 3c, in which the motion of a point $B\in Q_1$ is constrained to approximate a straight-line trajectory defined in the local coordinate system of the second moving plane $Q_2$.

A fixed hinge point $A\left(X_A, Y_A\right)\in Q$;

A moving point $B\left(x_B,y_B\right)\in Q_1$ - defined in the local coordinate frame of $Q_1$;

A line $ax+by+c=0$ defined in the coordinates system of plane $Q_2$.

$$ax+by+c=0$$

(6)

Such that the trajectory of point $B$, when transformed through $Q_1$ motion relative to the global base and observed from $Q_2$, remains as close as possible to this specified straight line during the entire motion cycle.

This is the synthesis of an ABD-type open initial kinematic chain with mixed kinematic pairs:

The point $A\in Q$ acts as a revolute joint to the input plane $Q_1$;

The point $B\in Q_1$ is constrained to follow a line on plane $Q_2$, corresponding to a prismatic pair.

Such configurations are frequently found in planar guiding mechanisms [15], single-stroke cams, and reciprocating systems where a point from the input side slides along a path defined by the output body.

The position of point $B\in Q_1$, given in local coordinates $x_B$,$y_B$, transforms to the global $OXY$ frame as:

$$\left[\begin{array}{c}{x}_{Bi}\\ y_{Bi}\\ 1\end{array}\right]=\left[\begin{array}{ccc}cos{\phi }_i& -sin{\phi }_i& X_A\\ sin{\phi }_i& cos{\phi }_i& Y_A\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_B\\ y_B\\ 1\end{array}\right]$$

The straight line in $Q_2$, defined by $ax+by+c=0$ in local coordinates, transforms into the global frame for each position $i$ of $Q_2$, using:

$$a\left[\left(x_{Bi}-X_{Di}\right)cos{\psi }_i+\left(y_{Bi}-Y_{Di}\right)sin\psi \right]+b\left[-\left(x_{Bi}-X_{Di}\right)sin{\psi }_i+\left(y_{Bi}-Y_{Di}\right)cos{\psi }_i\right]+c\approx 0$$

This expression computes the signed distance from point $Bi$ global position to the line defined on $Q_2$ at position $i$.

To ensure closeness of point $B$ to the desired trajectory across all N positions, we minimize the squared sum of line deviations:

$$S=\sum _{i=1}^N{\left[a·{\xi }_i+b·{\eta }_i-c\right]}^2$$

$${\xi }_i=\left(x_{Bi}-X_{Di}\right)cos{\psi }_i+\left(y_{Bi}-Y_{Di}\right)sin\psi ; {\eta }_i=-\left(x_{Bi}-X_{Di}\right)sin{\psi }_i+\left(y_{Bi}-Y_{Di}\right)cos{\psi }_i$$

(7)

Unknowns (7 total):

$X_A, Y_A$ – are the location of the fixed base pivot $A\in Q$;

$x_B,y_B$ – are the position of the moving point $B\in Q_1$;

$a, b, c$ – are the coefficients defining the straight-line constraint in $Q_2$.

This synthesis problem ensures that the motion of a point $B\in Q_1$, when projected onto plane $Q_2$, approximates a known line [6,16]. Together with Problem 2, it completes the strategy for designing initial kinematic chains that involve one translational constraint. Both Problem 2 (ACD) and Problem 3 (ABD) fall under the broader class of synthesis problems with translational pairs, where the guiding geometry is defined in a moving frame and the mechanism must conform to hybrid (rotational and prismatic) link constraints.

4. Weighted Difference for Synthesis with Rotational Pairs

4.1. First Form of the Weighted Difference (Optimizing over $X_A,Y_A,R$)

We want the length between the moving points $B_i$ and $C_i$ to be approximately constant:

$${{∆}_{qi}=\left|B_i{C}_i\right|}^2-R^2{=\left(x_{Ci}-x_{Bi}\right)}^2+{\left(y_{Ci}-y_{Bi}\right)}^2-R^2$$

(8)

This expression is a function of the seven unknown parameters: $X_A, Y_A$, $x_B,y_B$, ${x\prime }_C,{y\prime }_C$, $R$. To calculate ${∆}_{qi}$, we first need to compute $x_{Bi},y_{Bi}$ and $x_{Ci},y_{Ci}$ the positions of points $B$ and $C$ in the global coordinate system.

We define:

$${∆}_{qi}^{\left(1\right)}{=\left(X_{Ai}-X_A\right)}^2+{\left(Y_{Ai}-Y_A\right)}^2-R^2\mathrm{for} i=1,\dots ,N$$

(9)

where $X_{Ai}$, $Y_{Ai}$ is the composite transformed difference vector from $B_i$ to $C_i$, expressed as:

$$\left[\begin{array}{c}{X}_{Ai}\\ Y_{Ai}\\ 1\end{array}\right]=-\left[\begin{array}{ccc}cos{\phi }_i& -sin{\phi }_i& 0\\ sin{\phi }_i& cos{\phi }_i& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_B\\ y_B\\ 1\end{array}\right]+\left[\begin{array}{ccc}cos{\psi }_i& -sin{\psi }_i& X_{Di}\\ sin{\psi }_i& cos{\psi }_i& Y_{Di}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x\prime }_C\\ {y\prime }_C\\ 1\end{array}\right]$$

(10)

This gives the relative vector $\overrightarrow{C_i}$ - $\overrightarrow{B_i}$ mapped into global coordinates, allowing you to separate the search for optimal $\left(X_A, Y_A,R\right)$ while treating the rest as known or intermediate.

We breaking the total error ${∆}_{qi}$ into different function forms (likely for simplification in optimization):

Form 1: ${∆}_{qi}^{\left(1\right)}=f\left(X_A, Y_A,R\right)$

Form 2: ${∆}_{qi}^{\left(2\right)}=f\left(x_B,y_B,R\right)$

Form 3: ${∆}_{qi}^{\left(3\right)}=f\left({x\prime }_C,{y\prime }_C,R\right)$

Each form treats a different group of parameters as variables while holding the rest fixed or transformed into intermediate quantities. This can help in iterative numerical solutions, where we minimize over subsets of parameters alternately.

You define the sum of squared weighted deviations for $N$ given positions:

$$S^{\left(1\right)}=\sum _{i=1}^N{\left({∆}_{qi}^{\left(1\right)}\right)}^2=\sum _{i=1}^N{\left[{\left(X_{Ai}-X_A\right)}^2+{\left(Y_{Ai}-Y_A\right)}^2-R^2\right]}^2$$

(11)

where $X_{Ai}$, $Y_{Ai}$ are functions of known data and the transformation matrices; $X_A$, $Y_A$, $R$ are unknowns to be found.

To minimize $S$, take partial derivatives with respect to each unknown and set them to zero:

$${\displaystyle \frac{\partial S}{\partial X_A}}=0; {\displaystyle \frac{\partial S}{\partial Y_A}}=0; {\displaystyle \frac{\partial S}{\partial R}}=0$$

(12)

From your derivation, using the chain rule and applying the structure of ${∆}_{qi}^{\left(1\right)}$, we arrive at:

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

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

(13)

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

To avoid the meaningless solution $R=0$, assume $R\ne 0$. Then:

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

(14)

This simplifies the previous system. Substituting this result back in, the system becomes:

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

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

(15)

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

These are your final minimization conditions.

Step-by-step substitution from Eq. (9) into system Eq. (15):

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

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

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

(16)

We perform a powerful transformation that simplifies the nonlinear system into a linear system of equations by introducing a new variable:

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

This substitution transforms the previously nonlinear problem into a linear least-squares formulation, enabling direct matrix solution [17]. Let’s now clearly walk through we transformation and express everything in compact, structured matrix form.

Let’s write the system you obtained in matrix form:

$$\left[\begin{array}{ccc}\sum X_{Ai}^2& \sum X_{Ai}{Y}_{Ai}& \sum X_{Ai}\\ \sum X_{Ai}{Y}_{Ai}& \sum Y_{Ai}^2& \sum Y_{Ai}\\ \sum X_{Ai}& \sum Y_{Ai}& N\end{array}\right]\left[\begin{array}{c}{X}_A\\ Y_A\\ H_1\end{array}\right]={\displaystyle \frac{1}{2}}\left[\begin{array}{c}\sum R_{Ai}^2{X}_{Ai}\\ \sum R_{Ai}^2{Y}_{Ai}\\ \sum R_{Ai}^2\end{array}\right]$$

(17)

where $R_{Ai}^2=X_{Ai}^2+Y_{Ai}^2$ This system is linear and can be solved using any standard method (Gaussian elimination, LU decomposition, or matrix inversion), provided $R\ne 0$, which guarantees that the Jacobian of the transformation is non-zero:

$$‖{\displaystyle \frac{\partial \left(X_A,Y_A,H_1\right)}{\partial \left(X_A,Y_A,R\right)}}‖=R$$

As we stated for the first form:

$$\left(X_A,Y_A,H_1\right)={\displaystyle \frac{1}{D_1}}\left(D_{XA},D_{YA},D_{H1}\right), \mathrm{for} D_1\ne 0$$

(18)

Then, from:

$$R=\sqrt{X_A^2+Y_A^2+2H_1}$$

This remains valid and analytically elegant, assuming non-degeneracy of the determinant.

4.2. Second Form of the Weighted Difference (Optimizing over $x_B,y_B,R$)

In this second formulation, the goal is to express the squared deviation ${∆}_{qi}^{\left(2\right)}$ in terms of: unknowns $x_B,y_B,R$; given or pre-specified $X_A$, $Y_A$, ${x\prime }_C,{y\prime }_C$:

$${∆}_{qi}^{\left(2\right)}{=\left(x_{Bi}-x_B\right)}^2+{\left(y_{Bi}-y_B\right)}^2-R^2\mathrm{for} i=1,\dots ,N$$

(19)

where $x_{Bi}$, $y_{Bi}$ are known and computed as in your previous transformation.

This form interprets the deviation of the trajectory of $B_i\in Q_1$ from a constant-radius circle centered at $\left(x_B,y_B\right)$, while $B_i$ is obtained from the known transformation of points $C\prime \in Q_2$ and the offset between moving frames [18].

Transformation of coordinates:

$$\left[\begin{array}{c}{x}_{Bi}\\ y_{Bi}\\ 1\end{array}\right]=\underset{Rotation by {\phi }_i}{\underbrace{\left[\begin{array}{ccc}cos{\phi }_i& -sin{\phi }_i& 0\\ sin{\phi }_i& cos{\phi }_i& 0\\ 0& 0& 1\end{array}\right]}}\left[\begin{array}{c}{X_{Di}-X}_A\\ Y_{Di}-Y_A\\ 1\end{array}\right]+\underset{Relative orientation Q_2\to Q_1}{\underbrace{\left[\begin{array}{ccc}cos\left({\psi }_i-{\phi }_i\right)& -sin\left({\psi }_i-{\phi }_i\right)& 0\\ sin\left({\psi }_i-{\phi }_i\right)& cos\left({\psi }_i-{\phi }_i\right)& 0\\ 0& 0& 1\end{array}\right]}}\left[\begin{array}{c}{x\prime }_C\\ {y\prime }_C\\ 1\end{array}\right]$$

In the current second form, the goal is to minimize the sum:

$$S^{\left(2\right)}=\sum _{i=1}^N{\left({∆}_{qi}^{\left(2\right)}\right)}^2=\sum _{i=1}^N{\left[{\left(x_{Bi}-x_B\right)}^2+{\left(y_{Bi}-y_B\right)}^2-R^2\right]}^2$$

Let:

$$H_2={\displaystyle \frac{1}{2}}\left(R^2-x_B^2-y_B^2\right)⟹R=\sqrt{x_B^2+y_B^2+2H_2}$$

Then, the necessary conditions:

$${\displaystyle \frac{\partial S}{\partial x_B}}=0; {\displaystyle \frac{\partial S}{\partial y_B}}=0; {\displaystyle \frac{\partial S}{\partial R}}=0$$

(20)

lead to the linear system:

$$\left[\begin{array}{ccc}\sum x_{Bi}^2& \sum x_{Bi}{y}_{Bi}& \sum x_{Bi}\\ \sum x_{Bi}{y}_{Bi}& \sum y_{Bi}^2& \sum y_{Bi}\\ \sum x_{Bi}& \sum y_{Bi}& N\end{array}\right]\left[\begin{array}{c}{x}_B\\ y_B\\ H_2\end{array}\right]={\displaystyle \frac{1}{2}}\left[\begin{array}{c}\sum R_{Bi}^2{x}_{Bi}\\ \sum R_{Bi}^2{y}_{Bi}\\ \sum R_{Bi}^2\end{array}\right]$$

(21)

where $R_{Bi}^2=x_{Bi}^2+y_{Bi}^2$.

4.3. Third Form of the Weighted Difference (Optimizing over ${x\prime }_C,{y\prime }_C,R$)

In this third formulation, the goal is to express the squared deviation ${∆}_{qi}^{\left(3\right)}$ in terms of: unknowns ${x\prime }_C,{y\prime }_C,R$; given or pre-specified $X_A$, $Y_A$, $x_B,y_B$:

$${∆}_{qi}^{\left(3\right)}{=\left(x_{Ci}-{x\prime }_C\right)}^2+{\left(y_{Ci}-{y\prime }_C\right)}^2-R^2\mathrm{for} i=1,\dots ,N$$

(22)

where $x_{Сi}$, $y_{Сi}$ are known now derived from transformation involving $Q_2$.

Transformation for $x_{Сi}$, $y_{Сi}$:

$$\left[\begin{array}{c}{x}_{Ci}\\ y_{Ci}\\ 1\end{array}\right]=\underset{Rotation by {\psi }_i}{-\underbrace{\left[\begin{array}{ccc}cos{\psi }_i& sin{\psi }_i& 0\\ -sin{\psi }_i& cos{\psi }_i& 0\\ 0& 0& 1\end{array}\right]}}\left[\begin{array}{c}{X}_A-X_{Di}\\ Y_A-Y_{Di}\\ 1\end{array}\right]+\underset{Relative orientation Q_2\to Q_1}{\underbrace{\left[\begin{array}{ccc}cos\left({\phi }_i-{\psi }_i\right)& -sin\left({\phi }_i-{\psi }_i\right)& 0\\ sin\left({\phi }_i-{\psi }_i\right)& cos\left({\phi }_i-{\psi }_i\right)& 0\\ 0& 0& 1\end{array}\right]}}\left[\begin{array}{c}{x}_B\\ y_B\\ 1\end{array}\right]$$

This equation computes the position of point $C$ in the global coordinate system using relative motion from $Q₂$ point of view.

In the current second form, the goal is to minimize the sum:

$$S^{\left(3\right)}=\sum _{i=1}^N{\left({∆}_{qi}^{\left(3\right)}\right)}^2=\sum _{i=1}^N{\left[{\left(x_{Ci}-{x\prime }_C\right)}^2+{\left(y_{Ci}-{y\prime }_C\right)}^2-R^2\right]}^2$$

Let:

$$H_3={\displaystyle \frac{1}{2}}\left(R^2-{x\prime }_C^2-{y\prime }_C^2\right)⟹R=\sqrt{{x\prime }_C^2+{y\prime }_C^2+2H_3}$$

Then the minimization conditions:

$${\displaystyle \frac{\partial S}{\partial {x\prime }_C}}=0; {\displaystyle \frac{\partial S}{\partial {y\prime }_C}}=0; {\displaystyle \frac{\partial S}{\partial R}}=0$$

(23)

lead to the linear system:

$$\left[\begin{array}{ccc}\sum x_{Ci}^2& \sum x_{Ci}{y}_{Ci}& \sum x_{Ci}\\ \sum x_{Ci}{y}_{Ci}& \sum y_{Ci}^2& \sum y_{Ci}\\ \sum x_{Ci}& \sum y_{Ci}& N\end{array}\right]\left[\begin{array}{c}{x}_C\\ y_C\\ H_3\end{array}\right]={\displaystyle \frac{1}{2}}\left[\begin{array}{c}\sum R_{Ci}^2{x}_{Ci}\\ \sum R_{Ci}^2{y}_{Ci}\\ \sum R_{Ci}^2\end{array}\right]$$

(24)

where $R_{Ci}^2=x_{Ci}^2+y_{Ci}^2$.

Each system is linear and solvable using matrix methods or Cramer’s rule (as we noted, if determinant $D\ne 0$).

4.4. Reduction of Full 7 Variable System to 4 Variables

We are minimizing:

$$S=\sum _{i=1}^N{\left({∆}_{qi}\right)}^2=\sum _{i=1}^N{\left[{\left|B_i{C}_i\right|}^2-R^2\right]}^2$$

Earlier, you broke this down into:

First form ${∆}_{qi}^{\left(1\right)}$ expressed via $X_A,Y_A,R$;

Second form ${∆}_{qi}^{\left(2\right)}$ via $x_B,y_B,R$;

Third form ${∆}_{qi}^{\left(3\right)}$ via ${x\prime }_C,{y\prime }_C,R$.

We now derive the necessary conditions for minimum with respect to:

$${\displaystyle \frac{\partial S}{\partial j}}=0, j=X_A,Y_A,R,x_B,y_B,{x\prime }_C,{y\prime }_C$$

(25)

With:

$${\displaystyle \frac{\partial {∆}_{qi}^{\left(2\right)}}{\partial x_B}}=-2{\tilde{X}}_{Bi}, {\displaystyle \frac{\partial {∆}_{qi}^{\left(2\right)}}{\partial y_B}}=-2{\tilde{Y}}_{Bi}, {\displaystyle \frac{\partial {∆}_{qi}^{\left(3\right)}}{\partial x_C}}=-2{\tilde{X}}_{Ci}, {\displaystyle \frac{\partial {∆}_{qi}^{\left(3\right)}}{\partial y_C}}=-2{\tilde{Y}}_{Ci}$$

And:

$${∆}_{qi}^{\left(1\right)}{=\left(X_{Ai}-X_A\right)}^2+{\left(Y_{Ai}-Y_A\right)}^2-R^2=X_{Ai}{X}_A+Y_{Ai}{Y}_A+H_1-{\displaystyle \frac{1}{2}}{R}_{Ai}^2$$

Then:

$$\sum _{i=1}^N{∆}_{qi}^{\left(1\right)}*{\displaystyle \frac{\partial {∆}_{qi}^{\left(2\right)}}{\partial x_B}}=0⟹\sum _{i=1}^N\left(X_{Ai}{X}_A+Y_{Ai}{Y}_A+H_1-{\displaystyle \frac{1}{2}}{R}_{Ai}^2\right)\left(-2{\tilde{X}}_{Bi}\right)=0$$

Or rearranged:

$$\sum _{i=1}^N\left(X_{Ai}{X}_A+Y_{Ai}{Y}_A+H_1-{\displaystyle \frac{1}{2}}{R}_{Ai}^2\right)X_{Bi}=0$$

$$\sum _{i=1}^N\left(X_{Ai}{X}_A+Y_{Ai}{Y}_A+H_1-{\displaystyle \frac{1}{2}}{R}_{Ai}^2\right)Y_{Bi}=0$$

$$\sum _{i=1}^N\left(X_{Ai}{X}_A+Y_{Ai}{Y}_A+H_1-{\displaystyle \frac{1}{2}}{R}_{Ai}^2\right)X_{Ci}=0$$

$$\sum _{i=1}^N\left(X_{Ai}{X}_A+Y_{Ai}{Y}_A+H_1-{\displaystyle \frac{1}{2}}{R}_{Ai}^2\right)Y_{Ci}=0$$

(26)

Using Cramer’s rule:

$${\mathit{X}}_{\mathit{A}}={\displaystyle \frac{{\mathit{D}}_{\mathit{X}\mathit{A}}}{{\mathit{D}}_1}}, {\mathit{Y}}_{\mathit{A}}={\displaystyle \frac{{\mathit{D}}_{\mathit{Y}\mathit{A}}}{{\mathit{D}}_1}}, {\mathit{H}}_1={\displaystyle \frac{{\mathit{D}}_{\mathit{H}1}}{{\mathit{D}}_1}},$$

we substitute into the four equations and obtain:

$$\sum _{\mathit{i}=1}^{\mathit{N}}\left({\displaystyle \frac{{\mathit{X}}_{\mathit{A}\mathit{i}}{\mathit{D}}_{\mathit{X}\mathit{A}}+{\mathit{Y}}_{\mathit{A}\mathit{i}}{\mathit{D}}_{\mathit{Y}\mathit{A}}+{\mathit{D}}_{\mathit{H}1}}{{\mathit{D}}_1}}-{\displaystyle \frac{1}{2}}{\mathit{R}}_{\mathit{A}\mathit{i}}^2\right){\mathit{X}}_{\mathit{B}\mathit{i}}=0$$

$$\sum _{\mathit{i}=1}^{\mathit{N}}\left({\displaystyle \frac{{\mathit{X}}_{\mathit{A}\mathit{i}}{\mathit{D}}_{\mathit{X}\mathit{A}}+{\mathit{Y}}_{\mathit{A}\mathit{i}}{\mathit{D}}_{\mathit{Y}\mathit{A}}+{\mathit{D}}_{\mathit{H}1}}{{\mathit{D}}_1}}-{\displaystyle \frac{1}{2}}{\mathit{R}}_{\mathit{A}\mathit{i}}^2\right){\mathit{Y}}_{\mathit{B}\mathit{i}}=0$$

$$\sum _{\mathit{i}=1}^{\mathit{N}}\left({\displaystyle \frac{{\mathit{X}}_{\mathit{A}\mathit{i}}{\mathit{D}}_{\mathit{X}\mathit{A}}+{\mathit{Y}}_{\mathit{A}\mathit{i}}{\mathit{D}}_{\mathit{Y}\mathit{A}}+{\mathit{D}}_{\mathit{H}1}}{{\mathit{D}}_1}}-{\displaystyle \frac{1}{2}}{\mathit{R}}_{\mathit{A}\mathit{i}}^2\right){\mathit{X}}_{\mathit{C}\mathit{i}}=0$$

$$\sum _{\mathit{i}=1}^{\mathit{N}}\left({\displaystyle \frac{{\mathit{X}}_{\mathit{A}\mathit{i}}{\mathit{D}}_{\mathit{X}\mathit{A}}+{\mathit{Y}}_{\mathit{A}\mathit{i}}{\mathit{D}}_{\mathit{Y}\mathit{A}}+{\mathit{D}}_{\mathit{H}1}}{{\mathit{D}}_1}}-{\displaystyle \frac{1}{2}}{\mathit{R}}_{\mathit{A}\mathit{i}}^2\right){\mathit{Y}}_{\mathit{C}\mathit{i}}=0$$

(27)

This is the final reduced system in four unknowns: ${\mathit{x}}_{\mathit{B}},{\mathit{y}}_{\mathit{B}}$, ${\mathit{x}\prime }_{\mathit{C}},{\mathit{y}\prime }_{\mathit{C}}$, where ${\mathit{X}}_{\mathit{A}\mathit{i}}$, ${\mathit{Y}}_{\mathit{A}\mathit{i}}$, ${\mathit{R}}_{\mathit{A}\mathit{i}}^2$, ${\mathit{X}}_{\mathit{B}\mathit{i}}$, ${\mathit{Y}}_{\mathit{B}\mathit{i}}$, ${\mathit{X}}_{\mathit{C}\mathit{i}}$, ${\mathit{Y}}_{\mathit{C}\mathit{i}}$ are known, and ${\mathit{D}}_{\mathit{X}\mathit{A}}$, ${\mathit{D}}_{\mathit{Y}\mathit{A}}$, ${\mathit{D}}_{\mathit{H}1}$ are determinantal expressions from system Eq. (17).

We correctly analyzed the structure of system Eq. (28), derived from the necessary conditions:

$$\sum _{\mathit{i}=1}^{\mathit{N}}{∆}_{\mathit{q}\mathit{i}}^{\left(1\right)}{\displaystyle \frac{\partial {∆}_{\mathit{q}\mathit{i}}^{\left(2\right)}}{\partial {\mathit{x}}_{\mathit{B}}}}=0, \sum _{\mathit{i}=1}^{\mathit{N}}{∆}_{\mathit{q}\mathit{i}}^{\left(1\right)}{\displaystyle \frac{\partial {∆}_{\mathit{q}\mathit{i}}^{\left(2\right)}}{\partial {\mathit{y}}_{\mathit{B}}}}=0, \sum _{\mathit{i}=1}^{\mathit{N}}{∆}_{\mathit{q}\mathit{i}}^{\left(1\right)}{\displaystyle \frac{\partial {∆}_{\mathit{q}\mathit{i}}^{\left(3\right)}}{\partial {\mathit{x}\prime }_{\mathit{C}}}}=0, \sum _{\mathit{i}=1}^{\mathit{N}}{∆}_{\mathit{q}\mathit{i}}^{\left(1\right)}{\displaystyle \frac{\partial {∆}_{\mathit{q}\mathit{i}}^{\left(3\right)}}{\partial {\mathit{y}\prime }_{\mathit{C}}}}=0$$

(28)

This results in a system of four nonlinear algebraic equations of degree 9 in the unknowns: ${\mathit{x}}_{\mathit{B}},{\mathit{y}}_{\mathit{B}}$, ${\mathit{x}\prime }_{\mathit{C}},{\mathit{y}\prime }_{\mathit{C}}$.

Thus, the total number of complex roots is: $9^4=6561$. And the maximum number of real roots is approximately: $5^4=625$.

Each of these 625 real roots corresponds to: A local minimum, maximum, saddle, or ridge point of the surface: $\mathit{S}=\mathit{S}\left({\mathit{x}}_{\mathit{B}},{\mathit{y}}_{\mathit{B}},{\mathit{x}\prime }_{\mathit{C}},{\mathit{y}\prime }_{\mathit{C}}\right)$. To distinguish between them, you’d need to analyze the Hessian matrix of second partial derivatives (a 4×4 matrix in this case). The direct algebraic approach to solving this nonlinear system is computationally intractable and inefficient for practical synthesis [19], despite being mathematically rigorous.

A numerical algorithm for minimizing the function $\mathit{S}$, using the fact that it’s a smooth, differentiable function of four variables. Thus, while system Eq. (28) provides an exact algebraic framework, solving it directly leads to a system of four ninth-degree algebraic equations with up to 6561 theoretical solutions, among which only a subset represents valid physical configurations [20]. Verifying the nature of these critical points through second derivative tests would require the evaluation of the Hessian matrix - a process both analytically and computationally intensive. For this reason, the direct algebraic method is impractical. Instead, a more efficient and robust approach involves minimizing the function $\mathit{S}\left({\mathit{x}}_{\mathit{B}},{\mathit{y}}_{\mathit{B}},{\mathit{x}\prime }_{\mathit{C}},{\mathit{y}\prime }_{\mathit{C}}\right)$ numerically, using iterative optimization techniques, which will be presented below based on the following theorem.

Theorem 1. Let the minimal values of the sum of squared differences $\mathit{S}$ for $\mathit{N}$ given positions be defined as follows:

$${\mathit{S}}_1^{\left(1\right)}=\mathit{m}\mathit{i}\mathit{n}\mathit{S}\left[{\mathit{X}}_{\mathit{A}},{\mathit{Y}}_{\mathit{A}},\mathit{R},{\mathit{x}}_{\mathit{B}}^{\left(0\right)},{\mathit{y}}_{\mathit{B}}^{\left(0\right)},{\mathit{x}\prime }_{\mathit{C}}^{\left(0\right)},{\mathit{y}\prime }_{\mathit{C}}^{\left(0\right)}\right];$$

$${\mathit{S}}_2^{\left(1\right)}=\mathit{m}\mathit{i}\mathit{n}\mathit{S}\left[{\mathit{X}}_{\mathit{A}}^{\left(1\right)},{\mathit{Y}}_{\mathit{A}}^{\left(1\right)},\mathit{R},{\mathit{x}}_{\mathit{B}},{\mathit{y}}_{\mathit{B}},{\mathit{x}\prime }_{\mathit{C}}^{\left(0\right)},{\mathit{y}\prime }_{\mathit{C}}^{\left(0\right)}\right];$$

$${\mathit{S}}_3^{\left(1\right)}=\mathit{m}\mathit{i}\mathit{n}\mathit{S}\left[{\mathit{X}}_{\mathit{A}}^{\left(1\right)},{\mathit{Y}}_{\mathit{A}}^{\left(1\right)},\mathit{R},{\mathit{x}}_{\mathit{B}}^{\left(1\right)},{\mathit{y}}_{\mathit{B}}^{\left(1\right)},{\mathit{x}\prime }_{\mathit{C}},{\mathit{y}\prime }_{\mathit{C}}\right].$$

Let:

$\left({\mathit{x}}_{\mathit{B}}^{\left(0\right)},{\mathit{y}}_{\mathit{B}}^{\left(0\right)}\right)\in {\mathit{Q}}_1$ and $\left({\mathit{x}\prime }_{\mathit{C}}^{\left(0\right)},{\mathit{y}\prime }_{\mathit{C}}^{\left(0\right)}\right)\in {\mathit{Q}}_2$ be the initial guesses, $\left({\mathit{X}}_{\mathit{A}}^{\left(1\right)},{\mathit{Y}}_{\mathit{A}}^{\left(1\right)},{\mathit{R}}_1^{\left(1\right)}\right)$ be the solution minimizing $\mathit{S}$ over ${\mathit{X}}_{\mathit{A}},{\mathit{Y}}_{\mathit{A}},\mathit{R}$ (from system Eq. (17)) for fixed ${\mathit{B}}^{\left(0\right)}$, ${\mathit{C}}^{\left(0\right)}$;

$\left({\mathit{x}}_{\mathit{B}}^{\left(1\right)},{\mathit{y}}_{\mathit{B}}^{\left(1\right)},{\mathit{R}}_2^{\left(1\right)}\right)$ be the solution minimizing $\mathit{S}$ over ${\mathit{x}}_{\mathit{B}},{\mathit{y}}_{\mathit{B}},\mathit{R}$ (from system Eq. (21)) for fixed ${\mathit{A}}^{\left(1\right)}$, ${\mathit{C}}^{\left(0\right)}$;

$\left({\mathit{x}\prime }_{\mathit{C}}^{\left(1\right)},{\mathit{y}\prime }_{\mathit{C}}^{\left(1\right)},{\mathit{R}}_3^{\left(1\right)}\right)$ be the solution minimizing $\mathit{S}$ over ${\mathit{x}\prime }_{\mathit{C}},{\mathit{y}\prime }_{\mathit{C}},\mathit{R}$ (from system Eq. (24)) for fixed ${\mathit{A}}^{\left(1\right)}$, ${\mathit{B}}^{\left(1\right)}$.

Then the following inequality holds:

$${\mathit{S}}_3^{\left(1\right)}\le {\mathit{S}}_2^{\left(1\right)}\le {\mathit{S}}_1^{\left(1\right)}.$$

Proof of Theorem 1. Stage 1: For fixed initial values $\left({\mathbf{x}}_{\mathbf{B}}^{\left(0\right)},{\mathbf{y}}_{\mathbf{B}}^{\left(0\right)}\right)$, $\left({\mathbf{x}\prime }_{\mathbf{C}}^{\left(0\right)},{\mathbf{y}\prime }_{\mathbf{C}}^{\left(0\right)}\right)$, system (17) has a unique solution $\left({\mathbf{X}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{Y}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{R}}_1^{\left(1\right)}\right)$ minimizing $\mathbf{S}$ over ${\mathbf{X}}_{\mathbf{A}},{\mathbf{Y}}_{\mathbf{A}},\mathbf{R}$. Hence:

$${\mathbf{S}}_1^{\left(1\right)}=\mathbf{m}\mathbf{i}\mathbf{n}\mathbf{S}\left[{\mathbf{X}}_{\mathbf{A}},{\mathbf{Y}}_{\mathbf{A}},\mathbf{R},{\mathbf{x}}_{\mathbf{B}}^{\left(0\right)},{\mathbf{y}}_{\mathbf{B}}^{\left(0\right)},{\mathbf{x}\prime }_{\mathbf{C}}^{\left(0\right)},{\mathbf{y}\prime }_{\mathbf{C}}^{\left(0\right)}\right]=\mathbf{S}\left[{\mathbf{X}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{Y}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{R}}_1^{\left(1\right)},{\mathbf{x}}_{\mathbf{B}}^{\left(0\right)},{\mathbf{y}}_{\mathbf{B}}^{\left(0\right)},{\mathbf{x}\prime }_{\mathbf{C}}^{\left(0\right)},{\mathbf{y}\prime }_{\mathbf{C}}^{\left(0\right)}\right]$$

Stage 2: Fixing $\left({\mathbf{X}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{Y}}_{\mathbf{A}}^{\left(1\right)}\right)$ and $\left({\mathbf{x}\prime }_{\mathbf{C}}^{\left(0\right)},{\mathbf{y}\prime }_{\mathbf{C}}^{\left(0\right)}\right)$, system Eq. (21) gives a unique solution $\left({\mathbf{x}}_{\mathbf{B}}^{\left(1\right)},{\mathbf{y}}_{\mathbf{B}}^{\left(1\right)},{\mathbf{R}}_2^{\left(1\right)}\right)$ minimizing $\mathbf{S}$ over ${\mathbf{x}}_{\mathbf{B}},{\mathbf{y}}_{\mathbf{B}},\mathbf{R}$:

$${\mathbf{S}}_2^{\left(1\right)}=\mathbf{m}\mathbf{i}\mathbf{n}\mathbf{S}\left[{\mathbf{X}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{Y}}_{\mathbf{A}}^{\left(1\right)},\mathbf{R},{\mathbf{x}}_{\mathbf{B}},{\mathbf{y}}_{\mathbf{B}},{\mathbf{x}\prime }_{\mathbf{C}}^{\left(0\right)},{\mathbf{y}\prime }_{\mathbf{C}}^{\left(0\right)}\right]=\mathbf{S}\left[{\mathbf{X}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{Y}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{R}}_2^{\left(1\right)},{\mathbf{x}}_{\mathbf{B}}^{\left(1\right)},{\mathbf{y}}_{\mathbf{B}}^{\left(1\right)},{\mathbf{x}\prime }_{\mathbf{C}}^{\left(0\right)},{\mathbf{y}\prime }_{\mathbf{C}}^{\left(0\right)}\right]$$

Since this is a minimization over a larger space than in ${\mathbf{S}}_1^{\left(1\right)}$ (i.e., more variables allowed to vary), it follows that:

$${\mathit{S}}_2^{\left(1\right)}\le {\mathit{S}}_1^{\left(1\right)}$$

Stage 3: Fixing $\left({\mathbf{X}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{Y}}_{\mathbf{A}}^{\left(1\right)}\right)$ and $\left({\mathbf{x}}_{\mathbf{B}}^{\left(1\right)},{\mathbf{y}}_{\mathbf{B}}^{\left(1\right)}\right)$, system Eq. (24) provides a unique solution $\left({\mathit{x}\prime }_{\mathit{C}}^{\left(1\right)},{\mathit{y}\prime }_{\mathit{C}}^{\left(1\right)},{\mathit{R}}_3^{\left(1\right)}\right)$ minimizing $\mathbf{S}$ over ${\mathit{x}\prime }_{\mathit{C}},{\mathit{y}\prime }_{\mathit{C}},\mathit{R}$:

$${\mathit{S}}_3^{\left(1\right)}=\mathit{m}\mathit{i}\mathit{n}\mathit{S}\left[{\mathit{X}}_{\mathit{A}}^{\left(1\right)},{\mathit{Y}}_{\mathit{A}}^{\left(1\right)},\mathit{R},{\mathit{x}}_{\mathit{B}}^{\left(1\right)},{\mathit{y}}_{\mathit{B}}^{\left(1\right)},{\mathit{x}\prime }_{\mathit{C}},{\mathit{y}\prime }_{\mathit{C}}\right]=\mathit{S}\left[{\mathbf{X}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{Y}}_{\mathbf{A}}^{\left(1\right)},{\mathbf{R}}_3^{\left(1\right)},{\mathbf{x}}_{\mathbf{B}}^{\left(1\right)},{\mathbf{y}}_{\mathbf{B}}^{\left(1\right)},{\mathbf{x}\prime }_{\mathbf{C}}^{\left(1\right)},{\mathbf{y}\prime }_{\mathbf{C}}^{\left(1\right)}\right]$$

Again, since the minimization is over more variables than in ${\mathbf{S}}_2^{\left(1\right)}$, we have: ${\mathit{S}}_3^{\left(1\right)}\le {\mathit{S}}_2^{\left(1\right)}$. Hence, chaining both inequalities yields the final result: ${\mathit{S}}_3^{\left(1\right)}\le {\mathit{S}}_2^{\left(1\right)}\le {\mathit{S}}_1^{\left(1\right)}$.

This theorem provides the theoretical foundation for a stepwise iterative minimization algorithm:

• Start with initial guesses ${\mathit{B}}^{\left(0\right)}$, ${\mathit{C}}^{\left(0\right)}$;
• Solve system Eq. (17) → get ${\mathit{A}}^{\left(1\right)}$, ${\mathbf{R}}_1^{\left(1\right)}$;
• Solve system Eq. (21) → get ${\mathit{B}}^{\left(1\right)}$, ${\mathbf{R}}_2^{\left(1\right)}$;
• Solve system Eq. (24) → get ${\mathit{C}}^{\left(1\right)}$, ${\mathbf{R}}_3^{\left(1\right)}$;
• Repeat if desired for refinement.

This avoids solving the full 7-variable nonlinear system at once and guarantees non-increasing sequence of errors: ${\mathit{S}}^{(\mathit{k}+1)}\le {\mathit{S}}^{\left(\mathit{k}\right)}$.

5. Algorithm for Minimizing the Function S

This algorithm is shown in Figure 4 to iteratively minimize the weighted difference function $\mathit{S}$, determining the optimal positions of the hinges $\mathit{A}\in \mathit{Q}$, $\mathit{B}\in {\mathit{Q}}_1$ and $\mathit{C}\in {\mathit{Q}}_2$ for constructing a four-link initial kinematic chain ABCD with rotational pairs.

• Set initial position
  Select arbitrary initial points ${\mathit{B}}^{\left(0\right)}\left({\mathbf{x}}_{\mathbf{B}}^{\left(0\right)},{\mathbf{y}}_{\mathbf{B}}^{\left(0\right)}\right)\in {\mathit{Q}}_1$;
  Select arbitrary initial points ${\mathit{C}}^{\left(0\right)}\left({\mathit{x}\prime }_{\mathit{C}}^{\left(0\right)},{\mathit{y}\prime }_{\mathit{C}}^{\left(0\right)}\right)\in {\mathit{Q}}_2$.
• Verify that the determinant $\mathit{D}\ne 0$ (i.e., the linear system Eq. (17) is solvable).
  Main iterative process
  Repeat the following steps for $\mathit{i}=\mathrm{0,1},2,\dots ,:$
  Step 1: Update A
• Solve system Eq. (17) using ${\mathit{B}}^{\left(\mathit{i}\right)}$, ${\mathit{C}}^{\left(\mathit{i}\right)}$ to compute:

  $${\mathit{A}}^{(\mathit{i}+1)}\left({\mathit{X}}_{\mathit{A}}^{(\mathit{i}+1)},{\mathit{Y}}_{\mathit{A}}^{(\mathit{i}+1)}\right), {\mathit{R}}_1^{(\mathit{i}+1)}$$
• Ensure ${\mathit{D}}_1\ne 0$.
  Step 2: Update B
• Solve system Eq. (21) using ${\mathit{A}}^{(\mathit{i}+1)}$, ${\mathit{C}}^{\left(\mathit{i}\right)}$ to compute:

  $${\mathit{B}}^{(\mathit{i}+1)}\left({\mathit{x}}_{\mathit{B}}^{(\mathit{i}+1)},{\mathit{y}}_{\mathit{B}}^{(\mathit{i}+1)}\right), {\mathit{R}}_2^{(\mathit{i}+1)}$$
• Ensure ${\mathit{D}}_2\ne 0$.
  Step 3: Update C
• Solve system Eq. (24) using ${\mathit{A}}^{(\mathit{i}+1)}$, ${\mathit{B}}^{(\mathit{i}+1)}$ to compute:

  $${\mathit{C}}^{(\mathit{i}+1)}\left({\mathit{x}\prime }_{\mathit{C}}^{(\mathit{i}+1)},{\mathit{y}\prime }_{\mathit{C}}^{(\mathit{i}+1)}\right), {\mathit{R}}_3^{(\mathit{i}+1)}$$
• Ensure ${\mathit{D}}_3\ne 0$.

Stopping Criterion

The iterative process continues until the following conditions are met:

$$\left|{\mathit{R}}^{\left(\mathit{k}\right)}-{\mathit{R}}^{(\mathit{k}-1)}\right|\le \mathit{\epsilon };$$

$$\left|{\mathit{X}}_{\mathit{A}}^{\left(\mathit{k}\right)}-{\mathit{X}}_{\mathit{A}}^{(\mathit{k}-1)}\right|\le \mathit{\epsilon };$$

$$\left|{\mathit{Y}}_{\mathit{A}}^{\left(\mathit{k}\right)}-{\mathit{Y}}_{\mathit{A}}^{(\mathit{k}-1)}\right|\le \mathit{\epsilon }.$$

where $\mathit{\epsilon }$ is a specified small tolerance value.

Convergence Guarantee

This algorithm guarantees a non-increasing sequence of function values:

$${\mathit{S}}_1^{\left(1\right)}, {\mathit{S}}_2^{\left(1\right)}, {\mathit{S}}_3^{\left(1\right)}, {\mathit{S}}_1^{\left(2\right)}, {\mathit{S}}_2^{\left(2\right)}, {\mathit{S}}_3^{\left(2\right)},\dots \to {\mathit{S}}^{\mathit{*}}$$

Converging to a local minimum ${\mathit{S}}^{\mathit{*}}$ of the function $\mathit{S}$, as supported by the previously stated theorem.

6. Results

The iterative synthesis algorithm for four-link kinematic chains exhibits a clear pattern of structured convergence, as evidenced by the progressive reduction in the objective function $S$ over successive synthesis cycles, as can be seen in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15 of the algorithm [21]. The graph highlights the sequence $S_1^{\left(i\right)}$, $S_2^{\left(i\right)}$, $S_3^{\left(i\right)}$ for each cycle $i$, corresponding to the staged minimization process: solving for the base point $A$, the hinge point $B$, and the endpoint $C$, respectively. Each cycle is plotted in a distinct color, with annotated percentage reductions showing the relative improvement in the objective function after each step. The following observations are noteworthy:

The first synthesis cycle shows significant improvement, with reductions of: 25% from ${\mathit{S}}_1^{\left(1\right)}=0.20$ to ${\mathit{S}}_2^{\left(1\right)}=0.15$, 33.3% from ${\mathit{S}}_2^{\left(1\right)}=0.15$ to ${\mathit{S}}_3^{\left(1\right)}=0.10$. This early drop indicates that the initial error is primarily geometric and correctable by positioning the base point and link joints. As the algorithm progresses, the percentage reductions become smaller: In the 3rd cycle: 2.4% and 1.9%, In the 4th cycle: less than 0.1%. This behavior confirms that the algorithm approaches a stable local minimum, consistent with the convergence conditions $\left|{\mathit{R}}^{\left(\mathit{k}\right)}-{\mathit{R}}^{(\mathit{k}-1)}\right|\le \mathit{\epsilon }$; $\left|{\mathit{X}}_{\mathit{A}}^{\left(\mathit{k}\right)}-{\mathit{X}}_{\mathit{A}}^{(\mathit{k}-1)}\right|\le \mathit{\epsilon }$; and $\left|{\mathit{Y}}_{\mathit{A}}^{\left(\mathit{k}\right)}-{\mathit{Y}}_{\mathit{A}}^{(\mathit{k}-1)}\right|\le \mathit{\epsilon }$. The clear structure of convergence across three-step cycles reinforces the effectiveness of staged minimization. Optimizing over subsets of variables (first $\mathit{A}$, then $\mathit{B}$, then $\mathit{C}$) ensures that each variable is refined with respect to the updated configuration of the others, yielding a more stable and targeted convergence path [22]. The curve smoothly approaches the limiting value ${\mathit{S}}^{\mathit{*}}\approx 0.040$, with no oscillations or reversals. This is a strong indication of numerical stability, proper problem conditioning (as previously verified with ${\mathit{D}}_1\ne 0$, ${\mathit{D}}_2\ne 0$, ${\mathit{D}}_3\ne 0$, and convergence robustness.

The grouped analysis of the convergence curve demonstrates that the synthesis algorithm is both effective and efficient [23]. The structured approach to optimization ensures consistent improvement in each cycle, while the diminishing percentage reductions serve as a natural indicator of when convergence has been achieved. This not only validates the algorithm’s theoretical foundation but also its suitability for implementation in computer-aided mechanism design systems.

7. Discussion

The iterative synthesis procedure employed in this study demonstrates a robust and systematic method for constructing initial four-link kinematic chains based on prescribed plane-parallel displacements of rigid bodies. The algorithm, grounded in approximate kinematic geometry and staged minimization of the weighted difference function $\mathit{S}$, provides significant insights into the geometric and numerical properties of linkage configuration. One of the most significant observations is the consistent and rapid convergence of the objective function $\mathit{S}$ across successive iterations. The sequence ${\mathit{S}}_1^{\left(1\right)}$, ${\mathit{S}}_2^{\left(1\right)}$, ${\mathit{S}}_3^{\left(1\right)}$, ${\mathit{S}}_1^{\left(2\right)}$, ${\mathit{S}}_2^{\left(2\right)}$, ${\mathit{S}}_3^{\left(2\right)},\dots {\mathit{S}}_{\mathit{k}}^{\left(\mathit{i}\right)}\to {\mathit{S}}^{\mathit{*}}$ reveals a monotonic decrease in error, validating the effectiveness of the superposition and transformation approach. The convergence criteria involving changes in radius $\mathit{R}$, and base point coordinates ${\mathit{X}}_{\mathit{A}},{\mathit{Y}}_{\mathit{A}}$ fall within the prescribed tolerance $\mathit{\epsilon }$ after only a few cycles, underscoring the numerical stability of the method. Additionally, the determinant evaluations ${\mathit{D}}_1$, ${\mathit{D}}_2$, and ${\mathit{D}}_3$ for the linear systems corresponding to subproblems (solving for points $\mathit{A}$, $\mathit{B}$, and $\mathit{C}$) were consistently non-zero. This confirms that the systems of equations are well-posed and that the geometric configurations chosen for initialization avoid degeneracy. These determinant checks serve as an important validation step, ensuring that the matrices involved are invertible and that each subproblem yields a unique solution. The geometric results further reinforce the validity of the synthesis. Points ${\mathit{B}}^{\left(\mathit{i}\right)}$ and ${\mathit{C}}^{\left(\mathit{i}\right)}$, when plotted relative to the updated point ${\mathit{A}}^{(\mathit{i}+1)}$, fall on arcs of circles with radius ${\mathit{R}}^{(\mathit{i}+1)}$, maintaining the expected constraint geometry of a four-bar linkage with revolute joints. These visual confirmations complement the numerical convergence and indicate that the synthesized linkage can reliably reproduce the desired relative motions between moving planes ${\mathit{Q}}_1$ and ${\mathit{Q}}_2$. Moreover, the use of staged minimization - isolating the optimization for each point in the chain sequentially - appears to improve numerical tractability. Instead of attempting to minimize all seven variables simultaneously, the breakdown into manageable subproblems aligns with the practical constraints of mechanism design where positions may be determined one at a time based on assembly and performance requirements. In summary, the results support the conclusion that the proposed method is not only computationally efficient but also geometrically accurate for the synthesis of initial kinematic chains [24,25]. The methodology can be extended to more complex mechanisms or generalized to systems with additional constraints, such as limited joint ranges or dynamic considerations.

8. Conclusions

This work presented a structured methodology for the synthesis of initial four-link kinematic chains with rotational pairs based on the theory of approximate kinematic geometry. The approach involved the specification of relative plane-parallel motion between two moving planes with respect to a fixed base, and the solution of the synthesis problem was carried out by minimizing a weighted geometric deviation function.

Key conclusions from this study are as follows:

• Stable and convergent iterative algorithm:

The proposed iterative algorithm, which minimizes the objective function $S$ in three stages (for points $A$, $B$, and $C$), converges reliably to a local minimum. The observed convergence behavior is smooth, with successive values of $S$ decreasing monotonically, and convergence criteria based on the changes in $R$, $X_A$, and $Y_A$ are satisfied within a few iterations.

2.

Well-conditioned linear systems:

The synthesis process is supported by solving three linear systems (corresponding to systems Eq. (17), Eq. (21), and Eq. (24)), and determinant testing confirmed that the systems are non-singular $D_1\ne 0$, $D_2\ne 0$, $D_3\ne 0$ for the provided input configurations. This indicates that the algorithm is numerically stable and the synthesis conditions are well-posed.

3.

Geometrical accuracy of link placement:

The positions of the synthesized points $B^{\left(i\right)}$ and $C^{\left(i\right)}$ relative to the base $A^{(i+1)}$ lie on circles of radius $R^{(i+1)}$, verifying the geometrical correctness of the linkage design. The visualizations show that the synthesized four-link mechanism satisfies the constraints of consistent angular displacements and relative orientations.

4.

Modular and extendable framework:

The methodology allows for modular implementation, with each synthesis step being independently solvable and verifiable. This makes it especially suitable for computer-aided mechanism design applications and opens up possibilities for extension to mechanisms with translational joints or higher degrees of freedom.

5.

Effective error minimization strategy:

The use of staged error minimization - optimizing each subset of variables sequentially - reduces computational complexity while maintaining synthesis accuracy. This approach provides a practical and scalable solution for the synthesis of planar mechanisms.

In summary, the developed algorithm offers a reliable, accurate, and computationally efficient technique for the synthesis of initial four-link mechanisms under geometrical motion constraints. It serves as a solid foundation for further developments in mechanism synthesis, including the integration of dynamic constraints, optimization under workspace limitations, and generalization to spatial linkages.