Accurate Suspension Force Modeling and Its Control System Design Based on the Consideration of Degree-of-Freedom Interaction

Abstract

In this study, an accurate suspension force modeling method for the magnetic bearings of flywheel batteries considering degree-of-freedom (DOF) interactions and their control system is proposed to solve the problem that the traditional flywheel battery suspension force model does not consider DOF interactions, which makes the control system control effect poor. Firstly, according to the structural characteristics of the flywheel battery used, a suspension force model is established for the radial and axial magnetic bearings, which are most seriously interfered with by the torsional degrees of freedom of the flywheel battery. Next, by proposing DOF interaction factors, the complex changes due to DOF interactions are cleverly summarized into several interaction factors applied to the fundamental model to achieve accurate suspension force modeling considering DOF interactions. To better adapt the established accurate model and ensure precise control of the flywheel battery system under various working conditions, the firefly algorithm is employed to optimize the BP neural network (FA-BPNN). This optimization regulates the control system’s parameters, enabling the achievement of optimal control parameters in different scenarios and enhancing control efficiency. Compared to the flywheel battery controlled using the fundamental model, the radial and axial displacements are reduced by more than 30 percent and 20 percent, respectively, in the uphill condition using the accurate model with FA-BPNN.

1. Introduction

A flywheel battery (flywheel energy storage system) can break through the limitations of traditional chemical batteries with the advantages of no pollution, high energy conversion efficiency, and a long service life. When applied to electric vehicles, the energy efficiency of electric vehicles can be significantly improved [1,2,3,4,5,6]. Vehicle-mounted flywheel batteries, due to the application background of the vehicle environment, will inevitably be subjected to vehicle driving conditions and road conditions; due to the complexity of the disturbance, flywheel batteries are very easily destabilized [7]. In order to improve the stability of flywheel batteries, a novel shaftless flywheel structure is proposed in [8] to increase the energy density of the flywheel while ensuring good control linearity and load carrying capacity through the adjustable stiffness of each magnetic bearing; Ref. [9] proposes a virtual inertia spindle vehicle-mounted flywheel battery structure. It is based on the shaftless structure. The integration is further improved, and the spatial arrangement of magnetic bearings is optimized to enhance overall stability. Additionally, magnetic circuit decoupling is carried out, and a new disk-shaped flywheel battery with high perturbation resistance is proposed [10]. While improving the stability through flywheel structure optimization, a rotor disturbance suppression method based on an improved linear extended state observer is proposed in [11]; in [12], a novel stabilization control method is proposed to reduce the synchronous vibration and gyroscopic effect of the flywheel; and in [13], a new robust control method is proposed to enhance the rotational mode suppression of the flywheel rotor in order to improve the stability of the flywheel operation through the control algorithm. However, at present, neither the flywheel battery structure topology design nor the control strategy considers the interaction between each DOF. The interaction between the DOF will directly affect the control efficiency of each DOF.

To resolve the interaction of DOF, the five new DOF magnetic bearing designed by [14] separate the radial and axial magnetic circuits, so that the radial and axial DOF cannot interfere with each other. Further, in [15] the interaction of the radial two DOF is reduced from the structural point of view by the design of a six-pole radial magnetic bearing with a combination of composite magnetic circuits. Meanwhile, in [16], a hybrid magnetic bearing for flywheel batteries is decoupled between axial and radial DOF using an improved neural network based on a mathematical model. It can be concluded from the above that existing measures to reduce DOF interference mainly focus on optimizing the magnetic circuit to address DOF interference during structural design and implementing decoupling control for DOF interference during the control process. However, for an existing structure, it is no longer possible to reduce the DOF interaction through structural design. Moreover, decoupling control also needs to be based on the suspension force model. Therefore, compared to directly modifying the suspension force model to obtain an accurate suspension force model that accounts for a DOF interaction, these two methods lack universality and simplicity. Meanwhile, the interaction of torsional freedom on axial and radial DOF has rarely been investigated. For vehicle-mounted magnetic suspension flywheel batteries, under certain driving conditions, such as uphill and downhill, the torsional DOF offset and radial DOF offset as well as axial offset will occur at the same time, and the influence of the torsional DOF will simultaneously affect the effect of the suspension force of the radial and axial magnetic bearings. Therefore, it is important to consider the interaction of DOF, especially the torsional DOF for the radial and axial magnetic bearings, when modeling the suspension force of flywheel batteries, and apply them to the control to improve stability.

The main uses of wheel battery suspension force modeling for flywheel batteries are currently the equivalent magnetic circuit method, the Maxwell tensor method, and the mathematical model established by considering eddy current effects. The equivalent magnetic circuit method [17,18,19,20] is the most classical way to model the suspension force of magnetic suspension bearings. It is widely applied in the modeling of active, hybrid, and constant-current-source-biased magnetic suspension bearings. It is also used in the design of corresponding control systems. In [21], a novel six-pole, five-degree-of-freedom AC hybrid magnetic suspension bearing equivalent magnetic circuit model is proposed for the parameter design optimization of magnetic bearings. In [22], a suspension force model based on the equivalent magnetic circuit method is used to achieve control for a novel hybrid magnetic suspension bearing. However, conditions such as magnetic leakage, edge effects, hysteresis losses, and magnetic saturation need to be neglected when using the equivalent magnetic circuit method, which leads to a reduction in the accuracy of the suspension model. Therefore, in [23] leakage flux is included in the suspension modeling process, but it does not take into account edge effects. In [24], by splitting the magnetic field, edge effects as well as leakage flux are added to the suspension force model, which greatly improves the accuracy of the model. Existing suspension force modeling focuses on improving the accuracy of flywheel static-state modeling. However, none of the existing modeling methods consider the offsets of various degrees of freedom during the actual operation of the flywheel and incorporate the interference of DOF into the suspension force model. Moreover, the interference of DOF has a great influence on the suspension force model, meaning that the existing modeling methods cannot be well applied to the modeling of the suspension force of vehicle-mounted flywheel batteries under complex working conditions. At the same time, the flywheel conditions were further categorized into two categories with high and low DOF interaction, considering the accurate modeling of DOF. For different working conditions, the optimal values of the control parameters of the adopted control system are not the same. If the control parameters can be adjusted in real time based on the established accurate model to realize the control of the optimal control parameters under working conditions, the control efficiency of the control system can be greatly increased. The use of neural networks can be realized by the vehicle’s sensor data to determine the vehicle’s working conditions and facilitate real-time adjustment of control parameters. But the traditional BP neural network easily falls into local optimization and the initial parameters have a great influence on the results [25,26,27]. Ref. [28] proposes a control strategy based on CNN + LSTM + ATTENTION, so that the control parameters of the flywheel battery can change with the changing working conditions, but it is not analyzed based on an accurate model.

To solve these problems, this paper proposes an accurate suspension force modeling and control system design based on the consideration of DOF interaction. A magnetic suspension flywheel battery with a multi-function air gap is taken as an example. Firstly, the variation in suspension forces under DOF is analyzed, and it is found that the interaction of the DOF in the flywheel structure used is mainly manifested in the influence of the torsional DOF on the radial and axial suspension forces. Therefore, taking the law of change in the air gap magnetic flux density of the radial and axial magnetic bearings during torsion as a reference, assumptions will be made regarding the laws of change in the statically biased magnetic flux density and the resultant suspension force caused by torsion. Additionally, assumptions will be made for the force–current stiffness during torsion based on the change in the air gap magnetic flux density after applying the control current during torsion. Based on this, the interaction factors will be summarized, and accurate modeling will be achieved on the basis of the fundamental model. Finally, taking advantage of the firefly algorithm’s strength in optimality searching, the firefly algorithm is combined with the BP neural network. The firefly algorithm is used to adjust the weights and thresholds of the BP neural network, making the neural network parameters optimal. Then, this is applied to judge the vehicle’s working condition and adjust the control parameters. The experimental results show that the established accurate model is highly superior to the traditional suspension force model when flywheel torsion occurs, and in the control system using the accurate model and the control parameter adjustment based on FA-BPNN, after comparing with the control system based on the traditional model and the control system without parameter adjustment, the offset of the flywheel in the case where the DOF interaction case with each other is greatly reduced.

2. The Effect of DOF Interaction and Fundamental Model of Magnetic Bearings

2.1. Variation in Suspension Forces Due to DOF Interaction

As shown in Figure 1, the flywheel battery consists of a 5-DOF magnetic bearing, an outer rotor brushless DC motor, a bowl-shaped flywheel, and a protective housing with sensors. The multi-function air gap for the overall prototype has three main functions: guiding the magnetic circuit to the axial air gap and minimizing the control air gap using auxiliary permanent magnets, providing a complete and isolated control magnetic circuit, and enabling axial and torsional stators to share permanent magnets for integration through its addition. The finite element model used is based on an actual prototype, and the final energy error of the finite element calculation is kept within 1% to ensure the accuracy of the simulation.

The working environment of a vehicle-mounted flywheel battery makes its offset multiple degrees of freedom; it is necessary to additionally analyze the influence of magnetic bearings of a certain degree of freedom on magnetic bearings that are not of their corresponding degree of freedom. Due to the magnetic circuit decoupling design of the vehicle-mounted flywheel battery, the influence of axial offset on the radial magnetic bearing and that of radial offset on the axial magnetic bearing are very small. In the case of torsion, the static suspension force of the flywheel changes because of the variations in the air gap and magnetic flux density. Due to the limiting effect of the auxiliary bearing, the maximum torsion angle is 0.5 deg. As can be seen from the finite element simulation in Figure 2, as the torsion angle increases, the maximum suspension force generated by the radially biased magnetic field is about −25 N. Similarly, the axial magnetic suspension force increases from the suspension force of the equilibrium gravity to approximately 650 N. This greatly affects the stability and safety of the flywheel. Therefore, the following analyses of DOF interaction will focus on the radial and axial suspension force changes due to the torsional DOF.

2.2. Fundamental Modeling of Radial and Axial Suspension Forces

The equivalent magnetic circuit method is used to model the fundamental radial suspension force first. In Figure 1, F_r represents the radial permanent magnet potential; R_A1 is the upper air gap reluctance of phase A in the three-phase radial magnetic bearing, whereas R_A2 is the radial lower air gap reluctance of phase A. The remaining two phases have the same structure. Moreover, R_m is the radial permanent magnet reluctance. F_ca, F_cb, and F_cc are the magnetic potentials generated by the three-phase magnetic bearing control coil. Further, $\varphi$_Jpr (J = A, B, C) is the A-\B-\C-phase biased flux, and $\varphi$_Jcr is the A-\B-\C-phase control flux. The magnetic flux is as follows:

$${\varphi }_{\mathrm{J}pr}={\displaystyle \frac{F_{pr}}{R_m+R_{\mathrm{J}1}+R_{\mathrm{J}2}}}; {\varphi }_{\mathrm{J}cr}={\displaystyle \frac{F_{c\mathrm{J}}}{\frac{3}{2}{R}_{\mathrm{J}2}}}.$$

(1)

The radial synthetic flux $\mathsf{\varphi }$_Jr and the radial suspension force F_J under a single pole are as follows:

$$\left\{\begin{array}{l}{\varphi }_{\mathrm{J}r}={\varphi }_{\mathrm{J}cr}+{\varphi }_{\mathrm{J}pr}, F_{\mathrm{J}}={\displaystyle \frac{{({\varphi }_{\mathrm{J}pr}+{\varphi }_{\mathrm{J}cr})}^2}{2{\mu }_0{S}_r}}\hfill \\ F_x=F_{\mathrm{A}}-{\displaystyle \frac{1}{2}}{F}_{\mathrm{B}}-{\displaystyle \frac{1}{2}}{F}_{\mathrm{C}}, F_y={\displaystyle \frac{\sqrt{3}}{2}}{F}_{\mathrm{B}}-{\displaystyle \frac{\sqrt{3}}{2}}{F}_{\mathrm{C}}.\hfill \end{array}\right.$$

(2)

The forces generated by the ABC poles are projected onto the x\y axis and divided into x-direction and y-direction forces. The x-direction is parallel to the A-phase axis and the y-direction is orthogonal to the A-phase axis. Figure 3 shows the results of Equation (2) compared with the simulation: Figure 3a shows the x-direction and Figure 3b shows the y-direction. When the current exceeds a certain value, the calculated suspension force will be larger than the finite element simulation value because magnetic leakage and magnetic saturation are neglected in the magnetic circuit calculation. Therefore, the modification of the model and the corresponding control system involves the modification and construction of a nonlinear model.

The axial magnetic circuit is also shown in Figure 1. F_mm and R_mm, respectively, represent the magnetomotive force and permeability generated by the main permanent magnet. F_ma and R_ma, respectively, denote the magnetomotive force and reluctance of the auxiliary permanent magnet. The reluctances R_n, R_z, and R_h are associated with the air gaps through which the axial and twisting biased flux traverse; specifically, these air gaps are the axial torsion shared air gap, the axial air gap, and the multi-function air gap in sequence. R_cz is the air gap reluctance through which the axial control flux travels. Considering that the axial biased suspension force is generated by the flux between the two air gaps and the axial control suspension force is generated by the flux in only one air gap, the equations for the axial flux as well as the suspension force are expressed as follows, where $\varphi$_cz is the axial control flux, $\varphi$_pz is the axial biased flux, $\varphi$_z is the axial flux, and F_z is the axial suspension force. $\varphi$_pz in turn consists of two parts, one is the flux $\varphi$_pzm and the other is the flux $\varphi$_pzs produced by the auxiliary permanent magnet.

$$\left\{\begin{array}{l}\begin{array}{l}{\varphi }_{pzm}={\displaystyle \frac{F_{mm}{R}_{mm}}{(R_z+R_n)R_{ma}+(R_{mm}+R_h)(R_z+R_n+R_{ma})}}\\ {\varphi }_{pza}={\displaystyle \frac{F_{ma}(R_{ma}+R_h)}{(R_z+R_n)(R_{mm}+R_h)+R_{ma}(R_z+R_n+R_{mm}+R_h)}}\end{array}\hfill \\ \begin{array}{l}{\varphi }_{pz}={\varphi }_{pzm}+{\varphi }_{pza},{\varphi }_{cz}={\displaystyle \frac{N_z{i}_z}{R_z+R_{cz}+R_h}},{\varphi }_z={\varphi }_{cz}+{\varphi }_{pz}\\ F_z={\displaystyle \frac{{({\varphi }_{cz}+{\varphi }_{pz})}^2}{2{\mu }_0{\mathrm{S}}_{{\mathrm{z}}_1}}}+{\displaystyle \frac{{{\varphi }_{pz}}^2}{2{\mu }_0{\mathrm{S}}_{\mathrm{n}}}}={\displaystyle \frac{{{\varphi }_z}^2}{2{\mu }_0{\mathrm{S}}_{{\mathrm{z}}_1}}}+{\displaystyle \frac{{{\varphi }_{pz}}^2}{2{\mu }_0{\mathrm{S}}_{\mathrm{n}}}}.\end{array}\hfill \end{array}\right.$$

(3)

3. Modeling of Radial and Axial Suspension Forces Considering the Torsion DOF Interaction

3.1. Characteristics of the Variation in the Magnetic Flux Density in the Torsion DOF Interaction

When downhill, the flywheel is twisted around the y-axis in the direction of the A-phase, where the A-phase is defined as the main pole of the torsion, and the two symmetrical BC-phases are defined as the secondary poles of the torsion. The top-view of the three-phase air gap changes obtained from the finite element simulation is shown in Figure 4a, and the B-phase magnetic flux density distribution in the torsion, radial offset, and non-torsion cases is given in Figure 4b. It can be seen from Figure 4a,b that in the case of torsion, the air gap between the flywheel and the magnetic bearing has undergone significant changes. Moreover, the change in magnetic flux density caused by torsion is quite different from that caused by radial displacement. The changes in the air gap and magnetic flux density of the axial magnetic bearing are similar to those of the radial magnetic bearing.

3.2. Radial Biased Flux Correction Under the DOF Interaction

The torsional DOF interaction factors for radial magnetic bearings are first extracted starting from the distortion of the radially biased magnetic flux density. As shown in Figure 5, as can be seen from the finite element simulation, as the torsion angle increases, the static biased magnetic flux density of the main poles of the torsion decreases continuously and the static biased magnetic flux density of the secondary poles of the torsion is constantly distorted, meaning that the farther away they are from the torsion main pole, the larger the magnetic flux density, which makes the overall magnetic flux density increase.

Since the cross-sectional area of the three-phase air gap is constant, the change in the mean value of the magnetic flux density can be further expressed as a change in the magnetic flux. The fluxes under torsion are as follows:

$$\left\{\begin{array}{l}{\varphi }_{np\mathrm{M}}=K_{n\mathrm{M}}{B}_{p\mathrm{M}}{S}_r=K_{n\mathrm{M}}{\varphi }_{p\mathrm{M}}=(0.997-0.176{\theta }_N){\varphi }_{p\mathrm{M}}\hfill \\ {\varphi }_{np\mathrm{S}}=K_{n\mathrm{S}}{B}_{p\mathrm{S}}{S}_r=K_{n\mathrm{S}}{\varphi }_{p\mathrm{S}}=(0.998-0.108{\theta }_N){\varphi }_{p\mathrm{S}},\hfill \end{array}\right.$$

(4)

where K_nM is the torsional DOF interaction factor for the static biased flux of the main pole of torsion, and K_nS is the torsional DOF interaction factor for the static biased flux of the secondary poles of torsion. B_pM is the average magnetic flux density of the torsion main poles, and B_pS is the average magnetic flux density of the torsion secondary poles. $\varphi$_pM is the static biased magnetic flux density of the torsion main poles after the correction, and $\varphi$_pS is the static biased magnetic flux density of the torsion secondary poles after correction.

For the torsion main pole, the average value of the pole can be summarized with the torsion of the law of change. However, for the secondary poles of the torsion, the distortion effect needs further consideration. Although the geometry of the two phases B and C is still, with the A-phase presenting a 120-degree arrangement, the actual effect of the distortion of the magnetic pulling force is not the same as the 120-degree arrangement. Therefore, for the torsion case of the torsion secondary poles, the force action change is further corrected by adjusting the magnetic suspension force’s angle of action using the torsion correction angle based on the 120-degree distribution.

Considering that the magnetic flux density of the secondary poles of the torsion at this point increases the further away they are from the main poles of the torsion, the torsion correction angle is defined as the angle between the two phases of BC and phase A increases as the torsion increases. As shown in Figure 6, such a change causes the synthetic suspension force F_n to increase. The radial synthetic suspension force F_nx, F_ny in this case, is as follows (F_n is in the reversed direction of F_x when A is the main pole of the torsion):

$$\left\{\begin{array}{l}\begin{array}{l}{F}_{nx}=F_{n\mathrm{A}}-\cos (120+{\theta }_n)F_{n\mathrm{B}}-\cos (120+{\theta }_n)F_{n\mathrm{C}}\\ =F_{n\mathrm{M}}-2\cos (120+{\theta }_n)F_{n\mathrm{D}}\end{array}\hfill \\ \begin{array}{l}{F}_{ny}=\sin (120+{\theta }_n)F_{n\mathrm{B}}-\sin (120+{\theta }_n)\\ =\sin (120+{\theta }_n)F_{n\mathrm{D}}-\sin (120+{\theta }_n)F_{n\mathrm{D}},\end{array}\hfill \end{array}\right.$$

(5)

where θ_n is the torsional correction angle and F_nA, F_nB, and F_nC are the magnetic suspension forces generated by each phase of the radial magnetic bearing during torsion. Combining K_nM as well as K_nS, F_nA, F_nB, and F_nC in this case yields the following:

$$\left\{\begin{array}{l}{F}_{n\mathrm{A}}=F_{pn\mathrm{A}}=F_{pn\mathrm{M}}={\displaystyle \frac{{(K_{n\mathrm{M}}{\varphi }_{\mathrm{Apr}})}^2}{2{\mu }_0{\mathrm{S}}_r}}\hfill \\ F_{n\mathrm{B}}=F_{n\mathrm{C}}=F_{pn\mathrm{B}}=F_{pn\mathrm{C}}=F_{pn\mathrm{S}}={\displaystyle \frac{{(K_{n\mathrm{S}}{\varphi }_{\mathrm{Bpr}})}^2}{2{\mu }_0{\mathrm{S}}_r}}={\displaystyle \frac{{(K_{n\mathrm{S}}{\varphi }_{\mathrm{Cpr}})}^2}{2{\mu }_0{\mathrm{S}}_r}},\hfill \end{array}\right.$$

(6)

where F_pna, F_pnb, and F_pnc are the per-phase magnetic suspension forces generated by the static biased flux in the torsional case. In order to obtain the θ_n, it is assumed that θ_n is accompanied by a change in the torsion angle θ_N that exhibits a linear change in the angle from 0, denoted as

$${\theta }_n=K_n\cdot {\theta }_N,$$

(7)

where K_n is the torsional DOF interaction factor for the radial synthetic force. Substituting a torsion of 0.1 degrees into the Equation (5) inverse calculation allows a K_n of 32.42 to be obtained.

3.3. Radial Control of Flux Variation and Radial Accurate Suspension Forces Considering Torsional DOF Interaction

Figure 7 shows the comparison from the finite element simulation of the interpole magnetic flux densities when the control current is increased from 0 A to 1 A for a torsion angle of 0 deg as well as for a torsion angle of 0.5 deg, and the current–magnetic densities have the same stiffness, although the shapes of the interpole densities are different. Therefore, the control magnetic flux of each phase need not be corrected, since all the magnetic flux density distortion in torsion is normalized to the torsion correction angle.

Adding the control flux to Equation (5) yields the following formula for the suspension force in torsion.

$$\left\{\begin{array}{l}{F}_{nx}=F_{n\mathrm{A}}-\cos (120+{\theta }_n)F_{n\mathrm{B}}-\cos (120+{\theta }_n)F_{n\mathrm{C}}\\ ={\displaystyle \frac{{(K_{n\mathrm{M}}{\varphi }_{\mathrm{M}pr}+{\varphi }_{\mathrm{A}cr})}^2}{2{\mu }_0{S}_r}}-\cos (120+{\theta }_n){\displaystyle \frac{{(K_{n\mathrm{D}}{\varphi }_{\mathrm{D}pr}+{\varphi }_{\mathrm{B}cr})}^2}{2{\mu }_0{S}_r}}-\cos (120+{\theta }_n){\displaystyle \frac{{(K_{n\mathrm{D}}{\varphi }_{\mathrm{D}pr}+{\varphi }_{\mathrm{C}cr})}^2}{2{\mu }_0{S}_r}}\\ F_{ny}=\sin (120+{\theta }_n)F_{n\mathrm{B}}-\sin (120+{\theta }_n)F_{n\mathrm{C}}\\ =\sin (120+{\theta }_n)({\displaystyle \frac{{(K_{n\mathrm{D}}{\varphi }_{\mathrm{D}pr}+{\varphi }_{\mathrm{B}cr})}^2}{2{\mu }_0{S}_r}}-{\displaystyle \frac{{(K_{n\mathrm{D}}{\varphi }_{\mathrm{D}pr}+{\varphi }_{\mathrm{C}cr})}^2}{2{\mu }_0{S}_r}}).\end{array}\right.$$

(8)

Figure 8a,b show the comparison between the exact suspension force model (solid colors) and the FEA model (rainbow colors) in the radial x-direction and y-direction, respectively.

When the torsion angle is small and the control current applied is low, the two are almost coincident, while when the torsion angle and the control current increase, the error between the two gradually increases. The maximum deviation in the x-direction, for example, is 7% for a maximum torsion angle of 0.5 degrees and a maximum control current of 1 A. The difference is 28% for the suspension model regardless of the interference between the DOF. This indicates that summarizing the interference between the DOFs in terms of interaction factors K_nm and K_ns and incorporating it into the model leads to a great improvement in the accuracy of the model. The rainbow-colored part is represented as a smooth surface, which also means that the force–current stiffness is unchanged.

3.4. Axial Accurate Suspension Forces Considering Torsional DOF Interaction

The torsional DOF interaction factors for the axial magnetic bearing are also analyzed starting from the biased flux. The interference of the torsional DOF for the axial magnetic bearing is also seen in Figure 9 to cause a large distortion of the air gap magnetic flux density. Since it is unnecessary to consider the synthesis of the three phases like in the radial direction, the torsional DOF interference enables the aberrations in the axial density to be represented by the change in the average value of the axial air gap magnetic flux density. The magnetic flux density between the two air gaps is fitted by quadratic term fitting and cubic term fitting. The axial static biased flux factors K_pz1 and K_pz2 are obtained by dividing them by the sum of the average values of the three-phase magnetic flux density of each air gap when there is no torsion. Then, the expression for the axial static biased magnetic flux density is obtained by applying it to Equation (9).

$$\left\{\begin{array}{l}{K}_{pz1}=0.882{\theta }_N{}^2-0.1618{\theta }_N+1\\ K_{pz2}=1.26{\theta }_N{}^3+(-0.315{\theta }_N{}^2)+0.102{\theta }_N+1\\ F_{pz}={\displaystyle \frac{{(K_{pz1}{\varphi }_{pz})}^2}{2{\mu }_0{\mathrm{S}}_{{\mathrm{z}}_1}}}+{\displaystyle \frac{{(K_{pz2}{\varphi }_{pz})}^2}{2{\mu }_0{\mathrm{S}}_{\mathrm{n}}}}.\end{array}\right.$$

(9)

With the axial control flux passing through the axial air gap only, just the change in magnetic flux density is analyzed after applying the control current in this air gap.

The current–magnetic flux density stiffness in the axial direction, shown in Figure 10b, is the same as that in the radial direction, acting on the distorted air gap but with the same change in the magnetic flux density at each point, so that the force–magnetic flux density stiffness curve can be obtained from the translation of the curve when no torsion occurs. The formula for the overall axial magnetic tension is as follows:

$$F_z={\displaystyle \frac{{(K_{pz1}({\varphi }_{pz}+{\varphi }_{cz}))}^2}{2{\mu }_0{\mathrm{S}}_{{\mathrm{z}}_1}}}+{\displaystyle \frac{{(K_{pz2}{\varphi }_{pz})}^2}{2{\mu }_0{\mathrm{S}}_{\mathrm{n}}}}.$$

(10)

The force–current stiffness curves at 0 deg, 0.3 deg, and 0.5 deg of torsion are selected to compare the results of finite element simulations as well as the accurate suspension force model calculations in Figure 10a. When the torsion angle reaches 0.5 deg, the effect of magnetic leakage due to distortion between the air gaps makes the calculated value of the suspension force model larger than the finite element simulation. However, the overall error between the two is small, indicating that it is feasible to summarize the effect of torsion on the axial direction as torsional factors (K_pz1 as well as K_pz2).

4. Establishment of FA-BPNN-Based Control Parameter Adjustment Model for Flywheel Battery System

4.1. Firefly Algorithm

The basic idea of the firefly algorithm (FA) is as follows: There is a certain range of randomly shiny firefly individuals, and such individuals will also be attracted by the light emitted by other individuals; the brighter the light, the greater the attraction. The brighter individuals randomly “move”, other individuals “gather” around the brighter individuals, and in the process of gathering and moving, any one of the firefly flight positions changes because of its real-time adjustment to the brightest light-emitting individual. Finally, each individual firefly will gather around the brightest individual.

Brightness I(r) is defined as follows:

$$I(r)=I_0{e}^{-\gamma r^2},$$

(11)

where I₀ is the original luminous brightness of the firefly’s location, which is related to the value of the objective function, and the higher its own luminous brightness, the higher the value of the objective function on behalf of; γ is the luminous brightness absorption coefficient; and r is the spacing between any two fireflies i and j at their locations X_i and X_j, also known as the Cartesian distance.

r is calculated by the following formula:

$$r=‖X_i-X_j‖=\sqrt{{\displaystyle \sum _{k=1}^d{(X_{i,k}-X_{j,k})}^2}},$$

(12)

where X_i,k, X_j,k are kth dimensional coordinate values of the spatial coordinates X_i, X_j of firefly i and firefly j. k = 1~d and d is the dimension of the problem being solved.

Define the attractiveness β(r), since the attractiveness of fireflies is proportional to the brightness, thus.

$$\beta (r)={\beta }_0{e}^{-\gamma r^2},$$

(13)

where β₀ is the attractiveness at r = 0. Firefly i will move in the direction of other fireflies j that are brighter than it:

$$X_i^{\prime }=X_i+{\beta }_0{e}^{-\gamma r^2}(X_j-X_i)+\alpha (rand-{\displaystyle \frac{1}{2}}),$$

(14)

where $X_i^{\prime }$ is the position of firefly i after moving; α is the step factor; rand is a random real number uniformly distributed in [0, 1]; the second term on the right side of the equation represents the attraction of fireflies to light intensity; and the third term is a random perturbation term.

The updated equation for the position of the brightest firefly is obtained through Equation (14) as follows.

$$X_{besti}^{\prime }=X_{besti}+\alpha (rand-{\displaystyle \frac{1}{2}}),$$

(15)

where X_besti is the position of the brightest firefly; $X_{besti}^{\prime }$ is the updated position of the brightest firefly.

4.2. Firefly-Optimized BP Neural Network for Control Parameter Adjustment

The displacement detected by the sensors in each direction of the vehicle flywheel battery varies as it undergoes different working conditions. The results detected by the five position sensors, x₁, x₂, y₁, y₂, and z₁ are used as input characteristic parameters in the road condition judgment process. Before establishing the firefly algorithm-optimized Bp neural network (FA-BP) for control parameter adjustment, the vehicle’s working conditions need to be judged first to determine the neural network’s output. There are five common vehicle working conditions: vehicle acceleration, vehicle deceleration, vehicle turning, vehicle uphill, and vehicle downhill. And each working condition has its corresponding PID control parameters (K_p, K_i, and K_d).

Based on the above, a control parameter adjustment model using the firefly algorithm to optimize the BP neural network was established, as shown in Figure 11.

The parameters of input and output are determined by the type of working conditions and the initialization of the population in the FA, and then the structure of the BP neural network is determined. Each individual firefly corresponds to the threshold and weight in the BP neural network. The firefly adaptation f_f is calculated as follows:

$$f_f=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{p=1}^N{(h_p-h_m)}^2}},$$

(16)

where h_p when the actual output values of all pth training vectors, p = 1 to N; h_m is the predicted value of all training vectors; and N is the number of trainings.

The BP neural network selects the optimal firefly individual by varying the firefly position, firefly brightness, and decision radius, and replaces them with initial weights and thresholds to obtain the optimal solution. Test set data from different working conditions are brought into the trained BP neural network and when the conditions are met, the output determines the type of working conditions and control parameters of the vehicle. The target control parameters for different working conditions are shown in the

Table 1. Control parameters.

Working Conditions	Acceleration			Deceleration			Turning
Characteristic parameters	x1 /x2	y1 /y2	z	x1 /x2	y1 /y2	z	x1 /x2	y1 /y2	z
Proportionality coefficient Kp	32.5 /5	32.5 /5	8	35.5 /2.3	35.5 /2.3	12	7.8 /1.45	29.2 /1.75	0
Integral coefficient Ki	0 /0	0 /0	0	0 /0	0 /0	0	0.074 /0.02	0.273 /0.03	0
Differentiation coefficient Kd	0.1 /0.02	0.1 /0.02	0.03	0.08 /0.01	0.08 /0.01	0.06	0.102 /0.03	0.122 /0.03	0
Working conditions	Uphill			Downhill			General situation
Characteristic parameters	x1 /x2	y1 /y2	z	x1 /x2	y1 /y2	z	x1 /x2	y1 /y2	z
Proportionality coefficient Kp	3.72 /1.12	1.37 /0.85	12.6	14.6 /5.72	14.6 /5.72	40	1.6 /0.85	1.6 /0.85	9.6
Integral coefficient Ki	0.024 /0.01	0.015 /0.01	0.02	0.009 /0.07	0.009 /0.07	0.28	0 /0	0 /0	0
Differentiation coefficient Kd	0.047 /0.02	0.03 /0.02	0.09	0.073 /0.02	0.073 /0.02	0.142	0.15 /0.1	0.15 /0.1	0.053

; the last one is the control parameter to be used without parameter adjustment.

Figure 12 shows the flow chart of FA-BPNN.

The population size of fireflies is set to 50. The maximum number of iterations is 500, and the error threshold is set at 0.01. The luminous brightness absorption coefficient γ is 1.0, the maximum attractiveness β is 1.0, and the step factor is 0.025. Thirty sets of feature parameters for each common working condition obtained from experimental measurements are selected, and a total of 150 sets of data for five common working conditions are used as training parameters to train the neural network. Based on the initial parameters to establish the FA-BP neural network control parameter adjustment model, the 150 groups of data are used to train the model to calculate the model error, which is the size of the individual fitness, and according to the size of the fitness of the firefly position replacement, in order to find the best firefly individual, the best firefly individual position is used as the optimal BP neural network parameter, based on which the relationship between the working condition characteristic parameters and the working conditions are determined, and the corresponding control parameters are obtained. When the error meets the conditions, the training is completed; otherwise, it is recalculated.

To verify the correctness of the parameter adjustment model, another fifty sets of working condition characteristic data were selected, and it can be seen from Figure 13 that the number of times that the working condition can be accurately judged is 47 times, where 1–5 represent five working conditions and their corresponding control parameters. Therefore, the accuracy of using FA-BPNN to adjust the control parameters reaches 94%. Meanwhile, the BP neural network optimized by the genetic algorithm (GA-BPNN) and the BP neural network (BPNN) was employed to conduct working condition judgments under the same set of working condition characteristic parameters. Here, Figure 13b presents the result of the judgment made by GA-BPNN. Out of 50 sets of data, 45 sets were judged correctly, with an accuracy rate of 90%. And Figure 13c also shows the judgment result of BPNN. Among the 50 sets of data, 42 sets were judged correctly, and its accuracy rate was 84%.

To verify the superiority of the firefly-optimized BP neural network algorithm (FA-BPNN), further comparative training was conducted with BPNN and GA-BPNN in the context of vehicle-mounted flywheel battery operating condition judgment. The training error curves are presented in Figure 14. By observing the error convergence of the three algorithms, it is found that the error decline rates of all three is very rapid at the beginning. Eventually, FA-BPNN reaches the target error at around 200 iterations, while GA-BPNN fails to reach the target error within the given 250 iterations. BPNN falls into a local optimum at around 250 iterations and remains unable to reach the target error even after continuous training up to 500 iterations.

5. Experiments

5.1. Platform and Control System

To verify the accuracy of the modified suspension force model and the performance of the vehicle-mounted flywheel battery by using the modified model, as shown in Figure 15a, we designed and constructed an experimental platform capable of simulating the driving state of the vehicle to be used to simulate the working state of the vehicle as well as the disturbance of the road condition, and to analyze the dynamic performance of the vehicle-mounted flywheel battery when it is in operation. Figure 16 shows the control diagram for the flywheel battery with a multi-functional air gap. It mainly includes three modules: the controller module, power drive module, and displacement detection module. The displacement sensor detects the real-time displacement of the flywheel in 5 DOF. In case of flywheel offset, the sensor transmits the eccentricity signal of the flywheel to the controller. The corresponding control signal is then output by the controller through its calculations. The control signal is also used as an input signal to the power amplifier, which outputs a drive signal to regulate the control current of the magnetic floating bearing coil. The control system adopts double closed-loop feedback control of the displacement and current, which enhances the real-time stability of the system.

For the control system that considers the DOF interference factor and adjusts control parameters via FA-BPNN, after obtaining appropriate control parameters through it, as its radial and axial control is related to the flywheel’s torsion, the torsion is taken into account in the control. Then, the original force–current stiffness is corrected by applying the corresponding interference factor in the suspension force model, and the control current is calculated based on the corrected force–current stiffness. On this basis, the three-phase inverters corresponding to torsion control and radial control drive the control currents under the action of the control signals output from the PID controller to regulate the displacement of the flywheel in the radial and torsion directions, respectively. In axial control, the axial drive control current is output from the axial switching power amplifier, and the control signal output from the PID controller adjusts the axial displacement of the flywheel.

5.2. Performance Tests

In order to verify the accuracy of the model, it is first necessary to compare the force–current stiffnesses of the accurate model with that of the fundamental model as well as the experimental model, and Figure 17a shows the comparison in the radial x-direction and Figure 17b in the axial direction. The fundamental model without considering the DOF is the one that only applies when the torsion angle is 0 (in this case the accurate model is equal to the fundamental model), in which case it has a better fit with the experimental results. However, when the flywheel is in a 0.5 deg torsion, the fundamental model does not take into account the interaction between the degrees of freedom. Therefore, the fundamental model cannot be readily adjusted with the change in torsion angle, resulting in a large difference between the suspension force and the actual situation, and the control effect is greatly weakened if the fundamental model is used. In contrast, the suspension force of the radial and axial accurate model can be adjusted in real time with the change in torsion angle, so its force–displacement stiffness and force–current stiffness are closer to the experimental results.

To further validate the superiority of the accurate modeling and its control system in real working conditions, the uphill condition is chosen as the experimental condition. The experimental conditions are as follows: the vehicle-mounted flywheel battery is controlled using both the fundamental model without control parameters adjusted and the accurate modeling of control parameters adjusted with FA-BPNN, keeping phase A as the main poles of the torsion, and rushing up the 5 deg and 10 deg slopes along the x-direction at a speed of 1 m/s. As shown in Figure 15b, the uphill angle is adjusted in real time by an adjustable-angle experiment with a ramp.

The flywheel always remains suspended due to the presence of biased magnetic flux. The sensor provides real-time feedback on the position of the flywheel to quickly adjust the control current output. From Figure 18, under the influence of the working conditions uphill along the x-direction, the flywheel creates an x-direction offset in the radial direction and remains there for a period before the vehicle completes the uphill climb. The y-direction offset of the flywheel is negligible because it goes uphill in the x-direction. When comparing Figure 18a with Figure 18b and Figure 18c with Figure 18d, it can be found that the radial offset of the control based on the suspension force model considering the DOF interaction factor and control parameters adjusted with FA-BPNN is reduced by 44% at a slope of 5 deg and 35% at a slope of 10 deg compared to the fundamental model without control parameters adjusted.

Also, comparing Figure 19a with Figure 19b, and Figure 19c with Figure 19d, the degree of axial offset using the corrected model is reduced by 28% at 5 deg uphill and 24% at 10 deg uphill. This proves that the accurate modeling of control parameters adjusted with FA-BPNN has better control efficiency.

6. Conclusions

To enhance the stability of the vehicle-mounted flywheel battery system so that it can be applied to vehicle-mounted working conditions more safely and efficiently, this paper innovatively summarizes the complex interferences among degrees of freedom as interference factors to achieve precise modeling on the basis of the original model. Subsequently, the precise model is combined with the Firefly Optimization-Back Propagation Neural Network (FA-BPNN). FA-BPNN performs excellently in multi-input and multi-output situations. It can accurately obtain control parameters and cooperate with the precise model. On the one hand, an improved precise model is obtained by summarizing and applying the complex situation of magnetic field changes under the interaction of degrees of freedom to the basic model; on the other hand, FA-BPNN is used to adjust the control parameters in real time. The combination of the two significantly improves the stability of the flywheel battery system under common working conditions. The comparison with the finite element simulation proves the accuracy of the model, and the designed experiments demonstrate the superiority of the model and the control system. In the uphill condition, the radial and axial offsets are reduced by more than 30% and 20%, respectively, compared with the control based on the basic model without adjusted control parameters.

When different flywheel batteries are selected, considering that the interference law of the degree-of-freedom interference on the magnetic bearing does not change with the change of the flywheel battery, summarizing the degree-of-freedom interference as an interference factor and applying it to the original suspension force model can still quickly and accurately improve the accuracy of the model. As for FA-BPNN, although the working condition characteristic parameters may be different for different flywheel batteries, as long as the training data sets are complete, high-efficiency control can still be achieved. That is to say, the study proposed in this paper has universality.

The structural characteristics of the flywheel selected in this paper lead to an interaction between the DOF reflected in the torsional DOF, which leads to a large change in the axial and radial suspension forces. However, for other flywheel structures and magnetic bearings, the interaction of the DOF may be more pronounced in the other DOF of the magnetic bearings, and therefore the accurate modeling approach will be different. However, the modeling ideas presented in this paper can still be used for further analysis in specific cases. At the same time, the working conditions of the vehicle in this paper only selected a few of the most representative conditions; if more conditions under the control parameters can be summarized and added to the screening mechanism of the neural network, it can achieve more accurate control.