Identifying optimal trajectory parameters in robotic finishing operations using minimum number of physical experiments

Abstract

In this paper, we have presented a method to identify optimal trajectory parameters in process applications like robotic finishing, while minimizing the required number of physical experiments. Our method optimizes the task objective and meets the task performance constraints. The algorithm makes decisions based on the uncertainty in the surrogate models of task performances or black-box constraints. The method intelligently samples the parameter space to select a point for evaluation from the sampled set by determining its probability to be optimum within the set. The iterative method rapidly converges to the optimal point with a small number of experiments. We have proved that our method will converge faster than any other method in the expected sense. We have considered synthetic test problems to benchmark our method against other methods. We have validated our approach through physical experiments of robotic scrubbing and sanding. This algorithm is general enough to be applicable to mathematical constraint satisfaction problems where the objective function is known and the constraint functions are unknown.

Keywords

Robotic surface finishing parameter optimization constraint satisfaction problem optimization under black-box constraints

1. Introduction

Robots are traditionally used in assembly, machine tending, and material handling tasks. There is increasing interest in using robots for process applications. In these applications, the robot is required to make changes to the part being processed. Robotic surface finishing tasks such as cleaning, polishing, sanding, machining, and painting are some examples of process applications. The underlying physics-based model of the process application is known in some cases. Analytical optimization methods or simulation-based approaches can be used in these cases to identify and optimize operation parameters to achieve desired task performance. However, physics-based models are not known a priori for many applications. This is often the case when a new material, part, or tool is under consideration. This problem also arises when there is significant uncertainty in the state of the work-piece or the tool.

In this paper, we have focused on process applications where the underlying physics-based process models are not known a priori. A large number of physical experiments need to be conducted to build a complete physics-based model. This can be done for large volume production. Building a complete model first and then using it for parameter optimization is not practical for low volume production. Moreover, the complete model might not be reused due to the non-repetitive nature of tasks in small volume manufacturing. Our interest is in devising a method that can identify the optimal operation and trajectory parameters with a small number of experiments. This approach will have to focus on identifying optimal operation parameters instead of building an accurate physics-based model.

In process applications, the tasks are performed by robots manipulating a tool, which is in physical contact with a work-piece. The motion parameters (direction, speed, etc.) and the tool interaction parameters (stiffness, deflection, force, etc.) govern the trajectory followed by the robot. Our interest is in identifying the optimal set of parameters for the robot trajectory using the physical experimentation results.

There are three main considerations in processing tasks. The task objective is generally the first consideration. Most tasks require minimizing time and/or cost. We can express the task completion time and costs as a function of trajectory parameters. The task performance constraints are usually the second consideration. For example, a specific tolerance for matching operations, certain surface roughness, or minimum level of cleaning might need to be achieved. Usually, the mapping between the trajectory parameters and task performance variables is very complex. In this paper, we assume that such a model is not known in advance. The third consideration is the constraint on part damage. It comes from the fact that the robot will be applying force on the part using a tool. For example, we want to prevent breaking or scratching the part during the cleaning or polishing process. We can mathematically transform these three considerations into a constrained optimization problem with a known objective function and black-box constraint functions.

In this paper, we have represented the black-box performance constraints using surrogate models based on Gaussian Processes. The uncertainties in the surrogate models are taken into account to identify the optimal parameters. This enables us to calculate the probability of candidate points to be optimal and predict the parameter values that will have the maximum probability to be optimal. The algorithm iteratively proceeds by reducing uncertainties near the predicted optimal point. This allows it to converge quickly to the optimal point without the need for conducting an extensive number of experiments at other portions of the search space.

In our earlier work [1], we developed an algorithm for optimizing task objective considering black-box task performance constraints. We used a probability thresholding based search method that did not consider minimizing number of experiments. In our recent work [2], we generalized the method to minimize the number of experiments to identify optimal trajectory parameters. However, we did not have a way to minimize computation time for higher dimensions. This paper builds on our previous work [2] with the following new contributions:

1.
Theoretical analysis and proof of faster convergence rate over other methods.
2.
Candidate points pruning algorithm to reduce computation time in higher dimensions.
3.
More synthetic test problems to benchmarkagainst other methods.
4.
Illustration of the algorithm’s performance on a fine resolution fixed 2D dense grid.
5.
Validation through physical experiments ofrobotic sanding.

2. Related work

The problems addressed in the robotic finishing tasks that are considered in this paper are similar to many manufacturing processing applications like polishing [3], grinding [4], milling [5], cutting [6], spray coating [7], welding [8, 9], and assembly [10, 11, 12]. In all these processes, there is significant interaction between the tool tip and the work-piece. Usually, a large number of process parameters have varying degrees of influence on task performance. Model-driven optimization methods can perform well in simulation when the physical model of the system is available. However, a common challenge faced in all these problems is a lack of physics-based analytical models to derive functional relationships between process parameters and input-output behavior [6]. Only the input-output behavior of the system can be known through physical experiments. Therefore, most researchers in the field resorted to data-driven approaches for parameter optimization in different process applications. These data-driven approaches can be broadly classified into:(1) Design of experiments (DOE) approaches [13] (e.g., Taguchi method, factorial design, and response surface design methodology), (2) Surrogate model for task performance based approaches, and (3) Meta-heuristic search approaches (e.g., genetic algorithms, simulated annealing, Tabu search).

The Design Of Experiments (DOE or experimental design) is a method where a task is designed to identify the influence of certain variables on the outcome of an event. We can tune them to achieve a desired outcome by understanding the underlying influence of different parameters. This approach relies on conducting experiments by combinatorially selecting the sets of different parameter values. This method is applicable in scenarios where the outcome of the event is measurable from the experimental data. On the other hand, surrogate model based approaches are useful when the outcome of an event is not directly measurable by conducting experiments or it is not feasible to conduct required experiments. In this approach, a surrogate model, that can approximate the outcome of the real event, is used for experimentation and decision making. The surrogate model needs to be much cheaper to evaluate compared to the real event. It also needs to follow the input-output behavior of the actual model as closely as possible. Meta-heuristic approaches are heuristics designed to guide the search process in an optimization problem. These methods can find a sufficiently good solution without having sufficient knowledge about the entire search space. Meta-heuristics algorithms focus on efficiently exploring the search space to find near-optimal solution without heavy computational complexity. In contrast with DOE and surrogate model based approaches, meta-heuristic algorithms are not problem specific.

2.1 Design of Experiments (DOE) based approaches

Design of experiments (DOE) is a statistical method to determine which parameters are important in a process and the values of these parameters that give rise to optimal performance [13]. It involves fractional factorial experiments to identify influential parameters, full factorial experiments to find the optimal parameter set, and evaluation of the distribution of the objective metrics by running a number of experiments using the optimized parameter set.

Many researchers applied DOE approaches to grinding and polishing tasks [14, 3, 15, 16]. The objective in these tasks was either to minimize surface roughness or maximize material removal rate. Different process parameter sets were considered: (1) Abrasive size, contact force, belt linear velocity, and feed rate,(2) Depth of cut, wheel speed, and work speed,(3) Grinding time, rotating velocity, belt material, and workpiece material. Other machining processes where DOE was applied include cryogenic machining of stainless steel [17], machining of KDP crystals [18], and electro-chemical polishing [19, 20, 21].

DOE was applied in different works [11, 22, 23] to robotic assembly of a transmission torque converter. The goal was to optimize assembly cycle time and first-time-through-rate. Process parameters consisted of a starting point, assembly stage, insertion distance, timeout limit, max number of trials, searching force, rotation speed, rotation angle, force amplitude, and force period and so on. These process-related parameters were converted into robot force control parameters and corresponding robot motions to perform a particular assembly task.

Elangovan et al. [24] applied Taguchi’s method, which offers a distribution-free and orthogonal array design [25], to ultrasonic welding. The goal was to maximize weld strength as a function of parameters like welding pressure, weld time, and amplitude of the vibration. DOE was also applied to spot welding [26] and laser sintering of metallic powder [27].

If the grid resolution becomes finer in the parameter space or the number of parameters increases, then DOE-based methods become intractable. This is the limitation of complete DOE-based approach.

2.2 Surrogate model based approaches

Surrogate model based approaches have been used by Cheng and Chen [12] to optimize parameters in force control based robotic assembly. Their online optimization algorithm (GPRBOA) is based on surrogate models using a Gaussian Process (GP) [28] and Bayesian Optimization Algorithm (BOA) [29]. Their approach was presented to be more effective compared to DOE-based approaches. A tight peg-in-hole insertion problem (40 $\mu$ m clearance) was considered as the assembly process. Insertion force, search force, and search speed were the trajectory parameters under consideration. The GPRBOA algorithm was later improved by Wu et al. [30] using Orthogonal Exploration.

Another approach to surrogate models are metamodels [31] that approximate computationally expensive simulation response functions. The surrogate models presented in this paper are similar to local metamodels that are used within an iterative optimization strategy and updated as the optimization progresses. The metamodels are obtained by conducting experiments at different parameter settings and fitting parametric probability models to observed values of system response. However, we use non-parametric probability estimation models in our approach.

2.3 Meta-heuristics based approaches

Genetic Algorithms (GA) have been used in multi-constraint optimization of grinding Ni-based alloys [4] and force-based robotic assembly [32]. Non-traditional optimization algorithms like artificial bee colony(ABC), particle swarm optimization (PSO), and simulated annealing (SA) have been used for optimization of multi-pass milling operations [5]. Teaching-learning based optimization [33, 34] has been used for machining processes like abrasive water jet machining process, grinding, milling, and turning.

2.4 Constrained optimization under black-box constraints

Many process applications require optimization of a task objective with constraints on the task performance. The physics-based models are not available for the task performances. We consider such problems in this paper. Therefore, we frame the problem as an optimization problem under black-box constraints. We have benchmarked our approach presented in this paper against two other established optimization methods. These methods are designed to handle black-box constraints.

The first method is known as NOMAD [35]. This method is based on the Mesh Adaptive Direct Search (MADS) algorithm, which was developed by Abramson, Audet, and Dennis [36, 37]. The objective of NOMAD is to find the best feasible solution by conducting a small number of experiments. The authors used Kriging-based surrogate models for unknown constraint functions in one version of NOMAD.

The second method utilizes a stochastic Radial Basis Function (RBF) based algorithm. This was developed by Regis [38] for large scale optimization problems. His algorithm considers both objective and constraint functions to be black-box. His method was tested on a number of synthetic test problems and was shown to outperform several other algorithms including a pattern search algorithm, a genetic algorithm, a derivative free trust region based algorithm (COBYLA), a sequential quadratic programming algorithm, and NOMAD.

3. Problem formulation

We have adapted the problem formulation that we introduced in our earlier work [2]. We want to find the solution to the constrained optimization problem by conducting a small number of experiments. We follow the following steps to formulate the problem:

1.
Define the desired task performance targets. These

will be expressed as constraints in our approach. For example, we may want surface roughness to meet or go below a certain given threshold.
2.
Identify robot trajectory parameters (operation parameters) of interest. The task performances are unknown functions of these parameters.
3.
Define the limits of the operation parameters. The robot capabilities and safety considerations will dictate the range of the robot trajectory parameters.
4.
Formulate the task objective function in terms of the robot trajectory parameters. For example, we may be interested in performing the task at the minimum force and the highest possible speed. We express the task objective as a single function. Various sub-objectives are combined in that single function with proper weights.

The problem is mathematically formulated in Eq. (3).

The task objective function, $f(x)$ is known. The task performance functions or constraint functions, $g_{i}(x)$ , are black-box, i.e. they are known by their input-output behavior only.

$\displaystyle\min f(x)$ $\displaystyle\text{s.t.}\ x\in\mathbb{R}^{d},\ a_{j}\leqslant x_{j}\leqslant b% _{j}$ $\displaystyle a,b\in\mathbb{R}^{d},j=1,2,\ldots,d$ (1) $\displaystyle g_{i}(x)\geqslant 0,\ i=1,2,\ldots,n$ $\displaystyle f(x)\!:\ \textit{known},\ g_{i}(x)\!:\ \textit{black-box}$

A two-dimensional synthetic example of the problem is illustrated in Fig. 1. We are assuming that the unknown or black-box performance constraint ( $g(x)$ ) is a normalized 2D Schwefel function with non-uniform bias along the axes. This function is represented by the curved surface in Fig. 1a. We want to find a point which satisfies the task performance constraint, $g(x)\geqslant 0.65$ and minimizes the task objective function $f(x)=x_{1}+x_{2}$ . The points on the curved surface which are above the flat plane are the feasible points. We want to locate the point among the feasible points that has minimum value for $f(x)=x_{1}+x_{2}$ . Figure 1b shows the feasible points on 2D plane. The value of the objective function at these points are represented by colors. Darker blue represents lower value of the objective function. The optimal point is indicated by the red marker.

Figure 1.
A 2D example of parameter optimization problem under black-box constraints. Here the objective is to minimize $f(x)=x_{1}+x_{2}$ . The black-box constraint function, $g(x)$ , is a normalized 2D schwefel function with different bias along each axis. The constraint is $g(x)\geqslant 0.65$ . (a.) The 2D curved surface represents the constraint function $g(x)$ . The flat plane represents the boundary $g(x)=0.65$ . Points on the surface of Schwefel function above the plane are feasible points. The optimal point is pointed with an arrow. (b.) Feasible points on a 2D plane. The values of the objective function at these points are represented with color. Darker blue represents lower values of the objective function. The red marker indicates the optimal point.

4. Approach

Our method can be summarized in the following steps:

10. 1.
Select points in the task parameter space for initial exploratory experiments and conduct experiments at these points. These points can be selected by uniformly sampling the parameter space with a coarse resolution.
2.
Fit surrogate models for the task performance constraints using the experimental data.
3.
Identify candidate points from which the next point for experiment will be selected.
4.
Evaluate the objective function at the candidate points.
5.
Estimate the probability of meeting the task performance constraints for the candidate points using the surrogate models.
6.
For each of the candidate points, calculate its probability to be the optimal point among all the candidate points.
7.
Select the point from the candidate points with the maximum probability to be optimum.
8.
Conduct an experiment at the selected point. Generate new candidate points in the neighborhood of the picked point and remove this point from the list of candidate points. Use a finer sampling resolution for the new candidate points than the resolution that was associated with the experimented point.
9.
Update the surrogate models with the performance data from the new experiment.
10.
Go to step 4.

4.1 Initial exploratory experiments

Our method identifies optimal parameters to meet the desired task performance constraints and optimizes the given task objective function. We use the uncertainty in the task performance surrogate models to select points in the parameter space to conduct experiment. The initial points for exploratory experiments are uniformly sampled from the parameter space at a coarse resolution.

Let $\mathbb{X}_{0}$ and $\mathbb{X}$ be the set of the initial exploration points and set of tested points respectively. We add a point to $\mathbb{X}$ after conducting experiment on it. $\mathbb{X}$ is updated as, $\mathbb{X}\leftarrow\{\emptyset\}\cup\mathbb{X}_{0}$ after completing the initial exploratory experiments. We keep track of the associated sampling resolution of each point. We select the initial sampling resolution to achieve full factorial design [39] of the exploratory experiments. Let, $d$ be the dimension of the parameter space. The coarse sampling resolution, $r_{\textit{initial}}$ , is selected to have $3^{d}$ or $2^{d}+1$ full factorial design for the initial exploratory experiments. The values we used for $r_{\textit{initial}}$ in our experiments are summarized in Table 1. In $2^{d}$ and $3^{d}$ full factorial designs the combinations of [‘low’, ‘high’] and [‘low’, ‘medium’, ‘high’] values along each axis are considered respectively. We are calling $2^{d}+1$ full factorial design as the design where combinations of [‘low’, ‘high’] values along each axis are considered and in addition the mid point of the space is considered. For example, if $d=2$ then the points in $2^{d}+1$ full factorial design will be the followings [‘low’, ‘low’], [‘low’, ‘high’], [‘high’, ‘low’], [‘high’, ‘high’], [‘medium’, ‘medium’]. $2^{d}$ full factorial design only covers the boundary points of the space. We combined the mid point of the space with $2^{d}$ full factorial design to train the surrogate models with a better picture of the landscape.

For parameter space with dimensions 6, 7, 8, 9, and 10, the $2^{d}$ full factorial design will lead to 64, 128, 256, 512, and 1024 initial exploratory points respectively. We want to find the optimal solution by conducting a small number of experiments. Conducting exploratory experiments at every sampled point becomes infeasible for higher dimensional parameter spaces. In such cases we pursue $2^{d-m}+1$ fractional factorial design. We consider randomized fractional factorial design, i.e. we select the $2^{d-m}$ points as a random subset of $2^{d}$ points from the full factorial design. We also include the mid point of the entire space. In the experiments presented in this paper, we decided not to conduct more than 33 initial exploratory experiments. Therefore, for $d\geqslant 5$ , we selected $m$ such that $d-m=5$ .

Table 1
Parameters values used in our design of exploratory experiments

Dimension of	$n_{c}$	$r_{\textit{initial}}$	Design of exploratory
parameter space, $d$			experiments
$d<5$	3	0.5	$3^{d}+1$ full factorial
$5\leqslant d<9$	2	1.0	$2^{d-m}+1$ fractional
			factorial ( $d-m=5$ )
$d\geqslant 9$	1	1.0	$2^{d-m}+1$ fractional
			factorial ( $d-m=5$ )

4.2 Candidate points for refining experiments after the exploratory experiments

We refer to the experiments after the exploratory experiments as refining experiments. This is because we refine our surrogate models at each experiment and actively select the optimal point in expected sense for the next experiment. At each round of refining experiments, we select one point from the set of candidate points $\mathbb{C}_{x}$ . The candidate points are uniformly sampled from the neighborhood of points in $\mathbb{X}$ (set of already tested or experimented points). Along each axis around the tested point, we consider $n_{c}$ neighbor points in both positive and negative direction of the axis at a finer resolution. Table 1 lists values of $n_{c}$ and $r_{\textit{initial}}$ we used in our experiments. Suppose we conducted an experiment at a point, $x$ , which was originally sampled at a resolution of $r_{x}=0.5$ . Then we will sample points in the neighborhood of this tested point at a resolution of $r_{x}/2.0=0.25$ and add them to the set of candidate points, $\mathbb{C}_{x}$ . If we conduct an experiment at one of these new candidate points, then we will use a sampling resolution of 0.125 to generate candidate points around the new experimented point. Figure 2 illustrates examples of tested points and candidate points in 2D parameter space. The red points are the tested points and the black points are the candidate points. Figure 2a shows tested and candidate points after the initial explorations. Suppose, the candidate point marked by the blue ring was selected for the first refining experiment. Candidate points will be now generated around this point once an experiment is conducted. Figure 2b shows updated candidate points and tested points after the first refining experiment.

Figure 2.

Tested and candidate points in a 2D parameter space. Here $n_{c}=3$ . Left: Tested points and candidate points after exploratory experiments. Suppose the candidate point marked by the ring was selected for the first refining experiment. Right: Tested points and candidate points after the first refining experiment.

Our algorithm determines the point to conduct a refining experiment based on its probability to be optimum among the candidates. The probability computation follows through Eqs (4.3)–(4.3). It relies on the uncertainty of outcome at the candidate points. Therefore by design our algorithm avoids regions of the space where the uncertainty of the outcome is higher. Now, if the uncertainty around the global optimum is higher initially, the algorithm will take a larger number of iterations to converge. However, having uniformly sampled initial exploration points and generating candidate points around every tested point helps the algorithm to come out of a local optimum and go towards the global optimum as the surrogate models are refined in every iteration. At every iteration we identify the point $x_{\textit{best}}\in\mathbb{X}$ which minimizes the objective function and meets all the constraints. When $x_{\textit{best}}$ is updated at an iteration, the algorithm generates candidate points at finer resolution around the new $x_{\textit{best}}$ . This helps to find a better solution locally.

4.3 Selecting point for refining experiments

We build surrogate models for the black-box performance constraint functions from the experimental data. We are using Gaussian Processes (GP) for the surrogate models. A distinct GP is trained for each performance constraint. Using GPR (Gaussian Process Regression) [28], we estimate the means and standard deviations in performances for the candidate points. From the mean and standard deviation we can calculate the probability of meeting the performance constraints, $\forall x\in\mathbb{C}_{x}$ , using the following equations.

$\displaystyle P(g_{i}(x)\geqslant 0)=\frac{1}{2}\left[1-\textit{erf}\left(% \frac{-\mu_{i}(x)}{\sigma_{i}(x)\sqrt{2}}\right)\right],$ $\displaystyle\quad\text{ where }\textit{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{% x}e^{-t^{2}}\,dt$ (2)

$\mu_{i}(x)$ is the mean of predicted $g_{i}(x)$ and $\sigma_{i}(x)$ is the corresponding standard deviation.

We compute the probability of meeting all the constraints $\forall x\in\mathbb{C}_{x}$ . For example, if we have $n$ constraints, then the probability of $x_{i}$ meeting all the constraints,

$\displaystyle P_{i}^{c}=P(g_{1}(x_{i})\geqslant 0)\times P(g_{2}(x_{i})% \geqslant 0)\times$ $\displaystyle\quad∼{}\ldots\times P(g_{n-1}(x_{i})\geqslant 0)\times P(g_{n}(x% _{i})\geqslant 0)\!$ (3)

Let, $\mathbb{C}_{x}=\{x_{1},x_{2},\ldots,x_{m}\}$ , where $m$ is the number of candidate points. We keep the points in $\mathbb{C}_{x}$ sorted (low to high) in terms of the value of the objective function, i.e. $f(x_{1})\leqslant f(x_{2})\leqslant\ldots\leqslant f(x_{m})$ . If any subset of $\mathbb{C}_{x}$ happens to have the same value of objective function, then we sort them (high to low) based on their probability of meeting all the constraints. For example, if $f(x_{i})\leqslant f(x_{j})=f(x_{k})=f(x_{l})\leqslant f(x_{m})$ and $P_{k}^{c}\geqslant P_{j}^{c}\geqslant P_{l}^{c}$ , then the points get sorted as $\{x_{i},x_{k},x_{j},x_{l},x_{m}\}$ . Sorting in this way is important because there might be regions in the parameter space where the objective function is a level set, but the constraint functions may vary in the same regions.

We determine the probability of every element of $\mathbb{C}_{x}$ to be optimum among elements of $\mathbb{C}_{x}$ . Since the points in $\mathbb{C}_{x}$ are sorted in terms of the value of objective function, the probability of $x_{1}$ to be optimum is same as its probability to meet all the constraints. For $x_{2}$ , the probability to be optimum will be product of its probability to meet the constraints and probability that $x_{1}$ is not optimum.

Let, $P_{i}$ be the probability of $x_{i}\in\mathbb{C}_{x}$ to be optimum. Then,

$\displaystyle P_{1}=P_{1}^{c}$ $\displaystyle P_{2}=P_{2}^{c}\times(1-P_{1})$ $\displaystyle P_{3}=P_{3}^{c}\times(1-P_{1})\times(1-P_{2})$ $\displaystyle\ldots$ (4) $\displaystyle P_{i}=P_{i}^{c}\times\prod_{j=1}^{i-1}(1-P_{j})$ $\displaystyle\ldots$ $\displaystyle P_{m}=P_{m}^{c}\times\prod_{j=1}^{m-1}(1-P_{j})$

At each iteration we pick the point that has the maximum probability to be the optimal point and conduct one refining experiment. We repeat the process until we converge to a point that meets all performance constraints and optimizes the objective functions. After each refining experiment we update the surrogate models and set of candidate points. We remove a point from the set of candidate points after conducting a refining experiment on it.

We have described our constrained optimization method in Algorithm 1. At each iteration, Algorithm 1 makes call to Algorithms 2 and 3. Algorithm 2 performs uncertainty-based decision making to select the point for the refining experiment from the set of candidate points. We keep on accumulating more and more candidate points in each iteration. The computation time increases with the number of candidate points. Algorithm 3 prunes candidate points based on probabilistic decision making to reduce computation time.

Algorithm 1 takes the following inputs: (1) Dimension of the parameter space, $d$ , (2) Lower and upper bounds on the parameter values, $a$ and $b$ , where $a,b\in\mathbb{R}^{d}$ s.t. $a_{i}\leqslant x_{i}\leqslant b_{i}$ , $\forall x\in S$ , $i=1,2,\ldots,d$ . $S\subset\mathbb{R}^{d}$ is the parameter space, (3) Sampling resolution for initial exploration points, $r_{\textit{initial}}$ . Algorithm 2 takes the following inputs: (1) Set of candidate points, $\mathbb{C}_{x}$ , (2) Array of probabilities for each candidate point to meet each constraint, $G_{C_{x}}$ , and (3) Array of objective function values for each candidate point, $F_{C_{x}}$ . Algorithm 3 takes the following inputs: (1) Set of candidate points, $\mathbb{C}_{x}$ and (2) Array of probabilities for each candidate point to meet all the constraints, $P_{c}$ .

4.4 Computational complexity

Let, $M$ be the number of candidate points, and $N$ be the number of tested points. The computational complexity of the probability computation for the candidate points is $O(M)$ . We have used Gaussian Process Regression (GPR) for uncertainty estimation of the outcome at each iteration. GPR has $O(N^{3})$ time complexity and $O(N^{2})$ memory complexity, where $N$ is the size of the training sample. Therefore the computational complexity of our method at each iteration is $O(n\times N^{3})$ , where $n$ is the number of constraints.

4.5 Optimality in continuous space vs discrete space

We have designed our method to find the optimal solution up to the resolution of the finest grid in the discrete search space. This solution may not coincide with the true global optimum in the continuous space. The user can set the resolution of the finest grid in the discrete search space. Therefore, our method can provide the optimum solution up to the grid resolution required for specific application.

[] [1] Initialize $\mathbb{X}=\{\emptyset\}$ $\mathbb{X}_{0}:=$ {Set of uniformly sampled points in the parameter space}. The sampling resolution is $r_{\textit{initial}}$ . $\mathbb{X}\leftarrow\mathbb{X}\cup\mathbb{X}_{0}$ . Also, $\forall x\in\mathbb{X}$ , keep record of associated sampling resolution. $\mathbb{C}_{x}:=\{$ Set of uniformly sampled points in the neighborhood of each $x\in\mathbb{X}\}$ . The sampling resolution is $r_{\textit{initial}}/2$ . Evaluate $f(x)$ and $g_{i}(x)$ , by conducting experiments $\forall x\in\mathbb{X}$ . $i=1,2,\ldots,n$ . Find $x_{\textit{best}}\in\mathbb{X}$ , s.t. $g_{i}(x_{\textit{best}})\geqslant 0$ , $i=1,2,\ldots,n$ and $f(x_{\textit{best}})=minf(x),\forall x\in\mathbb{X}$ . If such $x_{\textit{best}}$ is found, then generate more candidate points at finer resolution around $x_{\textit{best}}$ . Use $r_{\textit{initial}}/2.0$ as sampling resolution. Train $n$ GP models, one for each of the constraint functions, using values of $g_{i}(x),\forall x\in\mathbb{X}$ . $i=1,2,\ldots,n$ ; where $n$ is the number of constraints. Evaluate $f(x),\forall x\in\mathbb{C}_{x}$ . Let $F_{C_{x}}$ be the array that stores these evaluated values. Each element of $F_{C_{x}}$ is associated with a specific element of $\mathbb{C}_{x}$ . Using GPR, estimate the mean ( $\mu_{i}(x)$ ) and standard deviation ( $\sigma_{i}(x)$ ) of $g_{i}(x)$ , $\forall x\in\mathbb{C}_{x}$ . $i=1,2,\ldots,n$ . Determine probability of meeting each constraint using Eq. (2), $\forall x\in\mathbb{C}_{x}$ . Let $G_{i}$ be the array that stores values of $P(g_{i}(x)\geqslant 0),\forall x\in\mathbb{C}_{x}$ , $i=1,2,\ldots,n$ . Each element of $G_{i}$ is associated with a specific element of $\mathbb{C}_{x}$ . Let $G_{C_{x}}$ be the array that holds $[G_{1},G_{2},\ldots,G_{n}]$ . $x_{\textit{test}},P_{c}\leftarrow$ PickTestPoint( $\mathbb{C}_{x},G_{C_{x}},F_{C_{x}}$ ). Generate new candidate points at the neighborhood of $x_{\textit{test}}$ . The sampling resolution will be $r_{x_{\textit{test}}}/2$ , where $r_{x_{\textit{test}}}$ is the resolution associated with $x_{\textit{test}}$ . Add these new candidate points to $\mathbb{C}_{x}$ and remove $x_{\textit{test}}$ from $\mathbb{C}_{x}$ . $\mathbb{X}\leftarrow\mathbb{X}\cup\{x_{\textit{test}}\}$ . $\mathbb{C}_{x}\leftarrow$ PruneCandidates( $\mathbb{C}_{x},P_{c}$ ). Evaluate $f(x_{\textit{test}})$ and $g_{i}(x_{\textit{test}})$ , $i=1,2,\ldots,n$ , by conducting experiment at $x_{\textit{test}}$ . If $g_{i}(x_{\textit{test}})\geqslant 0$ , $i=1,2,\ldots,n$ and $f(x_{\textit{test}})\leqslant f(x_{\textit{best}})$ , then $x_{\textit{best}}\leftarrow x_{\textit{test}}$ . Generate more candidate points at finer resolution around $x_{\textit{best}}$ . Use $r_{x_{\textit{test}}}/4$ as sampling resolution. Terminate if $x_{\textit{best}}$ has converged, else go to step 8. ParameterOptimization $(d,a,b,r_{\textit{initial}})$

[] [1] $\forall x\in\mathbb{C}_{x}$ , Compute the probability of meeting all the constraints. For example, the probability of $x_{i}$ meeting all the constraints, $P_{i}^{c}=G_{1}(x_{i})\times G_{2}(x_{i})\times\ldots G_{n}(x_{i})$ . We described in Algorithm 1 that $G_{C_{x}}=[G_{1},G_{2},\ldots,G_{n}]$ . Let, $P_{c}=[P^{c}_{1},P^{c}_{2},\ldots,P^{c}_{i},\ldots,P^{c}_{m}]$ be the array of probability of candidate points meeting all the constraints. Sort $F_{C_{x}}$ (low to high). Rearrange the elements of $\mathbb{C}_{x}$ to match their index with associated element in $F_{C_{x}}$ . Let, the rearranged set be $\mathbb{C}_{x}=\{x_{1},x_{2},\ldots,x_{m}\}$ , where $m=|\mathbb{C}_{x}|$ . If a subset of $\mathbb{C}_{x}$ has the same value for the objective function, then sort (high to low) the points in that subset in terms of their probability to meet all the constraints ( $P_{i}^{c}$ ). Calculate the probability of each element of $\mathbb{C}_{x}$ to be optimum among all the elements of $\mathbb{C}_{x}$ using Eq. (4). Let, $P_{i}$ be the probability of $x_{i}$ to be optimum. Let, $P_{\max}=\max\{P_{1},P_{2},\ldots,P_{m}\}$ , where $m=|\mathbb{C}_{x}|$ . Let, $x^{*}\in\mathbb{C}_{x}$ correspond to $P_{\max}$ . Return $x^{*}$ , $P_{c}$ . PickTestPoint $(\mathbb{C}_{x},G_{C_{x}},F_{C_{x}})$

[] [1] Initialize $\epsilon$ (Threshold on probability of meeting all the constraints). Initialize $\mu$ (Threshold on fraction of candidate points to prune). Initialize constants $k_{1},k_{2}$ . Find $T_{x}$ s.t. $\forall x\in T_{x}$ , $P_{x}<\epsilon$ , where $P_{x}$ is the probability of meeting all the constraints for $x$ . For $x=x_{i},P_{x}=P_{c}[i]$ . WHILE ( $|T_{x}|\ /\ |C_{x}|>\mu$ ): IF ( $\epsilon<k_{1}$ ): $\mu\leftarrow\mu/2.0$ ELSE: $\epsilon\leftarrow\epsilon/2.0$ Find $T_{x}$ s.t. $\forall x\in T_{x}$ , $P_{x}<\epsilon$ IF ( $\mu<k_{2}$ ): BREAK $\forall x\in T_{x}$ , remove $x$ from $\mathbb{C}_{x}$ Return $\mathbb{C}_{x}$ PruneCandidates( $\mathbb{C}_{x},P_{c}$ )

5. Theoretical analysis of convergence

We show that the strategy used in our algorithm to solve problem Eq. (3) requires the minimum number of experiments in the expected sense.

.

Let $\mathbb{C}_{x}=\{x_{1},x_{2},\ldots,x_{m}\}$ be the set of $m$ available candidate test points at an instance. Let $P_{i}$ be the probability of $x_{i}\in\mathbb{C}_{x}$ to be optimum, where $\sum_{i=1}^{m}P_{i}=1$ . Let $x^{*}\in\mathbb{C}_{x}$ be the candidate point with the highest probability of being the optimum. That is, $x^{*}=x_{j}$ with $P_{j}\geqslant P_{i}∼{}\forall i\neq j$ , where $i,j\in\{1,2,\ldots,m\}$ . Now, to solve problem Eq. (3) at any instance, if $x^{*}$ is chosen for the experiment (to evaluate the black-box function), then the expected number of trials required for convergence to solution will be less than or equal to any scheme where $x(\neq x^{*})\in\mathbb{C}_{x}$ is chosen for conducting the experiment.

Proof Let us sort the points in the set $\mathbb{C}_{x}$ from high to low w.r.t. their probability of being optimum:

$\displaystyle\mbox{If}∼{}P_{1}\geqslant P_{2}\geqslant\ldots\geqslant P_{m}$ (5)

then the sorted sequence of points is

$\displaystyle\{x_{1},x_{2},\ldots,x_{m}\}.$ (6)

Let $T$ be the number of trials required till convergence for the sequence in Eq. (6). If we use the selection scheme described in Theorem 1, then the expected number of trials $E(T)$ to converge to solution at an instance is

$\displaystyle E(T)=\sum_{i=1}^{m}i\times P_{i}$ (7)

Let us consider the sequence $\{x_{2},x_{1},x_{3},\ldots,x_{m}\}$ . Let $T^{\prime}$ be the number of trials required till convergence for this sequence. We can come up with expected number of trials for this sequence $E(T^{\prime})$ as following,

$\displaystyle E(T^{\prime})=1\times P_{2}+2\times P_{1}+\sum_{i=3}^{m}i\times P% _{i}=(2{-}1)\times P_{2}{+}(1{+}1)\times P_{1}{+}\sum_{i=3}^{m}i\times P_{i}=% \sum_{i=1}^{m}i\times P_{i}+P_{1}-P_{2}=E(T)+P_{1}-P_{2}\geqslant E(T),∼{}% \because P_{1}\geqslant P_{2}∼{}\mbox{from}∼{}\text{(5)}$ (8)

Let us consider another sequence where any two elements of the sequence in Eq. (6) is swapped. The expected number of trials $E(T^{\prime})$ is given by

$\displaystyle E(T^{\prime})=1\times P_{1}+2\times P_{2}+\ldots+q\times P_{r}+% \ldots+r\times P_{q}+\ldots+m\times P_{m}$ (9)

Now, $r>q$ and $P_{q}\geqslant P_{r}$ (from Eq. (5)). Therefore,

$\displaystyle q\times P_{r}+r\times P_{q}$ $\displaystyle\mspace{-6.0mu }=(r-(r-q))\times P_{r}+((r-q)+q)\times P_{q}$ $\displaystyle\mspace{-6.0mu }=r\times P_{r}+q\times P_{q}+(r-q)\times(P_{q}-P_% {r})$ (10)

Using Eq. (5) in Eq. (9),

$\displaystyle E(T^{\prime})=E(T)+(r-q)\times(P_{q}-P_{r})\geqslant E(T)$ (11)

Therefore, from Eqs (8) and (11), for any sequence $\{x_{1}^{\prime},x_{2}^{\prime},\ldots x_{m}^{\prime}\}$ $\neq$ $\{x_{1},x_{2},\ldots,x_{m}\}$ we have $E(T^{\prime})$ $\geqslant E(T)$ . Hence, Theorem 1 is proved.

In our implementation we have used the probability estimate described in Eq. (4.3). If a better design for probability estimate is available then it can be incorporated in our algorithm.

Now, the probability map over the set of candidate points will be updated after each black-box evaluation. We have proven Theorem 1 for each instance or trial. Therefore, in the expected sense, we will always converge faster compared to any other method if we conduct an experiment at the candidate point with highest probability to be the optimum at each iteration.

Figure 3.

Comparison of our method with [35, 38] on synthetic test problems (TP1–TP8). The plots show the mean objective function value at best feasible point at every iteration after the exploratory experiments. Our method, Regis’ method, and NOMAD are denoted by green, red, and blue lines respectively.

Figure 4.

Comparison of our method with [35, 38] on synthetic test problems (TP9–TP16). The plots show the mean objective function value at best feasible point at every iteration after the exploratory experiments. Our method, Regis’ method, and NOMAD are denoted by green, red, and blue lines respectively.

Table 2

Summary of synthetic test problems

Test problem	Dimension of	Inequality	Region defined by parameter space boundary
	parameter space	constraints
TP1	2	2	$[0,500]\times[0,500]$
TP2	2	2	$[13,100]\times[0,100]$
TP3	2	1	$[-1,1]^{2}$
TP4	3	1	$[0,10]^{3}$
TP5	4	3	$[0,1]^{2}\times[0,50]\times[0,240]$
TP6	4	7	$[0.125,10]\times[0.1,10]^{3}$
TP7	4	5	$[0,1200]^{2}\times[-0.55,0.55]^{2}$
TP8	5	3	$[-2.3,2.3]^{2}\times[-3.2,3.2]^{3}$
TP9	5	6	$[78,102]\times[33,45]\times[27,45]^{3}]$
TP10	6	6	$[0,5]\times[0,4]\times[1,5]\times[0,6]\times[1,5]\times[0,10]$
TP11	7	4	$[-10,10]^{7}$
TP12	7	11	$[2.6,3.6]\times[0.7,0.8]\times[17,28]\times[7.3,8.3]^{2}\times[2.9,3.9]\times[% 5,5.5]$
TP13	8	6	$[10^{2},10^{4}]\times[10^{3},10^{4}]^{2}\times[10,10^{3}]^{5}$
TP14	8	6	$[10^{2},10^{4}]\times[10^{3},10^{4}]^{2}\times[10,10^{3}]^{5}$
TP15	10	2	$[0,10]^{10}$
TP16	10	8	$[-10,10]^{10}$

6. Benchmarking on synthetic problems

Physical experiments can be expensive. We conducted experiments on synthetic test problems before physical experiments to validate our method. We benchmarked our method against two other established methods for solving optimization problems under black-box constraints.

Figure 5.

Results on dense uniform grid for Schwefel function. The color map illustrates value of objective function at every grid point. The Tested Points, Next Test Point, Best feasible Point among Tested Points, and True Optimum are plotted on top of the color map, at different iterations of refining experiments.

6.1 Synthetic test problems

We have selected some well established synthetic test problems for benchmarking from [40, 41, 42, 43, 44, 45, 46]. We replaced the equality constraints with inequality constraints. We considered the objective function to be known and the constraint functions to be black-box in all the problems. The description of our selected problems are in the Appendix. The dimensions of the parameter space and number of constraints for the synthetic test problems are summarized in Table 2.

6.2 Optimization methods for benchmarking

We have chosen methods for benchmarking that try to solve the problems of optimization under black-box constraints with a reduced number of experiments. We considered the NOMAD algorithm [35] and the Stochastic RBF (Radial Basis Function) based algorithm developed by Regis [38] for benchmarking our method. The NOMAD algorithm considers mesh adaptive direct search for optimization. This algorithm does not require to start from a feasible point.

RBF was used as surrogate models for the black-box functions in Regis’ algorithm. His algorithm requires at least one feasible point in the set of initial experimentation points. Regis benchmarked his method against six optimization methods, including NOMAD [35]. He considered the objective and constraint functions to be black-box and tested his method on the synthetic problems. His method converged to solution with fewer iterations compared to the other methods. However, none of these methods involve decision making based on uncertainty in the surrogate model.

We modified Regis’ method to solve for known objective functions to compare with our algorithm. For fair comparison, we considered the same initial exploration points for all the three methods.

6.3 Results on synthetic test problems

Comparison of our method against Regis’ algorithm and NOMAD algorithm is illustrated in Figs 3 and 4. The figures illustrate the results on synthetic test problems. From the plots we can see the convergence of the mean of objective function values at the best feasible point in each iteration. The plots only show the points after the initial exploratory experiments.

These problems were of dimension 2 to 10. There were $2^{5}+1=33$ initial exploration points for TP7–TP16. For TP1–TP3 there were $3^{2}=9$ , for TP4 $2^{3}+1=9$ and for TP5–TP6 $2^{4}+1=17$ initial exploration points. The method of selecting points for exploratory experiments were described in Section 4. We can see from the plots that our method took a lower number of experiments to converge to solution compared to the other methods on the synthetic test problems.

At every iteration, Regis’ method looks for points for future experiments in a local neighborhood of the current best feasible solution. Also, neither NOMAD nor Regis’ method consider uncertainty in the surrogate models while selecting points for experiment. Our algorithm stores candidate points in the neighborhood of all the tested or experimented points. This helps our algorithm to select points from a set of candidates sampled over the entire parameter space. This also helps to come out of a local optimum as the surrogate models keep refining themselves. Also, we consider the uncertainty in the models and select the points for experiments based on their expected probability to be optimum among all the candidate points. All these factors enable our algorithm to convergence faster compared to the other two methods.

Figure 6.

Results on dense uniform grid for Schwefel function. The color map illustrates the value of estimated constraint function at every grid point, at different iterations of refining experiments. The Tested Points, Next Test Point, Best feasible Point among Tested Points, and True Optimum are plotted on top of the color map at different iterations of refining experiments.

Figure 7.

Results on dense uniform grid for Schwefel function. The color map illustrates the value of probability of meeting constraints at every grid point, at different iterations of refining experiments. The Tested Points, Next Test Point, Best feasible Point among Tested Points, and True optimum are plotted on top of the color map, at different iterations of refining experiments.

Figure 8.

Results on dense uniform grid for Schwefel function. The color map illustrates the value of estimated probability to be optimum at every grid point, at different iterations of refining experiments. The Tested Points, Next Test Point, Best feasible Point among Tested Points, and True Optimum are plotted on top of the color map, at different iterations of refining experiments.

7. Results on dense uniform grid for Schwefel function

To illustrate the progression of our algorithm we conducted experiments on uniformly sampled fine resolution grid for the synthetic problem of 2D Schwefel function Eq. (12). Here we considered all the points on the grid to be candidate points, instead of hierarchically generating candidate points as described in Algorithm 1.

Figures 5–8 illustrate the value of objective function, estimated constraint function, probability of meeting constraints, and estimated probability to be optimum at the grid points.

The objective function is known, therefore it remains constant throughout the refining experiments. The estimated constraint function, probability of meeting constraints, and estimated probability to be optimum is updated at every refining experiment.

We have also plotted the tested points, next test point, best feasible point among Tested Points, and the true optimum on top of the color maps, to illustrate the progression of the algorithm at different iterations of refining experiments.

Our method does not focus on improving the probability map uniformly on the entire parameter space. From these figures we can see that our method improves the probability map around the promising regions where it believes the true optimum is located. This approach helps to find the optimum with a low number of experiments.

8. Results on physical experiments

We conducted robotic cleaning and sanding experiments to validate our method. The details on the experiments are described in the subsections.

8.1 Robotic scrubbing

In practice, the foreign particles are separated from any target surface by applying cleaning fluid or abrasive actions. Both methods are required in some cases. Our focus was on cleaning hard stains that require scrubbing.

We explored effects of the following robot trajectory parameters on cleaning performance: (1) Tool speed $v$ , (2) Force $f_{n}$ applied by the tool in the normal direction to the surface, (3) Stiffnesses $\kappa_{x}$ and $\kappa_{y}$ in the $x$ and $y$ directions of the tool reference frame, (4) Frequency $f_{x}$ and $f_{y}$ of the forced oscillation in the $x$ and $y$ direction with respect to the tool reference frame, and(5) Amplitude of overlaid force $A_{x}$ and $A_{y}$ for force oscillation. Changes in these parameters change the scrubbing trajectory from nominal trajectory. Figure 9a–c illustrates an example of a nominal trajectory and two actual trajectories generated by two different settings of trajectory parameters.

During robotic scrubbing, the applied force on the surface needs to be low to avoid damaging the part. The user may also want the cleaning process to complete quickly. In our experiments, we defined the objective functions as a weighted sum of normal force on the surface and tool speed. We defined the constraint to be the desired cleaning performance. It is very difficult to model the cleaning performance as a function of the trajectory parameters. Therefore we considered the cleaning performance to be a black-box function of the trajectory parameters.

Let $Y(\theta)$ represent the cleaning performance and $Y_{\rm th}$ be the performance constraint, where $\theta$ $=$ [ $v$ , $f_{n}$ , $\kappa_{x}$ , $f_{x}$ , $A_{x}$ ]. The optimization problem in robotic cleaning can be formulated as following:

Find $\theta^{*}$ s.t. $Y(\theta^{*})\geqslant Y_{\rm th}$ , $\theta^{*}=\mbox{arg}\underset{\theta}{\min}∼{}f(\theta)$ , where $f(\theta)=-(w_{1}\times v+w_{2}\times(1-f_{n}))$ and $w_{1}$ and $w_{2}$ are positive constant weights.

8.1.1 Experimental setup

A KUKA iiwa 7 manipulator was used in our robotic scrubbing experiments. We operated the arm in impedance control mode. Acrylic paint was used as a surrogate for stain and a plastic board was used as the surface to clean. The plastic board was subdivided in smaller blocks of stain. The robot performed cleaning with different sets of trajectory parameters in different blocks.

The experimental setup and an on-going experiment is illustrated in Fig. 9d and e. In our earlier work on parameter optimization [1], we had used small plastic tiles for our robotic cleaning experiments. Changing the slides and painting them individually for each experiment was time consuming. In our current setup the robot performed cleaning on 16 cells before the board needed to be rotated. Cleaning the entire surface from a fixed relative pose between robot and board was not possible due to kinematic constraints of the robot. This approach saved both time and effort as we had to paint only 3 boards and rotate or change the board for 5 times. In our previous setup we had to change the tile at each iteration accumulating a large setup change time.

Figure 9.

(a) Example of a nominal trajectory, (b, c) Actual trajectory originating from the nominal trajectory for different set of robotic operation parameters, (d) The experimental setup for robotic cleaning by scrubbing, (e) An on-going cleaning experiment, (f) An infeasible cleaning performance, (g) A feasible cleaning performance with high objective function value, (h) The feasible cleaning performance with the best found objective function value.

Figure 10.

Convergence to solution in robotic cleaning. Iteration 1–33 are exploration experiments. Iteration 34 to 80 are refining experiments.

8.1.2 Perception

The k-means clustering [47, 48] algorithm was used in our experiments for stain detection. The metric for cleaning performance was the difference in number of stain pixels before and after cleaning.

8.1.3 Result

The cleaning performance threshold was set to be $Y_{\rm th}=0.96$ in our robotic cleaning experiments (we wanted 96% of the stain pixels to be removed). A total of 80 experiments were conducted, including 33 initial exploratory experiments. We found the best feasible solution at the $63^{\rm rd}$ experiment. Examples of some cleaning results from the refining experiments are shown in Fig. 9f–h. We found the optimal parameters to be: $v=0.3$ , $f_{n}=11.25$ N, $\kappa_{x}=3500$ N/m, $f_{x}=3$ Hz, $A_{x}=5$ N. The graph of the best feasible objective function value is illustrated in Fig. 10.

8.2 Bi-manual robotic sanding

Sanding is a process application under the category of surface finishing. For our robotic sanding experiments we have used a rotary tool (Dremel 3000). Our objective was to achieve a desired surface roughness (Ra) value.

We explored effects of the following robot trajectory parameters on sanding performance: (1) Tool speed $v$ , (2) Force $f_{n}$ applied by the tool in the normal direction to the surface, (3) Stiffnesses $\kappa_{z}$ in the normal direction of the tool reference frame, and (4) Stiffnesses $\kappa_{x}$ and $\kappa_{y}$ in the $x$ and $y$ directions of the tool reference frame. In our experiments we considered $\kappa_{x}=\kappa_{y}$ .

Our focus was to minimize force and maximize velocity. Also, we wanted to reduce the Cartesian stiffness in the z-direction to aid tool manipulation on curved surface. We defined the objective functions as a weighted sum of normal force on the surface, tool speed, and stiffness in normal direction of tool. We defined the constraint to be the desired surface roughness. Developing a physics-based model for surface roughness in terms of robot trajectory parameters is not feasible. Therefore surface roughness was the black-box function in these experiments.

Let $Y(\theta)$ represent the surface roughness and $Y_{\rm th}$ be the performance constraint, where $\theta$ $=$ [ $v$ , $f_{n}$ , $\kappa_{z}$ , $\kappa_{x}$ ]. The optimization problem in robotic sanding can be formulated as following:

Find $\theta^{*}$ s.t. $Y(\theta^{*})\geqslant Y_{\rm th}$ , $\theta^{*}=\mbox{arg}\underset{\theta}{\min}f(\theta)$ , where $f(\theta)=-w_{1}\times v+w_{2}\times f_{n}+w_{3}\times\kappa_{z}$ ; $w_{1}$ , $w_{2}$ , and $w_{3}$ are positive constant weights.

8.2.1 Experimental setup

We have used one KUKA iiwa 7 manipulator and one Epson C3 manipulator in the robotic sanding experiments. In the bi-manual operation, the iiwa arm manipulated the tool while the C3 arm held the part. The holding arm changed the pose of the part so that the robot with the tool can reach all the regions on the surface. The advantage of bi-manual operation of such kind is that no human intervention is needed to change the setup. Also, this eliminates the need of a fixed fixture.

Figure 11.

An ongoing robotic sanding experiment.

Figure 12.

Some sanding performance measurement from refining experiments.

Figure 13.

Convergence to solution in robotic sanding. Iteration 1–17 are exploration experiments. Iteration 18 to 124 are refining experiments.

We operated the iiwa arm in impedance control mode. Additively manufactured ABS plastic parts were used in the sanding experiments. The surface roughness of the raw part was 90 $\mu$ m. The tool manipulating robot performed sanding with different sets of trajectory parameters in different regions on the part surface. We considered multiple passes (six passes) of sanding over the surface. Sandpaper of grit size 800 and 1200 were used in the first three passes and last three passes respectively. This selection was based on the opinion of an expert human user. The experimental setup and an ongoing experiment are illustrated in Fig. 11.

8.2.2 Perception

We have used a surface roughness measurement device (Mitutoyo SJ-410) to evaluate the sanding performance at each iteration. Figure 12 illustrates the surface roughness measurement process.

8.2.3 Result

The surface roughness threshold was set to be $Y_{\rm th}=$ 10 $\mu$ m in our robotic sanding experiments. A total of 124 experiments were conducted, including 17 initial exploratory experiments. We found the best feasible solution at the $44^{\rm th}$ experiment. Examples of some sanding results from the refining experiments are shown in Fig. 12. We found the optimal parameters to be: $v=0.01$ , $f_{n}=0.5$ N, $\kappa_{x}=\kappa_{y}=800$ N/m, $\kappa_{z}=400$ N/m. The graph of the best feasible objective function value is illustrated in Fig. 13.

8.3 Experimentation time

Robotic finishing of the entire surface of a part with fairly complex geometry or large surface area usually requires hours to complete. The total time to determine the optimal trajectory parameters needs to be small. We have designed our method such that each experiment will be done on a very small surface patch. Finishing a small surface patch can be done within few seconds. In our cleaning and sanding experiments, each trial took less than 5 seconds to complete. Therefore the entire process of optimization was less than 10 minutes in both cases.

9. Conclusions

In this paper, we described an approach to identify optimal trajectory parameters for robots. Our method optimizes a task objective function and satisfies the task performance constraints. Our iterative algorithm makes decisions based on the uncertainty in the parameter space. This method is applicable to general mathematical problems of constraint satisfaction where the objective function is known and the constraint functions are black-box. We have proven in this paper that our method will converge faster to solution in the expected sense compared to any other method.

We have benchmarked our approach against two established optimization methods on a number of synthetic test problems. Our algorithm converged to a solution faster in all the benchmarking problems. The probability based estimates generated by the surrogate models in our approach can be used as the fitness of the candidate points. Therefore, our method can also be benchmarked in the future with meta-heuristic based approaches that are typically utilized for hyper-parameter optimization in machine learning algorithms [49].

We validated our method by conducting robotic surface finishing (scrubbing and sanding) experiments. Constrained optimization problems in engineering can be solved efficiently using our developed method.

Currently, we have verified our approach on problems with two to ten dimensional parameter spaces and with one to eleven constraints. The benchmarking problems had single known objective functions and multiple black-box constraint functions. We plan to verify our method on higher dimensional problem with multiple black-box objective functions in the future. Our algorithm is designed to find the optimal solution at the finest resolution of discrete parameter space. We plan to extend our algorithm to find optimal solutions in a continuous parameter space in the future.

Footnotes

Acknowledgments

This work is supported in part by National Science Foundation Grant #1634431. Opinions expressed are those of the authors and do not necessarily reflect opinions of the sponsors. and Greg Bradach from Mitutoyo America Corporation for their help with surface roughness measurement using Mitutoyo SJ-410.

Appendix

Table 3

List of symbols

Symbol	Meaning
$d$	Dimension of parameter space
$S$	Parameter space, $S\subset\mathbb{R}^{d}$
$x$	A point in the parameter space (d-dimensional vector)
$a, b$	Lower and upper bounds for the parameter space. $a,b\in\mathbb{R}^{d}$ , $a^{i}\leqslant x^{i}\leqslant b^{i},i=1,\ldots,d$
$f(x)$	Known objective function
$g(x)$	Black-box constraint
$\mathbb{X}_{0}$	Set of initial exploration points
$\mathbb{X}$	Set of tested/evaluated points
$\mathbb{C}_{x}$	Set of candidate points
$n_{c}$	Number of candidate points at (a fixed resolution) in the neighborhood of a tested point
$r_{\textit{initial}}$	Coarse sampling resolution for initial experiments
$r_{x}$	Sampling resolution in the neighborhood of $x\in\mathbb{X}$
$x_{\textit{best}}$	$x_{\textit{best}}\in\mathbb{X}:g_{i}(x_{\textit{best}})\geqslant 0,i=1,\ldots,n$ and $x_{\textit{best}}=\text{argmin}f(x),x\in\mathbb{X}$
$P_{i}^{c}$	Probability of $x_{i}\in\mathbb{C}_{x}$ to satisfy all the constraints
$P_{i}$	Probability of $x_{i}\in\mathbb{C}_{x}$ to be the optimum in $\mathbb{C}_{x}$
$G_{\mathbb{C}_{x}}$	Array of probabilities for each candidate point to meet each constraint
$F_{\mathbb{C}_{x}}$	Array of objective function value for each candidate
$P_{c}$	Array of probabilities for each candidate point to meet all constraints
$Y$	Finishing performance
$Y^{\rm th}$	Desired finishing performance or threshold on finishing performance
$\theta$	Trajectory parameters ( $1\times d$ vector)
$v$	Tool speed
$f_{n}$	Normal force applied by the tool on the surface
$\kappa_{x}$	Stiffness of manipulator along $x$ -axis of tool
$\kappa_{y}$	Stiffness of manipulator along $y$ -axis of tool
$\kappa_{z}$	Stiffness of manipulator along $z$ -axis of tool
$A_{x}$	Amplitude of overlaid force-oscillation along $x$ -axis of tool
$A_{y}$	Amplitude of overlaid force-oscillation along $y$ -axis of tool
$w$	Weight vector for designing objective function in the physical experiments

References

Kabir

Langsfeld

Zhuang

Kaipa

Gupta

. Automated learning of operation parameters for robotic cleaning by mechanical scrubbing. In ASME Manufacturing Science and Engineering Conference, June 2016.

Kabir

Langsfeld

Zhuang

Kaipa

Gupta

. A systematic approach for minimizing physical experiments to identify optimal trajectory parameters for robots. In IEEE International Conference on Robotics and Automation, May 2017.

Ukar

Lamikiz

Lacalle

LNL

Liebana

Etayo

. Laser polishing parameter optimization for die and moulds surface finishing. In ASME International Manufacturing Science and Engineering Conference, October 2008.

Pavel

Wang

Srivastava

. Multi-constraint optimization for grinding nickel-based alloys. In ASME International Manufacturing Science and Engineering Conference, June 2013.

Rao

Pawar

. Parameter optimization of a multi-pass milling process using non-traditional optimization algorithms. Applied Soft Computing2010; 10(2): 445-456.

Mukherjee

Ray

. A review of optimization techniques in metal cutting processes. Computers and Industrial Engineering2006; 50(1-2): 15-34.

Liansheng

Xuemei

Aiping

. The study on the method of film thickness control by coating spray robot for ocean engineering. In ASME International Conference on Ocean, Offshore and Arctic Engineering, June 2015.

Longhurst

Strauss

Cook

. The identification of the key enablers for force control of robotic friction stir welding. In ASME Journal of Manufacturing Science and Engineering2011; 133: 031008.

Asadi

Goldak

Weck

. Welding distortion can be mitigated if welding current and traveling speed vary optimized along a weld path. In ASME Pressure Vessels and Piping Conference, July 2014.

10.

Yamanobe

Fujii

Maeda

Arai

Watanabe

Kato

Sato

Hatanaka

. Optimization of damping control parameters for cycle time reduction in clutch assembly. In IEEE/RSJ International Conference on Intelligent Robots and Systems Aug 2005; 3251-3256.

11.

Zhang

Gravel

Zhang

Wang

Bell

. Robotic force control assembly parameter optimization for adaptive production. In IEEE International Conference on Robotics and Automation (ICRA) May 2011; 464-469.

12.

Cheng

Chen

. Online parameter optimization in robotic force controlled assembly processes. In IEEE International Conference on Robotics and Automation (ICRA) May 2014; 3465-3470.

13.

Ilzarbe

Alvarez

Viles

Tanco

. Practical applications of design of experiments in the field of engineering: A bibliographical review. Quality and Reliability Engineering International2008; 24(4): 417-428.

14.

Zhao

Shi

Lin

Duan

Sun

Zhang

. Surface roughness prediction and parameters optimization in grinding and polishing process for ibr of aero-engine. The International Journal of Advanced Manufacturing Technology2014; 74(5): 653-663.

15.

Ren

Cabaravdic

Zhang

Kuhlenkötter

. A local process model for simulation of robotic belt grinding. International Journal of Machine Tools and Manufacture2007; 47(6): 962-970.

16.

Alagumurthi

Palaniradja

Soundararajan

. Optimization of grinding process through design of experiment (doe) – a comparative study. Materials and Manufacturing Processes2006; 21(1): 19-21.

17.

Kumar

KVBSK

Choudhury

. Investigation of tool wear and cutting force in cryogenic machining using design of experiments. Journal of Materials Processing Technology2008; 203(1-3): 95-101.

18.

Lahaye

Chomont

Dumont

Duchesne

Chabassier

. Using a design of experiment method to improve kdp crystal machining process. International Society for Optics and Photonics, 1999.

19.

Wang

Zhao

Luo

. Optimization of design of experiment for chemical mechanical polishing of a 12-inch wafer. Microelectronic Engineering2013; 112: 5-9.

20.

Liao

Shie

Yang

. Applications of taguchi and design of experiments methods in optimization of chemical mechanical polishing process parameters. The International Journal of Advanced Manufacturing Technology2008; 38(7): 674-682.

21.

Tam

Loh

Mah

CPA

Loh

. Electrochemical polishing of biomedical titanium orifice rings. Journal of Materials Processing Technology1992; 35(1): 83-91.

22.

Gravel

Zhang

Bell

Zhang

. Objective metric study for doe-based parameter optimization in robotic torque converter assembly. In 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems Oct 2009; 3832-3837.

23.

Zhang

Bell

Zhang

Wang

Martinez

. On-pendant robotic assembly parameter optimization. In 2008 7th World Congress on Intelligent Control and Automation June 2008; 547-552.

24.

Elangovan

Prakasan

Jaiganesh

. Optimization of ultrasonic welding parameters for copper to copper joints using design of experiments. The International Journal of Advanced Manufacturing Technology2010; 51(1): 163-171.

25.

Bement

. Taguchi techniques for quality engineering. Technometrics1989; 31(2): 253-255.

26.

Rowlands

Antony

. Application of design of experiments to a spot welding process. Assembly Automation2003; 23(3): 273-279.

27.

Liao

Shie

. Optimization on selective laser sintering of metallic powder via design of experiments method. Rapid Prototyping Journal2007; 13(3): 156-162.

28.

Rasmussen

. Gaussian processes for machine learning. 2006.

29.

Pelikan

Goldberg

Cantú-Paz

. Boa: The bayesian optimization algorithm. In Proceedings of the 1st; Annual Conference on Genetic and Evolutionary Computation, GECCO’99, San Francisco, CA, USA, Morgan Kaufmann Publishers Inc1999; 1: 525-532.

30.

. Improving efficiency with orthogonal exploration for online robotic assembly parameter optimization. In IEEE International Conference on Robotics and Biomimetics (ROBIO) Dec 2015; 958-963.

31.

Barton

Meckesheimer

. Chapter 18 metamodel-based simulation optimization. In Henderson SG, Nelson BL, (eds). Simulation, volume 13 of Handbooks in Operations Research and Management Science, Elsevier2006; 13: 535-574.

32.

Marvel

Newman

Gravel

Zhang

Wang

Fuhlbrigge

. Automated learning for parameter optimization of robotic assembly tasks utilizing genetic algorithms. In 2008 IEEE International Conference on Robotics and Biomimetics Feb 2009; 179-184.

33.

Pawar

Rao

. Parameter optimization of machining processes using teaching – learning-based optimization algorithm. The International Journal of Advanced Manufacturing Technology2013; 67(5): 995-1006.

34.

Rao

Kalyankar

. Parameter optimization of machining processes using a new optimization algorithm. Materials and Manufacturing Processes2012; 27(9): 978-985.

35.

Le Digabel

. Algorithm 909: NOMAD: Nonlinear optimization with the MADS algorithm. ACM Transactions on Mathematical Software2011; 37(4): 1-15.

36.

Audet

Dennis Jr

. Mesh adaptive direct search algorithms for constrained optimization. SIAM Journal on Optimization2006; 17(1): 188-217.

37.

Abramson

Audet

. Convergence of mesh adaptive direct search to second-order stationary points. SIAM Journal on Optimization2006; 17(2): 606-619.

38.

Regis

. Stochastic radial basis function algorithms for large-scale optimization involving expensive black-box objective and constraint functions. Computers & Operations Research2011; 38(5): 837-853.

39.

Montgomery

. Design and Analysis of Experiments. John Wiley & Sons, Inc, 2013.

40.

Michalewicz

Schoenauer

. Evolutionary algorithms for constrained parameter optimization problems. Evolutionary Computation1996; 4(1): 1-32.

41.

Karaboga

Akay

. A modified artificial bee colony (abc) algorithm for constrained optimization problems. Applied Soft Computing2011; 11(3): 3021-3031.

42.

Coello

Montes

. Constraint-handling in genetic algorithms through the use of dominance-based tournament selection. Advanced Engineering Informatics2002; 16(3): 193-203.

43.

Dvorkin

Goldschmit

Storti

. Constrained optimization problems in mechanical engineering design using a real-coded steady-state genetic algorithm. Mecánica Computacional2010; XXIX: 9287-9303.

44.

Floudas

Pardalos

. A collection of test problems for constrained global optimization algorithms, volume 455. Springer Science & Business Media, 1990.

45.

Kramer

. An augmented lagrange multiplier based method for mixed integer discrete continuous optimization and its applications to mechanical design. Journal of Mechanical Design1994; 116: 405.

46.

Deb

Datta

. A fast and accurate solution of constrained optimization problems using a hybrid bi-objective and penalty function approach. In IEEE Congress on Evolutionary Computation July 2010; 1-8.

47.

Arthur

Vassilvitskii

. k-means++: The advantages of careful seeding. In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, Society for Industrial and Applied Mathematics2007; 1027-1035.

48.

Hartigan

. Clustering algorithms. 1975.

49.

Bergstra

Bengio

. Random search for hyper-parameter optimization. J Mach Learn Res February 2012; 13: 281-305.

Identifying optimal trajectory parameters in robotic finishing operations using minimum number of physical experiments

Abstract

Keywords

1. Introduction

2.1 Design of Experiments (DOE) based approaches

2.2 Surrogate model based approaches

2.3 Meta-heuristics based approaches

2.4 Constrained optimization under black-box constraints

3. Problem formulation

Table 1 Parameters values used in our design of exploratory experiments

4.5 Optimality in continuous space vs discrete space

5. Theoretical analysis of convergence

.

6.2 Optimization methods for benchmarking

6.3 Results on synthetic test problems

8. Results on physical experiments

8.1 Robotic scrubbing

8.1.1 Experimental setup

8.1.3 Result

8.2 Bi-manual robotic sanding

8.2.1 Experimental setup

8.2.3 Result

8.3 Experimentation time

9. Conclusions

Footnotes

Acknowledgments

Appendix

References

Table 1
Parameters values used in our design of exploratory experiments