Intelligent Feedrate Optimization using an Uncertainty-aware Digital Twin within a Model Predictive Control Framework

1. Introduction

Quality and productivity are two important and frequently competing factors in manufacturing. As a result, manufacturers strive to maximize productivity while adhering to quality constraints. In practice, achieving this goal has involved a trial-and-error approach. However, there is a growing demand for self-optimizing intelligent manufacturing machines that have the capability to autonomously optimize their speed while ensuring desired quality levels, eliminating the need for extensive trial-and-error [14].

Motion-induced servo error, which is one of the major sources of quality degradation in manufacturing machines, can result from the limited bandwidth of feedback controllers, flexible structures, nonlinear friction, and backlash. Another source of servo error is process force, such as cutting force. Given that motion- and process-force-induced servo errors tend to increase with higher speeds, there is a keen interest in maximizing the speed of motion while respecting the tolerances on motion- and/or process-induced servo errors.

Extensive research has been conducted in the field of feedrate optimization with the objective of maximizing the feedrate while respecting servo error constraints. The majority of feedrate optimization or feedrate scheduling techniques primarily focus on maximizing the feedrate while considering kinematic limits such as speed, acceleration, and jerk [1,25,31,34,35,39]. However, the existing studies in [1,25,31,34,35,39] do not incorporate dynamic constraints such as servo error and cutting force, resulting in the need for a cautious selection of kinematic limits to indirectly meet the requirements on dynamic constraints. This indirect approach is necessary due to the complex relationship between kinematics and dynamic constraints, which often leads to sub-optimal feedrates.

In order to directly enforce dynamic constraints, certain feedrate optimization methods incorporate limits on the servo error, in addition to kinematics, by using steady-state [28] or static [20,40] approximations of servo models associated with motion velocity and acceleration. However, their limited ability to directly incorporate dynamic aspects of servo error, such as dynamic servo error pre-compensation, hinders their accuracy and effectiveness in optimizing feedrate.

To directly incorporate dynamic components via physics-based models, numerous feedrate scheduling methods for CNC machines maximize feedrate in each NC block while keeping cutting force under desired levels via mechanistic force models [9,10,15,16,27,29,38,42]. Some feedrate scheduling techniques maximize feedrate while regulating machining error due to tool deflection [2,4,18,24,30,37,43] or force-induced servo error [21,41] under desired tolerance in CNC machine tools. A few works in feedrate optimization [12,13] constrain motion-induced error via linear physics-based models of servo dynamics. However, the works in [2,4,9,10,12,13,15,16,18,21,24,27,29,30,37,38,41,42,43] are unable to effectively constrain actual cutting force or servo error in situations where uncertainties arise from nonlinear dynamics or disturbances that are not incorporated in the physics-based models. As a result, their capability to maximize feedrate while adhering to dynamic constraints is severely restricted.

There is increasing interest in the utilization of digital twins in manufacturing. A digital twin is a virtual representation, parallel to a physical system, built on a bi-directional link between simulation and actual data collection [14]. To effectively optimize feedrate with existence of uncertainties, digital twins can be used to provide more-accurate predictions of process or servo dynamics for feedrate optimization using data-driven models updated via sensor measurements. Model predictive control (MPC) is a framework often used for optimizing feedrate using digital twins. In MPC [23], predictions made using physics-based and data acquired from sensors are used to optimize an objective in a receding horizon fashion. In the context of feedrate optimization, a linear hybrid model was augmented with a periodic internal model in a MPC framework [11] to effectively predict and constrain servo errors due to motion and cutting forces. Luenberger state observers [44] were used in feedrate optimization to correct the initial system states of servo dynamic models in real-time for accurate contour error constraint enforcement, where the objective function was based on a tunable index of how far away an unattainable target position is from the current position. Similarly, a two-stage MPC approach was proposed in [33], where the first stage performed feedrate optimization with contour error constraints using a linear data-driven model, and the second stage further pre-compensated the contour error using artificial neural networks (ANNs). However, the data-driven models in [33] were trained offline through numerous experiments, which is time-consuming and may not be effective in predicting contour error in the presence of in-situ uncertainties that are not included in the training data. Moreover, the two-stage optimizations in [33] and the mixed objective functions used to balance quality and productivity in [33,44] require trial-and-error to tune the objective function weights to determine an acceptable quality level. Therefore, it is better to enforce quality requirements as contour error constraints that must be met, as is often the case in practice. Furthermore, the existing MPC methods in [11,33,44] do not quantify or exploit the uncertainty of the prediction in their feedrate optimization. Hence, they may not effectively adhere to constraints in the presence of high uncertainty due to a lack of training data, system variability and sudden changes in operating conditions.

To quantify uncertainties and impose robustness, studies exist on maximizing feedrate while regulating spindle power, where the spindle power is modeled using Gaussian process regression (GPR) [26]. The spindle power constraint is derived from a stochastic constraint with a fixed confidence level to safely optimize feedrate in uncertain environments. However, GPR in [26] is updated cycle-by-cycle, where the first cycle is initialized with a conservative feedrate profile followed by numerous sequential cycles for convergence. This requires optimizing a highly non-linear GPR objective [32] to estimate hyperparameters, which renders the process less adaptable to real-time control. Moreover, it does not constrain the servo error nor account for the propagation of model uncertainties to the servo error. This oversight is critical in achieving desired part accuracy in feedrate optimization, especially for toolpaths with high curvatures that can create significant structural vibration.

To address the shortcomings of the existing works, this paper proposes an intelligent feedrate optimization with contour error constraints using an uncertainty-aware digital twin, by making the following original contributions:

• It uses a novel physics-informed data-driven digital twin to predict the contour error and its uncertainty, where the digital twin incorporates the known uncertainty in the physics-based models and learns the unknown uncertainty by correcting the data-driven model on-the-fly using sensor feedback.
• It formulates an intelligent feedrate optimization framework capable of maximizing feedrate while accurately constraining contour error under desired tolerance and stringency, based on the uncertainty estimation from the digital twin using a model predictive control framework.
• It demonstrates the effectiveness of the proposed method by validating its performance in simulations and experiments using a desktop CNC machine tool and 3D printer.

The outline of the paper is as follows: Section 2 introduces our framework for intelligent feedrate optimization using an uncertainty-aware digital twin. Section 3 describes the methodology of the proposed digital twin for predicting the contour error distribution. Section 4 provides a formulation of the feedrate optimization with contour error and stringency constraints (i.e., the tolerance for constraint violation under uncertainty). Section 5 numerically validates the proposed method via a desktop CNC machine tol. Section 6 experimentally validates the proposed method via a desktop CNC machine tool and 3D printer. Finally, Section 7 concludes the paper and discusses future work.

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

The framework for the proposed intelligent feedrate optimization with an uncertainty-aware digital twin is depicted in Figure 1. First, a manufacturer submits a part together with the desired dimensions and contour error tolerance to an intelligent manufacturing machine. Then, the goal of the machine is to autonomously produce the part as quickly as possible while respecting the given error tolerance. The machine is equipped with a digital twin that predicts the contour error, which the machine can exploit for feedrate optimization with contour error constraints. Hence, the proposed framework is based on model predictive control.

However, several uncertainties exist in the physical system. Some portions are known from available data or expert knowledge, while others, such as nonlinear dynamics, may be unknown. If not considered, the known and unknown uncertainties cause a violation in the contour error tolerance, hence the part quality, as illustrated in Figure 2(a).

However, given that uncertainty exists in enforcing tolerance constraints, manufacturers have different levels of stringency in enforcing constraints. For example, a manufacturer may want at least 99% of the produced parts to satisfy the constraints. Therefore, stringency reflects a manufacturer’s tolerance for quality constraint violations under uncertainty. In this paper, we propose that, instead of relying on trial-and-error, manufacturers can impose a desired stringency $\eta$% on the given tolerance by incorporating the uncertainty of the contour error prediction as shown in Figure 2(b). Imposing the stringency represents constraining the worst case out of $\eta$% of the entire variation of contour error, so that $\eta$% of the manufactured parts adhere to the imposed kinematic and tolerance constraints under the given uncertainty.

To do so, the digital twin uses the known uncertainty from the physics-based models and trains a data-driven model using the machine’s sensor measurements to learn the unknown uncertainty. The digital twin then predicts and quantifies the uncertainty of the contour error, which is used in feedrate optimization with desired tolerance and stringency on the contour error. Together with the uncertainty-aware digital twin, the feedrate optimization determines the fastest feedrate to run the machine while respecting the limits for the contour errors (and the kinematic limits of the machine) in a robust way. The measured sensor output is compared with the predicted output and used to adjust the digital twin and optimization algorithm in the next iteration of feedrate optimization.

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

This section first describes the methodology of the contour error prediction using a deterministic digital twin in Section 3.1, and extends it to prediction of contour error and its distribution using an uncertainty-aware digital twin in Section 3.2.

3.1. Overview of Contour Error Prediction Using a Deterministic Digital Twin

A flowchart of the intelligent feedrate optimization using a deterministic digital twin based on the previous work in [11] is depicted in Figure 3. Note that the internal model in [11] is removed. Also note that Figure 3, as well as the discussions that follow in this section, focus on only the x-axis for the sake of brevity. However, the exact same process can be applied to the y-axis and other independent axes of a Cartesian configured manufacturing machine. Small batches (i.e., look-ahead windows) ${\mathit{x}}_d^j$ of a desired position trajectory ${\mathit{x}}_d$ are fed into an intelligent feedrate optimizer to produce the optimized motion command, ${\mathit{x}}_c^j$, where the superscript $j\in \{0,1,2,\cdots \}$ represents the batch index. Here, each batch has a window length of $n_w$ and is defined on the time domain $t\in \{0,T_s,2T_s,3T_s,\cdots \}$ which represents discrete time at sampling interval $T_s$, as illustrated in Figure 4. The optimized motion commands are sent to the servo system ${\mathit{H}}_x$ to produce the actual position ${\mathit{x}}^j$. The servo system is composed of a servo error pre-compensation ${\mathit{C}}_x$ followed by machine dynamics ${\mathit{G}}_x$, i.e., ${\mathit{H}}_x={\mathit{G}}_x{\mathit{C}}_x$.

A key requirement for the intelligent feedrate optimization is accurate prediction of the servo error, which is achieved using linear regression built upon the hybrid model presented in [3]. The hybrid model takes input ${\mathit{x}}_c^j$ and predicts the actual position ${\widehat{\mathit{x}}}^j$ using a stable, nominal (or representative) physics-based model ${\widehat{\mathit{H}}}_x^{\ast }$. The predictions ${\widehat{\mathit{x}}}^j$ do not capture the effects of unmodeled dynamics and external disturbances. Therefore, the prediction error of ${\widehat{\mathit{H}}}_x^{\ast }$ (delayed by one batch) is computed as ${\mathit{e}}_x^{j-1}={\mathit{x}}^{j-1}-{\widehat{\mathit{x}}}^{j-1}$ and combined with ${\widehat{\mathit{x}}}^j$ and ${\widehat{\mathit{x}}}^{j-1}$ to be fed into a data driven model to generate an improved prediction ${\tilde{\mathit{x}}}^j$ which is used for constraining contour errors in the feedrate optimization.

Each element $\tilde{x}\left(t\right)$ of $\tilde{\mathit{x}}$ is modeled as

$$\begin{array}{c}\hfill \tilde{x}\left(t\right)={\mathit{\beta }}^{\mathrm{T}}{\mathit{\psi }}_t\end{array}$$

(1)

where ${\mathit{\psi }}_t$ is the deterministic feature vector and $\mathit{\beta }$ is the weight vector that is learned using linear regression. Note that the superscript j, i.e., batch index, has been removed from $\tilde{x}\left(t\right)$ in Eq. (1) to directly define $\tilde{\mathit{x}}$ in time domain. The sub elements of ${\mathit{\psi }}_t$ are given by

$$\begin{array}{cc}\hfill {\mathit{\psi }}_t=[\underset{:={\mathit{\psi }}_{t1}}{\underbrace{1}},& \underset{:={\mathit{\psi }}_{t2}}{\underbrace{\widehat{x}(t-n_2{T}_s),\cdots ,\widehat{x}\left(t\right),}}\hfill \\ \hfill & \underset{:={\mathit{\psi }}_{t3}}{\underbrace{e_x(t-n_3{T}_s),\cdots ,e_x(t-T_s)}}{]}^{\mathrm{T}}\hfill \end{array}$$

(2)

The sub-elements ${\mathit{\psi }}_{t1}$, ${\mathit{\psi }}_{t2}$ and ${\mathit{\psi }}_{t3}$ were contained in the hybrid model of [3]. They respectively represent a bias term, the past $n_2$ and current time steps of $\widehat{\mathit{x}}$, and the past $n_3$ time steps of ${\mathit{e}}_x$, where $n_2$ and $n_3$ are user defined.

The weight $\mathit{\beta }$ is updated recursively in each timestep within a given window, where the final weight is carried over to the next window for prediction, i.e., ${\tilde{\mathit{x}}}^j$ is predicted based on weight $\mathit{\beta }$ from the previous batch $j-1$. The algorithm for learning $\mathit{\beta }$ is described as follows. For $t=0$, $\mathit{\beta }$ and its covariance matrix $\mathit{P}$ are initialized using ridge regression with a regularization factor $\lambda$ as

$$\begin{array}{cc}\hfill & \mathit{\beta }={(\lambda \mathrm{I}+{\mathit{\psi }}_t{\mathit{\psi }}_t^T)}^{-1}{\mathit{\psi }}_{t}x\left(t\right)\hfill \\ \hfill & \mathit{P}={(\lambda \mathrm{I}+{\mathit{\psi }}_t{\mathit{\psi }}_t^T)}^{-1}\hfill \end{array}$$

(3)

For the rest of the timesteps $t\in \{T_s,2T_s,...\}$, a recursive least-squares is used to correct $\mathit{\beta }$ and $\mathit{P}$ using a forgetting factor $f_0$ as

$$\begin{array}{cc}\hfill & \mathit{\beta }\leftarrow \mathit{\beta }+\mathit{k}(x\left(t\right)-{\mathit{\beta }}^{\mathrm{T}}{\mathit{\psi }}_t)\hfill \\ \hfill & \mathit{P}\leftarrow \frac{1}{f_0}(\mathit{P}-\mathit{k}{\mathit{\psi }}_t^{\mathrm{T}}\mathit{P})\hfill \\ \hfill & \mathrm{where} \mathit{k}=\mathit{P}{\mathit{\psi }}_t{(f_0+{\mathit{\psi }}_t^{\mathrm{T}}\mathit{P}{\mathit{\psi }}_t)}^{-1}\hfill \end{array}$$

(4)

Using the final weight in batch $j-1$ to substitute for $\mathit{\beta }$, ${\tilde{\mathit{x}}}^j$ can be predicted using the feature vector ${\mathit{\psi }}_t$ formulated by Eq. (2). Since the past sensor data ${\mathit{x}}^{j-1}$ is provided up to $t=(j{n}_w-1)T_s$, for entries in batch j that have unavailable terms in ${\mathit{\psi }}_{t3}$, ${\mathit{e}}_x$ for all batches is approximated using the predicted values of $\tilde{\mathit{x}}$, i.e.,

$${\mathit{e}}_x=\mathit{x}-\widehat{\mathit{x}}\approx \tilde{\mathit{x}}-\widehat{\mathit{x}}$$

(5)

The same procedure is applied to the y-axis to predict $\tilde{\mathit{y}}$. Lastly, the contour error $\tilde{\epsilon }$ can be estimated from the predicted axis tracking errors ${\tilde{\epsilon }}_x={\mathit{x}}_d-\tilde{\mathit{x}}$ and ${\tilde{\epsilon }}_y={\mathit{y}}_d-\tilde{\mathit{y}}$, using a linear approximation [5] as

$$\tilde{\epsilon }=-sin\left(\mathit{\theta }\right){\tilde{\epsilon }}_x+cos\left(\mathit{\theta }\right){\tilde{\epsilon }}_y$$

(6)

where $\mathit{\theta }$ is inclination angle of the curve (${\mathit{x}}_d,{\mathit{y}}_d$).

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 ${\mathit{\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.

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)={\mathit{\beta }}^{\mathrm{T}}{\mathit{\psi }}_t+ϵ$$

(7)

where ${\mathit{\psi }}_t\sim N({\mathit{\mu }}_{{\mathit{\psi }}_t},{\Sigma }_{{\mathit{\psi }}_t})$ is the feature vector defined as an uncorrelated Gaussian random variable with mean ${\mathit{\mu }}_{{\mathit{\psi }}_t}$ and variance ${\Sigma }_{{\mathit{\psi }}_t}$ derived from the uncertainty distribution of ${\mathit{\psi }}_t$ previously defined in Eq. (2). $\mathit{\beta }\sim N({\mathit{\mu }}_{\mathit{\beta }},{\Sigma }_{\mathit{\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({\mathit{\mu }}_{{\tilde{x}}^j},{\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 ${\mathit{\psi }}_t$ in Eq. (7) to enable efficient training of $\mathit{\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\Delta \omega$, of which $\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({\mathit{\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\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 ${\mathit{\psi }}_t\sim N({\mathit{\mu }}_{{\psi }_t},{\Sigma }_{{\psi }_t})$, where

$$\begin{array}{cc}\hfill {\mathit{\mu }}_{{\psi }_t}=& [\underset{:={\mathit{\mu }}_{{\psi }_{t1}}}{\underbrace{1,}}\underset{:={\mathit{\mu }}_{{\psi }_{t2}}}{\underbrace{{\mu }_{\widehat{x}(t-n_2{T}_s),}\cdots ,{\mu }_{\widehat{x}\left(t\right)},}}\hfill \\ \hfill & \underset{:={\mathit{\mu }}_{{\psi }_{t3}}}{\underbrace{x(t-n_3{T}_s)-{\mu }_{\widehat{x}(t-n_3{T}_s)},\cdots ,x(t-T_s)-{\mu }_{\widehat{x}(t-T_s)}}}{]}^{\mathrm{T}}\hfill \\ \hfill {\Sigma }_{{\psi }_t}=& \mathrm{diag}(\underset{:={\Sigma }_{{\psi }_{t1}}}{\underbrace{0,}}\underset{:={\Sigma }_{{\psi }_{t2}}}{\underbrace{{\sigma }_{\widehat{x}(t-n_2{T}_s),}^2\cdots ,{\sigma }_{\widehat{x}\left(t\right)}^2,}}\hfill \\ \hfill & \underset{:={\Sigma }_{{\psi }_{t3}}}{\underbrace{{\sigma }_{\widehat{x}(t-n_3{T}_s)}^2,\cdots ,{\sigma }_{\widehat{x}(t-T_s)}^2}}{)}^{\mathrm{T}}\hfill \end{array}$$

(13)

The prior over $\mathit{\beta }\sim N({\mathit{\mu }}_{\beta }$, ${\Sigma }_{\beta })$ is initialized using the nominal physics-based model ${\mathit{\mu }}_{{\widehat{h}}_x}$ from Eq. (11), as described in Appendix A. Now given $\mathit{\beta }\sim N({\mathit{\mu }}_{\beta },{\Sigma }_{\beta })$ and the feature uncertainty ${\mathit{\psi }}_t\sim N({\mathit{\mu }}_{{\psi }_t},{\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(\mathit{\beta }\right|x\left(t\right),{\mathit{\mu }}_{{\psi }_t},{\Sigma }_{{\psi }_t})& \propto p\left(x\left(t\right)\right|\mathit{\beta },{\mathit{\mu }}_{{\psi }_t},{\Sigma }_{{\psi }_t})·p\left(\mathit{\beta }\right)\hfill \\ \hfill & \propto \int p(x\left(t\right),{\mathit{\psi }}_t|\mathit{\beta },{\mathit{\mu }}_{{\psi }_t},{\Sigma }_{{\psi }_t})d{\mathit{\psi }}_t·p\left(\mathit{\beta }\right)\hfill \\ \hfill & \propto \int p\left(x\left(t\right)\right|{\mathit{\psi }}_t,\mathit{\beta })p\left({\mathit{\psi }}_t\right|{\mathit{\mu }}_{{\psi }_t},{\Sigma }_{{\psi }_t})d{\mathit{\psi }}_t\hfill \\ \hfill & ·p\left(\mathit{\beta }\right)\hfill \\ \hfill & \propto p\left(x\right(t\left)\right|\mathit{\beta })·p(\mathit{\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|\mathit{\beta })$ and $p\left(\mathit{\beta }\right)$ are given as as

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

(15)

As shown above, a central feature of our approach is that it accounts for input uncertainty through ${\mathit{\psi }}_t\sim N({\mathit{\mu }}_{{\psi }_t},{\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 [8] where we write log of $p\left(x\right(t\left)\right|\mathit{\beta })$ as a quadratic function of $\mathit{\beta }$. This is shown below:

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

(16)

where the local point $\overline{\mathit{\beta }}$ is ${\mathit{\mu }}_{\mathit{\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({\mathit{\mu }}_{\beta },{\Sigma }_{\beta })$, now Eq. (14) can re-written in a closed-form solution as

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

(17)

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

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

(18)

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

Note that, if the feature ${\mathit{\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 ${\mathit{\psi }}_t$ in the following Section 5. In the case of deterministic ${\mathit{\psi }}_t$, the posterior $\mathit{\beta }\sim N({\mathit{\mu }}_{\beta },{\Sigma }_{\beta })$ can be estimated as a closed-form solution using Bayes rule as

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

(20)

The same procedures are applied to learn ${\mathit{\beta }}_y\sim N({\mathit{\mu }}_{{\beta }_y},{\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 ${\mathit{\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(\theta \left(t\right)\right)(x_d\left(t\right)-{\mathit{\mu }}_{\beta }^{\mathrm{T}}{\mathit{\psi }}_t)\hfill \\ \hfill & +cos\left(\theta \left(t\right)\right)(y_d\left(t\right)-{\mathit{\mu }}_{{\beta }_y}^{\mathrm{T}}{\mathit{\psi }}_{ty})\hfill \\ \hfill & {\sigma }_{\tilde{\epsilon }\left(t\right)}^2=sin{\left(\theta \left(t\right)\right)}^2({\mathit{\psi }}_t^{\mathrm{T}}{\Sigma }_{\beta }{\mathit{\psi }}_t+{\sigma }_{ϵ}^2)\hfill \\ \hfill & +cos{\left(\theta \left(t\right)\right)}^2({\mathit{\psi }}_{ty}^{\mathrm{T}}{\Sigma }_{{\beta }_y}{\mathit{\psi }}_{ty}+{\sigma }_{{ϵ}_y}^2)\hfill \end{array}$$

(21)

4. Methodology for Intelligent Feedrate Optimization with Contour Error Constraints

The feedrate optimization with contour error constraints using the quantified uncertainty from the digital twin is formulated in accordance with authors’ previous work [12] using a model predictive control framework. Taking the x-axis, for example, a desired trajectory ${\mathit{X}}_d=f\left(p\right)$ is parametrized with respect to a normalized, monotonically increasing path variable $\mathit{p}$, which is a vectorized form of p. Then, $X_d\left(t\right)$ is linearized as $x_d\left(t\right)$ with respect to $p\left(t\right)$ using an estimated linearization point $\overline{p}\left(t\right)$ as

$$x_d\left(t\right)=-{\left.\frac{\partial f}{\partial p}\right|}_{p=\overline{p}\left(t\right)}·(p\left(t\right)-\overline{p}\left(t\right))+f\left(\overline{p}\left(t\right)\right)$$

(22)

The procedure for computing the optimal ${\mathit{p}}^j$ (corresponding to the optimal feedrate) using the uncertainty-aware digital twin is as follows. The path variable ${\mathit{p}}^j$ is maximized under monotonicity, maximum feedrate, and axis-acceleration constraints as

$$\begin{array}{c}\hfill max{\mathit{1}}^T{\mathit{p}}^j\\ \hfill \mathrm{s}.\mathrm{t}.& p(t-1)\le p\left(t\right)\le 1\hfill \\ \hfill & D\left[{\mathit{p}}^j\right]\le {\mathit{V}}_{max}{T}_s\hfill \\ \hfill & \left|D^2\left[{\mathit{x}}_d^j\right]\right|\le {\mathit{A}}_{max}{T}_s^2\hfill \end{array}$$

(23)

where $\mathit{1}$ is a ones-vector, D is a difference operator, and ${\mathit{V}}_{max}$ and ${\mathit{A}}_{max}$ are the vectorized representations of feedrate and acceleration limits, respectively. In addition, kinematic and dynamic continuity between adjacent windows is enforced. The process described above for the x-axis can be applied to the y-axis.

The feedrate optimization constrains the contour error under a given tolerance and stringency, using the posterior predictive distribution from Section 3.2. To do so, we show that ${\mathit{\mu }}_{\widetilde{ϵ}\left(t\right)}$ and ${\sigma }_{\widetilde{ϵ}\left(t\right)}$ are linear in terms of ${\mathit{x}}_d^j$, by showing that the only alterable feature in ${\mathit{\psi }}_t$, which is the last term in ${\mathit{\psi }}_{t2}$ (i.e., $\widehat{\mathit{x}}\left(t\right)$), is linear in ${\mathit{x}}_d^j$.

Let ${\mathbf{\Phi }}_x\in {\mathbb{R}}^{n_h\times n_h}$ be the matrix (lifted domain) representation of ${\mathit{\mu }}_{{\widehat{h}}_x}$ truncated by length $n_h$. The last $n_w$ rows in ${\mathbf{\Phi }}_x$ can further be decomposed into two parts: its first $n_h-n_w$ columns ${\mathbf{\Phi }}_{x,p}$ and its last $n_w$ columns ${\mathbf{\Phi }}_{x,c}$ as

$${\mathbf{\Phi }}_x=\left[\begin{array}{cc}⋮& ⋮\\ {\mathbf{\Phi }}_{x,p}& {\mathbf{\Phi }}_{x,c}\end{array}\right]$$

(24)

If ${\mathit{x}}_{c,p}$ represents the last $n_h-n_w$ elements of the ${\mathit{x}}_c$ at past timesteps, $\widehat{x}\left(t\right)$ can be re-written as

$$\begin{array}{cc}\hfill & {\widehat{\mathit{x}}}^j={\mathbf{\Phi }}_{x,c}{\mathit{x}}_d^j+{\mathbf{\Phi }}_{x,p}{\mathit{x}}_{c,p}\hfill \\ \hfill & \therefore \widehat{x}\left(t\right)=\underset{:={\mathit{T}}_x}{\underbrace{{\mathit{M}}_t{\mathbf{\Phi }}_{x,c}}}{\mathit{x}}_d^j+\underset{:={\mathit{T}}_{0x}}{\underbrace{{\mathit{M}}_t{\mathbf{\Phi }}_{x,p}{\mathit{x}}_{cp}}}\hfill \end{array}$$

(25)

where ${\mathit{M}}_t$ is a selection matrix that picks the entry at timestep t. Similarly, for y-axis, the alterable term in ${\mathit{\psi }}_{ty}$ can be derived to be linear in terms of ${\mathit{y}}_d^j$, by using a similar notation as Eq. (25), i.e., $\widehat{y}\left(t\right)={\mathit{T}}_y{\mathit{y}}_d^j+{\mathit{T}}_{0y}$.

Then, the worst-case out of the $\eta$ [%] variations of distribution of the contour error $\tilde{\epsilon }\left(t\right)$, where $\eta$ is a user-defined stringency, is bounded by the tolerance $E_{max}$ as a stochastic constraint by

$$p(\tilde{\epsilon }\left(t\right)\le E_{max})\ge \eta$$

(26)

For the sake of brevity, the negative stochastic contour error constraint $p(\tilde{\epsilon }\left(t\right)\ge -E_{max})\ge \eta$ is omitted. Then, inversion the both sides of Eq. (26) becomes

$${\mu }_{\tilde{\epsilon }\left(t\right)}\le E_{max}-{\mathrm{\Phi }}^{-1}\left(\eta \right){\sigma }_{\tilde{\epsilon }\left(t\right)}$$

(27)

where $\mathrm{\Phi }$ is the cumulative density function of the distribution of $\tilde{\epsilon }\left(t\right)$, which is invertible because $\tilde{\epsilon }\left(t\right)$ follows a Gaussian distribution as was shown in Eq. (21).

Eq. (27) can be rearranged as a linear constraint in terms of ${\mathit{x}}_d^j$ and ${\mathit{y}}_d^j$ using Eq. (21), written as

$${\mathit{U}}_x{\mathit{x}}_d^j+{\mathit{U}}_y{\mathit{y}}_d^j\le {\mathit{U}}_0$$

(28)

where the derivation of coefficients ${\mathit{U}}_x$, ${\mathit{U}}_y$ and ${\mathit{U}}_0$ is described in Appendix C. Finally, the contour error constraint is also linear with respect to the decision variable ${\mathit{p}}^j$ using the relationship in Eq. (22).

The methodology of feedrate optimization described in this section can be broadly considered as model predictive control [22], because it: (1) optimizes manipulatable inputs, e.g., desired trajectory, over a finite, receding horizon using (2) predictions of the dynamical system’s behavior through a model that is updated via feedback.

5. Numerical Validation of the Intelligent Feedrate Optimization Using Uncertainty-Aware Digital Twin

This section validates the importance of the uncertainty quantification of the proposed physics-informed data-driven (PIDD) uncertainty-aware digital twin in feedrate optimization, by comparing the method with the following approaches:

• Conservative method, which is defined as the benchmark generated using a trapezoidal acceleration profile [7] with kinematic limits tuned by trial-and-error to achieve the contour error tolerance with $\eta$% stringency, by allowing up to $(100-\eta )$% RMS violation normalized by $E_{max}$ defined in Section 4
• Physics-based (PB) method, which predicts the output position and its uncertainty using only the known uncertainty obtained from the set of physics-based models ${\left\{{\widehat{\mathit{H}}}_x^i\right\}}_{i=1}^{N_H}$ and ${\left\{{\widehat{\mathit{H}}}_y^i\right\}}_{i=1}^{N_H}$
• Data-driven (DD) method, which predicts the output position and its uncertainty by learning the unknown uncertainty without incorporating any known uncertainties, i.e., the prior ${\mathit{\mu }}_{{\beta }_0}$, ${\Sigma }_{\beta 0}$, ${\mathit{\mu }}_{{\beta }_{y}0}$ and ${\Sigma }_{{\beta }_{y}0}$ are initialized as zero at the 0-th batch, and $\mathit{\beta }$ and ${\mathit{\beta }}_y$ are learned via Bayesian linear regression for deterministic features in Eq. (19). Note that both the PB and DD methods are subsets of the proposed uncertainty-aware digital twin. However, we have separated them out to highlight the effect of introducing uncertainty in both the PB and the DD models through the PIDD method used in the uncertainty aware digital twin

A Nomad 3 three-axis desktop CNC machine tool is chosen as the simulated system, where its setup is shown in Figure 6. To analyze its known uncertainties, the position commands are generated and commanded by dSPACE DS1007 real-time control board running at 500 Hz sampling rate, connected to DRV8825 stepper motor drivers for the x-, y- and z-axes stepper motors. ADXL335 accelerometers are attached on the x- and y-axis gantries to measure the x and y-axes acceleration. The known uncertainties are identified by measuring FRFs, of which the input is a swept sine acceleration command to the stepper motors, and the output is the relative acceleration between the x- and y-axis using the accelerometers. The operating condition under which the FRFs are measured is varied by modifying the input acceleration amplitude at discrete values: 2 m/${\mathrm{s}}^2$, 3 m/${\mathrm{s}}^2$ and 4 m/${\mathrm{s}}^2$, and 3 FRFs are measured per each acceleration amplitude to collect a total of $N_H=9$ FRFs per axes.

The set of FRFs ${\left\{{\widehat{\mathit{H}}}_x^i\right\}}_{i=1}^9$ of the x- and y-axis of the printer are shown in Figure 7 and Figure 8, respectively. The uncertainties in ${\widehat{\mathit{H}}}_x$ are then propagated to ${\widehat{\mathit{h}}}_x\sim N({\mathit{\mu }}_{{\widehat{h}}_x},{\Sigma }_{{\widehat{h}}_x})$ to initialize ${\mathit{\mu }}_{\beta }$ and ${\Sigma }_{\beta }$ and construct ${\mathit{\mu }}_{{\psi }_t}$ and ${\Sigma }_{{\psi }_t}$ in the physics-informed data-driven digital twin. To validate the hypothesis that the FRF coefficients of ${\left\{{\widehat{\mathit{H}}}_x\right\}}_{i=1}^9$ and ${\left\{{\widehat{\mathit{H}}}_y\right\}}_{i=1}^9$ belong to Gaussian distributions, Lilliefors test [17] is performed on every frequency ${\omega }_k$ for $k\in \{1,2,...,N_H\}$ of $a\left({\omega }_k\right)$ and $b\left({\omega }_k\right)$ in x- and y-axis to compute the test decisions at the 5% significance level and p-values. Figure 9 and Figure 10 show the FRF coefficients $a\left({\omega }_k\right)$ and $b\left({\omega }_k\right)$ for x- and y-axis, respectively, in the upper plots. The lower plots of Figure 9 and Figure 10 show p-values and Lilliefors test results, where 0 and 1 represent acceptance and rejection of the hypothesis, respectively. The figures imply that Gaussian hypothesis for both $a\left({\omega }_k\right)$ and $b\left({\omega }_k\right)$ is accepted at 90% and 88% of the frequencies for the x- and y-axis, respectively. One way to improve the reliability of Lillifors test result is to gather more data of FRFs. However, for the sake of simplicity and computational efficiency, the Gaussian hypothesis will be assumed valid for all frequencies in the following sections.

The output position $\mathit{x}$ is simulated as the sum of motion-induced position ${\mathit{x}}_m$ and force-induced position ${\mathit{x}}_f$, as

$$\begin{array}{cc}\hfill & \mathit{x}={\mathit{x}}_m+{\mathit{x}}_f={\widehat{\mathit{h}}}_x\ast {\mathit{x}}_c+{\mathit{x}}_f\hfill \\ \hfill & \mathrm{where}{\widehat{\mathit{h}}}_x\sim N({\mathit{\mu }}_{{\widehat{h}}_x},{\Sigma }_{{\widehat{h}}_x})\mathrm{is}\mathrm{sampled}\mathrm{at}\mathrm{every}t\hfill \\ \hfill & \mathrm{and}{x}_f\left(t\right)=A_{f}sin{\omega }_{f}t\hfill \end{array}$$

(29)

where $A_f$ = 0.2 and ${\omega }_f$ = 733 rad/s (7000 rpm) are chosen.

The butterfly trajectory [36] with its contour of the toolpath on the x- and y-axis shown in Figure 11 is selected. For the DD and PIDD methods, $n_w$ = 10, $n_2$ = 3, $n_3$ = 10 and ${\sigma }_{ϵ}$ = 0.01 are selected. For the DD and PIDD methods, stringency $\eta$ = 95% is selected. $V_{max}$ = 30 mm/s, $A_{max}$ = 5 m/${\mathrm{s}}^2$, and contour error limit of $E_{max}$ = 0.4 mm are selected for the feedrate optimization. The tolerance violation $\gamma$, which will be analyzed for each method, is defined as

$$\gamma \left(t\right)=\{\begin{array}{cc}\left|\epsilon \left(t\right)\right|-E_{max}\hfill & \mathrm{if} \left|\epsilon \left(t\right)\right|>E_{max}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}$$

(30)

Figure 12 shows the optimized feedrate, acceleration, contour error, tolerance violation and prediction error of all methods. The cycle times and RMS of tolerance violation $\gamma$ are summarized in Table 1. The PB method is the worst in prediction performance because it is not aware of the unknown uncertainties caused by the force-induced servo error, and hence results in the highest RMS tolerance violation. The DD method improves adherence to the tolerance by learning the unknown uncertainties over time. However, DD method initially suffers from significant prediction error due to its unawareness of known uncertainties. The proposed PIDD method with $\eta$ = 95% enables restriction of the contour error under the desired stringency by incorporating known uncertainties and learning unknown uncertainties the quickest, which enables it to conservatively stay below the error limit most of the time. Overall, the PIDD method is able to reduce cycle time by 19.3% compared to the conservative approach while maintaining a similar tolerance violation level. To demonstrate the effect of the selection of stringency, Figure 13 compares the commanded feedrate, acceleration, contour error, tolerance violation and prediction error of the PIDD methods using $\eta$ = 95% and 98%. It is observed that tuning $\eta$ to a higher level has the effect of making the optimized feedrate more conservative and reducing the error violation.

Note that the proposed PIDD method is not perfect in satisfying the contour error constraints. One reason is that the prediction error is not perfectly zero, and the stringency constraints can only ensure that the worst-case out of $\eta$% of contour error distribution is within the tolerance. This issue can be mitigated by increasing $\eta$, which will entail more conservative feedrate. Another reason might be due to the nonlinear effects neglected by linearization of the contour error constraint in Eq. (33) and sub-optimal learning in $\mathit{\beta }$ introduced by Laplace’s approximation in Eq. (16). These problems can be addressed by applying nonlinear optimization and non-Gaussian Bayesian regression, at the price of higher computational cost.

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.

	Cons.	PB	DD	PIDD ($\mathit{\eta }$=95%)	PIDD ($\mathit{\eta }$=98%)
RMS of $\gamma$ [$\mathit{\mu }$m]	2.2	6.6	3.0	1.7	0.8
Cycle time [s]	8.89	4.49	5.17	6.47	7.17

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.

	Cons.	PB	DD	PIDD ($\mathit{\eta }$=95%)	PIDD ($\mathit{\eta }$=98%)
RMS of $\gamma$ [$\mathit{\mu }$m]	2.2	6.6	3.0	1.7	0.8
Cycle time [s]	8.89	4.49	5.17	6.47	7.17

6. Experimental Validation

For validation of the proposed approach, two experimental setups are used. The first set of experiments, described in Section 6.1, is carried out on an Ender 3 Pro desktop 3D printer, and the second set of experiments, described in Section 6.2, is carried out on a Nomad 3 desktop CNC machine tool used in Section 5. Demonstration of the proposed method on two experimental setups helps to show its versatility.

6.1. Desktop 3D Printer

6.1.1. Experimental Setup

The experimental setup using an Ender 3 Pro desktop 3D printer is shown in Figure 14. 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, z and e- axes stepper motors. The measured x and y- axes accelerations from ADXL335 accelerometers are fed back to dSPACE 1007. To deduce the x and y axes displacement from acceleration measurements, a Luenberger state observer [6] is used. The observer gains are chosen such that the dynamics of the observer error (i.e., discrepancy between estimated position using the nominal physics-based model ${\mathit{\mu }}_{{\widehat{h}}_x}$ and observed position) obtains global asymptotic convergence with an observer frequency $f_o$ = 15 Hz.

6.1.2. Experimental Results

This section validates the proposed approach experimentally using the desktop 3D printer, by comparing its performance with conservative, PB and DD methods. The butterfly toolpath in Figure 11 is used to optimize the feedrate for air-printing (i.e., no material extrusion) and actual printing of the 3D printer. The known uncertainties of x- and y-axis of the printer are incorporated from FRFs in Figure 15 and Figure 16. Similar to Section 5, Lilliefors test is performed on $a\left({\omega }_k\right)$ and $b\left({\omega }_k\right)$ of the FRFs in the x- and y-axis to validate the Gaussian assumption at the 5% significance level and p-values, where the hypothesis acceptance rates are computed as 92% and 91% out of all frequencies for the x- and y-axis, respectively. For the DD and PIDD methods, $n_w$ = 30, $n_2$ = 10, $n_3$ = 30 and ${\sigma }_{ϵ}$ = 0.01 are used. For the PIDD method, the desired stringency is selected as $\eta$ = 95%. For feedrate optimization, $V_{max}$ = 70 mm/s, $A_{max}$ = 3 m/${\mathrm{s}}^2$ and $E_{max}$ = 0.1 mm are chosen.

Figure 17 shows the profiles of the optimized feedrate, acceleration, contour error and prediction error of x- and $y-$ axis using the conservative, PB, DD and PIDD methods. The RMS prediction errors, cycle times and RMS tolerance violation of all methods are reported in Table 2. The PB approach cannot predict the unknown uncertainties, and hence shows the most significant violation in the contour error. The DD method mitigates the violation by learning the unknown uncertainties, and the PIDD method further improves the accuracy by staying the closest to the tolerance with the desired stringency. As a result, the PIDD method completes the motion 17.8% faster than, while yielding similar contour error tolerance satisfaction as the conservative one.

To further validate our findings, a 3D-augmentation of the trajectory in Figure 11 with z-height 8 mm is printed using the same printer. Conservative, PB, DD and PIDD methods are applied at each layer of the print. Figure 18 shows the top and side views of the printed butterflies using the four methods, as well as their print times. The physics-based and data-driven prints show vibration marks in the side view, while the proposed and conservative prints are able to achieve vibration-free surface quality. Overall, the proposed method is able to achieve 15.51% print time reduction compared to the conservative approach while achieving similar print quality.

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 [19]. 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 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

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

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 [μm]	N/A	17.8	17.5	10.2
	y Pred. Error [μm]	N/A	25.4	15.7	7.2
	Cycle time [s]	40.2	20.0	21.1	24.9
	RMS of $\gamma$ [μm]	5.0	16.7	16.5	2.3
Cutting	x Pred. Error [μm]	N/A	47.2	17.3	16.7
	y Pred. Error [μm]	N/A	34.2	17.7	18.1
	Cycle time [s]	39.4	20.0	23.0	27.2
	RMS of $\gamma$ [μ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

This paper presents an MPC framework and methodology for the intelligent feedrate optimization using an uncertainty-aware digital twin. Its key contributions are summarized as follows.

• 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.

As a limitation, the proposed method has made several assumptions in the methodology to enable efficient computation, such as Gaussian distribution of frequency response function for computing the known uncertainty and the linearization of contour error constraints for solving the feedrate optimization as a sequential linear programming problem. The future work will explore more sophisticated (non-Gaussian) uncertainty distributions and nonlinear contour error constraints to improve the accuracy of the digital twin, at the expense of higher computational cost and non-closed form solutions. Furthermore, additional forms of learning to the uncertainty-aware digital twin, such as part-to-part or machine-to-machine learning, will be investigated to improve prediction accuracy.