Use of filtered basis splines to compensate servo-induced motion errors

This application is a U.S. National Phase Application under 35 U.S.C. 371 of International Application No. PCT/US2016/044491 filed on Jul. 28, 2016 and published in English as WO 2017/019869 A1 on Feb. 2, 2017. This application claims the benefit of U.S. Provisional Application No. 62/197,682, filed on Jul. 28, 2015. The entire disclosures of the above applications are incorporated herein by reference.

This invention was made with government support under CMMI1350202 awarded by the National Science Foundation. The Government has certain rights in the invention.

The present disclosure relates to motion errors and, more particularly, relates to the use of filtered basis functions derived from a uniform or non-uniform rational B-spline curve to compensate for servo-induced motion errors.

This section provides background information related to the present disclosure which is not necessarily prior art. This section provides a general summary of the disclosure, and is not a comprehensive disclosure of its full scope or all of its features.

The present teachings provide an approach for minimizing tracking errors in dynamic systems (including those with non-minimum phase (NMP)) zeros by using filtered basis splines (B-splines). The output of the tracking controller is represented as a linear combination of B-splines having unknown coefficients. The basis functions are forward filtered using the dynamics of the system and their coefficients selected to minimize the errors in tracking a given trajectory. Among many other benefits, the use of B-splines allows: (i) constraints to be imposed on the input and outputs of the system as linear functions of the spline coefficients; (ii) non-zero initial conditions to be imposed; (iii) tracking of higher order derivatives; and (iv) the use of limited look ahead capabilities in the filtering process.

Further areas of applicability will become apparent from the description provided herein. The description and specific examples in this summary are intended for purposes of illustration only and are not intended to limit the scope of the present disclosure.

The drawings described herein are for illustrative purposes only of selected embodiments and not all possible implementations, and are not intended to limit the scope of the present disclosure.

FIG. 1 shows the dynamics of a CNC feed axis including unwanted vibration modes;

FIG. 2 is a flow chart of the proposed approach

FIG. 3 shows the experimental setup i.e., a biaxial stage with a flexible fixture;

FIG. 4 shows the measured and modeled frequency response functions for the X and Y axis of the experimental setup;

FIG. 5 shows the desired path used for experiments;

FIG. 6 shows the X and Y command signals for uncompensated and proposed method;

FIG. 7 shows the Fast Fourier Transform (FFT) of the measured acceleration signals of the fixture for the uncompensated and proposed trajectories;

FIG. 8 shows the X and Y axis position tracking errors for the uncompensated and proposed methods; and

FIG. 9 shows the contour errors for uncompensated and proposed trajectory;

FIG. 10 shows the desired path used for simulations to analyze the limited look ahead capability

FIG. 11 shows the tracking errors with and without limited look ahead for the proposed method

Corresponding reference numerals indicate corresponding parts throughout the several views of the drawings.

Example embodiments will now be described more fully with reference to the accompanying drawings. Example embodiments are provided so that this disclosure will be thorough, and will fully convey the scope to those who are skilled in the art. Numerous specific details are set forth such as examples of specific components, devices, and methods, to provide a thorough understanding of embodiments of the present disclosure. It will be apparent to those skilled in the art that specific details need not be employed, that example embodiments may be embodied in many different forms and that neither should be construed to limit the scope of the disclosure. In some example embodiments, well-known processes, well-known device structures, and well-known technologies are not described in detail.

The terminology used herein is for the purpose of describing particular example embodiments only and is not intended to be limiting. As used herein, the singular forms “a,” “an,” and “the” may be intended to include the plural forms as well, unless the context clearly indicates otherwise. The terms “comprises,” “comprising,” “including,” and “having,” are inclusive and therefore specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof. The method steps, processes, and operations described herein are not to be construed as necessarily requiring their performance in the particular order discussed or illustrated, unless specifically identified as an order of performance. It is also to be understood that additional or alternative steps may be employed.

Tracking of Dynamic Systems Including Input or Output Constraints

The goal of tracking control is to force the output of a system to follow a desired trajectory as closely as possible. Consider a linear time-invariant (LTI) dynamic system with transfer function G, as shown in FIG. 1 . The goal of tracking control is to ensure that the actual trajectory, χ, follows the desired trajectory, χ d , as closely as possible, while satisfying constraints imposed on χ and/or χ d (or their time derivatives).

Furthermore, consider a specific example where G represents the dynamics of the feed axis of a CNC machine with unwanted vibration modes, given by

where M represents the number of modes, while ∝ i and β i are the residue constants, ζ i is the damping ratio and ω n,i is the natural frequency of the i th mode in rad/s. Using tracking control, it is desired to generate χ dm *, a modified and optimized version of desired position trajectory χ d , such that, when applied to the feed axis in place of χ d , minimizes the tracking errors caused by G while satisfying the kinematic limits of the machine axis.

FIG. 2 shows a block diagram of the approach according to the present teachings for generating χ dm * using filtered B-splines. Let χ d =[χ d (0)χ d (1) . . . χ d (E)] T represent E+1 discrete time steps of the x-component of a toolpath interpolated by a CNC (i.e., the desired trajectory). We assume that the CNC has look ahead capabilities such that the E+1 time steps of χ d are known in advance. Furthermore, we assume that the modified but un-optimized motion command, χ dm =[χ dm (0)χ dm (1) . . . χ dm (E)] T , is parameterized in the CNC using B-splines such that

For a B-spline of degree m, having n+1≤E+1 control points, p 0 , p 1 , . . . , p n , and knot vector, [g 0 g 1 . . . g m+n+1 ] T , N j,m (ξ) are real-valued basis functions given by

where j=0,1, . . . , n and ξ∈[0,1] is the spline parameter, representing normalized time, which is discretized in Eq. (2) into E+1 uniformly spaced points, ξ 0 ,ξ 1 , . . . , ξ E . The knot vector can be selected to be uniform, such that

Notice that if χ dm given in Eq. (2) is applied as the motion command to the CNC feed axis of FIG. 1 , in place of χ d , the control point vector, p, is not affected by G, because p is not a function of time; however, the basis functions matrix, N, which is a function of time, is filtered by G. Therefore, the resulting position χ=[χ(0)χ(1) . . . χ(E)] T can be approximated by

χ= ÑP   (5)

e=χ d −χ=χ d −Ñ p   (6)

One useful mathematical property of B-splines is that the derivative of a B-spline curve is also a B-spline curve that can be expressed linearly in terms of the control points of the original curve. Hence, the r th derivative of χ dm can be written as

where N r is the basis function matrix for the r th derivative of position; its elements can be calculated readily from those of N. Therefore, to achieve the stated objective of the present teachings, χ dm *=Np* can be generated using optimal control points, p*, calculated to minimize the objective function given by

where L r represents the kinematic limit associated with the r th derivative of χ dm (e.g., its velocity, acceleration and jerk limits), while ∥·∥ ∞ denotes the infinity norm (i.e., the maximum absolute value of the elements) of a vector.

Note that constraints can be similarly placed on the output, χ, and its derivatives by swapping N r in Eq.(8) with Ñ r (i.e., the filtered version of N r ).

Notice that, if the kinematic constraints are not considered, the solution to Eq. (8) reduces to a traditional least-squares fitting problem where

p *=(( Ñ T Ñ ) −1 Ñ T )χ d   (9)

However, with the constraints included, Eq. (8) becomes a least squares problem with inequality constraints, which has a guaranteed solution as long as the constraints are feasible. The solution is unique if the matrix Ñ is full rank (i.e., if the filtered basis functions are linearly independent). The present method achieves trajectory parameterization, constraint enforcement and optimization simultaneously using industry-standard B-splines and basic least squares techniques, making it easy to use. Notice also that both χ d and χ dm * have the same length (i.e., E+1), meaning that the proposed method does not sacrifice motion speed. Though the method has been discussed for a single axis (i.e., χ), it is readily applicable to multiple axes of a CNC machine by treating each axis separately; therefore, it is very systematic and scalable.

Handling of Non-Zero Initial Conditions

In tracking control, the trajectory to be tracked could have non-zero initial conditions associated with the trajectory itself, as well as with its higher-order derivatives. This is realized through the linear relationship between the control points and the higher order derivatives in Eq. (7). For example, the first derivative of χ dm at initial time ξ 0 (i.e., when ξ=ξ 0 ) is expressed as

x dm ′(ξ 0 )= m ( n−m+ 1)( p 1 −p 0 )  (10)

Higher order derivatives of Eq. (7) can be obtained in a similar manner. Note that the value of the first derivative at ξ=ξ 0 depends on the first two control points. The value of second derivative will depend on the first three control points and so on. Consequently, an initial condition on the r th derivative of a given trajectory constraints the first r+1 control points associated with the trajectory (where r<m). Hence, the initial conditions of x d (which are the same as those of x dm ) can be expressed as

Qp=z   (11)

where z is a vector of initial conditions and Q is a matrix comprising the coefficients obtained from Eqs. (7) and (10), and similar equations derived for higher order derivatives. Accordingly, Eqs. (9) and (11) can be combined to obtain

where is a vector of Lagrange multipliers. The solution of Eq. (12) (in a manner similar to Eq. (9)) results in minimum tracking error while satisfying the conditions given by Eq. (11).

Enhanced Tracking of Higher Order Derivatives

In certain tracking applications, e.g., position tracking in robots and CNC machines, it may be of interest to also maintain close tracking of higher order derivatives of the desired trajectory, e.g., velocity and acceleration. The B-spline framework is very useful in this regard.

Filtering N 1 ,N 2 , . . . N r using G gives, Ñ 1 , Ñ 2 , . . . , Ñ r , respectively. Consequently, accurate tracking of χ d and its first r derivatives can be achieved by solving the composite equation given by

where α 1 , α 2 , . . . α r are user-defined weights attached to each derivative of χ d to indicate its importance relative to position tracking. The solution of Eq.(14) can be obtained using least squares in a manner similar to Eq.(9).

Limited Look Ahead

In explaining the proposed method so far, it has been assumed that all E+1 points of the trajectory, x dm , is known, and that the optimal coefficients p* have been solved using all E+1 points of x dm . Basically, complete (total) look ahead capability has been assumed.

However, in practice x dm may not be known in its entirety; instead it may be available in smaller batches (subsets) which sequentially combine to form x dm in its entirety. Even when all of x dm is available, it may be computationally prohibitive (e.g., for real time implementation) to use all of x dm at once for generating p*.

In such cases, the application of the proposed method can be performed sequentially (recursively) using smaller subsets of x dm , thus provided limited look ahead as opposed to complete look ahead capability.

One way of doing this is to generate matrices N and Ñ as described above and then partition them into smaller submatrices correctly sized to match any subset of x dm that is available. A corresponding subset of the optimal coefficient vector p* is then solved for using the smaller submatrices. This process is repeated sequentially for subsequent subsets of x dm (as they become available).

Continuity between solutions from each batch of x dm can be assured by creating some overlap between the submatrices and coefficients of successive subsets of x dm .

Experimental Examples

Experimental results are presented in this section to compare the proposed technique to the uncompensated case. The biaxial (X-Y) linear motor driven stage (Aerotech ALS 25010) shown in FIG. 3 is used for the experiments. The stage is controlled using a traditional P/PI feedback controller, augmented with velocity and acceleration feedforward. The controller is implemented on a dSPACE 1103 real-time control board with 10 kHz sampling frequency. The velocity and acceleration limits for both of its axes are L 1 =100 mm/s and L 2 =8000 mm/s 2 , respectively. As shown in the FIG. 3 , the stage is equipped with a fixture consisting of a block mounted on a rod. The block is assumed to represent an apparatus (e.g., a tool, workpiece or measurement device) whose x and y positions are expected to track their respective desired trajectories, x d and y d , accurately, in spite of inherent structural flexibilities. The acceleration of the fixture is measured using two unidirectional accelerometers (PCB Piezotronics 393B05) shown in FIG. 3 .

FIG. 4 shows the frequency response function (FRF) of the dynamics G of each axis of the stage, generated by applying swept sine acceleration commands to the stage and measuring the corresponding accelerations of the fixture (i.e., and) using the accelerometers. The dynamics of each axis exhibits two dominant modes, which are very obvious in FIG. 4 .

For each axis, the lower mode is a mode from the stage while the higher mode is from the flexible fixture (whose natural frequencies are different in the x and y directions because the mounting rod has a rectangular cross section). However, not obvious in FIG. 4 are two much less dominant modes in each FRF. In other words, there are four unwanted modes (M=4) in each axis' dynamics, whose parameters are identified using the rational fraction polynomial method, and are summarized in Tab. 1. FIG. 4 shows an excellent fit between the modeled and actual FRFs for the X-axis. The model for the Y-axis is not as accurate as that of the X-axis because the two dominant modes of the Y axis are relatively more damped and closer in resonance frequency, making it harder to curve fit them. However, the fitting errors of the Y axis provide an opportunity to test the performance of the proposed method in the presence of some modeling errors, which are likely to occur in practice.

FIG. 5 shows the butterfly curve used as the desired path for the experiments. FIG. 6 shows the position trajectories, x d and y d , (i.e., the uncompensated motion commands) of the desired path, as well as its velocity and acceleration profiles, as functions of time. Notice that the duration of the desired trajectories is 1 s (i.e., E=104, based on 10 kHz sampling frequency). The uncompensated trajectories are compared with the corresponding modified trajectories obtained from the proposed method. The proposed method is executed with B-splines of order m=5, having n+1=51 control points. Notice from FIG. 6 that the modified trajectories generated by the proposed method automatically satisfy the imposed kinematic constraints and do not have any delays (i.e., they last exactly as long as the desired trajectories).

The position trajectories from the uncompensated and the proposed methods are applied to the X and Y axes of the stage as motion commands. FIG. 7 shows the Fast Fourier Transform (FFT) of the acceleration signals (measured using the accelerometers). The X-axis acceleration energy spectral densities are 6.38×106 mm 2 /s 4 and 1.92×106 mm2/s4, while the Y-axis acceleration energy spectral densities are 7.32×106 mm2/s4 and 1.81×106 mm2/s4 for the uncompensated and the proposed methods, respectively. This means that the proposed method generates less vibration than the uncompensated. To demonstrate this fact, the tracking errors, xd-x and yd-y, for all three cases are plotted in FIG. 8 , based on position signals (x and y) derived from the measured acceleration signals using an observer. The X-axis RMS tracking errors are 0.3014 mm and 0.0852 mm, while the Y-axis RMS tracking errors are 1.0165 mm and 0.1198 mm for the uncompensated and the proposed methods, respectively. FIG. 9 compares the contour errors of the three methods. Their RMS values are 0.5508 mm and 0.0880 mm for the uncompensated and the proposed methods, respectively, indicating the ability of the proposed method to indirectly reduce contour errors via minimizing tracking errors.

Examples on Limited Look Ahead Simulation

FIG. 11 shows the positioning performance tracking a butterfly shape trajectory (shown in FIG. 10 ) with and without limited look ahead, and the statistics are provided in Table 2. The tracking performances are identical however 96% of the computational time is saved with the limited look ahead implementation.