Adjustable Speed Control and Damping Analysis of Torsional Vibrations in VSD Compressor Systems

: This paper proposes a model-based two-degree-of-freedom (2DOF) speed control for a medium voltage (MV) variable speed drive (VSD) connected to a centrifugal compressor (CC) train. Torsional mode excitations in the drive shaft due to converter switching behaviour are considered. An effective description of the harmonics transfer is proposed. The tuning strategy aims to optimize the tracking behaviour of the step and ramp command, taking care of critical speed excitations. The stability of the closed-loop dynamics against time delay and drive parameter variations are studied by means of Nyquist diagrams and time-domain simulations. A descriptive method for the process damping behaviour is proposed. The control strategy is evaluated through simulations as well as an experimental setup, based on a hardware in the loop (HIL) in a master–slave conﬁguration.


Introduction
In the oil and gas industry, voltage source drives (VSDs) are key elements from liquefied natural gas trains to subsea scenarios [1], because they adjust their speed according to the centrifugal compressor (CC) train requirements, maximizing overall efficiency and making them suitable as prime movers [2,3].However, there may be harmful interactions in the coupling shaft, such as torsional mode excitations and fatigue cycles.The intrinsic nature of switching-based drives produces a pulsating torque transfer on the transmission components that can match some torsional natural frequencies (TNFs), implying an oscillating behavior of the state variables [4,5].In fact, crack and/or lifetime reduction can happen if the torque constraints are not satisfied [6].Hence, resonant mechanical models, such as a dual-inertia elastic systems, are commonly used to design a control law [7][8][9].
Looking to the industrial feasibility of a digital signal processor (DSP), this paper proposes a 2DOF controller, based on classical PI cascade control architecture, with the aim to obtain good reference tracking and load torque rejection for each DOF.We additionally aim to achieve a MV-VSD that follows a CC speed-profile in the presence of TNFs excitation.The selection of this control architecture is due to the necessity of implementing this control on existing industrial controllers without changing the implemented control strategy.Model-based pole placement rules are used to design the feedback controller, and aimed to reject the load torque changes.A feedforward action that optimizes the speed tracking is properly selected according to the actual reference profile, avoiding excessive coupling torque vibrations.More advance controls or estimators [10] could be used to improve the control performances, but these implementations are out of the scope of this paper.
Usually, a CC profile involves dynamic ramp references with different slopes (e.g., the pressure cycle does not vary instantaneously), and more than only a step-command of the load torque (e.g., valve opening, starting another CC in parallel) [11].Hence, the case-study analysed in the paper considers a mixed-speed profile and a quadratic load torque, both coming from a high-level CC pressure map.
The paper also proposes a damping analysis, because a closed-loop strategy indirectly affects the process damping ζ P at each TNF, influencing the superimposition of the mechanical damping ζ mech and the electrical damping ζ elet .Delayed control algorithm may increase/sustain potential failure conditions [12,13].Thus, a descriptive approach to consider the time delay influence onto ζ P is given.The simulated control strategy is then experimentally validated through a hardware-in-the-loop (HIL) in a master-slave configuration.The outline of the paper is as follows.A description of the case study requirements and the related mechanical modeling are given in Sections 2 and 3. Section 4 presents the AC-DC-AC stage and the harmonics transfer, followed by the combined control design in Sections 5 and 6.Section 7 addresses robustness issues and process damping description.In Section 8, the control design is compared and verified by means of HIL test.

MV Field Case Application
This paper presents a case study of an 8 MW industrial application, which consists of a VSD directly connected to a CC train through a coupling shaft.For reasons of simplicity, this paper considers the CC as the only load stage.The VSD is based on an induction motor and a multilevel conversion unit, which comprises an Active Neutral Point Clamped 5-level (5-level ANPC) Voltage Source Inverter (VSI) driven by a 36-pulse diode-bridge rectifier.The VSD control system has to guarantee the command tracking of the CC speed profile, which implies a certain pressure injection from the load stage to the next (outer pressure/volume controls).A complete scheme is shown in Figure 1, while the parameters are given in Table 1.Let us consider for this analysis that the system should follow the speed profile shown in Figure 2. Two steady-state values are considered for testing purpose: 10% and 20% of the rated speed Ω rat .Furthermore, let us assume a non-linear load torque characteristic m load (t) affecting the CC: it follows m load (t) = f v Ω 2 load (t) + f s .In particular, considering m load (t) as a disturbance, Figure 2 also includes a critical event given by a sudden load change, e.g., a downstream valve opening at t = 50 s.

Drive train Mechanical Modeling
In a MV field, the mechanical design of the drive-train provides TNFs between 0-100 Hz; such range mainly includes the first TNF and second TNF, see [3,10].This paper focuses on a dual-inertia model, which correctly describes the first torsional mode associated to the first TNF.
The angular positions of the motor and the load are θ mot , and θ load , respectively.The motor torque is m mot , the load torque is m load and the coupling torque is m coupling .The angular speeds of the motor and the load are Ω mot and Ω load , respectively.Using the model (1), the transfer function G 31 (s) from Ω mot to m mot is defined in (2): where J tot = J mot + J load and the numerator N mot (s) and the denominator D(s) are given in (3).
Equation ( 2) can be rewritten as in (4), where the main mechanical parameters that characterize the elasticity of the drive shaft are identified: the resonance parameters ζ res , ω res , and the anti-resonance parameters ζ ant , ω ant .
In (4), ω res is associated with the poles of the denominator and ζ res is the damping coefficient; ω ant is associated to the zero of numerator and ζ ant is the damping coefficient.
These two eigenvalues for the system under analysis define the anti-resonance frequency at ω ant = 2π6.6 rad/s, ζ ant = 4.3 × 10 −4 and the first Torsional Natural Frequency at ω res = 2π17 rad/s, ζ res = 1.1 × 10 −3 .The first torsional mode associated with the first TNF at ω res is shown in a 3D plot in Figure 3.

Analysis of the VSD Vibrating Behavior
Although the motor drive pulsating torque is the main source of TNFs excitation, the vibration behaviour in terms of amplitude and oscillation decay is related to the VSD integration in the closed-loop control strategy.A complete analysis of the torsional response needs to study:

•
The AC-DC-AC topology and the torque spectrum;

•
The control loops design and the oscillation characteristics.
An electrical scheme overview of the motor drive is shown in Figure 4.As discussed in [4,11], the AC-DC-AC stage produces integer harmonics and interharmonics in torque (referred to DC-link coupling).The overall excitation content transferred to the shaft line is computed in (5), as a function of the motor speed and supply frequency.m, k are integer multipliers according to the following cases: where f mot = n pp Ω mot /60.Experimentally, the motor torque running the VSD at the rated speed Ω rat in no-load condition leads to a motor torque as in Figure 5. Performing a spectrum analysis of the air-gap torque pulsations, the highest excitations components n k are the sixth, twelfth, and eighteenth, with a torque magnitude close to the 2% of the rated torque; the other components have a negligible magnitude.Hence, the pulsating torque is mainly defined by the VSI switching behaviour, kf mot with k = 6, 12, 18.

Campbell Diagram and First Analysis
The torsional behaviour, the torque content and the application speed range information can be merged in a Campbell diagram: a 2D graphical representation shown in Figure 6.The intersection points between the TNF horizontal line and the excitations lines transferred to the shaft f excit (m, ±, k) identify the critical speeds Ω (i) crit .They can be computed for the first TNF, as in (6).According to Table 2, Ω crit (1) , Ω crit (2) , Ω crit (3) , are excited at quite low speeds: this can be a problem either during a start-up or for a speed reference in their neighborhood.
Critical Speeds Ω (i) crit (pu) (3) 0.056 (2) 0.028 (1) 0.019 Furthermore, due to the VSD closed-loop integration, the actual damping factor at each Ω (i)  crit is a combination of the mechanical one given by the shaft proprieties ζ mech and the electrical damping ζ elet mainly given by the electrical motor damping ζ mot [14].The control strategy defines the contributions of each term, and it is possible to analyze this effect with a more general definition of an overall process damping coefficient ζ P , as shown in Figure 7 and given in (7).One way to describe this coefficient is through the closed-loop transfer function of the process, where the control architecture adds additional dynamics besides G mot (s) (4).
For the speed control a 2DOF piece of architecture that includes a feedback controller R(s) aimed to improve the load torque rejection, a feedforward controller C(s) aimed to optimize the command-tracking performance [7,15] has been selected.
A full scheme is reported in Figure 8.The load speed Ω load is controlled by closing the feedback loop on the motor speed Ωmot that came out from an estimator E(s), seen as an intermediate stage between Ω mot and Ωmot , G load (s) is described through G mot and G 41 = Ω load /m mot , as given in (8).
where N load (s) = sC el + K el it is possible to define a cross-coupling term G d (s), as shown in (9), in order to consider the effect of the m load torque on motor speed Ω mot : where N d (s) = J mot s 2 + sC el + K el .

Torque Control Loop
The torque control loop transfer function T(s) considers a current control aimed to adjust the voltage input of a modulator that fed the VSI.A simple modeling splits T(s) into two transfer functions: the power stage P(s) and the harmonic stage H(s).Assuming a high bandwidth for the current controller, the VSI mainly influences high-frequency dynamics, then the power stage can be modeled as in (10): where T p = 5 ms according to the cycle time of the switching states for the system under analysis [16][17][18][19].Considering the spectrum analysis in Figure 5, the harmonic content of the motor torque m mot is represented through a critical speeds function H(s) as follows: where A is the average magnitude of the harmonic components (i.e., 2% from Figure 5).

Feedback Controller Design
The control scheme architecture is shown in Figure 8.The selected PI-feedback controller is R(s) = k p + k i /s where k p and k i are the proportional and integral gains, respectively.Closing the feedback loop on the motor speed Ω mot, the open-loop transfer function L(s) is given by ( 12): Hence, the command-tracking transfer function from A model-based pole placement rule is chosen as control design solution.This approach is based on the following assumptions: undamped condition (C el = 0), ideal actuation for T(s) = 1 and E(s) = 1.In particular, (13) can be rewritten in an equivalent form such as (14): where there are four complex closed-loop poles p 1,2 and p 3,4 defined by the coefficients (ω 1 , ζ 1 ) and (ω 2 , ζ 2 ) that are the function of the previous system parameters {J mot , J load , K el } and controller parameters {k p , k i }.Comparing ( 12) and ( 13), it is possible to write the algebraic system (15): As shown in Figure 9, to design the controller, a classical pole placement could be used assuming identical damping [15].However, a model-based design rule can be tailored to the first known TNF dynamics [20,21].Assuming p 1,2 and p 3,4 are placed into separate s-plane regions, one of them is selected as the dominant one, and then the other resonant poles are computed in order to assure enough damping.As R(s) has only two tunable parameters {k p , k i }, the poles cannot be arbitrarily placed [8], so an analytical procedure to design the controller could be: 1.
From (15), k p and k i are written as a function of the pole pair p 1,2 , based on {J mot , J load , K el , ω 1 , ζ 1 } as in (16): Both k p and k i parameters are graphically represented as two controller parameter surfaces, as shown in Figure 10.

2.
Choosing the pole-pair p 1,2 as the dominant ones, as it follows ω 1 < ω 2 .A fast closedloop dynamic is guaranteed by settling the frequency boundaries according to ω ant as given in (17).In addition, ζ 1 should be high in order to have R(s) robustly stable on load torque rejection.The coefficients of p 1,2 ⇒ (ω 1 , ζ 1 ) are chosen through the gains counter plots in Figure 10, then k p , and k i are selected, focusing on the linear parts of the surface.
The resonant coefficients of p 3,4 ⇒ (ω 2 , ζ 2 ) are a function of the previous parameters, as given in (18): This strategy leads to the coefficients of Table 3 and to the bode diagrams of the transfer function L(s) and F com (s), as shown in Figure 11.The performances of the designed controller R(s) are evaluated through two step responses: a speed step at constant torque, and a torque step at constant speed, as shown in Figure 12.

Feedforward Controller Design
It is possible to improve the control performances adding a feedforward action C(s) into the speed control.In this way, the reference-tracking transfer function F* com (s) becomes (19): In order to reduce the overshoot in the response to the load speed Ω load shown in Figure 11, the related feed forward control transfer function C(s) can be computed by imposing (20): where α is a design parameter that defines the new dynamics of the command-tracking response in terms of rising time.This leads to the transfer function ( 21): Once again, the speed step response at constant torque Ω load = 0.1, m load = 0 is evaluated.The benefits of a 2DOF strategy compared with the original regulator R(s) are shown in Figure 13, especially for α = 1 it is possible to mitigate the overshoot, while a higher α decreases the rising time with higher overshoot.However, the feed forward action C(s) should be designed, minimizing the tracking error of the closed-loop e tr (s) ( 22): According to the CC speed profile, a useful criterion is to evaluate the steady-state error e st (s) through the final value theorem in the following cases:

•
When a speed step reference Ω ref (s) = h/s is applied, where h represents the step size (e.g., h = 1): • When a speed ramp reference Ω ref (s) = h/s 2 is applied: For the ramp speed reference, to obtain e st (s) = 0 the order of C(s) should be in a higher order.Keeping a relative zero-grade, Equation (20) can be rewritten as in (25), where it is used for a second order type: This leads to a new C*(s) defined in ( 26): where the new closed-loop transfer function differs from the previous one, due to the double pole at s = −αhs.In Figure 14, the different choices of α for a speed ramp of a 1% slope are shown.Different actions, summarized by C opt (s) = F{C(s), C*(s)}, could be applied according to the given reference.Evaluating the actual slope of the speed reference Ω ref , a simple algorithm can be executed to select a proper feedforward action, or to reduce the ramp slope.A pseudocode example is given in Algorithm 1.This criterion can be extended to different reference profiles by solving an optimization problem aimed to setpoint tracking, e.g., using an H 2 approach: min + − e tr (s) 2  2 where . 2 stands for the 2-norm.
A ramp command excites the first TNF by running through the critical speeds Ω crit (i).The applied torque m coupling leads to shaft failures if the maximum oscillations amplitude M m (t) exceeds the crack boundary [10].From Table 1, it follows that this maximum value should respect (27): According to (27), a slope limit of 1% is followed for the system under analysis.If the slope is greater than 1%, according to Algorithm 1, the reference will be modified in order to have a lower ramp, even if this modification changes the settling time as a side effect.It is not useful to change the slope only where there are localized critical speeds, due to the introduction of some additional ramps with higher slopes (e.g., 1.23%) that aim to follow the reference.The comparison between these two strategies is shown in Figure 15.

Nyquist Analysis
Implementing the control algorithm in an industrial DSP leads to time delays due to continuous/discrete-time conversion and integration methods [13,[20][21][22][23]. Latencies modify the damping action on the first TNF due to the decreasing effect of the overall closed-loop stability.Time delays can be roughly split into torque-actuation delay T t , located in the inner torque control transfer function T(s), and speed-computation delay T m , located in the estimator transfer function E(s).They both can be modeled as first order transfer functions (28) (Figure 8).
T(s) : The behavior of the system can be analyzed through a Nyquist analysis performed on the open-loop transfer function L(s) (12).This is useful to compute a critical time delay T cr that gives the maximum latency after the stability is lost.Defining a total system delay as T d = T t + T m , from (12) this follows (29): . L (T d ) (s) = L(s)e −sT d for T d = 0, . . ., T cr , . . .(29) Increasing T d , the vibrating behavior changes in terms of higher amplitude and timedecay, both for critical speeds crossing and step response, as shown in Figure 16.We identified T cr ∼ = 10 ms as the critical value for the system under analysis; if T d < T cr the system results are stable, but the decaying behavior depends on T d .This sensitivity to T cr is fundamental to choose a proper sampling time for the discretization onto a DSP.

Process Damping and Characteristics
A novel descriptive method to analyze the process damping ζ P of ( 7) could be based on the sensitivity to time-delay of reference tracking transfer function F* com (s) [24][25][26][27][28].The denominator D F (s) of this transfer function has further terms in addition to G mot (s), with the closed-loop poles characterized as in (30): If some mechanical parts are changed (e.g., higher CC size = J load , bigger motor = J mot , lifetime of the shaft = K el , C el ), they lead to different {k p ,k i } until the worst case is reached: a (ω 2 , ζ 2 ) close to (ω 1 , ζ 1 ).Indirectly, the mechanical characteristics can be a positive/negative help to damp the first TNF through ζ mech .Similarly, the intrinsic behavior of the motor can add an electromagnetic damping ζ elet .Such dynamics can be model as a resonating transfer function given in (31), with an electrical resonance at ω elet .
ζ elet and ω elet can be computed from an FEM analysis of the electrical side only [2,26,29].Usually, in CC trains ζ elet is negligible compared with ζ mech , which is higher due to large J mot , J load .In addition, ω elet is located at higher frequency, on the contrary to ω mech .However, even without changes in the electrical or mechanical parts, ζ P changes because the control strategy affects the impact of ζ mech , ζ elet .Considering the latency effect of T d , from (19), it follows (32): As shown in Figure 17, increasing the time-delay, the magnitude of the resonance is increased, and this produces a higher amplitude of the first TNF oscillations.As a first approximation, the resonating behavior at s = j ω res can be described as in (33), where µ = 1 for simplicity., as in (34).
where ζ P (T d ) is computed using a polynomial curve fitting from samples of (34), changing T d = 0, . . ., T cr inside the stable region, as shown in Figure 18.This approach allows us to consider the electrical and mechanical design, in addition to the control topology.For example, changing the motor inertia to J mot * = 10 J mot , this follows ω* res = ω res and the stability boundary increases from T cr to T* cr .As shown in Figure 18, with the same amount of delay, ζ* P (T d ) is less sensitive to T d .

HIL Experimental Setup
The experimental setup consists of a HIL in a master-slave configuration.A master board is used to implement the control algorithm, and dynamic models for VSI (the DClink is assumed constant), motor and estimator units.The control algorithm is partially translated from MATLAB/Simulink into C-code, and completed with an assembly-code, then compiled onto a Motorola DSP 24-bit word.A slave board is used to implement the mechanical model, process boundaries and load a torque profile.The boards allow four sampling times T s from 25 µs to 2 ms, used for backward-discretization.The most sensitive operations are computed at the highest T s .The communications are based on fibre-optic links.This HIL provides a rapid prototyping platform to test the different speed profiles earlier analyzed in simulation.An overview is shown in Figure 19.
The CC reference profile given in Figure 2 is evaluated with different C(s) and C*(s), as shown in Figure 20.The simulation agrees with the HIL results: there are similar torsional vibrations in terms of amplitude and time-decay both on critical speed crossings and step responses.The torque m coupling constraint is respected due to a speed ramp of 0.5% and 1%, respectively, validating the early analyzed approach.The process damping is ζ P (7 ms) ∼ = 0.3.For effective simulations, an additional white noise generator should be added because the HIL presents a greater noise than actual simulation, i.e., due to hardware implementation delays.Comparison between HIL and simulation results are reported in Table 4.
In addition to time-domain simulations, the power spectra of the coupling torque is evaluated through a waterfall diagram.The diagrams in Figure 20 exhibit the spectrum variation over time, by plotting a series of frequency-domain spectrums calculated during a consecutive time window of 15 s.This method is accurate if it is considered as a steady-state signal.By applying a Hanning window in a range 55-70 s, it follows a noisier waterfall diagram for the HIL than the other, but there are clear first TNF peak lines for both.Step command M m ∼ = 3.5% Noise level % ∼ =1.6%

Conclusions
This paper presents a complete design procedure and performance analysis of a motor drive control, in order to minimize oscillations and improve system reliability.A modelbased pole placement design method and an adjustable feedforward action of a 2DOF speed control have been presented.Sensitivity to time delays, noise and parameter changes were analyzed by means of Nyquist diagrams, the process damping descriptive method, timedomain simulations and waterfall diagrams.Simulations and HIL experiments showed that the designed controller achieved good reference tracking for both step and ramp commands, as well as a robust load torque rejection to valve opening for an industrial compressor system.

Figure 4 .
Figure 4. Exploit of the electrical topology of a VSD.

Figure 6 .
Figure 6.Campbell diagram of the analyzed system.

Figure 7 .
Figure 7. Damping influences at critical speeds: critical frequencies representation in the Campbell and Bode diagrams.
(s)  and leads to the transfer function (13):

Figure 11 .
Figure 11.Bode diagram with k p = 10 and k i = 6.5.(a) Open-loop transfer function L(s) Bode diagram; (b) command-tracking transfer function F com (s) Bode diagram.Red lines represent rigid behavior.

Figure 12 .
Figure 12.Step responses with R(s).(a) Speed reference step changing response; (b) load torque step changing response.

Figure 13 .
Figure 13.First order transient comparison according α in terms of speed response and torque response.

Figure 14 .
Figure 14.Second order transient comparison according to α, in terms of speed response and torque response.

Figure 15 .
Figure 15.Energy consideration for ramp speed reference.(a) Speed and torque system behavior; (b) speed reference possible modifications.

Figure 16 .
Figure 16.System response (torque and speed) and stability analysis with different T d : (a) T d = 4 ms; (b) T d = 10 ms; (c) T d = 11 ms.

F 2 +
2ζ P ω res s + ω 2 res s=jω res , (33) Exploiting the right side of (33) follows a description of ζ P (T d ) as a function of T d through F (T d ) com (s) s=jω res

Figure 18 .
Figure 18.Process damping description as a function of T d .

Table 2 .
System critical speed.