Journal of Control Science and Engineering

Volume 2012, Article ID 478907, 12 pages

http://dx.doi.org/10.1155/2012/478907

## Detection and Reduction of Middle-Frequency Resonance for Industrial Servo with Self-Tuning Lowpass Filter

Department of Control Science and Engineering, Huazhong University of Science and Technology, Wuhan, Hubei 430074, China

Received 2 March 2012; Revised 8 May 2012; Accepted 8 May 2012

Academic Editor: Zhiyong Chen

Copyright © 2012 Wen-Yu Wang and An-Wen Shen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

A novel method for middle frequency resonance detection and reduction is proposed for speed control in industrial servo systems. Defects of traditional resonance reduction method based on adaptive notch filter in middle frequency range are analyzed. And the main reason is summarized as the difference between the resonance frequency and the oscillation frequency. A self-tuning low-pass filter is introduced in the speed feedback path, whose corner frequency is determined by FFT results and several self-tuning rules. With the proposed method the effective range of the adaptive filter is extended across the middle frequency range. Simulation and Experiment results show that the frequency detection is accurate and resonances during the speed steady states and dynamics are successfully reduced.

#### 1. Introduction

Servo drives are used in a wide range of industrial applications including metal cutting, packaging, textiles, web-handling, automated assembly, and printing. Closed loop speed controllers with proportional-integral control laws are widely adopted in servo drives. Such controllers must be configured with high gains to achieve high performance. However, mechanical resonances are easily caused by high-speed controller gains at the presence of torsional spring connection between the motor and load [1].

Several techniques have been studied to reduce mechanical resonances. Lowering the gain of the servo driver’s speed loop is the simplest. This technique usually limits the band width of the speed loop to 1/3 of the resonance frequency or even smaller, providing very poor speed performance. Increasing the stiffness of the mechanic coupling so that the resonance frequency becomes too high to excite continuous oscillation is another choice. However, redesigning or installing new mechanical components is difficult and costly. A third scheme to reduce resonances is to introduce an additional feedback loop through the measurement of load velocity signal [2] or shaft torque [3, 4]. This loop actively dampens the resonant frequencies experienced by the load. The measured signal passes through an additional filter and is fed directly to the current command to produce a motor torque to cancel oscillations. However, measurement of the load velocity or torque output is not available in certain applications. Additional sensors on the load are also costly. The methods above are either incompetent in performance or inapplicable in hardware configuration so the following methods are more popular in recent years.

The most reported methods are based on state variables feedback, where unmeasured variables are obtained through observers, and the controller is designed with pole assignment [5–7]. However, practical problems exist such as online tuning of the observer gains, and bandwidth of the filters. With these drawbacks, this method are seldom adopted in applications with oscillation frequency higher than 200 Hz [8–10]. Another well-studied cure of mechanical resonances is adaptive notch filter based on oscillation frequency detection [11]. The structure and design are much simpler than that of the state variables feedback method, resulting in high performance in many applications with high oscillation frequency. However, the effectiveness of the notch filters relies on the accuracy of frequency detection [12, 13]. In certain cases, the results may not be accurate enough, that is, when the resonance frequency is around the cross-over frequency of the speed loop, which is within the middle-frequency range of 200~500 Hz for industrial servo. As a result, there is a blank area where both the methods above have difficulties.

In this paper, defects of traditional adaptive notch filters in middle-frequency range are analyzed, and an improved resonance frequency detection and reduction method based on self-tuning lowpass filter is proposed. Compared with traditional methods using only an adaptive notch filter, an additional self-tuning lowpass filter is introduced in the speed feedback path. The parameter of the lowpass filter is updated by FFT results and several self-tuning rules, which is shown in Figure 1. With the proposed method, the effective range of the traditional adaptive method is extended across the middle-frequency range. The blank area between state variables feedback methods and traditional adaptive notch filter methods is filled up.

#### 2. Problem Description and Existing Solution

##### 2.1. Modelling of the Compliant Coupled Two-Mass System

The torsional spring connection consists of components such as shaft joints, gearboxes, lead screws, and belt pulley sets. The mechanical stiffness of these components is limited. If the inertia of the transmission components is small compared to the motor and load, the stiffness of the components can be treated as a single, composite, equivalent spring constant that interconnects motor and load. The mechanical drive train can be modeled as a compliant coupled two-mass system. As shown in Figure 2, where is the motor torque output, and are the moments of inertia of the motor and load respectively, and is the equivalent spring constant.

A detailed block diagram of the torsional spring connected mechanism is shown in Figure 3. The equivalent spring constant is shown as providing torque to the load in proportion to the difference of motor and load positions. Also, to represent loss producing properties, a mechanical damping term is shown producing torque in proportion to velocity differences via cross-coupled viscous damping coefficient .

The transfer function from torque output of the motor to motor velocity is

The transfer function above can be viewed as a pure inertia term on the left and a dual quadratic function on the right caused by the torsional spring connection.

If the viscous damping coefficient is relatively low, both the numerator and denominator of the dual quadratic function are lightly damped. This produces a very low gain at the antiresonant frequency where the numerator is minimized and a very high gain at the resonance frequency where the denominator is minimized. The values of these frequencies are as follows:

The dual quadratic function produces severe change in both amplitude and phase of the baseline speed control system (speed control system with rigid connected load). Especially around the resonant frequency, the peak in amplitude reduces the gain margin and could easily cause instability, which is shown in Figure 4. Methods of resonance reduction rely on modifying the effects of the dual quadratic term.

##### 2.2. Defects of the Traditional Adaptive Notch Filter

Infinite impulse response notch filter is widely used for its large amplitude reduction at the resonance frequency and small phase lag at lower frequencies. A typical notch filter transfer function is as follows:

With the notch filter in the forward path, the open loop gain is eliminated at the notch frequency. With proper selection of , the gain peak caused by the torsional spring connected load could be cancelled out. The notch band width is often set low to avoid deteriorating performances around notch frequency. As a result the effectiveness of the notch filter depends on the precision of resonance frequency identification. In certain cases, frequency error of tens of hertz could render the notch filter ineffective [10]. As the mechanic characters vary from machine to machine, adaptive notch filter is proposed. In this method, FFT is performed using speed feedback signal and the peak frequency in the amplitude spectrum is treated as the notch frequency, which is shown in Figure 5.

When continuous resonance occurs in a servo system, the output of the speed loop oscillates between the positive and negative saturation values. The speed error signal is almost sinusoidal and the high-frequency noises are efficiently damped for the lowpass character of the current inner loop. The speed loop with torsional load can be analyzed with describing function method. The speed controller is treated as a proportional controller with saturation. The current inner loop and speed feedback are approximated as second-order and first-order lowpass filters, respectively.

The describing function of proportional controller with saturation is
where *a* is half the width of the linear zone and *k* is the linear gain. is the argument of the describing function. It also stands for the amplitude of the sinusoidal input signal of the nonlinear system.

Place equal to , derivative of is as follows: when > *a, so **, * is an increasing function of , and and are decreasing functions. , is a ray on negative real axis. Therefore, continuous oscillation occurs at the intersection of the Nyquist locus and the ray , that is, the −180° phase cross-over frequency of the open-loop Bode diagram [14].

The open loop phase lag increases rapidly with frequency and produces a −180° phase cross-over if the spring constant is relatively large and the damping coefficient is small. As a result the oscillation frequency is very close to the resonance frequency and the notch filter automatically designed with the FFT result works well. However, in the following cases, the FFT result may not be accurate enough.

*Case 1 (Higher Middle-Frequency Resonance with Large Damping Coefficient). *In this case the resonance frequency is slightly lower than the −180° phase cross-over frequency, where the phase lag increases slowly as a result of the large damping coefficient. That results in the difference between the cross-over frequency and the resonance frequency. There will be continuous oscillation at the cross-over frequency if the speed controller gain is large enough. Traditional adaptive notch filter will set the notch frequency the same as the cross-over frequency, which will be higher than the actual resonance frequency. For example, as shown in Figure 6, the notch/cross-over/oscillation frequency is about 500 Hz, while the resonance frequency is about 400 Hz. In that case, the traditional adaptive notch filter may fail in suppressing the amplitude peak, or even deteriorates the oscillation, which is shown in Figure 7.

*Case 2 (Lower Middle-Frequency Resonance). *In this case, the resonance frequency is much lower than the −180° phase cross-over frequency, leaving continuous oscillation almost impossible, which is shown in Figure 8. However, there are jitters during the speed steps during the speed controller output saturation. The jitters fade immediately the controller recovers from saturation when the speed reaches its reference, which is shown in Figure 9. Sampling of the speed error signal is usually performed at the same rate with speed-loop calculation, and a total amount of over 512 samples are necessary for precision of the FFT result. That means at least 100 milliseconds is necessary for one single FFT, while the speed jitters fade too quickly to be sampled. In that case, the traditional adaptive notch filter could neither find an effective notch frequency nor suppress the speed jitters during the speed steps.

#### 3. The Proposed Method Based on Self-Tuning Lowpass Filter

In the cases above, FFT failed to provide accurate notch frequency for the notch filter because of the difference in the oscillation frequency and the actual resonance frequency. The difference could be smaller if the phase cross-over frequency can be lowered deliberately. The effective method of lowering the cross-over frequency is adding a lowpass filter in the control loop. The method should meet the requirements below.(1)The method should generate enough phase lag in order to bring the stable system into continuous oscillation so that Case 2 can be treated as Case 1 afterward. (2)The final oscillation frequency should be close enough with the resonance frequency.

For requirement 1, the filter’s phase lag should increase fast enough around the desired frequency. An IIR Butterworth filter is applicable for its nonlinear phase-frequency characteristic. The order of the filter should be high enough, but not too high to bring excessive amplitude attenuation that the system remains stable, as in some applications with low load inertia ratio. The corner frequency should be relatively low enough, but not too low to induce low-frequency oscillation, whose frequency is 1/3 or lower of the resonance frequency.

A second-order lowpass filter with corner frequency the same as the speed loop cross-over frequency meets requirement 1. The largest possible phase margin of the speed loop is 90° for the load inertia which is modeled as pure integral. A second-order lowpass filter lowers the gain margin by 90° at its corner frequency, guaranteeing a phase cross-over around the resonance frequency, which is shown in Figures 10 and 11.

For requirement 2, self-tuning rules could be applied to the second-order lowpass filter to further reduce the frequency difference. The corner frequency can be updated with the oscillation frequency so that the cross-over frequency is further lowered. Since the phase lag increases rapidly around the resonance frequency, the phase oscillation frequency changes more slowly as the frequency difference becomes smaller, as shown in Figure 12. In most cases, the frequency difference will diminish after the first update and stay below the resolution of FFT with further updates which is shown in Figure 13. Generally, up to two updates are enough for precision of resonance frequency identification. More updates could possibly take extra time without a better result, thus it is not favourable.

The proposed resonance detection method based on self-tuning lowpass filter is summarized as follows: The first FFT is performed with the speed error signal once the detection begins. The notch filter is designed straight with the oscillation frequency if it is much higher than the baseline speed loop cross-over frequency. If the oscillation frequency is around the cross-over frequency or no effective resonance peak is detected, a second-order lowpass filter is added to the speed feedback with the baseline speed loop cross-over frequency as its corner frequency. Another FFT is performed and the corner frequency is updated with the new oscillation frequency. The third FFT is performed calculating the final oscillation frequency, which the notch filter is finally designed with. The system block diagram of the speed control loop with the proposed method is shown in Figure 1.

#### 4. Simulation of the Proposed Method

A preliminary simulation analysis has been performed to verify the proposed method using MATLAB/Simulink. The simulation model of the servo driver has been set up in accordance with Figure 1. The current loop is modeled as a second-order lowpass filter and the compliant coupled system is modeled with transfer functions in (1) for simplicity.

Parameters of the simulation are listed in Table 1.

Two groups of simulations have been performed in accordance with the cases in Section 2 with different spring constants, viscous damping coefficient and ratio of load/motor inertia. and indicate the oscillation frequencies after each parameter update.

Group 1, viscous damping coefficient was 0.03 and ratio of load/motor inertia was 1/1, corresponding to Case 1 in Section 2 where the resonance frequency and the damping coefficient are both relatively high. Two different spring constants were applied with the same speed controller parameter 0.85 A/rad/s and 6 ms. The initial oscillation frequency was much higher than the resonance frequency before the introduction of the lowpass filter. The final oscillation frequency was very close to the calculated resonance frequency after two updates. Results of simulation group 1 are shown in Table 2.

The detail of the initial oscillation in test number 1 is shown in Figure 14. The frequency of the detected oscillation is 346 Hz, much higher than the actual resonance frequency 302 Hz.

After the notch frequency detection, the lowpass filter is removed. A notch filter described in Section 2 is introduced in test number 1 using the result resonance frequency 302 Hz and notch band width 0.2. The oscillations faded quickly after the introduction of the notch filter. As shown in Figure 15, the jitters during the speed dynamics are also reduced effectively.

Group 2, viscous damping coefficient was 0.03 and ratio of load/motor inertia was 3/1, corresponding to Case 2 in Section 2 where the resonance frequency is relatively low. Three different spring constants were applied with the same speed controller parameter as 1.2 A/rad/s and as 5 ms. No continuous oscillation was detected before the introduction of the lowpass filter. The final oscillation frequency is very close to the calculated resonance frequency after two updates. Results of simulations group 2 is listed in Table 3.

The detail of the speed step response in test number 1 is shown in Figure 16. The system remains stable because the resonance frequency is much lower than the speed loop base-line cross-over frequency. But the speed jitters were quite obvious during the speed step.

After the notch frequency detection, the lowpass filter was removed. A notch filter described in Section 2 was introduced in test number 1 using the detected resonance frequency 196 Hz and notch stop band width 0.2. Most of the jitters during the speed dynamics are reduced effectively, which is shown in Figure 17.

#### 5. Experiment of the Proposed Method

##### 5.1. Hardware Configuration

The experiment was performed on a test platform consisting of a permanent magnet servo motor (PMSM) servo system and a torsional spring connection. The servo driver was controlled with a Texas Instrument 32 bit fixed-point DSP (TMS320LF28035). The DSP performed all the control algorithms under a clock frequency of 60 MHz. A surface mount PMSM with incremental encoder was used as the driving motor and an identical one as load. The two motors were connected with KTR-ROTEX torsional flexible coupling. The spring constant could be adjusted by altering the depth of the two coupling claws. The load/motor inertia ratio could be adjusted with additional inertia plates. Figure 18 shows the configuration in Experiment 2, where the inertia plate was installed and the two coupling claws were adjusted to only about 1/3 of the full depth. Parameters of the test platform hardware are listed in Table 4.

##### 5.2. Software Configuration

The proposed method was implemented on the basis of a classic double loop speed control system, which is shown in Figure 1. The inner current loop was implemented with field orient control (FOC), which was most commonly used in PMSM control. Since this paper focuses on the performance of the speed loop only, further details of the current loop are not necessary. The speed controller with PI control law was performed at an interval of 3 ms. Its output saturation value was set to 6 A, the same as the rated current of the motor. The motor speed feedback signal was measured by computing differential of the rotor-position at the same interval as the speed controller. The speed and current data were transmitted to PC through synchronous serial interface. FFT was performed every 512 speed control cycles with 512 samples of speed error. The parameters of the lowpass filter were updated once the peak frequency in the amplitude spectrum was identified. In frequencies 2 and 3, two updates were performed before the final identification of the resonance identification. The notch filter was designed with the identified notch frequency and bandwidth 0.1.

The baseline speed loop cross-over frequency was verified through open loop frequency sweep test. In this test, the speed feedback was disconnected and the motor runs under sinusoidal speed references without rigid connected loads. The frequency of reference speed increased until the phase difference between the reference speed and the motor speed reached 180 degree. Parameters of the test platform software are listed in Table 5.

##### 5.3. Speed Profile Selection, Controller Parameter Tuning

In industry applications the trapezoidal profile and s-curve are widely used for they are less aggressive and also energy saving than step reference. However, in order to test the effectiveness of the proposed method, a harsh experiment condition is necessary. As a result, aggressive step speed references were adopted in all the experiments below. Parameters of the controller were also tuned for maximum performance until generating instability. The actual resonance frequency was verified by a simple test in which we applied an impact on the motor axis with hammer and the speed oscillation was analyzed with FFT. The −3 dB bandwidth was tested using sinusoidal speed references with amplitude 5% of the rated motor speed.

##### 5.4. Experiment Results

Three experiments were performed, with the first two corresponding to Case 1 and the third corresponding with Case 2 in Section 2.2. Experiment 1 provided the performance of the system without notch filter. Both Experiments 2 and 3 performed the entire process of resonance identification and reduction. Results of the experiments are as follows.

*Experiment 1. *The load/motor inertia was 1/1, and the spring constant of the torsional flexible coupling was adjusted to around 1800 Nm/rad, resulting in a resonance frequency of 302 Hz. The system had no notch filter. The speed loop proportional gain was raised as high as possible without generating instability ( A/rad/s). The speed loop integral gain was then raised until a step command generated 25% overshoot ( ms). The step response is shown in Figure 19. Note that the resonance is fading during the saturation of the speed controller (before 20 ms). However, when the controller recovered from saturation, tendency of instability appeared (20 to 60 ms). Settling time was about 60 ms. The −3 dB bandwidth was measured as 30 Hz.

*Experiment 2. *The experiment was performed again with the same mechanical configuration in Experiment 1. The speed loop proportional gain was raised until a continuous oscillation was detected. As it is shown in Figure 20, the frequency of the detected oscillation was 344 Hz. This frequency was much higher than the resonance frequency 302 Hz. The frequency difference of 42 Hz would easily make traditional adaptive notch filter method ineffective. The lowpass filter was introduced at 0.35 s. Its corner frequency was set the same as the first detected oscillation (344 Hz). The lowpass filter was updated at 0.55 s, lowering the oscillation frequency to 304 Hz. The lowpass filter was removed and the notch filter was introduced at 0.7 s. The oscillations were immediately suppressed. Then, the proportional gain was further raised as high as 0.85 A/rad/s and was raised as high as 6 ms. With this configuration, the speed overshoot under a step command was 10% and settling time was about 20 ms, which is shown in Figure 21. The −3 dB bandwidth was measured as 75 Hz. The waveforms of the entire process are shown in Figure 22. Note that the amplitude of speed oscillation became larger at 0.35 s, which showed that the updated oscillation frequency was closer to the resonance peak.

With the same parameters as test number 1 in simulation group 1, we can see that Figures 20 and 21 are quite similar with Figures 14 and 15 in Section 4. The continuous oscillation in both simulation and experiment are effectively suppressed with the proposed method. The main difference in that the noises in Figure 21 are much bigger than that in Figure 15 as a result of the limited resolution of the encoder and the aggressive controller parameter.

*Experiment 3. *An inertia plate was added to the load motor, making the load-motor inertia ratio 3/1. The spring constant of the torsional flexible coupling was adjusted to around 1200 Nm/rad, resulting in a resonance frequency of 198 Hz. The proportional gain was raised to an aggressive value of 1.2 A/rad/s and was raised as high as 5 ms. The system remained stable because the resonance frequency was much lower than the speed loop baseline cross-over frequency. But the speed jitters were quite obvious during the first speed step, which is shown in Figure 23. Note that the speed jitters faded quickly within 20 ms. The traditional adaptive notch filter method would have difficulties in detecting these jitters. The lowpass filter was introduced at 0.4 s. Its corner frequency was set the same as the speed loop baseline cross-over frequency, inducing a continuous oscillation. The lowpass filter was updated at 0.55 s, lowering the oscillation frequency to 195 Hz. The lowpass filter was removed and the notch filter was introduced at 0.7 s. The speed jitters were successfully suppressed during the next speed step, which is shown in Figure 24. The waveforms of the entire process are shown in Figure 25. The −3 dB bandwidth was measured as 39 Hz. Note that the amplitude of speed oscillation became larger at 0.55 s, which showed that the updated oscillation frequency was closer to the resonance peak.

With the same parameters as test number 1 in simulation group 2, we can see that Figures 23 and 24 are quite similar with Figures 16 and 17. The speed jitters in both simulation and experiment are effectively suppressed with the proposed method. The main difference in that the noises in Figure 24 are much bigger than that in Figure 17 as a result of the limited resolution of the encoder and the aggressive controller parameter.

#### 6. Conclusion

In this paper, an improved resonance frequency detection and reduction method based on self-tuning lowpass filter is proposed. Introduction of the self-tuning lowpass filter minimize, the difference between the oscillation frequency and the resonance frequency. The effective range of the adaptive notch filter is extended well below the speed loop cross-over frequency even with considerable damping coefficient. The middle-frequency blank area between state variables feedback methods and traditional adaptive notch filter methods is filled up. Simulation and experimental results fully confirm the validity of the proposed method. The method is simple and effective. Neither additional hardware nor complex parameter tuning is required. Resonances during the speed steady states and dynamics are successfully reduced.

#### Acknowledgments

The authors wish to acknowledge the motivation and advisement given by Motion Control Laboratory of the Huazhong University of Science and Technology. The authors also would like to extend their thanks to The Hangzhou Riding Control Technology Ltd. for the financial support.

#### References

- G. Zhang and J. Furusho, “Speed control of two-inertia system by PI/PID control,”
*IEEE Transactions on Industrial Electronics*, vol. 47, no. 3, pp. 603–609, 2000. View at Publisher · View at Google Scholar · View at Scopus - S. N. Vukosavic and M. R. Stojic, “Suppression of torsional oscillations in a high-performance speed servo drive,”
*IEEE Transactions on Industrial Electronics*, vol. 45, no. 1, pp. 108–117, 1998. View at Google Scholar · View at Scopus - K. Sugiura and Y. Hori, “Vibration suppression in 2- and 3-mass system based on the feedback of imperfect derivative of the estimated torsional torque,”
*IEEE Transactions on Industrial Electronics*, vol. 43, no. 1, pp. 56–64, 1996. View at Google Scholar · View at Scopus - K. Szabat and T. Orlowska-Kowalska, “Comparative analysis of different PI/PID control structures for two-mass system,” in
*Proceedings of the 7th International Conference on Electric and Electronics Equipment*, pp. 97–102, Brasov, Romania, 2004. - Y. Hori, H. Sawada, and Y. Chun, “Slow resonance ratio control for vibration suppression and disturbance rejection in torsional system,”
*IEEE Transactions on Industrial Electronics*, vol. 46, no. 1, pp. 162–168, 1999. View at Google Scholar · View at Scopus - T. Ohmae, T. Mastuda, M. Kanno, and K. Saito, “Microprocessor-based motor speed controller using fast-response state observer for reduction of torsional vibration,”
*IEEE Transactions on Industry Applications*, vol. 23, no. 5, pp. 863–871, 1987. View at Google Scholar - Y. Hori, H. Iseki, and K. Sugiura, “Basic consideration of vibration suppression and disturbance rejection control of multi-inertia system using SFLAC (State feedback and load acceleration control),”
*IEEE Transactions on Industry Applications*, vol. 30, no. 4, pp. 889–896, 1994. View at Publisher · View at Google Scholar · View at Scopus - T. Orlowska-Kowalska and K. Szabat, “Sensitivity analysis of state variable estimators for two-mass drive system,”
*Acta Electrotechnica et Informatica*, vol. 4, no. 1, 2004. View at Google Scholar - T. M. O'Sullivan, C. M. Bingham, and N. Schofield, “High-performance control of dual-inertia servo-drive systems using low-cost integrated SAW torque transducers,”
*IEEE Transactions on Industrial Electronics*, vol. 53, no. 4, pp. 1226–1237, 2006. View at Publisher · View at Google Scholar · View at Scopus - G. Ellis and R. D. Lorenz, “Resonant load control methods for industrial servo drives,” in
*Proceedings of the 35th IAS Annual Meeting and Wolrd Conference on Industrial Applications of Electrical Energy*, pp. 1438–1445, October 2000. View at Scopus - P. Schmidt and T. Rehm, “Notch filter tuning for resonant frequency reduction in dual inertia systems,” in
*Proceedings of the 34th IAS Annual Meeting on IEEE Industry Applications Conference*, pp. 1730–1734, October 1999. View at Scopus - C. I. Kang and C. H. Kim, “An adaptive notch filter for suppressing mechanical resonance in high track density disk drives,”
*Microsystem Technologies*, vol. 11, no. 8-10, pp. 638–652, 2005. View at Publisher · View at Google Scholar · View at Scopus - H. Wang, D. H. Lee, Z. G. Lee, and J. W. Ahn, “Vibration rejection scheme of servo drive system with adaptive notch filter,” in
*Proceedings of the 37th IEEE Power Electronics Specialists Conference (PESC '06)*, pp. 1–6, June 2006. View at Publisher · View at Google Scholar · View at Scopus - A. Gelb and W. E. Vander Velde,
*Multiple-Input Describing Functions and Nonlinear System Design*, McGraw Hill, 1968.