# Research of Mechanical Resonance Analysis and Suppression Control Method of the Servo Drive System

Authors: Wenli Li; Yongkang Liu; Shuaishuai Ge; Daming Liao
Journal: Shock and Vibration
 (2021)

---

## Abstract

Transmission mechanisms of the servo drive system are not a pure rigid body, and the existence of the elastic transmission mechanisms will make the system generate mechanical resonance. Aiming at mechanical resonance of the servo drive system, the resonance generation mechanism is analyzed, the four-mass model considering the time-varying meshing stiffness of the gear is established, and the influence of different stiffness parameters on the mechanical resonance of the system is researched in this paper. The composite controller of Model Predictive Control (MPC) with Notch Filter is used to simulate the mechanical resonance suppression of the four-mass servo system considering time-varying meshing stiffness, and it is compared with the mechanical resonance suppression method based on Model Predictive Control. The simulation results show that when the step speed is 200 r/min, the overshoot is reduced from 11.6 r/min to 1.1 r/min, which is reduced by 90.5%. Under the impact load condition, from 10 Nm to 30 Nm, overshoot is reduced from 34.3 r/min to 12.8 r/min, reduced by 62%, and torque oscillation is reduced by 81.5%. Therefore, the composite controller of Model Predictive Control with Notch Filter can suppress the mechanical resonance problem effectively, caused by elastic transmission, and improve the robustness of servo drive system.


---

## Body

## 1. Introduction

Servo drive system is a common electromechanical system, which has the relationship of multiphysical processes and multicomplex parameter coupling. It is widely used in industrial fields such as industrial robots and hydraulic machine tools [1, 2]. During the operation of the servo motor, the torque ripple is transmitted to load through the transmission shaft, gear, and other transmission mechanisms, which will lead to vibration of the drive motor and transmission mechanisms. Owing to the existence of backlash and other factors, vibration will also be produced and transmitted to the drive motor, resulting in the occurrence of coupling vibration and affecting the life of the motor in the process of transmission. Therefore, it is particularly important to deeply analyze and suppress the vibration mechanism of servo drive system.In early period, the mechanical vibration suppression control method is relatively straightforward, which is to avoid the sudden change of speed. The elastic effect of damping by adding a term is proportional to the differential speed (difference between connecting rod speed and motor speed) simply [3, 4]. However, for reasons of efficiency and cost, it is rarely applied to the industrial field [5].At present, mechanical vibration suppression methods mainly focus on passive suppression and active suppression. The passive suppression method does not need to import external energy and change the control system structure. Vibration reduction of servo system is by damping energy dissipation device, vibration isolation, or filtering. The frequently used methods include passive damping vibration reduction, insertion of filter, and improvement of the inertia ratio of motor to load [6–9]. In the passive suppression method, it is difficult to improve the inertia ratio of motor-load. Low-pass filter will bring a certain phase delay, which is easy to make the dynamic characteristics of the system to deteriorate. The insertion of Notch Filter is easy, and it does not require additional feedback links in the system and can achieve better vibration suppression effect [10]. The active suppression method relies on external given command and system response indicators and changes the controller parameters and structure to weaken the mechanical vibration of servo system [11]. Benefiting from the improvement of basic science and computer science, active suppression methods have been developed rapidly, mainly including state observer, active disturbance rejection control, sliding mode control, fuzzy control, and Model Predictive Control [12–15]. In the active suppression method, fuzzy control is easy to simplify complex information and result in the degradation of the precision of the system. Sliding mode control is greatly affected by the switching amplitude of control variables, which is not suitable for the precise control system.Model Predictive Control can obtain accurate models through analysis techniques and identification methods, which is an effective solution to nonlinear and constrained problems [16]. Model Predictive Control can realize the control and simultaneous optimization of multidimensional variables. The speed regulator and current regulator of the cascade structure are combined into a noncascaded multi-input and multioutput Model Predictive Controller, which can realize the current at the same time in one calculation cycle. The control system combines the three regulators of the original classic control structure into one controller and realizes simultaneous control of current and speed, which significantly simplifies the system structure and eliminates the current hysteresis link that needs to be considered when designing a speed regulator. The noncascaded structure can further improve its dynamic response while retaining the optimized performance of Model Predictive Control and obtain better control performance.The actual flexible connection part of the servo drive system is very complex. In order to facilitate the analysis and research, the existing research usually assumes that the transmission mechanisms of the flexible link do not have elastic deformation, and the mass is evenly distributed on a long connection shaft, so as to establish a dual inertia system model [17, 18]. The research of mechanical vibration mechanism of double inertia system is not deep enough, and the vibration of transmission device is not considered in detail. In this paper, starting from the coupling mechanism of mechanical and electrical, taking the four-mass servo drive system considering the time-varying meshing stiffness of gear as the research object, the genesis and characteristics of vibration are expounded. On the basis of this, the use of Model Predictive Control does not require accurate mathematical models and can deal with various problems with constraints. The advantages of Notch Filter that can rapidly attenuate the import signal at a specific frequency point. A composite controller combining Model Predictive Control with Notch Filter is used to suppress the mechanical resonance generated by the servo drive system and reduce the damage caused by the mechanical resonance during the operation of the four-mass servo drive system.


## 2. Model of the Servo System and Vibration Analysis

As shown in Figure1, the servo drive system is a complex nonlinear system composed of servo motor, transmission system, and load. Considered on the influence of some nonlinear factors of mechanical system, the servo drive system can be regarded as a four-mass system composed of rigid and flexible components. The concentrated mass method is used to establish the servo drive system model considering the single stage spur gear transmission mechanisms.Figure 1
Servo drive system.


### 2.1. Dynamic Model of the Servo System

There are not great differences in magnitude of inertia between motor, gear drive, and load in servo system. Considered on the torsional stiffness and damping of the drive shaft and the moment of inertia of the motor rotorJe, gear transmission device J1, J2 and load JL comprehensively, the motor rotor, load, and gear transmission device can be regarded as four parts of moment of inertia. Based on Lagrange dynamics analysis and Hooke's law, do the setup as shown in Figure 2 of the servo drive system dynamics model, as shown inwhere θe is the rotation angle of the equivalent motor shaft; θ1 and θ2 are, respectively, the rotation angles of the main and driven gears in the gear pair; θL is the load rotation angle; Ts1 and Ts2 are, respectively, the torque borne by the drive shaft; C1 and C2 are, respectively, the torsional damping values of input and output shafts; K1 and K2 are, respectively, the torsional stiffness of the input and output shafts; Te and TL are the electromagnetic torque and load torque of the motor, respectively; R1 and R2 are the radius of the gear base circle; and Wd is the gear dynamic meshing force.Figure 2
Electromechanical coupling dynamics model.(1)Jeθ¨e=Te−Ts1,J1θ¨1=Ts1−R1Wd,J2θ¨2=R2Wd−Ts2,JLθ¨L=Ts2−TL,Ts1=C1θ˙e−θ˙1+K1θe−θ1,Ts2=C2θ˙2−θ˙L+K2θ2−θL,


### 2.2. Torsion Dynamics Model of Gear

Only considering the torsional vibration of the gear, ignoring the lateral and axial elastic deformation of the transmission shaft and the elastic deformation of the supporting system, the torsional vibration analysis model of the single stage spur gear transmission is established as shown in Figure3.Figure 3
Torsional vibration model of gear transmission.The torsional dynamics analysis model of the gear is shown in(2)J1θ¨1+R1CmR1θ˙1−R2θ˙2−e˙t+KmR1fR1θ1−R2θ2−et=T1,J2θ¨2−R2CmR1θ˙1−R2θ˙2−e˙t−KmR2fR1θ1−R2θ2−et=−T2,where Cm is the meshing damping of the gear pair; Km is the time-varying meshing stiffness in the gear transmission process; e(t) is the comprehensive error of the gear meshing; and T1 and T2 are the main and driven gears to bear the torque.According to the aforementioned gear dynamic model, the dynamic meshing force acting on thei-th tooth is shown in(3)Wd=CmR1θ˙1−R2θ˙2−e˙t+KmfR1θ1−R2θ2−et.


#### 2.2.1. Calculation of Time-Varying Meshing Stiffness

In the process of gear meshing, single and double tooth pairs of gears alternate with the meshing periodically, which results in periodic changes in the comprehensive meshing stiffness of the gears, which causes dynamic excitation of the gear system. The load required to produce a deflection of 1μm on an average of 1 mm tooth width is called the meshing stiffness of the gear teeth [19]. Therefore, the establishment of a more practical stiffness model and the study of meshing stiffness are of great significance. At present, the meshing stiffness model is mainly obtained by the following two methods:(1)
Formula Methods. This method mainly includes Ishikawa formula and its improved algorithm, potential energy method and improved energy method, and Weber formula.(2)
Finite Element Method. The finite element method is the most widely used method currently. An empirical formula that can approximate the stiffness can be obtained through the calculated finite element analysis results.The potential energy method [20] uses the idea of variable increment infinite approximation to calculate the position of each meshing point in the gear meshing process, mapping each meshing point position to a single tooth profile, and using the material mechanics cantilever beam theory to calculate the time-varying meshing stiffness of the gear, which has high calculation efficiency. The result is accurate and includes other advantages. In this paper, the potential energy method is used to solve the time-varying meshing stiffness of gears. During the gear meshing process, the force on the gear teeth can be equivalent to the force on the cantilever beam. The equivalent force model of the spur gear is shown in Figure 4. In the figure, F is the normal load received by the spur gear, and α is the pressure angle. h is half the tooth thickness, where the load F is applied; d is the distance from the root circle to the addendum circle.Figure 4
Equivalent model of the spur gear.According to the knowledge of mechanics, the deformation of the gear teeth under the action of the normal loadF on the spur gear can be equivalently decomposed into bending deformation, shear deformation, and radial deformation along the tooth height direction. Based on the theory of material mechanics, the bending stiffness Kb, the shear stiffness Ks, the radial compression stiffness Kr, the stiffness Kd produced by the deformation of the gear teeth, and the stiffness Kh produced by the Hertz contact deformation can be obtained. Therefore, the meshing stiffness of single tooth pair of spur gear pair is expressed as(4)1Km=1∑i=121/Kbi+1/Ksi+1/Kri+1/Kdi+1/Khi,where(5)1kb=∫α1−α231+cosα1α2−αsinα−cosα2α2−αcosα2ELsinα+α2−αcosα3dα+∫0Rb−Rdd+x1cosα1−hsinαEIx1dx1Z<42∫α1−α331+cosα1α2−αsinα−cosα2α2−αcosα2ELsinα+α2−αcosα3dαZ≥42,1ks=∫−α1α21.21+μα2−αcosαcos2α1ELsinα+α2−αcosαdα+∫0Rb−Rd1.2cosα12GAx1dx1Z<42∫−α1α31.21+μα2−αcosαcos2α1ELsinα+α2−αcosαdαZ≥42,1kd=cos2αELLμfSf2+MμfSf+P1+Qtan2α,1kr=∫−α1α2α2−αcosαsin2α12ELsinα+α2−αcosαdα+∫0Rb−Rdsinα12EAx1dx1Z<42∫−α1α3α2−αcosαsin2α12ELsinα+α2−αcosαdαZ≥42,where i = 1, 2 are the main and driven gears, respectively, α is the pressure angle, G is the shear modulus of elasticity, and Ax = 2hxL is the gear cross-sectional area.Based on the distribution of the aforementioned gear pairs in a meshing period, MATLAB programming can be used to obtain the single-double-tooth meshing comprehensive stiffness of the spur gear pair in a period as shown in Figure5.Figure 5
Time-varying mesh stiffness of gears.


### 2.3. Mechanical Resonance Analysis

Servo drive system has many vibration characteristics during operation. This paper mainly discusses the characteristic vibration of torsional vibration generated by the electromechanical system.The mechanical rotor is an elastic body, so in the process of rotation, various rotating parts will form different degrees of torsional vibration. Torsional vibration may originate from the instability of the motor output or load, or it may originate from the internal excitation of the transmission system. These external and internal disturbances lead to forced vibration or free vibration of the mechanical system.During the operation of the servo drive system, the drive motor directly drives the input shaft to rotate, and the shaft and gear pair are transmitted to the load. During the rotation process, each component will produce a certain degree of elastic deformation, forming an angular deviation, which makes the mechanical transmission system torsionally vibrate. This torsional vibration is mainly reflected in the fluctuation of angle, angular velocity, and angular acceleration.


#### 2.3.1. Mechanical Resonance Analysis with Fixed Stiffness

With analysis of time-varying meshing stiffness, the average meshing stiffness of the gear is 3.8 × 108 N/m. Therefore, the meshing stiffness of the gear pair is equivalent to a fixed stiffness of 3.8 × 108 N/m. The Bode diagram shown in Figure 6 of the four-mass servo system can be obtained using Equation (1).Figure 6
Bode diagram of the system with fixed stiffness.The four-mass servo system has three resonant frequency points, which are the input shaft, output shaft, and gear pair in the system, and all introduce the corresponding mechanical resonance to the system. The transfer function between the system motor speed and electromagnetic torque can be expressed as(6)ωesTes=1Jes∏n=13s2+2ξanωan+ωan2s2+2ξbnωbn+ωbn2,where ωbn and ωan are the resonant frequency and antiresonant frequency of the nth resonant link, respectively; and ξbn and ξan are the resonant damping coefficient and antiresonant damping coefficient of the nth resonant link, respectively.The transfer function of four-mass system that brings mechanical resonance is given by Equation (6):(7)Gs=∏n=13s2+2ξanωan+ωan2s2+2ξbnωbn+ωbn2.The analysis of Equation (6) shows that the resonance frequency of the system is related to the stiffness coefficient of the transmission mechanism and the equivalent inertia of each part of the system.In order to investigate the speed control performance of the load, assume that the meshing fixed stiffnesses are 3 × 107, 3 × 108, and 3 × 109 N/m, and the rest of the system parameters remain unchanged. The Bode diagram of the input to the output with the four-mass servo system can be obtained as shown in Figure 7. It can be seen that the change of the fixed stiffness value will change the resonant frequency of the system; that is, the resonant frequency value of the system will increase with the increase of the fixed stiffness.Figure 7
Bode diagram of the servo system with different stiffnesses.


#### 2.3.2. Analysis of Mechanical Resonance with Time-Varying Meshing Stiffness

In order to further investigate the influence of the time-varying meshing stiffness on the mechanical resonance of the four-mass servo system, the change range of the resonance frequency within the range of the time-varying meshing stiffness is analyzed. Here, the different stiffness values (2.7 × 108,3.8 × 108 and 4.9 × 108 N/m) during the gear meshing process are taken. The Bode diagram of the input electromagnetic torque to the output load speed can be obtained as shown in Figure 8. It can be seen that the mechanical resonance frequency range brought about by the time-varying meshing stiffness is from 5966 Hz to 7638 Hz.Figure 8
Bode diagram of electromagnetic torque to load speed.According to Equation (1), a four-mass servo system model considering the time-varying meshing stiffness of gear is established. The servo system model established is compared and analyzed under the conditions of time-varying meshing stiffness and constant stiffness. When a speed command of 1000 r/min is given, the time domain characteristics of the torque transmitted at the load are obtained in Figure 9, and when considering the time-varying stiffness parameters, the system load transfer torque will produce greater fluctuations than considering the constant stiffness parameters. From the resonance analysis of the system under constant stiffness and variable stiffness parameters, it can be concluded that it is necessary to consider the gear meshing stiffness as a time-varying meshing stiffness when we research four-mass servo system.Figure 9
Transmit torqueTs2 change of the servo system under constant stiffness and variable stiffness.


## 2.1. Dynamic Model of the Servo System

There are not great differences in magnitude of inertia between motor, gear drive, and load in servo system. Considered on the torsional stiffness and damping of the drive shaft and the moment of inertia of the motor rotorJe, gear transmission device J1, J2 and load JL comprehensively, the motor rotor, load, and gear transmission device can be regarded as four parts of moment of inertia. Based on Lagrange dynamics analysis and Hooke's law, do the setup as shown in Figure 2 of the servo drive system dynamics model, as shown inwhere θe is the rotation angle of the equivalent motor shaft; θ1 and θ2 are, respectively, the rotation angles of the main and driven gears in the gear pair; θL is the load rotation angle; Ts1 and Ts2 are, respectively, the torque borne by the drive shaft; C1 and C2 are, respectively, the torsional damping values of input and output shafts; K1 and K2 are, respectively, the torsional stiffness of the input and output shafts; Te and TL are the electromagnetic torque and load torque of the motor, respectively; R1 and R2 are the radius of the gear base circle; and Wd is the gear dynamic meshing force.Figure 2
Electromechanical coupling dynamics model.(1)Jeθ¨e=Te−Ts1,J1θ¨1=Ts1−R1Wd,J2θ¨2=R2Wd−Ts2,JLθ¨L=Ts2−TL,Ts1=C1θ˙e−θ˙1+K1θe−θ1,Ts2=C2θ˙2−θ˙L+K2θ2−θL,


## 2.2. Torsion Dynamics Model of Gear

Only considering the torsional vibration of the gear, ignoring the lateral and axial elastic deformation of the transmission shaft and the elastic deformation of the supporting system, the torsional vibration analysis model of the single stage spur gear transmission is established as shown in Figure3.Figure 3
Torsional vibration model of gear transmission.The torsional dynamics analysis model of the gear is shown in(2)J1θ¨1+R1CmR1θ˙1−R2θ˙2−e˙t+KmR1fR1θ1−R2θ2−et=T1,J2θ¨2−R2CmR1θ˙1−R2θ˙2−e˙t−KmR2fR1θ1−R2θ2−et=−T2,where Cm is the meshing damping of the gear pair; Km is the time-varying meshing stiffness in the gear transmission process; e(t) is the comprehensive error of the gear meshing; and T1 and T2 are the main and driven gears to bear the torque.According to the aforementioned gear dynamic model, the dynamic meshing force acting on thei-th tooth is shown in(3)Wd=CmR1θ˙1−R2θ˙2−e˙t+KmfR1θ1−R2θ2−et.


### 2.2.1. Calculation of Time-Varying Meshing Stiffness

In the process of gear meshing, single and double tooth pairs of gears alternate with the meshing periodically, which results in periodic changes in the comprehensive meshing stiffness of the gears, which causes dynamic excitation of the gear system. The load required to produce a deflection of 1μm on an average of 1 mm tooth width is called the meshing stiffness of the gear teeth [19]. Therefore, the establishment of a more practical stiffness model and the study of meshing stiffness are of great significance. At present, the meshing stiffness model is mainly obtained by the following two methods:(1)
Formula Methods. This method mainly includes Ishikawa formula and its improved algorithm, potential energy method and improved energy method, and Weber formula.(2)
Finite Element Method. The finite element method is the most widely used method currently. An empirical formula that can approximate the stiffness can be obtained through the calculated finite element analysis results.The potential energy method [20] uses the idea of variable increment infinite approximation to calculate the position of each meshing point in the gear meshing process, mapping each meshing point position to a single tooth profile, and using the material mechanics cantilever beam theory to calculate the time-varying meshing stiffness of the gear, which has high calculation efficiency. The result is accurate and includes other advantages. In this paper, the potential energy method is used to solve the time-varying meshing stiffness of gears. During the gear meshing process, the force on the gear teeth can be equivalent to the force on the cantilever beam. The equivalent force model of the spur gear is shown in Figure 4. In the figure, F is the normal load received by the spur gear, and α is the pressure angle. h is half the tooth thickness, where the load F is applied; d is the distance from the root circle to the addendum circle.Figure 4
Equivalent model of the spur gear.According to the knowledge of mechanics, the deformation of the gear teeth under the action of the normal loadF on the spur gear can be equivalently decomposed into bending deformation, shear deformation, and radial deformation along the tooth height direction. Based on the theory of material mechanics, the bending stiffness Kb, the shear stiffness Ks, the radial compression stiffness Kr, the stiffness Kd produced by the deformation of the gear teeth, and the stiffness Kh produced by the Hertz contact deformation can be obtained. Therefore, the meshing stiffness of single tooth pair of spur gear pair is expressed as(4)1Km=1∑i=121/Kbi+1/Ksi+1/Kri+1/Kdi+1/Khi,where(5)1kb=∫α1−α231+cosα1α2−αsinα−cosα2α2−αcosα2ELsinα+α2−αcosα3dα+∫0Rb−Rdd+x1cosα1−hsinαEIx1dx1Z<42∫α1−α331+cosα1α2−αsinα−cosα2α2−αcosα2ELsinα+α2−αcosα3dαZ≥42,1ks=∫−α1α21.21+μα2−αcosαcos2α1ELsinα+α2−αcosαdα+∫0Rb−Rd1.2cosα12GAx1dx1Z<42∫−α1α31.21+μα2−αcosαcos2α1ELsinα+α2−αcosαdαZ≥42,1kd=cos2αELLμfSf2+MμfSf+P1+Qtan2α,1kr=∫−α1α2α2−αcosαsin2α12ELsinα+α2−αcosαdα+∫0Rb−Rdsinα12EAx1dx1Z<42∫−α1α3α2−αcosαsin2α12ELsinα+α2−αcosαdαZ≥42,where i = 1, 2 are the main and driven gears, respectively, α is the pressure angle, G is the shear modulus of elasticity, and Ax = 2hxL is the gear cross-sectional area.Based on the distribution of the aforementioned gear pairs in a meshing period, MATLAB programming can be used to obtain the single-double-tooth meshing comprehensive stiffness of the spur gear pair in a period as shown in Figure5.Figure 5
Time-varying mesh stiffness of gears.


## 2.2.1. Calculation of Time-Varying Meshing Stiffness

In the process of gear meshing, single and double tooth pairs of gears alternate with the meshing periodically, which results in periodic changes in the comprehensive meshing stiffness of the gears, which causes dynamic excitation of the gear system. The load required to produce a deflection of 1μm on an average of 1 mm tooth width is called the meshing stiffness of the gear teeth [19]. Therefore, the establishment of a more practical stiffness model and the study of meshing stiffness are of great significance. At present, the meshing stiffness model is mainly obtained by the following two methods:(1)
Formula Methods. This method mainly includes Ishikawa formula and its improved algorithm, potential energy method and improved energy method, and Weber formula.(2)
Finite Element Method. The finite element method is the most widely used method currently. An empirical formula that can approximate the stiffness can be obtained through the calculated finite element analysis results.The potential energy method [20] uses the idea of variable increment infinite approximation to calculate the position of each meshing point in the gear meshing process, mapping each meshing point position to a single tooth profile, and using the material mechanics cantilever beam theory to calculate the time-varying meshing stiffness of the gear, which has high calculation efficiency. The result is accurate and includes other advantages. In this paper, the potential energy method is used to solve the time-varying meshing stiffness of gears. During the gear meshing process, the force on the gear teeth can be equivalent to the force on the cantilever beam. The equivalent force model of the spur gear is shown in Figure 4. In the figure, F is the normal load received by the spur gear, and α is the pressure angle. h is half the tooth thickness, where the load F is applied; d is the distance from the root circle to the addendum circle.Figure 4
Equivalent model of the spur gear.According to the knowledge of mechanics, the deformation of the gear teeth under the action of the normal loadF on the spur gear can be equivalently decomposed into bending deformation, shear deformation, and radial deformation along the tooth height direction. Based on the theory of material mechanics, the bending stiffness Kb, the shear stiffness Ks, the radial compression stiffness Kr, the stiffness Kd produced by the deformation of the gear teeth, and the stiffness Kh produced by the Hertz contact deformation can be obtained. Therefore, the meshing stiffness of single tooth pair of spur gear pair is expressed as(4)1Km=1∑i=121/Kbi+1/Ksi+1/Kri+1/Kdi+1/Khi,where(5)1kb=∫α1−α231+cosα1α2−αsinα−cosα2α2−αcosα2ELsinα+α2−αcosα3dα+∫0Rb−Rdd+x1cosα1−hsinαEIx1dx1Z<42∫α1−α331+cosα1α2−αsinα−cosα2α2−αcosα2ELsinα+α2−αcosα3dαZ≥42,1ks=∫−α1α21.21+μα2−αcosαcos2α1ELsinα+α2−αcosαdα+∫0Rb−Rd1.2cosα12GAx1dx1Z<42∫−α1α31.21+μα2−αcosαcos2α1ELsinα+α2−αcosαdαZ≥42,1kd=cos2αELLμfSf2+MμfSf+P1+Qtan2α,1kr=∫−α1α2α2−αcosαsin2α12ELsinα+α2−αcosαdα+∫0Rb−Rdsinα12EAx1dx1Z<42∫−α1α3α2−αcosαsin2α12ELsinα+α2−αcosαdαZ≥42,where i = 1, 2 are the main and driven gears, respectively, α is the pressure angle, G is the shear modulus of elasticity, and Ax = 2hxL is the gear cross-sectional area.Based on the distribution of the aforementioned gear pairs in a meshing period, MATLAB programming can be used to obtain the single-double-tooth meshing comprehensive stiffness of the spur gear pair in a period as shown in Figure5.Figure 5
Time-varying mesh stiffness of gears.


## 2.3. Mechanical Resonance Analysis

Servo drive system has many vibration characteristics during operation. This paper mainly discusses the characteristic vibration of torsional vibration generated by the electromechanical system.The mechanical rotor is an elastic body, so in the process of rotation, various rotating parts will form different degrees of torsional vibration. Torsional vibration may originate from the instability of the motor output or load, or it may originate from the internal excitation of the transmission system. These external and internal disturbances lead to forced vibration or free vibration of the mechanical system.During the operation of the servo drive system, the drive motor directly drives the input shaft to rotate, and the shaft and gear pair are transmitted to the load. During the rotation process, each component will produce a certain degree of elastic deformation, forming an angular deviation, which makes the mechanical transmission system torsionally vibrate. This torsional vibration is mainly reflected in the fluctuation of angle, angular velocity, and angular acceleration.


### 2.3.1. Mechanical Resonance Analysis with Fixed Stiffness

With analysis of time-varying meshing stiffness, the average meshing stiffness of the gear is 3.8 × 108 N/m. Therefore, the meshing stiffness of the gear pair is equivalent to a fixed stiffness of 3.8 × 108 N/m. The Bode diagram shown in Figure 6 of the four-mass servo system can be obtained using Equation (1).Figure 6
Bode diagram of the system with fixed stiffness.The four-mass servo system has three resonant frequency points, which are the input shaft, output shaft, and gear pair in the system, and all introduce the corresponding mechanical resonance to the system. The transfer function between the system motor speed and electromagnetic torque can be expressed as(6)ωesTes=1Jes∏n=13s2+2ξanωan+ωan2s2+2ξbnωbn+ωbn2,where ωbn and ωan are the resonant frequency and antiresonant frequency of the nth resonant link, respectively; and ξbn and ξan are the resonant damping coefficient and antiresonant damping coefficient of the nth resonant link, respectively.The transfer function of four-mass system that brings mechanical resonance is given by Equation (6):(7)Gs=∏n=13s2+2ξanωan+ωan2s2+2ξbnωbn+ωbn2.The analysis of Equation (6) shows that the resonance frequency of the system is related to the stiffness coefficient of the transmission mechanism and the equivalent inertia of each part of the system.In order to investigate the speed control performance of the load, assume that the meshing fixed stiffnesses are 3 × 107, 3 × 108, and 3 × 109 N/m, and the rest of the system parameters remain unchanged. The Bode diagram of the input to the output with the four-mass servo system can be obtained as shown in Figure 7. It can be seen that the change of the fixed stiffness value will change the resonant frequency of the system; that is, the resonant frequency value of the system will increase with the increase of the fixed stiffness.Figure 7
Bode diagram of the servo system with different stiffnesses.


### 2.3.2. Analysis of Mechanical Resonance with Time-Varying Meshing Stiffness

In order to further investigate the influence of the time-varying meshing stiffness on the mechanical resonance of the four-mass servo system, the change range of the resonance frequency within the range of the time-varying meshing stiffness is analyzed. Here, the different stiffness values (2.7 × 108,3.8 × 108 and 4.9 × 108 N/m) during the gear meshing process are taken. The Bode diagram of the input electromagnetic torque to the output load speed can be obtained as shown in Figure 8. It can be seen that the mechanical resonance frequency range brought about by the time-varying meshing stiffness is from 5966 Hz to 7638 Hz.Figure 8
Bode diagram of electromagnetic torque to load speed.According to Equation (1), a four-mass servo system model considering the time-varying meshing stiffness of gear is established. The servo system model established is compared and analyzed under the conditions of time-varying meshing stiffness and constant stiffness. When a speed command of 1000 r/min is given, the time domain characteristics of the torque transmitted at the load are obtained in Figure 9, and when considering the time-varying stiffness parameters, the system load transfer torque will produce greater fluctuations than considering the constant stiffness parameters. From the resonance analysis of the system under constant stiffness and variable stiffness parameters, it can be concluded that it is necessary to consider the gear meshing stiffness as a time-varying meshing stiffness when we research four-mass servo system.Figure 9
Transmit torqueTs2 change of the servo system under constant stiffness and variable stiffness.


## 2.3.1. Mechanical Resonance Analysis with Fixed Stiffness

With analysis of time-varying meshing stiffness, the average meshing stiffness of the gear is 3.8 × 108 N/m. Therefore, the meshing stiffness of the gear pair is equivalent to a fixed stiffness of 3.8 × 108 N/m. The Bode diagram shown in Figure 6 of the four-mass servo system can be obtained using Equation (1).Figure 6
Bode diagram of the system with fixed stiffness.The four-mass servo system has three resonant frequency points, which are the input shaft, output shaft, and gear pair in the system, and all introduce the corresponding mechanical resonance to the system. The transfer function between the system motor speed and electromagnetic torque can be expressed as(6)ωesTes=1Jes∏n=13s2+2ξanωan+ωan2s2+2ξbnωbn+ωbn2,where ωbn and ωan are the resonant frequency and antiresonant frequency of the nth resonant link, respectively; and ξbn and ξan are the resonant damping coefficient and antiresonant damping coefficient of the nth resonant link, respectively.The transfer function of four-mass system that brings mechanical resonance is given by Equation (6):(7)Gs=∏n=13s2+2ξanωan+ωan2s2+2ξbnωbn+ωbn2.The analysis of Equation (6) shows that the resonance frequency of the system is related to the stiffness coefficient of the transmission mechanism and the equivalent inertia of each part of the system.In order to investigate the speed control performance of the load, assume that the meshing fixed stiffnesses are 3 × 107, 3 × 108, and 3 × 109 N/m, and the rest of the system parameters remain unchanged. The Bode diagram of the input to the output with the four-mass servo system can be obtained as shown in Figure 7. It can be seen that the change of the fixed stiffness value will change the resonant frequency of the system; that is, the resonant frequency value of the system will increase with the increase of the fixed stiffness.Figure 7
Bode diagram of the servo system with different stiffnesses.


## 2.3.2. Analysis of Mechanical Resonance with Time-Varying Meshing Stiffness

In order to further investigate the influence of the time-varying meshing stiffness on the mechanical resonance of the four-mass servo system, the change range of the resonance frequency within the range of the time-varying meshing stiffness is analyzed. Here, the different stiffness values (2.7 × 108,3.8 × 108 and 4.9 × 108 N/m) during the gear meshing process are taken. The Bode diagram of the input electromagnetic torque to the output load speed can be obtained as shown in Figure 8. It can be seen that the mechanical resonance frequency range brought about by the time-varying meshing stiffness is from 5966 Hz to 7638 Hz.Figure 8
Bode diagram of electromagnetic torque to load speed.According to Equation (1), a four-mass servo system model considering the time-varying meshing stiffness of gear is established. The servo system model established is compared and analyzed under the conditions of time-varying meshing stiffness and constant stiffness. When a speed command of 1000 r/min is given, the time domain characteristics of the torque transmitted at the load are obtained in Figure 9, and when considering the time-varying stiffness parameters, the system load transfer torque will produce greater fluctuations than considering the constant stiffness parameters. From the resonance analysis of the system under constant stiffness and variable stiffness parameters, it can be concluded that it is necessary to consider the gear meshing stiffness as a time-varying meshing stiffness when we research four-mass servo system.Figure 9
Transmit torqueTs2 change of the servo system under constant stiffness and variable stiffness.


## 3. Design of the Composite Controller

In the servo drive system, the drive motor, gear transmission, and the load are connected by a transmission shaft. The transmission shaft has certain elasticity. When the transmission shaft is torsionally deformed, a certain torque will be generated. This torque can be regarded as load torque or driving torque [21]. The electromagnetic torque Te provided by the motor and the torque when the transmission shaft is twisted and deformed together act on the motor shaft. Most transmission devices can be regarded as elastic elements, and all have their own resonant frequency. The existence of elastic elements will bring in one or even several resonance points in the system, which will aggravate the occurrence of mechanical resonance. The structure diagram of the composite controller system is shown in Figure 10.Figure 10
Block diagram of the control system.


### 3.1. Model Predictive Control

Most of the control methods of servo drive systems are based on cascade structure. This structure has multiple PI regulators with cascaded control loops, which are used for three-loop control of current loop, speed loop, and position loop. The cascaded vector control structure is relatively simple and widely used in the industry, but it also has its shortcomings.The first problem is that the vector control structure of the classic cascaded permanent magnet synchronous motor is complicated, and multiple control loops greatly limit the dynamic response of the system, and the dynamic response of the small inertia system cannot be improved; the other problem is the control of a separate control loop, the controller, runs independently, which means that the overall control performance optimization and effective constraint processing are almost impossible. When the microcontroller adjusts the complex control target, the controller cannot meet all the conditions.Model Predictive Control is, as shown in Figure11, a form of control, in which the control action is obtained by solving online, at each sampling instant, a finite horizon optimal control problem, in which the initial state is the current state of the plant. Optimization yields a finite control sequence, and the first control action in this sequence is applied to the plant. MPC differs from conventional control, in which the control law is precomputed of line [16]. Compared with PID and Fuzzy control algorithms, Model Predictive Control can handle complex coupling problems such as multi-inputs, outputs, and multivariable coupling [22].Figure 11
Structure diagram of Model Predictive Control.


#### 3.1.1. Establishment of the Prediction Model

Considering that the gear meshing stiffness is a time-varying meshing stiffness, the continuous system state space equation of the four-mass servo system can be obtained using Equation (1):(8)x˙t=Axt+But+Cdtyt=Ext,where Select x=θ˙eθ˙1θ˙2θ˙Lθeθ1θ2θLT is the state variable of system; u=Te are the control inputs; and d=TL is a measurable amount of external interference, and then the coefficient matrix is as follows:(9)A=−Je−1C1J1−1C1001000Je−1C1A22A32001000A23A33JL−1C2001000J2−1C2−JL−1C20001−Je−1K1J1−1K1000000Je−1K1A26A360000000A270A37J2−1K2JL−1K2−JL−1K200000000,A22=−J1−1C1+Cm⋅R12;A23=J1−1⋅Cm⋅R1⋅R2,A26=−J1−1K1+Km⋅R12;A27=J1−1⋅Km⋅R1⋅R2,A32=J2−1⋅Cm⋅R1⋅R2;A33=−J2−1C2+Cm⋅R22,A36=J2−1⋅Km⋅R1⋅R2;A37=−J2−1K2+Km⋅R22,B=Je−10000000T,C=000−JL−10000T,E=00010000.The forward Euler method [23] is used to discretize Equation (8), and the state space equation of the discrete system is obtained as follows:(10)xk+1=Akxk+Bkuk+Ckdkyk=Ekxk,where k is the current sampling time and k + 1 is the next sampling time.


#### 3.1.2. Design of Objective Function

Letηk=xkuk−1dk−1T ; then Equation (10) can be rewritten as follows:(11)ηk+1=A⌢kηk+B⌢kΔuk+C⌢kΔdkγk=E⌢kηk,where A⌢k=AkBkCk0m⋅nIm0001, B⌢k=BkIm0, C⌢k=Ck01T, E⌢k=Ek00, m is the dimension of control quantity, and n is the dimension of state quantity.Set the predictive horizon asNp, and control horizon as Np, Nc<Np. According to Equation (11), the outputs of future prediction of the system can be deduced as follows:(12)Yk=Φηk+ΨΔU+ΓΔD,where Yk=γk+1γk+2⋯γk+NpT; Φ=E⌢kA⌢kE⌢kA⌢k2⋯E⌢kA⌢kNpT; and(13)Ψ=E⌢kB⌢k00⋯0E⌢kA⌢kB⌢kE⌢kB⌢k0⋯0⋮E⌢kA⌢kNc−1B⌢k⋮E⌢kA⌢kNc−2B⌢k⋮⋯⋱⋯⋮E⌢kB⌢k⋮⋮⋮⋱⋮E⌢kA⌢kNP−1B⌢kE⌢kA⌢kNP−2B⌢kE⌢kA⌢kNP−3B⌢k⋯E⌢kA⌢kNP−Nc−1B⌢k,ΔU=ΔukΔuk+1⋯Δuk+NcT,Γ=E⌢kC⌢k00⋯0E⌢kA⌢kC⌢kE⌢kC⌢k0⋯0⋮E⌢kA⌢kNc−1C⌢k⋮E⌢kA⌢kNc−2C⌢k⋮⋯⋱⋯⋮E⌢kC⌢k⋮⋮⋮⋱⋮E⌢kA⌢kNP−1C⌢kE⌢kA⌢kNP−2C⌢kE⌢kA⌢kNP−3C⌢k⋯E⌢kA⌢kNP−Nc−1C⌢k,ΔD=ΔdkΔdk+1⋯Δdk+NcT,The system control objective is to achieve accurate tracking of load speed without exceeding the maximum allowable electromagnetic torque. The optimization objective function is(14)Jk=∑i=1NpYk+i|k−Yrefk+i|k2Q+∑i=0Nc−1ΔUk+i2R,where Yrefk+i|k is the reference value of the system output variable, k+i|k represents the system output at time k to predict the value at time k + 1,where i=1,2,…,NpΔUk+i is the time k+i system control input increment, namely, the electromagnetic torque increment, where i=0,1,…,Nc−1; Q and R are weight matrices.


#### 3.1.3. Design of Constraint Conditions

Constraints are imposed on the electromagnetic torque, electromagnetic torque increment, and system output of the control input. The constraints of the control input and its increment are(15)umin≤uk+i≤umaxΔumin≤Δuk+i≤Δumax,where umin and umax are the minimum and maximum values of the electromagnetic torque, respectively, and Δumin and Δumax are the minimum and maximum values of the electromagnetic torque increment, respectively.The output constraint condition of system is(16)γmin≤γk+i≤γmax,where γmin and γmax are the minimum value and maximum value of system output, respectively.


#### 3.1.4. Optimization Solution

The system optimization problem is that, in the prediction time domain, Equation (14) reaches the minimum under the constraints of Equations (15) and (16). That is, the following problems need to be solved in each control cycle:(17)minΔUJηk,ΔUks.t.Umin≤SΔU+Ut≤ΔUmaxΔUmin≤ΔU≤ΔUmaxγmin≤γk+i≤γmax,where Umin and Umax are the minimum value set and the maximum value set of the control quantity, respectively; ΔUmin and ΔUmax are the minimum value set and the maximum value set of control increment, respectively; and S=K⊗Im; Ut=1Nc⊗uk−1, where 1Nc is the column vector of NC row, uk−1 is the actual control quantity at time k-1, ⊗ is the Kronecker product, and(18)K=10⋯⋯0110⋯0111⋱0⋮⋮⋱⋱011⋯11Nc×Nc.By solving the quadratic programming problem, the increment of control input in the control time domain can be obtained as(19)ΔU=ΔukΔuk+1⋯Δuk+Nc−1T.Take out the first itemΔuk of the control input increment, and the system control input at the current time can be obtained as(20)uk=uk−1+Δuk.


### 3.2. Notch Filter

There are many uncertainties in the actual operation process, so there will be some errors between the predicted output estimated by the prediction model and the actual output of the system. Due to the nonlinear factor of the time-varying meshing stiffness of the gear pair, the speed deviation (the absolute value of the difference between the actual speed and the given speed) will be large, and the corresponding gain will increase accordingly. At this time, the extracted rotational speed error signal should be processed to reduce the influence of rotational speed fluctuation on the system. Some of Filters can do this well.Unlike low-pass filters and high-pass filters, the Notch Filter is a circuit that allows only specific frequencies to pass, while effectively suppressing signals at other frequencies. Because of its selectivity to signals, it is widely used in electronic design. It can be divided into two-parameter Notch Filters and three-parameter Notch Filters. Compared with Low-pass Filters and biquad Filters, there are fewer parameters to be determined, and the notch center frequency can be accurately set by Fast Fourier Transform (FFT) to detect the resonant frequency.The design of the two-parameter Notch Filter is simple, but it only reflects the notch frequency and the width of the notch and cannot control the depth of the notch, so it is easy to attenuate the gain of the system at the notch frequency to no response [24]. In order to improve the performance of the Notch Filter and make up for the deficiencies of the two-parameter Notch Filter, this paper adopts a relatively complete three-parameter Notch Filter.The transfer function of the three-parameter Notch Filter is shown in(21)Gs=s2+2zωns+ωn2s2+2pωns+ωn2.where ωn is the center frequency of the notch and z and p are the damping ratio coefficients of the notch, respectively.Figure12(a) shows the influence of different parameters on the Notch Filter. It can be seen from Figure 12(b) that different systems can choose different notch center frequencies, and unfitting notch depth parameters will also cause frequency signals near the notch center frequency to be filtered out, resulting in system instability. The attenuation amplitude at the notch center frequency ωn is 20lgp/z dB, and the bandwidth is zωn. It can be concluded that, compared with the two-parameter Notch Filter, the amplitude attenuation of the three-parameter Notch Filter can be controlled, which can better suppress the mechanical resonance at the resonance point.Figure 12
Parametric characteristics of three-parameter Notch Filter. (a) Different notch depths. (b)Different notch center frequencies.
(a)(b)


### 3.3. Resonance Suppression Control

By combining the above analysis, mechanical resonance suppression structure block diagram of the four-mass servo drive system is shown in Figure13. Aiming at the mechanical resonance suppression of time-varying meshing stiffness in this paper, a three-parameter Notch Filter is connected in series to speed feedback loop. Extract the resonant frequency of four-mass servo drive system as the notch center frequency of the Notch Filter by using FFT. The speed will generate the larger gain value at the frequency design, and the difference from the reference value will become a negative gain. The positive gain of transmission mechanisms at the resonant frequency is cancelled, thereby achieving suppression of mechanical resonance. According to Bode diagram of the system, there are three resonance frequency points in the system. The results of FFT analysis of the speed feedback are shown in Figure 14. The resonance frequency near 640 Hz corresponds to the larger positive gain value, which has the greatest impact on the performance of four-mass servo drive system. The results show that it is consistent with Bode diagram analysis. Therefore, it is the basis of the mechanical resonance analysis above, and the mechanical resonance at this resonant frequency point can be suppressed.Figure 13
Block diagram of the servo drive system resonant suppression.Figure 14
Frequency domain analysis of load feedback speed.


## 3.1. Model Predictive Control

Most of the control methods of servo drive systems are based on cascade structure. This structure has multiple PI regulators with cascaded control loops, which are used for three-loop control of current loop, speed loop, and position loop. The cascaded vector control structure is relatively simple and widely used in the industry, but it also has its shortcomings.The first problem is that the vector control structure of the classic cascaded permanent magnet synchronous motor is complicated, and multiple control loops greatly limit the dynamic response of the system, and the dynamic response of the small inertia system cannot be improved; the other problem is the control of a separate control loop, the controller, runs independently, which means that the overall control performance optimization and effective constraint processing are almost impossible. When the microcontroller adjusts the complex control target, the controller cannot meet all the conditions.Model Predictive Control is, as shown in Figure11, a form of control, in which the control action is obtained by solving online, at each sampling instant, a finite horizon optimal control problem, in which the initial state is the current state of the plant. Optimization yields a finite control sequence, and the first control action in this sequence is applied to the plant. MPC differs from conventional control, in which the control law is precomputed of line [16]. Compared with PID and Fuzzy control algorithms, Model Predictive Control can handle complex coupling problems such as multi-inputs, outputs, and multivariable coupling [22].Figure 11
Structure diagram of Model Predictive Control.


### 3.1.1. Establishment of the Prediction Model

Considering that the gear meshing stiffness is a time-varying meshing stiffness, the continuous system state space equation of the four-mass servo system can be obtained using Equation (1):(8)x˙t=Axt+But+Cdtyt=Ext,where Select x=θ˙eθ˙1θ˙2θ˙Lθeθ1θ2θLT is the state variable of system; u=Te are the control inputs; and d=TL is a measurable amount of external interference, and then the coefficient matrix is as follows:(9)A=−Je−1C1J1−1C1001000Je−1C1A22A32001000A23A33JL−1C2001000J2−1C2−JL−1C20001−Je−1K1J1−1K1000000Je−1K1A26A360000000A270A37J2−1K2JL−1K2−JL−1K200000000,A22=−J1−1C1+Cm⋅R12;A23=J1−1⋅Cm⋅R1⋅R2,A26=−J1−1K1+Km⋅R12;A27=J1−1⋅Km⋅R1⋅R2,A32=J2−1⋅Cm⋅R1⋅R2;A33=−J2−1C2+Cm⋅R22,A36=J2−1⋅Km⋅R1⋅R2;A37=−J2−1K2+Km⋅R22,B=Je−10000000T,C=000−JL−10000T,E=00010000.The forward Euler method [23] is used to discretize Equation (8), and the state space equation of the discrete system is obtained as follows:(10)xk+1=Akxk+Bkuk+Ckdkyk=Ekxk,where k is the current sampling time and k + 1 is the next sampling time.


### 3.1.2. Design of Objective Function

Letηk=xkuk−1dk−1T ; then Equation (10) can be rewritten as follows:(11)ηk+1=A⌢kηk+B⌢kΔuk+C⌢kΔdkγk=E⌢kηk,where A⌢k=AkBkCk0m⋅nIm0001, B⌢k=BkIm0, C⌢k=Ck01T, E⌢k=Ek00, m is the dimension of control quantity, and n is the dimension of state quantity.Set the predictive horizon asNp, and control horizon as Np, Nc<Np. According to Equation (11), the outputs of future prediction of the system can be deduced as follows:(12)Yk=Φηk+ΨΔU+ΓΔD,where Yk=γk+1γk+2⋯γk+NpT; Φ=E⌢kA⌢kE⌢kA⌢k2⋯E⌢kA⌢kNpT; and(13)Ψ=E⌢kB⌢k00⋯0E⌢kA⌢kB⌢kE⌢kB⌢k0⋯0⋮E⌢kA⌢kNc−1B⌢k⋮E⌢kA⌢kNc−2B⌢k⋮⋯⋱⋯⋮E⌢kB⌢k⋮⋮⋮⋱⋮E⌢kA⌢kNP−1B⌢kE⌢kA⌢kNP−2B⌢kE⌢kA⌢kNP−3B⌢k⋯E⌢kA⌢kNP−Nc−1B⌢k,ΔU=ΔukΔuk+1⋯Δuk+NcT,Γ=E⌢kC⌢k00⋯0E⌢kA⌢kC⌢kE⌢kC⌢k0⋯0⋮E⌢kA⌢kNc−1C⌢k⋮E⌢kA⌢kNc−2C⌢k⋮⋯⋱⋯⋮E⌢kC⌢k⋮⋮⋮⋱⋮E⌢kA⌢kNP−1C⌢kE⌢kA⌢kNP−2C⌢kE⌢kA⌢kNP−3C⌢k⋯E⌢kA⌢kNP−Nc−1C⌢k,ΔD=ΔdkΔdk+1⋯Δdk+NcT,The system control objective is to achieve accurate tracking of load speed without exceeding the maximum allowable electromagnetic torque. The optimization objective function is(14)Jk=∑i=1NpYk+i|k−Yrefk+i|k2Q+∑i=0Nc−1ΔUk+i2R,where Yrefk+i|k is the reference value of the system output variable, k+i|k represents the system output at time k to predict the value at time k + 1,where i=1,2,…,NpΔUk+i is the time k+i system control input increment, namely, the electromagnetic torque increment, where i=0,1,…,Nc−1; Q and R are weight matrices.


### 3.1.3. Design of Constraint Conditions

Constraints are imposed on the electromagnetic torque, electromagnetic torque increment, and system output of the control input. The constraints of the control input and its increment are(15)umin≤uk+i≤umaxΔumin≤Δuk+i≤Δumax,where umin and umax are the minimum and maximum values of the electromagnetic torque, respectively, and Δumin and Δumax are the minimum and maximum values of the electromagnetic torque increment, respectively.The output constraint condition of system is(16)γmin≤γk+i≤γmax,where γmin and γmax are the minimum value and maximum value of system output, respectively.


### 3.1.4. Optimization Solution

The system optimization problem is that, in the prediction time domain, Equation (14) reaches the minimum under the constraints of Equations (15) and (16). That is, the following problems need to be solved in each control cycle:(17)minΔUJηk,ΔUks.t.Umin≤SΔU+Ut≤ΔUmaxΔUmin≤ΔU≤ΔUmaxγmin≤γk+i≤γmax,where Umin and Umax are the minimum value set and the maximum value set of the control quantity, respectively; ΔUmin and ΔUmax are the minimum value set and the maximum value set of control increment, respectively; and S=K⊗Im; Ut=1Nc⊗uk−1, where 1Nc is the column vector of NC row, uk−1 is the actual control quantity at time k-1, ⊗ is the Kronecker product, and(18)K=10⋯⋯0110⋯0111⋱0⋮⋮⋱⋱011⋯11Nc×Nc.By solving the quadratic programming problem, the increment of control input in the control time domain can be obtained as(19)ΔU=ΔukΔuk+1⋯Δuk+Nc−1T.Take out the first itemΔuk of the control input increment, and the system control input at the current time can be obtained as(20)uk=uk−1+Δuk.


## 3.1.1. Establishment of the Prediction Model

Considering that the gear meshing stiffness is a time-varying meshing stiffness, the continuous system state space equation of the four-mass servo system can be obtained using Equation (1):(8)x˙t=Axt+But+Cdtyt=Ext,where Select x=θ˙eθ˙1θ˙2θ˙Lθeθ1θ2θLT is the state variable of system; u=Te are the control inputs; and d=TL is a measurable amount of external interference, and then the coefficient matrix is as follows:(9)A=−Je−1C1J1−1C1001000Je−1C1A22A32001000A23A33JL−1C2001000J2−1C2−JL−1C20001−Je−1K1J1−1K1000000Je−1K1A26A360000000A270A37J2−1K2JL−1K2−JL−1K200000000,A22=−J1−1C1+Cm⋅R12;A23=J1−1⋅Cm⋅R1⋅R2,A26=−J1−1K1+Km⋅R12;A27=J1−1⋅Km⋅R1⋅R2,A32=J2−1⋅Cm⋅R1⋅R2;A33=−J2−1C2+Cm⋅R22,A36=J2−1⋅Km⋅R1⋅R2;A37=−J2−1K2+Km⋅R22,B=Je−10000000T,C=000−JL−10000T,E=00010000.The forward Euler method [23] is used to discretize Equation (8), and the state space equation of the discrete system is obtained as follows:(10)xk+1=Akxk+Bkuk+Ckdkyk=Ekxk,where k is the current sampling time and k + 1 is the next sampling time.


## 3.1.2. Design of Objective Function

Letηk=xkuk−1dk−1T ; then Equation (10) can be rewritten as follows:(11)ηk+1=A⌢kηk+B⌢kΔuk+C⌢kΔdkγk=E⌢kηk,where A⌢k=AkBkCk0m⋅nIm0001, B⌢k=BkIm0, C⌢k=Ck01T, E⌢k=Ek00, m is the dimension of control quantity, and n is the dimension of state quantity.Set the predictive horizon asNp, and control horizon as Np, Nc<Np. According to Equation (11), the outputs of future prediction of the system can be deduced as follows:(12)Yk=Φηk+ΨΔU+ΓΔD,where Yk=γk+1γk+2⋯γk+NpT; Φ=E⌢kA⌢kE⌢kA⌢k2⋯E⌢kA⌢kNpT; and(13)Ψ=E⌢kB⌢k00⋯0E⌢kA⌢kB⌢kE⌢kB⌢k0⋯0⋮E⌢kA⌢kNc−1B⌢k⋮E⌢kA⌢kNc−2B⌢k⋮⋯⋱⋯⋮E⌢kB⌢k⋮⋮⋮⋱⋮E⌢kA⌢kNP−1B⌢kE⌢kA⌢kNP−2B⌢kE⌢kA⌢kNP−3B⌢k⋯E⌢kA⌢kNP−Nc−1B⌢k,ΔU=ΔukΔuk+1⋯Δuk+NcT,Γ=E⌢kC⌢k00⋯0E⌢kA⌢kC⌢kE⌢kC⌢k0⋯0⋮E⌢kA⌢kNc−1C⌢k⋮E⌢kA⌢kNc−2C⌢k⋮⋯⋱⋯⋮E⌢kC⌢k⋮⋮⋮⋱⋮E⌢kA⌢kNP−1C⌢kE⌢kA⌢kNP−2C⌢kE⌢kA⌢kNP−3C⌢k⋯E⌢kA⌢kNP−Nc−1C⌢k,ΔD=ΔdkΔdk+1⋯Δdk+NcT,The system control objective is to achieve accurate tracking of load speed without exceeding the maximum allowable electromagnetic torque. The optimization objective function is(14)Jk=∑i=1NpYk+i|k−Yrefk+i|k2Q+∑i=0Nc−1ΔUk+i2R,where Yrefk+i|k is the reference value of the system output variable, k+i|k represents the system output at time k to predict the value at time k + 1,where i=1,2,…,NpΔUk+i is the time k+i system control input increment, namely, the electromagnetic torque increment, where i=0,1,…,Nc−1; Q and R are weight matrices.


## 3.1.3. Design of Constraint Conditions

Constraints are imposed on the electromagnetic torque, electromagnetic torque increment, and system output of the control input. The constraints of the control input and its increment are(15)umin≤uk+i≤umaxΔumin≤Δuk+i≤Δumax,where umin and umax are the minimum and maximum values of the electromagnetic torque, respectively, and Δumin and Δumax are the minimum and maximum values of the electromagnetic torque increment, respectively.The output constraint condition of system is(16)γmin≤γk+i≤γmax,where γmin and γmax are the minimum value and maximum value of system output, respectively.


## 3.1.4. Optimization Solution

The system optimization problem is that, in the prediction time domain, Equation (14) reaches the minimum under the constraints of Equations (15) and (16). That is, the following problems need to be solved in each control cycle:(17)minΔUJηk,ΔUks.t.Umin≤SΔU+Ut≤ΔUmaxΔUmin≤ΔU≤ΔUmaxγmin≤γk+i≤γmax,where Umin and Umax are the minimum value set and the maximum value set of the control quantity, respectively; ΔUmin and ΔUmax are the minimum value set and the maximum value set of control increment, respectively; and S=K⊗Im; Ut=1Nc⊗uk−1, where 1Nc is the column vector of NC row, uk−1 is the actual control quantity at time k-1, ⊗ is the Kronecker product, and(18)K=10⋯⋯0110⋯0111⋱0⋮⋮⋱⋱011⋯11Nc×Nc.By solving the quadratic programming problem, the increment of control input in the control time domain can be obtained as(19)ΔU=ΔukΔuk+1⋯Δuk+Nc−1T.Take out the first itemΔuk of the control input increment, and the system control input at the current time can be obtained as(20)uk=uk−1+Δuk.


## 3.2. Notch Filter

There are many uncertainties in the actual operation process, so there will be some errors between the predicted output estimated by the prediction model and the actual output of the system. Due to the nonlinear factor of the time-varying meshing stiffness of the gear pair, the speed deviation (the absolute value of the difference between the actual speed and the given speed) will be large, and the corresponding gain will increase accordingly. At this time, the extracted rotational speed error signal should be processed to reduce the influence of rotational speed fluctuation on the system. Some of Filters can do this well.Unlike low-pass filters and high-pass filters, the Notch Filter is a circuit that allows only specific frequencies to pass, while effectively suppressing signals at other frequencies. Because of its selectivity to signals, it is widely used in electronic design. It can be divided into two-parameter Notch Filters and three-parameter Notch Filters. Compared with Low-pass Filters and biquad Filters, there are fewer parameters to be determined, and the notch center frequency can be accurately set by Fast Fourier Transform (FFT) to detect the resonant frequency.The design of the two-parameter Notch Filter is simple, but it only reflects the notch frequency and the width of the notch and cannot control the depth of the notch, so it is easy to attenuate the gain of the system at the notch frequency to no response [24]. In order to improve the performance of the Notch Filter and make up for the deficiencies of the two-parameter Notch Filter, this paper adopts a relatively complete three-parameter Notch Filter.The transfer function of the three-parameter Notch Filter is shown in(21)Gs=s2+2zωns+ωn2s2+2pωns+ωn2.where ωn is the center frequency of the notch and z and p are the damping ratio coefficients of the notch, respectively.Figure12(a) shows the influence of different parameters on the Notch Filter. It can be seen from Figure 12(b) that different systems can choose different notch center frequencies, and unfitting notch depth parameters will also cause frequency signals near the notch center frequency to be filtered out, resulting in system instability. The attenuation amplitude at the notch center frequency ωn is 20lgp/z dB, and the bandwidth is zωn. It can be concluded that, compared with the two-parameter Notch Filter, the amplitude attenuation of the three-parameter Notch Filter can be controlled, which can better suppress the mechanical resonance at the resonance point.Figure 12
Parametric characteristics of three-parameter Notch Filter. (a) Different notch depths. (b)Different notch center frequencies.
(a)(b)


## 3.3. Resonance Suppression Control

By combining the above analysis, mechanical resonance suppression structure block diagram of the four-mass servo drive system is shown in Figure13. Aiming at the mechanical resonance suppression of time-varying meshing stiffness in this paper, a three-parameter Notch Filter is connected in series to speed feedback loop. Extract the resonant frequency of four-mass servo drive system as the notch center frequency of the Notch Filter by using FFT. The speed will generate the larger gain value at the frequency design, and the difference from the reference value will become a negative gain. The positive gain of transmission mechanisms at the resonant frequency is cancelled, thereby achieving suppression of mechanical resonance. According to Bode diagram of the system, there are three resonance frequency points in the system. The results of FFT analysis of the speed feedback are shown in Figure 14. The resonance frequency near 640 Hz corresponds to the larger positive gain value, which has the greatest impact on the performance of four-mass servo drive system. The results show that it is consistent with Bode diagram analysis. Therefore, it is the basis of the mechanical resonance analysis above, and the mechanical resonance at this resonant frequency point can be suppressed.Figure 13
Block diagram of the servo drive system resonant suppression.Figure 14
Frequency domain analysis of load feedback speed.


## 4. Simulation Results and Discussion

Compared with the control effect based on MPC, the mechanical resonance suppression control performance of MPC + Notch Filter composite controller for four-mass servo drive system considering time-varying meshing stiffness is proved. The simulation model of four-mass servo drive system and resonance suppression controller are established in Matlab/Simulink. The main parameters of MPC controller are predictive horizon = 25, control horizon = 5, weighting parameterQ = 0.1, R = 1, and the electromagnetic torque range of the system control input is ±60 Nm. The notch center frequency of the Notch Filter is 637.9 Hz, the notch depth is 20 dB, and notch width is 0.5. The Bode diagram after resonance suppression of system is shown in Figure 15.Figure 15
Bode diagram of the system after resonance suppression.The parameters of servo drive system used in the simulation in this paper are shown in Table1.Table 1
Parameters of the servo system.
Parameter namesValuesRotational inertia of motor/loadJe/JL(kg·m2)6.15 × 10−3/5.75 × 10−3Rotational inertia of gearJ1/J2 (kg·m2)2.7 × 10−4/2.7 × 10−3Damping coefficient of shaftC1/C21Module of gearm (mm)3Pressure angle of gearα (°)20Tooth numbers of gearZ1/Z220/40Overlap ratio of gearε1.6Reference center distance of geara (mm)90Tooth widthW (mm)20Yang 's modulusE (PA)1.9 × 1011Poisson ratioυ0.33Torsional stiffnessK1/K2 (N·m/rad)105/105


### 4.1. Response Characteristic Analysis of Step Speed

The reference speeds of 0–200r/min and 200–1000r/min were selected for simulation experiments to simulate the speed response of servo drive system at different working conditions. The speed response curves of the two control methods are shown in Figure16. As shown in the Figure 15(a) and 15(b), with different reference speeds, the MPC + Notch Filter composite controller has small overshoot and great speed following performance. When using MPC method, the oscillation of transfer torque Ts2 is greater than that of MPC + Notch Filter composite controller. Analysis of the control effect of two control methods is given in Table 2. The speed response time of MPC method is faster than MPC + Notch Filter composite controller, because increasing the filter is equivalent to adding an additional delay link, resulting in a longer response time of the system. The servo drive system requires high accuracy of speed tremble and position offset, but the response time does not affect the control effect of four-mass servo drive system in this paper.Figure 16
Dynamic characteristic curve of step speed. (a) Speed response curve of 0–200r/min. (b) Speed response curve of 200–1000r/min. (c) Transfer torque response of load under 200–1000r/min.
(a)(b)(c)Table 2
Comparison of speed response characteristics.
Control methodMPCMPC + Notch FilterMPCMPC + Notch FilterSpeed0～200r/min200～1000r/minOvershoot11.6r/min1.1r/min17r/min3r/minResponse time0.04s0.06s0.15s0.18s


### 4.2. Dynamic Characteristic Analysis of Impact Load

The operation process of servo drive system is often affected by the impact load. Therefore, it is necessary to analyze the dynamic characteristics of the two control methods with impact load. The reference speed is still set to 1000 r/min. The impact load is applied at 0.5 s, and torque increases from 10 Nm to 30 Nm and then recovers to 10 Nm, lasting for 0.2 s. The process is shown in Figure17.Figure 17
Schematic diagram of the impact load.The simulation results are shown in Figure18. With the impact load, in order to balance the load, electromagnetic torque suddenly increased and then decreased. The whole process showed two oscillation attenuation phenomena. Table 3 lists some contrast effects with two control methods. The speed response curves of four-mass servo drive system controlled by Model Predictive Control with Notch Filter composite controller also have the advantages of small overshoot and great robustness.Figure 18
Dynamic response characteristics with impact load. (a) Curve of speed response. (b) Curve of transmission torqueTs2.
(a)(b)Table 3
Dynamic characteristics of impact load.
Control methodMPCMPC + Notch FilterMPCMPC + Notch FilterTorque10～30 Nm30～10 NmOvershoot−34.3 r/min−12.8 r/min+32 r/min+10.2 r/minTorque ripple0.65 Nm0.12 N/m0.68 Nm0.09 Nm


## 4.1. Response Characteristic Analysis of Step Speed

The reference speeds of 0–200r/min and 200–1000r/min were selected for simulation experiments to simulate the speed response of servo drive system at different working conditions. The speed response curves of the two control methods are shown in Figure16. As shown in the Figure 15(a) and 15(b), with different reference speeds, the MPC + Notch Filter composite controller has small overshoot and great speed following performance. When using MPC method, the oscillation of transfer torque Ts2 is greater than that of MPC + Notch Filter composite controller. Analysis of the control effect of two control methods is given in Table 2. The speed response time of MPC method is faster than MPC + Notch Filter composite controller, because increasing the filter is equivalent to adding an additional delay link, resulting in a longer response time of the system. The servo drive system requires high accuracy of speed tremble and position offset, but the response time does not affect the control effect of four-mass servo drive system in this paper.Figure 16
Dynamic characteristic curve of step speed. (a) Speed response curve of 0–200r/min. (b) Speed response curve of 200–1000r/min. (c) Transfer torque response of load under 200–1000r/min.
(a)(b)(c)Table 2
Comparison of speed response characteristics.
Control methodMPCMPC + Notch FilterMPCMPC + Notch FilterSpeed0～200r/min200～1000r/minOvershoot11.6r/min1.1r/min17r/min3r/minResponse time0.04s0.06s0.15s0.18s


## 4.2. Dynamic Characteristic Analysis of Impact Load

The operation process of servo drive system is often affected by the impact load. Therefore, it is necessary to analyze the dynamic characteristics of the two control methods with impact load. The reference speed is still set to 1000 r/min. The impact load is applied at 0.5 s, and torque increases from 10 Nm to 30 Nm and then recovers to 10 Nm, lasting for 0.2 s. The process is shown in Figure17.Figure 17
Schematic diagram of the impact load.The simulation results are shown in Figure18. With the impact load, in order to balance the load, electromagnetic torque suddenly increased and then decreased. The whole process showed two oscillation attenuation phenomena. Table 3 lists some contrast effects with two control methods. The speed response curves of four-mass servo drive system controlled by Model Predictive Control with Notch Filter composite controller also have the advantages of small overshoot and great robustness.Figure 18
Dynamic response characteristics with impact load. (a) Curve of speed response. (b) Curve of transmission torqueTs2.
(a)(b)Table 3
Dynamic characteristics of impact load.
Control methodMPCMPC + Notch FilterMPCMPC + Notch FilterTorque10～30 Nm30～10 NmOvershoot−34.3 r/min−12.8 r/min+32 r/min+10.2 r/minTorque ripple0.65 Nm0.12 N/m0.68 Nm0.09 Nm


## 5. Conclusions

In this paper, a four-mass servo drive system model considering time-varying meshing stiffness of the gear is proposed. The time-varying meshing stiffness will cause a large fluctuation of the transmission torque of the system and bring adverse effects on the dynamic characteristics of system. A mechanical resonance suppression control strategy based on Model Predictive Control with Notch Filter is presented, which is applied to the mechanical resonance suppression of four-mass servo drive system. The control effect of MPC and MPC + Notch Filter under different working conditions (step speed and impact load) is analyzed. The simulation results of this paper show the following:(1)
In the process of step speed response, the composite controller method of Model Predictive Control with Notch Filter can effectively suppress the adverse effects of mechanical resonance on the four-mass servo drive system. When the speed is from 0r/min to 200 r/min, the speed overshoot of the composite controller is 10.5r/min less than that of the MPC controller, and the decrease is 90.5%. When the step speed is from 200r/min to 1000r/min, the speed overshoot of the composite controller is reduced by 14r/min compared with the MPC controller, and the decrease is 82.4%. The speed overshoot is reduced, and the control effect of speed oscillation is better than that of the Model Predictive Control.(2)
In the process from 10 Nm to 30 Nm, the speed overshoot controlled by the composite controller is 12.8 r/min, and the average transmission torque ripple is 0.12 Nm; the speed overshoot controlled by the MPC controller is 34.3 r/min; the average transmission torque ripple is 0.65 Nm. In the process of returning from 30 Nm to 10 Nm, the speed overshoot controlled by the composite controller is 10.2 r/min, and the average transmission torque ripple is 0.09 Nm; the speed overshoot controlled by the MPC controller is 32 r/min; the average transmission torque ripple is 0.68 Nm. The composite controller of Model Predictive Control with Notch Filter can suppress the oscillation of the transmission torque efficiently and have great anti-interference performance and dynamic characteristics.



---

*Source: 5627734-2021-12-22.xml*