Multi-objective optimization design of parallel manipulators using a neural network and principal component analysis

In this work, a multi-objective optimization design method is proposed based on principal component analysis (PCA) and a neural network to obtain a mechanism's optimal comprehensive performance. First, multi-objective optimization mathematical modeling, including design parameters, objective functions, and constraint functions, is established. Second, the sample data are obtained through the design of the experiment (DOE) and are then standardized to eliminate the adverse effects of a non-uniform dimension of objective functions. Third, the first k principal components are established for p performance indices (k<p) using the variance-based PCA method, and then the factor analysis method is employed to define its physical meaning. Fourth, the overall comprehensive performance evaluation index is established by objectively determining weight factors. Finally, the computational cost of the modeling is improved by combining the neural network and a particle swarm optimization (PSO) algorithm. Dimensional synthesis of a Sprint (3RPS) parallel manipulator (PM) is taken as a case study to implement the proposed method, and the optimization results are verified by a comprehensive performance comparison of robots before and after optimization.

1 Introduction

Parallel manipulators (PMs) have become a research hotspot because of their excellent dynamic performance, simple inverse kinematics, and high stiffness (Yang et al., 2022). Sprint (3PRS) (Chen et al., 2014), Tricept (3UPS-UP) (Farooq et al., 2021), and delta (Wang et al., 2014) robots have been successfully commercialized, and the 3RPS PM is used as a parallel module in the five-axis machining center for blade surface machining (Arabshahi and Novinzadeh, 2015). Comprehensive performance of PMs is crucial for machining precision. Therefore, multi-objective optimization is necessary to achieve optimal comprehensive performance.

PMs have smaller workspaces compared with serial robots. Therefore, maximizing the workspace volume, especially the regular workspace volume, is a major task of a PM's optimization design (Yang et al., 2019). However, poor performance, e.g., kinematic, stiffness, and dynamic performance, may be caused while the workspace volume is maximized (Shi et al., 2014). This problem can be fixed by considering other performance indices as constraint functions or by carrying out a multi-objective optimization design.

Compared with kinematic analysis, dynamic analysis is more complicated because the overall stiffness and mass matrices are difficult to describe due to the multiple closed-loop characteristics of PMs (Shi et al., 2014). Limited discrete-node function values in the workspace are adopted to evaluate the objective performance (global performance index, GPI) due to the difficulty in obtaining its analytical expressions in the workspace. It significantly increases the performance index evaluation's computational cost (Lian et al., 2017). Therefore, effective multi-objective optimization design is challenging.

Since multiple performance indices are often presented as competitive coupling relations, this makes it difficult to obtain an optimal solution with a definite physical meaning. The comprehensive objective method is most commonly used to fix this problem, and it defines the weighted sum of objective functions as a comprehensive index (Wang and Zhang, 2017; Zhang and Nelson, 2011). Xu et al. (2017) optimized the hybrid mechanism's comprehensive performance through a comprehensive objective method wherein the workspace volume and stiffness performance were considered as objective functions. Huang et al. (2018) adopted this method to optimize reconfigurable a PM's workspace volume, global dexterity, and stiffness indices. Various combinations of weighting factors are provided to decision-makers. However, the weight coefficients in this method are often determined artificially, and the dimension and unit of the objective function are inconsistent, which leads to the limitation of the application of this method (Fan et al., 2022).

Yang et al. (2019) derived a globally optimal comprehensive performance by proposing a multi-objective optimization game algorithm. Globally optimal solutions are derived based on their weight factors rather than an objective function weight sum. The method is simple but needs to manually set the weight factors and calculate the objective function's values at each discrete node in advance. Discrete node density decides the precision. As the density of discrete nodes increases, the calculation cost increases exponentially. For instance, for a PM with 3 degrees of freedom (DOFs), a population n_p, and generations n_g, if the node density in the workspace is increased to 2 times that of the original, the computational cost will be n_pn_g2³ times that of the original (Yang et al., 2019; Wu, 2017).

The Pareto front method has become a research hotspot of multi-objective optimal design methods because it can obtain multiple sets of non-dominated solutions (Altuzarra et al., 2011; Ling et al., 2022). Cui et al. (2015) combined a genetic algorithm (GA) and insight software to optimize the kinematics and kinetic flexibility of the tunnel tube assembly system's 3-DOF PMs and obtained multiple sets of non-dominated solutions. Qi et al. (2018), Sun and Lian (2018), and Sun et al. (2019) defined the sum of minimum squares (i.e., the cooperative balance point) from the non-dominated solutions as the global optimal solution. However, the physical explanation for this method is not clearly defined.

Establishing workspaces' analytic formulas is challenging, especially for regular workspaces, so the GPIs often obtained in terms of numerical methods include the Monte Carlo method (Wang et al., 2015) and the discrete-node method (Xu et al., 2018). High computational costs greatly reduce the optimized design's computational efficiency. The response surface model (RSM), including a back-propagation (BP) neural network (Zhang and Gao, 2008; Gao et al., 2010), multivariate regression (Yang et al., 2022; Lian et al., 2019), and Gaussian regression (Vinh et al., 2019), is often used to fix this problem, wherein the BP neural network has attracted extensive attention due to its strong learning and nonlinear fitting capabilities. Therefore, the BP neural network model is adopted to establish the RSM of the GPIs and the design parameters in this work.

The main contribution of this work is to propose a multi-objective optimization design method in terms of factor analysis, BP neural network, principal component analysis (PCA), and particle swarm optimization (PSO) algorithm (Sun and Lian, 2018; Wang et al., 2017) to achieve the optimal comprehensive performance of the mechanism. The proposed method has a clear physical explanation and a low computational expense.

The remainder of the work is arranged as follows. The structure description of the 3RPS PM is presented in Sect. 2. Objective functions of the 3RPS PM are presented in Sect. 3. The multi-objective optimization design of the 3RPS PM is introduced in Sect. 4, and the optimization results are discussed in detail. Finally, conclusions are drawn in Sect. 5.

2 Procedure of the multi-objective optimization design method

The basic idea of the multi-objective optimization design for the PMs proposed in this work is to first use PCA to construct an overall comprehensive performance evaluation index with objectively determined weight factors and clear physical significance. Then a response surface model between objective functions and design parameters is constructed based on a neural network to solve the problem of the high computational cost of the objective functions in the workspace. Figure 1 shows the procedure of the proposed multi-objective optimization design method of PMs based on a neural network and PCA. The steps are described as follows.

Figure 1 Procedure of the multi-objective optimization design method of PMs based on the neural network and PCA.

Step 1. Establish the multi-objective optimization mathematical model, including design parameters, objective functions, and constraint functions, wherein the limitations of the design parameters, actuation lengths, and joint angles should be given.

Step 2. Obtain the sample data based on the design of the experiment: Latin hypercube sampling (LHS) is recommended in this work, and then it is necessary to standardize the original sample data to eliminate adverse influences caused by non-uniform dimensions of objective functions. The Z-score method is adopted in this work.

Step 3. Establish the correlation matrix of multiple performance indices and the first k principal components according to the cumulative contribution rate (the ratio of the sum of the first k eigenvalues of the correlation matrix to the sum of all the eigenvalues). The i principal components can be established based on the corresponding eigenvectors of the correlation matrix. Finally, the overall comprehensive performance index is defined as the sum of the product of the first k eigenvalues and the corresponding principal components.

Step 4. Establish the RSM between GPIs and design parameters based on the neural network to avoid the high computational cost of the objective functions in the workspace.

Step 5. Combine the overall comprehensive performance index and RSMs to carry out the optimization design of PMs to achieve optimal design parameters with a clear physical explanation.

3 Structure description of the 3RPS PM

Figure 2 shows a spatial 3RPS PM with 3 DOFs. Three identical branches connect the base and mobile platform at points B_i and A_i with revolute joints and spherical joints, respectively. An upper rod and a lower rod connected by a prismatic joint constitute limbs. Inputs to the mechanism refer to the three prismatic joints. Coordinate frames $B_i\left\{x_i,y_i,z_i\right\}$, $o\left\{x,y,z\right\}$, and $O\left\{X,Y,Z\right\}$ are attached to the limb, moving platform, and base, respectively (see Fig. 2). Both ΔA₁A₂A₃ and ΔB₁B₂B₃ are equilateral triangles with $\left|o{A}_{\mathrm{1}}\right|=\left|o{A}_{\mathrm{2}}\right|=\left|o{A}_{\mathrm{3}}\right|=r_{\mathrm{2}}$ and $\left|O{B}_{\mathrm{1}}\right|=\left|O{B}_{\mathrm{2}}\right|=\left|O{B}_{\mathrm{3}}\right|=r_{\mathrm{1}}$. Each of the revolute joint axes is perpendicular to vector OB_i and lies on the X–Y plane. Geometric and physical parameters before optimization are given as moving platform height h=0.05 m, r₁=1.7 m, elasticity modulus E=200 GPa, r₂=0.8 m, shear modulus G=77 GPa, and d_i=0.4 m (diameter of the link). The 3RPS PM's multi-objective optimization aims to select appropriate r₁ and r₂ values to achieve maximal comprehensive performance.

Figure 2 3RPS PM.

4 Objective functions of the 3RPS PM

The 3RPS PM is to be used as the parallel module of the high-speed machining equipment, wherein dynamic performance, motion or force transmission performance, and operation workspaces strongly influence the mechanism's comprehensive performance. Therefore, the regular workspace volume, kinematic performance, and natural frequency are considered the three objective functions in this paper.

4.1 Workspace evaluation

A PM's workspace is smaller than those of serial mechanisms. Thus, PM optimization should primarily maximize workspace volume. A regular workspace volume (Yang et al., 2019) is recommended as the workspace evaluation index due to some performance indices being poor at the boundary of the reachable workspace, and the reachable workspace does not apply to trajectory planning (Babu et al., 2016). The method proposed in Yang et al. (2019) is adopted to calculate the regular workspace volume, the reachable workspace is first classified into n layers with an equal interval height of d_z, and then the inscribed circle of the workspace boundary in each layer is calculated. The 3RPS PM's regular workspace volume is expressed as follows:

$$\begin{array}{}\text{(1)}& V_{\mathrm{g}}=\sum _{i=\mathrm{1}}^m{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(A_i+A_{i+\mathrm{1}}\right)\times d_z,\end{array}$$

where $A_i=\mathrm{2}\mathit{\pi }{\mathit{\rho }}_i^{\mathrm{2}}$ is the inscribed circle's area in the ith layer. ρ_i is the radius. The minimum radius of inscribed circles is set to 50 mm in the work to guarantee the necessary operation workspace of the 3RPS PM.

4.2 Kinematic performance evaluation

The kinematic condition number is often adopted as the kinematic index of PMs to measure their dexterity, manipulability, and singularity. The main limitation of the condition number index is the non-uniform dimension of the Jacobian matrix for the PMs with mixed DOFs. Characteristic length is often adopted to fix this problem. However, due to the subjectivity of characteristic length selection and the fact that only some elements of the Jacobian matrix are involved in characteristic length calculation, the physical explanation for the characteristic length is ambiguous (Zang et al., 2019; Wang et al., 2019). The motion or force transmission index (Wu et al., 2023; Liu et al., 2012) is recommended in this paper due to its advantages of being frame-free and dimensionless.

The ith limb's twist system is derived through screw theory:

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}{\mathit{S}}_{i\mathrm{1}}=\left[\begin{array}{cc}{\mathbf{R}}_i{\mathit{e}}_{\mathrm{2}}& \mathit{O}{\mathit{B}}_i\times {\mathbf{R}}_i{\mathit{e}}_{\mathrm{2}}\end{array}\right],{\mathit{S}}_{i\mathrm{2}}=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{l}}_i\end{array}\right],\\ {\mathit{S}}_{i\mathrm{3}}=\left[\begin{array}{cc}{\mathit{e}}_{\mathrm{1}}& \mathit{O}{\mathit{A}}_i\times {\mathit{e}}_{\mathrm{1}}\end{array}\right],{\mathit{S}}_{i\mathrm{4}}=\left[\begin{array}{cc}{\mathit{e}}_{\mathrm{2}}& \mathit{O}{\mathit{A}}_i\times {\mathit{e}}_{\mathrm{2}}\end{array}\right],\\ {\mathit{S}}_{i\mathrm{5}}=\left[\begin{array}{cc}{\mathit{e}}_{\mathrm{3}}& \mathit{O}{\mathit{A}}_i\times {\mathit{e}}_{\mathrm{3}}\end{array}\right],\end{array}\right.\end{array}$$

where S_ij is the jth twist of the ith limb, l_i is the unit vector along the vector B_iA_i. R_i represents a rotation matrix from limb coordinate frame {B_i} to global coordinate frame {O}. e_i in Eq. (2) is presented in ${\mathit{e}}_{\mathrm{1}}={\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}$, ${\mathit{e}}_{\mathrm{2}}={\left[\begin{array}{ccc}\mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}$, and ${\mathit{e}}_{\mathrm{3}}={\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}$.

The transmission wrench of the ith limb can be obtained based on its reciprocal product with the twist (S_ij) of the passive joints being zero:

$$\begin{array}{}\text{(3)}& {\mathit{S}}_{Ti}=\left[\begin{array}{cc}{\mathit{l}}_i& \mathit{O}{\mathit{B}}_i\times {\mathit{l}}_i\end{array}\right].\end{array}$$

The ith limb's input twist screw is

$$\begin{array}{}\text{(4)}& {\mathit{S}}_{Ai}=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{l}}_i\end{array}\right].\end{array}$$

All of the inputs except the one in the ith limb are locked to derive the instantaneous 1-DOF output twist screw S_Oi:

$$\begin{array}{}\text{(5)}& {\mathit{S}}_{Oi}\circ {\mathit{U}}_i=\mathrm{0},\end{array}$$

where U_i is the mechanism's constraint wrench system.

Accordingly, the ith limb's transmission and input indices are expressed as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\lambda }}_i={\displaystyle \frac{\left|{\mathit{S}}_{Ti}\circ {\mathit{S}}_{Ai}\right|}{{\left|{\mathit{S}}_{Ti}\circ {\mathit{S}}_{Ai}\right|}_{max}}},{\mathit{\eta }}_i={\displaystyle \frac{\left|{\mathit{S}}_{Ti}\circ {\mathit{S}}_{Oi}\right|}{{\left|{\mathit{S}}_{Ti}\circ {\mathit{S}}_{Oi}\right|}_{max}}}.\end{array}$$

Both of the output and input indices are considered to obtain the local transmission index (LTI):

$$\begin{array}{}\text{(7)}& \text{LTI}=min\left\{{\mathit{\lambda }}_i,{\mathit{\eta }}_i\right\}.\end{array}$$

The global transmission index (GTI) is obtained to determine the global kinematic performance in a PM's regular workspace:

$$\begin{array}{}\text{(8)}& \text{GTI}={\displaystyle \frac{{\int }_{V_r}\text{LTI}\mathrm{d}{V}_r}{{\int }_{V_r}\mathrm{d}{V}_r}}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\text{LTI}}_i,\end{array}$$

where N is the number of discrete grid nodes in the prescribed workspace.

4.3 Dynamic performance evaluation

The natural frequency is the frequency at which a system tends to oscillate in the absence of external forces or damping. The first natural frequency, i.e., the fundamental frequency (Ganesh and Rao, 2020), is recommended in this paper due to its comprehensive consideration of the stiffness, mass, and damping of PMs.

The global independent generalized displacement coordinates (IGDCs) (Yang et al., 2021) are recommended as the elastodynamic modeling method of PMs in this paper due to the Lagrangian multipliers and kinematic constraint equations not being required.

The ith link's elastodynamic equation expressed in the limb coordinate frames is given as follows considering that the link is discretized into one element (the number of elements has little influence on the calculation accuracy of the fundamental frequency; Yang et al., 2021):

$$\begin{array}{}\text{(9)}& {\mathbf{M}}_{\mathrm{e}}^L{\ddot{\mathit{u}}}_i+{\mathbf{K}}_{\mathrm{e}}^L{\mathit{u}}_i={}^L{\mathit{f}}_i,\end{array}$$

where ${}^L{\mathit{u}}_i=\left[\begin{array}{cccc}{}^L{\mathit{d}}_{Bi},& {}^L{\mathit{\phi }}_{Bi},& {}^L{\mathit{d}}_{Ai},& {}^L{\mathit{\phi }}_{Ai}\end{array}\right]$ is the ith link's displacement column vector. ^Ld_Bi(^Lφ_Bi) is point B_i's linear displacement coordinate. ^Ld_Ai(^Lφ_Ai) is point A_i's angular displacement coordinate. ${}^L{\mathit{f}}_i={\left[\begin{array}{cccc}{}^L{\mathit{f}}_{Bi}^{\mathrm{T}},& {}^L{\mathit{m}}_{Bi}^{\mathrm{T}},& {}^L{\mathit{f}}_{Ai}^{\mathrm{T}},& {}^L{\mathit{m}}_{Ai}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$, and ^Lf_Bi(^Lm_Bi) and ^Lf_Ai(^Lm_Ai) are the forces (moments) exerted on points B_i and A_i, respectively. K_e is the element stiffness matrix. M_e is the element mass matrix. L in the upper-left corner indicates that the vector is expressed in the local frame system.

Equation (9) can be further presented in the global coordinate system:

$$\begin{array}{}\text{(10)}& {\mathbf{M}}_i{\ddot{\mathit{U}}}_i+{\mathbf{K}}_i{\mathit{U}}_i={\mathit{F}}_i,\end{array}$$

where ${\mathbf{M}}_i={\mathbf{T}}_{\mathbf{i}}{\mathbf{M}}_{\mathrm{e}}{{\mathbf{T}}_{\mathbf{i}}}^{\mathrm{T}}$, ${\mathbf{K}}_i={\mathbf{T}}_{\mathbf{i}}{\mathbf{K}}_{\mathrm{e}}{{\mathbf{T}}_{\mathbf{i}}}^{\mathrm{T}}$, and F_i=T_if_i. ${\mathbf{T}}_{\mathbf{i}}=\text{diag}\left[\begin{array}{cccc}{\mathbf{R}}_{\mathbf{i}},& {\mathbf{R}}_{\mathbf{i}},& {\mathbf{R}}_{\mathbf{i}},& {\mathbf{R}}_{\mathbf{i}}\end{array}\right]$.

With a rigid assumption of the moving platform, its elastodynamic equation in the global coordinate system is

$$\begin{array}{}\text{(11)}& {\mathbf{M}}_{\mathrm{p}}{\ddot{\mathit{U}}}_{\mathrm{p}}={\mathit{F}}_{\mathrm{p}},\end{array}$$

where ${\mathbf{M}}_{\mathrm{p}}={\mathbf{T}}_{\mathbf{p}}{}^L{\mathbf{M}}_{\mathrm{p}}{\mathit{T}}_{\mathrm{p}}^{\mathrm{T}}$. ${\mathbf{T}}_{\mathbf{p}}=\text{diag}\left[\begin{array}{c}\mathbf{R},\mathbf{R}\end{array}\right]$ and R is the rotation matrix from frames {o} to {O}. ^LM_p is the moving platform's mass matrix. ${\mathit{U}}_{\mathrm{p}}={\left[{\mathit{d}}_{\mathrm{p}}^{\mathrm{T}},{\mathit{\phi }}_{\mathrm{p}}^{\mathrm{T}}\right]}^{\mathrm{T}}$, d_p, and φ_p, respectively, denote point o's linear and angular displacement coordinates in the global coordinate frame. ${\mathit{F}}_{\mathrm{p}}={\left[{\mathit{f}}_{\mathrm{p}}^{\mathrm{T}},{\mathit{m}}_{\mathrm{p}}^{\mathrm{T}}\right]}^{\mathrm{T}}$, and m_p and f_p are the moment and force, respectively, exerted on the moving platform.

The moving platform and link B_iA_i are connected at point A_i through a spherical joint. The compatibility equation between point o's displacement coordinates and point A_i's linear displacement coordinates of the ith link can be obtained according to multipoint constraint (MPC) theory (Yang et al., 2019).

$$\begin{array}{}\text{(12)}& {\mathit{d}}_{\mathit{A}\mathit{i}}=\left[\begin{array}{cc}{\mathbf{E}}_{\mathbf{3}}& [{\mathit{A}}_i\mathit{o}\times ]\end{array}\right]{\mathit{U}}_{\mathit{p}}\end{array}$$

Revolute joints are used to connect the base to link B_iA_i based on MPC theory. The boundary conditions in the local coordinate frame are as follows:

$$\begin{array}{}\text{(13)}& \begin{array}{ll}\left\{\begin{array}{l}{}^L{\mathit{d}}_{B_i}={\mathrm{0}}_{\mathrm{3}\times \mathrm{1}}\\ {}^L{\mathit{\phi }}_{B{x}_i}={}^L{\mathit{\phi }}_{B{z}_i}=\mathrm{0}\end{array}\right.& (i=\mathrm{1},\mathrm{2},\mathrm{3}),\end{array}\end{array}$$

where ${}^L{\mathit{\phi }}_{B{x}_i}$ and ${}^L{\mathit{\phi }}_{B{z}_i}$ are point B_i's angular displacement coordinates along the x_i and z_i axes, respectively.

Boundary conditions in Eq. (13) can be further expressed based on Ganesh and Rao (2020).

$$\begin{array}{}\text{(14)}& \begin{array}{ll}\left\{\begin{array}{l}{\mathit{d}}_{B_i}={\mathrm{0}}_{\mathrm{3}\times \mathrm{1}}\\ {\mathit{\phi }}_{\mathit{B}\mathit{x}}=\frac{R_i(\mathrm{1},\mathrm{2})}{R_i(\mathrm{2},\mathrm{2})}{\mathrm{\Phi }}_{By}\\ {\mathit{\phi }}_{\mathit{B}\mathit{z}}=\mathrm{0}\end{array}\right.& (i=\mathrm{1},\mathrm{2},\mathrm{3})\end{array}\end{array}$$

Accordingly, the global IGDC of the manipulator can be extracted as follows.

$$\begin{array}{}\text{(15)}& \mathit{U}={\left[\begin{array}{cccccccc}{\mathit{\phi }}_{By\mathrm{1}}& {{\mathit{\phi }}_{\mathit{A}\mathbf{1}}}^{\mathrm{T}}& {\mathit{\phi }}_{By\mathrm{2}}& {{\mathit{\phi }}_{\mathit{A}\mathbf{2}}}^{\mathrm{T}}& {\mathit{\phi }}_{By\mathrm{3}}& {{\mathit{\phi }}_{\mathit{A}\mathbf{3}}}^{\mathrm{T}}& {{\mathit{d}}_{\mathit{p}}}^{\mathrm{T}}& {{\mathit{\phi }}_{\mathit{p}}}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}\end{array}$$

Finally, the elastodynamic control equation of the 3RPS PM is given by

$$\begin{array}{}\text{(16)}& \mathbf{M}\ddot{\mathit{U}}+\mathbf{K}\mathit{U}=\mathbf{W},\end{array}$$

where $\mathbf{M}=\sum _{i=\mathrm{1}}^{\mathrm{3}}{{\mathbf{N}}_{\mathbf{i}}}^{\mathrm{T}}{\mathbf{M}}_{\mathbf{i}}{\mathbf{N}}_{\mathbf{i}}+{{\mathbf{N}}_{\mathbf{p}}}^{\mathrm{T}}{\mathbf{M}}_{\mathbf{p}}{\mathbf{N}}_{\mathbf{p}}$, $\mathbf{K}=\sum _{i=\mathrm{1}}^{\mathrm{3}}{{\mathbf{N}}_{\mathbf{i}}}^{\mathrm{T}}{\mathbf{K}}_{\mathbf{i}}{\mathbf{N}}_{\mathbf{i}}$, $\mathbf{W}=\sum _{i=\mathrm{1}}^{\mathrm{3}}{{\mathbf{N}}_{\mathbf{i}}}^{\mathrm{T}}{\mathit{F}}_i+{{\mathbf{N}}_{\mathbf{p}}}^{\mathrm{T}}{\mathit{F}}_{\mathrm{p}}$, U_i=N_iU, N_i is the mapping matrix from the generalized displacement coordinate U_i to the global IGDC U.

It can be seen that the Lagrangian multipliers or the simultaneous kinematic constraint equations are not required in Eq. (16): all the kinematic constraints of the mechanism were considered in the global IGDC and the mapping matrix N_i, which is beneficial for solving dynamic equations.

The angular frequency in radian per second is derived by

$$\begin{array}{}\text{(17)}& \text{det}\left(\mathbf{K}-{\mathit{\omega }}_i\mathbf{M}\right)=\mathrm{0}.\end{array}$$

The fundamental frequency is

$$\begin{array}{}\text{(18)}& f_{\mathrm{1}}={\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}}}{\mathrm{2}\mathit{\pi }}}.\end{array}$$

The PM's global dynamic index is

$$\begin{array}{}\text{(19)}& \text{GDI}={\displaystyle \frac{{\int }_{V_{\mathrm{g}}}{f}_{\mathrm{1}}\mathrm{d}{V}_{\mathrm{g}}}{{\int }_{V_{\mathrm{g}}}\mathrm{d}{V}_{\mathrm{g}}}}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{f}_{i\mathrm{1}}.\end{array}$$

5 Multi-objective optimization analysis using the PCA and the neural network

The multi-objective optimization mathematical modeling of the 3RPS PM (Babu et al., 2013) is

$$\begin{array}{}\text{(20)}& \left\{\begin{array}{l}max\left\{V_r,\text{GTI},\text{GDI}\right\},\\ \mathrm{0.5}\le r_{\mathrm{1}}\le \mathrm{2.5},\mathrm{0.5}\le r_{\mathrm{2}}\le \mathrm{2.5},\\ r_{\mathrm{1}}+r_{\mathrm{2}}=\mathrm{2.5},\\ L_{min}\le L_i\le L_{max},{\mathit{\theta }}_i\le {\mathit{\theta }}_{max},{\mathit{\phi }}_i\le {\mathit{\phi }}_{max},\end{array}\right.\end{array}$$

where $L_{min}=\mathrm{1}$ m and $L_{max}=\mathrm{4}$ m. L_i is link B_iA_i's length. φ_i and θ_i represent angles among the moving platform, base, and link, respectively, and ${\mathit{\theta }}_{max}={\mathit{\phi }}_{max}=\mathrm{50}$^∘.

The PCA recombines the m dependent original indices into a new set of p (p<m) independent principal components (F_i) and determines an overall comprehensive index with a deterministic weight factor based on the variance-based method, which satisfies $\text{cov}(F_i,F_j)=\mathrm{0}$. The physical meanings of the principal components are then explained by using the factor analysis method.

Consider the non-uniform dimensions of objective functions: it is necessary to convert original sample data $X(X=[X_{\mathrm{1}},X_{\mathrm{2}},\mathrm{\dots },X_{\mathrm{p}}],X_i=[X_{\mathrm{1}i},X_{\mathrm{2}i},\mathrm{\dots },X_{ni}{]}^{\mathrm{T}})$ to the standardized X^* before calculating the principal components. The Z-score method is adopted here.

$$\begin{array}{}\text{(21)}& {\mathit{X}}^{*}=\left[\begin{array}{cccc}{\mathit{X}}_{\mathrm{1}}^{*}& {\mathit{X}}_{\mathrm{2}}^{*}& \mathrm{\cdots }& {\mathit{X}}_{\mathrm{p}}^{*}\end{array}\right],\end{array}$$

where ${\mathit{X}}_{ij}^{*}=\frac{X_{ij}-{\mathit{\mu }}_j}{{\mathit{\sigma }}_j}$, ${\mathit{X}}_i^{*}={\left[X_{\mathrm{1}i}^{*},X_{\mathrm{2}i}^{*},\mathrm{\dots },X_{ni}^{*}\right]}^{\mathrm{T}}$, and ${\mathit{\sigma }}_j=\sqrt{\text{var}\left({\mathit{X}}_j\right)}$, μ_j=E(X_j).

Population characteristics are often represented by sample data in terms of the LHS method (Lian et al., 2017). Table 1 shows 276 sets of sample data designed based on the LHS method. Objective functions' standardized sample data of the 3RPS PM are presented through Eq. (21) (see Table 2). The correlation matrix R is established to reflect correlations of objective functions:

$$\begin{array}{}\text{(22)}& \mathbf{R}=\left[\begin{array}{cccc}{r}_{\mathrm{11}}& r_{\mathrm{12}}& \mathrm{\cdots }& r_{\mathrm{1}p}\\ \mathrm{⋮}& & & \\ r_{p\mathrm{1}}& r_{p\mathrm{2}}& & \mathrm{\cdots }& r_{pp}\end{array}\right],\end{array}$$

where $r_{ij}=\text{cov}(X_i,X_j)=\frac{\sum _{k=\mathrm{1}}^n(X_{ki}-{\mathit{\mu }}_i)(X_{kj}-{\mathit{\mu }}_j)}{n-\mathrm{1}}$ $(i,j=\mathrm{1},\mathrm{2},\mathrm{\dots },p)$.

Table 1 Objective functions' raw sample data of the 3RPS PM.

Table 2 Objective functions' standardized sample data of the 3RPS PM.

The correlation matrix can be defined in Eq. (23) for standardized data X^* as follows:

$$\begin{array}{}\text{(23)}& \mathit{R}={\displaystyle \frac{{\mathit{X}}^{*\mathrm{T}}{\mathit{X}}^{*}}{n-\mathrm{1}}}.\end{array}$$

The correlation matrix of the objective functions obtained by Eq. (23) is shown in Table 3. Dynamic performance and motion or force transmission performance are negatively correlated with the workspace volume. Dynamic performance and kinematic performance present positive correlations, which is consistent with the results in Chen and Yang (2021).

Table 3 Correlation matrix of objective functions.

The first k principal components' cumulative contribution rate and the ith principal component are defined by

$$\begin{array}{}\text{(24)}& \left\{\begin{array}{l}{F}_i=a_{\mathrm{1}i}{x}_{\mathrm{1}}^{*}+a_{\mathrm{2}i}{x}_{\mathrm{2}}^{*}+\mathrm{\dots }+a_{pi}{x}_{\mathrm{p}}^{*},\\ \mathit{\delta }=\frac{\sum _{i=\mathrm{1}}^k{\mathit{\lambda }}_i}{\sum _{i=\mathrm{1}}^p{\mathit{\lambda }}_i},\end{array}\right.\end{array}$$

where λ_i and $a_i=(a_{\mathrm{1}i},a_{\mathrm{2}i},\mathrm{\dots },\text{and}{a}_{pi}{)}^{\mathrm{T}}$ are the ith eigenvalue and eigenvector of the correlation matrix, respectively. δ is the first k principal components' cumulative contribution rate.

The first two eigenvalues' cumulative contribution rate is 92.65 %, which is greater than 85 % (Kuroda et al., 2011) (see Table 4). Objection function information is included in the first two principal components (see Eq. 25).

$$\begin{array}{}\text{(25)}& \left\{\begin{array}{l}{F}_{\mathrm{1}}=-\mathrm{0.3281}{V}_r^{*}+\mathrm{0.6663}{\text{GTI}}^{*}+\mathrm{0.6697}{\text{GDI}}^{*}\\ F_{\mathrm{2}}=\mathrm{0.9445}{V}_r^{*}+\mathrm{0.2439}{\text{GTI}}^{*}+\mathrm{0.2200}{\text{GDI}}^{*}\end{array}\right.\end{array}$$

The physical explanation of principal components is given by introducing the factor analysis method. The factor-loading matrix is first given by

$$\begin{array}{}\text{(26)}& \mathbf{A}=\left[\begin{array}{ccc}{\mathit{a}}_{\mathrm{1}}\sqrt{{\mathit{\lambda }}_{\mathrm{1}}},& \mathrm{\dots },& {\mathit{a}}_m\sqrt{{\mathit{\lambda }}_m}\end{array}\right].\end{array}$$

The factor-loading matrix should be further simplified to explain principal components until each objective function has a large loading in only one principal component. The maximum orthogonal rotation method of square differences is usually adopted (Mardian et al., 1979):

$$\begin{array}{}\text{(27)}& \mathbf{B}=\mathbf{A}\sum _{i=\mathrm{1}}^{m-\mathrm{1}}\sum _{j=i+\mathrm{1}}^m{\mathbf{\Omega }}_{jk},\end{array}$$

where Ω is the orthogonal matrix.

Table 4 Principal component cumulative contribution rate and correlation matrix eigenvalues and eigenvectors.

Repeat Eq. (28) until the allowable error is larger than the change in B's relative variances.

The kinematic performance and the dynamic performance have larger factor loadings for the first principal component (see Table 5). Therefore, the kinematic and dynamic performance factor is named for the first principal component and presents a PM's dynamic performance and kinematic performance. The workspace volume has a larger factor loading on the second principal component; thus, the workspace factor is named for the second principal component, which reflects the mechanism's workspace volume. So far, the physical explanation of the principal components has been fully defined through the factor analysis method.

Table 5 Factor-loading matrix after orthogonal transformation.

Accordingly, the 3RPS PM's overall comprehensive performance evaluation index is

$$\begin{array}{}\text{(28)}& F=\sum _{i=\mathrm{1}}^k{\mathit{\lambda }}_i{F}_i=\mathrm{0.2250}{V}_r^{*}+\mathrm{1.4747}{\text{GTI}}^{*}+\mathrm{1.4597}{\text{GDI}}^{*}.\end{array}$$

Generally, the GPI is obtained by calculating the variance or mean of the function value of nodes in the regular workspace. Computational costs increase exponentially with increased discrete points. The high computational intensity reduces the optimization design efficiency. The RSM between design parameters and GPIs is established using the BP neural network in this work to improve the computational intensity for its powerful nonlinear fitting effect. The BP neural network's hidden layer sizes are set to 30. The recommended acceptance levels of RMSE and R² in Agarwal and Renaud (2004) are below 0.2 and above 0.9, respectively. Table 6 shows the optimal accuracy assessment for the BP-based mapping model through five tests that verify the accuracy and effectiveness of the BP-based RSM.

Table 6 Accuracy assessment for the BP-based mapping model.

The population- and evolution-based PSO intelligence algorithm determines their search directions based on their velocities and replaces the GA's mutation and crossover. The PSO algorithm's convergence efficiency was proven in Wang et al. (2017). Thus, the PSO algorithm is used in this work. Table 7 shows the PSO algorithm's setting that is recommended in Sun et al. (2018). The 59 iterations are performed to converge the objective function to 3.5945, and the optimal design parameters are r₁=1.4340 and r₂=1.0660 (see Fig. 3). The calculation time of the proposed algorithm is 192 s on a computer 3.00 GHz CPU, while the PSO algorithm is 6285 s and 96.94 % of the computational cost is saved. Optimized dynamic performance, kinematic performance, and workspace volume increase by 44.57 %, 6.74 %, and 283.69 %, respectively; after optimization, the overall comprehensive index F increases by 118.92 % (see Table 8). The multi-objective optimization design results verified the effectiveness of the proposed method.

Table 7 Parameter setting of the PSO optimization algorithm.

Figure 3 Convergence of the overall comprehensive evaluation index using the PSO algorithm.

Table 8 Comparison of performance indices before and after optimizing the 3RPS PM.

Dynamic performance, kinematic performance, and regular workspace have been enhanced as shown in Table 8. The distributions of the performance indices are shown in Figs. 4–5, the maximum radius of the workspace has been increased from 0.096 to 0.156 m, the maximum LTI has been increased from 0.95 to 0.99, and the maximum first natural frequency has been increased from 28.42 to 66.57 Hz, which is consistent with the results in Table 8. It is worth noting that, since the mechanism has a relatively excellent kinematic performance before optimization, the kinematic performance has not been significantly improved after optimization.

Figure 4 Distribution of the LTI in the regular workspace. (a) Before optimization. (b) After optimization.

Figure 5 Distribution of f₁ in the regular workspace. (a) Before optimization. (b) After optimization.

6 Conclusions

The results of multi-objective optimization design significantly affect the comprehensive performance of PMs. Obtaining the mechanism's optimal comprehensive performance with clear physical significance is challenging because of weight factors' subjectivity and multiple performance indices' non-uniform dimensions. Meanwhile, the high computational costs of the GPI undermine PMs' multi-objective optimization design efficiency. In this paper, the multi-objective optimization design method based on the PCA, factor analysis, PSO algorithm, and neural network was proposed to fix these problems, and the multi-objective optimization of the 3RPS PM was taken as a case study to implement the proposed method.

Considering the fundamental frequency, the motion or force transmission performance, and the regular workspace volume of the 3RPS PM to be objective functions, the results of PCA and factor analysis showed that the first principal component reflected the mechanism's dynamic and kinematic performances. The second one reflected the mechanism's regular workspace volume. The computational cost of the proposed approach was reduced by 96.95 % compared with the PSO algorithm. Meanwhile, the mechanism's overall comprehensive performance increased by 118.92 %, and the dynamic performance, kinematic performance, and workspace volume increased by 44.57 %, 6.74 %, and 283.69 %, respectively.

Although neural networks have a strong fitting ability, they have not revealed the inherent mechanism between the objective functions and the design parameters. In our future work, we will study the analytical formula between the objective function and the design parameters.