Optimal parameters identification and sensitivity study for abrasive waterjet milling model

Abstract

In this paper we present the work related to the parameters identification for abrasive waterjet milling (AWJM) model that appears as an ill-posed inverse problem. The necessity of studying this problem comes from the industrial milling applications where the possibility to predict and model the final surface with high accuracy is one of the primary tasks in the absence of any knowledge of the model parameters that should be used. The adjoint approach based on corresponding Lagrangian gives the opportunity to find out the unknowns of the AWJM model and their optimal values that could be used to reproduce the required trench profile. Due to the complexity of the nonlinear problem and the large number of the model parameters, we use an automatic differentiation software tool. This approach also gives us the ability to distribute the research on more complex cases and consider different types of model errors and 3D time dependent model with variations of the jet feed speed. This approach gives us a good opportunity to identify the optimal model parameters and predict the surface profile both with self-generated data and measurements obtained from the real production. Considering different types of model errors allows us to receive the results acceptable in manufacturing and to expect the proper identification of unknowns.

Keywords:

• Inverse problems
• abrasive waterjet
• PDE constrained optimization
• Tikhonov regularization
• automatic differentiation
• numerical analysis
• parameters estimation

AMS Subject Classifications:

• 35B30
• 35Q93
• 35R30
• 49Kxx
• 65K10

1. Introduction

One of the crucial problems in the variety of natural, industrial and economical phenomena which could be modeled by (systems of) partial differential equations (PDEs), is the identification of model parameters by suitable comparison between the experimental measurements referred to the real systems and the predictions of the mathematical models.

In order to solve the direct problem, it is necessary to know in advance all the involved parameters such as coefficients or sources characteristics, but usually different model parameters are unknown or inaccessible and have to be identified only from experimental measurements. However, parameter identification is generally an ill-posed inverse problem [1,2], even if the PDEs are linear under some considerations due to the measurement noises and modeling errors. This kind of problems can be overcome by using so-called regularization methods [3–5], studied by many authors [6–9], and parameters identification could be reformulated in a stable case as a minimization problem with a data mismatch and a regularization term.

In this paper we present the mathematical method to identify unknown parameters of the proposed generic abrasive waterjet milling (AWJM) model which was previously described and studied in [10–13] and was developed according to the industrial needs for waterjet footprints prediction.

Assuming that model parameters and source terms are known, one can find the shape of the trench profile. In this context this is the forward problem that in our case involves a nonlinear PDE. Some studies of linear AWJM inverse problems were previously presented in [14].

Our goal is to determine the model parameters inaccessible from the experiments to predict the surface construction before performing the manufacturing simulations. Following this requirement we pose the minimization problem, that could be presented as searching the minimum of a cost function – the difference between the observations and the corresponding model solution. This functional can be minimized using a numerical optimization method based on the limited memory Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm [15–17] underlying the N2QN1 minimization package from the INRIA MODULOPT library [18]. This is a gradient descent algorithm which is an iterative process requiring the gradient of the cost function. In our work the gradient vector is obtained numerically using the automatic differentiation (AD) software TAPENADE [19], which can be interpreted as the Lagrange multiplier of the model equations in terms of adjoint problem [20,21].

By adding noise to the artificial data, we show that in fact the parameter identification problem is highly unstable and strictly depends on input measurements. Furthermore, we demonstrate that Tikhonov regularization could be effectively used to deal with the presence of data noise and to improve the identification correctness.

The paper is organized as follows. In the first section we present the proposed mathematical model of the problem, explain the adjoint approach based on corresponding Lagrangian that was used to get the gradient of the cost function, needed in the implementation of the gradient descent L-BFGS based algorithm. We also provide a short list of gradient descent algorithms with some explanations. Section 3 represents the parameters identification for linearized 2D cases, based on self-generated input data as well as on real experiment measurements. Section 4 contains the results corresponding to the sensitivity study of the model in case of various measurement errors and Section 5 finalizes this paper with some conclusions and perspectives of future research.

2. Mathematical model

2.1. Proposed model

The model of jet footprints is defined as a nonlinear PDE that describes the forming of the surface for any jet feed speed and milling direction regardless of the target material types. We consider the time interval [0, T] and denote by $\mathrm{\Omega }$ a bounded domain of ${\mathbb{R}}^2$.

The AWJM model introduced in [10–12] is described as:(1) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\frac{\mathit{E}(x,y){\mathrm{e}}^{a\mathit{Z}}}{{\left(1+{\left|\mathrm{\nabla }\mathit{Z}\right|}^2\right)}^{k/2}}\text{in}\mathrm{\Omega }\times [0,T],\end{array}$$(1)

where

• $(x,y)\in \mathrm{\Omega }$

• $\mathit{Z}(x,y,t)\in {\mathbb{H}}_0^1(\mathrm{\Omega },{\mathbb{R}}^{+})$ is the parametrization of the surface,

• $a,k\in {\mathbb{R}}^{+}$ are the model parameters,

• $\mathit{E}(x,y)\in {\mathbb{L}}^2(\mathrm{\Omega })$ is the Etching rate function which is also model parameter,

with initial and boundary conditions:$$\left\{\begin{array}{c}\mathit{Z}(x,y)=0\text{on}\mathit{\partial }\mathrm{\Omega },\hfill \\ \mathit{Z}(x,y)={\mathit{Z}}_0\text{at}t=0.\hfill \end{array}\right.$$

Further the AWJM model could be rewritten in general as:(2) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\mathit{F}(\mathit{Z},t,\mathit{u}),\end{array}$$(2)

where $\mathit{u}=\{a,k,\mathit{E}\}$ is a set of model parameters.

The proposed model above explains the impact of an abrasive waterjet of radius a on the flat surface at the initial time moment, which changes its shape during the process and movement of the jet in the y-direction under the action of Etching rate function $\mathit{E}(x,y)$. The etching rate function, waterjets radius and change of the trenches profile accordingly to the depth and distance between the surface and the jets nozzle have influence on the intensity of the jet impact, and on the modification of the surface, which is described here as a function $\mathit{Z}(x,y,t)$. To model more specifically some of these processes, the exponential component and the gradient of the surface parametrization were introduced. Due to the dependency of Etching rate function $\mathit{E}$ on the properties of the workpiece material and machine parameters as pressure, velocity and abrasives mass flow, it should be calibrated using the experimental data.

Schematically AWJM process and jet footprint at some time moment $t\in [0,T]$ could be presented as on the Figure 1.

Figure 1. Schematic of the AWJM process and jet footprint.

2.2. Cost function and adjoint model

For a given set of parameters $\mathit{u}=\{a,k,\mathit{E}\}$ we formulate a minimization problem in order to find ${\mathit{u}}^{\ast }$ such that(3) $$\begin{array}{c}\hfill \mathit{J}({\mathit{u}}^{\ast })=\underset{\mathit{u}}{inf}\mathit{J}(\mathit{u}),\end{array}$$(3)

where $\mathit{J}(\mathit{u})$ is the cost function measuring the difference between the model solution and experimental observations(4) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\underset{\mathrm{\Omega }}{\int }\Vert \mathit{Z}(x,y,T)-{\mathit{Z}}_{\exp }{(x,y)\Vert }^2\mathrm{d}x\mathrm{d}y+\mathit{\alpha }{\Vert \mathit{u}-{\mathit{u}}_{\mathrm{b}}\Vert }^2,\end{array}$$(4)

under the constraint that $\mathit{Z}$ is the solution of our model with inputs $\mathit{u}\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}\times {\mathbb{L}}^2(\mathrm{\Omega },{\mathbb{R}}^{+})$.

Here ${\mathit{Z}}_{\exp }$ are the experimental measurements, ${\mathit{u}}_{\mathrm{b}}$ is an a priori estimation of the set of parameters $\mathit{u}$, and $\mathit{\alpha }>0$ is the Tikhonov regularization coefficient, that plays an important role in bringing the adjoint problem to a well-posed form and helping to identify suitable unknown parameters.

We have already introduced the solution of the parameters identification problem as a minimization problem of the mismatch between solution and experiment measurements, and in order to overcome the instability and reduce the inaccuracy, regularization theory suggests to augment the data mismatch by a regularization term, which is known as Tikhonov regularization [4].

In order to find an optimal ${\mathit{u}}^{\ast }$ one needs the gradient of the cost function $\mathrm{\nabla }\mathit{J}(\mathit{u})$ for the minimization process based on iterative gradient descent principles as quasi-Newton algorithms.

The minimization problem (3(3) $$\begin{array}{c}\hfill \mathit{J}({\mathit{u}}^{\ast })=\underset{\mathit{u}}{inf}\mathit{J}(\mathit{u}),\end{array}$$(3) ) could be regarded as an optimal control problem, and the approach based on Lagrangian multipliers could be used to get the solution of it [20,21]. The solution of the optimal control problem can therefore be regarded as the search of critical ‘points’ of the Lagrangian functional $\mathit{L}(\mathit{u},\mathit{Z},\mathit{P})$ associated to the constrained minimization, that we introduce as following:(5) $$\begin{array}{c}\hfill \mathit{L}(\mathit{u},\mathit{Z},\mathit{P})=\mathit{J}(\mathit{u})+\underset{0}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\int }\mathit{P}(t,x,y)\left(\frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}-\mathit{F}(\mathit{Z},t,\mathit{u})\right)\mathrm{d}x\mathrm{d}y\mathrm{d}t,\end{array}$$(5)

where $\mathit{P}$ is the Lagrangian multiplier associated to the constraint that $\mathit{Z}$ is a solution of (1(1) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\frac{\mathit{E}(x,y){\mathrm{e}}^{a\mathit{Z}}}{{\left(1+{\left|\mathrm{\nabla }\mathit{Z}\right|}^2\right)}^{k/2}}\text{in}\mathrm{\Omega }\times [0,T],\end{array}$$(1) ) from the same space ${\mathbb{H}}_0^1(\mathrm{\Omega })$.

We get the optimal system(6) $$\begin{array}{cc}\hfill \frac{\mathit{\partial }\mathit{L}}{\mathit{\partial }\mathit{P}}& =0,\hfill \end{array}$$(6) (7) $$\begin{array}{cc}\hfill \frac{\mathit{\partial }\mathit{L}}{\mathit{\partial }\mathit{Z}}& =0,\hfill \end{array}$$(7) (8) $$\begin{array}{cc}\hfill \frac{\mathit{\partial }\mathit{L}}{\mathit{\partial }\mathit{u}}& =0,\hfill \end{array}$$(8)

from where the statement (6(6) $$\begin{array}{cc}\hfill \frac{\mathit{\partial }\mathit{L}}{\mathit{\partial }\mathit{P}}& =0,\hfill \end{array}$$(6) ) says that $\frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\mathit{F}(\mathit{Z},t,\mathit{u})$, where $\mathit{Z}$ is solution of the direct model. We can obtain the adjoint equation, where $\mathit{P}$ should be also a solution from the expression (7(7) $$\begin{array}{cc}\hfill \frac{\mathit{\partial }\mathit{L}}{\mathit{\partial }\mathit{Z}}& =0,\hfill \end{array}$$(7) ). And from (8(8) $$\begin{array}{cc}\hfill \frac{\mathit{\partial }\mathit{L}}{\mathit{\partial }\mathit{u}}& =0,\hfill \end{array}$$(8) ) we have the gradient of the cost function from the initial adjoint state $\mathrm{\nabla }\mathit{J}(\mathit{u})=\mathit{P}(0)$, if $\mathit{Z}$ is a solution of the direct problem (2(2) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\mathit{F}(\mathit{Z},t,\mathit{u}),\end{array}$$(2) ), corresponding to the set of parameters $\mathit{u}$.

From Equations (5(5) $$\begin{array}{c}\hfill \mathit{L}(\mathit{u},\mathit{Z},\mathit{P})=\mathit{J}(\mathit{u})+\underset{0}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\int }\mathit{P}(t,x,y)\left(\frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}-\mathit{F}(\mathit{Z},t,\mathit{u})\right)\mathrm{d}x\mathrm{d}y\mathrm{d}t,\end{array}$$(5) ) and (7(7) $$\begin{array}{cc}\hfill \frac{\mathit{\partial }\mathit{L}}{\mathit{\partial }\mathit{Z}}& =0,\hfill \end{array}$$(7) ) we get by integration by parts:(9) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{P}}{\mathit{\partial }t}=-{\left(\frac{\mathit{\partial }\mathit{F}}{\mathit{\partial }\mathit{Z}}\right)}^{\mathrm{T}}\mathit{P}+\left(\mathit{Z}-{\mathit{Z}}_{\exp }\right)\mathit{\delta }(t-T),\end{array}$$(9) $$\begin{array}{c}\hfill \mathit{P}(x,y,T)=0,\\ \hfill \mathit{P}(x,y,t){\mid }_{\mathit{\partial }\mathrm{\Omega }}=0.\end{array}$$

To find out the optimal model parameters we have to consider the discrete system, not the continuous one. Indeed we need to minimize the discrete cost function and it requires the gradient of the discrete cost function. For that the discretized adjoint statement could be obtained from the discretized forms of the Lagrangian and its discrete derivatives, but more efficient way is to use AD software (i.e. TAPENADE) which bases on the principles explained above. It computes and provides the gradient of the discretized cost function which is further used to solve the minimization problem (3(3) $$\begin{array}{c}\hfill \mathit{J}({\mathit{u}}^{\ast })=\underset{\mathit{u}}{inf}\mathit{J}(\mathit{u}),\end{array}$$(3) ). Then all the optimal model parameters can be computed from only the trench profile, in other words from one experiment done on AWJ machine.

2.3. Minimization process

The most common minimization techniques are gradient descent algorithms which are iterative processes converging to the minimum. The minimization of the cost function could be realized in such a way, where the descent step is computed by quasi-Newton type algorithms. Particularly we used the L-BFGS based algorithms implemented in N2QN1 minimization package from ‘MODULOPT’ library [18].

Starting from the initial given position ${\mathit{u}}_0$ (which is defined from some background estimation ${\mathit{u}}_0={\mathit{u}}_{\mathrm{b}}$), one has to compute the cost function $\mathit{J}({\mathit{u}}_j)$ on each iteration and its gradient $\mathrm{\nabla }\mathit{J}({\mathit{u}}_j)$ to shift to the next discretized step and update the solution:$${\mathit{u}}_{j+1}={\mathit{u}}_j-{\mathit{\rho }}_j\mathrm{\nabla }\mathit{J}({\mathit{u}}_j),$$

where ${\mathit{\rho }}_j>0$ is a descent step of minimization, and $j>0$ is the iteration number.

This is a first order scheme which could be justified from the first order Taylor expansion of the cost function $\mathit{J}(\mathit{u})$ for any two iterations ${\mathit{u}}_j$ and ${\mathit{u}}_{j+1}$:$$\mathit{J}({\mathit{u}}_{j+1})=\mathit{J}({\mathit{u}}_j-{\mathit{\rho }}_j\mathrm{\nabla }\mathit{J}({\mathit{u}}_j))\simeq \mathit{J}({\mathit{u}}_j)-{\mathit{\rho }}_j{\Vert \mathrm{\nabla }\mathit{J}({\mathit{u}}_j)\Vert }^2\le \mathit{J}({\mathit{u}}_j).$$

The classical Newton type of minimization algorithms uses $-{\mathit{\rho }}_j{\mathit{H}}_j^{-1}\mathrm{\nabla }\mathit{J}({\mathit{u}}_j)$ as the direction of descent, where ${\mathit{H}}_j={\mathrm{\nabla }}^2\mathit{J}({\mathit{u}}_j)$ is the Hessian of the cost function $\mathit{J}$ at any iteration j. Due to the size of the problem one cannot compute the Hessian directly nor its inverse.

In quasi-Newton algorithms the Hessian does not need to be computed and inverted, but it can be replaced by a symmetric positive approximation $({\mathit{Q}}_j)$ to ${\mathit{H}}^{-1}$.

The L-BFGS [16] is a BFGS algorithm in the quasi-Newton’s family methods which uses an approximation of the inverse Hessian, but with limited amount of memory: it stores only m last values of ${\mathit{u}}_j$ and $\mathrm{\nabla }\mathit{J}({\mathit{u}}_j)$. The update formulas for L-BFGS algorithm are presented in [15].

3. Model parameters identification

3.1. Identification based on self-generated 2D profiles

From Equation (1(1) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\frac{\mathit{E}(x,y){\mathrm{e}}^{a\mathit{Z}}}{{\left(1+{\left|\mathrm{\nabla }\mathit{Z}\right|}^2\right)}^{k/2}}\text{in}\mathrm{\Omega }\times [0,T],\end{array}$$(1) ) we get the AWJM model corresponding to steady problems for non-moving jets, with initial and boundary conditions:(10) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\frac{\mathit{E}(x){\mathrm{e}}^{a\mathit{Z}}}{{\left(1+{{\mathit{Z}}_x}^2\right)}^{k/2}},\end{array}$$(10) $$\begin{array}{c}\hfill \mathit{Z}{\mid }_{t=0}=0,\\ \hfill \mathit{Z}{\mid }_{\mathit{\partial }{\mathrm{\Omega }}_1}=0.\end{array}$$

In this case we define a symmetric domain ${\mathrm{\Omega }}_1=\{x:x\in [-x_1;x_1]\}$, where $x_1$ always depends on the actual experimental parameters or measurements, due to the changes of the jets radius. Parameter identification problems usually bring many difficulties due to model errors and measurement noise. First we assume that measurement errors ${\mathit{\epsilon }}_{\exp }$, which come from the noisy measurements, are random variables with a Gaussian probability density function and a zero mean. To generate the noisy trench profile we introduce the following model:(11) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\frac{\mathit{E}(x){\mathrm{e}}^{a\mathit{Z}}}{{\left(1+{{\mathit{Z}}_x}^2\right)}^{k/2}}+\mathit{\lambda }{\mathit{\epsilon }}_{\exp },\end{array}$$(11)

where $\mathit{\lambda }$ is the factor corresponding to the percentage of the applied calibrated uncorrelated noise, ${\mathit{\epsilon }}_{\exp }\sim \mathcal{N}(0,1)$.

3.2. Discretization parameters for self-generated data

For the first numerical experiments, ‘pseudo-experimental’ surface profiles are generated with arbitrary values of model parameters $a=2$, $k=3$ and Etching rate function ${\mathit{E}}_0(x)$ on ${\mathrm{\Omega }}_1$, defined by setting $x_1=1$ mm. To discretize the domain ${\mathrm{\Omega }}_1$, a regular grid of 200 points with a step $\mathrm{\Delta }x=0.01$ mm is used and the time period is taken as unit $t\in [0,1]$ with $\mathrm{\Delta }t=\frac{\mathrm{\Delta }{x}^2}{4}$. To solve our minimization problem numerically the central difference method is used for the time integration.

The identification of both model parameters a and k is realized in the conditions of noise applied with different levels (10 and 30%) to the input surface. A Tikhonov regularization term is used with different values of the factor $\mathit{\alpha }$. This allows us to neutralize problems in identification process regarding the ill-posedness and to accelerate the minimization process. Using various values of $\mathit{\alpha }$, we show the importance of the regularization.

Results and comparisons are shown on Figure 2. We can notice that the background estimation of the unknown model parameters accelerates the minimization process, but the higher the level of noise in the input data the bigger errors in the identification. Even with a high level of noise, identification works properly and allows us to determine the unknowns being far enough from the first estimations. The choice of the regularization coefficient plays a very important role in the identification problem.

Figure 2. Results of the identification of model parameters a and k under the conditions of applied noise with levels of $10\%$ and $30\%$, using different Tikhonov regularization coefficients.

For the identification of AWJM model parameters a and k, one can see on Figure 2 that the use of right regularization factors strongly improves the speed and accuracy of the identification also with high levels of noise. In this particular case $\mathit{\alpha }=0.1$ gives the opportunity to identify initial model parameters within 8 and 15 iterations respectively, but with the value $\mathit{\alpha }=0.01$ it takes more time to approach the required precision. The use of too large or small values may prevent to identify the parameters at all.

The first part of the experiments only includes the identification process corresponding to the standard cost function (4(4) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\underset{\mathrm{\Omega }}{\int }\Vert \mathit{Z}(x,y,T)-{\mathit{Z}}_{\exp }{(x,y)\Vert }^2\mathrm{d}x\mathrm{d}y+\mathit{\alpha }{\Vert \mathit{u}-{\mathit{u}}_{\mathrm{b}}\Vert }^2,\end{array}$$(4) ), but in order to obtain a smooth solution ($E(x,y)\in {\mathbb{H}}^1(\mathrm{\Omega })$ instead of ${\mathbb{L}}^2(\mathrm{\Omega })$), we now change the regularization term to the one with the gradient of the Etching rate function:(12) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\underset{\mathrm{\Omega }}{\int }\Vert \mathit{Z}(x,T)-{\mathit{Z}}_{\exp }{(x)\Vert }^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2.\end{array}$$(12)

We also need to change the configuration of the minimization process as well, and note that in order to identify the function $\mathit{E}$, another value of Tikhonov regularization term $\mathit{\alpha }={10}^{-6}$ is now the closest to the optimal one. It was obtained by L-curve method which was first applied by Lawson and Hanson [22] and more recently by Hensen and O’Leary [23].

The goal of the next simulations is to correctly identify the Etching rate function $\mathit{E}$. It was used to generate the surface profile in the perfect conditions without any noises or errors accounted in the model. The other values of the model parameters $a=2,k=3$ were fixed. Results of the identification are presented on Figure 3. The point of identifying the Etching rate function is to focus on the most influential term of AWJM model, that increases the number of unknowns by several orders.

Figure 3. Results of the identification of Etching rate function E starting from different initial conditions.

The results shown on Figure 3 confirm the possibility to identify the Etching rate function in the ideal situation even with a wrong first estimation. Here the total number of unknowns is 200 which is the number of mesh points. In case of a really bad first estimation (Figure 3(b)) the minimization process takes more or less the same time and stops after 19 iterations by reaching the defined accuracy of ${10}^{-4}$, when in the first case, shown on Figure 3(a), it takes 16 iterations.

More interesting is the case when the Etching rate function is not so smooth and simple, or even more is unknown, and the trench profile has some specialties caused ratherish by measurement errors.

3.3. Real experiment measurements

Unlike the previous cases, the trench profile is now the only input for our problem. We present the results of identification of the Etching rate function which should be used in the direct AWJM model to be able to reproduce the required trench profile. These measurements were obtained from experiments done in collaboration with STEEP project partners. This data corresponds to experiments with a jet feed speed of 2000 mm/min, and nozzle diameter of AWJ machine of 0.5 mm. We use this value as a background estimation of model parameter a. The number of mesh points in this case equals 228 with the step $\mathrm{\Delta }x=0.0048$ mm and ${\mathrm{\Omega }}_1$ is defined by setting $x_1=0.55$ mm. The minimization process is realized using N2QN1 minimizer from the ‘MODULOPT library’ with the regularization factor $\mathit{\alpha }={10}^{-5}$ that was obtained numerically by L-curve method. In lack of any knowledge about the Etching rate function we started the identification from the zero assumption ${\mathit{E}}_0=\mathbf{0}$.

Figure 4. Behavior of the cost function and its gradient in the minimization process, corresponding to the real measurement data.

Figure 4(a) and (b) shows respectively the behavior of the cost function and of the norm of its gradient during the minimization process. It shows that the minimizer N2QN1 works well. Using 15 minimization iterations, the Etching rate function is identified with the accuracy of $4\times {10}^{-3}$.

Figure 5. Results of the identification of Etching rate function and corresponding trench profile in comparison with experimental measurements.

The identified Etching rate function from the experimental data allows us to reproduce numerically the trench profile with an error in terms of $L^2$ norm smaller than 3% (Figure 5). The used mathematical model (1(1) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\frac{\mathit{E}(x,y){\mathrm{e}}^{a\mathit{Z}}}{{\left(1+{\left|\mathrm{\nabla }\mathit{Z}\right|}^2\right)}^{k/2}}\text{in}\mathrm{\Omega }\times [0,T],\end{array}$$(1) ) has some limitations itself and cannot cover some effects on the edges of the trench which in practice is explained by redeposition of the surface material. These secondary effects appearing as the result of high power of the waterjet impact should be studied separately and are not considered in our work.

4. Numerical results of sensitivity study

One of the most important things in the identification problems is the sensitivity study which could give us the opportunity to provide correct AWJM model parameters to reproduce required shape of the trench. The point of this section is to study the capacity of the approach, and to observe the possibilities to identify model parameters even with high level of measurement errors which are always present in the provided input data. Based on AWJM model (11(11) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\frac{\mathit{E}(x){\mathrm{e}}^{a\mathit{Z}}}{{\left(1+{{\mathit{Z}}_x}^2\right)}^{k/2}}+\mathit{\lambda }{\mathit{\epsilon }}_{\exp },\end{array}$$(11) ) we generate different trench profiles with predefined values of $\mathit{E}={\mathit{E}}_0,a=a_0,k=k_0$ respectively to different levels of uncorrelated noise up to 40%. We use these noisy profiles as the only input of the method and we try to identify the Etching rate function that could satisfy the requirements. As this part of the study is based again on self-generated input data, we use the same numerical parameters and initial assumptions as in Section 3.3. We also assume that initial Etching rate function ${\mathit{E}}_0$ has symmetrical gaps corresponding to the edges of the trench.

Tikhonov regularization has a very strong influence. It allows us to vary the identification results between accuracy and smoothness. We can then choose the optimal values that can be acceptable and suitable in the real experiments and lead to reconstruction of the surface profile with high precision at the same time.

The identification of the model parameter $\mathit{E}$ is based only on the experimental measurements which include the measurement errors. We simulate such data (Figure 7) by adding a Gaussian white noise to the initial surface profile, generated by use of the Etching rate function named ‘Original’ on Figures 6 and 8. The goal of this study is to identify the unknown $\mathit{E}$ the use of which in the AWJM model (10(10) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{Z}}{\mathit{\partial }t}=\frac{\mathit{E}(x){\mathrm{e}}^{a\mathit{Z}}}{{\left(1+{{\mathit{Z}}_x}^2\right)}^{k/2}},\end{array}$$(10) ) will form the closest trench to the initial one, named ‘Target’ on Figures 6, 8, 9.

In order to obtain a smooth solution, we again use the cost function (12(12) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\underset{\mathrm{\Omega }}{\int }\Vert \mathit{Z}(x,T)-{\mathit{Z}}_{\exp }{(x)\Vert }^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2.\end{array}$$(12) ) with the regularization term on the gradient. This will ensure the absence of high oscillations (note that the regularization coefficient may be newly estimated through a L-curve method).

Figure 6. Results of numerical experiments in the identification the Etching rate functions for AWJM model and reconstruction of the surface shapes based on single trench profiles for measurement noises of 5, 15 and 30% respectively.

Figure 6(a)–(c) demonstrates the results of identification of the Etching rate function $\mathit{E}$ in case of noise with levels of 5, 15 and 30%, and comparison of reproduced trenches with original profiles. One could also notice that even with considering a very high level of noise in the measurements, it is still possible to identify the model parameter $\mathit{E}$. Taking these founded values, we can then model and forecast the shape of the surface, even if the form and view of the founded function $\mathit{E}$ is not perfectly matching.

Assume now that we have two or more different measurements of exactly the same experiment. The difference between them is only in random noises that are present in measurements and influence the available input (Figure 7).

Figure 7. Simulation of the 2 independent trench measurements with the level of noise = 15%.

In this case our cost function transforms to(13) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\frac{1}{2}\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-{\mathit{Z}}_{{\exp }_1}(x)∥}^2\mathrm{d}x+\frac{1}{2}\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-{\mathit{Z}}_{{\exp }_2}(x)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2,\end{array}$$(13)

and it leads to the following results for the identification process (Figure 8).

Figure 8. Results of numerical experiments in the identification the Etching rate functions for AWJM model and reconstruction of the surface shapes based on two different measurements of one trench profile in case the level of the measurement errors is 15%.

As it is possible to see on Figure 8 the obtained results are much more precise than in the previous case (Figure 6(b), only one trench profile), as the form of the identified function E becomes more symmetric and smooth. Another improvement is the symmetry of the solution obtained from very unclear measurements and increased accuracy of the trench reconstruction.

The use of more than one measurement of the only experimental data in theory will provide some kind of averaging of the trench profiles that should lead to smoothing and give higher opportunity to reconstruct the surface more precisely. Based on that assumption we introduce the superposition of the same two different measurements in the cost function as in the previous numerical case, that transforms the cost function as follows:(14) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-\left(\frac{{\mathit{Z}}_{{\exp }_1}(x)+{\mathit{Z}}_{{\exp }_2}(x)}{2}\right)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2.\end{array}$$(14)

Numerical results of this case (Figure 9(a)) show that the accuracy of the reconstruction of the trench profile is also higher in comparison with the case where only two trenches are used independently (Figure 8).

Figure 9. Results of numerical experiments in the identification the Etching rate functions for AWJM model and reconstruction of the surface shapes based on different variants of cost functions in case the level of the measurement errors is 15%.

Considering some improvement above, we now include one more additional experimental measurement, which is involved into the superposition of trenches to form the input for minimization process. First we introduce it independently from other measurements, so that the corresponding cost function is:(15) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\sum _{i=1}^3\frac{1}{3}\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-{\mathit{Z}}_{{\exp }_i}(x)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2,\end{array}$$(15)

Numerical results are displayed on the Figure 9(b), where one could see that mismatch between solution and ‘Target’ is still very low. Note that identified Etching rate function already keeps the smooth symmetric form.

Table 1. Comparison of the accuracy in the prediction of the trench profile, corresponding to different cases of the cost functions and different levels of applied noise.

Trenches	1%	2%	5%	10%	15%	20%	30%	40%
Single	$2.37\times {10}^{-3}$	$2.24\times {10}^{-3}$	$1.87\times {10}^{-2}$	$2.38\times {10}^{-2}$	$6.08\times {10}^{-2}$	$6.14\times {10}^{-2}$	$8.78\times {10}^{-2}$	0.272735
2 independent	$2.23\times {10}^{-3}$	$4.66\times {10}^{-3}$	$1.15\times {10}^{-2}$	$2.30\times {10}^{-2}$	$1.84\times {10}^{-2}$	$9.08\times {10}^{-2}$	$9.38\times {10}^{-2}$	$6.46\times {10}^{-2}$
2 superposed	$3.43\times {10}^{-3}$	$1.12\times {10}^{-2}$	$1.19\times {10}^{-2}$	$2.56\times {10}^{-2}$	$1.79\times {10}^{-2}$	$8.64\times {10}^{-2}$	$9.06\times {10}^{-2}$	$6.43\times {10}^{-2}$
3 independent	$1.68\times {10}^{-3}$	$3.66\times {10}^{-3}$	$1.18\times {10}^{-2}$	$1.84\times {10}^{-2}$	$1.92\times {10}^{-2}$	$8.94\times {10}^{-2}$	$9.53\times {10}^{-2}$	$6.50\times {10}^{-2}$
3 superposed	$1.68\times {10}^{-3}$	$3.65\times {10}^{-3}$	$1.18\times {10}^{-2}$	$2.08\times {10}^{-2}$	$1.92\times {10}^{-2}$	$8.97\times {10}^{-2}$	$9.86\times {10}^{-2}$	$6.39\times {10}^{-2}$

Finally all the three measurements used before can also be combined together into the formation of the input for minimization process. It means that next results, presented on Figure 9(c) rely on the superposition of the three different measurements of the same trench, leading to the following cost function:(16) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-\left(\frac{{\mathit{Z}}_{{\exp }_1}(x)+{\mathit{Z}}_{{\exp }_2}(x)+{\mathit{Z}}_{{\exp }_3}(x)}{3}\right)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2.\end{array}$$(16)

The difference between using two and three measurements is not very impressive due to random nature of the noise applied to the input, and could be strongly increased by involving hundreds of experimental observations to reduce the influence of the errors. But it is already possible to claim that use of three trench measurements instead of two gives more close and correct shape of the trench, and allows us to identify a more acceptable Etching rate function in the manufacturing.

Figure 10. Accuracy of matching the obtained solutions to original profiles for different cost functions considering different levels of measurement errors.

Figure 10 gives a good overview of the improvement in the identification achieved by introducing different approaches to define the cost function. Keeping the acceptable level of errors (less than 10%) in matching between experimental measurement and obtained solution, introduction of more complex cost functions (13(13) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\frac{1}{2}\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-{\mathit{Z}}_{{\exp }_1}(x)∥}^2\mathrm{d}x+\frac{1}{2}\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-{\mathit{Z}}_{{\exp }_2}(x)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2,\end{array}$$(13) ) and (15(15) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\sum _{i=1}^3\frac{1}{3}\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-{\mathit{Z}}_{{\exp }_i}(x)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2,\end{array}$$(15) ) gives us the opportunity to identify more applicable and convenient AWJM model parameters. In general the more trenches measurements are available the higher accuracy in the surface prediction can be achieved due to elimination of the measurement errors by averaging the input data.

Cost functions (13(13) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\frac{1}{2}\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-{\mathit{Z}}_{{\exp }_1}(x)∥}^2\mathrm{d}x+\frac{1}{2}\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-{\mathit{Z}}_{{\exp }_2}(x)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2,\end{array}$$(13) ) and (14(14) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-\left(\frac{{\mathit{Z}}_{{\exp }_1}(x)+{\mathit{Z}}_{{\exp }_2}(x)}{2}\right)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2.\end{array}$$(14) ) are not identical, but they theoretically have the same gradient. However numerical implementation gives different results and flexibility to find more suitable realization for each particular problem, what has to be taken into account. The same situation is for the pair of cost functions (15(15) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\sum _{i=1}^3\frac{1}{3}\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-{\mathit{Z}}_{{\exp }_i}(x)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2,\end{array}$$(15) ) and (16(16) $$\begin{array}{c}\hfill \mathit{J}(\mathit{u})=\underset{\mathrm{\Omega }}{\int }{∥\mathit{Z}(x,T)-\left(\frac{{\mathit{Z}}_{{\exp }_1}(x)+{\mathit{Z}}_{{\exp }_2}(x)+{\mathit{Z}}_{{\exp }_3}(x)}{3}\right)∥}^2\mathrm{d}x+\mathit{\alpha }{\Vert \mathrm{\nabla }\mathit{E}\Vert }^2.\end{array}$$(16) ). More detailed results of comparison between all proposed approaches are presented in Table 1.

One could notice that in most of the results the use of several trenches instead of only one can improve the accuracy in the parameters identification, leading to reducing the errors in the prediction of the surface profile up to 20% in cases of low level of noise. In cases of very noisy input data (see columns related to noise higher than 20%) the use of several measurements plays strong role in decreasing the mismatch in reconstruction of the trenches. This effect could be explained by nature of the distribution of the applied noise. Certainly, it should be noted that sometimes only one measurement is available, and it might be enough to obtain the model parameters required to reconstruct the profile. Moreover, superposition of the trenches involved in the identification process which can be interpreted as an average of the input experiments usually also leads to higher accuracy, but mostly on the ‘long distance’ – the quantity of available measurements.

5. Conclusion

The identification of model parameters, especially from the noisy data, is a very challenging problem in the real production. In this paper, we have presented the application of inverse problems theory, based on minimization problems, to Abrasive WaterJet machining. The identification of optimal model parameters indeed gives a chance to model and predict the trench profile for AWJ machining. We illustrated how an even small level of errors or noise can highly influence the identification results and lead to significant errors in the determination. Also we explained how the high level of noise could completely change the identification process and why it is necessary to keep the regularization terms to get a more precise structure and shape of the reconstructed trench.

We presented a way to estimate the required surface profile in the lack of knowledge of exact AWJ model parameters. Proposed AWJM model coupled with measurement errors which were taken into account can precisely model the surface shape by identifying optimal values of the model parameters even with very poor and inaccurate experimental observations.

In this paper we considered measurement errors that were included in the mathematical model, which we tried to compensate by adding Tikhonov regularization terms. Development of this study in the case of non-steady AWJM process in 3D, including the variation of the jet feed speed, leads to complication of the model and numerical realization of the identification process. All these aspects may be studied more precisely in a future work, along with different regularization techniques that might help keeping the identification problem well-posed.

Acknowledgements

The authors thank Mr. Pablo Lozano Torrubia from the University of Nottingham for his help in the experimental work.

Additional information

Funding

This work was supported by the EU-FP7-ITN [grant number 316560] for the works presented as a part of the STEEP ITN project.

Notes

No potential conflict of interest was reported by the authors.