- Peer review
- Editorial board
- Manuscript tracking
Correspondence: Hao Li (email@example.com)
In high-speed and high-precision machinery, trajectories with high-frequency harmonic content are one of the main sources of reduction of operational precision. Trajectories with high-frequency harmonic content generally demand even higher-harmonic actuating forces/torques due to the nonlinear dynamics of such systems, which may excite natural modes of vibration of the system and/or be beyond the dynamic response limitation of the actuation devices. In this paper, a global interpolation algorithm that uses the trajectory pattern method (TPM) for synthesizing low-harmonic trajectories is presented. The trajectory synthesis with the TPM is performed with a prescribed fundamental frequency and continuous jounce boundary condition, which would minimize the number of high-harmonic components in the required actuation forces/torques and avoid excitation of the system modes of vibration. The minimal curvature variation energy method, Lagrange multiplier method, and contour error control are used to obtain smooth kinematic profiles and satisfy the trajectory accuracy requirements. As an example, trajectory patterns that consist of a fundamental frequency sinusoidal time function and its first three harmonics are used to synthesize the desired trajectories for a selected dynamic system. The synthesized trajectories are shown to cause minimal system vibration during its operation. A comparison with a commonly used trajectory synthesis method clearly shows the superiority of the developed TPM-based approach in reducing vibration and demand on the actuator dynamic response, thereby allowing the system to operate at higher speeds and precision.
For the trajectory synthesis (Chu et al., 2020; Dai et al., 2020; Van Loock et al., 2015) of high-speed machinery, such as CNC (computerized numerical control) machinery, robot manipulators and other computer-controlled machines, polynomial-based curves (Analooee et al., 2020; Liang and Su, 2019; Shen et al., 2020) are still the most widely used trajectories. It has the following advantages: (1) numerically stable algorithms; (2) it is more natural for designing and representing shape in a computer. The coefficients of functions possess considerable geometric significance. This translates into intuitive design methods. However, it is important to note that when expressed in Fourier series, all synthesized trajectories – such as those based on polynomial (Analooee et al., 2020; Shen et al., 2020), Akima (Bica, 2014; Wang et al., 2014), or B spline (Du et al., 2018; Simba et al., 2016) curves – contain significant high-harmonic content. In addition, the nonlinearity of the machine system dynamics would require actuating forces/torques with even higher-harmonic content to follow the prescribed trajectories. The lack of control over the harmonic content of the required actuating forces/torques and their frequencies in currently available and published methods for trajectory synthesis leads to excitation of natural modes of vibration of high-speed machinery.
In order to avoid exciting the natural modes of the mechanical structure or servo control system, many researchers have focused on flexible feed-rate planning methods for high-speed and high-precision machinery. Erkorkmaz and Altintas (2001) provided a continuous feed motion by constructing trapezoidal acceleration profiles along the quintic spline trajectory. Lee and Choi (2015) obtained smoother kinematic profiles with a discretized sinusoidal jerk function. To improve the flexibility, Wang (2015) proposed a trigonometric velocity scheduling algorithm, which could generate smooth velocity, acceleration, and jerk curves of parametric interpolation. Huang and Zhu (2016) represent jerk profiles by the sine series for feed-rate planning of parametric interpolator. The coefficients of the sine series are determined by selecting the geometric sequence. To reduce the time spent by the smoothing and planning process, Li et al. (2019) developed a real-time and look-ahead interpolation algorithm with axial jerk-smooth profiles. In one step, the algorithm finishes the trajectory synthesis and velocity planning by using the trigonometric velocity planning method. From these methods mentioned above, it can be seen that the vibration reduction effect is not obvious any more along with the increase of flexibility further, since the frequency distribution of synthesized kinematic profiles is not considered. The sine series method constructed jerk profiles with sinusoidal functions with an appropriate fundamental frequency and all harmonics, which contained considerable high-harmonic content. Although other methods did motion planning by using one sinusoidal function with a specified fundamental frequency, their velocity or/and displacement functions exist as polynomial items, which also results in significant high-harmonic components. Biagiotti and Melchiorri (2012) and Tajima et al. (2018) utilized finite impulse response filters to remove high-harmonic content of synthesized trajectories. However, filtering introduces unavoidable delay and induces large contouring errors in multi-axis motion that must be compensated since the filtered components are needed to achieve the prescribed motion.
Rastegar and Fardanesh (1990) introduced the concept of trajectory patterns that are constructed with sinusoidal functions with an appropriate fundamental frequency. Rastegar and Feng (2011) applied it to synthesis of trajectories with low-harmonic content such that the harmonic content of the required actuation forces/torques is controlled and does not contain significant high-harmonic content. The synthesized trajectory is a unique combination of a fundamental frequency harmonic and its second harmonic, which results in the minimum number of harmonics in the required actuation forces/torques and minimizes excitation of the system modes of vibrations. But this method is only used for point-to-point motions with zero end point velocity, acceleration, and jerk. To improve motion efficiency and precision, this paper, therefore, extends the trajectory pattern method (TPM) to construct trajectories with continuous velocity, acceleration, jerk, and jounce at mid-points. The fundamental frequency of the synthesized trajectories can be selected such that together with all its harmonics that appear in the trajectory and the required actuating forces/torques would not excite the natural modes of vibration of the system. The actuating forces/torques can then be readily optimized and used in a model-based control scheme to achieve high-speed and high-precision operation of the system.
In this paper, a path planning algorithm is presented that uses harmonic-based trajectory patterns for high-speed machinery. In Sect. 2, a global interpolation method using the TPM with a specified fundamental frequency, motion periods, and chord error is introduced. Minimal curvature variation energy and Lagrange multiplier method are used to synthesize the trajectories. The synthesized trajectories and axis kinematic profiles are smooth and only contain the fundamental frequency and its first three harmonics, which are designed to minimize vibration of the high-speed machinery. In Sect. 3, the effectiveness of the proposed algorithm is illustrated through computer simulations studies. Conclusions of the present study are presented in Sect. 4.
In this section, a trajectory pattern method (TPM)-based global interpolation algorithm, which achieves smooth axis motion profiles with chord error constraint, is presented. In these profiles, trajectory patterns that are synthesized with a fundamental frequency and its first three harmonics are used to achieve minimum trajectory harmonic content (Rastegar and Feng, 2011).
The class of trajectory patterns used in the present study is harmonic based and is formed by a fundamental sinusoidal function and (n−1) number of its harmonics:
where ai and bi are trajectory parameters; t is the time, , and T is the motion period; f is the selected fundamental frequency of the trajectory pattern.
In the present study, the selected trajectory pattern used for trajectory synthesis consists of a fundamental frequency sinusoidal function and its first three harmonics. The displacement function is then expressed as
From Eq. (2), the velocity, acceleration, jerk, and jounce functions are then derived as
In the following sections, the present low-harmonic trajectory synthesis is presented by its application to a point-to-point trajectory synthesis. The method is however general and may be applied to other motion scenarios, for example for motions starting from rest and ending to start a constant velocity motion or from a constant velocity motion to a stop, a scenario that is very often encountered in high-speed and high-precision machine operations.
Figure 1 shows the flow diagram of the developed global interpolation algorithm. The trajectory synthesis includes five modules: an initialization module; a boundary conditions definition module; a jerk energy function definition module; a trajectory parameters calculation module; and a contour error control module. A Lagrange multiplier method is used to calculate the trajectory parameters with the indicated objective function and equality constraints. The above modules are described in more detail in the following sections.
Figure 2 presents point-to-point motion from the point P1 to the point Pn, which is divided into n−1 segments. In this example, the path segments are considered to line segments. The path of the point-to-point motion is therefore given by n number of data points, P1, P2 …Pn−1, Pn, and (n−1) number of segment trajectories s1(t1), s2(t2) …, are to be synthesized. Here, n is a positive integer and n≥2, and t1, t2 …tn−2, tn−1 are trajectory segment times.
The user is to provide the required contour error ε, the motion period for each trajectory segment, T1, T2 …Tn−2, Tn−1, with , i=1, …, n−1, and the fundamental frequency f of the trajectory pattern for motion from P1 to Pn, which must be selected such that together with its first three harmonics, the natural modes of vibration of the system are not excited.
To simplify the construction process and obtain continuous smooth axis motion profiles, the cycle time T is set equal to a quarter of the fundamental frequency period, i.e., .
For the trajectory segments, motion periods and corresponding fundamental frequencies f1, f2 …fn−2, fn−1 are adjusted to reduce the high-frequency harmonic content of the trajectory. In addition, for each trajectory segment, such as the segment PiPi+1, if Ti>T, we let ; otherwise let Ti=T and fi=f.
In Fig. 2, the trajectory segment curves are drawn in red and are shown to pass through the path segment ends P1, P2 …Pn−1, Pn with smooth kinematic profiles.
To achieve continuous smooth axis motion profiles, boundary conditions are necessary. In addition, since the motion is point to point, the velocities, accelerations, jerks, and jounces at the start point P1 and end point Pn−1 are set to zero. Then using Eqs. (2)–(6), the initial conditions of the first curve segment s1(t1) at t1=0 are obtained as
Similarly, the terminal conditions of the last curve segment at are expressed as
In Fig. 2, for mid-points, P2 …Pn−1, adjacent curve segments must have the same velocity, acceleration, jerk, and jounce at the nodes, that is
Then substituting Eqs. (2)–(6) into Eq. (17), we can get
where i=1, …, n−2.
The internal energy of curve segments can affect its shape and smoothness. The well-known examples of internal energy functions (Zhang et al., 2001) are the stretch energy, strain energy, and jerk energy. The stretch energy measures the length of a curve. The strain energy measures how much a curve is bent. And the jerk energy (Meier and Nowacki, 1987) measures the curvature variation of a curve. In this paper, we use the jerk energy to obtain continuous smooth axis kinematic profiles with low-harmonic content. The jerk energy is defined as
As can be seen in Fig. 2, we have the following (9n−9) number of trajectory parameters: a1,0, a1,1, a1,2, a1,3, a1,4, b1,1, b1,2, b1,3, b1,4 …, , , , , , , , and . In addition, Eqs. (7)–(17) define (6n−2) number of boundary conditions. To calculate the trajectory parameters, the Lagrange multiplier method is used in the present study. In this formulation, the jerk energy (Eq. 24) is used as the objective function to be minimized. The optimization problem is thereby defined as
where λi, ξi, , , , , , and are Lagrange multipliers. Now let
Then, the gradient as the first-order partial derivative of Eq. (25) with respect to X is set to zero, i.e.,
The trajectory parameters, a1,0, a1,1, a1,2, a1,3, a1,4, b1,1, b1,2, b1,3, b1,4 …, , , , , , , , , can then be obtained by solving Eq. (26). By substituting the calculated trajectory parameters into Eqs. (2)–(6), a continuous smooth trajectory with low-harmonic content is then obtained.
After segment trajectories are synthesized using the above steps, the following process is used in this study for bounding the contour error. In the present example, the path segments are presented as lines as shown in Fig. 2, and since the segment trajectories are synthesized with a fundamental frequency time function and its three harmonics, by checking contour error at a limited number of intervals in each segment the maximum contour error is accurately estimated. In addition, since the harmonic structure of each segment trajectory and the harmonic amplitudes and phases are known, a more exact analytic formulation is possible and will be presented in future publications.
In general, the path segments may not be lines and the maximum contour error must then be determined, for example, numerically by calculating the error at appropriate time or distance intervals or by various search or approximation methods. In the present example, the contour error is calculated at three equal intervals in each segment of the path and used for maximum contour error determination. Then if the contour error is larger than the prescribed threshold ε, the period of the fundamental frequency of the trajectory of that segment is then reduced as described in the following steps.
Step 1. For the ith trajectory segment PiPi+1, as shown in Fig. 3, let , and the maximum contour error .
Step 2. Set the times , , and . Calculate the corresponding position points si(ti,1), si(ti,2), and si(ti,3) for each indicated time from Eq. (2) with the calculated trajectory parameters. Then compute the distances d1, d2, and d3 to the line segment PiPi+1 (Fig. 3).
Step 3. Set , d2, d3). If , the maximum contour error of the ith segment curve si(ti) exceeds the prescribed threshold and the motion period Ti is reduced in proportion to the amount that the contour error threshold has been exceeded as
And let ; otherwise, contour error threshold has been met in the ith segment of the trajectory.
Step 4. For all synthesized trajectory segments, go through the above steps to control their contour errors. If contour error threshold has been exceeded in any segment, go to the Sect. 2.2.4 and recalculate trajectory parameters with the updated segment trajectory periods and fundamental frequencies.
In this section, simulations demonstrating the effectiveness and performance of the proposed low-harmonic trajectory synthesis algorithm are presented. The results are also compared with those obtained with one of the commonly used polynomial algorithms for synthesizing robot trajectories. Visual Studio 2013 was used to implement both algorithms, and results are plotted in MATLAB 2015.
To evaluate the proposed algorithm, a quintic polynomial algorithm (Perumalsamy et al., 2019), which interpolates data points using quintic polynomials with continuous displacement, velocity, acceleration, and jerk values, is used as the comparison algorithm.
As shown in Fig. 4, a letter Z-shaped pattern is selected as the path, which contains five data points, P1(0,3), P2(3,3), P3(1.5,1.5), P4(0,0), P5(3,0), and consists of four segments, P1P2, P2P3, P3P4, P4P5. The following variables were initialized: a fundamental frequency of f=20 Hz, the motion periods of segments, T1=0.68 s, T2=0.32 s, T3=0.32 s, T4=0.68 s, and a contour error limit of ε=0.25 mm. The synthesized trajectories synthesized by applying the two algorithms are shown in Fig. 4, and we can see that they are all smooth.
Figure 5 shows that both algorithms can obtain smooth axis displacement, velocity, and acceleration profiles, and the jerk of the proposed algorithm is also smooth because of satisfying continuous jounce conditions, while that of the comparison algorithm is not smooth.
To evaluate the proposed algorithm further, frequency spectrum distribution of the synthesized kinematic profiles is studied using Fourier integration. In the plots of Fig. 6, the axial displacements, velocities, accelerations, and jerks of the trajectory synthesized using the quintic polynomial algorithm are shown to contain a considerable number of high-harmonic components. As a result, the related actuators are required to produce high-frequency components, which for high-speed motions would in general be beyond their dynamic response and could excite natural modes of vibration of the system. As expected, the velocity, acceleration, and jerk profiles of the synthesized trajectory with the proposed algorithm are seen in Fig. 6 to only contain a finite number of low-harmonic components.
To explain how trajectories with high-harmonic components affect vibration of machinery, an undamped, 1-degree-of-freedom linear spring–mass system of Fig. 7 is used to model the vibration behavior of a machine. Let the spring have a spring rate and free length of k and l, respectively, and the mass of the system be m. The displacements of the mass and the system from their initial position are u and d, respectively. From the Newton's second law, the system equation of motion is seen to be
Let , and set the natural frequency of the spring–mass system to be
Then, Eq. (28) can be rewritten as
Now, take the x-axis motion profiles calculated using the aforementioned two algorithms as input to the dynamic system described by Eq. (30). The input motion is the x-axis acceleration profile, i.e., . Substituting into Eq. (30) and solving the resulting second-order ordinary differential equation, the vibrational motion of the mass m of the spring–mass system of Fig. 7 is determined.
As can be seen in Fig. 8, the displacement z of the proposed algorithm is smoother than that of the quintic polynomial algorithm. In addition, when the x-axis motion ends, the mass m comes to a stop with the proposed algorithm, while the mass of the quintic polynomial algorithm continues to vibrate with its natural mode. It is noted that the trajectories synthesized using the quintic polynomial algorithm contain a significant number of high-harmonic content, which would in general cause excitation of the system modes of vibration irrespective of their frequencies. The synthesized trajectories of the proposed algorithm, however, only contain a finite number of low-frequency harmonics in its velocity, acceleration, and jerk terms and even in their higher-order derivatives, which would minimize the number of high-harmonic components in the required actuation forces/torques. The fundamental frequency of the trajectory pattern can also be selected to avoid natural modes of vibration of the system. As a result, the system should be capable of operating with minimal vibration.
Table 1 shows that the motion times of both algorithms are almost same and that the computing time cost by the proposed algorithm (0.5051 s) is shorter than the quintic polynomial algorithm (0.7763 s). The contour errors with both algorithms are presented in Fig. 9. As can be seen, the maximum contour error of the quintic polynomial algorithm is also significantly larger than that of the proposed algorithm.
A trajectory pattern method (TPM)-based method and related algorithm has been developed for synthesis of low-harmonic trajectories for high-speed and high-precision machinery. The trajectory synthesis with the TPM is performed with a prescribed fundamental frequency and continuous jounce boundary condition, which would minimize the number of high-harmonic components in the required actuation forces/torques and avoid excitation of the system modes of vibration.
The developed TPM-based algorithm is shown to yield smooth axial kinematic profiles with low-harmonic content and thereby reduces demand on the dynamic response of the drive actuators. As a result, and with the capability of selecting the frequency of the fundamental harmonic of the synthesized trajectory and its harmonics present in the actuating forces/torques, a machine can operate at high speed and high precision.
In comparison with the previously published algorithms, the developed algorithm has the following advantages:
The trajectories are constructed with a sinusoidal function with a selected fundamental frequency, which can avoid excitation of the natural modes of vibration of the machinery system.
All segment displacement, velocity, acceleration, and jerk profiles are expressed by a sinusoidal time function consisting of a selected frequency fundamental harmonic and at least its first three harmonics, and the overall synthesized trajectories only contain a finite number of low harmonics, which minimizes the number of high-frequency harmonics of the required actuation forces/torques with the aforementioned operational advantages.
The developed algorithm can realize smooth axial kinematic control action with specified contour error limitation.
In the provided example, since the trajectory synthesis with the TPM is performed with a prescribed fundamental frequency and continuous jounce boundary condition, the trajectory pattern must include at least three harmonics of the fundamental frequency. More harmonics may, however, be used to provide more trajectory parameters to achieve, for example, higher precision, usually at the expense of operating speed.
In the provided example, a point-to-point trajectory was shown to be synthesized. The method and algorithm may readily be used to synthesize general trajectories with or without defined intermediate precision points or with more loosely defined precision points to, for example, minimize contour errors.
All the data used in this paper can be obtained from the corresponding author upon request.
JR proposed the conceptualization of the trajectory pattern method. HL designed the algorithm. BW implemented the comparison and proposed algorithms. HL and BW did the performance analysis. HL wrote the original draft. JR reviewed and edited it.
The authors declare that they have no conflict of interest.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors thank the reviewers for their critical and constructive review of the article.
This research was funded by the Scientific Research Startup Foundation for High-level Talents of Nanjing Institute of Technology (grant no. YKJ2019112) and the Shuang Chuang Doctoral Project of Jiangsu Province.
This paper was edited by Dario Richiedei and reviewed by two anonymous referees.
Analooee, A., Kazemi, R., and Azadi, S.: SCR-Normalize: A novel trajectory planning method based on explicit quintic polynomial curves, Proc. Inst. Mech. Eng. Pt. K-J. Multi-Body Dyn., 234, 650–674, https://doi.org/10.1177/1464419320924196, 2020.
Biagiotti, L. and Melchiorri, C.: FIR filters for online trajectory planning with time- and frequency-domain specifications, Control Eng. Pract., 20, 1385–1399, https://doi.org/10.1016/j.conengprac.2012.08.005, 2012.
Bica, A. M.: Optimizing at the end-points the Akima's interpolation method of smooth curve fitting, Comput. Aided Geom. D., 31, 245–257, https://doi.org/10.1016/j.cagd.2014.03.001, 2014.
Chu, C. H., Chen, H. Y., and Chang, C. H.: Continuity-preserving tool path generation for minimizing machining errors in five-axis CNC flank milling of ruled surfaces, J. Manuf. Syst., 55, 171–178, https://doi.org/10.1016/j.jmsy.2020.03.004, 2020.
Dai, Y., Liu, Z., Qi, Y., Zhang, H., You, B., and Gao, Y.: Spatial cellular robot in orbital truss collision-free path planning, Mech. Sci., 11, 233–250, https://doi.org/10.5194/ms-11-233-2020, 2020.
Du, X., Huang, J., and Zhu, L. M.: A locally optimal transition method with analytical calculation of transition length for computer numerical control machining of short line segments, Proc. Inst. Mech. Eng. Pt. B-J. Eng. Manuf., 232, 2409–2419, https://doi.org/10.1177/0954405417697351, 2018.
Erkorkmaz, K. and Altintas, Y.: High speed CNC system design. Part I: jerk limited trajectory generation and quintic spline interpolation, Int. J. Mach. Tool. Manu., 41, 1323–1345, https://doi.org/10.1016/S0890-6955(01)00002-5, 2001.
Huang, J. and Zhu, L. M.: Feedrate scheduling for interpolation of parametric tool path using the sine series representation of jerk profile, Proc. Inst. Mech. Eng. Pt. B-J. Eng. Manuf., 231, 2359–2371, https://doi.org/10.1177/0954405416629588, 2016.
Lee, A. Y. and Choi, Y. J.: Smooth trajectory planning methods using physical limits, Proc. Inst. Mech. Eng. Pt. C-J. Eng. Mech. Eng., 229, 2127–2143, https://doi.org/10.1177/0954406214553982, 2015.
Li, H., Wu, W. J., Rastegar, J., and Guo, A.: A real-time and look-ahead interpolation algorithm with axial jerk-smooth transition scheme for computer numerical control machining of micro-line segments, Proc. Inst. Mech. Eng. Pt. B-J. Eng. Manuf., 233, 2007–2019, https://doi.org/10.1177/0954405418809768, 2019.
Liang, X. and Su, T.: Quintic Pythagorean-Hodograph curves based trajectory planning for delta robot with a prescribed geometrical constraint, Appl. Sci., 9, 4491, https://doi.org/10.3390/app9214491, 2019.
Meier, H. and Nowacki, H.: Interpolating curves with gradual changes in curvature, Computer Aided Geom. D., 4, 297–305, https://doi.org/10.1016/0167-8396(87)90004-5, 1987.
Perumalsamy, G., Visweswaran, P., Jose, J., Winston, S. J., and Murugan, S.: Quintic Interpolation Joint Trajectory for the Path Planning of a Serial Two-Axis Robotic Arm for PFBR Steam Generator Inspection, in: 3rd International and 18th National Conference on Machines and Mechanisms, Mumbai, India, 13–15 December 2017, 637–648, https://doi.org/10.1007/978-981-10-8597-0_54, 2019.
Rastegar, J. and Fardanesh, B.: A “basis trajectory” approach to the inverse dynamics formulation of robot manipulators, in: Proceedings of the 1990 IEEE International Conference on Systems Engineering, Pittsburgh, USA, 9–11 August 1990, 308–311, https://doi.org/10.1109/ICSYSE.1990.203158, 1990.
Rastegar, J. and Feng, D.: On the minimum harmonic trajectory pattern for point-to-point motion, in: ASME 2011 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Washington, USA, 28–31 August 2011, 359–363, https://doi.org/10.1115/DETC2011-47724, 2011.
Shen, P. Y., Zhang, X. B., Fang, Y. C., and Yuan, M. X.: Real-time acceleration-continuous path-constrained trajectory planning with built-in tradeoff between cruise and time-optimal motions, IEEE T. Autom. Sci. Eng., 17, 1911–1924, https://doi.org/10.1109/TASE.2020.2980423, 2020.
Simba, K. R., Uchiyama, N., and Sano, S.: Real-time smooth trajectory generation for nonholonomic mobile robots using Bezier curves, Robots Comput.-Integr. Manuf., 41, 31–42, https://doi.org/10.1016/j.rcim.2016.02.002, 2016.
Tajima, S., Sencer, B., and Shamoto, E.: Accurate interpolation of machining tool-paths based on FIR filtering, Precis. Eng., 52, 332–344, https://doi.org/10.1016/j.precisioneng.2018.01.016, 2018.
Van Loock, W., Pipeleers, G., and Swevers, J.: B-spline parameterized optimal motion trajectories for robotic systems with guaranteed constraint satisfaction, Mech. Sci., 6, 163–171, https://doi.org/10.5194/ms-6-163-2015, 2015.
Wang, Y. S., Yang, D. S., and Liu, Y. Z.: A real-time look-ahead interpolation algorithm based on Akima curve fitting, Int. J. Mach. Tools Manuf., 85, 122–130, https://doi.org/10.1016/j.ijmachtools.2014.06.001, 2014.
Wang, Y. S., Yang, D. S., Gai, R. L., Wang, S. H., and Sun, S. J.: Design of trigonometric velocity scheduling algorithm based on pre-interpolation and look-ahead interpolation, Int. J. Mach. Tools Manuf., 96, 94–105, https://doi.org/10.1016/j.ijmachtools.2015.06.009, 2015.
Zhang, C. M., Zhang, P. F., and Cheng, F.: Fairing spline curves and surfaces by minimizing energy, Computer-Aided Des., 33, 913–923, https://doi.org/10.1016/S0010-4485(00)00114-7, 2001.