Data-driven Reference Trajectory Optimization for Precision Motion Systems

by   Samuel Balula, et al.
ETH Zurich

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


From Low to High Order Motion Planners: Safe Robot Navigation using Motion Prediction and Reference Governor

Safe navigation around obstacles is a fundamental challenge for highly d...

Data-driven learning of nonlocal models: from high-fidelity simulations to constitutive laws

We show that machine learning can improve the accuracy of simulations of...

A Holistic Motion Planning and Control Solution to Challenge a Professional Racecar Driver

We present a holistically designed three layer control architecture capa...

Transfer Learning for High-Precision Trajectory Tracking Through L_1 Adaptive Feedback and Iterative Learning

Robust and adaptive control strategies are needed when robots or automat...

Flexible-Joint Manipulator Trajectory Tracking with Learned Two-Stage Model employing One-Step Future Prediction

Flexible-joint manipulators are frequently used for increased safety dur...

Adaptive Feedforward Control For Reset Feedback Control Systems – Application in Precision Motion Control

This paper presents a novel adaptive feedforward controller design for r...

No More Differentiator in PID:Development of Nonlinear Lead for Precision Mechatronics

Industrial PID consists of three elements: Lag (integrator), Lead (Diffe...

1 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 [kim1998high], micro-assembly [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 trade-off between productivity and accuracy and, in practice, the control strategy is often the determining factor in 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 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:

  1. does not require modifications to the existing hardware and low-level controllers, making it suitable for application in production machines;

  2. does not require a physics-based model of the machine, relying exclusively on experimental data;

  3. can be applied to individual machines with diverse technical characteristics;

  4. 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 data-driven 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 non-homogeneous machines. This makes it more important to exploit data to ensure that machines can execute designs with minimal calibration.

1.1 State-of-the-art

1.1.1 Real time control

One method to improve the performance of PMS is by re-designing 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 inner-loop controllers.

1.1.2 Learning-based control

Another common approach in manufacturing and precision motion systems is learning-based 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 Run-to-Run 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 Pre-compensation

Pre-compensation 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 look-ahead or constraints, motivating the use of offline MPCC algorithms based on physics-based 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 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.

1.2 Contributions

In this paper, we propose a novel pre-compensation method based on data-driven modelling and trajectory optimization to improve simultaneously the productivity and performance of PMS. Our contributions are three-fold:

  1. A data-driven framework for PMS modelling and identification for control, including design of experiments for collecting training data, and machine learning architecture in a feed-forward form.

  2. A novel contour control pre-compensation method featuring time-optimal trajectory design and motion constraints.

  3. Experimentally validated improvement over the existing controller across the precision/productivity trade-off curve.

Our approach does not require any modification of the design or control software of the system (D1) and is fully data-driven (D2). Moreover, it enables adaptation to changing conditions or system degradation through retraining the underlying data-driven model (D3). 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. (D4). We build upon our previous work [balula2020reference] 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.

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 two-dimensional (2D) workspace . The position of the tooltip as a function of time is denoted by


The ultimate goal in precision motion planning is to have the tooltip trace a spatial target geometry


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

Figure 1: An example target geometry, featuring a segment of an Archimedean spiral, described in polar coordinates by with and . The two plots show alternative representations of the same shape. In the left panel it is shown in a plot, parametric in the path progress , while in the right panel each axis is a function of .

Most industrial precision motion systems come with a built-in low-level 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 low-level 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


where the function encapsulates the closed-loop dynamics of the machine, and

is a noise term, with zero-mean and a typical standard deviation of

, which quantifies the repeatibility of the system.

Ideally, the low-level 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 low-level 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 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 through preliminary experiments, then inverting it using optimization.

Figure 2: A low-level controller tracks the input trajectory , producing the output trajectory . The deviation measures the distance from to the target geometry as a function of the path progress . We desire that the deviation is smaller than the tolerance at all points.

Formally, let be the deviation function, defined as the projection of the output into the target geometry


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 trade-off 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 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.

2.2 Experimental Setup

(a) Experimental apparatus
(b) Letters test case
(c) airfoil test case
Figure 3: The experimental setup shown in panel 2(a)

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 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 . 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 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
Low-level control loop frequency
Repeatability (standard deviation)
Table 1: Software-limited values and operational specifications of the experimental setup.

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 high-lift 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 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.

2.4 Architecture overview

Figure 4: Scheme of the approach used to optimize the input trajectory. In the first stage optimization (10) a minimum time trajectory for tracking the target path is computed, using the continuous time linear model. The speed and acceleration profile from the first stage is used to sample the target geometry at constant sample rate of , creating a surrogate target discrete time trajectory . In the second stage optimization (11) a discrete time ANN-based nonlinear model at a sampling frequency of is used to find the best input trajectory to track , using the the first stage solution as initial guess. The parameter, found in the constraints (4.1) and (4.2) is singled out since the maximum acceleration of the output trajectory is found to be decisive for how much time it takes to track the target geometry. The optimized input trajectory parameterized in discrete time is given to the experimental apparatus, and the output trajectory recorded by the axis quadrature encoders at a sample rate of .

We approach the trajectory design problem in several stages. First, we design experiments to generate informative input-output data for the unknown mapping 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 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 two-stage optimization-based 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 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.

3 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 provided by the lower-fidelity linear model.

3.1 Low-fidelity linear model

The closed-loop 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


The model (5) uses a non-physical 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 High-fidelity 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 trade-off 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 2D-positioning 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


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.

Figure 5: Neural network structure. Each layer is fully connected to the next. The input layer has neurons, there are hidden layers, with neurons, and the output layer has

neuron. The activation function of the hidden layers is a

function with slope , and a bias term. There are total of weights to be adjusted in the training. Two independent networks with identical structure were trained, one for each axis.

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

Figure 6: 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 , shown in yellow background.

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 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 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 training-optimisation-retraining cycle can be executed recursively until the model prediction accuracy is acceptable.

Figure 7: Flow of data used to train the linear and nonlinear models. First experimental data is collected for randomly generated trajectories and regular shapes. This dataset is used to fit and train the linear and a first generation nonlinear models. The models are used to compute optimized input trajectories for the set of regular shapes using the method of Section 4 the additional data is used to extend the training the nonlinear model. This is done because the optimizer has a tendency to exploit regions of the nonlinear model that are not well fitted to reality.

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 /
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
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 2: Prediction error mean and standard deviation of the linear and non-linear models. Both models are evaluated in the Letters test case after optimisation, with different maximum acceleration values . The linear model is also evaluated on its training data set, and the non-linear model on its validation data set.
(a) Validation data
(b) Letters test case with
(c) Letters test case with
Figure 8: Histogram of the prediction error of the non-linear model. In panel (a) the model is evaluated in the validation data set, as exemplified in Figure 6. In panels (b) and (c) the model is evaluated in the Letters test case, with the input trajectory optimised with up to .
Figure 9: Prediction error of the nonlinear model as a function of the maximum acceleration of the input trajectory. The error is evaluated for the letters and airfoil test cases. Above the training data do not contain enough information to train the model properly.

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 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 of the trajectory in a computationally tractable manner. This is accomplished by using a variable time-discretization in the first stage, enabling inclusion of as a decision variable. The second stage then fixes the duration 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 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 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, as shown schematically in Figure 7.

4.1 First-stage 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


We use the identified linear model (5) to approximate the process dynamics. As our approach involves a fixed spatial discretization , the time-discretization 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 non-uniform 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 Runge-Kutta (RK4)[aboanber2008generalized] time-marching method to (5)


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


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


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 .

Figure 10: For every point of the target geometry a local coordinate frame is placed. is tangent to the target geometry at and pointing in the direction of increasing whereas is orthogonal to pointing to the port side. is the angle between the global axis and . The maximum allowed deviation is computed from . is the path distance between points and .

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^N-1 Δτ_k + 1N-2 ∑_k=2^N-1 ∥∂^2 ρ_k∥_2^2 z_k+1= f_lm(z_k, ρ_k, Δτ_k), k = 1, …, N-1 γ_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, …, N-1 ∥ ∂^2 γ_k ∥_∞≤  a_max, k = 2, …, N-1 , 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 non-negative, (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 component-wise limits on velocity and acceleration of the output.

4.2 Second-stage 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 high-fidelity 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


, 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 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 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 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 trade-off 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 second-stage (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) i9-9900K CPU and of RAM.

For the problem sizes solved in this paper, the first-stage optimization takes approximately 1 minute to complete (), and the second-stage 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 data-driven 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 low-level 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, non-optimized shape (in this case the letter “r”) is sampled at a constant progress speed


where is the total path length, which for a closed shape corresponds to the perimeter.

Figure 11: Letter “r” of the ETH Zürich Logo. Optimized with , 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 12: Letter “r” of the ETH Zürich Logo optimized with . The letters A-G refer to the points labelled in Figure 11. Top panel shows 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 panel shows the velocity of the optimized and non-optimized outputs projected onto the target geometry.

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 non-optimized 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 13-15 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.

Figure 13: Experimental results for the airfoil test case optimized with . L denotes the leading and T the trailing edge of the airfoil. The shape is traversed in the clockwise direction starting from the leading edge.
Figure 14: A close-up of the trailing edge of the airfoil test case with . The optimized trajectory is able to track the desired geometry more precisely than the baseline, though still not within the tolerance band. The optimized input trajectory is not constrained by the tolerance band.
Figure 15: Experimental airfoil test case with . The letters L and T refer to the leading and trailing edges labelled in Figure 13. The interpretation of the panels in analogous to Figure 12.

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.

5.2 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 (10)-(11), the trade-off between speed and accuracy is controlled by the parameter 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 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 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 and (rms) norms as well as the norm to quantify the deviation between the output and target geometry , both in simulation and experimentally


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 non-optimized 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 high-fidelity 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 non-optimized 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.

Figure 16: Precision-speed trade-off curve for the concatenation of the letters “u”, “r”, “c”, and “h” from the ETH Zürich logo test case. The top panel shows the normalised deviation of the optimized and non-optimized outputs from the target geometry. The bottom panel shows the total time to perform the trajectory and the 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 letters for a given acceleration, then retrieve the deviation for the corresponding time from the top plot.
Figure 17: Precision-speed trade-off curve for the airfoil test case. The interpretation of the plots is similar to Figure 16
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
Table 3: 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 .

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


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