MOPRP: Multiobjective optimization on preoperative robot placement for robot-assisted dental implant surgery

Abstract

In robot-assisted dental implant surgery, properly arranging the robot system are important for performance. Following this, a multiobjective optimization approach for preoperative robot placement (MOPRP) is proposed. Firstly, after kinematics modeling, the 9-dimensional planning of both the joint’s angles and the base’s position is simplified into 3-dimensional planning of merely the base’s position, building the design variables. The simplification is based on a multiple inverse kinematics (IK) selection strategy, which contains a collision-free principle and score-based selection. Secondly, three multiobjective indices, including the error range index, dexterity index, and joint variation index, are derived to evaluate placement quality. Subsequently, a Pareto front is generated through the multiobjective genetic algorithm, offering a group of competitive solutions. Lastly, by selecting a favorable solution from the Pareto front and checking its 3D preview, the placement planning can be done. Both numerical tests and experimental validations are conducted, demonstrating that the result approximately ranks in the 80%–90%, 99%, and 80%–90% percentile respectively for the three indices. The results show that MOPRP can provide enhanced dexterity, better accuracy, and less joint variation, with a pretty system layout, verifying its effectiveness and feasibility as a useful planning tool.

Introduction

Tooth loss, especially among the elderly, significantly impacts the daily life.^1,2 Dental implantation can restore most of the functions and is the most effective therapy.^3,4 In recent years, robot-assisted dental implant surgery has become a new solution, which has great potential to provide better surgery and relieve the shortage of dentists. ⁵ With these advantages, robot-assisted dental implant surgery is likely to become mainstream in the future.

The field of dental implant robots has attracted many researchers.^6,7 Sun et al. developed a system integrating a drilling robot, a measuring arm, and an imaging system.^8,9 Cheng et al. proposed a system using a UR5 robotic arm and studied stereo vision. ¹⁰ Tao et al. introduced a collaborative robot system combining a parallel mechanism and force-based guiding, ¹¹ with similar work by Yan et al. ¹² While these works have developed robotic systems, an integrated system alone is not sufficient. The system arrangement (e.g. the position of robot and the angles of robot joints before surgery) remains an important issue. Good arrangement can enhance manipulability and accuracy,^13,14 while poor arrangement can degrade performance and even cause singular workspace issues. ¹⁵

For medical robots, some work has been conducted to plan the system. Hayashibe et al. developed simulation software for preoperative setup, which helps doctors preview the motion. ¹⁶ Yu et al. proposed an automatic planning method through reinforcement learning for their celiac surgery robot, building a simplified model of the abdomen, and obtaining the best entrance posture. ¹⁷ A similar mathematical model of the surgical area is made use of in another work, where three characteristic parameters for the incision are obtained and used to decide the incisions. ¹⁸ Liang et al. developed a preoperative optimization method by evaluating the motion performance including workspace, collision safety, and dexterity for another celiac surgical robot, ¹⁹ and then optimized the preoperative setup using particle swarm optimization. Preda et al. developed a planning software to determine the depth and angle of surgical instruments during laparoscopic surgery. ²⁰ Badani et al. optimized the placement of the camera port to reduce the risk of collision. ²¹ Trejos et al. assessed the performance of an endoscopic robot by involving the global conditioning index and then solved the port location. ²² Yang et al. proposed the planning method with the evaluation indices of visibility, operability, and hand-eye coordination, adopting a multiobjective particle swarm optimization method to determine the suitable area for cholecystectomy. ²³ In addition, Qin et al. pay attention to the preoperative planning of the optical tracking system (OTS) by generating the osteotomy paths and finding the optimal poses of the OTS. ²⁴

In summary, for the preoperative setup of robots, the entrance to intervene in the body and the operation performance are usually being focused on. Other issues such as collision, dexterity, and workspace are common items to be considered for evaluation. Optimization methods including PSO or genetic algorithm (GA), along with their multiobjective versions, are often adopted to obtain an acceptable setup. While previous studies have made an array of progress, there is still some room for further research. Current literature seldom addresses the layout planning for dental/oral surgical robot systems, and clinical doctors still need a planning tool for assistance. Also, reducing motion errors and additional disturbance is rarely considered, which are the specified features of oral/dental surgery. The main contributions/novelty of our MOPRP can be summarized as follows:

• (1) While current literature for robotic dental/oral surgery rarely plans the robot system, our MOPRP arranges the system to enhance performance, serving as a useful tool for practitioners.

• (2) For the high accuracy required in dental implantation, MOPRP analyzes the motion error and optimizes joint configuration, reducing the error transmission from joints to the actuator, which is an issue seldom discussed currently. ²⁵

• (3) MOPRP reduces the disturbance to patients during feeding motion, as patients closely watch the robot whose excessive disturbance can cause stress and impact the surgery, especially since the patient’s head is not fixed. This is a characteristic of oral surgery and contrasts with laparoscopic surgery where patients are under general anesthesia.

• (4) The simplification of design variables in MOPRP is useful for various planning problems, and multiobjective optimization provides a solution set for complex tasks, ²⁶ making MOPRP not only fit the specific needs of dental/oral surgery but also its components have potential for similar applications.

The rest of the paper is organized as follows: Section “Task and main line” introduces the task and an overview of MOPRP. Section “Robot kinematics and IK selection” models the kinematics and provides the IK selection strategy, which is the foundation optimization. Section “Multiobjective indices establishment” analyzes three multiobjective indices for performance evaluation. Numerical tests, experiments, and discussions are carried out in Section “Result and discussion.” Conclusion is provided in Section “Conclusion.”

Task and main line

Planning task

To conduct MOPRP, we use our dental implant robot system,^27–30 involving a robotic arm, an implanting actuator, ²⁷ and a drilling actuator. ³⁰ As shown in Figure 1(a), the arm carries the actuators to perform the surgical operation. As shown in Figure 1(b), MOPRP determines the placement of the robotic arm, including the base’s position and the joint angles, including the position from the oral cavity to the robot base ${\mathit{P}}_{\mathrm{ob}}$ and the joint angles ${\theta }_1,{\theta }_2,\cdots ,{\theta }_6$ , making it a 9-dimensional optimization problem.

Figure 1.

System overview and planning task: (a) the dental implant robotic system and (b) the preoperative placement task for the robotic arm.

Main line of MOPRP

The main line consists of three stages: establishing design variables, building indices, and multiobjective solving, as shown in Figure 2. In stage 1, the 9-dimension planning including ${\mathit{P}}_{\mathrm{ob}}$ and ${\theta }_1,{\theta }_2,\cdots ,{\theta }_6$ is simplified into a 3-dimension planning of ${\mathit{P}}_{\mathrm{ob}}$ using a score-based selecting strategy which selects one representative inverse kinematics (IK) result for ${\mathit{P}}_{\mathrm{ob}}$ . In stage 2, three indices are derived, including the dexterity index that evaluates the distortion of the motion from joints to actuator (controllability), the error range index which evaluates the error transmission from joints to actuator, and the joint variation index which evaluates joint angle changes during feeding motion (Figure 1(a)). In stage 3, the multiobjective genetic algorithm (GA) will solve the planning, returning a solution set called the Pareto front, where we can select a favorable solution to finish the planning.

Figure 2.

Main line and data flow of MOPRP: (a) the main line consists of three stages and (b) detailed data flow of MOPRP.

Robot kinematics and IK selection

Kinematics modeling

Forward kinematics

The product of exponentials (POE) method is used for forward kinematics, where frames {Base}, {End}, and {Oral} are shown in Figure 1(b). Required parameters are listed in Table 1, where ${\mathit{\omega }}_i$ is the unit rotation axis of joint $i$ described in {Base} under zero state, and ${\mathit{r}}_i$ is on the rotation axis. Then the transformation ${\mathit{T}}_{\mathrm{be}}$ from {Base} to {End} can be calculated using equation (1).

$$\{\begin{array}{c}{\mathit{S}}_i=\left[\begin{array}{c}{\mathit{\omega }}_i\hfill \\ \multicolumn{1}{c}{-{\mathit{\omega }}_i\times {\mathit{r}}_i}\end{array}\right]\hfill \\ {\mathit{T}}_{\mathrm{be}}={\mathrm{e}}^{{\hat{\mathit{S}}}_1{\theta }_1}{\mathrm{e}}^{{\hat{\mathit{S}}}_2{\theta }_2}{\mathrm{e}}^{{\hat{\mathit{S}}}_3{\theta }_3}{\mathrm{e}}^{{\hat{\mathit{S}}}_4{\theta }_4}{\mathrm{e}}^{{\hat{\mathit{S}}}_5{\theta }_5}{\mathrm{e}}^{{\hat{\mathit{S}}}_6{\theta }_6}{\mathit{T}}_{\mathrm{be}}^{(0)}\hfill \end{array}$$

(1)

where ${\mathit{S}}_i$ is the axis for joint $i$ , ${\mathit{T}}_{\mathrm{be}}^{(0)}$ represents the zero state, $e^{{\hat{\mathit{S}}}_i{\theta }_i}$ is the transformation motivated by joint $i$ , and $\hat{\mathit{S}}=\left(\begin{array}{cc}\hat{\omega }\hfill & \mathit{\omega }\times \mathit{r}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right)$ is the matrix form of ${\mathit{S}}_i$ , with $\hat{\omega }$ is the skew-symmetric matrix form of ${\mathit{\omega }}_i$ .

Table 1.

The POE parameters for kinematics modeling.

Joint	${\mathit{\omega }}_i$	${\mathit{r}}_i$ (mm)
1	${\left(\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right)}^T$	${\left(\begin{array}{ccc}0\hfill & 0\hfill & 0\hfill \end{array}\right)}^T$
2	${\left(\begin{array}{ccc}0\hfill & -1\hfill & 0\hfill \end{array}\right)}^T$	${\left(\begin{array}{ccc}0\hfill & 0\hfill & 162.8\hfill \end{array}\right)}^T$
3	${\left(\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \end{array}\right)}^T$	${\left(\begin{array}{ccc}0\hfill & 0\hfill & 809.8\hfill \end{array}\right)}^T$
4	${\left(\begin{array}{ccc}0\hfill & -1\hfill & 0\hfill \end{array}\right)}^T$	${\left(\begin{array}{ccc}0\hfill & 0\hfill & 1410.3\hfill \end{array}\right)}^T$
5	${\left(\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right)}^T$	${\left(\begin{array}{ccc}0\hfill & -201.3\hfill & 0\hfill \end{array}\right)}^T$
6	${\left(\begin{array}{ccc}0\hfill & -1\hfill & 0\hfill \end{array}\right)}^T$	${\left(\begin{array}{ccc}0\hfill & 0\hfill & 1512.8\hfill \end{array}\right)}^T$

Jacobian matrix

The Jacobian matrix is the core for the error mapping index and dexterity index. The space Jacobian ${\mathit{J}}_s$ is shown in equation (2). The body Jacobian ${\mathit{J}}_b$ , shown in equation (3), which is calculated from ${\mathit{J}}_s$ , can be directly used for the two indices.

$$\{\begin{array}{c}{\mathit{J}}_s=\left(\begin{array}{cccc}{\mathit{S}}_1\hfill & \mathbf{A}{\mathbf{d}}_{{\mathit{T}}_{b1}}{\mathit{S}}_2\hfill & \cdots \hfill & \mathbf{A}{\mathbf{d}}_{{\mathit{T}}_{b5}}{\mathit{S}}_6\hfill \end{array}\right)\hfill \\ {\mathit{T}}_{\mathrm{bi}}={\mathrm{e}}^{{\hat{\mathit{S}}}_1{\theta }_1}{\mathrm{e}}^{{\hat{\mathit{S}}}_2{\theta }_2}\cdots {\mathrm{e}}^{{\hat{\mathit{S}}}_i{\theta }_i}\hfill \end{array}$$

(2)

where $\mathbf{A}{\mathbf{d}}_{{\mathit{T}}_{\mathrm{bi}}}$ is the adjoint transformation for ${\mathit{T}}_{\mathrm{bi}}$ , written as $\mathbf{A}{\mathbf{d}}_{{\mathit{T}}_{\mathrm{bi}}}=\left(\begin{array}{cc}{\mathit{R}}_{\mathrm{bi}}\hfill & 0\hfill \\ \multicolumn{1}{c}{{\hat{\mathit{t}}}_{\mathrm{bi}}{\mathit{R}}_{\mathrm{bi}}}& {\mathit{R}}_{\mathrm{bi}}\hfill \end{array}\right)$ with ${\mathit{R}}_{\mathrm{bi}}$ and ${\mathit{t}}_{\mathrm{bi}}$ being the rotation and translation from {Base} to joint $i$ .

$${\mathit{J}}_b={\left(\mathbf{A}{\mathbf{d}}_{{\mathit{T}}_{\mathrm{be}}}\right)}^{-1}{\mathit{J}}_s$$

(3)

where $\mathbf{A}{\mathbf{d}}_{{\mathit{T}}_{\mathrm{be}}}$ is the adjoint transformation for ${\mathit{T}}_{\mathrm{be}}$ .

Target position and posture

To perform implantation, {End} should coincide with {Oral}, the Y-axis and X-axis of {End} should respectively coincide with the Z-axis and the negative Y-axis of {Oral}. Since Joint 1 is a 360° joint, it is natural to fix the orientation of {Base}: the X-axis and Z-axis of {Base} respectively parallel with the X-axis and Y-axis of {Oral}, displayed in Figure 1(b). With these settings, the joint angles for a given ${\mathit{P}}_{\mathrm{ob}}$ can be obtained in equations (4)–(6).

$${\mathit{T}}_{\mathrm{bo}}=\left(\begin{array}{cc}{\mathit{R}}_{\mathrm{bo}}\hfill & -{\mathit{R}}_{\mathrm{bo}}{\mathit{P}}_{\mathrm{ob}}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right)=\left(\begin{array}{cc}{\mathit{R}}_{\mathrm{be}}{\mathit{R}}_{\mathrm{eo}}\hfill & -{\mathit{R}}_{\mathrm{be}}{\mathit{R}}_{\mathrm{oe}}^{\mathrm{T}}{\mathit{P}}_{\mathrm{ob}}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right)$$

(4)

$${\mathit{T}}_{\mathrm{oe}}=\left(\begin{array}{cccc}0\hfill & 0\hfill & -1\hfill & 0\hfill \\ \multicolumn{1}{c}{-1}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right)$$

(5)

$$\mathit{\theta }=\mathrm{ik}({\mathit{T}}_{\mathrm{bo}}{\mathit{T}}_{\mathrm{oe}})=\mathrm{ik}({\mathit{T}}_{\mathrm{be}}^{\mathrm{oral}})$$

(6)

where ${\mathit{T}}_{\mathrm{bo}}$ is the transformation from {Base} to {Oral}, ${\mathit{R}}_{\mathrm{bo}}$ is the rotation from {Base} to {Oral}, ${\mathit{R}}_{\mathrm{be}}$ is the rotation from {Base} to {End}, ${\mathit{T}}_{\mathrm{oe}}$ is the transformation from {Oral} to {End} whose rotation part is ${\mathit{R}}_{\mathrm{oe}}$ , $\mathit{\theta }$ is the IK result, $\mathrm{ik}$ is the IK function, and ${\mathit{T}}_{\mathrm{be}}^{\mathrm{oral}}$ is the target position/posture from {Base} to {End} when {End} is aligned with {Oral}.

Multiple IK results selection

The IK result in equation (6) usually contains multiple candidates such as eight groups and the objective function cannot be uniquely determined. Thus, a representative IK result need be selected for ${\mathit{P}}_{\mathrm{ob}}$ , making the objective function uniquely determined, which includes two parts: removing invalid results and score-based selection.

Removing invalid results

Firstly, undesired IK results can be filtered using collision-free principle which has two conditions shown in equation (7).

$$\{\begin{array}{c}\Vert {\theta }_2\Vert <{90}^{°}\hfill \\ \multicolumn{1}{c}{{\mathit{z}}_e^{(\mathrm{xy})}·{\mathit{p}}_{\mathrm{be}}^{(\mathrm{xy})}<0}\end{array}$$

(7)

where condition 1 prevents the arm from contacting the cabinet shown in Figure 3(a). In condition 2, ${\mathit{z}}_e^{(\mathrm{xy})}$ is the Z-axis of {End} and projected to the XY plane of {Base}, ${\mathit{p}}_{\mathrm{be}}^{(\mathrm{xy})}$ denotes the vector from {Base} to {End} and projected to the XY plane of {Base}. Condition 2 indicates that the actuators should not be located on the inner side of joint 5, because when rotating ${\theta }_6$ to switch actuators, a collision displayed in Figure 3(b) will occur.

Figure 3.

Collision-free principle for removing invalid IK results: (a) condition for joint 2 and (b) collision between actuators and arm should be avoided.

Score-based selection

After removing invalid IK results, a score-based strategy will select the representative from the remaining candidates. Now, assuming three indices (Section “Multiobjective indices establishment”) and $N$ remaining candidates: For each index, all candidates will be sorted by their index values and assigned a score, with the $k$ -th highest candidate being assigned $N-k$ scores. For instance, if $N=4$ , under one index, the best candidate scores 3, the second best candidate scores 2, and the worst one scores 0. After all indices are sorted, the candidate with the highest total score will be the representative IK result. The flowchart of the selection strategy is shown in Figure 4.

Figure 4.

Flowchart of score-based selection, which selects one representative from multiple IK results.

In the IK selection strategy, solving the inverse kinematics is done by the robot’s SDK function, returning all the groups of IK candidates, which are the input of the score-based selection strategy. That is to say, IK selection is the post-processing of IK solving. Moreover, the IK selection plays the role of simplifying constrained optimization, for which other methods like the penalty function method or Lagrange multiplier method may also be available. However other traditional methods will incorporate the multidimensional and complex kinematics constraint function without reducing the problem dimension or even importing additional dimensions, making the optimization solving too complex. In addition, when two IK candidates simultaneously score the highest, we directly choose one randomly, record its index number for subsequent processing, and the unselected candidate will be ignored. This may lose some useful information (discussed later), but it still provides a sufficiently competitive IK candidate, and fully meets the requirement of dimension reduction.

Multiobjective indices establishment

Dexterity index

The dexterity index $f_{\mathrm{dex}}$ is based on the condition number of ${\mathit{J}}_b$ , written as ${\sigma }_{\max }/{\sigma }_{\min }$ . If the maximum singular value ${\sigma }_{\max }$ is very large, a small joint movement may lead to a large movement of the actuator. If the minimum singular value ${\sigma }_{\min }$ is very small, a small actuator’s movement may require a huge joint’s movements, indicating a singular configuration. So the smaller the condition number, the better the dexterity. However, the different scales of linear and rotational motion make it less intuitive. Since ${\mathit{J}}_b=({\mathit{J}}_{b\omega };{\mathit{J}}_{\mathrm{bv}})$ assembled by rotational block ${\mathit{J}}_{b\omega }$ the linear block ${\mathit{J}}_{\mathrm{bv}}$ , we can calculate the condition number of ${\mathit{J}}_{b\omega }$ and ${\mathit{J}}_{\mathrm{bv}}$ respectively, and return the larger one:

$$\begin{array}{c}\multicolumn{1}{c}{f_{\mathrm{dex}}=max(\frac{{\sigma }_{max}^{(\omega )}}{{\sigma }_{min}^{(\omega )}},\frac{{\sigma }_{max}^{(v)}}{{\sigma }_{min}^{(v)}})}\end{array}$$

(8)

where ${\sigma }_{\max }^{(\omega )}$ and ${\sigma }_{\min }^{(\omega )}$ are the maximum and minimum singular values of ${\mathit{J}}_{b\omega }$ , reflecting the distortion for the rotational motion mapping. ${\sigma }_{\max }^{(v)}$ and ${\sigma }_{\min }^{(v)}$ are the counterparts of ${\mathit{J}}_{b\omega }$ , checking the linear motion. The worse one is used to be the dexterity index.

Error range index

The error range index calculates the range of the actuator’s tracking error under the assumption that the joint angles have a random error of $[-\delta ,\delta ]$ . The body Jacobian ${\mathit{J}}_b$ is used for error mapping:

$$\{\begin{array}{c}\mathrm{\Delta }{\mathit{X}}_{\theta }={\mathit{J}}_{b\omega }\left(\mathit{\theta }\right)\mathrm{\Delta }\mathit{\theta }\hfill \\ \mathrm{\Delta }{\mathit{X}}_x={\mathit{J}}_{\mathrm{bv}}\left(\mathit{\theta }\right)\mathrm{\Delta }\mathit{\theta }\hfill \end{array}$$

(9)

where $\mathrm{\Delta }{\mathit{X}}_{\theta }$ and $\mathrm{\Delta }{\mathit{X}}_x$ are respectively the posture and position error caused by the joint space error $\mathrm{\Delta }\mathit{\theta }$ . Both of them are described in {End}. Next, the actuator’s position error $e_x$ is expressed in equation (10).

$$\begin{array}{c}\multicolumn{1}{c}{e_x=\sqrt{{\left(\mathrm{\Delta }{\mathit{X}}_x\right)}^{\mathrm{T}}\mathrm{\Delta }{\mathit{X}}_x}}\\ =\sqrt{\mathrm{\Delta }{\mathit{\theta }}^T\left[{\mathit{J}}_{\mathrm{bv}}{\left(\mathit{\theta }\right)}^{\mathrm{T}}{\mathit{J}}_{\mathrm{bv}}\left(\mathit{\theta }\right)\right]\mathrm{\Delta }\mathit{\theta }}\hfill \end{array}$$

(10)

Assuming the joint’s error range to be $\mathrm{\Delta }{\theta }_i=[-\delta ,\delta ]$ , the actuator’s error range is the following optimization problem:

$$\begin{array}{c}max{e}_x^2\left(\mathrm{\Delta }\theta \right)\hfill \\ \mathrm{s}.\mathrm{t}.\mathrm{\Delta }{\theta }_i\in [-\delta ,\delta ]i=1,2,\cdots ,6\hfill \end{array}$$

(11)

Since getting the error range is finding the maximum value of a positive definite quadratic form $e_x^2=\mathrm{\Delta }{\mathit{\theta }}^T{\mathit{J}}_{\mathrm{bv}}^T{\mathit{J}}_{\mathrm{bv}}\mathrm{\Delta }\mathit{\theta }$ when $\mathrm{\Delta }\mathit{\theta }$ distributes in a 6-dimension hypercube, we can transform it into the diagonal standard form:

$$\{\begin{array}{c}{e_x}^2=\mathrm{\Delta }{\mathit{\theta }}^T\left[{\mathit{J}}_{\mathrm{bv}}{\left(\mathit{\theta }\right)}^{\mathrm{T}}{\mathit{J}}_{\mathrm{bv}}\left(\mathit{\theta }\right)\right]\mathrm{\Delta }\mathit{\theta }={\mathit{\alpha }}^T{\mathbf{\Lambda }}_v\mathit{\alpha }\hfill \\ \mathit{\alpha }={{\mathit{Q}}_v}^{\mathrm{T}}\mathrm{\Delta }\mathit{\theta }\hfill \\ {\mathbf{\Lambda }}_v=\mathrm{diag}\left(\begin{array}{cccc}{\lambda }_1\hfill & {\lambda }_2\hfill & \cdots \hfill & {\lambda }_6\hfill \end{array}\right)\hfill \\ {\mathit{Q}}_v={\left(\begin{array}{cccc}{\mathit{\gamma }}_1\hfill & {\mathit{\gamma }}_2\hfill & \cdots \hfill & {\mathit{\gamma }}_6\hfill \end{array}\right)}^{\mathrm{T}}\hfill \end{array}$$

(12)

where ${\mathbf{\Lambda }}_v$ is the diagonal matrix with ${\mathit{J}}_{\mathrm{bv}}^T{\mathit{J}}_{\mathrm{bv}}={\mathit{Q}}_v{\mathbf{\Lambda }}_v{\mathit{Q}}_v^T$ whose ${\mathit{Q}}_v$ is a orthogonal matrix. $\mathit{\alpha }$ is transformed from $\mathrm{\Delta }\mathit{\theta }$ by ${\mathit{Q}}_v$ . ${\lambda }_1,\cdots ,{\lambda }_6$ are the eigenvalues of ${\mathit{J}}_{\mathrm{bv}}^T{\mathit{J}}_{\mathrm{bv}}$ , and ${\mathit{\gamma }}_1,\cdots ,{\mathit{\gamma }}_6$ are eigenvectors. By the diagonal standard form, $e_x^2$ can be expressed as:

$$e_x^2={\lambda }_1{\alpha }_1^2+{\lambda }_2{\alpha }_2^2+\cdots +{\lambda }_6{\alpha }_6^2$$

(13)

Because ${\mathit{Q}}_v$ can be regarded as a rotation operation, the distributing region of $\mathit{\alpha }$ is still a 6-dimension hypercube. Then, the maximum value of equation (11) can be obtained:

$$max\left(e_x^2\right)=({\lambda }_1+{\lambda }_2+\cdots +{\lambda }_6){\delta }^2$$

(14)

Since the eigenvalues of a quadratic form remain unchanged within orthogonal transformation, the value ${\lambda }_1+\cdots +{\lambda }_6$ , which is the trace of ${\mathbf{\Lambda }}_v$ , is also the trace of ${\mathit{J}}_{\mathrm{bv}}^T{\mathit{J}}_{\mathrm{bv}}$ . Then, the position error range $e_x^{\max }$ can be obtained in equation (15). Similarly, the error range for the posture error $e_{\theta }^{\max }$ is as follows, which involves the rotational part ${\mathit{J}}_{b\omega }$ of ${\mathit{J}}_b$ , is written in equation (16).

$$e_x^{max}=\delta \sqrt{\mathrm{tr}\left({\mathit{J}}_{\mathrm{bv}}{\left(\mathit{\theta }\right)}^{\mathrm{T}}{\mathit{J}}_{\mathrm{bv}}\left(\mathit{\theta }\right)\right)}$$

(15)

$$e_{\theta }^{max}=\delta \sqrt{\mathrm{tr}\left({\mathit{J}}_{b\omega }{\left(\mathit{\theta }\right)}^{\mathrm{T}}{\mathit{J}}_{b\omega }\left(\mathit{\theta }\right)\right)}$$

(16)

where $\mathrm{tr}$ is the trace of a matrix. For equation (16), $e_{\theta }^{\max }$ is actually a constant: Taking the rotational part of equation (3), the relation between the rotational part ${\mathit{J}}_{s\omega }$ of the space Jacobian ${\mathit{J}}_s$ is:

$${\mathit{J}}_{b\omega }\left(\mathit{\theta }\right)={\mathit{R}}_{\mathrm{sb}}^{-1}{\mathit{J}}_{s\omega }\left(\mathit{\theta }\right)={\mathit{R}}_{\mathrm{sb}}^{\mathrm{T}}{\mathit{J}}_{s\omega }\left(\mathit{\theta }\right)$$

(17)

By substituting equation (17) into equation (16), it can be seen that $e_{\theta }^{\max }$ can also be represented by ${\mathit{J}}_{s\omega }$ :

$$\begin{array}{c}\multicolumn{1}{c}{e_{\theta }^{max}=\delta \sqrt{\mathrm{tr}\left({\left({\mathit{R}}_{\mathrm{sb}}^{\mathrm{T}}{\mathit{J}}_{s\omega }\left(\mathit{\theta }\right)\right)}^{\mathrm{T}}\left({\mathit{R}}_{\mathrm{sb}}^{\mathrm{T}}{\mathit{J}}_{s\omega }\left(\mathit{\theta }\right)\right)\right)}}\\ =\delta \sqrt{\mathrm{tr}\left({\mathit{J}}_{s\omega }^{\mathrm{T}}\left(\mathit{\theta }\right){\mathit{R}}_{\mathrm{sb}}{\mathit{R}}_{\mathrm{sb}}^{\mathrm{T}}{\mathit{J}}_{s\omega }\left(\mathit{\theta }\right)\right)}\hfill \\ =\delta \sqrt{\mathrm{tr}\left({\mathit{J}}_{s\omega }^{\mathrm{T}}\left(\mathit{\theta }\right){\mathit{J}}_{s\omega }\left(\mathit{\theta }\right)\right)}\hfill \end{array}$$

(18)

Then, extract the rotational part of equation (2), the calculation of ${\mathit{J}}_{s\omega }$ is:

$${\mathit{J}}_{s\omega }=\left[\begin{array}{cccc}{\mathit{R}}_0{\mathit{\omega }}_1\hfill & {\mathit{R}}_1{\mathit{\omega }}_2\hfill & \cdots \hfill & {\mathit{R}}_5{\mathit{\omega }}_n\hfill \end{array}\right]$$

(19)

where ${\mathit{R}}_0$ is the rotation from {Base} to {Base}, ${\mathit{R}}_1$ is the rotation from {Base} to Joint 1, ${\mathit{R}}_5$ is the rotation from {Base} to Joint 5. ${\mathit{\omega }}_i$ is the unit axis of joint $i$ described in {Base}. By substituting equation (19) into equation (18), then equation (20) can be obtained:

$$\begin{array}{c}\multicolumn{1}{c}{e_{\theta }^{max}=\delta \sqrt{\mathrm{tr}\left(\left(\begin{array}{ccc}{\left({\mathit{R}}_0{\mathit{\omega }}_1\right)}^{\mathrm{T}}\left({\mathit{R}}_0{\mathit{\omega }}_1\right)\hfill & \cdots \hfill & {\left({\mathit{R}}_0{\mathit{\omega }}_1\right)}^{\mathrm{T}}\left({\mathit{R}}_5{\mathit{\omega }}_6\right)\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\left({\mathbf{R}}_5{\mathit{\omega }}_6\right)}^{\mathrm{T}}\left({\mathit{R}}_0{\mathit{\omega }}_1\right)}& \cdots \hfill & {\left({\mathit{R}}_5{\mathit{\omega }}_6\right)}^{\mathrm{T}}\left({\mathit{R}}_5{\mathit{\omega }}_6\right)\hfill \end{array}\right)\right)}}\\ =\delta \sqrt{\sum _{i=1}^6{{\mathit{\omega }}_i}^{\mathrm{T}}{\mathit{\omega }}_i}\hfill \end{array}$$

(20)

From equation (20), it can be drawn that the error range of posture is a constant. Therefore, to evaluate the error range, the position error range $e_x^{\max }$ is sufficient:

$$f_{\mathrm{err}}=\delta \sqrt{\mathrm{tr}\left({\mathit{J}}_{\mathrm{bv}}{\left(\mathit{\theta }\right)}^{\mathrm{T}}{\mathit{J}}_{\mathrm{bv}}\left(\mathit{\theta }\right)\right)}$$

(21)

Joints variation index

The joints variation index $f_{\mathrm{var}}$ is used to evaluate the amplitude of the joint’s motion for feeding the actuator into the oral cavity, and is defined to be the average angle’s variation of the six joints:

$$f_{\mathrm{var}}=\frac{\mathrm{sum}\left(|\mathrm{ik}\left({\mathit{T}}_{\mathrm{be}}^{\mathrm{oral}}\right)-\mathrm{ik}\left({\mathit{T}}_{\mathrm{be}}^{\mathrm{init}}\right)|\right)}{6}$$

(22)

where $\mathrm{ik}({\mathit{T}}_{\mathrm{be}}^{\mathrm{oral}})$ is the selected IK joint angles for the target ${\mathit{T}}_{\mathrm{be}}^{\mathrm{oral}}$ aligned with {Oral}, ${\mathit{T}}_{\mathrm{be}}^{\mathrm{init}}$ is the position/posture of {End} before feeding to the oral cavity, for which the frame {End} should retreat along the negative Z-axis of {Oral} by 350 mm. $\mathrm{ik}({\mathit{T}}_{\mathrm{be}}^{\mathrm{init}})$ is the IK result for ${\mathit{T}}_{\mathrm{be}}^{\mathrm{init}}$ . $\mathrm{sum}$ is the sum of all the elements of a vector. Note that the two IK calculations in equation (22) are different. $\mathrm{ik}({\mathit{T}}_{\mathrm{be}}^{\mathrm{oral}})$ is selected by the strategy in Figure 4, and $\mathrm{ik}({\mathit{T}}_{\mathrm{be}}^{\mathrm{init}})$ is the closest result to $\mathrm{ik}({\mathit{T}}_{\mathrm{be}}^{\mathrm{oral}})$ among multiple IK results, not selected by the strategy in Figure 4. The initial and aligned actuator’s position is illustrated in Figure 5(a), the data flow for calculating $f_{\mathrm{var}}$ is shown in Figure 5(b).

Figure 5.

Joints variation index: (a) the actuator’s position before and after motion and (b) flowchart for calculating the index.

Moreover, the above calculation only considers the start and end points. For the feeding, the master-slave control based on a virtual fixture is employed, ²⁹ which follows the manual driving and the virtual fixture that corrects the manual commands by hindering undesired movements. Although master-slave driving has randomness, a good placement can reduce the joint’s swings significantly. As patients in oral surgery are usually awake and closely focused on the moving robot, smaller joint variation means less stress to the patients, and less disturbance to the doctor, reducing the risk of misoperation.

Result and discussion

Numerical tests

In this section, MOPRP is programmed through MATLAB, and the multiobjective genetic algorithm (GA) is adopted to obtain a solution set called the Pareto front, where all the results are non-dominated. In GA, the population size is 500, the iteration number is 500, and the joint error is set to be $\delta =0.02\mathrm{deg}$ (equation (21)). We run MOPRP under two examples, as shown in Figure 6. In the first example (see Figure 6(a)), the patient is lying on a bed and inclined by 40°; In the second example, the patient is sitting on a chair (see Figure 6(b)) and assumed to be vertical.

Figure 6.

Two examples for numerical tests: (a) the first example where a patient is lying and (b) the second example where a patient is sitting.

For example 1, the Pareto front is displayed in the design and objective space (Figure 7(a) and (b)). We select ${\mathit{P}}_{\mathrm{ob}}^{*}={(-683.11,-383.57,-196.00)}^T$ as the planning result. To extract the joint angles ${\mathit{\theta }}^{*}$ , the multiple IK results of ${\mathit{P}}_{\mathrm{ob}}^{*}$ is taken into the selection strategy (see Figure 4) to get its representative IK result. The system layout after the feeding motion is displayed in Figure 7(c), and another status before feeding the actuator, where the end tool is retreated along the Z-axis of {Oral} by 350 mm, is shown in Figure 7(d). The joint angles ${\mathit{\theta }}_{\mathrm{init}}^{*}$ should be the nearest IK result toward ${\mathit{\theta }}^{*}$ among multiple IK results. In addition, we randomly generate 50,000 valid solutions and evaluate the ranking of ${\mathit{P}}_{\mathrm{ob}}^{*}$ : $f_{\mathrm{err}}$ ranks 86.12% percentile, $f_{\mathrm{dex}}$ ranks 99.81% percentile, and $f_{\mathrm{var}}$ ranks 81.64% percentile. All the parameters of the solution are listed in Table 2.

Figure 7.

Result returned by MOPRP for example 1: (a) Pareto front in design space, (b) Pareto front in objective space, (c) system layout of selected solution after feeding the actuator, and (d) system layout of selected solution before feeding the actuator.

Table 2.

Detailed parameters of the selected solution for example 1.

Parameter	Value
Base’s position ${\mathit{P}}_{\mathrm{ob}}$ (mm)	${(-683.11,-383.57,-196.00)}^T$
Error range index $f_{\mathrm{err}}$ (mm)	0.40
Dexterity index $f_{\mathrm{dex}}$ (times)	1.94
Joint variation index $f_{\mathrm{var}}$ (deg)	17.46
Joint angles after feeding ${\mathit{\theta }}^{*}$	${(2.82,-18.77,125.62,144.39,87.18,-130.00)}^T$
Joint angles before feeding ${\mathit{\theta }}_{\mathrm{init}}^{*}$	${(25.88,1.00,116.06,115.06,64.12,-130.00)}^T$
Ranking of $f_{\mathrm{err}}$ (percentile)	86.12%
Ranking of $f_{\mathrm{dex}}$ (percentile)	99.81%
Ranking of $f_{\mathrm{var}}$ (percentile)	81.64%

For example 2, the Pareto front is shown in Figure 8(a) and (b). The layout after feeding motion is shown in Figure 8(c), and the layout before feeding is shown in Figure 8(d). The ranking of the selected solution are: $f_{\mathrm{err}}$ ranks 80.65% percentile, $f_{\mathrm{dex}}$ ranks 98.76% percentile, and $f_{\mathrm{var}}$ ranks 86.38% percentile. All the detailed parameters for the result are listed in Table 3.

Figure 8.

Result returned by MOPRP for example 2: (a) Pareto front in design space, (b) Pareto front in objective space, (c) system layout of selected solution after feeding the actuator, and (d) system layout of selected solution before feeding the actuator.

Table 3.

Detailed parameters of the selected solution for example 2.

Parameter	Value
Base’s position ${\mathit{P}}_{\mathrm{ob}}$ (mm)	${(-645.85,-336.04,-513.85)}^T$
Error range index $f_{\mathrm{err}}$ (mm)	0.41
Dexterity index $f_{\mathrm{dex}}$ (times)	2.12
Joint variation index $f_{\mathrm{var}}$ (deg)	17.17
Joint angles after feeding ${\mathit{\theta }}^{*}$	${(7.77,-3.07,125.02,-51.91,-82.23,0.00)}^T$
Joint angles before feeding ${\mathit{\theta }}_{\mathrm{init}}^{*}$	${(43.02,-5.55,122.28,-52.17,-46.98,0.00)}^T$
Ranking of $f_{\mathrm{err}}$ (percentile)	80.65%
Ranking of $f_{\mathrm{dex}}$ (percentile)	98.76%
Ranking of $f_{\mathrm{var}}$ (percentile)	86.38%

Experimental validation

The experiment validation uses the example from Section “Kinematics modeling,” including tests on error mapping, joint variation, and dexterity performance. The experimental setup, shown in Figure 9(a), includes a robot arm, a master handle, an oral model, and an actuator. During the experiments, we tested the selected solution (Table 3) and 20 random solutions, providing a systematic comparison.

Figure 9.

Working site for conducting experimental validation.

Result of error mapping test

To evaluate error mapping, the robot arm is commanded to run a circular path with a 25 mm radius and 10 mm/s velocity, as shown in Figure 10(a). The motion error is evaluated by the deformation from the desired to the executed path, as shown in Figure 10(b) to (e). All the 21 solutions are tested, acquiring the distribution of path error, as depicted in Figure 11. Although samples 5, 8, 10, and 11 perform better, the selected solution ranks well compared to all random solutions.

Figure 10.

Evaluate error mapping by executing a circle path: (a) working site, (b) desired and executed path for the selected solution, (c) desired and executed path for the 5th random solution, (d) counterparts for the 10th random solution, and (e) counterparts for the 15th random solution.

Figure 11.

Average value and distribution of the selected solution and 20 random sampled solutions in path error tests, the label “win” represents the selected solution.

Result of joint variation test

For this test, we use the master-slave control based on virtual fixture. When moving the master handle, the actuator is pushed into the oral cavity. For each of the 21 solutions, 10 sub-groups of feeding are performed to obtain statistical results. The working site before and after feeding the actuator is shown in Figure 12(a) and (b). The average joint variation curve for all 10 sub-groups of the selected solution is shown in Figure 12(c), and the actuator’s path for all 10 sub-groups is shown in Figure 12(d). Note that the joint’s variation is the overall traveled distance, rather than simply subtracting the starting angle from the ending angle, as written in equation (23):

$$\mathrm{\Delta }{\theta }_i={\int }_{t_b}^{t_e}\left|v_i\right|\mathrm{d}t$$

(23)

where $\mathrm{\Delta }{\theta }_i$ is the swing angle, $t_b$ and $t_e$ are the beginning and ending time, $v_i$ is the angular velocity of the joint.

Figure 12.

Evaluate joint variation under master-slave control: (a) working site before feeding the actuator, (b) working site after feeding the actuator, (c) increase process of average joint swing angel of all the 10 sub-groups for the selected solution, and (d) executed paths of all the 10 sub-groups for the selected solution.

The average joint variation and its distribution for all 21 solutions are shown in Figure 13. Except for samples 9, 12, and 15, no samples rank better than the selected solution. For error mapping (Figure 11) where samples 5, 8, 10, and 11 had smaller motion errors, we see that these samples perform much worse in joint variation tests. Conversely, samples with smaller joint swings rank much poorer in path error tests. This indicates that while a random sample might perform slightly better in one test, it rarely excels in another, highlighting MOPRP’s advantage in providing a well-balanced solution.

Figure 13.

Average value and distribution of the selected solution and 20 random sampled solutions in joint variation tests, the label “win” represents the selected solution.

Furthermore, while the specific master-slave strategy can influence the result (joint swing angle) for a single test, it will not impact the statistical results of multiple tests under the same master-slave control strategy. There are various methods for master-slave control, such as absolute mapping, incremental mapping, Euler angle mapping, rotation vector mapping, and some virtual fixture methods. Under different master-slave strategies, the patterns of the robot motion process differ, thus affecting the result of a single experiment. However, because multiple tests under the same master-slave strategy can statistically demonstrate the result distribution, the whole tendency of repeated results will remain the same, therefore the specified master-slave method does not influence the entire joint swing angle test.

Result of dexterity test

To check dexterity, the General Isotropy Index (GII) is introduced, which is a more comprehensive evaluation. We randomly sample 100 points within a 50 × 50 × 50 mm, 50 × 50 × 50 deg space of the oral cavity, and calculate their 6 × 6 body Jacobian (equation (3)). Subsequently, we take the smallest and the largest singular value among all 100 Jacobian matrices and calculated their ratio to be the GII index, as shown in equation (24). The larger this index, the better the dexterity. The performance of all these 21 solutions in GII is plotted in Figure 14, and the selected solution definitely has the best ranking.

$$\mathrm{GII}=\frac{\underset{p_1\in W}{min}{\sigma }_{min}\left(\mathit{J}\left(p_1\right)\right)}{\underset{p_2\in W}{max}{\sigma }_{max}\left(\mathit{J}\left(p_2\right)\right)}$$

(24)

where W is the point set distributed in the oral space, J is the 6 × 6 Jacobian matrix, $p_1$ has the minimum singular value among all points of W, and $p_2$ has the maximum singular value among all points of W.

Figure 14.

Bar plots of GII values for all 21 solutions, the label “win” represents the selected solution.

Furthermore, there are two reasons why GII is not used as the optimizing objective. First, the linear and angular units are not unified, making it difficult to intuitively check the dexterity and in particular during the selection stage. Separating linear and angular motion in the dexterity index is more intuitive for assessing the motion isotropy. Second, sampling many points within the workspace, multiplied by the population of GA, will result in a huge computational burden and slow down the solving. Additionally, dexterity test is based on numerical calculations because it is fully characterized by the Jacobian matrix rather than direct measurement from robot motion.

Discussion

About indices and results

In MOPRP, three indices are introduced, including the dexterity index, error mapping index, and joint variation index. The reasons for using these indices are as follows. The dexterity index is responsible for checking the manipulability and avoiding too large/small joint-to-actuator motion transmission, thus it can improve the maneuverability and prevent malformed manipulating or uncontrollable situations from happening. Also, we separate the linear and angular parts of the Jacobian matrix and use the larger condition number as the index value. This processing can make users check the quality of the results more intuitively. The error mapping index is responsible for curbing the path tracking error during the implantation since the accuracy requirement is high. Under this requirement, a robot configuration that can transmit less error from the joint space to the actuator will be preferable. On the contrary, the joint error caused by servo control, trajectory interpolation, and mechanical gap will be magnified, lowering the quality of hole drilling and implant inserting. The joint variation index helps produce smooth feeding, reducing unnecessary impact or disturbance on patients and doctors. In other words, under this index, a good configuration can produce less swing of the arm, reducing possible stress of the patient, and making the doctor more focused on the operation. Thus, this easily neglected index can be of great practical importance. To sum up, the aspects of manipulability, accuracy, and operating experience are paid attention to, forming a comprehensive multiobjective evaluation.

As the indices $f_{\mathrm{dex}}$ , $f_{\mathrm{err}}$ , and $f_{\mathrm{var}}$ are respectively responsible for checking the manipulability, evaluating the path error, and producing a smooth feeding, they are inherently influenced and need to reach a trade-off. First, the index $f_{\mathrm{dex}}$ will limit both excessively large and small transmission ratios (from joint space to Cartesian space), meaning it emphasizes the balance. Then, the error range index $f_{\mathrm{err}}$ , which reflects the amplification of joint errors, prefers a smaller transmission ratio. At last, the joint variation index $f_{\mathrm{var}}$ prefers a larger transmission ratio, since a larger transmission ratio will require less joint motion, making the feeding less disturbing. Therefore, blindly optimizing $f_{\mathrm{err}}$ will lead to the worst $f_{\mathrm{dex}}$ and $f_{\mathrm{var}}$ , and merely optimizing $f_{\mathrm{var}}$ will deteriorate $f_{\mathrm{err}}$ , hence a balance needs to be reached to provide a good trade-off.

Following this, conflicts between the indices make it difficult to achieve ideal performance simultaneously. For instance, top rankings for $f_{\mathrm{err}}$ and $f_{\mathrm{var}}$ are hard to obtain together. Three clusters of the Pareto front are shown in Figures 7 and 8, each with different focuses. In Figure 8, cluster 1 has a small error range but more joint variation, cluster 2 emphasizes joint variation with worse dexterity, and cluster 3 has the best $f_{\mathrm{dex}}$ and $f_{\mathrm{var}}$ but a poorer error range. These clusters indicates the inherent influences between indices. The experimental results align with the numerical tests: the selected solution achieves about the 85% percentile in error mapping and joint variation experiments compared to 20 random samples, with the best ranking in dexterity among all samples. Although not the top in error mapping and joint variation, the selected solution offers the best balance. Some samples may perform better in one index but they always rank much lower in another index. The results also shows the inherent influences between indices, showing that only multiobjective optimization can yield a competitive result across all objectives.

About method details

In the method design, while the IK selection strategy plays an important role for reducing the problem dimensions, there are three potential limitations. The first case is more than one IK solutions achieve the highest score. In this case, we randomly select one, and other IK solutions will not be considered, losing some useful IK information. The second case involves a non 6-DOF arm. The IK solution for redundant arms is infinite, meaning the introducing of additional constraints or numerical iteration. For an underactuated arm, due to the loss of DOF, the actuator cannot freely satisfy a given oral-to-base position with a fixed posture, and thus additional base-to-oral posture adjustment should be solved. The third case is the potential unfairness of the score assignment principle. Because it is hard to precisely determine the weight of each index, so we treat the three indices equally. But, the equal treatment may cause some potential unfairness in that the actual importance differs between them.

Despite these potential limitations, most collaborative arms usually have 6 DOFs. Moreover, since half of IK solutions are removed by collision conditions, the possibility of tied first-place scores is also low. In addition, although some possible imbalance between indices may cause unfairness, the winner IK candidate is reliable enough and of high efficiency. Thus, these potential limitations in this DOF simplification process is slight and entirely acceptable.

About practical implementation

Numerical tests and experiments show two layout scenes: aligned status after feeding the actuator (Figures 7(c) and 8(c)) and initial status before feeding the actuator (Figures 7(d) and 8(d)). While optimization focuses on the post-feeding status, the initial status is more important for practice as it can be directly arranged by the doctor.

The layout implementation involves an adjustable patient bed/chair, an integrated cart, and a stereo-vision system. First, the patient bed height should be adjusted via a remote controller to meet the selected solution in the Z direction, making more Pareto front solutions accessible. Second, for oral surgery robots, the robot arm, actuator, and the control system are typically integrated into a wheeled cart (see Figure 7), which can be easily moved and locked, making it matches the selected solution in the X and Y direction. Third, a stereo vision system should be involved, which can detect a vision marker worn on the patient, determining the oral cavity’s position. After hand-eye calibration, the relative position between the vision system and the robot base can be known. The doctor can then adjust the bed/chair and move the cart until the measured position matches the selected solution. Thus, for a fully deployed surgical robot system, the selected solution can be easily implemented.

Conclusion

In this study, we propose a preoperative placement planning approach called MOPRP, aimed at improving motion quality by planning the system layout before robot-assisted dental implant surgery. The approach is conducted by establishing design variables, deriving multiobjective indices, and optimization solving. During the establishment of design variables, we designed a multiple inverse kinematics (IK) selection method that simplifies the original 9-dimensional optimization problem into a 3-dimensional optimization by removing invalid IK results and employing a score-based strategy. Additionally, we derived three indices to evaluate the actuator’s dexterity, error range, and joint variation performance when certain joint angles are given and found that only the position error range can be optimized, and the posture error range remains unchanged. By running the multiobjective genetic algorithm, we obtained a series of competitive solutions, the Pareto front, from which a solution is selected and validated by evaluating its ranking and checking its system layout. Based on the selected solution, we conducted a series of experiments on the real robot, confirming the effectiveness and feasibility of the proposed MOPRP method.

In the future, we plan to develop 3D software that packs the MOPRP core solver, result selection, percentile ranking, and 3D previewing, which can facilitate on-site use and guide the layout arrangement for practitioners. Also, we will further test MOPRP by comparing our MOPRP with other similar works in terms of accuracy and efficiency, offering more detailed performance data.