1. Introduction

2. Framework for Feedrate Optimization with Uncertainty-aware Digital Twin

Figure 2. (a) Need for a tolerance range due to violation of error tolerance in the presence of uncertainties, (b) Proposed method of feedrate optimization with desired tolerance stringency using quantified uncertainty.

Figure 2. (a) Need for a tolerance range due to violation of error tolerance in the presence of uncertainties, (b) Proposed method of feedrate optimization with desired tolerance stringency using quantified uncertainty.

3. Contour error prediction with a real-time uncertainty-aware digital twin

3.2. Prediction and uncertainty quantification of contour error using uncertainty-aware digital twin

The accuracy of the predictions of the physics-based and data-driven servo models in Section 3.1 can be improved by incorporating the known uncertainty from the physics-based models. To do so, a digital twin that uses physics-informed data-driven servo model is proposed, which is updated on-the-fly via Bayesian approach that is capable of considering uncertainty in the feature vector ${\mathsf{\psi }}_t$. A flowchart of the proposed intelligent feedrate optimization using the uncertainty-aware digital twin is given in Figure 5. The key difference between Figure 3 and Figure 5 is that known uncertainty is included in ${\widehat{\mathit{H}}}_x$, and unknown uncertainty embedded in ${\mathit{x}}^j$ is learned using the data-driven model and used to predict ${\widehat{\mathit{x}}}^j$ with quantified uncertainty.

Figure 5. Flowchart of intelligent feedrate optimization using an uncertainty-aware digital twin (with physics-based and data-driven servo models, y-axis omitted).

Figure 5. Flowchart of intelligent feedrate optimization using an uncertainty-aware digital twin (with physics-based and data-driven servo models, y-axis omitted).

Each element $\tilde{x}\left(t\right)$ of the digital twin’s prediction of the output position ${\tilde{\mathit{x}}}^j$ is modeled as

$$\tilde{x}\left(t\right)={\mathsf{\beta }}^{\mathrm{T}}{\mathsf{\psi }}_t+ϵ$$

(7)

where ${\mathsf{\psi }}_t\sim N({\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})$ is the feature vector defined as an uncorrelated Gaussian random variable with mean ${\mathsf{\mu }}_{{\psi }_t}$ and variance ${\mathsf{\Sigma }}_{{\psi }_t}$ derived from the uncertainty distribution of ${\mathsf{\psi }}_t$ previously defined in Eq. (2). $\mathsf{\beta }\sim N({\mathsf{\mu }}_{\beta },{\mathsf{\Sigma }}_{\beta })$ is the weight vector defined as a correlated Gaussian random variable learned via Bayesian linear regression, and $ϵ\sim N(0,{\sigma }_{ϵ}^2)$ is the unobserved Gaussian noise. Unlike the deterministic point-estimation of ${\tilde{\mathit{x}}}^j$ from Section 3.1, a distribution ${\tilde{\mathit{x}}}^j\sim N({\mathsf{\mu }}_{{\tilde{x}}^j},{\mathsf{\Sigma }}_{{\tilde{x}}^j})$ is estimated in this section based on the uncertainties of the features and the weights.

This paper proposes that the known uncertainties of the physics-based models are embedded into the feature vector ${\mathsf{\psi }}_t$ in Eq. (7) to enable efficient training of $\mathsf{\beta }$. To do so, a set of $N_H$ stable physics-based models ${\left\{{\widehat{\mathit{H}}}_x^i\right\}}_{i=1}^{N_H}$ is obtained, where each model ${\widehat{\mathit{H}}}_x^i$ for $i\in \{1,2,...,N_H\}$ is identified in the form of the complex-valued frequency response function (FRF) of the physical system ${\mathit{H}}_x$ at discrete frequencies ${\omega }_k$ via experiments as

$${\widehat{\mathit{H}}}_x^i\left({\omega }_k\right)=a^i\left({\omega }_k\right)+b^i\left({\omega }_k\right)j$$

(8)

where ${\omega }_k=k\mathsf{\Delta }\omega$, of which $\mathsf{\Delta }\omega$ is the increment of frequencies, $k\in \{1,2,...,\frac{N_{\omega }}{2}-1\}$ where $N_{\omega }$ is the number of discrete negative and positive frequencies at which FRF is identified, and j is the unit imaginary number (which should not be confused with the batch index j used as a superscript elsewhere).

Then, the uncertainty in ${\widehat{\mathit{H}}}_x$ is propagated to the finite impulse response ${\widehat{\mathit{h}}}_x$ as follows. The discrete sets ${\left\{{\widehat{\mathit{H}}}_x^i\left({\omega }_k\right)\right\}}_{i=1}^{N_H}$ for each k will introduce discrete sets of their real and imaginary coefficients, namely ${\left\{a^i\left({\omega }_k\right)\right\}}_{i=1}^{N_H}$ and ${\left\{b^i\left({\omega }_k\right)\right\}}_{i=1}^{N_H}$. For computational efficiency, it is assumed that ${\left\{a^i\left({\omega }_k\right)\right\}}_{i=1}^{N_H}$ and ${\left\{b^i\left({\omega }_k\right)\right\}}_{i=1}^{N_H}$ are sampled from Gaussian distributions of $a\left({\omega }_k\right)$ and $b\left({\omega }_k\right)$, of which 99.73%, i.e., the 3-sigma range, lie within the minimum and maximum of the identified discrete sets. Then, $a\left({\omega }_k\right)\sim N({\mu }_{a\left({\omega }_k\right)},{\sigma }_{a\left({\omega }_k\right)}^2)$ can be approximated as

$$\begin{array}{cc}\hfill {\mu }_{a\left({\omega }_k\right)}& =\frac{max\left\{a\left({\omega }_k\right)\right\}+min\left\{a\left({\omega }_k\right)\right\}}{2}\hfill \\ \hfill {\sigma }_{a\left({\omega }_k\right)}& =\frac{max\left\{a\left({\omega }_k\right)\right\}-{\mu }_{a\left({\omega }_k\right)}}{3}\hfill \end{array}$$

(9)

and $b\left({\omega }_k\right)\sim N({\mu }_{b\left({\omega }_k\right)},{\sigma }_{b\left({\omega }_k\right)}^2)$ can be approximated using the same procedure as Eq. (9).

Then, the impulse response ${\widehat{\mathit{h}}}_x$ (i.e., ${\widehat{h}}_x\left[n\right]$ for $n\in \{1,2,...,N_h\}$) of the physics-based model with sampling time $T_s$ can be formulated using discrete inverse Fourier transform as

$$\widehat{h}\left[n\right]=\frac{1}{N_{\omega }}\sum _{k=\frac{-N_{\omega }}{2}}^{\frac{N_{\omega }}{2}-1}(a\left({\omega }_k\right)+b\left({\omega }_k\right)j)e^{\frac{2\pi j(n-1)(k-1)}{N_{\omega }}}$$

(10)

where $N_h=\frac{1}{T_s\mathsf{\Delta }\omega }$ is the maximum estimable length of the impulse response. Then, due to the linearity of Eq. (10), $\widehat{h}\left[n\right]$ follows a Gaussian distribution $\widehat{h}\left[n\right]\sim N({\mu }_{\widehat{h}\left[n\right]},{\sigma }_{\widehat{h}\left[n\right]}^2)$ as

$$\begin{array}{cc}\hfill & {\mu }_{\widehat{h}\left[n\right]}=\frac{1}{N_{\omega }}\sum _{k=1}^{N_{\omega }}({\mu }_{a\left({\omega }_k\right)}+{\mu }_{b\left({\omega }_k\right)}j)e^{\frac{2\pi j(n-1)(k-1)}{N_{\omega }}}\hfill \\ \hfill & {\sigma }_{\widehat{h}\left[n\right]}^2=\frac{1}{N_{\omega }}\sum _{k=1}^{N_{\omega }}({\sigma }_{a\left({\omega }_k\right)}^2-{\sigma }_{b\left({\omega }_k\right)}^2)e^{\frac{2\pi j(n-1)(k-1)}{N_{\omega }}}\hfill \end{array}$$

(11)

Next, the uncertainty in $\widehat{h}\left[n\right]$ can be propagated to $\widehat{x}\left(t\right)\sim N({\mu }_{\widehat{x}\left(t\right)},{\sigma }_{\widehat{x}\left(t\right)}^2)$ as

$$\begin{array}{cc}\hfill & {\mu }_{\widehat{x}\left(t\right)}=\sum _{i=1}^{N_h}{\mu }_{\widehat{h}\left[i\right]}{x}_d(t-i{T}_s)\hfill \\ \hfill & {\sigma }_{\widehat{x}\left(t\right)}^2=\sum _{i=1}^{N_h}{\sigma }_{\widehat{h}\left[i\right]}^2{x}_d{(t-i{T}_s)}^2\hfill \end{array}$$

(12)

Finally, the feature vector follows an uncorrelated multivariate normal distribution ${\mathsf{\psi }}_t\sim N({\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})$, where

(13)

The prior over $\mathsf{\beta }\sim N({\mathsf{\mu }}_{\beta }$, ${\mathsf{\Sigma }}_{\beta })$ is initialized using the nominal physics-based model ${\mu }_{{\widehat{h}}_x}$ from Eq. (11), as described in Appendix A. Now given $\mathsf{\beta }\sim N({\mathsf{\mu }}_{\beta },{\mathsf{\Sigma }}_{\beta })$ and the feature uncertainty ${\mathsf{\psi }}_t\sim N({\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})$ from Eq. (13), we propose to take a Bayesian route to obtain the posterior over $\mathit{\beta }$ as more data is collected in real time.

Mathematically, the posterior is given as:

$$\begin{array}{cc}\hfill p\left(\mathsf{\beta }\right|x\left(t\right),{\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})& \propto p\left(x\left(t\right)\right|\mathsf{\beta },{\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})·p\left(\mathsf{\beta }\right)\hfill \\ \hfill & \propto \int p(x\left(t\right),{\mathsf{\psi }}_t|\mathsf{\beta },{\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})d{\mathsf{\psi }}_t·p\left(\mathsf{\beta }\right)\hfill \\ \hfill & \propto \int p\left(x\left(t\right)\right|{\mathsf{\psi }}_t,\mathsf{\beta })p\left({\mathsf{\psi }}_t\right|{\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})d{\mathsf{\psi }}_t\hfill \\ \hfill & ·p\left(\mathsf{\beta }\right)\hfill \\ \hfill & \propto p\left(x\right(t\left)\right|\mathsf{\beta })·p(\mathsf{\beta })\hfill \end{array}$$

(14)

where $x\left(t\right)$ is the actual position on the x-axis at timestep t defined in Section 3.1. Here, $p\left(x\right(t\left)\right|\mathsf{\beta })$ and $p\left(\mathsf{\beta }\right)$ are given as as

$$\begin{array}{cc}\hfill & p\left(x\left(t\right)\right|\mathsf{\beta })=\exp (-\frac{1}{2}{(x\left(t\right)-{\mathsf{\beta }}^{\mathrm{T}}{\mathsf{\mu }}_{{\psi }_t})}^{\mathrm{T}}{({\mathsf{\beta }}^{\mathrm{T}}{\mathsf{\Sigma }}_{{\psi }_t}\mathsf{\beta }+{\sigma }_{ϵ}^2)}^{-1}\hfill \\ \hfill & (x\left(t\right)-{\mathsf{\beta }}^{\mathrm{T}}{\mathsf{\mu }}_{{\psi }_t}))\hfill \\ \hfill & p\left(\mathsf{\beta }\right)=\exp \left(-\frac{1}{2}{(\mathsf{\beta }-{\mathsf{\mu }}_{\beta })}^{\mathrm{T}}{\mathsf{\Sigma }}_{\beta }^{-1}(\mathsf{\beta }-{\mu }_{\mathsf{\beta }})\right)\hfill \end{array}$$

(15)

As shown above, a central feature of our approach is that it accounts for input uncertainty through ${\mathsf{\psi }}_t\sim N({\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})$. Unfortunately, the price to pay is the lack of a closed form solution due to the quadratic term in the variance. While monte carlo (MC) sampling approach such as Markov chain MC (i.e. MCMC) or Hamiltonian MC can be used, our goal is to obtain the posterior on the fly to enable real-time control.

To this end, we take a Gaussian Laplacian approximation [38] where we write log of $p\left(x\right(t\left)\right|\mathsf{\beta })$ as a quadratic function of $\mathsf{\beta }$. This is shown below:

$$\begin{array}{cc}\hfill lnp\left(x\right(t\left)\right|\mathsf{\beta })& \approx \underset{:=c_0}{\underbrace{lnp\left(x\left(t\right)\right|\mathsf{\beta }){|}_{\mathsf{\beta }=\overline{\mathsf{\beta }}}}}\hfill \\ \hfill & +\underset{:={\mathit{c}}_1}{\underbrace{\frac{\mathrm{d}lnp\left(x\right(t\left)\right|\mathsf{\beta })}{\mathrm{d}\mathsf{\beta }}{|}_{\mathsf{\beta }=\overline{\mathsf{\beta }}}^{\mathrm{T}}}}(\mathsf{\beta }-\overline{\mathsf{\beta }})\hfill \\ \hfill & +\frac{1}{2}{(\mathsf{\beta }-\overline{\mathsf{\beta }})}^{\mathrm{T}}\underset{:={\mathit{C}}_2}{\underbrace{\frac{{\mathrm{d}}^{2}lnp\left(x\left(t\right)\right|\mathsf{\beta })}{\mathrm{d}{\mathsf{\beta }}^2}{|}_{\mathsf{\beta }=\overline{\mathsf{\beta }}}}}(\mathsf{\beta }-\overline{\mathsf{\beta }})\hfill \end{array}$$

(16)

where the local point $\overline{\mathsf{\beta }}$ is ${\mu }_{\mathsf{\beta }}$ from the previous batch, i.e., the prior mean. A detailed derivation of the coefficients $c_0,{\mathit{c}}_1$ and ${\mathit{C}}_2$ in Eq. (16) is provided in Appendix B. Given that $p\left(x\left(t\right)\right|\mathit{\beta })\sim N({\mathsf{\mu }}_{\beta },{\mathsf{\Sigma }}_{\beta })$, now Eq. (14) can re-written in a closed-form solution as

$$\begin{array}{cc}\hfill p\left(\mathsf{\beta }\right|x\left(t\right),{\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})=\exp (& {\mathsf{\beta }}^{\mathrm{T}}\left(\frac{1}{2}{\mathit{C}}_2-\frac{1}{2}{\mathsf{\Sigma }}_{\beta }^{-1}\right)\mathsf{\beta }\hfill \\ \hfill & +\left({\mathit{c}}_1^{\mathrm{T}}+{\mathsf{\mu }}_{\beta }^{\mathrm{T}}{\mathsf{\Sigma }}_{\beta }^{-1}\right)\mathsf{\beta }+\mathrm{constant})\hfill \end{array}$$

(17)

Thus, the updated variance and mean of the Gaussian posterior $p\left(\mathsf{\beta }\right|x\left(t\right),{\mathsf{\mu }}_{{\psi }_t},{\mathsf{\Sigma }}_{{\psi }_t})$ are

(18)

where the posteriors ${\mathsf{\mu }}_{\beta }$ and ${\mathsf{\Sigma }}_{\beta }$ are used as priors for the next batch.

Note that, if the feature ${\mathsf{\psi }}_t$ is assumed to be deterministic (i.e., zero feature uncertainty), no known uncertainties will be included in the Bayesian regression, of which the case will be compared with the proposed method using normally distributed ${\mathsf{\psi }}_t$ in the following Section 5. In the case of deterministic ${\mathsf{\psi }}_t$, the posterior $\mathsf{\beta }\sim N({\mathsf{\mu }}_{\beta },{\mathsf{\Sigma }}_{\beta })$ can be estimated as a closed-form solution using Bayes rule as

$$\begin{array}{cc}\hfill & p\left(\mathsf{\beta }\right|{\mathsf{\mu }}_{\beta },{\mathsf{\Sigma }}_{\beta })\hfill \\ \hfill & =\exp (-\frac{1}{2}{\mathsf{\beta }}^{\mathrm{T}}({\sigma }_{ϵ}^{-2}{\mathsf{\psi }}_t^{\mathrm{T}}{\mathsf{\psi }}_t+{\mathsf{\Sigma }}_{\beta }^{-1})\mathsf{\beta }\hfill \\ \hfill & +{({\sigma }_{ϵ}^{-2}{\mathsf{\psi }}_t^{\mathrm{T}}x\left(t\right)+{\mathsf{\Sigma }}_{\beta }^{-1}{\mathsf{\mu }}_{\beta })}^{\mathrm{T}}\mathsf{\beta })\hfill \\ \hfill & \therefore {\mathsf{\Sigma }}_{\beta }\leftarrow {({\sigma }_{ϵ}^{-2}{\mathsf{\psi }}_t^{\mathrm{T}}{\mathsf{\psi }}_t+{\mathsf{\Sigma }}_{\beta }^{-1})}^{-1}\hfill \\ \hfill & {\mathsf{\mu }}_{\beta }\leftarrow {\mathsf{\Sigma }}_{\beta }({\sigma }_{ϵ}^{-2}{\mathsf{\psi }}_t^{\mathrm{T}}x\left(t\right)+{\mathsf{\Sigma }}_{\beta }^{-1}{\mathsf{\mu }}_{\beta })\hfill \end{array}$$

(19)

Lastly, using the trained weight distribution from Eq. (18), the posterior predictive distribution $\tilde{x}\left(t\right)\sim N({\mu }_{\tilde{x}\left(t\right)},{\sigma }_{\tilde{x}\left(t\right)}^2)$ can be written as

$$\begin{array}{cc}\hfill & {\mu }_{\tilde{x}\left(t\right)}={\mathsf{\mu }}_{\beta }^{\mathrm{T}}{\mathsf{\psi }}_t\hfill \\ \hfill & {\sigma }_{\tilde{x}\left(t\right)}^2={\mathsf{\psi }}_t^{\mathrm{T}}{\mathsf{\Sigma }}_{\beta }{\mathsf{\psi }}_t+{\sigma }_{ϵ}^2\hfill \end{array}$$

(20)

The same procedures are applied to learn ${\mathsf{\beta }}_y\sim N({\mathsf{\mu }}_{{\beta }_y},{\mathsf{\Sigma }}_{{\beta }_y})$ and predict $\tilde{y}\left(t\right)\sim N({\mu }_{\tilde{y}\left(t\right)},{\sigma }_{\tilde{y}\left(t\right)}^2)$ based on y-axis feature vector ${\mathsf{\psi }}_{ty}$ and unobserved noise ${ϵ}_y$. Using Eq. (6), the contour error distribution $\tilde{\epsilon }\left(t\right)\sim N({\mu }_{\tilde{\epsilon }\left(t\right)},{\sigma }_{\tilde{\epsilon }\left(t\right)}^2)$ can be predicted as

$$\begin{array}{cc}\hfill & {\mu }_{\tilde{\epsilon }\left(t\right)}=-sin\left(\mathsf{\theta }\left(t\right)\right)(x_d\left(t\right)-{\mathsf{\mu }}_{\beta }^{\mathrm{T}}{\mathsf{\psi }}_t)\hfill \\ \hfill & +cos\left(\mathsf{\theta }\left(t\right)\right)(y_d\left(t\right)-{\mathsf{\mu }}_{{\beta }_y}^{\mathrm{T}}{\mathsf{\psi }}_{ty})\hfill \\ \hfill & {\sigma }_{\tilde{\epsilon }\left(t\right)}^2=sin{\left(\mathsf{\theta }\left(t\right)\right)}^2({\mathsf{\psi }}_t^{\mathrm{T}}{\mathsf{\Sigma }}_{\beta }{\mathsf{\psi }}_t+{\sigma }_{ϵ}^2)\hfill \\ \hfill & +cos{\left(\mathsf{\theta }\left(t\right)\right)}^2({\mathsf{\psi }}_{ty}^{\mathrm{T}}{\mathsf{\Sigma }}_{{\beta }_y}{\mathsf{\psi }}_{ty}+{\sigma }_{{ϵ}_y}^2)\hfill \end{array}$$

(21)

4. Methodology for intelligent feedrate optimization with contour error constraints

5. Numerical validation of the intelligent feedrate optimization using uncertainty-aware digital twin

Figure 12. Feedrate, acceleration, contour error, tolerance violation, and prediction error using conservative (Cons.) physics-based (PB), data-driven (DD) and proposed (PIDD) methods with $\eta$ = 95% for numerical validation.

Figure 12. Feedrate, acceleration, contour error, tolerance violation, and prediction error using conservative (Cons.) physics-based (PB), data-driven (DD) and proposed (PIDD) methods with $\eta$ = 95% for numerical validation.

Figure 13. Feedrate, acceleration, contour error, tolerance violation, and prediction error using the proposed (PIDD) methods with $\eta$ = 95% (from Figure 12 and $\eta$ = 98% for numerical validation.

Figure 13. Feedrate, acceleration, contour error, tolerance violation, and prediction error using the proposed (PIDD) methods with $\eta$ = 95% (from Figure 12 and $\eta$ = 98% for numerical validation.

Table 1. Cycle times and RMS of tolerance violation $\gamma$ for conservative (Cons.), physics-based (PB), data-driven (DD) and proposed (PIDD) methods in Figure 12 and Figure 13.

6. Experimental validation

6.2. CNC machine tool

6.2.1. Experimental setup

For experimental validation, the same experimental setup with the machine tool in Figure 6 in used. The optimization algorithm is implemented on dSPACE 1007 real-time control board running at 500Hz sampling rate, connected to DRV8825 stepper motor drivers for x, y, and z- axes stepper motors. Renishaw RKLC20-S optical linear encoders are attached to the x and y- axes gantries to measure x and y- axes positions that are fed back to dSPACE 1007.

6.2.2. Experimental results

This section validates proposed feedrate optimization using the same set of methods for command generation, which are conservative, PB, DD and the proposed PIDD methods. The same desired butterfly trajectory in Figure 11 is used for air cutting and machining an aluminum workpiece with a 3.175 mm diameter flat-end mill and spindle speed of 7000 rpm. Kinematic limits are set as $V_{max}$ = 20 mm/s and $A_{max}$ = 0.5 m/$s^2$, and contour error bound is chosen as $E_{max}$ = 0.1 mm in the feedrate optimization; $n_w$ = 30, $n_2$ = 2, $n_3$ = 30 and ${\sigma }_{ϵ}$ = 0.01 are used in the DD and PIDD methods. The desired stringency is chosen as $\eta$ = 95% in the DD and PIDD method.

Figure 19 shows the profiles of optimized feedrate, acceleration, contour error, tolerance violation, and prediction error of x- and y-axis in air-cutting using the conservative, PB, DD and PIDD approaches. The PB method frequently violates the tolerance due to unmodeled dynamics, which is caused by the significant prediction error. The DD method slightly improves prediction accuracy, which is further improved in the PIDD method where the contour error is able to be constrained close to the limit using the desired stringency, similar to conservative approach. The proposed PIDD approach completes the motion 38.06% faster than the conservative method, while maintaining a similar level of tolerance adherence. The RMS prediction errors in x- and y-axis, cycle times and RMS tolerance violations of each method are summarized in Table 3.

Figure 20 shows the profiles of optimized feedrate, acceleration, contour error, tolerance violation, and prediction error of x- and y-axis in actual cutting using the conservative, PB, DD and PIDD approaches. Similar to air-cutting, the PB method is worst in constraining the contour error due to unmodeled dynamics and/or cutting force. The DD and PIDD methods reduce the prediction error compared to the PB method and are able to constrain the contour error close to the limit using the desired stringency. However, occasionally, PIDD shows worse performance than DD, which may be due to the difference between measured FRFs when the machine tool is not cutting (shown in Figure 7 and Figure 8) and the FRFs while it is cutting. Research has shown that there can be significant differences between FRFs measured without cutting and those measured while cutting [44]. A possible solution solve this issue is to measure the FRFs and compute known uncertainties during cutting using operational modal analysis, which may be complicated because of FRF’s variability on the operating condition of the cutting tool. Another solution is to discard the known uncertainties when they are inaccurate, i.e., use DD only when PIDD may have errors. Overall, the proposed PIDD approach completes the motion 29.02% faster than the conservative method while maintaining a similar level of tolerance adherence. The RMS prediction errors in x- and y-axis, cycle times and RMS tolerance violations of each method are summarized in Table 3.

Figure 19. Commanded feedrate, acceleration, contour error, tolerance violation, and prediction error (with its zoomed-in images) of conservative, physics-based (PB), data-driven (DD) and proposed (PIDD) approaches in air-cutting.

Figure 19. Commanded feedrate, acceleration, contour error, tolerance violation, and prediction error (with its zoomed-in images) of conservative, physics-based (PB), data-driven (DD) and proposed (PIDD) approaches in air-cutting.

Figure 20. Commanded feedrate, acceleration, contour error, tolerance violation, and prediction error (with its zoomed-in images) of conservative, physics-based (PB), data-driven (DD) and proposed (PIDD) approaches in actual cutting.

Figure 20. Commanded feedrate, acceleration, contour error, tolerance violation, and prediction error (with its zoomed-in images) of conservative, physics-based (PB), data-driven (DD) and proposed (PIDD) approaches in actual cutting.

Table 3. Comparison of RMS prediction errors, cycle times and RMS of tolerance violation $\gamma$ using conservative (Cons.), physics-based (PB), data-driven (DD) and proposed (PIDD) methods.

Table 3. Comparison of RMS prediction errors, cycle times and RMS of tolerance violation $\gamma$ using conservative (Cons.), physics-based (PB), data-driven (DD) and proposed (PIDD) methods.

		Cons.	PB	DD	PIDD
Air-cutting	x Pred. Error [$\mu$m]	N/A	17.8	17.5	10.2
	y Pred. Error [$\mu$m]	N/A	25.4	15.7	7.2
	Cycle time [s]	40.2	20.0	21.1	24.9
	RMS of $\gamma$ [$\mu$m]	5.0	16.7	16.5	2.3
Cutting	x Pred. Error [$\mu$m]	N/A	47.2	17.3	16.7
	y Pred. Error [$\mu$m]	N/A	34.2	17.7	18.1
	Cycle time [s]	39.4	20.0	23.0	27.2
	RMS of $\gamma$ [$\mu$m]	6.8	22.2	3.7	7.4

Figure 21 shows the machined surfaces and their zoomed-in images of upper left wing using the trajectories from Figure 20. The surface machined using PB method shows vibration marks, while the DD and PIDD methods mitigate the vibration and achieve similar quality to that of conservative approach.

7. Conclusion and future work

• A novel uncertainty-aware digital twin that predicts contour error is proposed. The digital twin is able to incorporate known uncertainty from physics-based models and learn unknown uncertainty using an online data-driven model to predict contour error’s distribution.
• For the first time, a feedrate optimization with constraints on kinematics and contour error using quantified uncertainty is introduced. The contour error’s uncertainty using digital twin enables the manufacturer to impose stringency constraints, which can replace trial-and-error approach of tuning the tolerance used in practice.
• We have demonstrated the effectiveness of the intelligent feedrate optimization using uncertainty-aware digital twin, to show up to 38% and 17% cycle time reduction using a desktop CNC machine tool and a desktop 3D printer, respectively, while achieving similar stringency of tolerance to that of the a conservative trial-and-error approach similar to those used in practice.
• The proposed intelligent feedrate optimization is expected to bring impact in achieving desired quality with higher productivity, using less trial-and-error. It is applicable to any machines that use feed drives, such as coordinate measurement machines (CMMs), and precision machine tools.

References

1. HC, M.; P, W.; K, E.; Y, K. Self-optimizing machining systems. CIRP Annals 2020, 69, 740–763. [CrossRef]
2. W, F.; XS, G.; CH, L.; K, Z.; Q, Z. Time-Optimal Interpolation for Five-Axis CNC Machining along Parametric Tool Path Based on Linear Programming. Int J Adv Manuf Technol 2013, 69, 1373–1388. [CrossRef]
3. K, E.; QGC, C.; Zhao MY, B.X.; XS, G. Linear Programming and Windowing Based Feedrate Optimization for Spline Toolpaths. CIRP Annals 2017, 66, 393–396. [CrossRef]
4. A, B.; J, D. Feedrate Optimization for Smooth Minimum-Time Trajectory Generation with Higher Order Constraints. Int J Adv Manuf Technol 2016, 82, 1029–1040. [CrossRef]
5. Q, Z.; SR, L. Efficient Computation of Smooth Minimum Time Trajectory for CNC Machining. Int J Adv Manuf Technol 2013, 68, 683–692. [CrossRef]
6. S, T.; B, S. Real-time trajectory generation for 5-axis machine tools with singularity avoidance. CIRP Annals 2020, 69, 349–352. [CrossRef]
7. S, T.; B, S. Online interpolation of 5-axis machining toolpaths with global blending. International Journal of Machine Tools and Manufacture 2022, 175, 103862. [CrossRef]
8. M, C.; J, X.; Y, S. Adaptive feedrate planning for continuous parametric tool path with confined contour error and axis jerks. The International Journal of Advanced Manufacturing Technology 2017, 89, 1113–1125. [CrossRef]
9. J, D.; JA, S. Optimal feed-rate scheduling for high-speed contouring. J. Manuf. Sci. Eng. 2007, 129, 63–76. [CrossRef]
10. W, F.; XS, G.; CH, L.; K, Z.; Q, Z. Time-optimal interpolation for five-axis CNC machining along parametric tool path based on linear programming. The International Journal of Advanced Manufacturing Technology 2013, 69, 1373–1388. [CrossRef]
11. H, E.; I, L.; B, O. Feedrate scheduling strategies for free-form surfaces. International Journal of Machine Tools and Manufacture 2006, 46, 747–757. [CrossRef]
12. L, Z.; J, F.; Y, W.; M, C. Feedrate scheduling strategy for free-form surface machining through an integrated geometric and mechanistic model. International Journal of Advanced Manufacturing Technology 2009, 40, 740–763. [CrossRef]
13. SD, M.; Y, A. Virtual cutting and optimization of three-axis milling processes. International Journal of Machine Tools and Manufacture 2008, 48, 1063–1071. [CrossRef]
14. ND, R.; BK, F.; RB, J. Efficient NC machining using off-line optimized feedrates and on-line adaptive control. 2002, Vol. 3641, pp. 181–191. [CrossRef]
15. H, E.; I, L.; M, K. Free-form surface machining and comparing feedrate scheduling strategies. Machining Science and Technology 2007, 11, 117–133. [CrossRef]
16. HS, P.; B, Q.; DV, D.; DY, P. Development of smart machining system for optimizing feedrates to minimize machining time. Journal of Computational Design and Engineering 2018, 5, 299–304. [CrossRef]
17. Y, A.; Kersting P, B.D.; E, B.; B, D.; I, L. Virtual process systems for part machining operations. CIRP Annals 2014, 63, 585–605. [CrossRef]
18. HU, L.; DW, C. An intelligent feedrate scheduling based on virtual machining. The International Journal of Advanced Manufacturing Technology 2003, 22, 873–882.
19. BK, F.; RB, J.; JG, H. Robust feedrate selection for 3-axis NC machining using discrete models. J. Manuf. Sci. Eng. 2001, 123, 214–224. [CrossRef]
20. Y, A.; SD, M. Virtual high performance milling. CIRP annals 2007, 56, 81–84. [CrossRef]
21. JH, K.; WS, Y.; DW, C. Off-line feed rate scheduling using virtual CNC based on an evaluation of cutting performance. Computer-Aided Design 2003, 35, 383–393. [CrossRef]
22. SD, M.; Y, A. Virtual simulation and optimization of milling applications—part II: optimization and feedrate scheduling. [CrossRef]
23. HY, F.; N, S. Integrated tool path and feed rate optimization for the finishing machining of 3D plane surfaces. International Journal of Machine Tools and Manufacture 2000, 40, 1557–1572. [CrossRef]
24. CM, B.; M, L.; P, B. A model-based adaptive controller for chatter mitigation and productivity enhancement in CNC milling machines. Robotics and Computer-Integrated Manufacturing 2016, 40, 34–43. [CrossRef]
25. P, Q.; M, W.; L, S. Feed Rate Variation Strategy for Semi-Conical Shell Workpiece in Ball Head End Milling Process. Applied Sciences 2020, 10, 9135. [CrossRef]
26. J, Y.; D, A.; Y, A. A feedrate scheduling algorithm to constrain tool tip position and tool orientation errors of five-axis CNC machining under cutting load disturbances. CIRP journal of Manufacturing Science and Technology 2018, 23, 78–90. [CrossRef]
27. Y, A.; J, Y.; ZM, K. Virtual prediction and constraint of contour errors induced by cutting force disturbances on multi-axis CNC machine tools. CIRP Annals 2019, 68, 377–380. [CrossRef]
28. H, K.; CE, O. Simultaneous servo error pre-compensation and feedrate optimization with tolerance constraints using linear programming. Int J Adv Manuf Technol 2020, 109, 809–821. [CrossRef]
29. H, K.; CE, O. Accurate and computationally efficient approach for simultaneous feedrate optimization and servo error pre-compensation of long toolpaths—with application to a 3D printer. The International Journal of Advanced Manufacturing Technology 2021, 115, 2069–2082. [CrossRef]
30. JB, R.; DQ, M.; M, D. Model predictive control: theory, computation, and design; Vol. 2, Madison, WI: Nob Hill Publishing.
31. H, K.; C, O. Intelligent feedrate optimization using a physics-based and data-driven digital twin. CIRP Annals. 2023, pp. 308–322. [CrossRef]
32. YC, C.; CW, C.; TC, T. Near time-optimal real-time path following under error tolerance and system constraints. CIRP Annals 2018, 140, 071004. [CrossRef]
33. S, B.; D, L.M.; A, R.; J, L. Data-driven Reference Trajectory Optimization for Precision Motion Systems. arXiv preprint 2205.15694 2022.
34. L, R.; I, L.; ED, K.; HC, M. Safe optimization for feedrate scheduling of power-constrained milling processes by using Gaussian processes. Procedia CIRP 2021, 99, 127–132. [CrossRef]
35. R, K.; G, R.; Z, S. Minimizing negative transfer of knowledge in multivariate gaussian processes: A scalable and regularized approach. IEEE Transactions on Pattern Analysis and Machine Intelligence 2020, 43, 3508–3522. [CrossRef]
36. CH, C.; M, D.; CE, O. A linear hybrid model for enhanced servo error pre-compensation of feed drives with unmodeled nonlinear dynamics. CIRP Annals 2021, 70, 301–304. [CrossRef]
37. D, L.; C, M.; MC, G. Model Predictive Contouring Control for Biaxial Systems. IEEE Transactions on Control Systems Technology 2012, 2, 552–559. [CrossRef]
38. DJ, M. Information theory, inference and learning algorithms; Cambridge university press, 2003.
39. JB, R. Model predictive control technology. 1999, Vol. 1, pp. 662–676.
40. DJ, G.; K, E. Accurate control of ball screw drives using pole-placement vibration damping and a novel trajectory prefilter. Precision Engineering 2013, 37, 308–322. [CrossRef]
41. HW, L. On the Kolmogorov-Smirnov test for normality with mean and variance unknown. Journal of the American Statistical Association 1967, 62, 339–402. [CrossRef]
42. S, Y.; AH, G.; X, L.; CE, O. Pre-Compensation of Servo Contour Errors Using a Model Predictive Control Framework. Int J Mach Tool Manu 2015, 98, 50–60. [CrossRef]
43. DG, L. An Introduction to Observers. IEEE Trans. On Automatic Control 1971, 16, 596–602. [CrossRef]
44. I, Z.; V, S. Estimation of machine-tool dynamic parameters during machining operation through operational modal analysis. International Journal of Machine Tools and Manufacture 2009, 49, 945–957. [CrossRef]