**2. Materials and Methods**

In this paper, a robust design model is presented through a methodological and systematic approach based on the design philosophy proposed by Genichi Taguchi. The integration of optimization processes and ANN design are methodological tools that allow the performance and generalization capacity in ANN models to be improved. In this study, an RDANN methodology was used, which was initially proposed for the reconstruction of spectra in the field of neutron dosimetry by means of ANNs [7].

The RDANN methodology was used to estimate the optimal structural parameters in a BPNN to calculate the inverse kinematics in a six-DOF robot, where the main objective was the development of accurate and robust ANN models. In other words, it was sought that the selection of the structural parameters of the proposed model allows us to obtain the best possible performance in the network.

#### *2.1. Kinematic Analysis*

The robot called Quetzal is based on an open source, 3D-printable, and low-cost manipulator [42]. The modeling and representation were carried out using the Denavit– Hartemberg (D–H) parameters to obtain a kinematic model through four basic transformations that are determined based on the geometric characteristics of the Quetzal manipulator to be analyzed [43]. The D–H parameters are shown in Table 1.


**Table 1.** D–H parameters of the Quetzal manipulator.

The basic transformations represent a sequence of rotations and translations, where the reference system of element *i* is related to the system of element *i* − 1. The transformation matrix is given by Equation (2).

$$\sigma\_i^{i-1} T = R\_Z(\theta\_i) D\_Z(d\_i) D\_{X\_r}(a\_{i-1}) R\_X(a\_{i-1}) \tag{2}$$

Carrying out the multiplication of the four matrices, Equation (3) is obtained:

$$\begin{aligned} \;\_i^{i-1}T = \begin{bmatrix} c\theta\_i & -s\theta\_i c\alpha\_{i-1} & s\theta\_i s\alpha\_{i-1} & a\_{i-1}c\theta\_i \\ s\theta\_i & c\theta\_i c\alpha\_{i-1} & -c\theta\_i s\alpha\_{i-1} & a\_{i-1}s\rho\_i \\ 0 & s\alpha\_{i-1} & c\alpha\_{i-1} & d\_i \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{aligned} \tag{3}$$

where *<sup>i</sup>*−<sup>1</sup> *<sup>i</sup> T* is the D–H transformation matrix from coordinate system *i* to *i* − 1. *RZ*(*θi*) is the rotation matrix representing a rotation *θi* around the *Z* axis, *DZ*(*di*) is the translation matrix representing a translation of *di* on the *Z* axis, *DX*,(*ai*−1) is the translation matrix representing a translation of *ai*−<sup>1</sup> on the *X* axis, *RX*(*αi*−1) is the rotation matrix representing a rotation *αi*−<sup>1</sup> around the *X* axis, *cθ<sup>i</sup>* is shorthand for cos *cos*(*θi*), *sθ<sup>i</sup>* is shorthand for *sin*(*θi*), etc. The transformation matrices of each of the joints are multiplied to obtain the initial position of the end effector in the base reference system, as shown in Equation (4).

$$\begin{aligned} \,^0\_6T &= \,^0\_1A \cdot \,^1\_2A \cdot \,^2\_3A \cdot \,^3\_4A \cdot \,^4\_5A \cdot \,^5\_6A = \begin{bmatrix} ^0\_6R\_{3 \times 3} & ^0\_6P\_{3 \times 1} \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} n\_x & o\_x & a\_x & p\_x \\ n\_y & o\_y & a\_y & p\_y \\ n\_z & o\_z & a\_z & p\_z \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{aligned} \tag{4}$$

Therefore, the equation of the forward kinematics of the Quetzal manipulator can be expressed as shown in Equation (5).

$$F\_{forward\\_k}(\theta\_1, \theta\_2, \theta\_3, \theta\_4, \theta\_5, \theta\_6) = \left(p\_{x\_\tau} p\_{y\_\tau} p\_{z\_\tau} n\_{x\_\tau} n\_{y\_\tau} n\_{z\_\tau} o\_{x\_\tau} o\_{y\_\tau} o\_{z\_\tau} a\_{x\_\tau} a\_{y\_\tau}, a\_z\right) \tag{5}$$

As shown in Equation (5), the position of the end effector of the manipulator can be obtained from the angular values of the six joints of the manipulator. However, in practice, it is necessary to obtain the angles at each of the joints through a given position, so it is necessary to calculate the inverse kinematics, which can be expressed as shown in Equation (6).

$$F\_{\text{invers}\ge \underline{k}}(p\_{\mathbf{x}\prime}p\_{y\prime}p\_{z\prime}n\_{\mathbf{x}\prime}n\_{y\prime}n\_{z\prime}o\_{\mathbf{x}\prime}o\_{y\prime}o\_{z\prime}a\_{\mathbf{x}\prime}a\_{y\prime}a\_{z\prime}) = (\theta\_{1\prime}\theta\_{2\prime}\theta\_{3\prime}\theta\_{4\prime}\theta\_{5\prime}\theta\_{6})\tag{6}$$

Solving (4) gives the orientation and position of the final effector in regard to the reference system, as shown in Equation (7), where the position vector [*p*] = {42.215, 0, 20.2} and the orientation vector [*noa*] = {0, 0, 1, 0, −1, 0, 1, 0, 0}.

$$T\_0^6 = \begin{bmatrix} 0 & 0 & 1 & 42.215 \\ 0 & -1 & 0 & 0 \\ 1 & 0 & 0 & 20.2 \\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{7}$$

The graphic representation of the initial position of the Quetzal robotic manipulator is shown in Figure 2 through a simulation carried out with the Robotics Toolbox for MATLAB software [44].

**Figure 2.** Representation of the initial position of the Quetzal manipulator.

### *2.2. Training and Testing Datasets*

The dataset was generated from the equations obtained in the forward kinematics of the Quetzal manipulator. The variables involved in the proposed dataset were the orientation vector [*noa*] = 5 *nx*, *ny*, *nz*, *ox*, *oy*, *oz*, *ax*, *ay*, *az* 6 , the position vector [*p*] = 5 *px*, *py*, *pz* 6 , and the vector of joint angles [*θ*] = {*θ*1, *θ*2, *θ*3, *θ*4, *θ*5, *θ*6}, for a total of 18 variables.

Table 2 shows the ranges of movement established for each of the joints in the workspace of the manipulator. Dataset *X* was generated with a spatial resolution of 25 × 25 × 25 × 25 × 25 × 25 × 25 in a six-dimensional matrix with 18 variables involved. The total size of the dataset was 4,394,531,250 data, occupying an approximate physical memory space of 32.74 Gb due to the data class of type *Single* − *precision* used in the study, where the data are stored as 4-byte (32-bit) floating point values.

**Table 2.** Range of motion and step size at each of the joints.


Figure 3 shows a graphical representation of the workspace using a 3D data scatter plot corresponding to the position vector [*p*] based only on the joints *θ*1, *θ*2, *θ*4, *θ*<sup>5</sup> and *θ*6, where it is possible to appreciate the workspace of the robotic manipulator without taking into account joint *θ*3. The illustrated workspace was generated from the forward kinematic equations as a function of the six joints.

**Figure 3.** Position vector [*p*] = 5 *px*, *py*, *pz* <sup>6</sup> in function of joints *<sup>θ</sup>*1, *<sup>θ</sup>*2, *<sup>θ</sup>*4, *<sup>θ</sup>*<sup>5</sup> and *<sup>θ</sup>*6.

In order to process the enormous volume of data in a conventional processor, a data reduction filter (DRF) based on linear systematic sampling (LSS) was applied to reduce the set to a size of 190.7 Kb in memory [45]. The data were processed on an 8-core AMD Ryzen 7 5000 series processor with a base clock of 1.8 GHz, 16 GB of RAM, and integrated Radeon 16-thread graphics with a maximum clock of 4.3 GHz [46]. The dataset and code are available in Supplementary Materials.

Figure 4 shows the scatter matrices corresponding to the position dataset before and after applying the FRD filter, where a reduction of 99.99% of the data was obtained with a size of 24,410 data. In Figure 4b, it is shown that the data maintain a constant and uniform distribution with respect to the dataset of Figure 4a.

The data were normalized with a mean of zero in the range from −1 to 1 using Equation (8), where *data* is the data to normalize, *min* is the minimum value of the dataset, *range* is a value established by the difference between the maximum value and the minimum value of the dataset, *h* and *l* are the maximum and the minimum desired values for normalization.

$$DataNorm = \left(\frac{data-min}{range} \times (h-l)\right) + l \tag{8}$$

#### *2.3. Robust Design Methodology*

Figure 5 shows the RDANN methodology based on the Taguchi philosophy that consists of four stages. The designer must know the problem and choose an appropriate network model, as well as the parameters involved in the design of the network for its optimization (planning stage). By implementing an OA and systematically training a small number of ANN models (experimentation stage), the response to be analyzed is determined using the S/N relationship of the Taguchi method (analysis stage). Finally, through a confirmation process, the best performance values of the model are obtained (confirmation stage).

**Figure 4.** Scatter matrix of the position dataset: (**a**) before applying the reduction filter; (**b**) after applying the data reduction filter.

**Figure 5.** Robust design methodology for the optimization of structural parameters.

The graphic representation of the BPNN used in this work is shown in Figure 6, with 12 input variables and 6 output variables that correspond to the position vector [*p*] = 5 *px*, *py*, *pz* 6 and orientation vector [*noa*] = 5 *nx*, *ny*, *nz*, *ox*, *oy*, *oz*, *ax*, *ay*, *az* 6 as input and the vector of joint angles [*θ*] = {*θ*1, *θ*2, *θ*3, *θ*4, *θ*5, *θ*6} as output.

**Figure 6.** BPNN network topology used in this study.

#### 2.3.1. Planning Stage

In this stage, the design variables, noise, and the objective function are identified. The objective function is defined according to the purpose and requirements of the system. In this work, the objective function is related to the prediction or classification errors between the calculator data and the data predicted by the ANN model during the testing stage. The performance at the output of the ANN or the mean square error (MSE) is used as the objective function and is described in the following equation:

$$MSE = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left(\theta\_i^{PREDICED} - \theta\_i^{ORIGINAL}\right)^2} \tag{9}$$

In this case, *N* represents the number of attempts, *θPREDICTED <sup>i</sup>* represents the set of joint values that are predicted by the BPANN, and *θORIGINAL <sup>i</sup>* represents the set of joint values.

The design variables correspond to those that can be controlled by the user, such as the number of neurons in the first layer, the number of neurons in the second layer, the momentum constant, and the learning rate. By contrast, the noise variables are commonly not directly controlled by the user in most cases, such as the initialization of synaptic weights that are generally assigned randomly, the size of the training sets versus test sets, and the random selection of training and test sets. According to the requirements of the problem, the user can choose the factors related to variation in the system during the optimization process. Four design variables and three noise variables were selected because they were directly involved with the performance of the ANN, as described below in Table 3 with their respective configuration levels.

In terms of variables, *A* is the number of neurons in the first hidden layer; *B* is the number of neurons in the second hidden layer; *C* is the momentum constant, which allows the stabilization of the updating of each of the synaptic weights taking into account the sign of the gradient; *D* is the learning rate, which allows us to define the cost that the gradient has in updating a weight because the increase or decrease in the synaptic weight is related to the magnitude of the proposed value, so it may or may not affect the convergence of the MSE, causing instability and divergence [47]. *X* is the initial set of weights, *Y* is the size in proportions of the dataset, and *Z* is the random selection from the training and testing set.


**Table 3.** Design and noise variables.

Once the variables and their respective levels were chosen, a suitable OA was chosen to carry out the training sessions. An OA is described as *Lr*(*Sc*), where *r* represents the number of rows, *c* represents the number of columns, and *s* represents the number of levels in each of the columns. In this experiment, the columns of the OA represent the parameters to be optimized, and the rows represent the tests carried out by combining the three proposed levels.

#### 2.3.2. Experimentation Stage

The success in this stage depends on an adequate choice of the OA because, in this process, a series of calculations are carried out in order to evaluate the interaction and the effects produced between the variables involved through a reduced number of experiments. For the implementation of a robust design, Taguchi suggests the use of a configuration in two crossed OAs with *L*<sup>9</sup> 34 and *L*<sup>4</sup> 32 , as shown below in Table 4.



#### 2.3.3. Analysis Stage

Through the S/N ratio, a quantitative evaluation is carried out, where the mean and the variation in the responses measured by the ANN with different design parameters are considered. The unit of measure is the decibel, and the formula is described as follows:

$$S/N = 10 \cdot \log\_{10}(MSD) \tag{10}$$

In this case, *MSD* is the root mean square deviation in the ANN performance. The best topology is considered when more signals and less noise are obtained; therefore, a high S/N ratio at this stage allows us to identify the best design values in the BPANN with the help of statistical analysis with the JMP software.

#### 2.3.4. Confirmation Stage

In this stage, the value of the robustness measure is obtained based on the specifications and optimal conditions of the design. A confirmation experiment is carried out using the optimal design conditions that were previously chosen, in order to verify if the calculated value is close to the value predicted by the BPANN.

#### **3. Results**

In this work, the RDANN methodology was used for the optimal selection of the structural parameters in a feed-forward backpropagation network, known as BPNN, to find the solution to the inverse kinematics in a Quetzal robot. For the BPNN training, the "resilient backpropagation" training algorithm and *mse* = 1*E*−<sup>4</sup> were selected. In accordance with the RDANN methodology, an OA corresponding to the design and noise variables, respectively, was implemented in configurations *L*<sup>9</sup> 34 and *L*<sup>4</sup> 32 to determine the response to the tests during the 36 training sessions carried out.

The results obtained after applying the RDANN methodology are presented in the next sections.

#### *3.1. Planning Stage*

Table 5 shows the design and noise variables with their respective assigned values for the different levels proposed during the experiment.


**Table 5.** Design and noise variables with their assigned levels.

The values for the three levels established in each of the tests regarding the number of neurons for the first hidden layer were *A* = 80, 100, and 120, respectively; for the number of neurons in the second hidden layer, they were *B* = 30, 60, and 90, respectively; for the constant momentum, they were *C* = 0.1, 0.2, and 0.3, respectively; and for the learning rate, they were *D* = 0.01, 0.1, and 0.2, respectively.

The values for the initial sets of weights *X* were randomly determined at all three levels. The values set in the proportions of the dataset for level 1 were *Y* = 70% training and 30% testing; for level 2, they were 80% and 20%, and for level 3, they were 90% and 30%, respectively; finally, the random selection of the training and testing set for level 1 was *Z* = *Training*1/*Test*; for level 2, it was *Training*2/*Test*2, and for level 3, it was *Training*3/*Test*3.
