Analysis and Control of Flywheel Energy Storage Systems

1. Introduction

Since a few years ago, electrical energy storage has been attractive as an effective use of electricity and coping with the momentary voltage drop. Above all, flywheel energy storage systems (FESS) using superconductor have advantages of long life, high energy density, and high efficiency (Subkhan & Komori, 2011), and is now considered as enabling technology for many applications, such as space satellites and hybrid electric vehicles (Samineni et al., 2006; Suvire & Mercado, 2012). Also, the contactless nature of magnetic bearings brings up low wear, absence of lubrication and mechanical maintenance, and wide range of work temperature (Bitterly, 1998; Beach & Christopher, 1998). Moreover, the closed-loop control of magnetic bearings enables active vibration suppression and on-line control of bearing stiffness (Cimuca et al., 2006; Park et al., 2008).

Active magnetic bearing is an open-loop unstable control problem. Therefore, an initial controller based on a rigid rotor model has to be introduced to levitate the rotor. In reality, the spinning rotor under the magnetic suspension may experience two kinds of whirl modes. The conical whirl mode gives rise to the gyroscopic forces to twist the rotor, thereby severely affecting stability of the rotor if not properly controlled (Okada et al, 1992; Williams et al., 1990). The translatory whirl mode constrains the rotor to synchronous motion in the radial direction so as to suppress the gyroscopic rotation, which has been extensively used in industry (Tomizuka et al, 1992; Tsao et al., 2000). The synchronization control has also been shown to be very capable in dealing with nonlinear uncertain models, and to be very effective in disturbance rejection for systems subject to synchronous motion. Until the advent of synchronization control, the prevalent use of the synchronization controller has been limited to stable mechanical systems and therefore is not readily applicable to magnetic systems which are unstable in nature and highly nonlinear (Yang & Chang, 1996).

In the past three decades the theory of optimal control has been well developed in nearly all aspects, such as stability, nonlinearity, and robustness (Summers et al., 2011; Rawlings et al., 2008; Mayne, et al., 2000). It is known that multivariable constrained control problems in state-space can be effectively handled using Linear Quadratic Gaussian (LQG). An application of the optimal control to synchronize multiple motion axes has been reported in (Zhu & Chen, 2001; Xiao & Zhu, 2006), where cross-coupling design of generalized predictive control was presented by compensating both the tracking error and the synchronous error. In this chapter, robust MPC control algorithms for the flywheel energy storage system with magnetically assisted bearings are developed. The controllers are derived through minimization of a modified cost function, in which the synchronization errors are embedded so as to reduce the synchronization errors in an optimal way.

2. Flywheel structure

Fig.1 illustrates the basic structure of a flywheel system with integrated magnetic bearings. The motor and generator with disk-type geometry are combined into a single electric machine, and the rotor is sandwiched between two stators. Each of the stators carries a set of three-phase copper winding to be fed with sinusoidal currents. Furthermore, both axial faces of the rotor contain rare-earth permanent magnets embedded beneath the surfaces. The radial magnetic bearing which consists of eight pairs of electromagnets is constructed around the circumference of hollow center. A combination of active and passive magnetic bearings allows the rotor to spin and remain in magnetic levitation.

The control of such a system normally includes two steps. First, the spinning speed and the axial displacement of the rotor are properly regulated (Zhang & Tseng, 2007). Second, a synchronization controller is introduced to suppress the gyroscopic rotation of the rotor caused by the outside disturbance and model uncertainty (Xiao et al., 2005).

Figure 1.

The flywheel energy storage system

3. System dynamics

Let $$x_c$$ and $$y_c$$ denote the displacements of the mass center of the rotor in the $$x$$ and $$y$$-directions, and $$\theta$$ and $$\phi$$ the roll angles of rotation about $$x$$-axis and $$y$$-axis, respectively. Note that $$\theta$$ and $$\phi$$ are assumed to be small since the air gap is very narrow within the magnetic bearings. It is also assumed that the rotor is rigid with its inertia perfectly balanced about the $$z$$-axis so that the flexibility and eccentricity of the rotor are not considered herein; thereby, the variation effects of tensor of inertia due to the roll motion of the rotor can be negligible.

The mass center of the rotor in the radial direction can be described by

$$\{\begin{array}{l}m{\ddot{x}}_c=f_{xa}+f_{xb}+e_x\\ m{\ddot{y}}_c=f_{ya}+f_{yb}+e_y\\ J_x\ddot{\theta }=f_{yb}{l}_b-f_{ya}{l}_a-J_{z}w\dot{\phi }+e_{\theta }\\ J_y\ddot{\phi }=f_{xa}{l}_a-f_{xb}{l}_b+J_{z}w\dot{\theta }+e_{\phi }\end{array}$$E1

where$$x_c=\frac{l_b}{l_a+l_b}{x}_a+\frac{l_a}{l_a+l_b}{x}_b$$，$$y_c=\frac{l_b}{l_a+l_b}{y}_a+\frac{l_a}{l_a+l_b}{y}_b$$, $$\theta =\frac{1}{l_a+l_b}\left(y_b-y_a\right)$$, $$\phi =\frac{1}{l_a+l_b}\left(x_a-x_b\right)$$, $$m$$ is the mass of the rotor, $$$$$$J_x$$, $$J_y$$and $$J_z$$ are the moments of inertia about $$x$$-axis, $$y$$-axis and $$z$$-axis respectively, $$w$$is the spinning rate about $$z$$-axis, $$f_{xa}$$, $$f_{xb}$$, $$f_{ya}$$and $$f_{yb}$$ are the magnetic forces along the radial directions, $$e_x$$, $$e_y$$, $$e_{\theta }$$and $$e_{\phi }$$ are the disturbances.

According to the Maxwell’s law, the magnetic forces$$f_{xa}$$, $$f_{xb}$$, $$f_{ya}$$and $$f_{yb}$$ have nonlinear relationships with the control currents and displacements of the rotor. Then, the magnetic forces at equilibriums can be linearized with Taylor’s method (Zhu et al., 2009),

$$\left[\begin{array}{c}{f}_{xa}\\ f_{xb}\\ f_{ya}\\ f_{yb}\end{array}\right]\approx K_p\left[\begin{array}{c}{x}_a\\ x_b\\ y_a\\ y_b\end{array}\right]+K_c\left[\begin{array}{c}{i}_{xa}\\ i_{xb}\\ i_{ya}\\ i_{yb}\end{array}\right]$$E2

where$$K_p=\frac{8G{I}_0^2}{h_0^3}$$, $$K_c=\frac{8G{I}_0\sin \beta }{h_0^2}$$, $$\beta$$ is a constant angle corresponding to the structure of electromagnets, $$I_0$$is the bias current, $$i_x$$and $$i_y$$ are the control currents near $$x$$-axis and $$y$$-axis, respectively, $$h_0$$is the nominal air gap at equilibrium, $$G$$is an electromagnet constant given by$$G=\frac{1}{4}{\mu }_0{A}_g{N}^2$$, $${\mu }_0$$is the air permeability, $$A_g$$is the cross-sectional area of air gap, and $$N$$ is the number of turns of the winding circuit.

Then, the state-space model of (1) is obtained,

$$\{\begin{array}{l}\dot{X}=A_{c}X+B_{c}u+B_x{d}_x\\ Z=C_{c}X+D_z{d}_z\end{array}$$E3

where $$X={\left[x_a,{\dot{x}}_a,x_b,{\dot{x}}_b,y_a,{\dot{y}}_a,y_b,{\dot{y}}_b\right]}^T$$is the state variable, $$u={\left[i_{xa},i_{xb},i_{ya},i_{yb}\right]}^T$$ is the forcing vector, $$Z={\left[x_a,x_b,y_a,y_b\right]}^T$$is the output vector, and $$T_1=T_2=T_3=T_4=\left[1,0\right]$$ are the output transition matrices, $$d_x$$and $$d_z$$ denote model uncertainties or system disturbances with appropriate matrices $$B_x$$ and$$D_z$$,

$$$$

$$A_c=\left[\begin{array}{cccccccc}0& 1& 0& 0& 0& 0& 0& 0\\ A_{c21}& 0& A_{c23}& 0& 0& A_{c26}& 0& A_{c28}\\ 0& 0& 0& 1& 0& 0& 0& 0\\ A_{c41}& 0& A_{c43}& 0& 0& A_{c46}& 0& A_{c48}\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& A_{c62}& 0& A_{c64}& A_{c65}& 0& A_{c67}& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\\ 0& A_{c82}& 0& A_{c84}& A_{c85}& 0& A_{c87}& 0\end{array}\right]$$E4

,$$B_c=\left[\begin{array}{cccc}0& 0& 0& 0\\ B_{c21}& B_{c22}& 0& 0\\ 0& 0& 0& 0\\ B_{c41}& B_{c42}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& B_{c63}& B_{c64}\\ 0& 0& 0& 0\\ 0& 0& B_{c83}& B_{c84}\end{array}\right]$$,

$$C_c=diag\left(T_1,T_2,\cdots ,T_4\right)$$E5

.

where

$$A_{c21}=\left(\frac{1}{m}+\frac{l_a^2}{J_y}\right)K_p$$, $$A_{c23}=\left(\frac{1}{m}-\frac{l_a{l}_b}{J_y}\right)K_p$$,

$$A_{c26}=-\frac{J_{z}w{l}_a}{J_y\left(l_a+l_b\right)}$$E6

, $$A_{c28}=\frac{J_{z}w{l}_a}{J_y\left(l_a+l_b\right)}$$, $$A_{c41}=\left(\frac{1}{m}-\frac{l_a{l}_b}{J_y}\right)K_p$$, $$A_{c43}=\left(\frac{1}{m}+\frac{l_b^2}{J_y}\right)K_p$$, $$A_{c46}=\frac{J_{z}w{l}_b}{J_y\left(l_a+l_b\right)}$$, $$A_{c48}=-\frac{J_{z}w{l}_b}{J_y\left(l_a+l_b\right)}$$, $$A_{c62}=\frac{J_{z}w{l}_a}{J_x\left(l_a+l_b\right)}$$, $$A_{c64}=-\frac{J_{z}w{l}_a}{J_x\left(l_a+l_b\right)}$$, $$A_{c65}=\left(\frac{1}{m}+\frac{l_a^2}{J_x}\right)K_p$$, $$A_{c67}=\left(\frac{1}{m}-\frac{l_a{l}_b}{J_x}\right)K_p$$, $$A_{c82}=-\frac{J_{z}w{l}_b}{J_x\left(l_a+l_b\right)}$$, $$A_{c84}=\frac{J_{z}w{l}_b}{J_x\left(l_a+l_b\right)}$$, $$A_{c85}=\left(\frac{1}{m}-\frac{l_a{l}_b}{J_x}\right)K_p$$, $$A_{c87}=\left(\frac{1}{m}+\frac{l_b^2}{J_x}\right)K_p$$, $$B_{c21}=\left(\frac{1}{m}+\frac{l_a^2}{J_y}\right)K_c$$, $$B_{c22}=\left(\frac{1}{m}-\frac{l_a{l}_b}{J_y}\right)K_c$$, $$B_{c41}=\left(\frac{1}{m}-\frac{l_a{l}_b}{J_y}\right)K_c$$, $$B_{c42}=\left(\frac{1}{m}+\frac{l_b^2}{J_y}\right)K_c$$, $$B_{c63}=\left(\frac{1}{m}+\frac{l_a^2}{J_x}\right)K_c$$, $$B_{c64}=\left(\frac{1}{m}-\frac{l_a{l}_b}{J_x}\right)K_c$$, $$B_{c83}=\left(\frac{1}{m}-\frac{l_a{l}_b}{J_x}\right)K_c$$,$$B_{c84}=\left(\frac{1}{m}+\frac{l_b^2}{J_x}\right)K_c$$.

During a closed-loop control phase, the position and rate of the shaft are constantly monitored by contactless sensors, and are processed in a controller, so that a control current to the coils of electromagnets which attract or repel the shaft is amplified and fed back.

4. Controller design

Let he discrete-time model of (3) be described by

$$\{\begin{array}{l}X\left(k+1\right)=AX\left(k\right)+Bu\left(k\right)\\ Z\left(k\right)=CX\left(k\right)\end{array}$$E7

where $$k$$ denotes the discrete time. Note that the disturbance term is ignored.

By introducing the following synchronization errors,

$$\delta \left(k\right)=LCX\left(k\right)$$E8

where

$$\delta \left(k\right)=\left[\begin{array}{c}{x}_a-x_b\\ x_b-y_a\\ y_a-y_b\\ y_b-x_a\end{array}\right],L=\left[\begin{array}{cccc}1& -1& 0& 0\\ 0& 1& -1& 0\\ 0& 0& 1& -1\\ -1& 0& 0& 1\end{array}\right],$$E9

it has the modified cost function,

$$\begin{array}{l}J\left(k\right)={\displaystyle \sum _{i=1}^{H_P}{\widehat{Z}}^T\left(k+i|k\right)\widehat{Z}\left(k+i|k\right)+\lambda {\displaystyle \sum _{i=1}^{H_c}{\widehat{u}}^T\left(k+i-1|k\right)\widehat{u}\left(k+i-1|k\right)}}\\ +\nu {\displaystyle \sum _{i=1}^{H_p}{\widehat{\delta }}^T\left(k+i|k\right)\widehat{\Delta }\left(k+i|k\right)}\end{array}$$E10

where $$\widehat{Z}\left(k+i|k\right)$$ is the future output vector, $$\widehat{u}\left(k+i-1|k\right)$$is the future control input vector, $$\widehat{\delta }\left(k+i|k\right)$$is the future synchronization errors, $$H_p$$is the prediction horizon, $$H_c$$is the control horizon, $$\lambda$$is the positive weighting factor used to adjust the control action, $$v$$is the non-negative weighting factor for the synchronization error.

Rewrite (6) as,

$$J\left(k\right)={\displaystyle \sum _{j=0}^{H_p-1}\Theta (\widehat{X},\widehat{u})+{\widehat{X}}^T\left(k+H_p|k\right)P_0{\widehat{X}}^T\left(k+H_p|k\right)}$$E11

where

$$\begin{array}{l}\Theta \left(\widehat{X},u\right)={\widehat{X}}^T\left(k+H_p-1+j|k\right)Q_j\widehat{X}(k+H_p-1+j|k)\\ +{\widehat{u}}^T\left(k+H_p-1+j|k\right)R_j\widehat{u}\left(k+H_p-1+j|k\right)\end{array}$$E12

$$P_0=C^{T}C+\nu C^T{L}^{T}LC$$E13

$$Q_j=\{\begin{array}{l}{C}^{T}C+\nu C^T{L}^{T}LCifj=0,1,\dots ,H_p-2\\ 0ifj=H_p-1\end{array}$$E14

$$R_j=\{\begin{array}{l}\infty I,ifj=0,1\dots ,H_p-H_c-1\\ \lambda I,ifj=H_p-H_c,\dots ,H_p-1\end{array},$$E15

where $$I$$ is the unit matrix with appropriate dimension.

Hence, minimization of the cost function (7) results in the synchronization control law,

$$u\left(k\right)=-K_{H_p-1}X\left(k\right)$$E16

where

$$K_{H_p-1}={(B^T{P}_{H_p-1}B+\lambda I)}^{-1}{B}^T{P}_{H_p-1}A$$E17

$$P_{j+1}=A^T{P}_{j}A-A^T{P}_{j}B{(B^T{P}_{j}B+R_j)}^{-1}{B}^T{P}_{j}A+Q_j$$E18

from the initial condition$$P_0$$.

Indeed, as receding horizon LQG control is a stationery feedback strategy, over an infinite interval, questions of stability naturally arise while solutions are slow to emerge. On the other hand, the stability of the proposed controller (12) can sometimes be guaranteed with finite horizons, even if there is no explicit terminal constraint. The finite horizon predictive control problem is normally associated with a time-varying RDE, which is related to the optimal value of the cost function. Attempts at producing stability result for MPC on the basics of its explicit input–output description have been remarkably unsuccessful, usually necessitating the abandonment of a specific control performance.

5. Stability analysis

Lemma 1. Consider the following ARE with an infinite-horizon linear quadratic control (Souza et al., 1996),

$$P=A^{T}PA-A^{T}PB{(B^{T}PB+R)}^{-1}{B}^{T}PA+Q$$E19

where

• $$\left[A,B\right]$$
• $$\left[A,Q^{1/2}\right]$$
• $$Q\ge 0\text{and}R>0.$$E20

Then

• there exists a unique, maximal, non-negative definite symmetric solution $$\overline{P}$$.

• $$\overline{P}$$$$i.e$$$$A-B{(B^T\overline{P}B+R)}^{-1}{B}^T\overline{P}A$$

Rewrite (15) as

$$P={(A-BK)}^{T}P(A-BK)+K^{T}RK+Q$$E21

In order to connect the RDE (14) to the ARE (15), the Fake Algebraic Riccati Technique (FART) is used as follows:

$$P_j={(A-B{K}_j)}^T{P}_j(A-B{K}_j)+K_j{}^T{R}_j{K}_j+{\overline{Q}}_j$$E22

where$${\overline{Q}}_j=Q_j-\left(P_{j+1}-P_j\right)$$. Clearly, while one has not altered the RDE in viewing it as a masquerading ARE, the immediate result from Lemma 1 and (17) can be obtained.

Theorem 1. Consider (17) with $${\overline{Q}}_j$$. If

• $$\left[A,B\right]$$
• $$\left[A,Q_j^{1/2}\right]$$
• $${\overline{Q}}_j\ge 0\text{and}{R}_j>0.$$E23

then $$P_j$$ is stabilizing,$$i.e$$. the closed-loop transition matrix

$${\overline{A}}_j=A-B{\left(B^T{P}_{j}B+R_j\right)}^{-1}{B}^T{P}_{j}A$$E24

has all its eigenvalues strictly within the unit circle.

Regarding the receding horizon strategy, only $$P_j$$ with $$j=H_p-1$$ will be applied. This leads to

$$P_{H_p-1}=A^T{P}_{H_p-1}A-A^T{P}_{H_p-1}B{\left(B^T{P}_{H_P-1}B+\lambda I\right)}^{-1}{B}^T{P}_{H_p-1}A+{\overline{Q}}_{H_P-1}$$E25

where$${\overline{Q}}_{H_p-1}=P_{H_p-1}-P_{H_p}$$. Then, the stability result of the control system can be given by the following theorem.

Theorem 2. Consider (19) with the weighting matrix$${\overline{Q}}_{H_p-1}$$. If

• $$\left[A,B\right]$$
• $$\left[A,{\overline{Q}}_{H_p-1}^{1/2}\right]$$
• $$P_{H_p-1}$$$$\lambda >0$$

then the controller (12) is stabilizing, i.e., the closed-loop transition matrix $${\overline{A}}_{H_p-1}=A-B{\left(B^T{P}_{H_p-1}B+\lambda I\right)}^{-1}{B}^T{P}_{H_p-1}A$$ has all its eigenvalues strictly within the unit circle.

Proof. The proof is completed by setting $$j=H_p-1$$ in Theorem 1.

It can be seen from the above theorem that the prediction horizon $$H_p$$ is a key parameter for stability, and an increasing $$H_p$$ is always favorable. This was the main motivation to extend the one-step-ahead control to long range predictive control. However, a stable linear feedback controller may not remain stable for a real system $$P\left(z\right)$$ with model uncertainty, which is normally related to stability robustness of the system. The most common specification of model uncertainty is norm-bounded, and the frequency response of a nominal model (3) can be obtained by evaluating:

$$\widehat{P}\left(z\right)=C{\left(zI-A\right)}^{-1}B$$E26

Then, the real system $$P\left(z\right)$$ is given by a ‘norm-bounded’ description:

$$P\left(z\right)=\widehat{P}\left(z\right)+{\Delta }_A,\text{for additive model uncertainties}$$E27

where $${\Delta }_A$$ is stable bounded operator, and $$P\left(z\right)$$ is often normalized in such a way that$$\Vert \Delta \Vert \le 1$$.

Because one does not know exactly what $$\Delta$$ is, various assumptions can be made about the nature of$$\Delta$$: nonlinear, linear time-varying, linear parameter-varying and linear time-invariant being the most common ones. Also, various norms can be used, and the most commonly used one is the ‘H-infinity’ norm $${\Vert \Delta \Vert }_{\infty }$$, which is defined as the worst-case ‘energy gain’ of an operator even for nonlinear systems. It then follows from the small-gain theorem that the feedback combination of this system with the uncertainty block $${\Delta }_A$$ will remain stable if

$$\overline{\sigma }\left[K_{H_p-1}\left(e^{j\omega T_s}\right)S\left(e^{j\omega T_s}\right)\right]{\Vert {\Delta }_A\Vert }_{\infty }<1$$E28

where $$\overline{\sigma }[\cdot ]$$ denotes the largest singular value, $$S\left(z\right)={\left[I+\widehat{P}\left(z\right)K_{H_p-1}\left(z\right)\right]}^{-1}$$is the sensitivity function. Note that (22) is only a sufficient condition for robust stability; if it is not satisfied, robust stability may nevertheless have been obtained. In practice, when tuning a controller, one can try to influence the frequency response properties in such a way as to make (22) hold.

6. Simulation study

Stability robustness with respect to variable control parameters will first be carried out. The $$y$$-axis of each graph indicates the maximum singular value of $$\overline{\sigma }\left[K_{H_p-1}\left(e^{j\omega T_s}\right)S\left(e^{j\omega T_s}\right)\right]$$, and the $$x$$-axis is the frequency range, $$\omega ={10}^{-2}~1$$Hz. Then, the performance of the proposed controller will be demonstrated in the presence of external disturbances and model uncertainties.

Consider the flywheel system with parameters given in (Zhu & Xiao, 2009), and assume that the rotor is spinning at a constant speed. As the eigenvalues of $$A_c$$ are:$$\pm 2.0353i$$, $$\pm 10.4i$$, $$\pm 149.3$$,$$\pm 149.3$$ , the open-loop continuous system is obviously unstable. With appropriate control parameters for the discrete-time model (sampling period $$T_s=0.008$$ s), such as$$H_p=6$$, $$H_c=1$$, $$\lambda =0.01$$, $$v=10$$, all of the eigenvalues of the closed-loop transition matrix $${\overline{A}}_{Hp-1}$$ are within the unit circle, which are: $$-0.782\pm 0.555i$$, $$-0.379$$, $$0.481$$, $$0.378$$, $$0.127\pm 0.195i$$and $$-0.027$$ respectively. In another word, the system can be stabilized with this feedback controller.

6.1. Stability robustness against control parameters

The prediction and control horizons are closely related to the stability of the closed-loop system. In the case of additive uncertainties, the maximum singular value $$\overline{\sigma }\left[K_{H_p-1}\left(e^{j\omega T_s}\right)S\left(e^{j\omega T_s}\right)\right]$$ against variation of prediction horizon is illustrated in Fig. 2, while$$H_c=1$$, $$\lambda =0.01$$and $$v=0$$ are set. It can be seen that a larger prediction horizon results in a smaller singular value, which means that the stability robustness of the control system can be improved. As a rule of thumb, $$H_p$$can be chosen according to$$H_p=\mathrm{int}\left(2{\omega }_s/{\omega }_b\right)$$, where $${\omega }_s$$ is the sampling frequency and $${\omega }_b$$ is the bandwidth of the process. Fig. 3 shows the singular value when the control horizon is varying. Clearly, a smaller control horizon $$H_c$$ may enhance the stability robustness of the control system. However, if the nominal model of the process is accurate enough, and the influence of model uncertainties is negligible, then $$H_c>1$$ is preferred for faster system responses.

Figure 2.

Maximum singular value $$\overline{\sigma }\left[K_{H_p-1}\left(e^{j\omega T_s}\right)S\left(e^{j\omega T_s}\right)\right]$$ against prediction horizon $$H_p$$

Figure 3.

Maximum singular value $$\overline{\sigma }\left[K_{H_p-1}\left(e^{j\omega T_s}\right)S\left(e^{j\omega T_s}\right)\right]$$ against control horizon $$H_c$$

The stability robustness bounds shown in Fig. 4 is obtained by varying$$\lambda$$, while$$H_p=6$$, $$H_c=1$$and $$v=0$$ are set. Clearly, a larger value of $$\lambda$$ can improve the stability robustness of the control system. This is because that the increasing $$\lambda$$ will reduce the control action and the influence of the model uncertainties on system stability will become less important. Consequently, the stability robustness can be enhanced. If$$\lambda \to \infty$$, the feedback action disappears and the closed loop is broken. In general, a larger $$\lambda$$ should be chosen when the system stability might be degraded due to significant model uncertainty. However, if the model uncertainty is insignificant, a smaller $$\lambda$$ would then be expected as the system response can be improved in this case, i.e., a decrease in the response time. In practice, a careful choice of $$\lambda$$ is necessary as it may have a large range of the values and is difficult to predetermine it.

The synchronization factor $$v$$is introduced to compensate the synchronization error of the rotor in radial direction. Fig. 5 shows that the influence of $$v$$ on stability robustness is not consistent over frequency. In particular, a lower value of $$v$$ can enhance the stability robustness at certain frequencies, but the performance will be degraded at higher frequencies. Another interesting observation is that the two boundaries for $$v=5$$ and $$v=10$$ are almost overlapping. It means that the stability robustness of the control system will not be affected if a further increase of $$v$$ is applied. In general, one can increase the prediction horizon and the synchronization control weighting factor so that the stability of the control system is maintained while the synchronization performance can be improved.

Figure 4.

Maximum singular value $$\overline{\sigma }\left[K_{H_p-1}\left(e^{j\omega T_s}\right)S\left(e^{j\omega T_s}\right)\right]$$ against weighting factor $$\lambda$$

Figure 5.

Maximum singular value $$\overline{\sigma }\left[K_{H_p-1}\left(e^{j\omega T_s}\right)S\left(e^{j\omega T_s}\right)\right]$$ against synchronization factor $$v$$

6.2. Disturbances on magnetic forces

In this simulation, force disturbances are introduced to the bearings of the rotor at different time instants, and amplitudes are 0.5N, -0.5N, 0.5N and -0.5N on xa-axis, xb-axis, ya-axis and yb-axis respectively. The duration of 0.2 seconds for each disturbance is assumed. Figs. 6-11 show the numerical results of the control algorithm when$$H_p=10$$, $$H_c=1$$, $$\lambda =0.01$$are set for the two cases: with$$v=0$$, and$$v=10$$. Clearly, without cross-coupling control action due to$$v=0$$, evident synchronization errors and a conical whirl mode during the transient responses are resulted. However, when $$v=10$$ is introduced, the synchronization performance can be improved significantly, especially in terms of the rolling angles, as shown in Fig. 11. Therefore, with adequately selected control parameters the improved synchronization performance as well as guaranteed stability of the FESS can be obtained, and in consequence, the whirling rotor in the presence of disturbances would be suppressed near the nominal position.

Figure 6.

Radial displacements of the rotor along x-axis

Figure 7.

Radial displacements of the rotor along y-axis

Figure 8.

Control currents to bearings along x-axis

Figure 9.

Control currents to bearings along y-axis

Figure 10.

Displacements of rotor mass center

Figure 11.

Rolling angles of rotor mass center

7. Conclusion

In this chapter, stability problem of magnetic bearings for a flywheel energy storage system has been formulated, and a synchronization design has been presented by incorporating cross-coupling technology into the optimal control architecture. The basic idea of the control strategy is to minimize a new cost function in which the synchronization errors are embedded, so that the gyro-dynamic rotation of the rotor can be effectively suppressed.

However, as optimal control, using receding horizon idea, is a feedback control, there is a risk that the resulting closed-loop system might be unstable. Then, stability of the control system based on the solution of the Riccati Difference Equation has also been analyzed, and some results are summarized. The illustrative example reveals that with adequately adjusted control parameters the resulting control system is very effective in recovering the unstable rotor and suppressing the coupling effects of the gyroscopic rotation at high spinning speeds as well as under external disturbances and model uncertainties.