1 Introduction

With the development of industrial technology, it is difficult for traditional mechanical bearings to meet the requirements of high-performance applications. However, due to the significant advantages of low friction, no wear and high speed, magnetic bearings have played an important role in advanced equipment fields such as aero-engines, large-scale mechanical equipment and marine gas turbines [1,2,3]. Therefore, the corresponding stability and safety problems for magnetic bearing system have gradually been widely concerned [4, 5].

In order to better suppress the vibration of magnetic bearing and improve the stability of magnetic bearing operation, the control algorithm needs to be carefully designed. Liu and his colleagues [6] solved the unbalanced vibration problem of active magnetic bearing rigid rotor system by automatic balancing method. Based on the traditional influence coefficient method for controlling the vibration of a flexible rotor, the optimal influence coefficient control method with weights for multinode unbalanced vibration of flexible electric spindle rotors was proposed by Qiao and others [7]. Marx and Nataraj [8] combined PD feedback and optimal feedforward control to effectively diminish the influence of pedestal simple harmonic motion on magnetic bearing system. Theodore et al. [9] designed a PID controller for linear magnetic bearing model, and tested the performance of the proposed PID controller according to the existing experimental data. Adaptive suppression for unbalanced vibration of magnetic suspension control components was analyzed by Xu and the team [10]. Miroslav and the co-workers [11] introduced the realization of magnetic bearing and its active control system based on PI and PID controllers. Wu et al. [12] designed a backstepping controller to improve the dynamical performance and robustness of three-pole radial hybrid magnetic bearing system. According to the characteristic equation and Routh-Hurwitz criterion, Sun et al. [13] proposed a PID parameter adjustment strategy of radial active magnetic bearing based on dynamical stiffness. In order to improve the stability and performance of active magnetic bearing system, the fractional PID controller was introduced by Jonathan and Pabitra [14].

Magnetic bearing system has very rich dynamical characteristics based on high order, nonlinearity and coupled terms. Therefore, some scholars were attracted to study its dynamical characteristics, in which the mathematical model of nonlinear magnetic bearing system was established [15, 16]. Inayat [17] reported the numerical investigation of magnetic bearing model and considered the nonlinear dynamical characteristics caused by the geometric coupling of magnetic actuators. Jiang et al. [18] conducted bifurcation analysis of active magnetic bearing system with flexible rotor, and explored the main factors affecting the dynamical characteristics of magnetic bearing system. Hou et al. [19] combined harmonic balance method, alternating frequency domain and time domain method to identify the bifurcation behavior of the rub-impact rotor system. Liu [20] established a new dynamical model to analyse the rotor bearing system with local fault, which included additional excitation area. The stability of rotor bearing system under the influences of deviation and parameter uncertainty was analyzed by Sun et al. [21]. Lin and the copartners [22] studied the nonlinear vibration characteristics and stability of double-disk rotor bearing system under multi-frequency excitation, in which the experimental results showed that the change of system parameters will lead to saddle node and Hopf bifurcation. Aldemir and the cooperators [23] analyzed the dynamical behaviors of a flexible rotor supported by hydrodynamic bearings with uncertain parameters. Based on the state feedback control and the injection of dither signals, Zhang [24] explored the nonlinear dynamics and chaos suppression for the magnetic bearing system. Sun and the collaborators [25] analyzed the nonlinear dynamical characteristics of active magnetic bearing system based on the element mapping method.

In fact, the magnetic bearing system is inevitably affected by various stochastic factors during its operation. Therefore, it is more practical to study the nonlinear stochastic magnetic bearing model. Leng et al. [26] used Monte Carlo numerical simulation to simulate the bifurcation and chaos characteristics of cracked rotor with white noise disturbance. Wang et al. [27] discussed Hopf bifurcation of magnetic bearing system with stochastic parameter disturbance, and proposed a hybrid feedback control strategy for the dynamical behavior. Do and the team [28] provided a Lyapunov-type theorem to study the fitness and exponential stability of dynamical systems driven by \(L\acute{e}vy\) process, and applied this theorem to the design of tracking controllers to realize the exponential stability of magnetic bearings. A multi-scale method was used by Ma et al. [29] to study the vibration characteristics of magnetic bearing system excited by narrow-band noise. The bifurcation characteristics of the active magnetic bearing-rotor system subjected to the external excitation were investigated analytically by Zhang et al. [30]. Yao et al. [31] discussed the nonlinear dynamics of the rotor sealing system with stochastic delay with elastic support.

At present, some scholars have done enough work on the deterministic magnetic bearing system, which was to design appropriate controllers to ensure the stability of the system. Others have studied the dynamical behaviors of magnetic bearing system influenced by stochastic factors. However, the researches about magnetic bearing system influenced by internal and external noise were unusually found which is more important for designing and controlling effectively from global operating system. Therefore, this paper takes this as the breakthrough point to explore the dynamical behavior evolution and control of magnetic bearing system under the influence of multisource stochastic factors. Briefly, the innovation of this paper is mainly reflected in the following:

(1) Considering the joint influence of additive and multiplicative noise on the system, we construct the mathematical model of magnetic bearing system under the multisource stochastic factors.

(2) The fractional PID controller is introduced, which has better effect on the stable operation of the system than the PID controller.

The rest of the paper is organized as follows: in Sect. 2, we establish the mathematical model of magnetic bearing system and analyze the related theories. Afterwards, we explore the influence of parameters on the system in Sect. 3. Then, the controlled magnetic bearing system is studied and simulated numerically in Sect. 4. In the end, a brief conclusion is given in Sect. 5.

2 Theoretical Analysis

There are two types of magnetic bearing-rotor system: rigid rotor system and flexible rotor system. Due to the high rigidity of rigid rotor system, it is assumed that the rotor system will not deform during the rotation in this paper. According to the actual situation, the model of magnetic bearing-rigid rotor system affected by multisource stochastic factors is as follows [32]:

$$\begin{aligned} \ddot{x} +2\mu \dot{x} +\omega ^{2}x-(\alpha _{1}x^{3}+\alpha _{2}x^{2}\dot{x}+\alpha _{3}x\dot{x}^{2})=W_{1}(t)+\dot{x}W_{2}(t), \end{aligned}$$
(1)

where x, \(\mu\), and \(\omega\) respectively represent the displacement, magnetic flux rate and natural frequency of the rotor in direction x. Additive noise \(W_{1}(t)\) and multiplicative noise \(W_{2}(t)\) are independent Gaussian white noises, and their power spectral density and autocorrelation function are

$$\begin{aligned} \Phi _{W_{i}W_{i}}(\omega _{0} ) =k_{ii}, R_{W_{i}W_{i}}(\tau )=2\pi k_{ii}\delta (\tau ), i=1, 2, \end{aligned}$$
(2)

where \(k_{ii}\) is the Gaussian white noise intensity, and \(\delta (\tau )\) is the Dirac function.

Theorem 1

The joint PDF of magnetic bearing-rigid rotor Eq. (1) driven by multiplicative noise and additive noise is shown as:

$$\begin{aligned} p(x, \dot{x})=\frac{2 \omega C}{\pi ^{2}} \exp \left( \frac{\left( x^{2}+\frac{\dot{x}^{2}}{\omega ^{2}}\right) \alpha _{2}}{6 k_{22} \pi }\right) \left( 3 k_{22} x^{2} \omega ^{2}+3 k_{22} \dot{x}^{2}+4 k_{11}\right) ^{D-1}, \end{aligned}$$
(3)

where C is a constant, and

$$\begin{aligned} D=\frac{-3 \pi k_{22} ^{2} \omega ^{2}+12 \mu k_{22} \omega ^{2}+2 k_{11} \alpha _{2}}{9 k_{22} ^{2} \omega ^{2} \pi } \end{aligned}$$

Proof

Firstly, use the amplitude envelope stochastic average method to solve the Eq. (1). And then the average formula is used to average the time on its quasi-period. Finally, the Van der pol transformation is performed on the magnetic bearing-rigid rotor system as follows:

$$\begin{aligned} x=A(t)cos\theta, \dot{x}=-A(t)\omega sin\theta, \theta =\omega t+\varphi (t), \end{aligned}$$
(4)

where A(t) is the amplitude process and \(\varphi (t)\) is the phase process. Simultaneous Eq. (4) and Eq. (1),

$$\begin{aligned} \begin{aligned} \dot{A}&=-A\omega \sin \theta \cos \theta -\frac{1}{\omega }\ddot{x}\sin \theta,\\ \dot{\varphi }&=-\omega \cos ^{2}\theta -\frac{1}{\omega A}\ddot{x}\cos \theta. \end{aligned} \end{aligned}$$
(5)

From Eq. (1), we get

$$\begin{aligned} \begin{aligned} \ddot{x}&=-\omega ^{2}A\cos \theta +2 \mu A\omega \sin \theta +\alpha _{1}A^{3}\cos ^{3}\theta -\alpha _{2}A^{3}\cos ^{2}\theta \omega \sin \theta \\&\quad +\alpha _{3}A^{3}\cos \theta \omega ^{2}\sin ^{2}\theta +W_{1}(t)-A\omega \sin \theta W_{2}(t). \end{aligned} \end{aligned}$$
(6)

Substituting Eq. (6) into Eq. (5), we have

$$\begin{aligned} & \begin{aligned} \dot{A}&=\frac{\sin \theta }{\omega }[-2\mu A\omega \sin \theta -\alpha _1A^3\cos ^3\theta +\alpha _2A^3\cos ^2\theta \omega \sin \theta \\&\quad -\alpha _3A^3\cos \theta \omega ^2\sin ^2\theta ]+\frac{\sin \theta }{\omega }[-W_1(t)+A\omega \sin \theta W_2(t)], \\ \end{aligned} \end{aligned}$$
(7)
$$\begin{aligned} & \begin{aligned} \dot{\varphi }&= \frac{\cos \theta }{\omega A}[-2\mu A\omega \sin \theta -\alpha _1A^3\cos ^3\theta +\alpha _2A^3\cos ^2\theta \omega \sin \theta \\&\quad -\alpha _3A^3\cos \theta \omega ^2\sin ^2\theta ]+\frac{\cos \theta }{\omega A}[-W_1(t)+A\omega \sin \theta W_2(t)]. \end{aligned} \end{aligned}$$
(8)

Then according to the standard stochastic average method and Eq. (8), we get

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{A}(t)=f_{1}(A, \varphi )+g_{11} W_{1}(t)+g_{12} W_{2}(t) \\ \dot{\varphi }=f_{2}(A, \varphi )+g_{21} W_{1}(t)+g_{22} W_{2}(t), \end{array}\right. } \end{aligned}$$
(9)

where

$$\begin{aligned} & \begin{aligned} f_1(A,\varphi )&=\frac{\sin \theta }{\omega }[-2\mu A\omega \sin \theta -\alpha _{1}A^{3}\cos ^{3}\theta \\&\quad +\alpha _2A^3\cos ^2\theta \omega \sin \theta -\alpha _3A^3\cos \theta \omega ^2\sin ^2\theta ], \\ \end{aligned} \end{aligned}$$
(10)
$$\begin{aligned} & \begin{aligned} f_{2}(A,\varphi )&= \frac{\cos \theta }{\omega A}[-2\mu A\omega \sin \theta -\alpha _1A^3\cos ^3\theta \\&\quad +\alpha _{2}A^{3}\cos ^{2}\theta \omega \sin \theta -\alpha _{3}A^{3}\cos \theta \omega ^{2}\sin ^{2}\theta ], \\ \end{aligned} \end{aligned}$$
(11)
$$\begin{aligned} & \begin{aligned} g_{11}&=-\frac{\sin \theta }{\omega },g_{12}=A\sin ^{2}\theta, \end{aligned} \end{aligned}$$
(12)
$$\begin{aligned} & \begin{aligned} g_{21}&=-\frac{\cos \theta }{\omega A},g_{22}=\cos \theta \sin \theta, \end{aligned} \end{aligned}$$
(13)

A(t) and \(\varphi (t)\) constitute a Markov diffusion process, and the system is quasi-periodic motion with a period of \(2\pi /\omega =T\). Therefore, the following time averages are performed on the quasi-period for Eq. (10) to Eq. (13)

$$\begin{aligned} \left\langle [\cdot ]\right\rangle _{t}=\frac{1}{T}\int _{0}^{T}\left[ \cdot \right] \textrm{d}t \end{aligned}$$

When solving Markov diffusion process, Ito equation about amplitude is obtained by Ito differential formula. Note that the equation of average amplitude A(t) does not contain phase \(\varphi (t)\). Consequently, the smooth Ito equation of amplitude is

$$\begin{aligned} \textrm{d}A=m(A)\textrm{d}t+\sigma (A)\textrm{d}B(t). \end{aligned}$$
(14)

The drift coefficient m(A) and diffusion coefficient \(\sigma (A)\) are obtained by the following time average formula, namely

$$\begin{aligned} m_{j}(X)&=\left\langle f_{j}(X_{t},t)\right\rangle _{t}+\sum _{l,s=1}^{m}\sum _{r=1}^{n_{1}}\int _{-\infty }^{0}\left\langle g_{rs}(X_{t+\tau },t+\tau )\frac{\partial }{\partial X_{r}}g_{jl}(X_{t},t)\right\rangle _{t}R_{ls}(\tau )\textrm{d}\tau, \end{aligned}$$
(15)
$$\begin{aligned} (\sigma \sigma ^{T})_{jk}&=\sum _{r=1}^{n_{1}}\sigma _{jr}(X)\sigma _{kr}(X)=\sum _{l,s=1}^{m}\int _{-\infty }^{\infty }\left\langle g_{jl}(X_{t},t)g_{ks}(X_{t+\tau },t+\tau )\right\rangle _{t}R_{ls}(\tau )\textrm{d}\tau. \end{aligned}$$
(16)

Through Eq. (15) and Eq. (16), we obtain

$$\begin{aligned} m(A)&=-\mu A+\frac{\alpha _{2}A^{3}}{8}+\frac{5\pi k_{22}A}{8}+\frac{\pi k_{11}}{2\omega ^{2}A}, \end{aligned}$$
(17)
$$\begin{aligned} \sigma ^{2}(A)&=\frac{3\pi k_{22}A^{2}}{4}+\frac{\pi k_{11}}{\omega ^{2}}. \end{aligned}$$
(18)

The corresponding FPK equation is

$$\begin{aligned} \frac{\partial }{\partial t}p+\frac{\partial }{\partial t}[m(A)p]-\frac{1}{2}\frac{\partial ^{2}}{\partial A^{2}}[\sigma ^{2}(A)p]=0. \end{aligned}$$

In order to solve the FPK equation, firstly, we need to solve the stationary PDF p(A). And p(A) satisfies the following simplified FPK equation

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}A}G=\frac{\textrm{d}}{\textrm{d}A}\{m(A)p(A)-\frac{1}{2}\frac{\textrm{d}}{\textrm{d}A}[\sigma ^{2}(A)p(A)]\}=0. \end{aligned}$$
(19)

By integrating Eq. (19), we get

$$\begin{aligned} G(A)=m(A)p(A)-\frac{1}{2}\frac{\textrm{d}}{\textrm{d}A}[\sigma ^{2}(A)p(A)]=G_{c}, \end{aligned}$$
(20)

As can be seen from Eq. (20), the probability G(A) is a constant. To ensure the existence of p(A), the most common case for the boundary is to take its probability flow as zero. Hence, from Eq. (20), we get a simplified solution

$$\begin{aligned} \begin{aligned} p(A)&=\frac{C}{\sigma ^{2}(A)}\exp [\int \frac{2m(A)}{\sigma ^{2}(A)}\textrm{d}A] \\&=\frac{4\omega ^{2}C}{3\pi k_{22}A^{2}\omega ^{2}+4\pi k_{11}}\exp \left[ \ln (A)+\frac{A^{2}\alpha _{2}}{6k_{22}\pi } \right. \\&\quad \left. -\frac{\ln (3k_{22}A^2\omega ^2+4k_{11})\left( -3\pi k_{22} ^2\omega ^2+12\mu k_{22}\omega ^2+2k_{11}\alpha _2\right) }{9k_{22} ^2\omega ^2\pi }\right]. \end{aligned} \end{aligned}$$
(21)

Then, by using the transformation of Eq. (4), the stationary response PDF about x and \(\dot{x}\) is obtained as follows:

$$\begin{aligned} p(x,\dot{x})=\frac{1}{2\pi \omega A}p(A), A=\sqrt{x^2+\frac{\dot{x}^2}{\omega ^2}}, \end{aligned}$$
(22)

Substitute Eq. (22) into Eq. (21), and the proof is completed. \(\square\)

3 Parameter Influence Analysis

Based on the stationary response PDF, shown in Eq. (22), the corresponding joint PDF and its cross-section at \(\dot{x} =0\) are obtained. And by controlling variables, we explore the effects of additive noise, multiplicative noise and magnetic permeability for the magnetic bearing-rigid rotor system.

Fig. 1
figure 1

Joint PDF of the system under different additive noise intensities and the cross-section comparison at \(\dot{x} =0\). (a) \(k_{11}=0.02\), (b) \(k_{11}=0.1\), (c) \(k_{11}=0.5\), (d) the cross-section comparison at \(\dot{x} =0\)

Full size image

Firstly, we analyze the influences of additive noise on the magnetic bearing-rigid rotor system. Let \(\alpha _{2}=-1.8591\), \(\omega =1\), \(\mu =1\), \(k_{22}=0.2\), the additive noise intensity \(k_{11}\) is chosen as 0.02, 0.1 and 0.5 respectively. We draw its joint stationary PDF and cross-section at \(\dot{x}=0\), as shown in Fig. 1.

In Fig. 1a, we can clearly see that when the additive noise intensity is 0.02, the PDF is crater-shaped. When the additive noise intensity is 0.1, the crater shape is not obvious, as shown in Fig. 1b. Besides with the increase of the additive noise intensity, it is not difficult to see that the crater of the PDF gradually decreases until it disappears and becomes a single peak. In addition, it can be clearly observed from Fig. 1d that with the increase of additive noise intensity, the peak value of the PDF go through the firstly increasing and then gradually decreasing. Correspondingly, the number of peaks gradually decreases from the classic double peaks until it becomes a single peak.

Fig. 2
figure 2

Joint PDF of the system under different multiplicative noise intensities and the cross-section comparison at \(\dot{x} =0\). (a) \(k_{22}=0.05\), (b) \(k_{22}=0.1\), (c) \(k_{22}=0.2\), (d) the cross-section comparison at \(\dot{x} =0\)

Full size image

Secondly, the influences of multiplicative noise intensity on the magnetic bearing-rigid rotor system are expounded. Fixed the additive noise intensity and the magnetic permeability intensity, let \(\alpha _{2}=-1.8591\), \(\omega =1\), \(\mu =1\), \(k_{11}=0.3\), the multiplicative noise intensity \(k_{22}\) is chosen as 0.02, 0.1 and 0.5 respectively. We draw the joint stationary PDF and the cross-section at \(\dot{x}=0\). The joint stationary PDF of the magnetic bearing system and the cross-section comparison at \(\dot{x} =0\) are displayed in Fig. 2. When the multiplicative noise intensity is selected as 0.05, 0.1 and 0.2 respectively, the joint PDF always presents a single peak. However, when the multiplicative noise increases from 0.05 to 0.1, the PDF increases significantly, while when the multiplicative noise increases from 0.1 to 0.2, the PDF increases slowly. Figure 2d also shows that with the increase of multiplicative noise intensity, the single peak height of the PDF is gradually raising, and the amplitude in the displacement direction of the system is increasing as well.

Fig. 3
figure 3

Joint PDF of the system under different permeability and the cross-section comparison at \(\dot{x} =0\). (a) \(\mu =0.1\), (b) \(\mu =0.5\), (c) \(\mu =1\), (d) the cross-section comparison at \(\dot{x} =0\)

Full size image

Finally, we discuss the influences of magnetic permeability on the magnetic bearing-rigid rotor system. Let \(\alpha _{2} =-1.8591\), \(\omega =1\), \(k_{11} =0.3\), \(k_{22} =0.3\), the magnetic permeability intensity \(\mu\) is chosen as 0.1, 0.5 and 1 respectively. Its joint stationary PDF and cross-section at \(\dot{x}=0\) are shown in Fig. 3.

As magnetic permeability intensity \(\mu\) increases from 0.1 to 1, the PDF always presents a single peak from the Fig. 3. And with the obvious increase of magnetic permeability, the peak height does not increase obviously, as shown in Fig. 3a, b and c. It can also be clearly seen from Fig. 3d that with the increase of magnetic permeability, the single peak height of the PDF is slowly increasing. In brief, Fig. 3d obviously shows that the increase of magnetic permeability increases the amplitude in the displacement direction and the vibration frequency.

4 Theoretical Analysis

Assume that the corresponding parameters of the magnetic bearing rigid rotor system are as follows:

$$\begin{aligned} \begin{aligned} \alpha _{1} = -0.9262, \alpha _{2} = -1.8591, \alpha _{3} = -0.0113, \mu = 1, \omega = 1, k_{11} = 0.3, k_{22} = 0.4, \end{aligned} \end{aligned}$$
Fig. 4
figure 4

Velocity, displacement time history diagram and phase diagram of a sample solution of initial system

Full size image

A sample solution of magnetic bearing system is randomly selected for numerical simulation. We obtain the phase diagram of the system, the time history of displacement and velocity, as shown in Fig. 4a, b and c, in which the initial value of the system is \((x,\dot{x} )=(0,0)\).

Figure 4a shows that the system presents an irregular motion state. Furthermore, it can be seen that the velocity oscillation is fierce, while the displacement fluctuation is relatively stable, as shown in Fig. 4b and c. Compared with the traditional controller, the fractional PID controller adds two adjustable parameters, which makes the design more flexible. In order to make the magnetic bearing system run more smoothly under the influence of stochastic factors, the fractional PID control is carried out for the system.

4.1 Controlled Magnetic Bearing-Rigid Rotor System

The equation of magnetic bearing-rigid rotor system with multisource stochastic noise and fractional \(PI^{\lambda } D^{\eta }\) controller is as follows:

$$\begin{aligned} \ddot{x}+2 \mu \dot{x}+\omega ^{2} x-\left( \alpha _{1} x^{3}+\alpha _{2} x^{2} \dot{x}+\alpha _{3} x \dot{x}^{2}\right) =\varepsilon u(x, \dot{x})+W_{1}(t)+\dot{x} W_{2}(t), \end{aligned}$$
(23)

where \(\mu\), \(\omega\), \(\alpha _{1}\), \(\alpha _{2}\), \(\alpha _{3}\), \(W_{1}(t)\) and \(W_{2}(t)\) are expressed similarly to the above Eq. (1). \(\varepsilon u(x, \dot{x})=\varepsilon k_{0} x+\varepsilon k_{1} I^{\lambda } x+\varepsilon k_{2} D^{\eta } x\) is a fractional \(PI^{\lambda } D^{\eta }\) controller, \(k_{0}\), \(k_{1}\) and \(k_{2}\) represent the corresponding control intensity, \(I^{\lambda } x\) and \(D^{\eta } x\) denote \(\lambda\)-order fractional integral and \(\eta\)-order fractional derivative defined by Riemann-Liouville, and \(\varepsilon\) represents the coefficient of the controller, satisfying \(0<\varepsilon <1\).

In the definition of Riemann-Liouville, let \(n=1\), \(t-\tau =s\), the fractional integral of \(\lambda\) order and the fractional derivative of \(\eta\) order can be expressed as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} D^{\eta }x(t)=\frac{1}{\Gamma (1-\eta )}\frac{\textrm{d}}{\textrm{d}t}\int _{0}^{t}\frac{x(t-s)}{s^{\eta }} \textrm{d}s, 0<\eta<1\\ I^{\lambda }x(t)=\frac{1}{\Gamma (\lambda )}\int _{0}^{t}\frac{x(t-s)}{s^{1-\lambda }}\textrm{d}s, 0<\lambda <1. \end{array}\right. } \end{aligned}$$

The Eq. (23) disturbed by weak noise and light damping is a quasi-conservative system, so the energy consumed by system damping is quite much as that generated by noise excitation, i.e., the response of Eq. (23) is quasi-periodic. Using the definition of Gamma function, \(\Gamma (\eta )=\int _{0}^{\infty }e^{-z}z^{\eta -1}dz\) and its property \(\Gamma (\eta )\Gamma (1-\eta )=\pi /(\sin \pi \eta )\), the fractional derivative can be expressed as the following double integral

$$\begin{aligned} D^{\eta }x(t)=\int _{0}^{\infty }\int _{0}^{t}\frac{\sin \pi \eta }{\pi }e^{-z}z^{\eta -1}\frac{\dot{x}(s)}{(t-s)^{\eta }} \textrm{d}s\textrm{d}z. \end{aligned}$$

In the above formula, let \(z=(t-s)y^{1/\eta }\) and define \(\psi (y,t)=(\sin \pi \eta /\pi \eta )\int _{0}^{t}e^{-(t-s)y^{1/\eta } }\dot{x}(s) \textrm{ds}\), thus the fractional derivative can be further converted into improper integral with \(\psi (y,t)\) as follows:

$$\begin{aligned} D^{\eta }x(t)=\int _{0}^{\infty }\int _{0}^{t}\frac{\sin \pi \eta }{\pi \eta }e^{-(t-s)y^{1/\eta }}\dot{x}(s)\textrm{d}s\textrm{d}y=\int _{0}^{\infty }\psi (y,t)\textrm{d}y. \end{aligned}$$
(24)

According to Eq. (24), take the partial derivative of a about time t based on Leibniz integral rule, and we have

$$\begin{aligned} \begin{aligned} \frac{\partial \psi (y,t)}{\partial t}=&\frac{\sin \pi \eta }{\pi \eta }e^{-(t-t)y^{\frac{1}{\eta }}}\dot{x}(t)\cdot 1 -\frac{\sin \pi \eta }{\pi \eta }e^{-(t-0)y^{\frac{1}{\eta }}}\dot{x}(0)\cdot 0 \\&+\int _{0}^{t}\frac{\partial }{\partial t} \left( \frac{\sin \pi \eta }{\pi \eta }e^{-(t-s)y^{\frac{1}{\eta }}}\dot{x}(s)\right) \textrm{d}s. \end{aligned} \end{aligned}$$

Simplified, we get

$$\begin{aligned} \frac{\partial \psi (y,t)}{\partial t}=\frac{\sin \pi \eta }{\pi \eta }\dot{x}(t)+\int _{0}^{t}\frac{\sin \pi \eta }{\pi \eta }\left( -y^{\frac{1}{\eta }}e^{-(t-s)y^{\frac{1}{\eta }}}\right) \dot{x}(s)\textrm{d}s, \end{aligned}$$

that is

$$\begin{aligned} \frac{\partial \psi (y,t)}{\partial t}=\frac{\sin \pi \eta }{\pi \eta }\dot{x}(t)-y^{\frac{1}{\eta }}\int _{0}^{t}\frac{\sin \pi \eta }{\pi \eta }e^{-(t-s)y^{\frac{1}{\eta }}}\dot{x}(s)\textrm{d}s. \end{aligned}$$

Noting that the integral term is the definition of \(\psi (y,t)\), so \(\psi (y,t)\) satisfies the following partial differential equation

$$\begin{aligned} \frac{\partial \psi (y,t)}{\partial t}=-y^{\frac{1}{\eta }}\psi (y,t)+\frac{\sin \pi \eta }{\pi \eta }\dot{x}(t). \end{aligned}$$
(25)

Eq. (25) is a linear inhomogeneous ordinary differential equation. With the help of constant variation method and Van der Pol transformation, we have

$$\begin{aligned} \psi (y,t)=\frac{\sin \pi \eta }{\pi \eta }A(t)\frac{\omega ^{2}\cos \theta (t)}{\omega ^{2}+y^{2/\eta }}-\frac{\sin \pi \eta }{2}A(t)\frac{\omega \sin ^{2}\theta (t)y^{1/\eta }}{\omega ^{2}+y^{1/\eta }}+ce^{ty^{-1/\eta }}, \end{aligned}$$
(26)

where c is derived from the initial condition \(\psi (y,0)=0\). Then substitute Eq. (26) into the Eq. (25), we get the slow variable expression of \(\eta\) order fractional derivative in the fractional controller as follows:

$$\begin{aligned} D^{\eta }x(t)=A\omega ^{\eta }\left( \cos \theta \cos \frac{\pi \eta }{2}-\sin \theta \sin \frac{\pi \eta }{2}\right), 0<\eta <1. \end{aligned}$$
(27)

According to the relationship between fractional derivative and fractional integral, i.e.

$$\begin{aligned} D^{\eta }x(t)=D^{(n)}I^{n-\eta }x(t)=I^{n-\eta }D^{(n)}x(t). \end{aligned}$$

As \(n=1\), fractional integral can be expressed as

$$\begin{aligned} D^{\eta }x(t)=D^{(1)}I^{1-\eta }x(t)=\frac{\textrm{d}}{\textrm{d}t}I^{1-\eta }x(t). \end{aligned}$$

Using the method of variable separation, it is obtained

$$\begin{aligned} I^{1-\eta }x(t)=A\omega ^{\eta -1}\left( \sin \theta \cos \frac{\pi \eta }{2}+\cos \theta \sin \frac{\pi \eta }{2}\right), \end{aligned}$$

Let\(1-\eta =\lambda\), we have

$$\begin{aligned} I^{\lambda }x(t)=A\omega ^{-\lambda } \left( \sin \theta \cos \frac{\pi (1-\lambda )}{2}+\cos \theta \sin \frac{\pi (1-\lambda )}{2}\right), 0<\lambda <1. \end{aligned}$$
(28)

By substituting Eq. (28) and Eq. (27) into the expression of fractional-order controller, we get

$$\begin{aligned} \begin{aligned} u(x,\dot{x})&=k_{0}x+k_{1}I^{\lambda }x+k_{2}D^{\eta }x \\&=k_{0}A\cos \theta +k_{1}\left[ a\omega ^{-\lambda }\left( \sin \theta \cos \frac{\pi (1-\lambda )}{2}+\cos \theta \sin \frac{\pi (1-\lambda )}{2}\right) \right] \\&\quad +k_2\left[ A\omega ^\eta \left( \cos \theta \cos \frac{\pi \eta }{2}-\sin \theta \sin \frac{\pi \eta }{2}\right) \right] \\&=\left( k_{0}+k_{1}\omega ^{-\lambda }\sin \frac{\pi (1-\lambda )}{2}+k_{2}\omega ^{\eta }\cos \frac{\pi \eta }{2}\right) A\cos \theta \\&+\left( k_1\omega ^{-\lambda -1}\cos \frac{\pi (1-\lambda )}{2}-k_2\omega ^{\eta -1}\sin \frac{\pi \eta }{2}\right) A\omega \sin \theta. \end{aligned} \end{aligned}$$

To facilitate the calculation, we write down the formula as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} m(\eta,\lambda )=k_{0}+k_{1}\omega ^{-\lambda } \sin \frac{\pi (1-\lambda )}{2}+k_{2}\omega ^{\eta }\cos \frac{\pi \eta }{2}\\ n(\eta,\lambda )=k_{1}\omega ^{-\lambda -1}\cos \frac{\pi (1-\lambda )}{2}-k_{2}\omega ^{\eta -1}\sin \frac{\pi \eta }{2}. \end{array}\right. } \end{aligned}$$

Then the fractional expression is

$$\begin{aligned} u(x,\dot{x})=m(\eta,\lambda )A\cos \theta +n(\eta,\lambda )A\omega \sin \theta. \end{aligned}$$
(29)

By introducing Eq. (29) into Eq. (23) and Van der Pol transformation, we obtain the controlled system as follows:

$$\begin{aligned} \begin{aligned}&\ddot{x}+2\mu \dot{x}+\omega ^{2}x-(\alpha _{1}x^{3}+\alpha _{2}x^{2}\dot{x} +\alpha _{3}x\dot{x}^{2})=\varepsilon x\left( k_{0}+k_{1}\omega ^{-\lambda }\sin \frac{\pi (1-\lambda )}{2}\right. \\&\quad \left. +k_{2}\omega ^{\eta }\cos \frac{\pi \eta }{2}\right) +\varepsilon \dot{x}\left( k_{1} \omega ^{-\lambda -1}\cos \frac{\pi (1-\lambda )}{2}-k_{2}\omega ^{\eta -1}\sin \frac{\pi \eta }{2}\right) \\&\quad +W_{1}(t)+\dot{x}W_{2}(t), \end{aligned} \end{aligned}$$
(30)

Finally, it is necessary to prove the equivalence of Eq. (23)and Eq. (30).

Fig. 5
figure 5

Response curves of the system derived from theoretical method and numerical method respectively

Full size image

The system response curves obtained by the two methods are shown in Fig. 5. According to numerical analysis, we find that the time history trend of displacement and velocity is basically consistent, which also proves that Eq. (23) and Eq. (30) are equivalent.

4.2 Controlled Magnetic Bearing-Rigid Rotor System

Let \(k_{0}=0.34\), \(k_{1}=0.2\), \(k_{2}=0.56\), adjust the fractional fraction \(\lambda =0.64\), the integral number \(\eta =0.53\) and the controller coefficient \(\varepsilon =0.45\). Then we carry out numerical simulation on the Eq. (29), and the results are shown in Fig. 6a, b and c.

Fig. 6
figure 6

Comparison of velocity, displacement time history and phase diagram of sample solutions before and after control in system (1)

Full size image

Figure 6 shows the corresponding velocity time history diagram in Fig. 6a, system phase diagram in Fig. 5b and displacement time history diagram in Fig. 6c before and after the magnetic bearing-rigid rotor system is controlled respectively. It can be clearly seen from Fig. 6a that the velocity oscillation is intense before the control, but becomes relatively mild after the control. The phase diagram trajectory of the system presents an irregular state at first, but gradually converges to the equilibrium point after being controlled, as shown in Fig. 6b. Figure 6c shows that the fluctuation of displacement is more stable than before. According to Fig. 6, it can be concluded that the fractional PID controller can control the magnetic bearing-rigid rotor system as expected. In Fig. 7, the system is compared with the ordinary PID control. It can be observed that the velocity time history diagram, displacement time history diagram and phase diagram of the system sample solution controlled by fractional PID are more concentrated in the fluctuation near the origin than the ordinary PID controller, and the control effect is better.

Fig. 7
figure 7

Comparison of velocity, displacement time history and phase diagram of the sample solution of the system by the ordinary PID controller and the fractional PID controller

Full size image

4.3 Controller Parameter Analysis

Fixing other parameters, we will further explore the changes of the velocity time history diagram and phase diagram of the sample solution of the controlled Eq. (29) with different fractional fractions \(\eta\), integral numbers \(\lambda\) and controller coefficients \(\varepsilon\).

Fig. 8
figure 8

Velocity time history diagram and phase diagram of sample solution of controlled system as \(\varepsilon =0.35\), \(\varepsilon =0.5\), \(\varepsilon =0.75\)

Full size image
Fig. 9
figure 9

Velocity time history diagram and phase diagram of sample solution of controlled system as \(\lambda =0.1\), \(\lambda =0.3\), \(\lambda =0.8\)

Full size image

According to Figs. 8, 9 and 10 when other parameter values are fixed, the velocity time history diagram and phase diagram of the corresponding sample solution will converge to a certain range by changing the fractional fraction \(\eta\), integral number \(\lambda\) and controller coefficient \(\varepsilon\) of the fractional controller. In Fig. 8, with the constant reduction of the controller coefficient, the vibration amplitude of the velocity time history diagram of the controlled magnetic bearing-rigid rotor system decreases, and the phase diagram trajectory further converges. In addition, with the increasing of fractional fraction \(\eta\) and integral number \(\lambda\), the time history diagram of velocity oscillates in a small range, and the phase diagram trajectory converges more obviously, as shown in Figs. 9 and 10. Therefore, fractional PID control has a good effect on the safe and stable operation of magnetic bearing-rigid rotor system.

Fig. 10
figure 10

Velocity time history diagram and phase diagram of sample solution of controlled system as \(\eta =0.27\), \(\eta =0.6\), \(\eta =0.75\)

Full size image

5 Conclusion

In this paper, we analyze the probability evolution law and control of the response PDF for the magnetic bearing system affected by multisource stochastic factors in complex environment. Firstly, the model of magnetic bearing-rigid rotor system under the influence of additive and multiplicative noise is established. Secondly, the FPK equation of the system response is obtained by using the stochastic average method of amplitude envelope. Then, through numerical simulation, we verify the evolution law of response PDF. Finally, we introduce the fractional PID controller to analyze the controlled magnetic bearing-rotor system.

Through the above theoretical and numerical analyses, we find that with the increase of additive noise intensity, the amplitude in the displacement direction of the magnetic bearing system decreases, while with the increase of multiplicative noise intensity and permeability, the amplitude of the magnetic bearing system in the displacement direction will gradually increase. We draw the conclusions that the instability of the magnetic bearing rotor system is mainly caused by the uncertainty of internal parameters. Therefore, improving the stability of magnetic bearing itself is the key to the safe operation of the system. Based on this, the magnetic bearing system can be better controlled by properly selecting the parameters of the fractional controller.