The Higher Education and Research forge

Home My Page Projects Code Snippets Project Openings Complex Surface Machining Optimization
Summary Activity SCM

SCM Repository

1 %\documentclass[review]{elsarticle}
2 \documentclass[a4paper,3p, twocolumn]{elsarticle}
4 \usepackage{soulutf8}
5 \usepackage[dvipsnames, svgnames]{xcolor}
6 \usepackage[ruled,vlined]{algorithm2e}
7 \usepackage[colorlinks=true]{hyperref}
8 \usepackage{amssymb}
9 \usepackage{tabularx}
10 \usepackage{amsmath}
12 \journal{Journal of \LaTeX\ Templates}
13 \bibliographystyle{elsarticle-num}
15 \def\thename{C4DO1D}
17 \begin{document}
18 \newenvironment{itemize*}%
19   {\begin{itemize}%
20     \setlength{\itemsep}{0pt}%
21     \setlength{\parskip}{0pt}}%
22   {\end{itemize}}
25 \begin{frontmatter}
26 \title{
27     \vspace*{-0.2cm}
28 \begin{tabularx}{\textwidth}{Xr}
29  \includegraphics{elsevier.png}
30  % elsevier.png: 898x992 px, 1200dpi, 1.90x2.10 cm, bb=0 0 54 60
31  &
32  \includegraphics{cirp.png} \\
33  % cirp.png: 560x414 px, 500dpi, 2.84x2.10 cm, bb=0 0 81 60
34  ~ & ~
35 \end{tabularx}
36 {\small 14th CIRP Conference on Intelligent Computation in Manufacturing Engineering, Gulf of Naples, Italy}
37 \begin{tabularx}{\textwidth}{Xr}
38  ~ & ~
39 \end{tabularx}
40 Toolpath planning optimization for end milling of free-form surfaces using a clustering algorithm
41 }
44 %% or include affiliations in footnotes:
45 \author[mymainaddress]{Mahfoud Herraz}
46 \author[mymainaddress]{Jean-Max Redonnet\corref{mycorrespondingauthor}}
47 \cortext[mycorrespondingauthor]{corresponding author}
48 \ead{jean-max.redonnet@univ-tlse3.fr}
49 \author[mysecondaryaddress]{Mohammed Sbihi}
50 \author[mysecondaryaddress]{Marcel Mongeau}
52 \address[mymainaddress]{Institut Cl\'ement Ader, Universit\'e Paul Sabatier, Toulouse, France}
53 \address[mysecondaryaddress]{ENAC, Universit\'e de Toulouse, France}
55 \begin{abstract}
56 A two-step approach to address the toolpath planning optimization problem is introduced.
57 The first step divides the surface into zones with similar local geometric properties in order to improve efficiency of toolpath planning algorithms. To do this, an unsupervised clustering algorithm (like {\em K-means}) is used on a mesh defined by isoparametric curves.
58 The second step applies a toolpath planning algorithm to each zone, according to its optimal machining direction. This optimal direction is calculated using a black-box optimization software. The final goal is to  enhance gradually the formulation of the optimization problem to obtain better and better results.
59 \end{abstract}
61 \begin{keyword}
62 end-milling \sep free-form surface \sep toolpath planning \sep clustering \sep black-box optimization
63 \end{keyword}
65 \end{frontmatter}
67 \section{Introduction}
68 Tool path planning is a very important issue for the end milling of free-form surfaces, due to the high machining cost commonly observed for manufacturing high added value parts such as molds and stamping dies. Therefore, optimizing the toolpath planning process may lead to significant gains in terms of machining costs.
70 In order to make progress in optimizing toolpath planning of free-form surface machining, a two-step approach, that we shall call \thename\footnote{\thename{} stands for Clustering 4D Optimization 1D}(pronounced ``CADOID''), is introduced in this paper. Because it is well established that toroidal cutter may lead to better results than its ball-end counterpart \cite{bedi_toroidal_1997}, this type of cutter have been chosen to carry on this study. Actually, a toroidal cutter provides better results than a ball-end one when machining along the steepest-slope direction, but when used perpendicularly to the steepest-slope direction, its performances are worse than those of a ball-end cutter. Because the steepest-slope direction may vary a lot across a free-form surface, partitioning this surface into several zones (each of which will then be machined along an appropriate direction) is considered as the best approach to improve efficiency of the toroidal-cutter choice \cite{senatore_corr_2012}.
72 \section{Partitioning the surface by a clustering algorithm}
73 In this section, a new method to define zones suitable for machining with a toroidal cutter is presented. This method relies on a clustering approach. Clustering algorithms are part of unsupervised machine learning algorithms. From an initial set of sample points, they provide a predefined number, $K$, of clusters gathering sample points according to a given metric. One of the most commonly used clustering algorithm is the {\em K-means} algorithm. 
75 \subsection{Formulating the clustering problem}
76 In order to apply the K-means algorithm to the surface partitioning problem, a set of sample points must first be defined. To do this, a mesh over the entire surface $\mathbf{S}(u, v) $ is defined by a set of $(u,\,v)$ isoparametric curves. This way, a set of elementary meshes is also defined. Taking the center of each elementary mesh as a sample point provides a suitable set of sample points ensuring the entire surface is covered (Figure \ref{fig:sample}). In other words, the sample points $\mathbf{S_i}$, $i \in \{1,\dots,n\}$, are defined such that $\mathbf{S_i}=\mathbf{S}(u_i,v_i)$, where $u_i$ and $v_i$ are parametric values discretizing the whole parametric domain.
77 \begin{figure}[ht]
78  \centering
79  \def\svgwidth{\linewidth}
80  \input{sampling.pdf_tex}
81  % sampling.pdf: 351x223 px, 72dpi, 12.38x7.87 cm, bb=0 0 351 223
82  \caption{Sampling-point definition}
83  \label{fig:sample}
84 \end{figure}
86 The second element required to implement the K-means algorithm is a metric space adapted to the sought objective. A metric space consists of two components: a vector of significant \textit{features}, and the definition of distance between two vectors: that we shall call the \textit{metric}. Actually, a metric provides a measure of the dissimilarity between a pair of sample points. Thus, two points that are similar, have a reduced dissimilarity distance, while two points that are very different have a greater distance. As a first naive approach, the dissimilarity can be evaluated directly through the Euclidean distance between the two 3D points $\mathbf{S_i}$ and $\mathbf{S_j}$. However, this distance does not totally represent the nearness of these sample points from a milling point of view. In fact, to be suitable for milling with a toroidal cutter, one expects the zones proposed by clustering algorithm to be \textit{connected} (in the mathematical sense, {\em i.e.}, roughly speaking, a zone must not show any discontinuity). Indeed, a discontinuity in a zone leads to a tool lifting that is harmful to the total machining time. Thus, some definition of geometric proximity of the sample points should intervene in the chosen metric. The parametric values $u_i$ and $v_i$ are well suited to give such an insight.
88 Moreover, in order to improve the efficiency of the toroidal cutting tool, the resulting zones are expected to feature rather uniform steepest slopes. In other words, the steepest-slope value and its orientation should vary as little as possible within a single zone. Thus, these two values should also be taken into account in the chosen metric.
90 To summarize, for each sample point $\mathbf{S_i}$, $i \in \{1,\dots,n\}$, we define the four-component feature vector $\mathbf{F_i} \in \mathbb{R}^4$ whose four components are (Figure \ref{fig:features}):
91 \begin{itemize}
92  \item the $u_i$ parametric value of the sample point $\mathbf{S_i}$
93  \item the $v_i$ parametric value of the sample point $\mathbf{S_i}$
94  \item the steepest-slope angle, noted $s_i$. This is the angle between the vector $\mathbf{n_i} = \mathbf{n}(u_i,v_i)$, normal to the surface at $\mathbf{S_i}$, and the horizontal plane $(\mathbf{X},\mathbf{Y})$.
95  \item the steepest-slope orientation, noted $\theta_i$. This is the angle between the projection of the vector $\mathbf{n_i}$ onto the horizontal plane $(\mathbf{X},\mathbf{Y})$, and the $\mathbf{X}$ axis.
96 \end{itemize}
97 \begin{figure}[ht]
98  \centering
99 % \def\svgwidth{\linewidth}
100  \input{features.pdf_tex}
101  % sampling.pdf: 351x223 px, 72dpi, 12.38x7.87 cm, bb=0 0 351 223
102  \caption{The four components ($u_i$, $v_i$, $s_i$ and $\theta_i$) of the feature vectors}
103  \label{fig:features}
104 \end{figure}
106 To obtain a fully defined metric space, a distance calculation method (the metric) must also be defined. The most natural choice is the Euclidean distance in $\mathbb{R}^4$, the space of the chosen parameters ($u$, $v$, $s$ and $\theta$). The dissimilarity between two sample points $\mathbf{S_i}$ and $\mathbf{S_j}$ is thereby defined as:
107 \begin{multline}\label{eq:dist}
108  d(i,j) = \|\mathbf{F_i}-\mathbf{F_j} \| = \left((u_i-u_j)^2+(v_i-v_j)^2 \right. \\
109  \left.+(s_i-s_j)^2+(\theta_i-\theta_j)^2\right)^\frac{1}{2}
110 \end{multline}
111 The clustering problem  aims at partitioning the given $n$  sample points (or, to be more precise, the $n$ feature vectors) into $K$ $(\leq n)$ clusters $Z_1, Z_2, \dots, Z_K$  so as to minimize the sum of intra-cluster distances (squared, simply to avoid discontinuity in the derivative of the objective function). More precisely, the objective is to find $(Z_1,\dots,Z_K) \in \mathcal{P}$ minimizing:
112 $$
113 \sum_{k=1}^K\sum_{i \in Z_k}\|\mathbf{F_i}-\mathbf{\bar{F}_k} \|^2
114 $$
115 where $\mathcal{P}$ denotes the set of all possibles \textit{K-partitions} (into K sets) of the sample-point index set $\{1,\dots,n\}$, and $\mathbf{\bar{F}_k}=\frac{1}{|Z_k|}\sum_{i \in Z_k}\mathbf{F_i}$ (called the \textit{centroid} of $Z_k$) is the barycenter of the features vectors of all sample points belonging to $Z_k$.
117 \subsection{Running the K-means algorithm on the 4D parameter vectors}
118 Once the sample-point dataset with their feature vectors and the metric defined, the K-means algorithm may be carried out. Let $K$ be the number of clusters chosen. The K-means algorithm  is described in Algorithm~\ref{alg:Kmeans}.
119 \begin{algorithm}
120  \DontPrintSemicolon
121  \SetAlgoLined
122  \tcc{Initialize centroids}
123  choose $K$ initial centroids $\mathbf{\bar{F}_1},\dots, \mathbf{\bar{F}_K}$ among the $n$ feature vectors\;
124  \Repeat{no more change occur}{
125     \tcc{allocate points to clusters}
126     \ForEach{sample point $\mathbf{S_i}$}
127     {attach $\mathbf{S_i}$ to the cluster $k$ corresponding to $\mathrm{argmin}_k \|\mathbf{F_i}- \bar{\mathbf{F}}_k\|$}
128    \tcc{cluster-centroid update}
129     \ForEach{cluster $k$}{recalculate centroid $\mathbf{\bar{F}_k}$ as the barycenter of cluster-$k$ feature vectors\;}
130  }
131  \caption{The K-means algorithm}\label{alg:Kmeans}
132 \end{algorithm}
134 In general, the resulting K-partition of the K-means algorithm depends on the particular choice of the initial centroids. Empirical tests showed however that in the case of machinable surfaces, the resulting zones are almost not sensitive to the choice of the initial centroids. This is due to the (relative) regularity of machinable surfaces and the large (enough) number of sampling points, so that outliers are avoided and the final result is less sensitive to initialization. In this study, an equally-spread angular repartition, in the parametric space, of initial centroids is adopted.
136 The running loop of the algorithm is divided into two parts. In the allocation step, each sample point $\mathbf{S_i}$ is attached to the closest centroid, according to the metric distance (\ref{eq:dist}), {\em i.e.} a sample point $\mathbf{S_i}$ is attached to the cluster $k$ corresponding to $\mathrm{argmin}_k \|\mathbf{F_i}- \bar{\mathbf{F}}_k\|$. In the update step, the centroids are recalculated to reflect the allocation changes. To do this, for each cluster $Z_k$, the 4D barycenter of the feature vectors $(u_i$, $v_i$, $s_i$, $\theta_i)$, $i\in Z_k$, of the sampling points belonging to the cluster $Z_k$, is calculated, which gives a new centroid,  $\bar{\mathbf{F}}_k$. Remark that this new centroid does not generally correspond to one of the $n$ given sample points. Finally the repeat loop ends when stabilization occurs, {\em i.e.} no more sample point is switching from a cluster to another one and the position of centroids remain constant.
138 Figure \ref{fig:choi3} shows an example of surface partitioning obtained on the surface used in \cite{keun_choi_tool_2007} with $K=3$ and a $200\times200$ tessellation grid.
139 \begin{figure}[ht]
140  \centering
141  \includegraphics[width=\columnwidth]{choi200x200.png}
142  % choi200x200.png: 1517x885 px, 300dpi, 12.84x7.49 cm, bb=0 0 364 212
143  \caption{Example of a surface partitioned into $K=3$ zones}
144  \label{fig:choi3}
145 \end{figure}
147 \subsection{Dividing disconnected zones}
148 Carrying out the K-means algorithm is fast, even for fine tessellation meshes (for example the partitioning scheme presented in Figure \ref{fig:choi3} is fully computed within only 684~ms on a 16 cores Intel(R) Xeon(R) 1.70GHz system). However, it may happen that the above procedure provides disconnected zones, like in the case depicted by Figure \ref{fig:ccs}, where the green zone is composed of two disconnected sub-zones. This green zone must then be further partitioned into two zones, to prevent the toolpath planning algorithm to provide a machining process involving a single run, without lifting the cutter, which may cause severe damages to the surface, and lead to unexpected results.
149 \begin{figure}[ht]
150  \centering
151  \includegraphics[width=\columnwidth]{ccsMap.png}
152  % ccsMap.png: 1218x1028 px, 300dpi, 10.31x8.70 cm, bb=0 0 292 247
153  \caption{Example of disconnected zone (in green)}
154  \label{fig:ccs}
155 \end{figure}
157 We therefore propose a connectivity-check procedure to ensure that the zones provided by the above K-means algorithm are connected. The connected-component search algorithm we are proposing to carry out this task goes as follows. Finding connected components is a well-known polynomial problem in graph theory, but to be able to use it with a partitioning scheme, an undirected graph corresponding to this partitioning scheme must first be defined. Each sample point corresponds to a vertex of the graph. Any two sample points belonging to a same cluster are connected with an edge if and only if they come from side-by-side elementary meshes. Once the graph corresponding to the partitioning scheme is defined, a connected-component search can be performed, {\em e.g.}, using the software package  JGraphT \cite{jgrapht}. As a result, any disconnected zone will simply be sub-divided into several connected zones. The time needed to perform this task is heavily dependant on the number of sample points. However, even if it is a bit longer than the K-means procedure itself, it is still short enough to be used interactively. For example, the connected-component search on the surface presented in Figure \ref{fig:ccs} requires less than 15~seconds of CPU time (on the same computer as above). As expected, for this particular case, the procedure results in four different zones displayed on Figure \ref{fig:ccsDone}.
158 \begin{figure}[ht]
159  \centering
160  \includegraphics[width=\columnwidth]{ccsDone.png}
161  % ccsDone.png: 1228x1026 px, 300dpi, 10.40x8.69 cm, bb=0 0 295 246
162  \caption{Result of the connected component search}
163  \label{fig:ccsDone}
164 \end{figure}
166 To summarize, using the proposed connected-component search strategy appears to be a fast an efficient choice to deal with the zone-partitioning issue.
168 \section{Optimizing machining direction within each zone}
169 Once the zone-partitioning scheme defined, the optimal milling process for machining the resulting $K$ zones still needs to be determined. For a given zone, let us call $\mathbf{X}$ the vector of the (numerous) variables that can have an impact on the machining time, and $\mathbf{P_i(\mathbf{X})}$, the $i^{th}$ point of the toolpath ($i=1 \ldots n(\mathbf{X})$). For each zone, the optimization problem can then be expressed as:
170 \begin{equation}\label{eq:optProbGlobal}
171 \begin{aligned}
172   \min_{\mathbf{X}} \quad & machining\_time(\mathbf{X}) \\
173   \textrm{s.t.} \quad & sh_i(\mathbf{X}) \leqslant sh_{max}, \quad i=1 \ldots n(\mathbf{X})
174 \end{aligned}
175 \end{equation}
176 where:
177 \begin{itemize}
178     \item $sh_i(\mathbf{X})$ is the scallop height at the vicinity of the point $\mathbf{P_i}(\mathbf{X})$
179     \item $sh_{max}$ is the maximum scallop height allowed for this surface.
180 \end{itemize}
182 From an optimization point of view, this is a fairly difficult problem as it is complicated to characterize the entire search space. Moreover, even for a given path, the objective function ($machining\_time$) cannot be calculated analytically ({\em i.e.} as an explicit mathematical formula in terms of optimization variables $\mathbf{X}$). Actually, the objective function is the result of a computer simulation. Therefore, the machining time has to be considered as a so-called \textit{black-box} objective function involving numerous variables and constraints. In a preliminary approach, some choices can be made to simplify the objective function. First of all, we restrict our study to 3-axis machines. Second a given toolpath generation strategy is chosen: the strategy of {\em zig-zag} parallel planes. This strategy is widely used in the industry as it ensures that the entire surface is covered and it is easy to implement. The paths are thereby defined as the intersection of surface of the part and parallel vertical planes. The distance between two adjacent planes (called step-over distance, noted $sod$)  is defined as the maximum distance such as the scallop-height constraint (\ref{eq:optProbGlobal}) is respected for the whole path. This choice facilitates the management of the quality criterion. Indeed, at a given plane, the position of the adjacent plane is easy to determine: it suffices to consider the worst point, in terms of scallop height, of the path defined by the given plane. This is enough to ensure that the quality criterion will be satisfied for the entire trajectory. From the optimization point of view, this strategy is very useful because the constraints are thereby satisfied by construction. Using a fast scallop height calculation method, such as in \cite{redonnet_study_2013}, the above-described toolpath planning process is rapid enough to be used as the objective function of an optimization procedure.
184 The only parameter that remains to be defined is, for each zone, the direction of the parallel planes used to run the machining strategy. For the zone $Z_k$, this angle-determination problem may be summarized as:
185 \begin{equation}\label{eq:optProbSimple}
186   \min_{\Theta_k} \quad machining\_time(\Theta_k)
187 \end{equation}
188 where $\Theta_k$ is the angle defining the machining direction of zone $Z_k$ with reference to the $\mathbf{X}$ axis.
190 This is a univariate (one-dimensional) optimization problem whose objective function is, again,  a black-box ({\em i.e.}, its analytic form is not known). Thus, one single evaluation of this function requires a full toolpath planning simulation. Therefore,  the optimization process may be time expensive, especially if running  the whole procedure within a time suitable for interactive usage is required. Such a problem can be solved using a simple optimization algorithm, like the popular Nelder-Mead algorithm, used in \cite{redonnet_optim_2018} to solve (\ref{eq:optProbSimple}) over the entire surface considered as a single zone. However, given the previous considerations on the computation time, in the present study, a more efficient black-box optimization approach is chosen. An extensive review of black-box optimization methods can be found in \cite{AuHa2017a}. The software package used in the present study to perform optimization is {\em NOMAD} \cite{Nomad, AuLeTr09a}. It is a LGPL-licence software that relies on the algorithms described in \cite{Le2011a, AuDe2006}. Using this package makes the resolution more flexible, and suitable for parallel computing.
192 In order to run the optimization procedure an initial guess of the angle $\Theta_k$ is advised. For this, we choose the $\theta$ value of the initial centroid of zone $Z_k$, as it may be interpreted as the overall steepest-slope direction for the zone $Z_k$. As the steepest-slope direction is the most effective machining direction when using a toroidal cutter, this choice should be close to the optimal direction.
194 Applying this optimization procedure on the surface studied in \cite{keun_choi_tool_2007} (see also Figure \ref{fig:choi3}) provides the toolpath planning depicted in figure \ref{fig:choiToolpath}.
195 \begin{figure}[ht]
196  \centering
197  \includegraphics[width=\columnwidth]{choi_toolpath.png}
198  % choi_toolpath.png: 1525x924 px, 300dpi, 12.91x7.82 cm, bb=0 0 366 222
199  \caption{Example of toolpath optimization}\label{fig:choiToolpath}
200 \end{figure}
202 To carry out the parallel-planes strategy, a cutter of outer radius 5~mm and torus radius 2~mm is used, and the maximum scallop height is set to 0.01~mm, a value commonly used in industry.
204 Using the optimized toolpath, the two symmetric zones of Figure \ref{fig:choiToolpath} need about 36~s to be machined, while the front zone requires about 28~s. Thus, the total machining time is around 100~s. This value should be considered keeping in mind the bulk size of the part is about $50\times75$~mm on the $(\mathbf{X},\mathbf{Y})$ plane.
206 \section{Discussion}
207 We perform empirical tests comparing K-means with other clustering algorithms, namely the {\em Hierarchical clustering} \cite{szekely_hierarchical_2005} and the {\em Rival penalized competitive learning} \cite{xu_rival_1993}. These tests reveals that the latter are more sensitive to algorithm parameters, such as learning coefficients and mesh size, and that the K-means algorithm provides the best results.
209 The choice of the particular metric used in  the clustering algorithm has a deep impact on the quality of the results. This choice deserves thereby attention in a deeper study.
211 The optimization process duration is heavily dependant on the implementation parameters chosen. It ranges from a few seconds of CPU time to several minutes, depending on: whether the surrogate-model option is activated or not, which stopping criterion is chosen, the mesh size, and so on. More work needs to be pursued to find out a good trade-off between CPU time and quality of the results.
213 Toolpath planning for milling of free-form surfaces is a field of research where very few authors release enough data to allow relevant comparisons between methods. Nevertheless, a comparison can be conducted between the \thename{} method introduced in this paper and the one proposed by Choi \& Banerjee in \cite{keun_choi_tool_2007} because numerical result values are published in their paper. However, this comparison should be considered with caution because the cutter used in their paper is a ball-end cutter. Remark also that the measurement units in their paper are inches, thus conversion to millimeters gives unusual values. The cutter radius used in~\cite{keun_choi_tool_2007} is $0.125~\mathrm{in} = 3.175~\mathrm{mm}$; thus we choose a toroidal cutter with same outer radius to carry out the comparison. The torus radius was set to $0.05~\mathrm{in} = 1.27~\mathrm{mm}$. Choi \& Banerjee present two experimental simulations: one with a maximum scallop height $sh_{max}=0.01~\mathrm{in} = 0.254~\mathrm{mm}$, and the other using $sh_{max} = 0.05~\mathrm{in} = 1.27~\mathrm{mm}$. Besides, the machining time is not precised in \cite{keun_choi_tool_2007}; thus, we only compare total toolpath lengths. The resulting values are summarized in table \ref{tab:cmp}.
214 \begin{table}[ht]
215 \centering
216 \begin{tabular}{|l|c|c|}
217     \hline
218     $sh_{max}$ & \thename & Choi \& Banerjee \\
219     \hline
220   0.254  & 1437  & 1717  \\
221   1.27  & 323 & 955 \\
222     \hline
223 \end{tabular}\\
224 \vspace{1pt}
225 \raggedleft
226 \begin{footnotesize}units: millimeters\end{footnotesize}
227 \caption{Comparison results with Choi \& Banerjee's method}\label{tab:cmp}
228 \end{table} 
230 Computation time is not specified in \cite{keun_choi_tool_2007}, but simply to give a rough idea, running the full \thename{} method (clustering and optimization) takes less than 1 min of CPU time on a 16 cores Intel(R) Xeon(R) 1.70GHz system.
232 These comparative results show a significant improvement of performances using the proposed clustering/black-box optimization method. But once again, they should be considered with caution because the surface is rather small and the maximum scallop height values are relatively high. It would be interesting to compare both methods on parts that are more representative of industrial ones, and using parameters commonly used in industry.
234 \section{Conclusion}
235 In the context of machining free-form surfaces, the toolpath optimization problem is a challenge. The \thename{} approach introduced in this paper provides the optimal solution (given the restricting choices made). First, a K-means algorithm is carried out to define machining zones. A zone segmentation verification is then performed, and finally, for each zone, a black-box optimization procedure is carried out to define an optimal machining direction.
237 This work is a first attempt to address the toolpath planning optimization problem using a clustering/black-box optimization approach. Future tracks of research include improvements of both the clustering step and the optimization step. In regard to this latter point, the choice of a black-box optimization software may be very helpful to integrate, step by step, new degrees of complexity into the optimisation model.
239 \bibliography{biblio}
241 \end{document}