Abstract
The analysis of dynamics of rotordriven artificial ventricle (VAD) was conducted in the work. A comparison of two types of control is given: a linearquadratic (LC) optimization and PIDregulator. It was shown that LC – control allows the pump rotor positioning with an accuracy of 0.2 mm at speeds ranging from 5,000 to 12,000 rev/min.
1. Introduction
The problems of heart failure are particularly acute in recent years. One of the variants of solution became circulatory support devices –ventricle assist devices (VADs). In accordance with articles [13], the axial VADs are preferred nowadays. The work is devoted to describing the VAD rotor dynamics: the theoretical bases were shown; the equation of motion was obtained; active magnetic bearings were considered. The question of the rotor movement control was considered: two types of control were selected for comparison, PID control and LCmanagement. The numerical simulation of the rotor stabilization was realized. The rotor positioning error had not to exceed 0.2 mm. The behavior of the rotor and the control response is examined in the speed range from 5000 rev/min to 12000 rev/min.
2. Statement
The symmetric homogeneous rigid rotor rotating along the longitudinal axis at a constant angular velocity in two radial active magnetic bearings AMP1 (A) and AMP2 (B) of axial pump VAD is considered. Specifications are given in the tables (Table 1).
Table 1Rotor technical characteristics
Nomination  Symbol  Value 
Diameter, m  $D$  15.6×10^{3} 
Length, m  $l$  24×10^{3} 
Mass, kg  $m$  12.42×10^{3} 
Positioning accuracy, mm  $\epsilon $  ≤ 0.2 
Rotation frequency, rev./min.  $\mathrm{\Omega}$  500012000 
Length from rotor center mass to bearings, m  $a$  9×10^{3} 
$b$  9×10^{3}  
Length from rotor center mass to censors, m  $c$  10×10^{3} 
$d$  10×10^{3} 
The assumptions for the equation of motion are the following [4]:
1) The rotor is symmetric and rigid.
2) Deviations from the reference position are small in comparison with the rotor dimensions.
3) The angular velocity $\mathrm{\Omega}$ of the rotor about its longitudinal axis $z$ is assumed tobe constant.
The inclinations and the angular motion around the rotor spin axis are described by the three socalled Cardan angles $\alpha $, $\beta $, $\gamma $. Linearization leads to characterizing the angles $\alpha $, $\beta $ as inclinations about the $X$ and $Y\text{.}$ The equations of motion are given for the variables $\mathbf{q}={\left\{\beta ,{x}_{S},\alpha ,{y}_{S}\right\}}^{T}$.
3. Model and equation of motion
The equations of motion follow from Lagrange’s equations:
with the kinetic energy $T$ and generalized forces ${Q}_{i}$, ${q}_{i}$ – the generalized coordinates. The kinetic energy $T$ is:
where ${\dot{x}}_{S}$, ${\dot{y}}_{S}$, ${\dot{z}}_{S}$ – components of center mass velocity, ${{\rm I}}_{S}$ – rotor inertia tensor $\mathrm{\Omega}$ – angular velocity vector.
Equation of motion of the rotor with active magnetic bearings [4]:
where $\mathbf{M}$ [4×4] – symmetric positive definite mass matrix, $\mathbf{G}\mathbf{}$[4×4]  skewsymmetric gyroscopic matrix, $\mathbf{K}$ [4×4] – stiffness matrix, ${\mathbf{B}}_{q}$ [4×4] – transformation matrix which relates the generalized coordinates of the rotor center mass to the rotor displacements in magnetic bearings, ${\mathbf{K}}_{i}$ [4×4] – matrix of current stiffness’s, $\mathbf{i}\left(t\right)$ [4×1] – vector of currents in magnets, ${\mathbf{F}}_{ext}$ [4×1] – vector of generalized external forces.
The rotor receives the load ${\mathbf{F}}_{ext}$ as the force of gravity and the moments from the hydrodynamic forces in the fluid flow:
where $\mu $ – the blood viscosity, $\mu =$ (34)×10^{3} Pa∙s at 37 °С, $C$ – the drag coefficient, $C=$ 0,910,85, $R$– radius, $l$– length, $\mathrm{\Omega}$ – rotation frequency. The impact of external influences on a person is taken into account in the form of vibration, expressed by harmonic functions, acting at the $x$ and $y$ axes: ${A}_{1}$, ${A}_{2}$ – the amplitudes of the oscillation of transport.
4. Control types
4.1. The decentralized PIDcontrol
The local control shown in Fig. 1 feeds each local sensor signal back to the corresponding bearing control current using the feedback gains [5]. The four output signals $\mathbf{i}\left(t\right)$ combine in the output vector $\mathbf{q}$:
with $\mathbf{P}$, $\mathbf{D}$, $\mathbf{I}$ – the diagonal matrixes of control coefficients. The equation of motion then takes the form:
where ${\mathbf{K}}_{C}={\mathbf{B}}_{q}{\mathbf{K}}_{i}\mathbf{P}\mathbf{C}$ and ${\mathbf{D}}_{C}={\mathbf{B}}_{q}{\mathit{K}}_{i}\mathbf{D}\mathbf{C}$ are stiffness and damping matrixes respectively and ${\mathrm{{\rm I}}}_{C}={\mathbf{B}}_{q}{\mathbf{K}}_{i}\mathbf{I}\mathbf{C}$ – matrix of integral components, $\mathbf{C}$ – transformation matrix.
4.2. The linear quadratic method
Linearquadratic regulator (LQR) – optimal control algorithm based on the idea of minimizing a certain functional [6]. The novelty of this method is based on varying parameters which can be chosen for different types of conditions. Besides, for this type of control the equation of motion will change, because it should be implemented by bearing coordinates:
In the standard form:
$\mathbf{A}=\left[\begin{array}{cc}0& \mathbf{E}\\ {\mathbf{M}}_{b}^{1}{\mathbf{B}}_{q}{\mathbf{K}}_{S}& {\mathbf{M}}_{b}^{1}{\mathbf{G}}_{b}\end{array}\right],\mathbf{B}=\left[\begin{array}{c}0\\ {\mathbf{M}}_{b}^{1}{\mathbf{B}}_{q}{\mathbf{K}}_{i}\end{array}\right],{\mathbf{C}}^{\mathrm{*}}=\left[\begin{array}{c}0\\ {\mathbf{M}}_{b}^{1}\end{array}\right].$
Formally control problem can be represented as follows [6]. It is required to find a control law $\mathbf{i}\left(t\right)$ that minimizes the objective function:
where $\mathbf{z}\left(t\right)={\left\{{\mathbf{q}}_{b}^{T}\left(t\right){\dot{\mathbf{q}}}_{b}^{T}\left(t\right)\right\}}^{T}$ – is the solution of the system Eq. (9). $\mathbf{Q}\in {R}^{8\times 8}$ is positive semidefinite and $\mathbf{R}\in {R}^{8\times 8}$ – positive definite matrixes, respectively, obtained in accordance with Bryson rule [10]. The control law will be found as:
The equation of motion becomes:
5. Modeling based on PID control
The study of rotor oscillations at fixed positional stiffness of bearings for different rotor speeds was conducted. Yellow line shows oscillations at a rotor speed of 5000 rev/min, the purple – at a speed of 10000 rev/min, green – 12000 rev/min.
For a fixed positional stiffness ${k}_{S}=$ 10^{4}^{}N∙m^{1} the results are illustrated with Fig. 1.
6. Modeling based on LQR method
The comparison of the rotor displacements was held at different speeds and constant stiffness. The value of the positional stiffness ${k}_{S}$ is 10^{5}^{}N∙m^{1}, of the current stiffness ${k}_{i}$_{}– 10 N∙А^{1}. The black line shows changes of displacements at 5000 rev/min, blue line – at 10000 rev/min, red line at – 12000 rev/min (Fig. 2).
Table 2Comparison of the displacement amplitude values for different rotation speeds and positional bearing stiffnesses
Rotation speed, rev/min  Max values of displacements ${k}_{S}=$ 10^{4 }N×m^{1}, m  Max values of displacements ${k}_{S}=$ 10^{5 }N×m^{1}, m  Max values of displacements ${k}_{S}=$ 10^{6 }N×m^{1}, m 
5000  3.9×10^{7}  3.6×10^{7}  2.54×10^{7} 
10000  4.8×10^{7}  3.8×10^{7}  2.6×10^{7} 
12000  5.21×10^{7}  3.93×10^{7}  2.71×10^{7} 
Fig. 1Displacements of the section centers A and B
Fig. 2Displacements of the section centers A and B
It can be concluded that with the increase of rotor speed, values of the oscillation amplitudes and of the control currents increase too. However, the obtained values are permissible. The values of currents in the magnets arranged along the $y$axis are greater than in magnets along the $x$axis. It is due to the fact that in addition to the influence of the hydrodynamic moments of the blood flow the force of gravity acts along the $y$axis.
The results obtained for the three types of control have been tabulated (Table 2) for the velocity of 10000 rev/min and 12000 rev/min respectively.
Table 3Comparison of rotor center displacements at different control types at 10000 rev/min. The value of the positional stiffness kS is 105N∙m1
PIDcontrol  LQRmethod  
${x}_{bAmax}$, m  3.7×10^{7}  5.4×10^{7} 
${y}_{bAmax}$, m  3.9×10^{7}  5.6×10^{7} 
${i}_{bxmax}$, mА  5  5.3 
${i}_{bymax}$, mА  8  5.5 
7. Conclusions
The simulation results show that LQR method meets the requirements of the rotor control the best way. It provides the position of the centers of rotor sections within the permissible error – 0.2 mm and allows to optimize the control on several criteria, in this case, criteria were: the position of the bearing section centers and control currents. The results of these studies can be used for real axial pump design.
References

Birks E. J. Left ventricular assist devices. Heart, Vol. 96, 2010, p. 6371.

Thunberg Christopher A., Gaitan Brantley Dollar, Arabia Francisco A., Cole Daniel J., Grigore Alina M. Ventricular assist devices today and tomorrow. Journal of Cardiothoracic and Vascular Anesthesia, Vol. 24, Issue 4, 2010, p. 656680.

Griffith B. P., Kormos R. L., Borovetz H. S., Litwak K., Antaki J. F., Poirier V. L., et al. HeartMate II left ventricular assist system: from concept to first clinical use. The Annals of Thoracic Surgery, Vol. 71, 2001, p. 1620.

Schweizer G., Maslen E. H. Chapter 7: Dynamics of the Rigid Rotorsa; Chapter 10: Dynamics of Flexible Rotors. Magnetic Bearings. Theory, Design and Application to Rotating Machinery, Springer, Berlin Heidelberg, 2009, p. 167189, p. 251297.

Franklin G. F., Powell J. D., EmamiNaeini A. Feedback Control of Dynamic Systems. 4th Edition. Prentice Hall, Upper Saddle River, NJ, 2002.

Barbaraci G., Virzì Mariotti G. Suboptimal control law for active magnetic bearings suspension. Journal of Control Engineering and Technology, Vol. 2, Issue 1, 2012, p. 110.
About this article
This work was supported by the Russian Foundation for Basic Research (Grant No. 152901085 ofi_m).