High Performance Motion Tracking Control - RPI - Rensselaer ...
High Performance Motion Tracking Control Benjamin Potsaid and John T. Wen Center for Automation Technologies Rensselaer Polytechnic Institute Troy, NY 12180. Emails: {potsaid,wen}@cat.rpi.edu
Abstract— Motion control in electronic manufacturing requires high speed and precision. Due to the large production volume, even a small improvement in the motion control performance can result in significant economic gains. This paper presents an experimental study of the feedforward controller design for a pointto-point motion control system. We systematically apply model identification, inverse dynamics control, iterative refinement (to address system nonlinearities), and adaptive least-mean-square to achieve high speed trajectory tracking. Experimental results are provided to show the efficacy of the proposed approach.
I. I NTRODUCTION Due to the volume and performance requirement, motion control in electronic manufacturing requires high speed and precision. To achieve this goal, feedforward compensation as well as feedback control are needed. In this paper, we present a systematic approach to improve the performance of point-to-point motion of a positioning system by designing a feedforward controller after a feedback controller has already been implemented. Since only a single positioner is considered, the system is single-input/single-output (SISO). The overall experimental setup is shown in Figure 1. The position is measured through a high resolution joint encoder. The data acquisition and closed loop control system is capable operating at sampling frequency up to 200KHz. Velocity and acceleration constraints are imposed on the motion to avoid internal saturation within the feedback controller and servo amplifier. Based on our experimental testbed and requirement in electronic manufacturing, we choose (for the vmax = 5m/sec, amax = 1000g (for the equivalent laser beam location). Our overall approach to the feedforward control design is summarized below: • Verify system linearity for small amplitude inputs. • Obtain the experimental frequency response. • Identify a transfer function model based on the experimental frequency response. • Use the identified model to design an inverse dynamics feedforward controller. • Iteratively refine the feedforward control signal based on the output tracking error and using the identified model to generate a gradient descent direction. • Use the iterative refinement result to train a finite impulse response (FIR) filter to generate a corrective input. • Adjust the FIR coefficients using the output tracking error based on a gradient type algorithm.
Fig. 1.
Positioning System Experimental Testbed
Experimental results are included throughout to show the effectiveness of this approach. II. I DENTIFICATION Linear time invariant (LTI) model identification based on input/output responses may be performed in either time or frequency domains. The time domain approach is adversely affected by high frequency noise. Therefore we decided to use the frequency domain subspace identification method [1], [2]. The experimental transfer function is first obtained by a sine sweep (with input amplitude chosen small enough to avoid saturation but large enough to increase signal-to-noise ratio). The maximum sweep frequency is 20KHz to avoid the noise frequency spectrum of around 50KHz. In system identification, it is important to consider possible input/output delays. If the experimental system contains a pure delay, direct application of LTI identification will often approximate the delay with non-minimum phase zeros and additional poles (Pade approximation). The approximation can be avoided by including a pure time delay in the identification procedure by time shifting the output signal relative to the input signal before performing the identification. Including the delay in this manner can result in identified models of lower order, better agreement between simulated and experimental
Fig. 3. 500µm Fig. 2. Frequency Gain/Phase and Time Step Response Comparison between Experimental Data and Identified Model with 16-sample Delay
responses (no artificial undershoot), and significantly less nonminimum phase behavior in the identified model. For the inverse dynamics controller, the non-minimum phase zeros can compromise the tracking performance, while the pure delay simply leads to a time shift in the response. Our identified plant model exhibited significant non-minimum phase zeros when identified with no delay included, while including a delay of 16 samples (at 5 µsec sampling period) showed minimal non-minimum phase zero contribution and good agreement between simulation and experiment. For the final identified model used for the control design, we use a 16-sample delay, which results in an 16th order identified model. The gain and phase comparison between the experiment frequency and time step responses are shown in Fig. 2. III. I NVERSE DYNAMICS C ONTROL AND I TERATIVE R EFINEMENT The next step is to use the identified model, G, to construct ˆ † , by replacing the unstable zeros an inverse dynamics filter, G by the stable mirror images and inverting the transfer function. We then apply the inverse dynamics filter to the desired ˆ † ydes . output ydes to generate the command input: u = G A trajectory generator is used to generate ydes (t) based on the position, velocity, and acceleration constraints. To avoid unbounded jerk, we use a half-sine profile for acceleration and deceleration. The experimental result of the 500µm move is shown in Figure 3. Large tracking errors near the entry into and inside the settling zone are observed. Clearly, using inverse dynamics filter alone does not meet the desired performance specification. The discrepancy between the experimental and simulated responses is due to the mismatch between the model and the physical system. To correct for this mismatch, we modify the ˆ † ydes + ∆u where ∆u is obtained using the input to u = G iterative refinement algorithm based on the output tracking error (described in Appendix I). Using the complete trajectory tracking error to iteratively update the command input is known as iterative learning control (ILC). This concept was originally introduced in [3] for robot tracking control. Most
Simulated vs. experimental results using inverse dynamics filter:
ILC algorithms (proportional, Newton, Secant) update the input at time t based on the previous input and tracking error also at time t [4], [5]. In [6], [7], the input is updated based on the gradient of the complete trajectory tracking error with respect to the input. A similar idea is also proposed in an ILC general setting in [8] using the Newton algorithm. For this study, we apply the gradient descent approach with the nominal LTI model, G, as the approximate gradient. The basic algorithm is summarized below: Algorithm 1: Given y ∗ := {y ∗ (ti ) : i ∈ 0, 1, . . . , N } and u0 := {u0 (ti ) : i ∈ 0, 1, . . . , N }. Set u = u0 . 1) Apply u to the physical system and obtain the output sequence y := {y(ti ) : i ∈ 0, 1, . . . , N }. 2) Update u by adding a corrective term ∆u = −αG∗ (y − y ∗ )
(1)
where G∗ is the adjoint of G, and α may be set as a sufficiently small constant or found by using a line search (which would require additional runs). 3) Iterate until y − y ∗ or k∆uk becomes sufficiently small. The key step in the above algorithm is the update equation (1). Let the state space parameters of G be (A, B, C, D). The adjoint G∗ is given by (−AT , −C T , B T , DT ) but it must propagate backward in time from the zero state. To implement G∗ ∆y, we first reverse ∆y backwards in time, filter it forward in time through the filter (AT , C T , B T , DT ), and then reverse the result in time again. The result of iterative refinement for the move length 500µm is shown in Figure 4. The desired trajectory is duplicated almost exactly after 6 iterations. IV. FIR A PPROXIMATION Excellent tracking performance is obtained by using iterative refinement, but the update process is non-causal, i.e., the output error of a complete run is needed to update the input at any given time. To allow for real-time implementation, we use the results from iterative refinement to train a filter mapping ydes to the corrective input, ∆u. Many filter parameterizations are possible; we decide on the FIR filter structure due to its
Fig. 4. Experimental results of applying iterative refinement to 500µm move
efficient real-time implementation and the ease for finding the filter coefficients. We parameterize the corrective input by the following FIR filter: ∗ ∗ ∗ ∆uk = w1 yk+n + w2 yk+n + . . . + wn yk−n , 1 1 −1 2 +1
(2)
where n = n1 + n2 is the order of the filter, n1 is the lookahead horizon. The filter coefficients, {wi , i = 1, . . . , n}, may be obtained through a least square fit to the (ydes , ∆u) data obtained from iterative refinement. The idea is to choose a sufficiently rich set of ydes so that the FIR filter is applicable to any specific ydes of interest. First re-write the FIR equation as y∗ yk∗1 −n2 +1 k1 +n1 w1 ∆uk1 .. . .. .. . . . . . = . , . .. . ∗ . . y k +n . . 1 1 . .. ∗ wn ∆uk2 y k2 +n1 | {z } | {z } | {z } U
Y
Fig. 5. Command input comparison between inverse dynamics, iterative refinement, and FIR filter: 500µm move
w
where (k1 , k2 ) (range of fit), n (order of FIR), and n1 (lookahead horizon) are parameters to be chosen. Note that Y is a Toeplitz matrix. The least square solution of filter coefficients w is given by w = Y †U (3)
Fig. 6. Output comparison between iterative refinement and FIR filter: 500µm move
generated by the iterative refinement. This presents two drawbacks: 1. As seen in the previous section, the FIR fit is not perfect which results in a larger output tracking error; 2. When the system changes over time, the performance of the FIR filter will likely degrade. To address these issues, we apply the adaptive least mean square (LMS) controller [9] to adjust the FIR filter coefficients based on the output tracking error. This approach is equivalent to the standard gradient parameter estimation algorithm as shown in Appendix II. The overall controller architecture is given in Figure 7.
where Y † is the Moore-Penrose pseudo-inverse of Y . After some experimentation, we decided on a 100-tap (n = 100) FIR filter with 1 look-ahead step (n1 = 1). For real-time implementation, the FIR runs at 20KHz (while the inverse dynamics filter runs at 200KHz). The iterative refinement data for move lengths 200µm to 700µm are used to fit for the FIR filter coefficients. The comparison between ∆u for the iterative refinement and FIR filter for the 500µm move case is shown in Figure 5, where the FIR reproduces the iterative refinement inputs very well. The output responses for the 500µm move is shown in Figure 6. Clearly, the use of FIR does result in slightly larger errors, especially in terms of the steady state values. V. A DAPTIVE FIR F ILTER The FIR approach described in the previous section requires fitting the filter coefficients to the ideal input/output response
Fig. 7.
Block diagram of adaptive LMS filter
We have experimentally determined the best adaptive LMS
architecture involves three FIR filters, with the corresponding inputs being the desired position, velocity, and acceleration: ∆uk = wxT y desk+n
x1
+ wvT y˙ desk+n
v1
+ waT y¨desk+n
(4) a1
where y desj
:= {ydes (ti+j ) : i = 0, 1, . . . , Nx }
y˙ desj
:= {y˙ des (ti+j ) : i = 0, 1, . . . , Nv }
y¨desj
:= {¨ ydes (ti+j ) : i = 0, 1, . . . , Na } .
the first two equations in (5) with wv initially set at zero, wa remains zero). The output tracking result is shown in Figure 10. The result is quite striking: the peak tracking error in the high velocity range is reduced by about 2/3! This is a direct result of using wv to reduce the velocity tracking error. The steady state error remains small.
The coefficient update using adaptive LMS is given by δwx δwv δwa
= −x (y − ydes )Gd y des = −v (y − ydes )Gd y˙ des = −a (y − ydes )Gd y¨des ,
(5)
where Gd is the identified discrete LTI system. Through some initial experimentation, we decide to use 100 taps for the position portion of FIR, and 4 taps each for the velocity and acceleration portion. The FIR filter and adaptive LMS are both implemented at 20KHz. The constant coefficient FIR (using the iterative refinement result to find the coefficients) performance is shown in Figure 8 (only wx is trained in (4), wv and wa are both zero). Note the presence of steady error as well as large tracking error in the high velocity region.
Fig. 8.
Lastly, we use all three portions (position, velocity, acceleration) of the FIR filter in (4) and update them using the adaptive LMS rule (5). The result is shown in Figure 11. There is some improvement in the 1000µm to 0µm transition region, but the improvement is inconclusive in other high acceleration regions. This may be due to the fact that the tracking error is already quite small (less than 2µm throughout the move). Thus the coefficient update is very slow. A possible reason is the presence of noise which limits the amount of achievable tracking performance.
Response with constant coefficient FIR filter
We first apply adaptive LMS to only the position part of the filter in (4) (i.e., wx is updated by the first equation in (5), wv and wa remain zero). The result after 330 repeated runs (to allow the coefficients to reach steady state) is shown in Figure 9. The steady state error is reduced since it strongly affects the update of wx . However, the transient error in the high velocity region remains.
Fig. 9.
Fig. 10. Response with adaptive FIR filter using position and velocity errors
Response with adaptive FIR filter using position error only
We next apply adaptive LMS to the position and velocity portions of the FIR filter (4) (i.e., wx and wv are updated by
Fig. 11. Response with adaptive FIR filter using position, velocity, and acceleration errors
VI. R ANDOM M OVE R ESULT To compare the performance of various approaches described so far, we use a random motion profile as a test case. An FIR filter is trained using the results of move lengths 200µm to 700µm at 100µm intervals. The desired trajectory is shown in Figure 12. The comparison between the tracking errors for the inverse dynamics filter and fixed-coefficnet FIR filter are shown in Figure 13. As expected, the tracking error is significantly and uniformly reduced in the FIR case. The addition of the adaptive LMS filters further improves the response, as shown in Figure 14. The adaptive filter clearly shows an improvement over the FIR filter with the xv and xva filters demonstrating the best performance. The quantitative comparison between these different controllers is summarized in Table I. Based on these results, we recommend the xv filter as the most promising candidate.
all units in µm inv dyn FIR adaptive x adaptive xv adaptive xva
min -10.88 -3.77 -2.91 -2.07 -1.99
max 13.11 2.56 2.57 2.38 2.15
RMS 2.96 1.04 0.52 0.37 0.46
mean 0.29 -0.68 -0.05 -0.03 0.21
stdev 2.94 0.79 0.52 0.37 0.41
TABLE I P ERFORMANCE COMPARISON FOR A RANDOM
MOVE PROFILE
Fig. 13.
Fig. 12.
Tracking error for random moves: inverse dynamics vs. FIR filter
Desired output trajectory
VII. C ONCLUSION This paper presents the experimental results of developing high performance feedforward control of a positioning system. Frequency domain subspace identification method with input delays is applied to the experimental frequency response. The identified model is first used in the inverse dynamics control law. Further performance enhancement is obtained by using iterative refinement using the identified model to provide a descent direction. The iterative refinement results are next used to train an FIR filter. Finally, the FIR filter coefficients are updated by using an adaptive LMS scheme. Tracking error for a random motion profile shows significant tracking error reduction by using the adaptive LMS filter. ACKNOWLEDGMENT This material is based upon work supported in part by the National Science Foundation under Grant CMS-0301827, and in part by the Center for Automation Technologies under a block grant from the New York State Office of Science, Technology, and Academic Research.
Fig. 14. Tracking error for random moves: comparison between FIR, x, xv, and xva filters
[5] J.X. Xu and Y. Tan. On the P-type and Newton-type ILC schemes for dynamic systems with non-affine-in-input factors. Automatica, 38(7):1237–1242, July 2002. [6] W. Cheng and J.T. Wen. A class of learning controllers with application to tracking control of a flexible beam. In Proc. 1993 IEEE Conference on Robotics and Automation, pages 411–416, Atlanta, GA, May 1993. [7] D.M. Gorinevsky. Direct learning of feedforward control for manipulator path tracking. In Proceedings of the 1992 IEEE International Symposium on Intelligent Control, pages 42–47, Glasgow, UK, August 1992. [8] K.E. Avrachenkov. Iterative learning control based on quasi-newton methods. In Conference on Decision and Control, pages 170–174, Tampa, FL, December 1998. [9] S. Haykin. Adaptive Filter Theory. Prentice-Hall, 3 edition, 1996. [10] G.C. Goodwin and K.S. Sin. Adaptive Filtering Prediction and Control. Prentice-Hall, 1984.
A PPENDIX I I TERATIVE R EFINEMENT A LGORITHM
R EFERENCES [1] T. McKelvey, H. Akcay, and L. Ljung. Subspace-based multivariable system identification from frequency response data. IEEE Trans. Automat. Contr., 41(7):960–979, July 1996. [2] P. Van Overschee and B. De Moor. Coninuous-time frequency domain subspace system identification. signal Processing, Special Issue on Subspace Methods, Part: System Identification, 52:179–194, 1996. [3] S. Arimoto. Mathematical theory of learning with applications to robot control. In K.S. Narendra, editor, Adaptive and Learning Systems: Theory and Applications. Plenum Press, New York, 1985. [4] J.X. Xu and Q. Ji. New ILC algorithms with improved convergence for a class of non-affine functions. In Conference on Decision and Control, pages 660–665, Tampa, FL, December 1998.
Consider a nonlinear dynamical system H : Lm 2 [0, T ] → that is initially at rest:
Lp2 [0, T ]
Define V =
y = H(u).
(6)
1 2 ky − ydes kLp . 2 2
(7)
Note that positive time dependent weighting may be used in the above norm with obvious modifications in the subsequent
update law. The variation of V is δV = hy − ydes , δyiLp
(8)
2
where h·, ·iL2 denotes the L2 inner product. From (6), δy is δy = ∇u Hδu
(9)
where ∇u H is the Fr´echet derivative of H. Suppose H is given by the following state space representation: x˙ = f (x, u), x(0) = x0 ;
y = h(x, u).
(10)
Eq. (19) needs to be integrated backward in time. Through a change of variable, z1 (t) = z(T − t), we can instead use the forward solution. The adjoint system in this case is then (with z1 (0) = 0 and δv1 (t) := δv(T − t)): z˙1 (t) = AT (T − t)z1 (t) + C T (T − t)δv1 (t) w(T − t) = B T (T − t)z1 (t) + DT (T − t)δv1 (t).(20) This can be implemented by reversing δv backward in time, filtering it through (20), and then reversing the result in time again.
Then ∇u H is the time varying linearized system about x and u (x is the state trajectory corresponding to u):
A PPENDIX II A DAPTIVE FIR A LGORITHM
∂f (x, u) ∂f (x, u) δx + δu, δx(0) = 0 ∂x ∂u ∂h(x, u) ∂h(x, u) δy = δx + δu, (11) ∂x ∂u where ∇u H is the operator that maps δu to δy. The gradient descent algorithm choose δu based on ∇u H ∗ :
where ydes is the desired output, w is the constant weighting vector, and u∗ is a feedforward control that depends on ydes only. Then the output is
δu = −α∇u H ∗ (y − ydes ).
y = (Hy des )T w + Hu∗ .
δ x˙ =
δy = −α(∇u H)(∇u H ∗ (y − ydes ))
2
δV = −α k∇u H ∗ (y − ydes )kLp . 2
uk+1 = uk − αk ∇u H (uk )(yk − ydes ),
ydes = (Hy des )T w∗ + Hu∗ .
(14)
Our goal is to choose an update law for w so that y → ydes even though w∗ is unknown. To simplify notation, let a = Hy des , and write (23) in discrete time as
(15)
where αk is found through a line search to ensure Vk+1 (V with uk+1 ) is strictly smaller than Vk . Now we consider the computation of ∇u H ∗ . Write ∇u H as a linear time varying system (same as (11)): δ x˙ = A(t)δx + B(t)δu, δy = C(t)δx + D(t)δu.
(16)
Let Φ(t, τ ) be the state transition matrix from τ to t. The output can now be written as Z t δy(t) = C(t) Φ(t, τ )B(τ )δu(τ ) dτ + D(t)δu(t). 0 ∗
To find ∇u H , we use the defining relationship hδv, ∇u Hδui = h∇u H ∗ δv, δui for all δv and δu in L2 [0, T ]. We then obtain Z T ∇u H ∗ δv(τ ) = B T (τ )ΦT (t, τ )C T (t)δv(t) dt+DT (τ )δv(τ ). τ
(17)
z(τ ) =
T
ΦT (t, τ )C T (t)δv(t) dt.
(22)
(13)
Provided that ∇u H is onto, V will converge to zero. For a discrete implementation, (12) is replaced by ∗
(21)
Assume that an optimal w∗ exists such that
which implies
Z
u = y Tdes w + u∗
(12)
If this can be implemented exactly, then from (9)
Define
Consider a SISO LTI system, y = Hu. Choose
(18)
τ
The state space representation of ∇u H ∗ is then given by the mapping from δv to w as follows: z(t) ˙ = −AT (t)z(t) − C T (t)δv(t), z(T ) = 0 w(t) = B T (t)z(t) + DT (t)δv(t). (19)
ydes (t) − Hu∗ (t) = aT (t − 1)w∗ ,
(23)
(24)
where the left hand side and a(t − 1) are known. The problem then becomes a standard least square parameter estimation problem [10, Sec.3.3]. Various algorithms, ranging from projection to least square with covariance resetting may be used. To reduce the real-time computation load, we use the gradient update law: w(t) = w(t − 1) − a(t − 1)(y(t) − y des (t)).
(25)
To analyze the convergence of this scheme, first compute the evolution of the parameter estimation error (δw := w − w∗ ): 2
2
kδw(t)k − kδw(t − 1)k = (−2 + aT (t − 1)a(t − 1))e2 (t). (26) where e(t) := ydes (t) − y(t). To ensure that the parameter error does not increase, we need choose sufficiently small (a(t) is assumed uniformly bounded). It follows that ∞ X
(aT (t − 1)a(t − 1) − 2)e2 (t) < ∞
(27)
t=1
which implies (aT (t − 1)a(t − 1) − 2)e2 (t) → 0 as t → ∞. If aT (t − 1)a(t − 1) does not converge to 2, then e(t) converges to zero as desired. Other type of parameter estimation algorithms with more relaxed convergence condition may be used (e.g., projection, least square with covariance reset, etc.). However, the realtime computation load would be more demanding, possibly preventing their implementation.
Abstract— Motion control in electronic manufacturing requires high speed and precision. Due to the large production volume, even a small improvement in the motion control performance can result in significant economic gains. This paper presents an experimental study of the feedforward controller design for a pointto-point motion control system. We systematically apply model identification, inverse dynamics control, iterative refinement (to address system nonlinearities), and adaptive least-mean-square to achieve high speed trajectory tracking. Experimental results are provided to show the efficacy of the proposed approach.
I. I NTRODUCTION Due to the volume and performance requirement, motion control in electronic manufacturing requires high speed and precision. To achieve this goal, feedforward compensation as well as feedback control are needed. In this paper, we present a systematic approach to improve the performance of point-to-point motion of a positioning system by designing a feedforward controller after a feedback controller has already been implemented. Since only a single positioner is considered, the system is single-input/single-output (SISO). The overall experimental setup is shown in Figure 1. The position is measured through a high resolution joint encoder. The data acquisition and closed loop control system is capable operating at sampling frequency up to 200KHz. Velocity and acceleration constraints are imposed on the motion to avoid internal saturation within the feedback controller and servo amplifier. Based on our experimental testbed and requirement in electronic manufacturing, we choose (for the vmax = 5m/sec, amax = 1000g (for the equivalent laser beam location). Our overall approach to the feedforward control design is summarized below: • Verify system linearity for small amplitude inputs. • Obtain the experimental frequency response. • Identify a transfer function model based on the experimental frequency response. • Use the identified model to design an inverse dynamics feedforward controller. • Iteratively refine the feedforward control signal based on the output tracking error and using the identified model to generate a gradient descent direction. • Use the iterative refinement result to train a finite impulse response (FIR) filter to generate a corrective input. • Adjust the FIR coefficients using the output tracking error based on a gradient type algorithm.
Fig. 1.
Positioning System Experimental Testbed
Experimental results are included throughout to show the effectiveness of this approach. II. I DENTIFICATION Linear time invariant (LTI) model identification based on input/output responses may be performed in either time or frequency domains. The time domain approach is adversely affected by high frequency noise. Therefore we decided to use the frequency domain subspace identification method [1], [2]. The experimental transfer function is first obtained by a sine sweep (with input amplitude chosen small enough to avoid saturation but large enough to increase signal-to-noise ratio). The maximum sweep frequency is 20KHz to avoid the noise frequency spectrum of around 50KHz. In system identification, it is important to consider possible input/output delays. If the experimental system contains a pure delay, direct application of LTI identification will often approximate the delay with non-minimum phase zeros and additional poles (Pade approximation). The approximation can be avoided by including a pure time delay in the identification procedure by time shifting the output signal relative to the input signal before performing the identification. Including the delay in this manner can result in identified models of lower order, better agreement between simulated and experimental
Fig. 3. 500µm Fig. 2. Frequency Gain/Phase and Time Step Response Comparison between Experimental Data and Identified Model with 16-sample Delay
responses (no artificial undershoot), and significantly less nonminimum phase behavior in the identified model. For the inverse dynamics controller, the non-minimum phase zeros can compromise the tracking performance, while the pure delay simply leads to a time shift in the response. Our identified plant model exhibited significant non-minimum phase zeros when identified with no delay included, while including a delay of 16 samples (at 5 µsec sampling period) showed minimal non-minimum phase zero contribution and good agreement between simulation and experiment. For the final identified model used for the control design, we use a 16-sample delay, which results in an 16th order identified model. The gain and phase comparison between the experiment frequency and time step responses are shown in Fig. 2. III. I NVERSE DYNAMICS C ONTROL AND I TERATIVE R EFINEMENT The next step is to use the identified model, G, to construct ˆ † , by replacing the unstable zeros an inverse dynamics filter, G by the stable mirror images and inverting the transfer function. We then apply the inverse dynamics filter to the desired ˆ † ydes . output ydes to generate the command input: u = G A trajectory generator is used to generate ydes (t) based on the position, velocity, and acceleration constraints. To avoid unbounded jerk, we use a half-sine profile for acceleration and deceleration. The experimental result of the 500µm move is shown in Figure 3. Large tracking errors near the entry into and inside the settling zone are observed. Clearly, using inverse dynamics filter alone does not meet the desired performance specification. The discrepancy between the experimental and simulated responses is due to the mismatch between the model and the physical system. To correct for this mismatch, we modify the ˆ † ydes + ∆u where ∆u is obtained using the input to u = G iterative refinement algorithm based on the output tracking error (described in Appendix I). Using the complete trajectory tracking error to iteratively update the command input is known as iterative learning control (ILC). This concept was originally introduced in [3] for robot tracking control. Most
Simulated vs. experimental results using inverse dynamics filter:
ILC algorithms (proportional, Newton, Secant) update the input at time t based on the previous input and tracking error also at time t [4], [5]. In [6], [7], the input is updated based on the gradient of the complete trajectory tracking error with respect to the input. A similar idea is also proposed in an ILC general setting in [8] using the Newton algorithm. For this study, we apply the gradient descent approach with the nominal LTI model, G, as the approximate gradient. The basic algorithm is summarized below: Algorithm 1: Given y ∗ := {y ∗ (ti ) : i ∈ 0, 1, . . . , N } and u0 := {u0 (ti ) : i ∈ 0, 1, . . . , N }. Set u = u0 . 1) Apply u to the physical system and obtain the output sequence y := {y(ti ) : i ∈ 0, 1, . . . , N }. 2) Update u by adding a corrective term ∆u = −αG∗ (y − y ∗ )
(1)
where G∗ is the adjoint of G, and α may be set as a sufficiently small constant or found by using a line search (which would require additional runs). 3) Iterate until y − y ∗ or k∆uk becomes sufficiently small. The key step in the above algorithm is the update equation (1). Let the state space parameters of G be (A, B, C, D). The adjoint G∗ is given by (−AT , −C T , B T , DT ) but it must propagate backward in time from the zero state. To implement G∗ ∆y, we first reverse ∆y backwards in time, filter it forward in time through the filter (AT , C T , B T , DT ), and then reverse the result in time again. The result of iterative refinement for the move length 500µm is shown in Figure 4. The desired trajectory is duplicated almost exactly after 6 iterations. IV. FIR A PPROXIMATION Excellent tracking performance is obtained by using iterative refinement, but the update process is non-causal, i.e., the output error of a complete run is needed to update the input at any given time. To allow for real-time implementation, we use the results from iterative refinement to train a filter mapping ydes to the corrective input, ∆u. Many filter parameterizations are possible; we decide on the FIR filter structure due to its
Fig. 4. Experimental results of applying iterative refinement to 500µm move
efficient real-time implementation and the ease for finding the filter coefficients. We parameterize the corrective input by the following FIR filter: ∗ ∗ ∗ ∆uk = w1 yk+n + w2 yk+n + . . . + wn yk−n , 1 1 −1 2 +1
(2)
where n = n1 + n2 is the order of the filter, n1 is the lookahead horizon. The filter coefficients, {wi , i = 1, . . . , n}, may be obtained through a least square fit to the (ydes , ∆u) data obtained from iterative refinement. The idea is to choose a sufficiently rich set of ydes so that the FIR filter is applicable to any specific ydes of interest. First re-write the FIR equation as y∗ yk∗1 −n2 +1 k1 +n1 w1 ∆uk1 .. . .. .. . . . . . = . , . .. . ∗ . . y k +n . . 1 1 . .. ∗ wn ∆uk2 y k2 +n1 | {z } | {z } | {z } U
Y
Fig. 5. Command input comparison between inverse dynamics, iterative refinement, and FIR filter: 500µm move
w
where (k1 , k2 ) (range of fit), n (order of FIR), and n1 (lookahead horizon) are parameters to be chosen. Note that Y is a Toeplitz matrix. The least square solution of filter coefficients w is given by w = Y †U (3)
Fig. 6. Output comparison between iterative refinement and FIR filter: 500µm move
generated by the iterative refinement. This presents two drawbacks: 1. As seen in the previous section, the FIR fit is not perfect which results in a larger output tracking error; 2. When the system changes over time, the performance of the FIR filter will likely degrade. To address these issues, we apply the adaptive least mean square (LMS) controller [9] to adjust the FIR filter coefficients based on the output tracking error. This approach is equivalent to the standard gradient parameter estimation algorithm as shown in Appendix II. The overall controller architecture is given in Figure 7.
where Y † is the Moore-Penrose pseudo-inverse of Y . After some experimentation, we decided on a 100-tap (n = 100) FIR filter with 1 look-ahead step (n1 = 1). For real-time implementation, the FIR runs at 20KHz (while the inverse dynamics filter runs at 200KHz). The iterative refinement data for move lengths 200µm to 700µm are used to fit for the FIR filter coefficients. The comparison between ∆u for the iterative refinement and FIR filter for the 500µm move case is shown in Figure 5, where the FIR reproduces the iterative refinement inputs very well. The output responses for the 500µm move is shown in Figure 6. Clearly, the use of FIR does result in slightly larger errors, especially in terms of the steady state values. V. A DAPTIVE FIR F ILTER The FIR approach described in the previous section requires fitting the filter coefficients to the ideal input/output response
Fig. 7.
Block diagram of adaptive LMS filter
We have experimentally determined the best adaptive LMS
architecture involves three FIR filters, with the corresponding inputs being the desired position, velocity, and acceleration: ∆uk = wxT y desk+n
x1
+ wvT y˙ desk+n
v1
+ waT y¨desk+n
(4) a1
where y desj
:= {ydes (ti+j ) : i = 0, 1, . . . , Nx }
y˙ desj
:= {y˙ des (ti+j ) : i = 0, 1, . . . , Nv }
y¨desj
:= {¨ ydes (ti+j ) : i = 0, 1, . . . , Na } .
the first two equations in (5) with wv initially set at zero, wa remains zero). The output tracking result is shown in Figure 10. The result is quite striking: the peak tracking error in the high velocity range is reduced by about 2/3! This is a direct result of using wv to reduce the velocity tracking error. The steady state error remains small.
The coefficient update using adaptive LMS is given by δwx δwv δwa
= −x (y − ydes )Gd y des = −v (y − ydes )Gd y˙ des = −a (y − ydes )Gd y¨des ,
(5)
where Gd is the identified discrete LTI system. Through some initial experimentation, we decide to use 100 taps for the position portion of FIR, and 4 taps each for the velocity and acceleration portion. The FIR filter and adaptive LMS are both implemented at 20KHz. The constant coefficient FIR (using the iterative refinement result to find the coefficients) performance is shown in Figure 8 (only wx is trained in (4), wv and wa are both zero). Note the presence of steady error as well as large tracking error in the high velocity region.
Fig. 8.
Lastly, we use all three portions (position, velocity, acceleration) of the FIR filter in (4) and update them using the adaptive LMS rule (5). The result is shown in Figure 11. There is some improvement in the 1000µm to 0µm transition region, but the improvement is inconclusive in other high acceleration regions. This may be due to the fact that the tracking error is already quite small (less than 2µm throughout the move). Thus the coefficient update is very slow. A possible reason is the presence of noise which limits the amount of achievable tracking performance.
Response with constant coefficient FIR filter
We first apply adaptive LMS to only the position part of the filter in (4) (i.e., wx is updated by the first equation in (5), wv and wa remain zero). The result after 330 repeated runs (to allow the coefficients to reach steady state) is shown in Figure 9. The steady state error is reduced since it strongly affects the update of wx . However, the transient error in the high velocity region remains.
Fig. 9.
Fig. 10. Response with adaptive FIR filter using position and velocity errors
Response with adaptive FIR filter using position error only
We next apply adaptive LMS to the position and velocity portions of the FIR filter (4) (i.e., wx and wv are updated by
Fig. 11. Response with adaptive FIR filter using position, velocity, and acceleration errors
VI. R ANDOM M OVE R ESULT To compare the performance of various approaches described so far, we use a random motion profile as a test case. An FIR filter is trained using the results of move lengths 200µm to 700µm at 100µm intervals. The desired trajectory is shown in Figure 12. The comparison between the tracking errors for the inverse dynamics filter and fixed-coefficnet FIR filter are shown in Figure 13. As expected, the tracking error is significantly and uniformly reduced in the FIR case. The addition of the adaptive LMS filters further improves the response, as shown in Figure 14. The adaptive filter clearly shows an improvement over the FIR filter with the xv and xva filters demonstrating the best performance. The quantitative comparison between these different controllers is summarized in Table I. Based on these results, we recommend the xv filter as the most promising candidate.
all units in µm inv dyn FIR adaptive x adaptive xv adaptive xva
min -10.88 -3.77 -2.91 -2.07 -1.99
max 13.11 2.56 2.57 2.38 2.15
RMS 2.96 1.04 0.52 0.37 0.46
mean 0.29 -0.68 -0.05 -0.03 0.21
stdev 2.94 0.79 0.52 0.37 0.41
TABLE I P ERFORMANCE COMPARISON FOR A RANDOM
MOVE PROFILE
Fig. 13.
Fig. 12.
Tracking error for random moves: inverse dynamics vs. FIR filter
Desired output trajectory
VII. C ONCLUSION This paper presents the experimental results of developing high performance feedforward control of a positioning system. Frequency domain subspace identification method with input delays is applied to the experimental frequency response. The identified model is first used in the inverse dynamics control law. Further performance enhancement is obtained by using iterative refinement using the identified model to provide a descent direction. The iterative refinement results are next used to train an FIR filter. Finally, the FIR filter coefficients are updated by using an adaptive LMS scheme. Tracking error for a random motion profile shows significant tracking error reduction by using the adaptive LMS filter. ACKNOWLEDGMENT This material is based upon work supported in part by the National Science Foundation under Grant CMS-0301827, and in part by the Center for Automation Technologies under a block grant from the New York State Office of Science, Technology, and Academic Research.
Fig. 14. Tracking error for random moves: comparison between FIR, x, xv, and xva filters
[5] J.X. Xu and Y. Tan. On the P-type and Newton-type ILC schemes for dynamic systems with non-affine-in-input factors. Automatica, 38(7):1237–1242, July 2002. [6] W. Cheng and J.T. Wen. A class of learning controllers with application to tracking control of a flexible beam. In Proc. 1993 IEEE Conference on Robotics and Automation, pages 411–416, Atlanta, GA, May 1993. [7] D.M. Gorinevsky. Direct learning of feedforward control for manipulator path tracking. In Proceedings of the 1992 IEEE International Symposium on Intelligent Control, pages 42–47, Glasgow, UK, August 1992. [8] K.E. Avrachenkov. Iterative learning control based on quasi-newton methods. In Conference on Decision and Control, pages 170–174, Tampa, FL, December 1998. [9] S. Haykin. Adaptive Filter Theory. Prentice-Hall, 3 edition, 1996. [10] G.C. Goodwin and K.S. Sin. Adaptive Filtering Prediction and Control. Prentice-Hall, 1984.
A PPENDIX I I TERATIVE R EFINEMENT A LGORITHM
R EFERENCES [1] T. McKelvey, H. Akcay, and L. Ljung. Subspace-based multivariable system identification from frequency response data. IEEE Trans. Automat. Contr., 41(7):960–979, July 1996. [2] P. Van Overschee and B. De Moor. Coninuous-time frequency domain subspace system identification. signal Processing, Special Issue on Subspace Methods, Part: System Identification, 52:179–194, 1996. [3] S. Arimoto. Mathematical theory of learning with applications to robot control. In K.S. Narendra, editor, Adaptive and Learning Systems: Theory and Applications. Plenum Press, New York, 1985. [4] J.X. Xu and Q. Ji. New ILC algorithms with improved convergence for a class of non-affine functions. In Conference on Decision and Control, pages 660–665, Tampa, FL, December 1998.
Consider a nonlinear dynamical system H : Lm 2 [0, T ] → that is initially at rest:
Lp2 [0, T ]
Define V =
y = H(u).
(6)
1 2 ky − ydes kLp . 2 2
(7)
Note that positive time dependent weighting may be used in the above norm with obvious modifications in the subsequent
update law. The variation of V is δV = hy − ydes , δyiLp
(8)
2
where h·, ·iL2 denotes the L2 inner product. From (6), δy is δy = ∇u Hδu
(9)
where ∇u H is the Fr´echet derivative of H. Suppose H is given by the following state space representation: x˙ = f (x, u), x(0) = x0 ;
y = h(x, u).
(10)
Eq. (19) needs to be integrated backward in time. Through a change of variable, z1 (t) = z(T − t), we can instead use the forward solution. The adjoint system in this case is then (with z1 (0) = 0 and δv1 (t) := δv(T − t)): z˙1 (t) = AT (T − t)z1 (t) + C T (T − t)δv1 (t) w(T − t) = B T (T − t)z1 (t) + DT (T − t)δv1 (t).(20) This can be implemented by reversing δv backward in time, filtering it through (20), and then reversing the result in time again.
Then ∇u H is the time varying linearized system about x and u (x is the state trajectory corresponding to u):
A PPENDIX II A DAPTIVE FIR A LGORITHM
∂f (x, u) ∂f (x, u) δx + δu, δx(0) = 0 ∂x ∂u ∂h(x, u) ∂h(x, u) δy = δx + δu, (11) ∂x ∂u where ∇u H is the operator that maps δu to δy. The gradient descent algorithm choose δu based on ∇u H ∗ :
where ydes is the desired output, w is the constant weighting vector, and u∗ is a feedforward control that depends on ydes only. Then the output is
δu = −α∇u H ∗ (y − ydes ).
y = (Hy des )T w + Hu∗ .
δ x˙ =
δy = −α(∇u H)(∇u H ∗ (y − ydes ))
2
δV = −α k∇u H ∗ (y − ydes )kLp . 2
uk+1 = uk − αk ∇u H (uk )(yk − ydes ),
ydes = (Hy des )T w∗ + Hu∗ .
(14)
Our goal is to choose an update law for w so that y → ydes even though w∗ is unknown. To simplify notation, let a = Hy des , and write (23) in discrete time as
(15)
where αk is found through a line search to ensure Vk+1 (V with uk+1 ) is strictly smaller than Vk . Now we consider the computation of ∇u H ∗ . Write ∇u H as a linear time varying system (same as (11)): δ x˙ = A(t)δx + B(t)δu, δy = C(t)δx + D(t)δu.
(16)
Let Φ(t, τ ) be the state transition matrix from τ to t. The output can now be written as Z t δy(t) = C(t) Φ(t, τ )B(τ )δu(τ ) dτ + D(t)δu(t). 0 ∗
To find ∇u H , we use the defining relationship hδv, ∇u Hδui = h∇u H ∗ δv, δui for all δv and δu in L2 [0, T ]. We then obtain Z T ∇u H ∗ δv(τ ) = B T (τ )ΦT (t, τ )C T (t)δv(t) dt+DT (τ )δv(τ ). τ
(17)
z(τ ) =
T
ΦT (t, τ )C T (t)δv(t) dt.
(22)
(13)
Provided that ∇u H is onto, V will converge to zero. For a discrete implementation, (12) is replaced by ∗
(21)
Assume that an optimal w∗ exists such that
which implies
Z
u = y Tdes w + u∗
(12)
If this can be implemented exactly, then from (9)
Define
Consider a SISO LTI system, y = Hu. Choose
(18)
τ
The state space representation of ∇u H ∗ is then given by the mapping from δv to w as follows: z(t) ˙ = −AT (t)z(t) − C T (t)δv(t), z(T ) = 0 w(t) = B T (t)z(t) + DT (t)δv(t). (19)
ydes (t) − Hu∗ (t) = aT (t − 1)w∗ ,
(23)
(24)
where the left hand side and a(t − 1) are known. The problem then becomes a standard least square parameter estimation problem [10, Sec.3.3]. Various algorithms, ranging from projection to least square with covariance resetting may be used. To reduce the real-time computation load, we use the gradient update law: w(t) = w(t − 1) − a(t − 1)(y(t) − y des (t)).
(25)
To analyze the convergence of this scheme, first compute the evolution of the parameter estimation error (δw := w − w∗ ): 2
2
kδw(t)k − kδw(t − 1)k = (−2 + aT (t − 1)a(t − 1))e2 (t). (26) where e(t) := ydes (t) − y(t). To ensure that the parameter error does not increase, we need choose sufficiently small (a(t) is assumed uniformly bounded). It follows that ∞ X
(aT (t − 1)a(t − 1) − 2)e2 (t) < ∞
(27)
t=1
which implies (aT (t − 1)a(t − 1) − 2)e2 (t) → 0 as t → ∞. If aT (t − 1)a(t − 1) does not converge to 2, then e(t) converges to zero as desired. Other type of parameter estimation algorithms with more relaxed convergence condition may be used (e.g., projection, least square with covariance reset, etc.). However, the realtime computation load would be more demanding, possibly preventing their implementation.
