Simply structured controllers for vibration suppression in long rotors

In this paper suppression of the transient flexural vibrational disturbances in long rotors, with fluid film bearings, is investigated. The rotor is described by a series of distributed shafts connected by the lumped discs, and the system is mounted on lumped fluid film bearings. Upon determination of the dynamic stiffness matrix of the system, the best approximate transfer function matrix description of the rotor, is determined. Initially vibration suppression by simple diagonal Proportional + Integral (PI) controllers is studied and via direct search optimisation techniques the PI parameters which exhibit fast vibration suppression is found. The resulted high integration rate, and low proportional gain PI controller, theoretically provided fast suppression time. However, it is shown that due to the strong coupling effect in the rotor system, and high rate of integration, the closed loop relative stability is weak, and feasibility of controller is questionable. Therefore, an alternative simple first order controller without integration action, that is named “attenuation filter “is suggested that can produce stronger stability and produces significant (not full) vibration suppression. The closed loop multivariable control of the rotor system comprising two vibration sensors and two magnetic actuators using such attenuation filter, is then simulated. The response to step disturbances, has provided 95% suppression with significantly fast response. It is concluded that although the attenuation filter may not provide 100% suppression, but it more reliable since the integration of the error, that results weak stability is avoided.


Introduction
It is nearly four decades that research on active vibration control is continuing to accommodate the needs for providing reliable aircraft engines and active anti-seismic structures. Although the technology has progressed significantly, applicable and commercially viable solutions, are few. The two major control strategies are electromagnetic excitation type and described in [1]. And wrapped piezoelectric patches on rotors that is still under research and development [2].
Since there are several sources that can cause excessive vibration in rotors, the suppression of unwanted vibration by unique controller will be sophisticated. Recently, a model free control (MFC) based on adaptive machine learning [3] is tested for such purpose. There are not simply structured controllers except in limited circumstances, where PD controllers can be used [4].
Traditional methods for vibration control are usually passive and, rely upon viscous dampers or dynamic vibration absorbers [5]. However, the synthesis of active controllers for vibrating systems, was investigated in the early seventies using modal control theory [6,7]. Such works were extended to include active control of flexible rotors using state variable feedback [8], where An approximation to behaviour uses discrete models with lumped masses.
Difficulties arise, when higher order descriptions of behaviour which are derived from finite element analysis of rotor-bearing system [9], are employed in these types of investigations. The complexity initiated, from assigning a large number of eigenvalues (fifty pairs for example) this results in high order, non-realisable controllers [10] especially when observer based controllers are designed.
In order to achieve a realisable control strategy for lumped of rotor-bearing systems, a vibration performance index in the frequency domain has been defined by Burrows [11]. The optimal location of the control force was computed to minimise this index and was confirmed using experiments with a magnetic control unit [12]. Implementing this type of control is often expensive and impractical in rotating machinery applications, since it requires moving actuators.
Practical control scheme known as Adaptive Force Balancing (AFB) can be implemented and described by Shafai et al. [13]. However this strategy is limited to removal of steady state vibration, thereby achieving "on-line" balancing. This has been applied to single variable and multivariable models of a rotor supported on magnetic bearings, that enable removal of steady state vibration of a given frequency between them.
In order to increase the range of frequencies, further theoretical and experimental studies was carried out by Kenospe et al. [14]. They proposed a multi-tasking computer program, that performs AFB and feedback loop computations, simultaneously, and computes the required control forces in the magnetic bearings.
Investigations into removing loads due to shock or transient vibration as given by sudden blade loss or unpredicted disturbances, is very limited due to its complexity. Palazzolo et al. [15] use optimal control theory for this purpose, to search for an optimal state feedback rather than seek to reject disturbances. Hori et al. [16] consider a simple torsional system, with State Feedback and Load Acceleration Control (SFLAC) to suppress the torsional vibration. Their simulation results do not indicate appropriate vibration rejection, since their state feedback strategy considers pole placement rather than disturbance rejection.
Hathout and Shafei [17] have considered a simple Jeffcott rotor model (a hypothetical distributed shaft with its mass concentrated in the middle), they show that a Hybrid Squeeze Film Damper (HSFD) coupled to a servomechanism with PI controller can reduce the rotor vibration. The results of their simulation do not indicate a significant degree of rejection. This arises because a single control element is applied to a two degree of freedom system.
In all previous investigations two major assumptions are ignored. First, the rotor ought to be modelled as a distributed-lumped parameter system, rather than purely lumped. Moreover, the rotor-bearing model can be made multivariable, with multiple sensors and control force actuators. Apart from that, stationary actuators are necessary prelude in order to implement practical control strategies in long rotor systems. The single moving actuators used in [11,12] is hard to implement in practice.
Aleyaasin et al. [18] have shown that both assumptions can be effectively combined into a single system. Although the actual distributed-lumped rotor model would be more complicated, by implementing dynamic stiffness matrix method, an appropriate approximation for describing a multivariable model can be obtained. Moreover in a review on various control strategies [19], it is proposed that the response of a low order PI controller designed by optimisation is much better than a high order H ∞ controller [20]. This paper sets out to explore the above issues in more detail.
In particular a rotor that is mounted on fluid film bearings, and is described by a series of distributed and lumped elements. Upon determination of the dynamic stiffness matrix of the system, the best multivariable transfer function matrix description of the rotor is computed. It is shown that vibration disturbances in the rotor could be suppressed by simply structured controllers. The first type investigated in this article is diagonal PI controllers, that can provide 100% suppression, but they suffer from stability weakness. To overcome this weakness, another first order controller named "attenuation filter" is also investigated and can produce 95% suppression with stronger stability.
To find parameters of the PI controller and also the attenuation filter, that ensures stability together with reducing the suppression time, a powerful pattern search algorithm based on the Hooke and Jeeves (HJ) technique is implemented. By considering the location of the dominant closed loop pole as an objective function, and the closed loop system stability as a constraint, the range of PI parameters and attenuation filter are determined, that can reduce the suppression time and also increase the suppression percentage. This is then applied to a system comprising two vibration sensors and two magnetic actuators. The response to a step disturbance, is shown to be quick. However, the examination of the full Nyquist plot of multivariable system shows that the relative stability for PI controller is very weak. Therefore, PI controllers that are useful in process industries, cannot be successful in vibration suppression particularly in long rotors.
The alternative controller, which is the combination of minimum effort controller and an attenuation filter, does not provide full disturbance rejection (100%), but can supress the vibration significantly i.e. (95%) in simulation. Since there is not integrator in the structure of alternative controller, the stability weakness, does not exist in attenuation filter. where G � lm (j ) in (1) is the frequency response function, in which the corresponding flexibility matrix is:

Multivariable models and transfer function estimation
The above partitioned flexibility matrix is derived in [18] based on dynamic stiffness matrix method (DSMM) in which each of the sub-matrices are described in appendix.
Frequency response data can be obtained from each transfer function, so that a vector of frequency response data, as a function of , could be written as.
From the multivariable frequency response matrix in Eqs. (3) and (4) estimates of the transfer functions can be obtained. An elementary algorithm for the computation of the estimated second-order transfer functions is given in [18].
The simple route for obtaining low order estimated transfer functions arises from the use of frequency domain identification techniques. A program in the MATLAB named MFD toolbox [21] is available, enabling a proper transfer function of any order to be estimated from frequency response data. This satisfies both the amplitude and phase of Eqs. (3) and (4). However, the estimation procedures are based on least squares techniques [22] and there is the possibility of achieving unstable transfer functions, even if the stability of the system is guaranteed, as explained by the procedure used in [18]. Therefore, the user should choose the lowest possible order transfer function, commensurate with steady state and transient fidelity, by matching data from Eqs. (3) and (4). Such transfer functions could be represented by the form:

Controller synthesis by optimisation
Upon the computation of an estimated transfer function matrix G(s) a control system could be determined, for the multivariable configuration, shown in Fig. 2 where the D(s) is the vibrational disturbance vector, E(s) is the error vector and Y(s) are the vibrational displacement, measured by the The necessary condition for steady state vibration rejection can be expressed as In order to connect the frequency and time domain, the final value theorem in multivariable version should be used [23]. When step vibrational disturbances are applied to Eq. (9), it results the following expression between time and frequency i.e. mapping from t = ∞ into s = 0: Substituting (10) into (11), then to achieve a non trivial solution for the system of linear equations in (11) we can conclude that [23]: To satisfy Eq. (12) the compensator K(s) should be constructed, in diagonal form with PI (Proportional + Integral) elements as follows: The parameters of each PI element could be determined to ensure stability and also to produce low disturbance rejection times. For this purpose direct search optimisation techniques are employed. The procedure consist of a series of iterations in which the closed loop transfer function matrix is computed using: The iteration procedure is based on choosing the initial gains k i and integral rates k p then by using the symbolic operations in MATLAB, the characteristic determinant of the system can be determined and represented as: The location of the dominant closed loop pole of the system can be determined as an objective function in the optimisation procedure, which could be shown by: Since all of the closed loop poles have negative real parts, the minus sign in (16) is considered to find the location of domain poles. An efficient search method is then required to continue the iteration until the dominant pole location Δ in (16) is minimised enabling the required low disturbance rejection time to be achieved. In order to determine realisable controllers, the following constraints should be imposed on the parameters of the PI elements: Since the objective function for this problem is complicated, and there are several searching variables, gradient based optimisation will be formidable task. Therefore, direct search optimisation techniques are the best route for this case. Direct search optimisation techniques do not form the fastest vehicle, but their convergence is not too sensitive to the initial starting values of the variables. An efficient and ingenious direct search technique was proposed by Hooke and Jeeves [24]. The method is based on a pattern search scheme which is described by the flow chart in Fig. 3. The search consists of a sequence of exploration steps about a base point, which if successful is followed by pattern moves. An exploration step is described by the flow chart in Fig. 4. To satisfy the constraints, a suitable penalty function [25] could be added to the main objective function.

Numerical example and application
A rotor bearing system which is shown in Fig. 5 consists of three lumped steel discs each with thickness t = 50 mm, diameter of 100 mm connected by four distributed shaft elements each with length of L = 400 mm and diameter This example is also considered by Burrows and Sahinkaya [11]. The rotor with angular velocity of Ω = 360 rad/s is supported by two identical fluid film bearings each with the length l = 16.9 mm. The radial clearance of the shaft bearing is c = 0.127 mm and the lubricant viscosity is μ = 0.015 N·s/m 2 . Using the above data the stiffness and damping of the bearings could be computed according to the procedure outlined in the "Appendix 2".
For the system in Fig. 5, in the case where three sensors are used to measure the lateral vibration of the lumped discs, while the control force is applied through the bearings, considering Eqs. (1) and (2) the system of equations for the flexural vibrations in (52) gives: The actual G′ ij found via distributed-lumped modelling is irrational (not polynomial) because it is related to flexibility matrix in Eq. (2). In order to find open loop poles it should be approximated by polynomial type G ij . In MFD toolbox a powerful identification algorithm [22] exists that approximates a G ij polynomial form with poles and zeros.
Using the elements of the multivariable frequency response matrix displayed in Eqs. (18) and (19) and the algorithm for the estimation of the closest low order transfer functions, which is described by Eq. (5), the transfer function matrix of the rotor bearing system can be computed. In the configuration where two sensors are located at the discs 1 and 2, while actuators are placed at bearings 1 and 2, the system is MIMO with two inputs and two outputs, hence: where the transfer functions in (20) are:  All the above estimations satisfy both the amplitude and phase angle equivalence. For example in Fig. 6, the estimated transfer functions (21) are very close to the actual transfer function, where both the amplitude and phase angle of the actual and estimated transfer functions are drawn.
These figures also show that when ω = 112 rad/s there is a peak in the frequency response curves. The transfer function matrix in [18] for the above system is different, since it is an approximated 2nd order system fit with poles only (no zeros in individual transfer function). However, herein the type of fitting above for the "same system" is more accurate and is carried out by Multivariable Frequency Domain (MFD) toolbox that provides the closest fit as possible such that both poles and zeros exist in the open loop system. The MFD toolbox identified G 11 as an order 3 system, from the general form expressed for any order in (5).
The G 11 is very close to (approximated value) of the G`1 1 which gives the displacement Y 1 resulted by unit force of Q 1 . This can be obtained from the following ration by using the elements of the flexibility matrix (inversion of stiffness matrix) in Eq. (2). The elements of the stiffness matrix can be found from the Eqs. (38) and (48) in "Appendix 1".
Attempts to find a diagonal compensator with PI elements in order to supress vibrational disturbances in the rotor system, described by the transfer functions (21) to (24) resulted in a low proportional gain together with a high integral rate. Elsewhere, in the gain space it was not possible to find PI parameters which would provide closed loop stability. In Fig. 7, a three dimensional graph shows the location of the dominant closed loop poles of the system versus the parameters of the PI elements. From this graph it can be concluded that, the dominant pole of the system could be located in the right hand side of the complex plane, with a high proportional gain.  Fig. 6 Frequency responses indicate actual G 11 and estimated G11 Fig. 7 Location of the dominant closed loop pole versus proportional gain and integral rate A controller matrix consisting of PI diagonal elements was calculated using the Hooke and Jeeves search method. Table 1, shows the search results including the initial PI parameters, number of iterations, and the final resulting PI parameters which are used to obtain the minimum rejection time when a step disturbance is imposed on the rotor system. The resulting PI controller can be expressed as: In Fig. 8, the block diagram of the multivariable closed loop configuration for the rotor-bearing system is shown, where two vibrational step disturbances are applied at positions 1 and 2, simultaneously. The physical explanation for the high integral rate used herein, is because of the existence of high frequency oscillations in the system model, which requires substantial control effort for disturbance rejection. The simulation is for the worst case, where double step disturbances are introduced. The responses are shown in Fig. 9 where a rapid suppression time of about 2 s. is achieved. It should be noted that the response would be faster if the disturbance originated from one position only. The suppression time in Fig. 9, is in transient phase and after that phase, the steady state vibration will remain if the unbalance force has not been removed or cancelled by active balancing or AFB.
However, the low proportional gain with high integral rate for stabilising this system is unusual and the frequency domain interpretation of this action requires further explanation. Unfortunately, due to strong coupling in the system, and also the huge growth of the vibration amplitude at the resonance frequency, the stability information cannot be assessed via Gershgorin bands. In the 2 × 2 MIMO system in this article, the radius is q 12 for q 11 band and q 21 for q 22 band and q = (GK) −1 .
In Fig. 10 these bands, which are based on the inverse Nyquist plot of the open loop system are shown for the case where feedback reversal is not employed. Both of the bands include the origin, and therefore stability cannot be investigated [23]. Step r1 Step r Step 2 Step 1

Fig. 8 Simulation block diagram of the multivariable closed loop control of rotor system
In order to investigate the stability condition Nyquist's stability theorem could be employed. This deals with the closed loop characteristic equation of the system: The full Nyquist plot will be drawn by using (26). The Gershgorian band in Fig. 10, shows that in vibrational systems (long rotors) "diagonal dominance" cannot be observed in the band. This the only tool illustrate the strong interaction in the system. The full Nyquist plot of the system is shown in Fig. 11, where the mapping starts at = −∞ on the positive real axis and goes up to the second quadrant at = 0 and makes a infinite semi-circular clockwise turn (which excludes the origin), and from the third quadrant progresses until it crosses the positive real axis again at = +∞ . This curve does not make any encirclement around the origin, and it can be concluded that there are no right hand plane closed loop poles, indicating closed loop stability. The plot also shows that small model perturbations may cause instability.
In Fig. 11 the design route for PI controller in which the full suppression is guaranteed, is based on the pole placement in frequency domain, and can be designated by objective function (16). This leads to low "suppression time" similar to "settling time". However, the huge control force that is resulted from error integration cannot remove the stability weakness. Therefore, an alternative controller will be investigated in the next section in which "error integration" is avoided.

Alternative attenuation filter
To overcome weak stability performance of PI controller and simultaneously achieve a significant (not full) amplitude suppression, an alternative controller is indicated in Fig. 12 is suggested. This includes a filter that is a simple time delay located in the forward path of block diagram. The filter is defined by the gain K and time delay T, which is multiplied by the gains k 1 and k 2 . Then the gain vector will be converted into following functions: The gains 5.2 and 5.3 in equation is the result of two step controller design that is suggested by the author. In the first step the minimum effort control strategy that is described in [26] is implemented.
The polynomials N(s) and D(s) in Eq. (31) can be determined in terms of the filter parameters K and T as follows:   The roots of N(s) give the poles of the closed loop system versus the filter parameters of K and T. An advantage of the minimum effort design procedure herein, can be realised by finding a range for the gains k 1 and k 2 for which the system is stable. In Fig. 13, the location of the dominant closed loop poles, which are designated by their negative real parts versus the filter parameters K and T, are shown.
This enables a direct search optimisation technique to be employed for obtaining the filter gain K and time delay T, for which the minimum vibrational amplitude occurs. The objective function of the search could be the sum of the absolute values of the steady state amplitudes, expressed by: Via the simulation of the closed loop control system, the minimum value of the objective function ∆ min can be achieved after several iterations.
The objective function for this problem cannot be expressed as a function of K and T and requires the system simulation at each trial, although only two searching variables are used. Direct search optimisation techniques still form the best solution approach for this case, since gradient based optimization techniques are a formidable task due to the complexity of objective function in Eq. (34).
Table2, shows the search results including, the initial gain and time delay, number of iterations and the final gain and time delay which are selected as the attenuation filter parameters.
Using this filter the corresponding values of the steady state amplitude of each output are: The simulation block diagram of the closed loop minimum effort control system with the attenuation filter is shown in Fig. 14, where two step disturbances are applied at the system outputs. A similar situation occurs in wide span rotors during transient phase of the sudden blade loss, when the disturbance occurs at multiple locations across the rotor system.
The simulation results are shown in Fig. 15 where the responses due to the unit step disturbances are plotted. The vertical axis, which shows vibration displacement, has the same scale as the disturbances.
Therefore, the curves indicate directly the reduction achieved by using this control scheme. According to (35) the average equation vibration reduction is about 95% which is remarkable. This reduction occurs rapidly and the amplitude of the oscillations are not visible. It should be noted that Fig. 15, does not indicate full disturbance rejection, as achieved in PI controller. However, a significant reduction in vibration disturbances is obtained with minimum effort control only, without using integral control elements.
For the attenuation filter herein the error integration is avoided, and the main objective functions (34) is "suppression level" similar to "steady state error". Moreover, the physical system limits the real part of the dominant pole and the important phenomena of "actuator saturation" should be considered. In this paper a basic assumption is that actuators are strong enough and will not be saturated.

Conclusions
This paper shows that in distributed-lumped rotor systems, vibrational disturbances can be suppressed by simply structured first order controllers. When the number of sensors is the same as the number of actuators in the control system, suppression could be achieved by using simple controllers either in form of diagonal PI compensators or nondiagonal attenuation filters. The PI controllers require high integral rate to obtain low suppression time, such that the closed loop system, provides significant control forces to supress the vibrational disturbances completely. However, the relative stability for PI controller, according to Fig. 11 is very weak.
It can be concluded that when a blade loss or a sudden shock occurs in a rotor system, the full vibrational disturbance suppression can be achieved in theory. Although, it is not reliable enough, to be implemented in practice. This is because a slight change in PI parameter can cause instability i.e. rotor -stator clash.
Thereafter, to overcome this an alternative controller strategy is used by combination of minimum effort controller with an attenuation filter used in Sect. 5 of this article. Although this alternative controller cannot produce full suppression (100%), but significant suppression (95%) can be achieved. Therefore, reliable (more stable) controller with 95% suppression is superior to unreliable (weakly stable) with 100% suppression and can be recommended in practice.

Appendix 1
The corresponding block diagram for the system is shown in Fig. 16. For the distributed part there are left and right hand terminations of eight parameters. These are respectively denoted by L and R subscripts and are: displacements The input to the lumped elements is the output from the distributed elements and, therefore, the left and right terminations for the lumped part are specified by the above eight parameters. The equations of motion for the lumped mass considers rigid body motion, while for the distributed element, these are derived from Euler beam model.
The following relationships exist for the input and output of each distributed, and lumped element, as give by detailed analysis in [18]. The governing equations for a distributed element are succinctly written in terms of Γ = s √ L 0 C 0 with L 0 = A and C 0 = 1 E I , where A is the cross sectional area, E is the modulus of elasticity and I is the moment of inertia of the cross section in bending.
where Φ (i) d is the 8 × 8 matrix: With the submatrices: For the lumped parts, the gyroscopic moments and rotary inertia of the discs are included which couples the in y and z coordinates. The corresponding governing equation is described by: where the Φ (i) l is taken from [27]: with: J p and J t are the polar and transverse moments of inertia of the lumped disc. In Eqs. (1) to (6) s is the Laplace variable, m is the mass of each lumped element, Ω is the rotational speed of the shaft and is the density of the shaft material. L is the length of each distributed element. The transfer matrix for the fluid film bearings could be The damping of the fluid film bearings is captured by a spring and dashpot, as indicated in Fig. 1  and for the second lumped bearing: From Fig. 1, when the rotor bearing system, consists of n distributed shaft elements connected by n lumped discs and is supported on two bearings at the left and right ends, using Eqs. (49)-(51), the equations for the flexural vibrations of the rotor could be expressed by a dynamic stiffness matrix form. When the external forces applied at ith lumped element, is designated by and the resulting displacement designated by . Then the equations of motion could be expressed in terms of the flexibility matrix as follows:

Appendix 2
For the computation of the stiffness and damping coefficients of a journal bearing the following data is required.
(1) length of the journal bearing l. Then the load carrying capacity of the bearing would be where: The variable ε in (54) is the eccentricity ratio.
The following successive steps are necessary: Step 1-using l, R, c, μ, and ω compute: Step 2-when the load carried by bearing P is known the eccentricity ration ε can be determined by finding the roots of the fourth order equation.
Step 3-compute the dimensionless parameters.
Step 4-compute the angle ψ and further parameters. (52) Step 5-compute the stiffness of the bearing from: and damping coefficients of the bearing from.