# Data-driven Reference Trajectory Optimization for Precision Motion Systems

Samuel Balula, Dominic Liao-McPherson, Alisa Rupenyan, and John Lygeros The
authors are with the ETH Zürich Automatic Control Laboratory, Physikstrasse 3,
8092 Zürich, Switzerland. E-mail: {sbalula, dliaomc, ralisa,
<EMAIL_ADDRESS>S. Balula and A. Rupenyan are also with Inspire AG. This
project has been funded by the Swiss Innovation Agency (Innosuisse, Grant
Number 46716) and the Swiss National Science Foundation through NCCR
Automation, a National Centre of Competence in Research (Grant Number 180545).

###### Abstract

We propose a data-driven optimization-based pre-compensation method to improve
the contour tracking performance of precision motion stages by modifying the
reference trajectory and without modifying any built-in low-level controllers.
The position of the precision motion stage is predicted with data-driven
models, a linear low-fidelity model is used to optimize traversal time, by
changing the path velocity and acceleration profiles then a non-linear high-
fidelity model is used to refine the previously found time-optimal solution.
We experimentally demonstrate that the proposed method is capable of
simultaneously improving the productivity and accuracy of a high precision
motion stage. Given the data-based nature of the models, the proposed method
can easily be adapted to a wide family of precision motion systems.

###### Index Terms:

Precision Engineering, Motion control, Pre-compensation, Data-driven modeling

##  I Introduction

A precision motion stage (PMS) is a robotic setup designed to precisely
position an end-effector, typically in the micro- to nanometer range. PMS are
used in lithography [1], micro-assembly [2], precision metrology [3], and
precision engineering [4], and many other applications [5, 6]. They are often
used as components in machine tools for the production of parts, such as metal
laser cutting or grinding processes where the goal is to trace a desired
target geometry as quickly and accurately as possible.

The accuracy of PMS systems is critical as the parts they produce must satisfy
tight engineering tolerances. At the same time, economic considerations
dictate rapid production of parts to boost productivity. Unfortunately, the
physical inertia of PMS systems dictates a trade-off between productivity and
accuracy and, in practice, the quality of the control system is critical for
managing this trade-off. Typical industrial precision systems achieve high-
precision using well-tuned “classical” feedback controllers, such as
Proportional-Integral-Derivative (PID), to track a given target trajectory.
These low-level controllers are designed by the manufacturer and users are
typically not able to modify/tune the control algorithm.

In this paper, we propose an optimization-based reference (input trajectory)
design (pre-compensation) methodology for simultaneously improving the
precision and productivity of PMS by leveraging the widespread availability of
system data. Our methodology is guided by the following requirements:  (i) it
does not require modifications to the existing hardware and low-level
controllers, making it suitable for application in production machines;  (ii)
it does not require a physics-based model of the machine, relying exclusively
on experimental data;  (iii) it can be applied to individual machines with
diverse technical characteristics; and  (iv) it can produce good results for
previously unseen parts at the first attempt.

These requirements are motivated by the requirements of “Industry 4.0” where
machines are no longer standalone setups and are instead connected in a
network, making data more easily available [7]. This data can be exploited to
improve performance; for example, it can be leveraged to create data-driven
models, avoiding the labour intensive process of developing and maintaining
first principles based ones [8]. Future manufacturing is also expected to be
more distributed and personalized (e.g., “Manufacturing as a service”) [9],
leading to smaller volumes for any given part, possibly produced across a pool
of non-homogeneous machines. This makes it more important to exploit data to
ensure that machines can execute designs with minimal calibration.

Several methods for improving PMS performance have been proposed in
literature. One method is to re-designing the feedback controller. A common
approach is model predictive contouring control (MPCC) [10], wherein a model
of contouring error is minimized in a receding horizon fashion, drawing on
methods from autonomous racing [11]. However these methods have high online
computational loads, rely on accurate process models, and require modifying
the inner-loop controllers.

Another common approach in manufacturing and precision motion systems is
learning-based control [12, 13]. These methods assume that the processes is
repeatable and adjust the input trajectories, using the error from the
previous repetition as feedback. The input trajectories can be modified
between iterations (Iterative Learning Control, and Run-to-Run control), or
periods (Repetitive Control). These methods can exploit models [14], include
constraints [15] , and are effective for repeatable processes. However, the
computed trajectory is specific for a particular part and machine, with
multiple trials required to achieve the desired accuracy. While this may be
acceptable in large production runs, it may be inefficient for smaller runs
e.g., in ”manufacturing as a service” applications.

The final major class of methods is pre-compensation, which aim to modify the
input trajectory offline to reduce the contouring error for arbitrary parts
[16]. A common approach is to drive a model of the contouring error to zero by
manipulating the input trajectory using a classical controller (e.g., a PID)
[17]. These methods cannot include look-ahead or constraints, motivating the
use of offline MPCC algorithms based on physics-based nonlinear [18],
identified linear models [19] or combined with a reference governor [20]. In
[21] the compensation is computed online using previous measurements and a
linear model. In [22] the contour error is predicted with an Artificial Neural
Network (ANN), and the pre-compensation computed with reinforcement learning.
This approach does not allow constraints to be included and uses only zeroth-
order information about the ANN model, despite the ready availability of ANN
derivatives. In [23, 24] a simultaneous feedrate optimization and error pre-
compensation method is proposed for a system described by a linear model. This
approach is analogous to the first stage optimization described in IV-A.

Therefore, to the best of our knowledge, there is a gap in the existing
literature on pre-compensation methods for improving PMS
performance/productivity that are data-driven, capable of generalizing to
novel geometries, and applicable to a broad variety of machines that can only
be accuratelly described by a nonlinear model.

In this paper we tackle the challenge of simultaneously improving the
performance/productivity of a PMS using only data from the system. Our main
contribution is a novel pre-compensation method based on data-driven modelling
and trajectory optimisation. Our approach does not require modification of the
design or control software of the system (i) and is fully data-driven (ii).
Moreover, it enables adaptation to changing conditions or system degradation
through retraining the underlying data-driven model (iii). Finally, the
approach is independent of the target geometry and significantly shifts the
precision/productivity trade-off curve111Relative to a method without pre-
compensation, our proposed method is 30-60% more accurate for new geometries
for same trajectory time budget or alternatively, trajectories can be traced
in 25-50% of the time for the same contouring accuracy. (iv).

We build upon our previous work [25] based on linear system identification.
Here we improve and extend this approach by incorporating a nonlinear ANN-
based model of system performance, which offers two orders of magnitude higher
precision, a refined optimization approach, and experimental validation of the
simulation results.

This paper is organized as follows: in Section II we precisely define the
problem tackled, while in Section III we describe the modelling structure, the
model fitting strategy, and the accuracy achieved. In Section IV we detail the
input trajectory design strategy proposed, specifying the optimization
problems. In Section V we show experimental validation for the method with a
set of test cases.

##  II Problem Formulation

Consider a tooltip, e.g., the laser in a laser cutting machine, which
traverses a two-dimensional (2D) workspace $\mathcal{W}\subset\mathbb{R}^{2}$.
The position of the tooltip as a function of time is denoted by
$\gamma(t)=[x(t)~{}~{}y(t)]^{\top}.$ The ultimate goal in precision motion
planning is to have the tooltip trace a spatial target geometry
$\Xi:[0,S]\to\mathcal{W}\,,$ where $S$ is the path length. The target geometry
$\Xi$ is assumed to be a continuous function parameterized by the path
progress variable $s$. Typical design specifications require the target
geometry to be traced with a precision of tens of micrometers; for example, in
industrial laser cutting systems a typical precision specification is
$20~{}\mathrm{\mu m}$

Most industrial precision motion systems come with a built-in low-level
controller designed by the manufacturer to track an input trajectory
$r:[0,T]\to\mathcal{W}$, where $T>0$ is the time necessary for the tooltip to
complete tracing the target geometry. When $r$ is given as an input to the
low-level controller, its actions and the dynamics of the machine give rise to
an output trajectory $\gamma:[0,T]\to\mathcal{W}$; note that for simplicity we
assume that both the input and the output trajectory lie in the workspace
$\mathcal{W}$. We denote the mapping between the commanded input trajectory
$r$ and the resulting output trajectory $\gamma$ by

$\gamma=F(r)+e,$ (1)

where the function $F:\mathcal{W}^{[0,T]}\to\mathcal{W}^{[0,T]}$ encapsulates
the closed-loop dynamics of the machine. The noise $e$ is zero-mean and
quantifies the repeatibility of the system, a typical standard deviation is
$\approx 2\,\mathrm{\mu m}$. Ideally, the low-level controller is well
designed so that $\gamma\approx r$ and $F$ is approximately the identity
function.

We are interested in the following inverse problem: Given the target geometry
$\Xi$, how do we generate an input trajectory $r$ for the machine, so that the
resulting output $\gamma=F(r)$ matches $\Xi$ as closely as possible. In
addition to the usual difficulties associated with inverse problems, for
precision motion systems the function $F$ is often poorly known, partly due to
lack of information about the low-level controller.

In standard practice, input trajectories are generated by traversing the path
at a constant speed $v>0$, i.e. setting $r(t)=\Xi(vt)$. However, due to sensor
noise, actuator limits, the inertia of the machine, and other factors we
observe experimentally a significant error between $\gamma$ and $r$,
especially when $v$ is high. This induces a natural trade-off between speed
and precision that is characteristic of precision motion systems. Our aim is
to improve this trade-off by learning the function $F$ through preliminary
experiments, then inverting it using optimization.

Figure 1: A low-level controller tracks the input trajectory $r$, producing
the output trajectory $\gamma$. The deviation $d(s,\gamma)$ measures the
distance from $\gamma$ to the target geometry $\Xi$ as a function of the path
progress $s$. We desire that the deviation is smaller than the tolerance
$\mu(s)$ at all points.

Formally, let $d:[0,S]\times\mathcal{W}^{[0,T]}\rightarrow\mathbb{R}$ be the
deviation function, defined as the projection of the output into the target
geometry $d(s,\gamma)=\min_{t\in[0,T]}\|\Xi(s)-\gamma(t)\|_{2}.$ Our goal is
then to solve the optimization problem [c1] rT(γ)+ λ∫_0^S d(s,γ) ds γ= F(r)
d(s,γ) ≤ μ(s), ∀s ∈[0,S] , where $T(\gamma)$ is the duration of the output
trajectory, $\lambda>0$ governs the trade-off between speed and accuracy, and
$\mu(s):[0,S]\rightarrow\mathbb{R}$ is a tolerance bound around the target
geometry; for the typical industrial laser cutting systems mentioned above,
one would use $\mu(s)=20~{}\mathrm{\mu m},~{}\forall s\in[0,S]$. The different
quantities are illustrated in Figure 1.

The optimization problem (1) is stated in a trajectory space and is not
computationally tractable as written. In the following sections we present and
experimentally validate a practical methodology for approximately solving (1)
using only data gathered from the process. Because the methodology is data-
driven and noninvasive (in the sense that it does not require any knowledge of
the underlying controllers but only the ability to run experiments on the
machine) it is applicable to a wide range of practical industrial
applications.

###  II-A Experimental Setup

(a) Experimental apparatus

(b) Letters test case

(c) airfoil test case

Figure 2: The experimental setup shown in panel 2(a) features two degrees of
freedom, $x$ and $y$, driven by linear actuators. The end effector is mounted
on the $y$ axis, which in turn is mounted on top of the $x$ axis. In panel
2(b) the contour of the letters in black is the Letters test case, and the
outline in panel 2(c) the airfoil test case, labeled with the leading (L) and
trailing (T) edges locations.

While our proposed methodology is applicable for general PMS, we use a
specific experimental setup as a running example throughout the paper and to
experimentally validate our approach. The experimental setup is shown in
Figure 2(a), and consists of an ETEL DXL-LM325 two stage motion machine. The
machine is instrumented with quadrature encoders that measure the position of
both axes and is actuated with ETEL LMS15-100 linear motors.

The system is equipped with a feedback controller that adjusts the motor
voltages to track a supplied input trajectory $r$. The controller is
implemented on a Speedgoat rapid prototyping unit using MATLAB/Simulink
software. For the purposes of this paper we consider this loop as a black-box
and consider only the desired input trajectory $r$ as an input. The machine
enforces limits on the acceleration, velocity, and position of the input
trajectory to prevent accidental damage. The main characteristics of the
machine are collected in Table I.

Property | Value | Unit
---|---|---
Maximum acceleration | $40$ | $\mathrm{m\,s^{-2}}$
Maximum velocity | $1.5$ | $\mathrm{m\,s^{-1}}$
Working area | $[-0.19,0.19]^{2}$ | $\mathrm{m^{2}}$
Low-level control loop frequency | $10$ | $\mathrm{kHz}$
Repeatability (standard deviation) | $2$ | $\mathrm{\mu m}$
TABLE I: Software-limited values and operational specifications of the
experimental setup.

###  II-B Test cases

We selected two test cases, a series of letters and an airfoil, as target
geometries to validate our proposed methodology. These geometries contain a
rich set of features that are representative of industrial applications. To
avoid “cherry picking”, and to demonstrate generalization to new geometries,
neither test case was used to generate training data or during the development
of the method.

The letters test case, as can be seen in Figure 2(b), consists of four letters
(“u”,“r”,“c”, and “h”) from the ETH Zürich logo. These letters were selected
due to the rich set of relevant and challenging features they contain,
including sharp corners, straight segments, and curves with various
curvatures. The limiting contour of each letter is used as the target
geometry.

The airfoil test case, shown in Figure 2(c), is based on a high-lift Eppler
e344 airfoil[26]. Geometries of this kind are typically laser cut out of e.g.,
balsa wood, and are used in the primary structural sub-assembly of an
aircraft. The main challenging features of the airfoil are a cut-out at the
leading edge (L) used to attach a structural spar, and the trailing edge (T)
whose geometry has a strong effect on the aerodynamics of the aircraft.

###  II-C Architecture overview

Figure 3: In the first stage optimization (6) a minimum time trajectory for
tracking the target path $\Xi(s)$ is computed. The speed and acceleration
profile from the first stage is used to sample the target geometry at constant
sample rate of $400~{}\mathrm{Hz}$, creating a surrogate target discrete time
trajectory $\bar{\Xi}$. In the second stage optimization (7) nonlinear model
is used to find the best input trajectory $r$ to track $\bar{\Xi}$. The
maximum allowable acceleration $\mathtt{a_{max}}$ governs the
accuracy/trajectory time trade-off. The optimized input trajectory $r$ is
given to the experimental apparatus.

We approach the trajectory design problem in several stages. First, we design
experiments to generate informative input-output data for the unknown mapping
$F$ and use the resulting data to identify two models for the system, a low-
fidelity linear one and a high-fidelity one based on an ANN. Throughout, we
incorporate machine learning best practices, such as independent training and
validation data sets, normalization, and diverse training data to avoid over
fitting and to ensure that the resulting models are able to generalize to
trajectories outside the training data set, enabling optimization of new
geometries.

Using these models, we propose a two-stage optimization-based approach for
offline optimization of references as schematized in Figure 3. In the first
stage, we compute an input trajectory by solving a contouring control problem
using the linear model of the system. This first stage solution yields a fast
trajectory that respects the problem constraints, by reducing the speed in
intricate features of the path, such as sharp corners, and accelerating
through long smooth segments. The output of the first-stage is used as the
initial condition for the second stage, which employs the high-fidelity neural
network model to correct errors from the first stage. The resulting input
trajectory is then given to the low-level controller as an input trajectory.

##  III Data driven modelling

As a first approximation, we model the system using a continuous-time linear
model which captures the dominant dynamics. The structure of this model allows
for efficient computations of a time-optimal solution in the first stage, but
lacks the high-precision required by the application. This shortcoming is
alleviated by a second high-fidelity neural network model, used to refine the
initial solution.

###  III-A Low-fidelity linear model

The closed-loop system maps trajectories to trajectories and is of the form
$\gamma=F(r)$. We approximate this mapping using a linear time invariant (LTI)
state space model of the following form

$\displaystyle\dot{z}(t)$ $\displaystyle=Az(t)+Br(t)\mathrm{,}$ (2a)
$\displaystyle\gamma(t)$ $\displaystyle=Cz(t)+Dr(t)\mathrm{.}$ (2b)

The model (2) uses a non-physical hidden state $z\in\mathbb{R}^{4}$; the
matrices $A$, $B$, $C$, and $D$ are identified from experimental data using
the `MATLAB` function `n4sid` from the system identification toolbox. The
dimension of the state was chosen empirically, noting that increasing the
model order further has diminishing returns.

###  III-B High-fidelity Nonlinear Model

To obtain a higher precision model, the 2D-positioning stage is modeled using
two independent causal ANN, one for each axis. Unlike (III-A) the nonlinear
model operates in discrete time, at a sample frequency of $400\,\mathrm{Hz}$
and it has the following form

$\gamma_{i}=f_{nlm}(r_{i-h},\dots,r_{i-1}),$ (3)

where $\gamma$ is the output, $f_{nlm}$ is the nonlinear function defined by
the ANN, $r$ is the input trajectory given to the closed loop system, and $h$
is the length of the input history considered for the prediction. The two ANN
are structurally identical, but are trained separately, leading to different
weights. Each takes as input the most recent $500\mathrm{ms}$ of the input
trajectory for the corresponding axis, subsampled at $400\mathrm{Hz}$, leading
to $200$ inputs. Each network has a single output, namely the predicted
position at the subsequent time step. The networks have $8$ fully connected
hidden layers of LeakyReLU [27] activation functions with a slope of $0.01$,
and a triangular structure. The length of the time history was selected to be
approximately $25$ times the timescale of the slowest mode of the identified
linear model. The sample rate was selected by analyzing how much distance can
be covered at the highest expected speeds, and comparing it to the desired
range of precisions. The number of hidden layers and their dimensions were
selected empirically, to achieve good performance with the lowest possible
complexity. The choice of feed-forward ANN was made to ensure that the model
prediction accuracy does not degrade over time due to compounding errors. We
also note that by training separate ANN for each axis we are implicitly
ignoring coupling effects between the two axes; this choice is supported by
experimental data that suggests coupling effects are negligible, as will be
demonstrated in Section III-D.

Figure 4: Sample of the data collected to train and test the nonlinear neural
network model. Six trajectories are shown, five randomly generated and one
regular shape (a spiral). The trajectories are contained in the workspace
$\mathcal{W}$, shown in yellow background.

###  III-C Training data generation

Choosing input trajectories that cover a wide range of operational regimes is
essential for generating informative data to train models that can generalize
to previously unseen trajectories.

####  III-C1 Linear Model

For the identification of the linear model we opted for randomly generated
trajectories, yielding $6\times 10^{6}$ points in total. By randomly changing
the velocity, acceleration, and jerk of the input trajectory we made sure that
the trajectories represent features seen in real parts. This included, for
example, curves of different radii, traversed at different speeds and
appearing in different locations within the workspace, as well as sharp
transitions in acceleration and jerk. The trajectories were not designed for a
specific part, but for a particular set of constraints on velocity and
acceleration.

####  III-C2 Nonlinear Model

Having sufficient high-quality data is necessary to ensure that the ANN is
accurate and generalizes well. We started by training a first generation ANN
with the random trajectories used to fit the linear model, and a set of varied
regular geometries including circles, polygons, and spirals. Each shape was
traversed at varied speeds, constant for each trial. To enrich the data set
further, the trained first generation ANN was used in the optimization
methodology of Section IV applied to the same regular geometries. This
resulted in optimized input trajectories and, after applying these to the
system, additional experimental data. This data was used to train the ANN
further obtaining a second generation ANN. This training-optimization-
retraining cycle can be executed recursively until the model prediction
accuracy is acceptable. These data sets lead to a targeted reduction in over
fitting in areas that are favored by the optimization. Overall, our strategy
is to include sufficient random (exploratory) and structured (exploitative)
data to ensure the ANN is both accurate and generalizes well in the areas of
interest.

A representative training data sample is shown in Figure 4 In section III-D we
present the prediction quality of a third generation ANN to generate the
results of Section V. The data was divided $80\%/20\%$ for training and
testing purposes. This results in a ratio of $\approx 34$ data points per
parameter of the ANN. About $60\%$ of the data comes from the randomly
generated trajectories, and the remaining from regular shapes both before and
after optimization. The experimental data is divided in segments of
$500~{}\mathrm{ms}$, subsampled to achieve an effective sample rate of
$400\mathrm{Hz}$, and scaled. For the loss function we use the mean square of
the prediction error. We use batches of $16\mathrm{k}$ segments which take
maximum advantage of the GPU computation power. Adam [28] was used for the
optimizer, with decaying learning rate. Modeling and training were done in
PyTorch [29] with NVIDIA CUDA support.

###  III-D Model Validation

|  | Training / | $\mathtt{a_{max}}$ | $\mathtt{a_{max}}$
---|---|---|---|---
|  | Validation | $1.0~{}\mathrm{m\,s^{-2}}$ | $3.0~{}\mathrm{m\,s^{-2}}$
Model | Axis | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
– | – | $\mathrm{\mu m}$ | $\mathrm{\mu m}$ | $\mathrm{\mu m}$ | $\mathrm{\mu m}$ | $\mathrm{\mu m}$ | $\mathrm{\mu m}$
Linear | x | 3.0 | 48.0 | 4.7 | 30.1 | 5.4 | 57.2
| y | -13.2 | 156.3 | 0.0 | 108.5 | 0.5 | 229.4
Non-linear | x | 0.0 | 3.2 | -0.3 | 8.2 | 1.0 | 11.1
| y | 0.0 | 6.0 | -2.0 | 13.4 | -2.3 | 18.7
TABLE II:  Prediction error mean $\mu$ and standard deviation $\sigma$ of the
linear and non-linear models. Both models are evaluated in the Letters test
case after optimization, with different maximum acceleration values
$\mathtt{a_{max}}$. The linear model is also evaluated on its training data
set, and the non-linear model on its validation data set. Figure 5: Prediction
error of the nonlinear model as a function of the maximum acceleration
$\mathtt{a_{max}}$ of the input trajectory. The error is evaluated for the
letters and airfoil test cases. Above $5~{}\mathrm{m\,s^{-2}}$ the training
data do not contain enough information to train the model properly.

The quality of the model predictions is evaluated over different data sets and
acceleration limits. The results are summarized in Table II.

The linear model training data set includes trajectories with a range of
accelerations up to $5~{}\mathrm{m\,s^{-2}}$. However, due to its simple
structure, the linear model cannot overfit the data as it does not have the
representative power to provide accurate predictions everywhere in its
training set. We observed that it performs better when the acceleration is
restricted to a lower value, even for trajectories not seen during the fit,
such as the letters test case. The linear model accuracy is insufficient for
the demanding requirements of the application.

Regarding the nonlinear model, the prediction accuracy is best for data with
the same characteristics as the one used for the training. When presented with
data coming from shapes which have not been used for the training of the
model, the prediction accuracy degrades. However, the accuracy is still
adequate for the application, as the standard deviations remain below
$20~{}\mathrm{\mu m}$ for each axis for accelerations below
$3~{}\mathrm{m\,s^{-2}}$.

Figure 5 shows how the standard deviation of the prediction evolves with the
maximum acceleration $\mathtt{a_{max}}$ that was set as a constraint in the
optimization of Section IV. For acceleration values where there is sufficient
training data, the model prediction error increases gradually with the
acceleration. If we try to use the model outside of the acceleration values it
was trained for, the predictions rapidly deteriorate.

##  IV Trajectory Optimization

With the the data-driven models in hand, we now discuss the two-stage
optimisation architecture used to compute input trajectories. The decision to
divide the optimization problem in two stages is motivated by the need to
include the duration $T$ of the trajectory in a computationally tractable
manner. This is accomplished by using a variable time-discretization in the
first stage, enabling inclusion of $T$ as a decision variable. The second
stage then fixes the duration $T$ and refines the accuracy of the first-stage
solution using the high-fidelity neural network model.

Our experiments suggest that treating the time discretization directly with
the nonlinear model is intractable. A key difficulty is that real machines
typically generate noisy output measurements in a sampled, time series format.
These cannot be readily used to train the high-fidelity continuous time model
needed for a computationally tractable variable-time discretization, due the
challenges associated with numerical differentiation of noisy data. Trying to
solve the inverse problem in one shot using the resulting model leads to low
quality solutions; often the solver fails to even return a feasible solution
in a reasonable amount of time.

On the other hand, skipping the second stage and inverting using only the
linear model is computationally tractable, but typically leads to low quality
solutions; this observation is supported by the model prediction errors shown
in Table II. Training a fixed sampling time nonlinear model and seeding with
the solution of the first stage optimization substantially reduces the overall
computation time and improves the quality of the solutions.

A key advantage of our methodology is that, unlike iterative learning control,
once the models are trained and we are given a target geometry, the input
trajectory is computed off-line, without the need for additional experiments.
This is especially important for small batch manufacturing where it is crucial
to minimize failed parts. However we can still use the information collected
from experimental runs to improve the model.

###  IV-A First-stage Optimization Formulation

The purpose of the first stage is to fix the duration of the trajectory, $T$,
by trading off speed and precision. We employ a contouring control approach
with a fixed spatial discretization. This involves sampling the target
geometry at $N$ equally spaced points in space
$\\{\Xi_{k}\doteq\Xi(s_{k}):k\in\\{0,1,\dots,N-1\\}\\}$, where
$s_{k}=k\frac{S}{N-1}=k\Delta s.$

We use the identified linear model (2) to approximate the process dynamics. As
our approach involves a fixed spatial discretization $\Delta s$, the time-
discretization must be variable. Starting with $t_{0}=0$, we denote by $t_{k}$
as the time at which we would like the trajectory to reach the point $\Xi_{k}$
and consider the non-uniform time steps
$\Delta\tau_{k}=t_{k+1}-t_{k},k=\\{0,\dots,N-1\\}$ as decision variables; note
that the time $T=t_{N}=\sum_{k=0}^{N-1}\Delta\tau_{k}$ at which the trajectory
is completed is implicitly also a decision variable. We then obtain a discrete
time model by applying the $4$th order Runge-Kutta (RK4) [30] time-marching
method to (2)

$\displaystyle z_{k+1}$
$\displaystyle=f_{\mathrm{lm}}(z_{k},r_{k},\Delta\tau_{k})$ (4a)
$\displaystyle\gamma_{k}$ $\displaystyle=Cz_{k}+Dr_{k},$ (4b)

where $f_{\mathrm{lm}}$ is the RK4 function, and $\Delta\tau$ will be treated
as decision variables.

To evaluate the contouring error $d_{k}$ we start with the distance
$\nu_{k}=\gamma_{k}-\Xi_{k}$ and introduce a local coordinate frame as shown
in Figure 6. The linear transformation

$T_{k}=\begin{bmatrix}\cos{\alpha_{k}}&\sin{\alpha_{k}}&\Xi_{k,x}\\\
-\sin{\alpha_{k}}&\cos{\alpha_{k}}&\Xi_{k,y}\\\ 0&0&1\end{bmatrix},$ (5)

where $\alpha_{k}$ is the rotation angle of the local frame, encodes a
translation and rotation between the local coordinate frame $k$ and the global
coordinate frame. With

$\begin{bmatrix}\gamma_{k}\\\ 1\end{bmatrix}=T_{k}\begin{bmatrix}\nu_{k}\\\
1\end{bmatrix},$ (6)

the output can be transformed from local coordinates $\nu_{k}$ to global
coordinates $\gamma_{k}$. Since time discretization is variable, we can always
ensure that the $x^{l}$ component of $\nu_{k}$ is zero by adding a constraint
to the optimization problem. The contouring error is then simply the $y^{l}$
component of $\nu_{k}$. The tolerance bound condition in (1) can then also be
computed for the corresponding discretization point $\mu_{k}=\mu(s_{k})$.

Figure 6: For every point $\Xi_{k}$ of the target geometry $\Xi(s)$ a local
coordinate frame is placed. $x^{l}$ is tangent to the target geometry at
$\Xi_{k}$ and pointing in the direction of increasing $s$ whereas $y^{l}$ is
orthogonal to $x^{l}$ pointing to the port side. $\alpha_{k}$ is the angle
between the global $x$ axis and $x^{l}$. The maximum allowed deviation
$\mu_{k}$ is computed from $\mu(s)$. $\Delta s$ is the path distance between
points $\Xi_{k}$ and $\Xi_{k+1}$.

The aim of the optimization problem is to minimize the total time needed to
complete the target geometry subject to the tolerance bounds, velocity and
acceleration limits. To this basic formulation, we introduce an additional
term for the input to promote smoothness; this is used to suppress high
frequency content on the input trajectory, which is filtered out by the linear
model. This results in the following optimization problem Δτ, r, γ, z, ν
∑_k=1^N-1 Δτ_k + 1N-2 ∑_k=2^N-1 ∥∂^2 r_k∥_2^2 z_k+1= f_lm(z_k, r_k, Δτ_k), k =
1, …, N-1 γ_k= C z_k + D r_k, k = 1, …, N Δτ_k ≥ 0, k = 1, …, N γ_k ∈ W, k =
1, …, N [γk1 ]= T_k [νk1 ], k = 1, …, N ν_k_1 = 0, k = 1, …, N — ν_k_2 — ≤
μ_k, k = 1, …, N ∥ ∂γ_k ∥_∞≤ v_max, k = 1, …, N-1 ∥ ∂^2 γ_k ∥_∞≤ a_max, k = 2,
…, N-1 , where $N$ is the number of discretization points, selected based on
the total path length and tolerance bound,
$\Delta\tau=\\{\Delta\tau_{k}\\}_{k=1}^{N-1}\subseteq\mathbb{R}$,
$\partial\gamma_{k}=\frac{\gamma_{k+1}-\gamma_{k}}{\Delta\tau_{k}}\text{ and
}\partial^{2}\gamma_{k}=\frac{\gamma_{k+1}-2\gamma_{k}+\gamma_{k-1}}{\Delta\tau_{k}\Delta\tau_{k+1}}$
are discrete derivative operators,
$\nu=\\{\nu_{k}\\}_{k=1}^{N}\subseteq\mathbb{R}^{2}$ is the output trajectory
in local coordinates,
$\gamma=\\{\gamma_{k}\\}_{k=1}^{N}\subseteq\mathbb{R}^{2}$ is the output
trajectory in global coordinates,
$r=\\{r_{k}\\}_{k=1}^{N-1}\subseteq\mathbb{R}^{2}$ is the input trajectory,
$f_{\mathrm{lm}}$ is the RK4 approximation of the linear model (4), with
identified matrices $C$ and $D$, $T_{k}$ is the transformation matrix between
global and local coordinates in (5), $\mu_{k}$ is the discretized maximum
allowed deviation from the target geometry, $\mathcal{W}$ are the limits of
the working space of the device, $\mathtt{v_{max}}$ is the maximum allowable
speed, and $\mathtt{a_{max}}$ is the acceleration limit.

The cost (IV-A) contains two terms. The first penalizes the total time to
perform the trajectory while the second penalizes the acceleration of the
input, promoting fast trajectories with a smooth input, (IV-A) and (IV-A)
imposes the RK4 discretization of the linear dynamics (4), (IV-A) requires the
variable time step to be non-negative, (IV-A) imposes that the output stays
within the working space of the device, (IV-A) the transformation between
local and global coordinates in (6), (IV-A) requires that the $x_{k}^{l}$
component of $\nu_{k}$ is zero, in which case the $y_{k}^{l}$ component is the
deviation, (IV-A) bounds the $y_{k}^{l}$ component of $\nu_{k}$ to be less
than the maximum allowed deviation, (IV-A) and (IV-A) impose component-wise
limits on velocity and acceleration of the output. This formulation is
analogous to [23] where the feedrate and contour error are optimized using a
linear model.

###  IV-B Second-stage Optimization Formulation

Experiments applying the trajectory $r$ generated by (6) directly to the
machine suggest that this generally results in large deviations from the
target geometry, as can be seen in Figure 10, which we attribute to the low
fidelity of the linear model. To improve performance, we refine the trajectory
generated by the first stage using the high-fidelity ANN model (3). In
principle, it is possible to skip the first stage and directly embed the ANN
model in an optimization problem similar to (6). However numerical experiments
suggest that this often leads to poorer performance even compared to the
linear model, as the solver tends to converge to poor local minima, or fails
to return a feasible solution. Here we take advantage of the output of (6) to
fix the time discretization and further improve precision through a second
optimization problem based on the nonlinear model. In practice this leads to a
more reliable and scalable overall method.

Given the structure of the nonlinear model, which is trained using time series
data, we sample the target geometry at $M$ equally spaced points in time

$\bar{\Xi}_{i}\doteq\Xi(s(t_{i})):i\in\\{0,1,\dots,M\\}$ (7)

$t_{i+1}=t_{i}+\Delta t$, where $\Delta t$ is the fixed sample rate of the
time series data used to train the nonlinear model. Given this fixed
discretization, the target geometry of the second-stage optimization problem
is simply to minimize the deviation from the target trajectory while
satisfying the tolerance, speed, and acceleration bounds. This leads to the
optimization problem [2]r, γ ∑_i=1^M ∥γ_i - ¯Ξ_i∥^2_2 γ_i = f_nlm(r_i-h, …,
r_i-1), i=1,…,M — γ_i - ¯Ξ_i — ≤ μ_i i=1,…,M γ_i ∈ W i=1,…,M — ∂ γ_i — ≤ v_max
i=1,…,M-1 — ∂^2 γ_i — ≤ a_max i=2,…,M-1 , where $M$ is the number of
discretization points,
$\gamma=\\{\gamma_{i}\\}_{i=1}^{M}\subseteq\mathbb{R}^{2}$ is the output,
$r=\\{r_{i}\\}_{i=1}^{M}\subseteq\mathbb{R}^{2}$ is the input trajectory,
$\bar{\Xi}=\\{\bar{\Xi}_{i}\\}_{i=1}^{M}\subseteq\mathbb{R}^{2}$ is the
sampled target geometry, $f_{\mathrm{nlm}}$ is the nonlinear ANN model (3),
$\mu_{i}$ is the tolerance bound evaluated for each discretization point,
$\mathcal{W}$ are the limits of the working space of the device,
$\mathtt{v_{max}}$ is the maximum allowable speed, and $\mathtt{a_{max}}$ is
the acceleration limit.

The cost (IV-B) penalizes deviations of the output from the target geometry.
The constraint (IV-B) imposes the nonlinear system dynamics modeled by the
ANN, (IV-B) constrains the output to be within the tolerance bound for each
discretization point, While (IV-B), (IV-B) and (IV-B) ensure that the output,
its velocity and acceleration stay within working space and operational limits
of the device.

In practice, to ensure feasibility, the tolerance constraint is replaced with
an exact $L_{1}$ penalty which is implemented using slack variables [31]. For
longer trajectories, due to memory constraints, the optimization problem can
be solved with a receding horizon strategy. In this case we use an horizon of
$11$ steps. The maximum acceleration limit $\mathtt{a_{max}}$ is the main
parameter determining the total trajectory time, tuning the solution along the
accuracy/trajectory time trade-off curve.

###  IV-C Implementation Details

The optimization problems defined in (6) and (7) are both nonlinear and are
solved using IPOPT [32], with the JuMP [33], and Casadi [34] interfaces in
Julia and Python respectively. The first stage optimization problem (6) is
initialized by taking as $r$ the target geometry traversed at a constant
speed. The solution of (6) is then used to initialize the second-stage (7).

The computer used throughout this paper runs Arch Linux, with Linux kernel
version 5.15, and it is equipped with a `GeForce RTX 2080 Ti` GPU with
$12~{}\mathrm{GB}$ of dedicated memory, an `Intel(R) Core(TM) i9-9900K` CPU
<EMAIL_ADDRESS>and $64~{}\mathrm{GB}$ of RAM.

For the problem sizes solved in this paper, the first-stage optimization takes
approximately 1 minute to complete ($N\approx 10^{3}$), and the second-stage
takes approximately between 1 to 4 hours depending on the trajectory time
($M\in\\{200,\dots,1000\\}$). Lower $\mathtt{a_{max}}$ results in slower
trajectories which take more time to compute due to the increased number of
function evaluations required. Training of the ANN model takes approximately
$24~{}\mathrm{h}$ per axis.

##  V Experimental Results

In this section, we apply our proposed data-driven optimization methodology
from Section IV to the experimental apparatus in II-A, using the test cases in
II-B. Our experiments demonstrate that the methodology is capable of improving
system performance relative to the baseline – the trajectory obtained with the
low-level controller following a not optimized input trajectory – by tracing
desired geometries faster and more precisely.

###  V-A Individual trajectories

We first focus on the letter “r” as a case study and consider a scenario with
an acceleration limit of $1~{}\mathrm{m\,s^{-2}}$; the individual results for
the other letters in the test case are comparable and will be discussed
collectively in Section V-B. Applying our method results in an optimized input
trajectory with a total time of $T=0.807~{}\mathrm{s}$. The experimental
result of this new input trajectory is compared with a baseline performed in
the same total time. For this the original, non-optimized shape (in this case
the letter “r”) is sampled at a constant progress speed
$r(t)=\Xi\left(t\frac{S}{T}\right)\mathrm{,}$ where $S$ is the total path
length, which for a closed shape corresponds to the perimeter.

Figure 7: Letter “r” of the ETH Zürich Logo. Optimized with
$\mathtt{a_{max}}=1~{}\mathrm{m\,s^{-2}}$, including a detail view of station
B. Note how the optimized input trajectory significantly deviates from the
target geometry near the corner in B to compensate for the dynamics of the
machine. This leads to significantly smaller error between the output and the
target geometry compared to using the state-of-practice (non-optimized) input
trajectory. Figure 8: A close-up of the trailing edge of the airfoil test case
with $\mathtt{a_{max}}=2~{}\mathrm{m\,s^{-2}}$. The optimized trajectory is
able to track the desired geometry more precisely than the baseline, though
still not within the $20~{}\mathrm{\mu m}$ tolerance band. The optimized input
trajectory is not constrained by the tolerance band.

(a) Letter “r” of the ETH Zürich logo optimized with
$\mathtt{a_{max}}=1~{}\mathrm{m\,s^{-2}}$. The letters A-G refer to the points
labeled in Figure 7.

(b) Airfoil test case, optimized with
$\mathtt{a_{max}}=2~{}\mathrm{m\,s^{-2}}$. The letters L and T refer to the
leading and trailing edges labeled in Figure 2(c).

Figure 9: Top panels show the deviation from the target geometry for the
optimized and non-optimized output, as well as the optimized input trajectory.
The non-optimized input trajectory by construction does not deviate from the
target geometry. Bottom panels shows the velocity of the optimized and non-
optimized outputs projected onto the target geometry.

Figure 7 compares the output and input trajectories for the optimized and
baseline cases in the $x-y$ plane while Figure 9(a) displays the same data as
a function of progress.

Analyzing Figures 7 and 9(a), we see that the optimized output effectively
stays within the selected $\mu=20~{}\mathrm{\mu m}$ while the baseline output
deviates by more than $100\mu m$. The optimized output trajectory speeds up in
areas with low curvature and slows down near corners or other intricate
features. Figure 9(a) also illustrates that the optimized input trajectory
deviates aggressively from the target geometry near areas where the baseline
output performs poorly. This can be interpreted as an attempt to “cancel out”
the error; The detail in Figure 7 illustrates this behavior in the $x-y$
plane.

Figure 9(b) display the same information for the airfoil test geometry with a
maximum acceleration of $2~{}\mathrm{m\,s^{-2}}$. The baseline and optimized
trajectories complete the part in $1.081~{}\mathrm{s}$. Similarly to the
letter “r” example, the optimized trajectory slows down near intricate
features, in this case the leading and trailing edges of the airfoil, and
accelerates between them. Figures 8 and 9(b) show that the optimized input
again seeks to “compensate” for deviations in the baseline trajectory.

To summarize, we observe that the optimizer exploits two main mechanisms to
improve performance of the machine relative to the baseline system with the
same trajectory time:  1. It reduces the speed of the input trajectory near
intricate geometric features and increases it in areas with low curvature.  2.
It “compensates” for errors in the baseline trajectory by moving the optimized
input trajectory in an opposing direction.  In effect, the optimizer exploits
information encapsulated in the data-driven models to push the limits of
system performance and adapt the input trajectory to the machine capabilities,
for the selected maximum acceleration value.

###  V-B Precision-accuracy Trade-off

In precision motion systems there is a fundamental trade-off between speed and
accuracy caused primarily by the inertia of the machine. In this section, we
demonstrate that our proposed methodology is able to improve system
performance by shifting this trade-off.

In our formulation (6)-(7), the trade-off between speed and accuracy is
controlled by the parameter $\mathtt{a_{max}}$ which limits the maximum
acceleration of the input trajectory. We generated trade-off curves for both
the letters and airfoil test geometries by conducting experiments for
different values of this parameter between $0.1~{}\mathrm{m\,s^{-2}}$ and
$3.3~{}\mathrm{m\,s^{-2}}$ for the letters and between
$0.1~{}\mathrm{m\,s^{-2}}$ and $10~{}\mathrm{m\,s^{-2}}$ for the airfoil, and
tested each resulting input trajectory experimentally. Similarly to the two
individual test cases presented above, the non-optimized input trajectories
are subsequently run with the same total time as the one found for the
optimized input trajectories.

We use the normalized
$L_{1}=\frac{1}{M}\sum_{i=1}^{M}\min_{s}\left|\left|\gamma_{i}-\Xi(s)\right|\right|_{2}$
and
$L_{2}=\sqrt{\frac{1}{M}\sum_{i=1}^{M}\min_{s}\left|\left|\gamma_{i}-\Xi(s)\right|\right|_{2}^{2}}$
(rms) norms as well as the
$L_{\infty}=\max_{i\in\\{1,\dots,M\\}}{\min_{s}\left|\left|\gamma_{i}-\Xi(s)\right|\right|_{2}}$
norm to quantify the deviation between the output $\gamma$ and target geometry
$\Xi$, both in simulation and experimentally.

The results are shown in Figure 10. In all cases the optimized experimental
output results in significantly more precise trajectories for a given part
completion time than the non-optimized trajectories. The input trajectories
optimized in the first stage exhibit higher deviations than the non-optimized
trajectories. This is to be expected given the prediction error of the linear
model, as quantified in Table II. Moreover, the simulated deviation is lower
than the experimental one due to model mismatch. Figure 10 also show that for
the same deviation values the optimized trajectories take less time to
complete. For example, in the letters test case the non-optimized version
takes $8.1~{}\mathrm{s}$ and achieves a deviation of $26~{}\mathrm{\mu m}$,
while the optimized version takes only $2.2~{}\mathrm{s}$ for an identical
deviation. In this case the optimized trajectory reduces the time needed to
trace the shape in 73%. For the airfoil test case a similar analysis yields a
reduction of 57% for a deviation of $47~{}\mathrm{\mu m}$. In Table III we
show the results for $5$ different test cases, each at two different
$\mathtt{a_{max}}$ scenarios. For all cases our method improves system
performance in $L_{1}$, $L_{2}$ and $L_{\infty}$ norms.

Figure 10: Precision-speed trade-off curve for the letters test case. The top panel shows the normalized $L_{2}$ deviation of the outputs from the target geometry. The bottom panel shows the total time to perform the trajectory and the $\mathtt{a_{max}}$ used in optimization. To determine the performance at a given acceleration value, one can use the lower panel to determine the time needed to complete the trajectory for a given acceleration, then retrieve the deviation for the corresponding time from the top plot. | $\mathtt{a_{max}}=1.0~{}\mathrm{m\,s^{-2}}$ | $\mathtt{a_{max}}=3.0~{}\mathrm{m\,s^{-2}}$
---|---|---
Shape | time | $L_{1}$ | $L_{2}$ | $L_{\infty}$ | time | $L_{1}$ | $L_{2}$ | $L_{\infty}$
– | $\mathrm{s}$ | % | % | % | $\mathrm{s}$ | % | % | %
Letter u | 1.052 | 52.7 | 60.0 | 63.0 | 0.618 | 36.4 | 38.7 | 46.1
Letter r | 0.807 | 58.9 | 64.3 | 75.1 | 0.473 | 32.8 | 29.7 | 42.9
Letter c | 0.807 | 80.0 | 78.5 | 74.4 | 0.473 | 40.2 | 38.6 | 54.3
Letter h | 0.971 | 26.6 | 32.8 | 31.0 | 0.568 | 22.6 | 27.0 | 43.1
airfoil | 1.322 | 73.5 | 69.0 | 64.9 | 1.024 | 75.0 | 63.0 | 28.1
TABLE III:  Percent improvement in accuracy after applying the proposed input
trajectory design method compared to the default control without input
trajectory optimization. Experimental results for different shapes and
$\mathtt{a_{max}}$.

###  V-C Effect of horizon length

Figure 11: Trajectory time and computation time of the first stage as function
of the horizon length using a receding horizon strategy. Both times are
normalized by the one-shot optimization case, the right-most points where the
horizon is equal to the number of points used in the discretization $N=128$.
In this study the tolerance band is set to $200\,\mathrm{\mu m}$.

In Figure 11 we study the effect of optimizing the reference trajectory with a
receding horizon strategy on the first stage optimization. We observe that the
trajectory time is the smallest for the one-shot optimization, and increases
with a reduction of the horizon length. The increase is more pronounced for
the Airfoil test case, which can be traced at higher speeds compared to the
Letters test case. The computation time is also the smallest for the one-shot
optimization. Short horizons are fast to compute, but require the largest
number of optimization problems to be solved. The maximum computation time is
reached for long horizons, where each optimization problem takes some time to
compute and many are required.

##  VI Conclusion

In this paper we proposed a method that improves the precision vs.
productivity trade-off of a PMS. We use models built exclusively from
experimental data, and only modify the input trajectory provided to the closed
loop control system. Experimental data obtained for shapes outside of the
training dataset corroborates simulation results and shows that the method can
significantly improve system performance, reliably shifting the precision vs.
productivity trade-off curve across a wide range of operating conditions. This
is accomplished by exploiting variable speeds (possible since the full
trajectory is optimized in one go), and error compensation while respecting
system operational constraints.

Future work will focus on reducing the computational load of the offline input
trajectory optimization and further increasing the ANN complexity to improve
the prediction accuracy. In another direction, the ANN can be combined with an
ILC loop. In an ILC setting the number of trials can be reduced by using the
ANN gradient information, either independently or with an initial solution
provided by the method we proposed here, taking advantage of the strengths of
both methods.

## Acknowledgment

The authors would like to thank Dr. Natanael Lanz for his support with the
operation of the experimental apparatus.

## References

  * [1] W.-j. Kim and D. L. Trumper, “High-precision magnetic levitation stage for photolithography,” Precision engineering, vol. 22, no. 2, pp. 66–77, 1998\.
  * [2] J. J. Gorman and N. G. Dagalakis, “Force control of linear motor stages for microassembly,” in ASME International Mechanical Engineering Congress and Exposition, vol. 37130, pp. 615–623, 2003.
  * [3] T. Kim and J. J. Gorman, “A 3-dof mems motion stage for scanning tunneling microscopy,” in 2022 IEEE 35th International Conference on Micro Electro Mechanical Systems Conference (MEMS), pp. 470–472, IEEE, 2022.
  * [4] M. Khosravi, C. Koenig, M. Maier, R. S. Smith, J. Lygeros, and A. Rupenyan, “Safety-aware cascade controller tuning using constrained bayesian optimization,” IEEE Transactions on Industrial Electronics, pp. 1–1, 2022\.
  * [5] K. K. Tan, T. H. Lee, and S. Huang, Precision Motion Control: Design and Implementation.  Springer Science & Business Media, 2007.
  * [6] P. Ouyang, R. Tjiptoprodjo, W. Zhang, and G. Yang, “Micro-motion devices technology: The state of arts review,” The International Journal of Advanced Manufacturing Technology, vol. 38, no. 5, pp. 463–478, 2008.
  * [7] M. Ghobakhloo, “The future of manufacturing industry: a strategic roadmap toward industry 4.0,” Journal of Manufacturing Technology Management, 2018\.
  * [8] N. Kouraytem, X. Li, W. Tan, B. Kappes, and A. D. Spear, “Modeling process–structure–property relationships in metal additive manufacturing: a review on physics-driven versus data-driven approaches,” Journal of Physics: Materials, vol. 4, no. 3, p. 032002, 2021.
  * [9] J. S. Srai, M. Kumar, G. Graham, W. Phillips, J. Tooze, S. Ford, P. Beecher, B. Raj, M. Gregory, M. K. Tiwari, et al., “Distributed manufacturing: scope, challenges and opportunities,” International Journal of Production Research, vol. 54, no. 23, pp. 6917–6935, 2016.
  * [10] D. Lam, C. Manzie, and M. Good, “Model predictive contouring control,” in 49th IEEE Conference on Decision and Control (CDC), pp. 6137–6142, IEEE, 2010.
  * [11] J. L. Vázquez, M. Brühlmeier, A. Liniger, A. Rupenyan, and J. Lygeros, “Optimization-based hierarchical motion planning for autonomous racing,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 2397–2403, 2020.
  * [12] Y. Wang, F. Gao, and F. J. Doyle III, “Survey on Iterative Learning Control, Repetitive Control, and Run-to-Run control,” Journal of Process Control, vol. 19, no. 10, pp. 1589–1600, 2009.
  * [13] E. C. Balta, K. Barton, D. M. Tilbury, A. Rupenyan, and J. Lygeros, “Learning-based repetitive precision motion control with mismatch compensation,” arXiv preprint arXiv:2111.10246, 2021.
  * [14] K. L. Barton and A. G. Alleyne, “A norm optimal approach to time-varying ILC with application to a multi-axis robotic testbed,” IEEE Transactions on Control Systems Technology, vol. 19, no. 1, pp. 166–180, 2010.
  * [15] D. Liao-McPherson, E. C. Balta, A. Rupenyan, and J. Lygeros, “On robustness in optimization-based constrained Iterative Learning Control,” arXiv preprint arXiv:2203.05291, 2022.
  * [16] K. Zhang, A. Yuen, and Y. Altintas, “Pre-compensation of contour errors in five-axis CNC machine tools,” International Journal of Machine Tools and Manufacture, vol. 74, pp. 1–11, 2013.
  * [17] C.-C. Lo and C.-Y. Hsiao, “CNC machine tool interpolator with path compensation for repeated contour machining,” Computer-Aided Design, vol. 30, no. 1, pp. 55–62, 1998.
  * [18] T. Haas, S. Weikert, and K. Wegener, “MPCC-based set point optimisation for machine tools,” International Journal of Automation Technology, vol. 13, no. 3, pp. 407–418, 2019.
  * [19] X. Yang, R. Seethaler, C. Zhan, D. Lu, and W. Zhao, “A model predictive contouring error precompensation method,” IEEE Transactions on Industrial Electronics, vol. 67, no. 5, pp. 4036–4045, 2019.
  * [20] S. Di Cairano, A. Goldsmith, U. V. Kalabić, and S. A. Bortoff, “Cascaded reference governor–MPC for motion control of two-stage manufacturing machines,” IEEE Transactions on Control Systems Technology, vol. 27, no. 5, pp. 2030–2044, 2018.
  * [21] Z. Wang, C. Hu, Y. Zhu, and L. Zhu, “Prediction-model-based contouring error iterative precompensation scheme for precision multiaxis motion systems,” IEEE/ASME Transactions on Mechatronics, vol. 26, no. 5, pp. 2274–2284, 2021\.
  * [22] Y. Jiang, J. Chen, H. Zhou, J. Yang, P. Hu, and J. Wang, “Contour error modeling and compensation of cnc machining based on deep learning and reinforcement learning,” The International Journal of Advanced Manufacturing Technology, vol. 118, no. 1, pp. 551–570, 2022.
  * [23] H. Kim and C. E. Okwudire, “Simultaneous servo error pre-compensation and feedrate optimization with tolerance constraints using linear programming,” The International Journal of Advanced Manufacturing Technology, vol. 109, no. 3, pp. 809–821, 2020.
  * [24] H. Kim and C. E. Okwudire, “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, vol. 115, no. 7, pp. 2069–2082, 2021.
  * [25] S. Balula, A. Liniger, A. Rupenyan, and J. Lygeros, “Reference design for closed loop system optimization,” in 2020 European Control Conference (ECC), pp. 650–655, IEEE, 2020.
  * [26] R. Eppler, “Airfoil data,” in Airfoil Design and Data, pp. 163–512, Springer, 1990.
  * [27] G. Masetti and F. D. Giandomenico, “Analyzing forward robustness of feedforward deep neural networks with LeakyReLU activation function through symbolic propagation,” in Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pp. 460–474, Springer, 2020.
  * [28] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [29] A. Paszke et al., “Pytorch: An imperative style, high-performance deep learning library,” in Advances in Neural Information Processing Systems 32, pp. 8024–8035, Curran Associates, Inc., 2019.
  * [30] A. Aboanber and Y. Hamada, “Generalized Runge–Kutta method for two-and three-dimensional space–time diffusion equations with a variable time step,” Annals of Nuclear Energy, vol. 35, no. 6, pp. 1024–1040, 2008.
  * [31] G. Goodwin, M. M. Seron, and J. A. De Doná, Constrained Control and Estimation: An Optimisation Approach.  Springer Science & Business Media, 2006.
  * [32] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematical programming, vol. 106, no. 1, pp. 25–57, 2006.
  * [33] I. Dunning, J. Huchette, and M. Lubin, “JuMP: A modeling language for mathematical optimization,” SIAM Review, vol. 59, no. 2, pp. 295–320, 2017\.
  * [34] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi – A software framework for nonlinear optimization and optimal control,” Mathematical Programming Computation, vol. 11, no. 1, pp. 1–36, 2019.

|  Samuel Balula received his B.Sc. and M.Sc. in Engineering Physics from
Instituto Superior Técnico, University of Lisbon, Portugal in 2014 and 2016
respectively. Between 2017-2018 he was head of research at a tech startup. He
is currently a PhD student in the Automatic Control Laboratory at ETH Zürich,
Switzerland. His research interest include optimal non-linear control,
trajectory planning, data-driven control and optimal experimental design with
applications in robotics and manufacturing.
---|---
|  Dominic Liao-McPherson received his BASc. in Engineering Science from the
University of Toronto in 2015 and his PhD. in Aerospace Engineering and
Scientific Computing from the University of Michigan, Ann Arbor, in 2020. He
is currently a postdoctoral researcher at the ETH Zürich Automatic Control
Laboratory. His research interests include real-time optimization, game-
theoretic multi-agent control and numerical algorithms with applications in
robotics, energy, and manufacturing.
---|---
|  Alisa Rupenyan received the BSc in Engineering physics and the MSc degree
in Laser physics from the University of Sofia, Sofia, Bulgaria, in 2004 and
2005, and the Ph.D. degree from the Department of Physics and Astronomy of
Vrije Universiteit Amsterdam, Amsterdam, The Netherlands in 2010. Between
2011-2014 she was a postdoctoral fellow at ETH Zürich, Switzerland, and
between 2014-2018 a lead scientist in a robotic start-up. She is currently a
group leader in automation at inspire, the technology transfer unit at ETH
Zürich, and a senior scientist and PI at the Automatic Control Laboratory
there. Her research interests are on the intersection between machine
learning, control, and optimization for industrial applications and robotics.
---|---
|  John Lygeros received a B.Eng. degree in 1990 and an M.Sc. degree in 1991
from Imperial College, London, U.K. and a Ph.D. degree in 1996 at the
University of California, Berkeley. After research appointments at M.I.T.,
U.C. Berkeley and SRI International, he joined the University of Cambridge in
2000 as a University Lecturer. Between March 2003 and July 2006 he was an
Assistant Professor at the Department of Electrical and Computer Engineering,
University of Patras, Greece. In July 2006 he joined the Automatic Control
Laboratory at ETH Zurich where he is currently serving as the Professor for
Computation and Control and the Head of the laboratory. His research interests
include modelling, analysis, and control of large scale systems, with
applications to biochemical networks, energy systems, transportation, and
industrial processes. John Lygeros is a Fellow of IEEE, and a member of IET
and the Technical Chamber of Greece. Since 2013 he is serving as the Vice-
President Finances and a Council Member of the International Federation of
Automatic Control and since 2020 as the Director of the National Center of
Competence in Research “Dependable Ubiquitous Automation” (NCCR Automation).
---|---


11institutetext<EMAIL_ADDRESS>
# Learning From How Humans Correct

Tong Guo

###### Abstract

In industry NLP application, our manually labeled data has a certain number of
noisy data. We present a simple method to find the noisy data and re-label
them manually, meanwhile we collect the correction information. Then we
present novel method to incorporate the human correction information into deep
learning model. Human know how to correct noisy data. So the correction
information can be inject into deep learning model. We do the experiment on
our own text classification dataset, which is manually labeled, because we re-
label the noisy data in our dataset for our industry application. The
experiment result shows that our method improve the classification accuracy
from 91.7% to 92.5%. The 91.7% accuracy is trained on the corrected dataset,
which improve the baseline from 83.3% to 91.7%.

###### Keywords:

Deep Learning Text Classification NLP

##  1 Introduction

In recent years, deep learning [2] and BERT-based [1] model have shown
significant improvement on almost all the NLP tasks. However, past methods did
not inject human correction information into the deep learning model. Human
interact with the environment and learn from the feedback from environment to
correct the their own error or mistake. Our method try to solve the problem
that let deep learning model imitate how human correct.

In order to solve the problem we present a learning framework. The framework
mainly works for our industry application dataset, because the problem starts
from our industry application. In order to solve the text classification
problem in our industry application, we first label a dataset. Then we can
find the noisy data in the dataset and re-label them. The re-labeling step
collects the human correction information.

To the best of our knowledge, this is the first study exploring the
improvement of injecting human correction information into deep model for
natural language understanding. Our key contribution are 4 folds:

1\. Based on our dataset, we first present the simple step to find the noisy
data and re-label the noisy data.

2\. We present the method to inject the human correction information into BERT
for text classification.

3\. The experiment result shows our framework gets gain of 0.8% accuracy
against the strong baseline.

4\. Our learning framework can apply to a broad set of deep learning industry
applications whose dataset is manually labeled.

Figure 1: Our learning framework

##  2 Relate work

BERT [1] is constructed by the multi-layer transformer encoder [10], which
produces contextual token representations that have been pre-trained from
unlabeled text and fine-tuned for the supervised downstream tasks. BERT
achieved state-of-the-art results on many sentence-level tasks from the GLUE
benchmark [3]. There are two steps in BERT’s framework: pre-training and fine-
tuning. During pre-training, the model is trained on unlabeled data by using
masked language model task and next sentence prediction task. Apart from
output layers, the same architectures are used in both pre-training and fine-
tuning. The same pre-trained model parameters are used to initialize models
for different down-stream tasks.

Our method is different to semi-supervised learning. Semi-supervised learning
solve the problem that making best use of a large amount of unlabeled data.
These works include UDA [6], Mixmatch [7], Fixmatch [8], Remixmatch [9]. These
works do not have the human correction information. Our work has a clear goal
that is to learn how humans correct their mistake. [11] studies the PN
learning problem of binary classification. But our idea can apply to a broad
set of deep learning industry applications.

The work[12] of OpenAI also use model predictions for human as references to
label. The work[12] use the new human-labeled data from model predictions to
train the reward model. Our work do not have the reward model of [12], and I
think reinforcement learning from human feedback(RLHF) is same to re-label
here. Also the new human-labeled data from model predictions do not conflict
to our method, because the new data can be simply merged to the whole dataset
(of [12]’s policy model), and then to re-label/re-correct. Correcting the
[12]’s policy by reward model is same to correct all the related data/labels
in training dataset. The detail pattern-based/substring-based human-correct
method is illustrated at [13].

##  3 Our method

In this section, we describe our method in detail. Our learning framework is
shown in Fig 1. The framework includes 5 steps:

Step 1, in order to solve the industry text classification problem. We label
2,040,000 data and split them into 2,000,000 training data and 40,000 test
data. The 2,040,000 data are sampled from our application database whose data
size is 500 million.

Step 2, we train / fine-tune the BERT model on the 2,000,000 training data. We
named the result model of this step Model-A. Note that Model-A should not
overfit the training dataset.

Step 3, we use Model-A to predict for all the 2,040,000 data. Then we find
300,000 data whose predicted label and human label are not equal. We consider
it is the noisy data. In detail, there are 294,120 noisy data in the training
dataset and 5,880 noisy data in the test dataset.

Step 4, we manually re-label 300,000 data and merge back to the 2,040,000
data. Then we get the merged 2,040,000 data. During the re-labeling, the last
label by human and the model-A’s predicted label are listed as the references
for people. But labeling people also make their own decision.

Step 5, we add the Model-A’s predicted one-hot label as the addition input for
training / fine-tuning a new model. The detail encoding method for the
predicted label is described in the next section. We named the result model of
this step Model-C. The Model-A’s predicted one-hot labels represent the
before-corrected information. The training ground truth for Model-C contains
the 294,120 corrected human label, which represent the corrected information.
So Model-C is learning how to correct and learning the text classification
task in the same time.

Figure 2: The encoding detail for Model-A, which is corresponding to Fig 1.
Figure 3: The encoding detail for Model-C, which is corresponding to Fig 1.

##  4 The Model

We use BERT as our model. The training steps in our method belongs to the
fine-tuning step in BERT. We follow the BERT convention to encode the input
text. The encoding detail for Model-A is shown in Fig 2. The encoding detail
for Model-C is shown in Fig 3.

##  5 Experiments

In this section we describe detail of experiment parameters and show the
experiment result. The detail result is shown in Table 2. The data size in our
experiment is shown in Table 1.

In fine-tuning, we use Adam [4] with learning rate of 1e-5 and use a dropout
[5] probability of 0.1 on all layers. We use BERT-Base (12 layer, 768 hidden
size) as our pre-trained model.

Table 1: The data size for text classification Data Size | Description
---|---
500,000,000 | All the data in our application database.
2,040,000 | The data we label in step 1 of Fig 1.
2,000,000 | The the training dataset we split from the 2,040,000 data.
300,000 | All the noisy data we select from the 2,040,000 data to re-label.
294,120 | The noisy data in the training dataset.
40,000 | The test dataset we split from the 2,040,000 data.
5,880 | The noisy data in the testing dataset.
Table 2: The experiment result of text classification. The test dataset is the 40,000 data. The accuracy reported by human is 5000 data that sampled from the 500 million data. Model-A, Model-B and Model-C are corresponding to Fig 1. Model In Fig 1 | Test Dataset Accuracy | Human Evaluate Accuracy
---|---|---
Model-A | 83.3% | 88.0%
Model-B | 91.7% | 97.2%
Model-C | 92.5% | 97.7%

##  6 Analysis

In step 4 of Fig 1, the manually re-labeling can correct the noisy 300,000
data. Because the selected 300,000 data is unfitting to the Model-A, the re-
labeling’s ’worst’ result is the origin last human’s label.

The step 5 in Fig 1 is the core contribution of our work. In step 5, the
predicted label of Model-A contains the before-human-corrected information.
The ground truth for Model-C contains the after-human-corrected information.
So the model is learning the human correction.

We could use the before-corrected human label (i.e., the ground truth in step
1 of Fig 1) as the input for Model-C. But this way can not apply to real
industry inference. Because we can not get the before-corrected human label as
the input of Model-C in real industry application. In real industry
application, we use the Model-A to predict one-hot label as input for Model-C.

Human evaluation accuracy is higher than the test dataset accuracy, because we
randomly sampled 5000 data from the 500 million data. The sampled 5000 data
represents the great majority of the 500 million data.

Why re-label method work? Because deep learning is statistic-based. Take
classification as example. (In a broad sense, all the machine learning tasks
can be viewed as classification.) If there are three very similar data
(data-1/data-2/data-3) in total, which labels are class-A/class-A/class-B.
Then the trained model will probably predict class-A for data-3. We assume
that data-3 is wrong labeled to class-B by human , because more people label
its similar data to class-A.

##  7 Conclusion

Human interact with environment and learn from the feedback from environment
to correct human’s own error. Base on the human’s correction idea, we design a
learning framework to inject the information of human’s correction into the
deep model. The experiment result shows our idea works. For further
improvement, we will try to re-label the data which top-1 predicted score and
top-2 predicted score are very close. It means that the model can hardly
classify these data and the dataset contains noisy data here. Our idea can
apply to a broad set of deep learning industry applications.

## References

  * [1] Devlin J, Chang M W, Lee K, et al. Bert: Pre-training of deep bidirectional transformers for language understanding[J]. arXiv preprint arXiv:1810.04805, 2018.
  * [2] Krizhevsky A, Sutskever I, Hinton G E. Imagenet classification with deep convolutional neural networks[J]. Advances in neural information processing systems, 2012, 25: 1097-1105.
  * [3] Wang A, Singh A, Michael J, et al. GLUE: A multi-task benchmark and analysis platform for natural language understanding[J]. arXiv preprint arXiv:1804.07461, 2018.
  * [4] Kingma D P, Ba J. Adam: A method for stochastic optimization[J]. arXiv preprint arXiv:1412.6980, 2014.
  * [5] Srivastava N, Hinton G, Krizhevsky A, et al. Dropout: a simple way to prevent neural networks from overfitting[J]. The journal of machine learning research, 2014, 15(1): 1929-1958.
  * [6] Xie Q, Dai Z, Hovy E, et al. Unsupervised data augmentation for consistency training[J]. arXiv preprint arXiv:1904.12848, 2019.
  * [7] Berthelot D, Carlini N, Goodfellow I, et al. Mixmatch: A holistic approach to semi-supervised learning[J]. arXiv preprint arXiv:1905.02249, 2019.
  * [8] Sohn K, Berthelot D, Li C L, et al. Fixmatch: Simplifying semi-supervised learning with consistency and confidence[J]. arXiv preprint arXiv:2001.07685, 2020.
  * [9] Berthelot D, Carlini N, Cubuk E D, et al. Remixmatch: Semi-supervised learning with distribution alignment and augmentation anchoring[J]. arXiv preprint arXiv:1911.09785, 2019.
  * [10] Vaswani A, Shazeer N, Parmar N, et al. Attention is all you need[C]//Advances in neural information processing systems. 2017: 5998-6008.
  * [11] Northcutt C G, Wu T, Chuang I L. Learning with confident examples: Rank pruning for robust classification with noisy labels[J]. arXiv preprint arXiv:1705.01936, 2017.
  * [12] Ouyang L, Wu J, Jiang X, et al. Training language models to follow instructions with human feedback[J]. arXiv preprint arXiv:2203.02155, 2022.
  * [13] Guo, Tong (2022): Re-Label By Data Pattern Is All You Need For Knowledge Driven Deep Learning. TechRxiv. Preprint. https://doi.org/10.36227/techrxiv.20485917.v7