**1. Introduction**

Over the past two decades, the number of biomedical research studies, including the study of the foundations of the pathogenesis of infectious pathology, carried out all over the world, has increased [1]. The level of such research is largely determined by the availability of a certified collection of biological material. This applies equally to the development of countermeasures in the face of the threat to the global community associated with the COVID-19 [2] pandemic. The need to reduce the number of production processes with the use of manual labor is also significant for preventing contamination of equipment and work surfaces by infectious agents and, consequently, minimizing the risk of personnel infection. The more steps that can be performed with high reproducibility, the lower the level of random error there will be in the final data. As robotic hardware and computer control systems have advanced and the price has reduced, more steps in the sample preparation process have become amenable to automation and more end users are able to justify the capital investment relative to the anticipated labor savings [3].

The term aliquot comes from the word "aliquot", which means the exact measured fraction of the sample taken for analysis, which retains the properties of the main sample. Thus, the aliquot process involves dividing the whole sample of biomaterial into subfractions. The main developers of robotic equipment for aliquoting biological samples are Dornier-LTF (Germany) [4], TECAN (Switzerland) [5], and Hamilton [6]. The laboratory equipment they produce is as follows: robotic digging station PIRO manufactured by Dornier-LTF (Germany), robotic laboratory station (Freedom EVO® series) manufactured by TECAN, and automatic dispensing system (robot dispenser) Microlab STAR manufactured by Hamilton. They are designed to perform the same type of technological processes

**Citation:** Malyshev, D.; Rybak, L.; Carbone, G.; Semenenko, T.; Nozdracheva, A. Optimal Design of a Parallel Manipulator for Aliquoting of Biomaterials Considering Workspace and Singularity Zones. *Appl. Sci.* **2022**, *12*, 2070. https:// doi.org/10.3390/app12042070

Academic Editor: Alessandro Gasparetto

Received: 10 November 2021 Accepted: 10 February 2022 Published: 16 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and include the main elements: a robotic module with dosing multichannel pipettes and removable tips; a robotic module with a gripper for transferring racks; and auxiliary devices (pumps, barcode readers, thermoshakers, UV emitters, detectors of the position of objects of the working surface, etc.), the number of which depends on the purpose of using the robotic station. These automated stations have a typical architecture and successfully solve the problems of aliquoting biomaterial using tubes of various sizes. For liquid dispensing, intelligent systems are used that record the moments when the pipette touches the liquid surface, the absence of liquid in the test tube, or the presence of colloidal clots that impede pipetting.

However, their significant disadvantage is the inability to aliquot serum from test tubes with different liquid phases. For the operation of the equipment, only samples that have undergone primary manual processing in the form of separating serum from a blood clot can be used.

At the same time, it is important to design a robot, taking into account the optimization of its parameters in order to achieve the best combination of design parameters. Optimal design of robotic systems is an actual and important topic [7–12].

#### **2. Building a Model of a Robotic System**

The work [13] shows the effectiveness of automated aliquoting compared to manual labor. Research in the field of application of robotic systems for liquid dosing and automation of laboratory processes for sample preparation is carried out by many scientists. Active research is being conducted at the Institute of Automation and Center for Life Science Automation, University of Rostock (Rostock, Germany). In works [14,15], the use of a cognitive two-handed robot, working in conjunction with an operator and remotely controlled, is proposed. Pipetting of liquid is carried out with an automatic hand-type pipette, widely used in laboratory practice. Each robot arm has 7 degrees of freedom and can perform any manipulations within the workspace. The disadvantage of this solution is the low speed of the end-effector and the low accuracy of its positioning. So, in [16], an integrated system based on robots is considered for automating a multistage system of laboratory research, including the movement of equipment and biomaterial samples. The possibility of integrating devices controlled through electronic interfaces such as RS232, USB, Firewire, or Ethernet into a common network is shown, which allows centralized control of the technological process.

The article discusses a new architecture and design of a robotic system, which must meet the following requirements to ensure the technological process of aliquoting:


**Figure 1.** Round bottom blood tubes.

The flow chart of the aliquoting process is shown in Figure 2.

**Figure 2.** Aliquoting process.

Figure 3 shows a 3D model of a robotic system, which includes a parallel DeLi manipulator with 4 degrees of freedom, made based on the Delta mechanism and equipped with an excavating head for performing basic operations for aliquoting biomaterial, and a collaborative robot (Uni) with a serial structure for performing transportation operations, racks with test tubes and plates, and microtubes. Uni robot is used for the stage when liquid is divided into fractions. The DeLi robot should take only one upper fraction from a test tube and distribute it into smaller tubes, while not allowing the lower fraction (sediment) to enter. After taking the desired fraction, the DeLi robot can move at any speed available

to this type of robot. The amount of the withdrawn fraction is determined using the developed technical vision system, which is not described in this work. The DeLi manipulator dispenses only one fraction of blood; therefore, the velocity for the DeLi manipulator is not limited. The Delta parallel robot topology has been chosen due to the fact that it has a better performance in terms of dynamics and accuracy as compared with serial robots. Furthermore, Delta architectures are better suited for moving along curvilinear trajectories that may arise in the proposed specific application due to different heights of test tubes and racks. This is justified by the features of the robot structure, which has low inertia and increased characteristics of speeds and accelerations, which are important when moving along circular trajectories, including those involving reverse movement of drives. This makes the Delta architecture better suited also in comparison with portal architecture.

**Figure 3.** Three-dimensional model of the robotic system: 1—body, 2—DeLi manipulator, 3— dispensing head, 4—dispenser tip, 5—robot manipulator, 6—base of the robot manipulator, 7— workspace, 8—rack with test tubes, 9—tray for consumables, 10—trolley.

It is necessary to determine the boundaries of its workspace to select the geometric and design parameters of the robotic system. Determination of the workspace of parallel robots is much more complex than for serial robots, due to the peculiarities of structure, kinematics, and dynamics. There are various methods for determining the workspace of translational robots. A multiobjective optimum design procedure to 3 degrees of freedom (DOF) parallel robots with regards to three optimality criteria: workspace boundary, transmission quality index, and stiffness, is presented in article [17]. A kinematic optimization was performed to maximize the workspace of the parallel robot. Genetic algorithms were applied to optimize the objective function. Article [18] discusses the kinematic and geometric aspects of the 3-DOF translational orthoglide robot. New solutions are proposed to solve the inverse and forward kinematics and conduct a detailed analysis of the workspace and features, taking into account specific joint limit constraints. In [19], such methods of optimizing the workspace as the genetic algorithm and the maximum surrounding workspace are considered. In [20], the cylindrical algebraic decomposition is mentioned, which is used to calculate the robot's workspace in the projection space (x, y, z), taking into account some joint constraints. In [21], an approach to calculating the workspace of a parallel robot Delta is proposed based on the forward kinematics model. In paper [22], the proposed method

is implemented for a spatial 3-DOF parallel robot, known as the Tripteron. The mechanical interference, including interference of links, interference of links with obstacles, and interference of the end-effector with obstacles, are investigated by using new geometrical reasoning. For this purpose, a new geometric method is proposed which is based on the segmen<sup>t</sup> to segmen<sup>t</sup> intersection test. This method can be well extended to a wide range of robotic mechanical systems, including, among others, parallel robots. This method is also applied in [23] for the Delta robot. However, the method considered in this article has disadvantages; in particular, the authors propose to determine the interference of segments on the auxiliary plane and not the distance between the nearest points. This makes it impossible to identify such interference of the links, in which there is no interference of the axes. The workspace for a 3-PUU Delta mechanism with translational drive joints is considered in the article [24]. The analysis of the achievable workspace can be performed using the inverse kinematics method; the workspace is expressed with three-dimensional dot clouds of points. The article [25] presents a workspace, a shared space, and features of a family of Delta-like parallel robots using algebraic tools. In articles [26,27], attention is paid to the analysis of the influence of singularity zones and interference of links. However, the articles did not consider the various configurations of the Delta robot or the influence of singularity zones on the amount of workspace.

Although manipulators of the Delta type are well studied, the development of new methods for analyzing the workspace that are more efficient in terms of increasing the accuracy of approximation and versatility of their application for various versions of their architecture is an important task. The task of determining the workspace can be solved with a given accuracy using interval analysis methods [28–30]. One of the disadvantages of deterministic interval analysis methods is a large number of elements in the covering set. Nevertheless, the article's authors proposed an approach to the transformation of the covering set, which was considered earlier in [31]. It is used to transform the covering set before rendering the stage. The interval analysis method is applied to approximate the set of solutions to a system of nonlinear inequalities that determine the constraints on the parameters of the DeLi robot links (Figure 4). The manipulator has 4 degrees of freedom and includes three RUU kinematic chains. In each of the chains, a rotary drive joint *Ai* is used to connect to the base, and there are two universal joints— *Ci* to connect to the working platform and *Bi* to connect two links together. The base and work platform are equilateral triangles. As an end-effector, we will consider point P—the center of the working platform.

We write the boundaries of the workspace in the form of a system of inequalities:

$$\begin{cases} \theta\_i - \theta\_{\max} \le 0 \\ \theta\_{\min} - \theta\_i \le 0 \end{cases} \tag{1}$$

where *θi* are the angles of rotation in the drive joints, and *θmin* and *θmax* are the minimum and maximum angles, respectively. Based on the solution of the inverse kinematics described in [32], we obtain

$$\theta\_{\bar{i}} = 2 \tan^{-1} \left( \frac{-F\_{\bar{i}} \pm \sqrt{E\_{\bar{i}}^2 + F\_{\bar{i}}^2 + G\_{\bar{i}}^2}}{G\_{\bar{i}} - E\_{\bar{i}}} \right) \tag{2}$$

where *Ei*, *Fi*, and *Gi* are defined as:

$$E\_1 = 2d\left(y\_P + \frac{a - 2c}{2\sqrt{3}}\right), F\_1 = 2z\_P d\,, \tag{3}$$

$$\mathbf{G}\_1 = \mathbf{x}\_P^2 + y\_P^2 + z\_P^2 + \left(\frac{a - 2c}{2\sqrt{3}}\right)^2 + d^2 + 2y\_P \left(\frac{a - 2c}{2\sqrt{3}}\right) - e^2 \tag{4}$$

$$E\_2 = -d\left(\sqrt{3}\left(x\_P + \frac{2c - a}{4}\right) + y\_P + \frac{2c - a}{4\sqrt{3}}\right), \; F\_2 = 2z\_P d,\tag{5}$$

$$G\_2 = x\_P^2 + y\_P^2 + z\_P^2 + \left(\frac{2c - a}{4}\right)^2 + \left(\frac{2c - a}{4\sqrt{3}}\right)^2 + d^2 + 2\left(x\_P\left(\frac{2c - a}{4}\right) + y\_P\left(\frac{2c - a}{4\sqrt{3}}\right)\right) - c^2 \tag{6}$$

$$E\_3 = d\left(\sqrt{3}\left(xp - \frac{2c - a}{4}\right) - yp - \frac{2c - a}{4\sqrt{3}}\right), \text{  $P\_3 = 2zpd, $ \tag{7}$$

$$\mathbf{G}\_3 = \mathbf{x}\_P^2 + y\_P^2 + z\_P^2 + \left(\frac{2c - a}{4}\right)^2 + \left(\frac{2c - a}{4\sqrt{3}}\right)^2 + d^2 + 2\left(y\_P\left(\frac{2c - a}{4\sqrt{3}}\right) - \mathbf{x}\_P\left(\frac{2c - a}{4}\right)\right) - c^2 \tag{8}$$

where *a* is the side length of a regular triangle of the base, *c* is the side length of a regular triangle of the movable platform, *d* is the distance between *A* and *B* joints, *e* is the distance between *B* and *C* joints, and *xP*, *yP*, *zP* are the coordinates of point P.

**Figure 4.** The structure of the DeLi manipulator.

There is an ambiguity in the solution of the inverse kinematics, when for the same coordinates of the point P there are two possible angles of rotation in the drive rotary joint for each of the three kinematic chains. An example of two possible positions of the kinematic chain *AiBiCi* is shown in Figure 5. The first position corresponds to [−] in Equation (2), the second position corresponds to [+].

Thus, there are eight possible solutions to the inverse kinematics, and, therefore, for each of them, a workspace can be determined. In this case, one should take into account the presence of singularity zones of the mechanism when hit, in which the dynamic loads on the links increase significantly, and the robot loses controllability. We consider the method proposed in [33,34], based on the analysis of the Jacobi matrix, to determine the singularity zones. The determinant of the Jacobi matrix has the form:

$$\det(f\_A) = \begin{bmatrix} \frac{\partial \Theta\_1}{\partial x\_P} & \frac{\partial \Theta\_1}{\partial y\_P} & \frac{\partial \Theta\_1}{\partial z\_P} \\ \frac{\partial \Theta\_2}{\partial x\_P} & \frac{\partial \Theta\_2}{\partial y\_P} & \frac{\partial \Theta\_2}{\partial z\_P} \\ \frac{\partial \Theta\_3}{\partial x\_P} & \frac{\partial \Theta\_3}{\partial y\_P} & \frac{\partial \Theta\_3}{\partial z\_P} \end{bmatrix},\tag{9}$$

where Θ*i* are determined by Equation (2).

**Figure 5.** The ambiguity of the solution of the inverse kinematics.

In view of the cumbersomeness of the formulas obtained for each of the elements of the determinant, we present only the first of them:

$$\frac{\partial \Theta\_1}{\partial \mathbf{x}\_P} = \frac{2\left(\frac{2\mathbf{x}\_p \left(2d\mathbf{z}\_p + s\_1\right)}{\left(s\_2 - 2d\left(y\_r - \left(\frac{2c-s}{2\sqrt{3}}\right)\right)\right)^2} + \frac{2\mathbf{x}\_p s\_2}{\mathbf{x}10\left(s\_2 - 2d\left(y\_r - \left(\frac{2c-s}{2\sqrt{3}}\right)\right)\right)}\right)}{\left(\frac{\left(2d\mathbf{z}\_p + s\_1\right)}{2\mathbf{x}\_p s\_2}\right)^2 + 1}\,\tag{10}$$

where

$$s\_1s\_1 = \sqrt{4d^2z\_p^2 + 4d^2\left(y\_p - \left(\frac{2c-a}{2\sqrt{3}}\right)\right)^2 - s\_2^2} \\ s\_2 = \left(\frac{2c-a}{2\sqrt{3}}\right)^2 - 2y\_p\left(\frac{2c-a}{2\sqrt{3}}\right) - c^2 + d^2 + x\_p^2 + y\_p^2 + z\_p^2$$

We write down the condition for the presence of singularity zones:

$$\det(f\_A) = 0\tag{11}$$

Each of the eight solutions of the inverse kinematics also corresponds to eight possible Jacobi matrices and, accordingly, their determinants. It is necessary to ensure the constant sign of the determinant of the Jacobi matrix to exclude singularity zones from the workspace of the robot. To ensure this condition, it is necessary to add one of the conditions to inequality system (1): det(*JA*) < 0 or det(*JA*) > 0, depending on the sign of the determinant. So, singularity surfaces separate the workspace into non-overlapping volumes. It is then necessary for the robot to remain within the same volume as it moves, in order to avoid crossing a singularity.

#### **3. Synthesis of an Algorithm for Determining the Workspace**

Taking into account Equations (1)–(5), an algorithm for approximating the workspace of the DeLi robot is synthesized. The workspace analysis method is a direct application of interval analysis methods, as described in Refs. [28–31].

The algorithm works with a system of inequalities written in a general form:

$$\begin{cases} & g\_1(\mathbf{x}) \le 0, \\ & \cdots, \\ & g\_m(\mathbf{x}) \le 0 \\ & a\_i \le x\_i \le b\_{i\prime} \ i = 1, \dots, n. \end{cases} \tag{12}$$

The initial box *Q*, which includes the entire set of solutions X, is determined by the constraints of the intervals *ai* ≤ *xi* ≤ *bi*, *i* = 1, ... , *n*. Consider an arbitrary box *B*. Let m

*m*(*B*) = max *i*=1,...,*<sup>m</sup>* min *x*∈*B gi*(*x*) and *M*(*B*) = max *i*=1,...,*<sup>m</sup>* max *x*∈*B gi*(*x*). If *m*(*B*) > 0, then B contains no possible points for system (12). The proposed algorithm excludes such boxes. If *M* (*B*) ≤ 0, then each point of the box is a possible solution. Therefore, it can be added to the coverage as an inner box. If a box cannot be excluded, it is divided into two smaller boxes if its diameter is not less than the specified precision *δ*. (Figure 6). The check is iteratively repeated for boxes formed after division until their size becomes less than precision *δ*.

**Figure 6.** Division of the original box *Q*: *Q*2 —into 2 parts, *Q*4 —into 4 parts, *Q*8 —into 8 parts.

The algorithm for determining the workspace, taking into account singularity zones, is similar and differs in a large number of inequalities in the system, as well as additional conditions when checking *m (B)* and *M (B)*, taking into account the strict inequality for the determinant of the Jacobian matrix.

Equation (11) is used instead of the system of inequalities in the algorithm for determining the singularity zones. To find the set of its solutions, we calculate *M (B)*, which is calculated as *M*(*B*) = min *i*∈1,*k* max *x*∈*B gj*(*x*), the condition for adding is that there is no inner coverage,andtheconditionforexcludingboxeswilltaketheform:

 *M*(*B*) < 0 V *m*(*B*) > 0 (13)

We use synthesized algorithms for numerical simulation. A program was written in the C++ programming language to implement the developed algorithms.

#### **4. Simulation Results of Workspace Taking into Account Singularity Zones**

The results of modeling the workspace were obtained for the following parameters of the DeLi manipulator: *a* = 450 mm, *c* = 200 mm, *d* = 150 mm, *e* = 230 mm. The visualization of the results of modeling the workspace is carried out by converting the covering set describing the workspace into a universal format of 3D models—an STL file. The constructed workspace without taking into account singularity zones is shown in Figure 7. Figure 7a shows the workspace in full, 7b—in section by the XOZ plane. In order to increase the performance, the processor multithreading is applied using the OpenMP package. The computation time for the approximation accuracy *δ* = 2 mm and the grid dimension for performing computations of the 32 × 32 × 32 functions was 67 s for each inverse kinematics case.

**Figure 7.** (**a**) Volumetric workspace, (**b**) in section in the XOZ plane.

The results of modeling the workspace taking into account singularity zones are shown in Figure 8. Figure 8a shows the regions corresponding to the singularity zones when using the "+" sign in inverse kinematics (2), Figure 8b—when using the "−" sign, Figure 8c— with all eight possible combinations of the "+" and "−" signs in the inverse kinematics, Figure 8d—for all combinations within the workspace. The average computation time for an approximation accuracy *δ* = 2 mm and a grid dimension of 64 × 64 × 64 on a personal computer was 32 s.

**Figure 8.** Areas corresponding to singularity zones: (**a**) for "+" in inverse kinematics, (**b**) for "−", (**c**) for all possible combinations, (**d**) inside the workspace.

The simulation results show the following conclusions. The inner part of the workspace has a positive sign of the determinant of the Jacobi matrix if two or three kinematic chains have a "−" sign in the inverse kinematics. The negative sign of the determinant is obtained

in the opposite cases, that is, if the "+" sign is selected for at least two chains. With this in mind, when determining the workspace, taking into account singularity zones for four possible solutions, we add the condition of the positive determinant, i.e., det (*JA*) > 0, and for the other four, the condition of negativity, i.e., det (*JA*) < 0. Table1showsthesimulationresultsforthesecases.

**Configuration No. Sign in Inverse Kinematics Sign of the Jacobi Determinant Volume of the Workspace,** *×***10<sup>4</sup> mm<sup>3</sup> Volume Reduction in%** θ1 θ2 θ3 **With Singularity Zones Without Singularity Zones** 1 −−− − 5140 4822 6.18 2 − − + − 5034 4285 14.88 3 − + − − 5034 4285 14.88 4 − + + + 4931 4705 4.58 5 + −− − 5034 4285 14.88 6 + − + + 4931 4705 4.58 7 ++ − + 4931 4705 4.58 8 + + + + 4831 4449 7.91

**Table 1.** Volumes of workspaces taking into account singularity zones.

As can be seen from the table, out of eight options for the inverse kinematics, four groups of solutions can be distinguished: 1—all kinematic chains with "+", 2—all with " −", 3—one chain with "+" and two with "−", and also 4—two with "+" and one with " −". The first two groups correspond to one option, the second two to three options each. Considering the symmetry of the robot structure, the volume of the resulting workspace is practically the same for all variants of the group. This is confirmed by the calculated values of the volume of the workspace. The table also shows that options 2, 3, and 5, as well as 4, 6, and 7 correspond to the same volume of the workspace. Figure 9 shows the workspaces corresponding to each of the groups.

In Figure 10, it can be seen that taking into account the singularity zones led to a decrease in the size of the workspace by 4.58–14.88%, depending on the combination of possible solutions to the inverse kinematics. In addition to singularity zones, the shape and volume of the workspace are directly influenced by the interferences of the robot links with the platforms and with each other.

**Figure 9.** Workspaces, taking into account singularity zones: (**a**)—for option 1, (**b**)—for 2, 3, and 5, (**c**)—for 4, 6, and 7, (**d**)—for 8.

**Figure 10.** A kinematic scheme with highlight of joint rotation angles for a limb of the proposed mechanism.

#### **5. Determination of the Interference of the Links of the Mechanism**

The interferences of the links of the mechanism can be divided into three groups:


The first group can be determined, taking into account the restrictions on the angles of rotation in the joints *Bi* and *Ci* (Figure 10):

$$\left\{ \begin{array}{c} \alpha\_{i} \in \left[ \mathfrak{a}\_{\text{min}}; \alpha\_{\text{max}} \right] \\ \beta\_{i} \in \left[ \beta\_{\text{min}}; \beta\_{\text{max}} \right] \\ \gamma\_{i}^{(j)} \in \left[ \gamma\_{\text{min}}; \gamma\_{\text{max}} \right] \end{array} \right\} \tag{14}$$

The angles can be determined using the formula for the cosines between vectors:

$$a\_i = \frac{(\mathbf{x}\_{Ai} - \mathbf{x}\_{Bi})(\mathbf{x}\_{Ci} - \mathbf{x}\_{Bi}) + (y\_{Ai} - y\_{Bi})(y\_{Ci} - y\_{Bi}) + (z\_{Ai} - z\_{Bi})(z\_{Ci} - z\_{Bi})}{de} \tag{15}$$

$$\beta\_{i} = \frac{\sqrt{3}((\mathbf{x}\_{Bi} - \mathbf{x}\_{Ci})(\mathbf{x}\_{P} - \mathbf{x}\_{Ci}) + (y\_{Bi} - y\_{Ci})(y\_{P} - y\_{Ci}) + (z\_{Bi} - z\_{Ci})(z\_{P} - z\_{Ci}))}{\infty} \tag{16}$$

$$\gamma\_i^{(\text{Bi1})} = \frac{(\mathbf{x}\_{\text{Ci1}} - \mathbf{x}\_{\text{Bi1}})(\mathbf{x}\_{\text{Bi2}} - \mathbf{x}\_{\text{Bi1}}) + (y\_{\text{Ci1}} - y\_{\text{Bi1}})(y\_{\text{Bi2}} - y\_{\text{Bi1}}) + (z\_{\text{Ci1}} - z\_{\text{Bi1}})(z\_{\text{Bi2}} - z\_{\text{Bi1}})}{h\epsilon} \tag{17}$$

$$
\gamma\_i^{(B2)} = \gamma\_i^{(C1)} = -\gamma\_i^{(C2)} = -\gamma\_i^{(Bi1)} \tag{18}
$$

$$\begin{split} \text{where } & \mathbf{x}\_{A1} = z\_{A1} = z\_{A2} = z\_{A3} = \mathbf{x}\_{B1} = 0, \mathbf{x}\_{A2} = \frac{a}{4}, \mathbf{x}\_{A3} = -\frac{a}{4}, \mathbf{y}\_{A1} = -\frac{a}{2\sqrt{3}}, \mathbf{y}\_{A2} = \mathbf{y}\_{A3} = 0, \mathbf{x}\_{B1} = -\frac{a}{4\sqrt{3}}, \mathbf{y}\_{B2} = \mathbf{x}\_{B3} = \mathbf{y}\_{B2} = \mathbf{x}\_{B3} = \mathbf{y}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B1} = \mathbf{x}\_{B2} = \mathbf{x}\_{B2} = \mathbf{x}\_{B$$

The second group of interference includes possible interference of links with platforms *A*1, *A*2, *A*3 and *C*1, *C*2, *C*3. In this case, the interference with the platform *C*1, *C*2, *C*3 is excluded by condition (14). The interference of links *Ai*, *Bi* and *Bi*1, *Bi*2 with platform *A*1, *A*2, *A*3 is also excluded by condition (1). Let us add the conditions for the appearance of the remaining interference, namely links *Bi*1, *Ci*1 and *Bi*2, *Ci*2 with platform *A*1, *A*2, *A*3 (Figure 11). The first condition for the occurrence of interference is:

$$z\_{Bi} \ge 0\tag{19}$$

If condition (19) is satisfied for the joint *Bi*, then we calculate the coordinate points of the points *Di*1 and *Di*2, lying on the segments *Bi*1, *Ci*1 and *Bi*2, *Ci*2, respectively, which can be the interference points of the segments with the platform *A*1, *A*2, *A*3:

$$z\_{\rm Di1} = z\_{\rm Di2} = 0 \tag{20}$$

$$\mathbf{x}\_{Di1} = \mathbf{x}\_{Bi1} + \frac{-z\_{Bi1}(\mathbf{x}\_{Ci1} - \mathbf{x}\_{Bi1})}{z\_{Ci1} - z\_{Bi1}},\\ y\_{Di1} = y\_{Bi1} + \frac{-z\_{Bi1}(y\_{Ci1} - y\_{Bi1})}{z\_{Ci1} - z\_{Bi1}}\tag{21}$$

$$\mathbf{x}\_{\rm Di2} = \mathbf{x}\_{\rm Bi2} + \frac{-z\_{\rm Bi2}(\mathbf{x}\_{\rm CI2} - \mathbf{x}\_{\rm Bi2})}{z\_{\rm Ci2} - z\_{\rm Bi2}},\\\ y\_{\rm Di2} = y\_{\rm Bi2} + \frac{-z\_{\rm Bi2}(y\_{\rm Ci2} - y\_{\rm Bi2})}{z\_{\rm Ci2} - z\_{\rm Bi2}}\tag{22}$$

**Figure 11.** Interference of links with a fixed platform of the mechanism.

The platform is an equilateral triangle with points *A*1, *A*2, *A*3 being the midpoints of the sides. The coordinates of the vertices *<sup>A</sup>*1, *<sup>A</sup>*2, *A*3 are defined as:

$$z\_{A\_1'} = z\_{A\_2'} = z\_{A\_3'} = \mathbf{x}\_{A\_3'} = 0,\ \mathbf{x}\_{A\_1'} = -\frac{a}{2},\ \mathbf{x}\_{A\_2'} = \frac{a}{2},\ \mathbf{y}\_{A\_1'} = \mathbf{y}\_{A\_2'} = -\frac{a}{2\sqrt{3}},\tag{23}$$

$$y\_{A\_3'} = \frac{a}{\sqrt{3}}\tag{24}$$

As a result, it is required to check in the two-dimensional plane XOY whether the points *Di*1 and *Di*2 belong to the triangle *<sup>A</sup>*1*A*2*A*3. If at least one of the points is included in the triangle, then there is interference with the platform.

We define the third group of interference using an approach based on determining the minimum distance between the segments drawn between the centers of the joints of each of the links. In [22,23], a similar condition is used, but the approach has drawbacks. In particular, the authors propose to determine the interference of the segments on the auxiliary plane, and not the distance between the nearest points. This does not allow identifying such an interference of links in which there is no interference of the axes.

The approach proposed in the current work is as follows. We construct an auxiliary plane to determine the interference of the links of the mechanism. This plane is parallel to the axis of one of the links and to which the axis of the other link belongs. In this case, the condition for the absence of interference of the links will take the form:

$$
\sqrt{\mathbf{u}\_1^2 + \mathbf{u}\_2^2} > \mathbf{r}\_{link1} + \mathbf{r}\_{link2} \tag{25}
$$

where *u1* is the distance between the axis of a link that does not belong to the plane and the auxiliary plane, *u2* is the distance between the nearest points of the segments connecting the centers of the joints of each of the links when projecting a segmen<sup>t</sup> that does not belong to the auxiliary plane onto this plane, and <sup>r</sup>*link*1, <sup>r</sup>*link*2 are the radiuses of the links.

It is worth noting the special case when the links are parallel to each other and the construction of an auxiliary plane is not required.

The algorithm for determining the interference of links is as follows:



A more detailed approach to the definition of the third group of interference is presented in [35].

The proposed approach can be similarly applied to determine the interference between the links of the Uni and DeLi robots. In this case, it is possible to discretely separate the trajectories of the robots in the process of collaborative work of two robots. Then, it is required to determine the positions of the joints for each of the positions of the robots. Using the coordinates of the joints, one must check the intersection of all combinations of links using condition (25) and the above algorithm.

#### **6. Analysis of the Workspace Taking into Account the Interference of the Links**

The coverage obtained at the stage of defining the workspace is used to determine the positions of the end-effector inside the workspace, in which the interference of the links occurs. The coverage is divided into many boxes of equal size less than the approximation accuracy.

For the coordinate of the center of each of the boxes, conditions (14) and (25) are checked. If all the conditions are met, it is entered into a new set. The resulting set will contain boxes that correspond to the workspace without interference.

Workspace for configuration 1 according to Table 1 after excluding interference areas, as well as interference areas for *a* = 450 mm, *c* = 200 mm, *d* = 150 mm, *e* = 230 mm, *αmin* = *βmin* = *γmin* = 10◦, *αmax* = *βmax* = *γmax* = 350◦, D*link* = 10 mm are shown in Figure 12. Table 2 shows the change in the volume of the workspace, taking into account the interference, without interference of links, as well as volume reduction for various configurations by taking into account the interference.


**Table 2.** Volumes of workspaces taking into account link interference.

**Figure 12.** Workspace of the mechanism after excluding the areas of interference: (**a**) for configuration 1, (**b**) 2,3,5, (**c**) 4,6,7, (**d**) 8.

As can be seen from Table 2, with an increase in the number of kinematic chains, which corresponds to "+" in the inverse kinematics, the greater the percentage of the workspace that is excluded due to interference. This is due to the fact that in this case the kinematic chains are "bent" inward. As a result, the volume of the workspace for the first configuration is 3.05 times greater than the volume for the eighth configuration.

Visualization of the positions of the links at which they intersect occurs to verify the results of determining the interference. Figure 13 shows some of the link crossings that occur.

Figure 12 shows that the main inner part of the work area is round. Simulations were also performed for other link sizes (Table 3). The simulation results are shown in Figure 14. As can be seen from the figures, the workspace for various sizes has a round shape.

**Table 3.** Volumes of workspaces taking into account link interference.


**Figure 13.** Examples of identified link interference: (**a**) links B31C31 and B32C32 with a fixed platform and link A3B3 with link B32C32, (**b**) links B31C31 and B32C32 with a fixed platform, (**c**) link A2B2 with link B22C22.

**Figure 14.** Workspace of the mechanism after excluding the areas of interference: (**a**) for configuration 1, (**b**) 2, (**c**) 3.

#### **7. Geometric Parameters Optimization**

Next, we will optimize the geometric parameters of the mechanism. The volume of the workspace, taking into account singularity zones and the interference of links, is an optimization criterion. Determination of the optimal sizes of links a, c, d, and e was carried out in several stages. At each stage, the volume of the workspace was calculated taking into account the interference of the links for various combinations of sizes. The ranges of resizing and the iteration step are reduced with each stage to reduce the computational complexity.

Numerical values are given in Table 4. The last column contains the volume of the workspace, which is a criterion for excluding or including in the next stage a certain part of the size range. Constant dimensions for modeling: *h* = 100 mm, D*link* = 30 mm. The sum of the dimensions *a*, *c*, *d*, and *e* was taken to be constant and equal to 1850 mm. The amount was selected in accordance with the dimensions of the industrial robot ABB IRB 360: *a* ≈ 600 mm, *c* ≈ 100 mm, *d* ≈ 350 mm, *e* ≈ 100 mm. It should be noted that at the first stage, the computations were carried out for all configurations; however, the condition for the required volume of the workspace was fulfilled only for configurations 1 and 2/3/5. At the second stage, the volume condition was met only for the first configuration. In this regard, modeling was carried out only for the first configuration starting from the third stage.


**Table 4.** Stages of determining the optimal link lengths.

The maximum workspace is reached at *a* = 61 mm, *c* = 86 mm, *d* = 513 mm, and *e* = 1190 mm. For such a ratio of dimensions, the total volume of the workspace, taking into account singularity zones, is 4.597 m3, of which 2.714 m<sup>3</sup> is an area without interference of links. Thus, assuming the conditional size of the side of the base *a* = 1, the ratios of the sizes *c*, *d*, and *e* were 1.41, 8.41, and 19.51, respectively.

The dependence of the volume of the workspace for various configurations on the change in the parameters *d* for *a* = 61 mm and *c* = 86 mm, *e* = (1850 − *a* − *c* − *d*) is shown in Figure 15. The numerical values of the workspace volume for the plot are obtained by iterative execution of the developed program for various parameters.

**Figure 15.** Dependence of the volume of the workspace on the change in the length of the links *d* and *e*.

Figure 16 shows the dependence of the proportion of areas in which interference occurs in the total volume of the workspace depending on the parameters *d* and *e*. The graphs show that for the first configuration, the fraction of interference decreases with increasing size *d*. For the rest of the configurations, the percentage of interference for almost all sizes *d* is more than 90%.

Let us perform the second stage of optimization. It takes into account the requirements of the aliquoting process. The required workspace height *hc* is based on the requirements of tube and pipette sizes (350 mm). We can conclude that the manipulator workspace is round based on the previous analysis. Based on this, an optimization problem can be posed, which maximizes the radius *rc* of a cylinder of a given height *hc*, which can be inscribed into the manipulator's workspace. In this case, the sum of the dimensions *a*, *c*, *d*, and *e* of the manipulator should also be minimized. Therefore, the general criterion k can be written as:

$$k = \frac{a+c+d+e}{r\_c} \tag{26}$$

**Figure 16.** Dependence of the ratio of areas of interference to the workspace on the length of the links.

The optimization was performed in several stages, similar to the optimization performed above to maximize the workspace. Initial data for optimization: *a* ∈ [30; <sup>1500</sup>], *c* ∈ [30; <sup>1500</sup>], *d* ∈ [30; <sup>1500</sup>], *e* ∈ [30; <sup>1500</sup>], *hc* = 350 mm. As a result, the optimal value of the criterion *k* = 1.778 was obtained for the following dimensions: *a* = 61 mm, *c* = 87 mm, *d* = 478 mm, *e* = 1485 mm, *rc* = 1187, 29 mm. The simulation of the workspace was performed for the obtained dimensions. Figure 17 shows that the workspace completely includes the cylinder with the calculated radius.

**Figure 17.** Workspace for optimal parameters.

An additional simulation was performed to compare the results of the two optimizations. The volume of the workspace and the maximum cylinder diameter were calculated for both optimal results. The parameters obtained during the first optimization were changed while maintaining the ratio 1:1.41:8.41:19.51. The change is necessary due to the incorrectness of comparing results that have a different sum of parameters. The sum of the parameters for the second optimization was 2111 mm. Therefore, the dimensions of the first optimization were changed so that their sum was 2111 mm. The obtained initial data and simulation results are shown in Table 5.


**Table 5.** Stages of determining the optimal link lengths.

As can be seen from the table, the workspace volume was larger as expected for the parameters obtained as a result of the first optimization and the maximum radius of the inscribed cylinder for the second optimization. However, it should be noted that the workspace volume in the first case turned out to be only 0.47% larger, while the radius for the second case was larger by 9.34%. Therefore, in this case, the radius of the inscribed cylinder is the main criterion for the problem of aliquoting. So, we can conclude that second optimization results should be currently accepted as the final results. However, a formal optimization algorithm will be planned in the future for a further optimization of the obtained results.

The aliquoting process assumes that Uni and DeLi should have a common part of workspaces. The Uni kinematic scheme is shown in Figure 18.

**Figure 18.** The kinematic scheme of Uni robot.

Workspace of Uni is calculated using D-H parameters and transformation matrices, where D-H parameters are equal to *a*2 = 240 mm, *a*3 = 210 mm, *d*1 = 150 mm, *d*4 = 110 mm, *d*5 = 85 mm, and *d*6 = 81 mm. The end-effector positions were calculated in discrete steps using the forward kinematics. The calculated positions were entered into the workspace array. The resulting array was exported to the STL model. Likewise, a program was written in the C++ programming language to implement the developed algorithm. The Uni workspace is shown in Figure 19.

**Figure 19.** The calculated Uni workspace.

Taking into account the calculated workspace of the Uni robot, we perform a simulation of a robotic system. The workspace models of robots were exported to the CAD system. The relative position of the robots was determined taking into account their workspaces and the requirements for the absence of intersections of links of mechanisms. Simulation of collaborative work of robots is shown in Figure 20. The execution of operations should be consistent: the Uni manipulator delivers a rack with test tubes to the joint workspace; then the DeLi manipulator performs the aliquoting process. The simulation confirmed the absence of intersections during collaborative work.

**Figure 20.** CAD model of RS movement with the movement of robots: (**a**) Uni moves the rack into the joint workspace, (**b**) DeLi interacts with the rack within a joint workspace.

Figure 21 shows the intersection of the DeLi and Uni workspaces. The intersection area is sufficient to accommodate a rack with 96 test tubes. Based on this, it can be concluded that the workspace of the DeLi robot and the relative position of the robots meet the requirements of the aliquoting process. The above conditions are met if the robot Uni is fixed relative to the coordinate system DeLi at a point with coordinates (0, 830, −897).

**Figure 21.** The intersection of the workspaces of a Uni and DeLi robot.
