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 [
1].
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 [
2,
3,
4,
5,
6,
7]. However, the existing studies in [
2,
3,
4,
5,
6,
7] 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 [
8] or static [
9,
10] 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 [
11,
12,
13,
14,
15,
16,
17,
18]. Some feedrate scheduling techniques maximize feedrate while regulating machining error due to tool deflection [
19,
20,
21,
22,
23,
24,
25] or force-induced servo error [
26,
27] under desired tolerance in CNC machine tools. A few works in feedrate optimization [
28,
29] constrain motion-induced error via linear physics-based models of servo dynamics. However, the works in [
11,
12,
13,
14,
15,
16,
17,
18,
19,
20,
21,
22,
23,
24,
25,
26,
27,
28,
29] 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 [
1]. 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 [
30], 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 [
31] to effectively predict and constrain servo errors due to motion and cutting forces. Luenberger state observers [
32] 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 [
32,
33] 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 [
31,
32,
33] 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) [
34]. 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 [
34] 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 [
35] 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
% 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
% of the entire variation of contour error, so that
% 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.
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 [
28] using a model predictive control framework. Taking the
x-axis, for example, a desired trajectory
is parametrized with respect to a normalized, monotonically increasing path variable
, which is a vectorized form of
p. Then,
is linearized as
with respect to
using an estimated linearization point
as
The procedure for computing the optimal
(corresponding to the optimal feedrate) using the uncertainty-aware digital twin is as follows. The path variable
is maximized under monotonicity, maximum feedrate, and axis-acceleration constraints as
where
is a ones-vector,
D is a difference operator, and
and
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
and
are linear in terms of
, by showing that the only alterable feature in
, which is the last term in
(i.e.,
), is linear in
.
Let
be the matrix (lifted domain) representation of
truncated by length
. The last
rows in
can further be decomposed into two parts: its first
columns
and its last
columns
as
If
represents the last
elements of the
at past timesteps,
can be re-written as
where
is a selection matrix that picks the entry at timestep
t. Similarly, for
y-axis, the alterable term in
can be derived to be linear in terms of
, by using a similar notation as Eq. (
25), i.e.,
.
Then, the worst-case out of the
[%] variations of distribution of the contour error
, where
is a user-defined stringency, is bounded by the tolerance
as a stochastic constraint by
For the sake of brevity, the negative stochastic contour error constraint
is omitted. Then, inversion the both sides of Eq. (
26) becomes
where
is the cumulative density function of the distribution of
, which is invertible because
follows a Gaussian distribution as was shown in Eq. (
21).
Eq. (
27) can be rearranged as a linear constraint in terms of
and
using Eq. (
21), written as
where the derivation of coefficients
,
and
is described in
Appendix C. Finally, the contour error constraint is also linear with respect to the decision variable
using the relationship in Eq. (
22).
The methodology of feedrate optimization described in this section can be broadly considered as model predictive control [
39], 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:
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/
, 3 m/
and 4 m/
, and 3 FRFs are measured per each acceleration amplitude to collect a total of
FRFs per axes.
The set of FRFs
of the
x- and
y-axis of the printer are shown in
Figure 7 and
Figure 8, respectively. The uncertainties in
are then propagated to
to initialize
and
and construct
and
in the physics-informed data-driven digital twin. To validate the hypothesis that the FRF coefficients of
and
belong to Gaussian distributions, Lilliefors test [
41] is performed on every frequency
for
of
and
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
and
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
and
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
is simulated as the sum of motion-induced position
and force-induced position
, as
where
= 0.2 and
= 733 rad/s (7000 rpm) are chosen.
The butterfly trajectory [
42] with its contour of the toolpath on the
x- and
y-axis shown in
Figure 11 is selected. For the DD and PIDD methods,
= 10,
= 3,
= 10 and
= 0.01 are selected. For the DD and PIDD methods, stringency
= 95% is selected.
= 30 mm/s,
= 5 m/
, and contour error limit of
= 0.4 mm are selected for the feedrate optimization. The tolerance violation
, which will be analyzed for each method, is defined as
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
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
= 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
= 95% and 98%. It is observed that tuning
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
% of contour error distribution is within the tolerance. This issue can be mitigated by increasing
, 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. (
A3) and sub-optimal learning in
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 = 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 = 95% for numerical validation.
Figure 13.
Feedrate, acceleration, contour error, tolerance violation, and prediction error using the proposed (PIDD) methods with
= 95% (from
Figure 12 and
= 98% for numerical validation.
Figure 13.
Feedrate, acceleration, contour error, tolerance violation, and prediction error using the proposed (PIDD) methods with
= 95% (from
Figure 12 and
= 98% for numerical validation.
Table 1.
Cycle times and RMS of tolerance violation
for conservative (Cons.), physics-based (PB), data-driven (DD) and proposed (PIDD) methods in
Figure 12 and
Figure 13.
Table 1.
Cycle times and RMS of tolerance violation
for conservative (Cons.), physics-based (PB), data-driven (DD) and proposed (PIDD) methods in
Figure 12 and
Figure 13.
|
Cons. |
PB |
DD |
PIDD (=95%) |
PIDD (=98%) |
RMS of [m] |
2.2 |
6.6 |
3.0 |
1.7 |
0.8 |
Cycle time [s] |
8.89 |
4.49 |
5.17 |
6.47 |
7.17 |
Figure 1.
Diagram of intelligent feedrate optimization using an uncertainty-aware digital twin. A manufacturer provides a part tolerance and stringency (i.e., the tolerance for quality constraint violation under uncertainty). The intelligent machine leverages an uncertainty-aware digital twin to optimize feedrate while satisfying the tolerance and stringency requirements.
Figure 1.
Diagram of intelligent feedrate optimization using an uncertainty-aware digital twin. A manufacturer provides a part tolerance and stringency (i.e., the tolerance for quality constraint violation under uncertainty). The intelligent machine leverages an uncertainty-aware digital twin to optimize feedrate while satisfying the tolerance and stringency requirements.
Figure 3.
Flowchart of intelligent feedrate optimization using a deterministic digital twin with physics-based and data-driven servo models (y-axis omitted for simplicity).
Figure 3.
Flowchart of intelligent feedrate optimization using a deterministic digital twin with physics-based and data-driven servo models (y-axis omitted for simplicity).
Figure 4.
Batches and j defined on a discrete-time domain.
Figure 4.
Batches and j defined on a discrete-time domain.
Figure 7.
Frequency response functions of the Nomad 3’s x-axis showing the known uncertainty obtained under different input acceleration amplitudes.
Figure 7.
Frequency response functions of the Nomad 3’s x-axis showing the known uncertainty obtained under different input acceleration amplitudes.
Figure 8.
Frequency response functions of the Nomad 3’s y-axis showing the known uncertainty obtained under different input acceleration amplitudes.
Figure 8.
Frequency response functions of the Nomad 3’s y-axis showing the known uncertainty obtained under different input acceleration amplitudes.
Figure 9.
Upper plot: FRF coefficients and of Nomad 3 x-axis; lower plot: p-value and accept/reject results of Lilliefors test of FRF coefficients.
Figure 9.
Upper plot: FRF coefficients and of Nomad 3 x-axis; lower plot: p-value and accept/reject results of Lilliefors test of FRF coefficients.
Figure 10.
Upper plot: FRF coefficients and of Nomad 3 y-axis; lower plot: p-value and accept/reject results of Lilliefors test of FRF coefficients.
Figure 10.
Upper plot: FRF coefficients and of Nomad 3 y-axis; lower plot: p-value and accept/reject results of Lilliefors test of FRF coefficients.
Figure 11.
Desired toolpath.
Figure 11.
Desired toolpath.
Figure 15.
Frequency response functions of the Ender 3 Pro’s x-axis showing the known uncertainty obtained under different input acceleration amplitudes.
Figure 15.
Frequency response functions of the Ender 3 Pro’s x-axis showing the known uncertainty obtained under different input acceleration amplitudes.
Figure 16.
Frequency response functions of the Ender 3 Pro’s y-axis showing the known uncertainty obtained under different input acceleration amplitudes.
Figure 16.
Frequency response functions of the Ender 3 Pro’s y-axis showing the known uncertainty obtained under different input acceleration amplitudes.
Figure 17.
Commanded feedrate, acceleration, contour error, tolerance violation and prediction error using conservative (Cons.), physics-based, data-driven and proposed approaches in air-printing.
Figure 17.
Commanded feedrate, acceleration, contour error, tolerance violation and prediction error using conservative (Cons.), physics-based, data-driven and proposed approaches in air-printing.
Figure 18.
Top and side views of 3D-printed butterfly models using conservative, physics-based, data-driven and proposed approaches and their print times.
Figure 18.
Top and side views of 3D-printed butterfly models using conservative, physics-based, data-driven and proposed approaches and their print times.
Figure 21.
Machined parts and the zoomed-in images of upper left wing using conservative, physics-based (PB), data-driven (DD) and proposed (PIDD) approaches
Figure 21.
Machined parts and the zoomed-in images of upper left wing using conservative, physics-based (PB), data-driven (DD) and proposed (PIDD) approaches
Table 2.
Comparison of RMS prediction errors, cycle times and RMS of tolerance violation using conservative (Cons.), physics-based (Physics.), data-driven (Data.) and proposed methods in air-printing.
Table 2.
Comparison of RMS prediction errors, cycle times and RMS of tolerance violation using conservative (Cons.), physics-based (Physics.), data-driven (Data.) and proposed methods in air-printing.
|
Cons. |
Physics. |
Data. |
Proposed |
x Pred. Error [m] |
N/A |
37.4 |
22.1 |
18.1 |
y Pred. Error [m] |
N/A |
31.7 |
23.9 |
19.3 |
Cycle time [s] |
4.70 |
1.97 |
2.73 |
3.86 |
RMS of [m] |
1.8 |
5.5 |
3.9 |
1.9 |