1 Introduction
A precision motion stage (PMS) is a robotic setup designed to precisely position an endeffector, typically in the micro to nanometer range. PMS are used in lithography [kim1998high], microassembly [gorman2003force], precision metrology [kim20223], and precision engineering [yang2018two, Khosravi2022], and many other applications [tan2007precision, ouyang2008micro]. 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 inertia of PMS systems dictates a tradeoff between productivity and accuracy and, in practice, the control strategy is often the determining factor in this tradeoff. Typical industrial precision systems achieve highprecision using welltuned “classical” feedback controllers, such as ProportionalIntegralDerivative (PID), to track a given target trajectory. These lowlevel controllers are designed by the manufacturer and users are typically not able to modify/tune the control algorithm.
In this paper, we propose an optimizationbased reference (input trajectory) design 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 design principles:

does not require modifications to the existing hardware and lowlevel controllers, making it suitable for application in production machines;

does not require a physicsbased model of the machine, relying exclusively on experimental data;

can be applied to individual machines with diverse technical characteristics;

can produce good results for previously unseen parts at the first attempt.
These design principles 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 [ghobakhloo2018future]. This data can be exploited to improve performance; for example, it can be leveraged to create datadriven models, avoiding the labour intensive process of developing and maintaining first principles based ones [kouraytem2021modeling]. Future manufacturing is also expected to be more distributed and personalized (e.g., “Manufacturing as a service”) [srai2016distributed, tao2011cloud], leading to smaller volumes for any given part, possibly produced across a pool of nonhomogeneous machines. This makes it more important to exploit data to ensure that machines can execute designs with minimal calibration.
1.1 Stateoftheart
1.1.1 Real time control
One method to improve the performance of PMS is by redesigning the feedback controller. A common approach is model predictive contouring control (MPCC) [lam2010model, YANG201550], wherein a model of contouring error is minimized in a receding horizon fashion, drawing on methods from autonomous racing [Vazques2020]. However these methods have high online computational loads, rely on accurate process models, and require modifying the innerloop controllers.
1.1.2 Learningbased control
Another common approach in manufacturing and precision motion systems is learningbased control [wang2009survey, ahn2007iterative, balta2021learning, altin2014robust]. 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 RuntoRun control), or periods (Repetitive Control). These methods can exploit models [barton2010norm, lu2017nonlinear], include constraints [mishra2010optimization, liao2022robustness] 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.
1.1.3 Precompensation
Precompensation refers to a class of methods that aim to modify the input trajectory offline to reduce the contouring error for arbitrary parts [zhang2013pre]. 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) [lo1998cnc]. These methods cannot include lookahead or constraints, motivating the use of offline MPCC algorithms based on physicsbased nonlinear[haas2019mpcc], identified linear models [yang2019model] or combined with a reference governor [di2018cascaded]. However, these generate suboptimal input trajectories since the optimization uses a receding horizon, instead of optimizing the trajectory as a whole.
In [jiang2022contour]
the contour error is predicted with an Artificial Neural Network (ANN), and the precompensation computed with reinforcement learning. This approach does not allow constraints to be included and uses only zerothorder information about the ANN model, despite the ready availability of ANN derivatives.
1.2 Contributions
In this paper, we propose a novel precompensation method based on datadriven modelling and trajectory optimization to improve simultaneously the productivity and performance of PMS. Our contributions are threefold:

A datadriven framework for PMS modelling and identification for control, including design of experiments for collecting training data, and machine learning architecture in a feedforward form.

A novel contour control precompensation method featuring timeoptimal trajectory design and motion constraints.

Experimentally validated improvement over the existing controller across the precision/productivity tradeoff curve.
Our approach does not require any modification of the design or control software of the system (D1) and is fully datadriven (D2). Moreover, it enables adaptation to changing conditions or system degradation through retraining the underlying datadriven model (D3). Finally, the approach is independent of the target geometry and significantly shifts the precision/productivity tradeoff curve^{1}^{1}1Relative to a method without precompensation, our proposed method is 3060% more accurate for new geometries for same trajectory time budget or alternatively, trajectories can be traced in 2550% of the time for the same contouring accuracy. (D4). We build upon our previous work [balula2020reference] based on linear system identification. Here we improve and extend this approach by incorporating a nonlinear ANNbased model of system performance, which offers two orders of magnitude higher precision, a refined optimization approach, and experimental validation of the simulation results.
1.3 Organization
In section 2 we precisely define the problem tackled, while in section 3 we describe the modelling structure, the model fitting strategy, and the accuracy achieved. In section 4 we detail the input trajectory design strategy proposed, specifying the optimization problems. In section 5 we show experimental validation for the method with a set of test cases.
2 Problem Formulation
2.1 Target geometry tracking problem
Consider a tooltip, e.g., the laser in a laser cutting machine, which traverses a twodimensional (2D) workspace . The position of the tooltip as a function of time is denoted by
(1) 
The ultimate goal in precision motion planning is to have the tooltip trace a spatial target geometry
(2) 
where is the path length. The target geometry is assumed to be a continuous function parameterized by the path progress variable . An example is shown in Figure 1. 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
Most industrial precision motion systems come with a builtin lowlevel controller designed by the manufacturer to track an input trajectory , where is the time necessary for the tooltip to complete tracing the target geometry. When is given as an input to the lowlevel controller, its actions and the dynamics of the machine give rise to an output trajectory ; note that for simplicity we assume that both the input and the output trajectory lie in the workspace . We denote the mapping between the commanded input trajectory and the resulting output trajectory by
(3) 
where the function encapsulates the closedloop dynamics of the machine, and
is a noise term, with zeromean and a typical standard deviation of
, which quantifies the repeatibility of the system.Ideally, the lowlevel controller is well designed so that and is approximately the identity function.
We are interested in the following inverse problem: Given the target geometry , how do we generate an input trajectory for the machine, so that the resulting output matches as closely as possible. In addition to the usual difficulties associated with inverse problems, for precision motion systems the function is often poorly known, partly due to lack of information about the lowlevel controller.
In standard practice, input trajectories are generated by
traversing the path at a constant speed , i.e. setting . However, due to sensor noise, actuator limits, the inertia of the machine, and other factors we observe experimentally a significant error between and , especially when is high.
This induces a natural tradeoff between speed and precision that is characteristic of precision motion systems. Our aim is to improve this tradeoff by learning the function through preliminary experiments, then inverting it using optimization.
Formally, let be the deviation function, defined as the projection of the output into the target geometry
(4) 
Our goal is then to solve the optimization problem [c1] ρT(γ)+ λ∫_0^S d(s,γ) ds γ= F(ρ) d(s,γ) ≤ μ(s), ∀s ∈[0,S] , where is the duration of the output trajectory, governs the tradeoff between speed and accuracy, and is a tolerance bound around the target geometry; for the typical industrial laser cutting systems mentioned above, one would use . The different quantities are illustrated in Figure 2.
The optimization problem (4) 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 (4) using only data gathered from the process. Because the methodology is datadriven 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.
2.2 Experimental Setup
features two degrees of freedom,
and , driven by linear actuators. The end effector is mounted on the axis, which in turn is mounted on top of the 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, labelled 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 DXLLM325 two stage motion machine. The machine is instrumented with quadrature encoders that measure the position of both axes and is actuated with ETEL LMS15100 linear motors.
The system is equipped with a feedback controller that adjusts the motor voltages to track a supplied input trajectory . 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 blackbox and consider only the desired input trajectory 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 1.
Property  Value  Unit 

Maximum acceleration  
Maximum velocity  
Working area  
Lowlevel control loop frequency  
Repeatability (standard deviation) 
2.3 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 highlift Eppler e344 airfoil[eppler1990airfoil]. Geometries of this kind are typically laser cut out of e.g., balsa wood, and are used in the primary structural subassembly of an aircraft. The main challenging features of the airfoil are a cutout 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.
2.4 Architecture overview
We approach the trajectory design problem in several stages. First, we design experiments to generate informative inputoutput data for the unknown mapping and use the resulting data to identify two models for the system, a lowfidelity linear one and a highfidelity one based on an ANN. Throughout, we incorporate machine learning best practices, such as independent training and validation datasets, 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 dataset, enabling optimization of new geometries.
Using these models, we propose a twostage optimizationbased approach for offline optimization of references as schematized in Figure 4. 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 firststage is used as the initial condition for the second stage, which employs the highfidelity neural network model to correct errors from the first stage. The resulting input trajectory is then given to the lowlevel controller as an input trajectory.
3 Data driven modelling
As a first approximation, we model the system using a continuoustime linear model which captures the dominant dynamics. The structure of this model allows for efficient computations of a timeoptimal solution in the first stage, but lacks the highprecision required by the application. This shortcoming is alleviated by a second highfidelity neural network model, used to refine the initial solution provided by the lowerfidelity linear model.
3.1 Lowfidelity linear model
The closedloop system maps trajectories to trajectories and is of the form . We approximate this mapping using a linear time invariant (LTI) state space model of the following form
(5a)  
(5b) 
The model (5) uses a nonphysical hidden state ; the matrices , , , and 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. In particular dimensions higher than have practically the same prediction error, but with added complexity.
The linear model captures the dominant system dynamics (model validation results are presented in Section 3.4), allowing us to formulate and solve a tractable optimization problem that effectively trades off the total trajectory time and its accuracy, as described in Section 4.1. However, it fails to capture some of the nonlinearities inherent in the system and is not sufficiently accurate to achieve the precision required by the application. On the other hand, the computational burden of using only the more complex nonlinear model (described in Section 3.2) without the linear model would limit the method to unrealistic small problem instances.
3.2 Highfidelity Nonlinear Model
To obtain a higher precision model, we use an ANN to capture the nonlinearities of the system dynamics. ANN are much more expressive than linear models, but require more informative training data to avoid over fitting; the latter is not a drawback in our application, where a large amount of data can be collected at low cost.
We also investigated other alternatives, such as Gaussian Processes, but found that the ANN architecture below offered a more favourable tradeoff between accuracy and computation effort. They are well suited for our application, where a large amount of data can be collected at low cost.
The 2Dpositioning stage is modelled using two independent causal ANN, one for each axis. Unlike (3.1) the nonlinear model operates in discrete time, at a sample frequency of and it has the following form
(6) 
Where is the output, is the nonlinear function defined by the ANN, is the input trajectory given to the closed loop system, and is the length of the input history considered for the prediction. Note that we have introduced the variables in place of , and in place of to make explicit that and are parameterized in continuous time, while and are defined in discrete time.
The two ANN are structurally identical, but are trained separately, leading to different weights. Each takes as input the most recent of the input trajectory for the corresponding axis, subsampled at , leading to inputs. Each network has a single output, namely the predicted position at the subsequent time step. The networks have fully connected hidden layers of LeakyReLU [masetti2020analyzing] activation functions with a slope of , and a triangular structure as illustrated in Figure 5. The length of the time history was selected to be approximately times the timescale of the slowest mode of the identified linear model. The sample rate was selected by analizing 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 feedforward 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 3.4. Finally, note that the nonlinear and linear models are independent, in particular the predictions of the linear model are not used in the nonlinear model evaluation.
3.3 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. Classical system identification methods make use of a toolbox of excitation signals such as Pseudo Random Binary Signal, sinusoidal, or chirp signals. These inputs have the property of being persistently exciting for linear systems. Unfortunately none of these proved adequate for the identification of the nonlinear system in this paper.
3.3.1 Linear Model
For the identification of the linear model we opted for randomly generated trajectories, yielding 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.
3.3.2 Nonlinear Model
Having sufficient highquality 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 4 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 trainingoptimisationretraining 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 favoured 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. Our experiments suggest that there is no added benefit in refitting the linear model, since its small number of parameters is already identified well with the initial dataset.
A representative training data sample is shown in Figure 6, and the flow of data is summarized in Figure 7. In section 3.4 we present the prediction quality of a third generation ANN to generate the results of Section 5. The data was divided for training and testing purposes. This results in a ratio of data points per parameter of the ANN. About 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 , subsampled to achieve an effective sample rate of
, and scaled. For the loss function we use the mean square of the prediction error. We use batches of
segments which take maximum advantage of the GPU computation power. Adam [kingma2014adam] was used for the optimizer, with decaying learning rate. Modelling and training were done in PyTorch[NEURIPS2019_9015] with NVIDIA CUDA support.3.4 Model Validation
Training /  
Validation  
Model  Axis  
–  –  
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  
Nonlinear  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 



The quality of the model predictions is evaluated over different datasets and acceleration limits. The results are summarized in Table 2 and Figure 8.
The linear model training data set includes trajectories with a range of accelerations up to . However, due to its simple structure, the linear model cannot overfit the data, hence cannot accurately predict all points in its training data 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.
Regarding the nonlinear model, the prediction accuracy is best for data with the same characteristics as the one used for the training, as can be seen in Figure 7(a) where the error is evaluated for the validation dataset. The empirical error distribution has zero mean. When presented with data coming from shapes which have not been used for the training of the model (Figures 7(b) and 7(c)), the prediction accuracy degrades. However, the accuracy is still adequate for the application, as the standard deviations remain below for each axis for accelerations below as shown in Figure 7(c).
Figure 9 shows how the standard deviation of the prediction evolves with the maximum acceleration that was set as a constraint in the optimisation of Section 4. 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.
4 Trajectory Optimization
With the the datadriven models in hand, we now discuss the twostage 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 of the trajectory in a computationally tractable manner. This is accomplished by using a variable timediscretization in the first stage, enabling inclusion of as a decision variable. The second stage then fixes the duration and refines the accuracy of the firststage solution using the highfidelity 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 highfidelity continuous time model needed for a computationally tractable variabletime 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 2. 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 offline, 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, as shown schematically in Figure 7.
4.1 Firststage Optimization Formulation
The purpose of the first stage is to fix the duration of the trajectory, , by trading off speed and precision. We employ a contouring control approach with a fixed spatial discretization. This involves sampling the target geometry at equally spaced points in space , where
(7) 
We use the identified linear model (5) to approximate the process dynamics. As our approach involves a fixed spatial discretization , the timediscretization must be variable. Starting with , we denote by as the time at which we would like the trajectory to reach the point and consider the nonuniform time steps as decision variables; note that the time at which the trajectory is completed is implicitly also a decision variable. We then obtain a discrete time model by applying the th order RungeKutta (RK4)[aboanber2008generalized] timemarching method to (5)
(8a)  
(8b) 
Where is the RK4 function, and will be treated as decision variables.
To evaluate the contouring error we start with the distance and introduce a local coordinate frame as shown in Figure 10
. The linear transformation
(9) 
where is the rotation angle of the local frame, encodes a translation and rotation between the local coordinate frame and the global coordinate frame. With
(10) 
the output can be transformed from local coordinates to global coordinates . Since time discretization is variable, we can always ensure that the component of is zero by adding a constraint to the optimization problem. The contouring error is then simply the component of . The tolerance bound condition in (4) can then also be computed for the corresponding discretization point .
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 Δτ, ρ, γ, z, ν ∑_k=1^N1 Δτ_k + 1N2 ∑_k=2^N1 ∥∂^2 ρ_k∥_2^2 z_k+1= f_lm(z_k, ρ_k, Δτ_k), k = 1, …, N1 γ_k= C z_k + D ρ_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, …, N1 ∥ ∂^2 γ_k ∥_∞≤ a_max, k = 2, …, N1 , where is the number of discretization points, selected based on the total path length and tolerance bound, , are discrete derivative operators, is the output trajecotry in local coordinates, is the output trajectory in global coordinates, is the input trajectory, is the RK4 approximation of the linear model (8), with identified matrices and , is the transformation matrix between global and local coordinates in (9), is the discretized maximum allowed deviation from the target geometry, are the limits of the working space of the device, is the maximum allowable speed, and is the acceleration limit.
The cost (4.1) 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, (4.1) and (4.1) imposes the RK4 discretization of the linear dynamics (8), (4.1) requires the variable time step to be nonnegative, (4.1) imposes that the output stays within the working space of the device, (4.1) the transformation between local and global coordinates in (10), (4.1) requires that the component of is zero, in which case the component is the deviation, (4.1) bounds the component of to be less than the maximum allowed deviation, (4.1) and (4.1) impose componentwise limits on velocity and acceleration of the output.
4.2 Secondstage Optimization Formulation
Experiments applying the trajectory generated by (10) directly to the machine suggest that this generally results in large deviations from the target geometry, which we attribute to the low fidelity of the linear model; the data is omitted in the interest of space. To improve performance, we refine the trajectory generated by the first stage using the highfidelity ANN model (6). In principle, it is possible to skip the first stage and directly embed the ANN model in an optimization problem similar to (10). 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 (10) 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 equally spaced points in time
(11) 
, where 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 secondstage 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 r, ϕ ∑_i=1^M ∥ϕ_i  ¯Ξ_i∥^2_2 ϕ_i = f_nlm(r_ih, …, r_i1), i=1,…,M — ϕ_i  ¯Ξ_i — ≤ μ_i i=1,…,M ϕ_i ∈ W i=1,…,M — ∂ ϕ_i — ≤ v_max i=1,…,M1 — ∂^2 ϕ_i — ≤ a_max i=2,…,M1 , where is the number of discretization points, is the output, is the input trajectory, is the sampled target geometry, is the nonlinear ANN model (6), is the tolerance bound evaluated for each discretization point, are the limits of the working space of the device, is the maximum allowable speed, and is the acceleration limit.
The cost (4.2) penalizes deviations of the output from the target geometry, in this case, and contrary to the first stage, we do not include a smoothing term. Trajectories that correct the error, even with aggressive input trajectories (as long as they comply with the constraints) are acceptable. The constraint (4.2) imposes the nonlinear system dynamics modelled by the ANN, (4.2) constrains the output to be within the tolerance bound for each discretization point, While (4.2), (4.2) and (4.2) 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 penalty which is implemented using slack variables [goodwin2006constrained]. Note that the limits the maximum acceleration of the output trajectory and it is a decisive parameter determining the total trajectory time. It is the main tuning parameter of the method along the accuracy/trajectory time tradeoff curve.
4.3 Implementation Details
The optimization problems defined in (10) and (11) are both nonlinear problems, which were solved using IPOPT [Ipopt], with the JuMP [JuMP], and Casadi [Casadi] interfaces in Julia and Python respectively. The first stage optimization problem (10) is initialized by taking as the target geometry traversed at a constant speed. The solution of (10) is then used to initialize the secondstage (11).
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 of dedicated memory, an Intel(R) Core(TM) i99900K
CPU and of RAM.
For the problem sizes solved in this paper, the firststage optimization takes approximately 1 minute to complete (), and the secondstage takes approximately between 1 to 4 hours depending on the trajectory time (). Lower 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 per axis.
5 Experimental Results
In this section, we apply our proposed datadriven optimization methodology from Section 4 to the experimental apparatus in 2.2, using the test cases in 2.3. Our experiments demonstrate that the methodology is capable of improving system performance relative to the baseline – the trajectory obtained with the lowlevel controller following a not optimized input trajectory – by tracing desired geometries faster and more precisely.
5.1 Individual trajectories
We first focus on the letter “r” as a case study and consider a scenario with an acceleration limit of ; the individual results for the other letters in the test case are comparable and will be discussed collectively in Section 5.2. Applying our method results in an optimized input trajectory with a total time of . The experimental result of this new input trajectory is compared with a baseline performed in the same total time. For this the original, nonoptimized shape (in this case the letter “r”) is sampled at a constant progress speed
(12) 
where is the total path length, which for a closed shape corresponds to the perimeter.
Figure 11 compares the output and input trajectories for the optimized and baseline cases in the plane while Figure 12 displays the same data as a function of progress.
Analyzing Figures 11 and 12, we see that the optimized output effectively stays within the selected while the baseline output deviates by more than . The optimized output trajectory speeds up in areas with low curvature and slows down near corners or other intricate features. While the baseline input trajectory path speed is constant, the speed of the projection of the output onto the target geometry shows peaks around the corners. This is an artifact of the projection, given the way the nonoptimized output cuts through corners, as shown in the detail of Figure 11. Figure 12 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 11 illustrates this behaviour in the plane.
Figures 1315 display the same information for the airfoil test geometry with a maximum acceleration of . The baseline and optimized trajectories complete the part in . 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 14 and 15 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:

It reduces the speed of the input trajectory near intricate geometric features and increases it in areas with low curvature.

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 datadriven models to push the limits of system performance and adapt the input trajectory to the machine capabilities, for the selected maximum acceleration value.
5.2 Precisionaccuracy Tradeoff
In precision motion systems there is a fundamental tradeoff 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 tradeoff.
In our formulation (10)(11), the tradeoff between speed and accuracy is controlled by the parameter which limits the maximum acceleration of the input trajectory. We generated tradeoff curves for both the letters and airfoil test geometries by conducting experiments for different values of this parameter between and for the letters and between and for the airfoil, and tested each resulting input trajectory experimentally. Similarly to the two individual test cases presented above, the nonoptimized input trajectories are subsequently run with the same total time as the one found for the optimized input trajectories.
We use the normalized and (rms) norms as well as the norm to quantify the deviation between the output and target geometry , both in simulation and experimentally
(13) 
(14) 
(15) 
The results are shown in Figures 16 and 17. In all cases the optimized experimental output results in significantly more precise trajectories for a given part completion time. Moreover, the simulated deviation is lower than the experimental one due to model mismatch. The nonoptimized experimental output shows lower deviation than the simulation predicts in most cases. Model predictions are distributed approximately symmetrically around the experimental output as seen in Figure 8. However the error is measured between a) the prediction and the target geometry for the simulated case and b) the experimental output and target geometry for the experimental case. For a) the prediction error can be seen as the sum of the expected error, in this case measured experimentally and the model error, while in b) is only the experimental error. For acceleration values higher than , as seen in the left side of Figure 17, the model prediction exhibits high deviations. Since the model used for simulation is the same as the one used for optimization, in this region the designed trajectory has over exploited the model, resulting in an overly optimistic simulated deviation.
As seen in Figure 17, precision begins to degrade rapidly for high accelerations starting at around . This occurs because the training data for the highfidelity model does not contain sufficient data at those acceleration levels, leading to poor model accuracy, as shown in Figure 9. In contrast, the predicted (simulated) precision grows more slowly as it does not account for model mismatch.
Figures 16 and 17 also show that for the same deviation values the optimized trajectories take less time to complete. For example, in the letters test case the nonoptimized version takes and achieves a deviation of , while the optimized version takes only 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 .
In Table 3 we show the results for different test cases, each at two different scenarios. For all cases our method improves system performance in , and norms.
Shape  time  time  

–  %  %  %  %  %  %  
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 
6 Conclusion
In this paper we proposed a method that improves the precision vs. productivity tradeoff 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 tradeoff 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 invaluable support with the operation of the experimental apparatus.