Charging–Discharging Control Strategy for a Flywheel Array Energy Storage System Based on the Equal Incremental Principle

Abstract

The widely used flywheel energy storage (FES) system has such advantages as high power density, no environment pollution, a long service life, a wide operating temperature range, and unlimited charging–discharging times. The flywheel array energy storage system (FAESS), which includes the multiple standardized flywheel energy storage unit (FESU), is an effective solution for obtaining large capacity and high-power energy storage. In this paper, the strategy for coordinating and controlling the charging–discharging of the FAESS is studied in depth. Firstly, a deep analysis is conducted on the loss generated during the charging–discharging process of the FESU. The results indicate that the loss is related to the charging–discharging of power. To solve the problems of over-charging, over-discharging, and overcurrent caused by traditional charging–discharging control strategies, this paper proposes a charging–discharging coordination control strategy based on the equal incremental principle (EIP). This strategy aims to minimize the total loss and establish a mathematical model of optimal coordination control with the constraints of total charging–discharging power, rated power limit, over-charging, over-discharging, and overcurrent. Based on the EIP, the optimal distribution scheme of power charging–discharging is determined. Secondly, this paper gives the specific implementation scheme of the optimal coordinated control strategy. Lastly, the charging–discharging coordinated control strategy is verified by examples. The results show that the coordinated control strategy can effectively reduce the loss during the charging–discharging process and can prevent over-charging, over-discharging, and overcurrent of the system. Overall, it has a better control effect than the existing charging–discharging control strategies.

1. Introduction

To address the increasingly prominent environmental and energy problems, renewable energy sources, such as solar and wind energy, have developed rapidly in recent years. However, if renewable energy sources with characteristics of intermittence and randomness are connected to the power grid on a large scale, the voltage and frequency of the power grid will fluctuate drastically, resulting in poor continuity and stability of the power supply [1,2]. The configuration of a large capacity energy storage system in the power grid that utilizes charging and discharging functions can smooth out the power fluctuation in renewable energy, achieving decoupling control of both the power generation side and the power consumption side to improve the stability of the power grid [3]. At present, the energy storage form adopted in the power grid system is generally the storage battery, whose charging–discharging process is achieved by electrochemical reactions, with drawbacks such as slow speed, a limited recycling time, and a short service life. At the same time, the replacement of batteries will significantly increase operating costs [4]. In addition, wasted batteries pollute the environment.

Flywheel energy storage (FES) is a form of energy storage that uses a high-speed rotating flywheel rotor as a carrier to convert electrical energy into mechanical energy. It has the advantages of high power density, no environmental pollution, a long service life, and an almost unlimited charging–discharging time. Therefore, FES has been widely used in various fields, including renewable energy power generation, rail brake energy recovery, uninterrupted power supply (UPS), frequency modulation and voltage regulation in power systems, satellite energy storage, and attitude adjustment [5,6]. There are two main technical approaches used to acquire a larger storage capacity, higher charging–discharging power, and longer backup time. One is to develop a large-capacity single-unit FES system, and the other one is to connect multiple standardized flywheel energy storage units (FESUs) in parallel to form a flywheel array energy storage system (FAESS). The single large-capacity FES, however, not only increases the cost significantly but is also difficult to use due to technical limitations [3]. Comparatively, the FAESS can decrease both the cost and the time required for research and development; accordingly, it is a better solution for obtaining a high-capacity and high-power energy storage system.

Few studies have investigated the optimization of coordination control of the FAESS, both at home and abroad. In regard to charging, based on the average distribution algorithm, [7] proposes a distribution algorithm according to the chargeable capacity. As for discharging, discharging power distribution algorithms (average distribution, distribution by speed ratio, and distribution by residual energy ratio) are proposed in [7,8,9,10], and a preliminary simulation study was conducted. In terms of the FAESS used in wind farms, a layered distribution control strategy is presented in [11]: The upper layer is used to allocate energy to each FESU based on the proportional consistency algorithm, and the lower layer relies on the local controller in the FESU to adjust the rotating speed. This strategy, however, does not consider the effect of loss on the charging–discharging of the system. In [12], the FAESS is applied in the field of urban rail transit to absorb braking energy, and it achieves a good control effect; however, its power distribution scheme has not been studied in depth. The study presented in [13] advances a comprehensive electromechanical model of FES system, analyzes small-signal stability, and conducts time domain simulation research. The results show that the FAESS improves the stability of grid system, but the internal power distribution has not been studied in depth. In [14,15,16], an optimal control method aimed at FAESS is proposed for the grid-connected micro-grid. The upper optimization center establishes a corresponding charging–discharging optimized model according to the power vacancy and rotating speed of each flywheel, and solves the power parameter of each flywheel. The lower flywheel controller adopts a dual-module, dual-ring control method to control the rotating speed and output power of the flywheel. This method is still, essentially, a distribution by residual power ratio.

In general, the existing coordinated operation control of the FAESS mostly adopts the algorithms of average distribution, distribution by speed ratio, or residual energy ratio. These algorithms are likely to cause the distributed power to be greater than the available output power of the FESU, further causing overcurrent of the flywheel motor. At the same time, over-charging or over-discharging is likely to occur. In addition, the existing algorithms do not consider the internal loss factor of the system during the actual operation, resulting in low system operation efficiency. In order to solve these problems, further research on the FAESS is needed.

In this paper, the distribution of the charging–discharging power of the FAESS is studied in depth. Firstly, this paper conducts an in-depth analysis of the power loss of each component of the FESU, including the loss from the grid-side converter, the generator-side converter, and the motor. Then, with minimum loss as the optimization goal, an objective function of optimal control is established and solved by the equal incremental principle (EIP), with constraints on the total charging–discharging power, rated power limit, overcharging, over-discharging, and overcurrent. Secondly, based on theoretical research, a specific implementation scheme for the charging–discharging control strategy is proposed. Finally, the charging–discharging control strategy based on EIP is analyzed and verified by examples. The results demonstrate the effectiveness of the control strategy.

2. Loss Analysis of FESU

During charging, the FESU absorbs energy from the grid, whereas during discharge, it outputs energy into the grid. As shown in Figure 1, the converter generates switching and conduction loss during the power conversion process. During the electromechanical conversion process, the loss generated by the motor includes copper loss, iron loss (hysteresis loss, eddy current loss), and friction loss (bearing friction loss, wind resistance friction loss). When the FESU is connected to the grid for discharging, the loss distribution is as shown in Figure 2. The mechanical energy stored in the flywheel rotor is first converted into electrical energy through electromechanical conversion, and is then transferred into the power grid through the converter. Similar to that, during the charging process, the motor loss during the electromechanical conversion process includes copper loss, iron loss, and friction loss. During the process of power conversion, there are also switching and conduction losses.

2.1. Loss Analysis of the Converter in the FESU

The converter loss in the FESU includes the loss from the grid-side and the motor-side converters. As shown in Figure 3, taking the voltage source converter (VSC) as an example, the grid-side converter contains six insulated gate bipolar transistors (IGBTs) and six anti-parallel diodes, and the loss mainly includes switching and conduction loss. The switching loss in the grid-side converter is as shown in (1):

$$\begin{array}{l}{P}_{sw\_G}=6\times (P_{sw\_IGBT}+P_{sw\_D})\\ =\frac{6}{\pi }\times f_{sw}\times (E_{on\_test}+E_{off\_test}+E_{rec\_test})\times \frac{I_m}{I_{test}}\times \frac{U_{dc}}{V_{test}}\end{array}$$

(1)

where $P_{sw\_G}$ represents the switching loss of the grid-side converter; $P_{sw\_IGBT}$ represents the switching loss of a single IGBT; $P_{sw\_D}$ represents the switching loss of a single diode; $E_{on\_test}$ represents the switching-on loss of the IGBT; $E_{off\_test}$ represents the switching-off loss of the IGBT; $E_{rec\_test}$ represents the switching-off loss of the diode; $I_{test}$ represents the test current; $V_{test}$ represents the test voltage; $I_m$ represents the output current of the grid-side converter; $U_{dc}$ represents the DC bus voltage; and $f_{sw}$ represents the switching frequency of the grid-side converter.

Since the conduction loss of the negative half-cycle of the grid-side converter is the same as that of the positive half-cycle, we can only consider the conduction loss of the positive half-cycle. In the one-phase leg of the grid-side converter, the conduction loss of the IGBT is

$$P_{IGBT}=V_{T0}·I_m(\frac{1}{\pi }+\frac{1}{4}m\cos \phi )+R_{CE}·I_m{}^2·(\frac{1}{4}+\frac{2}{3\pi }m\cos \phi )$$

(2)

where $P_{IGBT}$ represents the conduction loss of the IGBT in the one-phase leg of the grid-side converter; $V_{T0}$ represents the IGBT threshold voltage; $I_m$ represents the output phase current in the grid-side converter; $m$ represents the modulation ratio; $\phi$ represents the power factor angle; and $R_{CE}$ represents the IGBT conduction resistance.

By calculating the diode conduction loss, we can acquire

$$P_D=V_{D0}·I_m(\frac{1}{\pi }-\frac{1}{4}m\cos \phi )+R_D·I_m{}^2·(\frac{1}{4}-\frac{2}{3\pi }m\cos \phi )$$

(3)

where $P_D$ represents the conduction loss of the diode in the one-phase leg of the grid-side converter; $V_{D0}$ represents the diode threshold voltage; $I_m$ represents the output phase current in the grid-side converter; $m$ represents the modulation ratio; $\phi$ represents the power factor angle; and $R_D$ represents the diode conduction resistance.

The total conduction loss of grid-side converter $P_{cond}$ is

$$\begin{array}{ll}{P}_{cond}& =3(P_{IGBT}+P_D)\\ & =3I_m·(\frac{V_{T0}+V_{D0}}{\pi }+m\cos \phi \frac{V_{T0}-V_{D0}}{4})+\\ & 3I_m{}^2·\left[\frac{R_{CE}+R_D}{4}+m\cos \phi \frac{2(R_{CE}-R_D)}{3\pi }\right]\end{array}$$

(4)

The absolute value of $m\cos \phi$ is less than 1, and

$$\left|m\cos \phi \frac{V_{T0}-V_{D0}}{4}\right|<<\frac{V_{T0}+V_{D0}}{\pi }$$

(5)

$$\left|m\cos \phi \frac{2(R_{CE}-R_D)}{3\pi }\right|<<\frac{R_{CE}+R_D}{4}$$

(6)

Therefore, the total conduction loss can be simplified to

$$P_{cond}=\frac{3(V_{T0}+V_{D0})}{\pi }{I}_m+\frac{3(R_{CE}+R_D)}{4}{I}_m{}^2.$$

(7)

According to (1) and (7), the total conduction loss in the grid-side $P_{G\_tot}$ is

$$P_{G\_tot}=b{I}_m+c{I}_m{}^2$$

(8)

$$b=\frac{6}{\pi }\times f_{sw}\times (E_{on\_test}+E_{off\_test}+E_{rec\_test})\times \frac{1}{I_{test}}\times \frac{U_{dc}}{V_{test}}+\frac{3(V_{T0}+V_{D0})}{\pi }$$

$$c=\frac{3(R_{CE}+R_D)}{4}.$$

Since the relation between the amplitude of output phase current $I_m$ and the active power directive value $P^{*}$ is

$$I_m=\frac{P^{*}}{\sqrt{3}{U}_N}\times \sqrt{2},$$

(9)

the total loss in the grid-side converter $P_{G\_tot}$ is

$$P_{G\_tot}=d\left|P^{\ast }\right|+f{P}^{{\ast }^2}$$

(10)

where

$$\begin{array}{l}d=b\times \frac{\sqrt{2}}{\sqrt{3}\ast U_N}\\ =\left[\frac{6}{\pi }\times f_{sw}\times (E_{on\_test}+E_{off\_test}+E_{rec\_test})\times \frac{1}{I_{test}}\times \frac{U_{dc}}{V_{test}}+\frac{3(V_{T0}+V_{D0})}{\pi }\right]\times \frac{\sqrt{2}}{\sqrt{3}{U}_N}\end{array}$$

$$f=c\times \frac{2}{3U_N{}^2}=\frac{3(R_{CE}+R_D)}{4}\times \frac{2}{3U_N{}^2}.$$

Taking the VSC converter as an example for the motor-side converter, the loss mechanism is exactly the same. In the dq coordinate system, the d axis current $i_d$ is set to be 0, and then the amplitude of the output phase current is the absolute value $\left|i_q\right|$ of the q axis current $i_q$. Similarly, the total loss in the motor-side converter $P_{M\_tot}$ is

$$P_{M\_tot}=b\left|i_q\right|+c{i}_q{}^2.$$

(11)

2.2. Motor Loss of FESU

As shown in Figure 2 and Figure 3, the motor loss in the FESU is mainly copper loss, iron loss, and friction loss. The cooper loss $P_{cu}$ in the motor is

$$P_{cu}=\frac{3}{2}{R}_s{i}_q{}^2=g{i}_q{}^2$$

(12)

where

$$g=3\raisebox{1ex}{$R_s$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$

where $R_s$ is the resistance of the stator in the motor.

The iron loss of the permanent magnet synchronous motor includes hysteresis loss and eddy current loss. The iron loss can be represented by the equivalent iron loss resistor $R_c$. The equivalent steady state circuit of the permanent magnet synchronous motor that includes the iron loss resistor $R_c$ is shown in Figure 4.

From Figure 4, we can acquire (13):

$$\{\begin{array}{l}{i}_{cd}=\frac{-{\omega }_e{L}_q{i}_{oq}}{R_c}=-k{L}_q{i}_{oq}\\ i_{cq}=\frac{{\omega }_e(L_d{i}_{od}+{\psi }_f)}{R_c}=k(L_d{i}_{od}+{\psi }_f)\end{array}$$

(13)

where $i_{cd}$ represents the d axis equivalent current of the equivalent iron loss current $i_c$; $i_{cq}$ represents the q axis equivalent current of the equivalent iron loss current $i_c$; ${\omega }_e$ represents the electric angular speed of the motor; ${\psi }_f$ represents the excitation flux linkage of the motor permanent magnet; $L_d$ represents the d axis equivalent inductance of the motor; $L_q$ represents the q axis equivalent inductance of the motor; $i_{od}$ represents the d axis equivalent current of the motor; and $i_{oq}$ represents the q axis equivalent current of the motor.

According to (13), the iron loss of motor $P_{Fe}$ is

$$P_{Fe}=\frac{3}{2}(i_{cd}^2+i_{cq}^2)R_c=\frac{3}{2}kp{\omega }_m[{(L_d{i}_{od}+{\psi }_f)}^2+{(L_q{i}_{oq})}^2].$$

(14)

From

$$\{\begin{array}{l}{i}_d=i_{cd}+i_{od}=-k{L}_q{i}_{oq}+i_{od}=0\\ i_q=i_{cq}+i_{oq}=k(L_d{i}_{od}+{\psi }_f)+i_{oq}\end{array}$$

(15)

we can acquire

$$\{\begin{array}{l}{i}_{od}=\frac{k{L}_q(i_q-k{\psi }_f)}{1+k^2{L}_d{L}_q}\\ i_{oq}=\frac{i_q-k{\psi }_f}{1+k^2{L}_d{L}_q}\end{array}.$$

(16)

Bringing (16) into (14), we can acquire

$$P_{Fe}=(k_1{i}_q{}^2+k_2{i}_q+k_3)\times {\omega }_m$$

(17)

where ${\omega }_m$ is the mechanical angular speed of the flywheel rotor and

$$\begin{array}{l}{k}_1=\frac{3}{2}\frac{kp(1+k^2{L}_d{}^2)L_q{}^2}{{(1+k^2{L}_d{L}_q)}^2}\\ k_2=\frac{3k^{2}p{\psi }_f(L_d-L_q)L_q}{{(1+k^2{L}_d{L}_q)}^2}\\ k_3=\frac{3}{2}\frac{kp(1+k^2{L}_q{}^2){\psi }_f{}^2}{{(1+k^2{L}_d{L}_q)}^2}\end{array}$$

The friction loss of the flywheel $P_{fr}$ can be represented by

$$P_{fr}=B\times {\omega }_m{}^2$$

(18)

where $B$ represents the viscous friction coefficient and ${\omega }_m$ represents the mechanical angular speed of the flywheel rotor.

2.3. Analysis of the Total Charging–Discharging Loss of the FESU

Figure 5 shows the energy flowing when the FESU is charging. $P_C$ is the absolute value of the directive value of active power during charging; $P_{dcin}$ is the DC input power of the generator-side converter; $P_{acin}$ is the AC output power of the generator-side converter; $P_M$ is the electromagnetic power of the flywheel motor; and $P_F$ is the actual storage power of the flywheel. During charging in the grid, the total loss $P_{lossin}$ of FESU is

$$P_{lossin}=\alpha P_C^2+\beta P_C+\gamma$$

(19)

where

$$\alpha =f+(c+g+k_1{\omega }_m)k_{\omega }{}^2$$

$$\beta =d+(b+k_2{\omega }_m)k_{\omega }$$

$$\gamma =k_3{\omega }_m+B{\omega }_m{}^2.$$

The energy flowing in the FESU during discharge is shown is Figure 6. $P_D$ is the absolute value of the directive value of active power during discharge; $P_{dcout}$ is the DC output power of the motor-side converter; $P_{acout}$ is the AC input power of the motor-side converter; $P_M$ is the electromagnetic power of the flywheel motor; and $P_F$ is the actual output power of the flywheel. During charging of the grid, the total loss $P_{lossout}$ of FESU is

$$P_{lossout}={\alpha }^{\prime }{P}_D{}^2+{\beta }^{\prime }{P}_D+{\gamma }^{\prime }$$

(20)

where

$${\alpha }^{\prime }=f+(c+g+k_1{\omega }_m){k^{\prime }}_{\omega }{}^2$$

$${\beta }^{\prime }=d+(-b+k_2{\omega }_m){k^{\prime }}_{\omega }$$

$${\gamma }^{\prime }=k_3{\omega }_m+B{\omega }_m{}^2.$$

It can be seen from the above analysis that when the FESU is charged or discharged, the generated loss is related to power charging or discharging. In the case that the total charging or discharging power is fixed, the charging or discharging power that distributes to each unit is different and so is the generated loss and the total loss in the system. Therefore, it is necessary to conduct in-depth research on the charging–discharging coordination control strategy of the FAESS.

3. Charging Control Strategy of the FAESS

3.1. Charging Objective Function for the FAESS

When the total charging power of the system is fixed, in order to minimize the total loss $P_{loss\_tot}$ during the charging operation, the objective function is as follows:

$$P_{loss\_tot}={\displaystyle \sum _{i=1}^N{P}_{loss\_i}}(P_{C\_i})$$

(21)

$$P_{loss\_i}(P_{C\_i})={\alpha }_i{P}_{C\_i}^2+{\beta }_i{P}_{C\_i}+{\gamma }_i(i=1,2,\dots ,N)$$

(22)

where $P_{loss\_tot}$ represents the total loss of the FAESS during charging; $P_{loss\_i}$ represents the loss of FESU i during charging; $P_{C\_i}$ represents the charging power reference value of FESU i; and

$${\alpha }_i=f+(c+g+k_1{\omega }_{m\_i})k_{\omega \_i}{}^2$$

$${\beta }_i=d+(b+k_2{\omega }_{m\_i})k_{\omega \_i}$$

$${\gamma }_i=k_3{\omega }_{m\_i}+B{\omega }_{m\_i}{}^2$$

$$k_{\omega \_i}=\frac{1-d}{b+h{\omega }_{m\_i}}.$$

3.2. Charging Constraints of FAESS

3.2.1. Equity Constraint

The sum of the total charging power of FESU is equal to the charging power reference value of the system $P_{C\_array}$, which is

$${\displaystyle \sum _{i=1}^N{P}_{C\_i}}=P_{C\_array}.$$

(23)

3.2.2. Inequity Constraints

(1) Inequity constraints of the limited amplitude of the rated power

The charging power of each FESU cannot exceed its rated power, which is

$$0\le P_{C\_i}\le P_{rated}(i=1,2,\dots ,N).$$

(24)

(2) Inequity constraints of overcharging prevention

The current rotating speed of FESU i is set as ${\omega }_{{}_{m\_i}}$, and the kinetic energy is set as $E_i$. After operating cycle $\Delta T$, the rotating speed is changed to ${\omega }_{{}_{m\_i\_next}}$ and the kinetic energy is changed to $E_{{}_{i\_next}}$. In order to prevent overcharging, the rotating speed of the next moment ${\omega }_{{}_{m\_i\_next}}$ should not exceed the maximum permitted rotating speed ${\omega }_{{}_{max}}$, as shown in Figure 7.

The kinetic energy of the flywheel energy storage unit at the next moment $E_{i\_next}$ should not exceed the maximum allowable kinetic energy $E_{max}$ (25).

$$E_{i\_next}=E_i+P_{F\_i}\times \Delta T\le E_{max}$$

(25)

where $P_{F\_i}$ is the current input power of i of the flywheel rotor.

Equation (25) can be represented by

$$P_{F\_i}\le P_{C\_i}-d{P}_{C\_i}-f{P}_{C\_i}^2-B{\omega }_{m\_i}^2.$$

(26)

According to (26),

$$E_i+P_{F\_i}\times \Delta T\le E_i+(P_{C\_i}-d{P}_{C\_i}-f{P}_{C\_i}^2-B{\omega }_{m\_i}^2)\times \Delta T.$$

(27)

It can be seen from (27) and (25) that if (28) is established, (25) is also established:

$$E_i+(P_{C\_i}-d{P}_{C\_i}-f{P}_{C\_i}^2-B{\omega }_{m\_i}^2)\times \Delta T\le E_{max}.$$

(28)

From (28), (29) is acquired:

$$\frac{1}{2}J{\omega }_{m\_i}^2+(P_{C\_i}-d{P}_{C\_i}-f{P}_{C\_i}^2-B{\omega }_{m\_i}^2)\times \Delta T\le \frac{1}{2}J{\omega }_{max}^2.$$

(29)

which is

$$f·\Delta T·P_{C\_i}^2+l·\Delta T·P_{C\_i}+\frac{1}{2}J{\omega }_{max}^2-\frac{1}{2}J{\omega }_{m\_i}^2+B·\Delta T·{\omega }_{m\_i}^2\ge 0$$

(30)

where $l=d-1$ and

$$A_i=\frac{1}{2}J{\omega }_{max}^2-\frac{1}{2}J{\omega }_{m\_i}^2+B·\Delta T·{\omega }_{m\_i}^2.$$

(31)

Then,

$$f·\Delta T·P_{C\_i}^2+l·\Delta T·P_{C\_i}+A_i\ge 0.$$

(32)

According to (32),

$$P_{C\_i}\le \frac{-l·\Delta T-\sqrt{{(l·\Delta T)}^2-4A_{i}f·\Delta T}}{2f·\Delta T}=P_{OCh\_i\_max}(i=1,2,\dots ,N)$$

(33)

where for $A_i$, as in (31), $l=d-1$.

Equation (33) is the constraint for preventing the overcharging of FESU i.

(3) Inequity constraints of preventing overcurrent in the motor

If the distributed charging power is excessive, overcurrent can occur in the flywheel motor when the FESU is connected to the grid and charging. The q axis current of flywheel motor $i$ is set as $i_{q\_i}$. After operating cycle $\Delta T$, the q axis current in the next moment will be $i_{q\_i\_next}$. Since the charging power is constant and the rotating speed is increasing within operating cycle $\Delta T$, there will be $i_{q\_i\_next}<i_{q\_i}$, as shown in Figure 8. In order to prevent the overcurrent of flywheel motor $i$, the q axis current $i_{q\_i}$ cannot exceed the maximum permitted q axis current $i_{q\_max}$. Thus, the q axis current will naturally not exceed $i_{q\_max}$, so (34) should be guaranteed:

$$i_{q\_i}\le i_{q\_max}.$$

(34)

According to (34),

$$i_{q\_i}=k_{\omega \_i}{P}_{C\_i}\le i_{q\_max}$$

(35)

and then,

$$P_{C\_i}\le i_{q\_max}/k_{\omega \_i}=P_{OCur\_i\_max}$$

(36)

where

$$k_{\omega \_i}=(1-d)/(b+h{\omega }_{m\_i})$$

$$h=\raisebox{1ex}{$3p{\psi }_f$}\!\left/ \!\raisebox{-1ex}{$2$}\right..$$

Equation (38) is the constraint for preventing the overcurrent of flywheel motor $i$.

According to (24), (33), and (36), the upper limit of charging power should be the minimum of these three values, which is

$$P_{C\_i\_max}=\min (P_{rated},P_{OCh\_i\_max},P_{OCur\_i\_max})(i=1,2,\dots ,N).$$

(37)

Comprehensively, when the FAESS is charging, the charging power of each unit should meet the following inequality constraints:

$$0\le P_{C\_i}\le P_{C\_i\_max}(i=1,2,\dots ,N).$$

(38)

3.3. Research of the Charging Control Method for FAESS

The charging power distribution of FAESS can be presented with the equality constraint of

$${\displaystyle \sum _i^N{P}_{C\_i}}=P_{C\_array}.$$

(39)

Let the objective function:

$$P_{loss\_tot}={\displaystyle \sum _{i=1}^N{P}_{loss\_i}(P_{C\_i})}$$

(40)

be the minimum, where $p_{loss\_i}(P_{C\_i})={\alpha }_i{P}^2{}_{C\_i}+{\beta }_i{P}_{C\_i}+{\gamma }_i$ represents the loss characteristics of each FESU during charging.

At the same time, the charging power of each unit $P_{C\_i}$ should also meet the following inequity constraints (38). When the minimum loss of the charging operation of the system is calculated by EIP, the LaGrange function shown in (41) should be constructed according to Objectives (39) and (40):

$$L=P_{loss\_tot}-\lambda ({\displaystyle \sum _{i=1}^N{P}_{C\_i}}-P_{C\_array})$$

(41)

where $\lambda$ represents the LaGrange multiplier.

The necessary condition for LaGrange function L to take the extreme value is

$$\{\begin{array}{l}\frac{\partial L}{\partial P_{C\_i}}=\frac{\partial P_{loss\_tot}}{\partial P_{C\_i}}-\lambda =0(i=1,2,\dots ,N)\\ {\displaystyle \sum _{i=1}^N{P}_{C\_i}}-P_{C\_array}=0\end{array}.$$

(42)

The charging loss $P_{loss\_i}(P_{C\_i})$ of each unit is only the function that relates to its charging power $P_{C\_i}$; therefore,

$$\frac{\partial p_{loss\_tot}}{\partial P_{C\_i}}=\frac{d{p}_{loss\_i}(P_{C\_i})}{d{P}_{C\_i}}=2{\alpha }_i{P}_{C\_i}+{\beta }_i={\lambda }_i(i=1,2,\dots ,N)$$

(43)

where ${\lambda }_i$ represents the incremental rate of consumption of each FESU.

Equation (42) can be rewritten as

$$\{\begin{array}{l}2{\alpha }_i{P}_{C\_i}+{\beta }_i=\lambda (i=1,2,\dots ,N)\\ {\displaystyle \sum _{i=1}^N{P}_{C\_i}}-P_{C\_array}=0\end{array}.$$

(44)

There are $N+1$ equations and $N+1$ variables in (44), which means the value of $P_{C\_i}$ solved by (44) is the value of distributed charging power with minimum loss.

When the charging power of FESU $i$ is lower than the lower limit or higher than the upper limit, the charging power of the system will be zero or $P_{C\_i\_max}$. Then, the residual charging power will be distributed to other FESUs.

3.4. Implementation Scheme of Charging Control for FAESS

After the calculation of the distributed power value corresponding to the minimum loss based on the equal incremental principle, the directive charging order will be sent to each underlying controller to realize coordinated charging control of each unit. As shown in Figure 9, the implementation consists of an array controller and N bottom controllers. The array controller runs the charging coordinated control strategy based on the EIP and outputs the most optimal charging power value corresponding to the minimum loss; each bottom controller will control the charging process of each unit by following the distributed charging power.

The main steps of the coordination control strategy based on EIP are as follows:

(1) The directive value of the total charging power $P^{*}{}_{array}$ and the rotating speed of each flywheel ${\omega }_{m\_i}$ are sampled and systematic parameters are invoked.

(2) The upper limit of charging power is calculated to prevent overcharging and overcurrent of FESU, $P_{OCh\_i\_max}$, and $P_{OCur\_i\_max}$.

(3) The upper charging power limit of each FESU, $P_{C\_i\_max}$, is determined according to (36).

(4) Coefficients, ${\alpha }_i$, ${\beta }_i$, and ${\gamma }_i$ are calculated using the quadratic equation of charging loss of FESU, $P_{{loss\_i(P}_{C\_i})}$.

(5) The minimum and maximum values of the incremental ratio of charging consumption of each FESU, ${\lambda }_i$, are determined:

$$\{\begin{array}{c}{\lambda }_{i\_min}={\beta }_i\\ {\lambda }_{max}==2{\alpha }_i{P}_{C\_i\_max}+{\beta }_i\end{array}.$$

(45)

(6) The minimum and maximum values of the incremental ratio of public consumption of the FAESS are determined:

$$\{\begin{array}{c}{\lambda }_{min}=\min ({\lambda }_{1\_min},{\lambda }_{2\_min},\cdots ,{\lambda }_{N\_min})\\ {\lambda }_{max}=\max ({\lambda }_{1\_max},{\lambda }_{2\_max},\cdots ,{\lambda }_{N\_max})\end{array}.$$

(46)

(7) The incremental ratio of public consumption is set as ${\lambda }^{\left(0\right)}={\lambda }_{max}$.

(8) The charging power $P_{C\_i}^{(0)}$ corresponding to ${\lambda }^{\left(0\right)}$ is calculated using

$$P_{C\_i}^{(0)}=\frac{{\lambda }^{(0)}-{\beta }_i}{2{\alpha }_i}(i=1,2,\dots ,N).$$

(47)

If $P_{C\_i}^{(0)}$ is over the limit, then the limit value is used.

(9) The deviation between the sum of charging power and the total charging power $P_{C\_array}$ corresponding to ${\lambda }^{\left(0\right)}$ is calculated:

$$\Delta P_C^{(0)}={\displaystyle \sum _{i=1}^N{P}_{C\_i}^{(0)}-P_{C\_array}}.$$

(48)

(10) Verification of whether the power deviation will satisfy the permitted range is carried out using

$$\Delta P_C^{(0)}\le \epsilon .$$

(49)

If (49) is satisfied, then $P_{C\_i}^{(0)}$ is the value of distributed charging power with minimum loss, then $P_i^{\ast }$ is given a negative value of $P_{C\_i}^{(0)}$ as the directive value of the charging power of the FESU.

(11) If it is not satisfied, and ${\lambda }^{\left(0\right)}>{\lambda }_{min}$, then

$${\lambda }^{\left(1\right)}>{\lambda }^{\left(0\right)}-\Delta \lambda .$$

(50)

The iterative computations start from step (8) until the power deviation satisfies the permitted range, or until the incremental ratio of public consumption is lower than or equal to ${\lambda }_{min}$.

The charging coordination control program based on EIP is shown in Figure 10.

4. Discharging Control Strategy of the FAESS

4.1. Discharging Objective Function for the FAESS

When the FAESS is discharging, the loss $P_{loss\_i}$ of FESU $i$ is related to its active discharging power $P_{D\_i}$. In the case that the total discharging power is fixed, the distributed power of FESU $i$ $P_{D\_i}$ is different, and so is its loss $P_{loss\_i}$ and the total loss of the system. With a fixed total discharging power, in order to minimize the total loss $P_{loss\_tot}$, the objective function is defined as follows:

$$p_{loss\_tot}={\displaystyle \sum _{i=1}^N{p}_{loss\_i}}(P_{D\_i})$$

(51)

$$P_{loss\_i}(P_{D\_i})={\alpha }_i^{\prime }{P}_{D\_i}^2+{\beta }_i^{\prime }{P}_{D\_i}+{\gamma }_i^{\prime }(i=1,2,\dots ,N)$$

(52)

where

$${\alpha }_i^{\prime }=f+(c+g+k_1{\omega }_{m\_i})k_{\omega \_i}^{\prime }{}^2$$

$${\beta }_i^{\prime }=d+(-b+k_2{\omega }_{m\_i})k_{\omega \_i}^{\prime }$$

$${\gamma }_i^{\prime }=k_3{\omega }_{m\_i}+B{\omega }_{m\_i}{}^2$$

$$k_{\omega \_i}^{\prime }=\frac{1+d}{b-h{\omega }_{m\_i}}.$$

4.2. Discharging Constraints of the FAESS

4.2.1. Equity Constraints

When the system is discharged, the equity constraint is that the sum of the discharging power of FESUs is equal to the total discharging power of the system $P_{D\_array}$, which is

$${\displaystyle \sum _{i=1}^N{P}_{D\_i}}=P_{D\_array}.$$

(53)

4.2.2. Inequity Constraints

(1) Inequity constraint of the limited amplitude of the rated power

The discharging power of each unit cannot exceed its rated power, which is

$$0\le P_{D\_i}\le P_{rated}(i=1,2,\dots ,N).$$

(54)

(2) Inequity constraint of preventing over-discharging

The current rotating speed of FESU i is set as ${\omega }_{m\_i}$ and the kinetic energy is set as $E_i$. After operating cycle $\Delta T$, the rotating speed is changed to ${\omega }_{{}_{m\_i\_next}}$ and the kinetic energy is changed to $E_{i\_next}$. In order to prevent over-discharging, the rotating speed of the next moment ${\omega }_{{}_{m\_i\_next}}$ should not be lower than the minimum permitted rotating speed ${\omega }_{{}_{min}}$, as shown in Figure 11.

The kinetic energy of FESU in the next moment $E_{i\_next}$ should be higher than its minimum limit $E_{min}$, and Equation (55) should be guaranteed:

$$E_{i\_next}=E_i-P_{F\_i}\times \Delta T\ge E_{min}.$$

(55)

In (55), $P_{F\_i}$ represents the current output power of flywheel rotor $i$:

$$P_{F\_i}\ge P_{D\_i}+d{P}_{D\_i}+f{P}_{D\_i}^2+B{\omega }_{m\_i}^2.$$

(56)

According to (55) and (56),

$$E_i\ge E_{min}+(P_{D\_i}+d{P}_{D\_i}+f{P}_{D\_i}^2+B{\omega }_{m\_i}^2)\times \Delta T$$

(57)

which is

$$\frac{1}{2}J{\omega }_{m\_i}^2\ge \frac{1}{2}J{\omega }_{min}^2+(P_{D\_i}+d{P}_{D\_i}+f{P}_{D\_i}^2+B{\omega }_{m\_i}^2)\times \Delta T.$$

(58)

From (58), (59) is acquired:

$$f\times \Delta T\times P_{D\_i}^2+e\times \Delta T\times P_{D\_i}+A_i^{\prime }\le 0,$$

(59)

and (60) can be acquired from (59):

$$P_{D\_i}\le \frac{-e\times \Delta T+\sqrt{{(e\times \Delta T)}^2-4A_i^{\prime }f\times \Delta T}}{2f\times \Delta T}=P_{ODch\_i\_max}$$

(60)

where

$$A_i^{\prime }=B\times \Delta T\times {\omega }_{m\_i}^2+\frac{1}{2}J{\omega }_{min}^2-\frac{1}{2}J{\omega }_{m\_i}^2$$

$$e=1+d.$$

Equation (60) is the constraint for the prevention of over-discharge of $i$ FESU.

(3) Inequity constraints of preventing overcurrent of motor

If the distributed charging power is excessive, overcurrent of the flywheel motor can occur when FESU is connected to the grid and discharging. The q axis current of flywheel motor $i$ as $i_{q\_i}$ is set. After operating cycle $\Delta T$, the q axis current in the next moment will be $i_{q\_i\_next}$. Since it is within operating cycle $\Delta T$, where the discharging power is constant and the rotating speed is decreasing, the absolute value of the q axis current at the next moment will be greater than the absolute value of the q axis current. In order to prevent overcurrent in flywheel motor $i$, the absolute value $\left|i_{q\_i\_next}\right|$ of the q axis current in the next moment cannot exceed the maximum permitted q axis current $i_{q\_max}$, as shown in Figure 12.

During discharging, there is

$$\frac{1}{2}J{\omega }_{m\_i\_next}^2\le E_i-P_{D\_i}\times \Delta T\Rightarrow \frac{1}{{\omega }_{m\_i\_next}}\ge \sqrt{\frac{J}{2(E_i-P_{D\_i}\times \Delta T)}}$$

(61)

and

$$h{\omega }_{m\_i\_next}\left|i_{q\_i\_next}\right|\ge P_{D\_i}\Rightarrow \left|i_{q\_i\_next}\right|\ge \frac{P_{D\_i}}{h{\omega }_{m\_i\_next}}.$$

(62)

From (61) and (62) we can acquire

$$\left|i_{q\_i\_next}\right|\ge \frac{P_{D\_i}}{h}\times \sqrt{\frac{J}{2(E_i-P_{D\_i}\times \Delta T)}}$$

(63)

and

$$i_{q\_max}\ge \left|i_{q\_i\_next}\right|.$$

(64)

Therefore,

$$i_{q\_max}\ge \frac{P_{D\_i}}{h}\times \sqrt{\frac{J}{2(E_i-P_{D\_i}\times \Delta T)}}$$

(65)

which is

$$J{P}_{D\_i}^2+2i^2{}_{q\_max}{h}^2\times \Delta T\times P_{D\_i}-i^2{}_{q\_max}{h}^{2}J{\omega }_{m\_i}^2\le 0.$$

(66)

From (66), we can acquire

$$P_{D\_i}\le \frac{-i^2{}_{q\_max}{h}^2\times \Delta T+\sqrt{i^4{}_{q\_max}{h}^4\Delta T^2+J^2{i}^2{}_{q\_max}{h}^2{\omega }_{m\_i}^2}}{J}=P_{OCur\_i\_max}^{\prime }(i=1,2,\dots ,N).$$

(67)

Equation (67) is the constraint for preventing the overcurrent of FESU $i$.

According to (54), (60), and (67), the upper limit of the discharging power of FESU $i$ should be the minimum of these three equations, which is

$$P_{D\_i\_max}=\min (P_{rated},P_{ODch\_i\_max},P_{OCur\_i\_max}^{\prime }).$$

(68)

Generally, when the FAESS is discharged, the discharging power of each unit should meet the following inequity constraints:

$$0\le P_{D\_i}\le P_{D\_i\_max}(i=1,2,\dots ,N).$$

(69)

4.3. Research of the Discharging Control Method for the FAESS

The discharging power distribution of the system can be presented as follows. There are the equality constraints:

$${\displaystyle \sum _{i=1}^N{P}_{D\_i}}=P_{D\_array},$$

(70)

and the objective function

$$P_{loss\_tot}={\displaystyle \sum _{i=1}^N{P}_{loss\_i}}(P_{D\_i})$$

(71)

is the minimum.

$P_{loss\_i}(P_{D\_i})=a^{\prime }{}_i{P}^2{}_{D\_i}+{\beta }^{\prime }{}_i{P}_{D\_i}+{\gamma }^{\prime }{}_i$ represents the loss characteristics of each FESU during discharge. At the same time, the discharging power of each unit $P_{D\_i}$ should satisfy Equation (69).

When the minimum discharging loss of the system is calculated by EIP, a LaGrange function should firstly be constructed according to Objectives (70) and (71).

$$L=p_{loss\_tot}-\lambda ({\displaystyle \sum _{i=1}^N{P}_{D\_i}}-P_{D\_array})$$

(72)

where $\lambda$ is the LaGrange multiplier.

The necessary conditions for LaGrange function L at the extremum are

$$\{\begin{array}{l}\frac{\partial L}{\partial P_{D\_i}}=\frac{\partial P_{loss\_tot}}{\partial P_{D\_i}}-\lambda =0(i=1,2,\dots ,N)\\ {\displaystyle \sum _{i=1}^N{P}_{D\_i}}-P_{D\_array}=0\end{array}.$$

(73)

Since the charging loss $P_{loss\_i}(P_{D\_i})$ of each unit is only the function that relates to the discharging power $P_{D\_i}$ of each unit,

$$\frac{\partial P_{loss\_tot}}{\partial P_{D\_i}}=\frac{d{P}_{loss\_i}(P_{D\_i})}{d{P}_{D\_i}}=2{\alpha }_i^{\prime }{P}_{D\_i}+{\beta }_i^{\prime }={\lambda }_i(i=1,2,\dots ,N)$$

(74)

where ${\lambda }_i$ represents the incremental rate of consumption of each FESU.

Equation (73) can be rewritten as

$$\{\begin{array}{l}2{\alpha }_i^{\prime }{P}_{D\_i}+{\beta }_i^{\prime }=\lambda (i=1,2,\dots ,N)\\ {\displaystyle \sum _{i=1}^N{P}_{D\_i}}-P_{D\_array}=0\end{array}.$$

(75)

There are $N+1$ equations and $N+1$ variables in (75), which means the $P_{D\_i}$ solved by (75) is the value of the distributed discharging power with the minimum loss.

When the discharging power of FESU $i$ based on EIP is lower than the lower limit or higher than the upper limit, the discharging power of the unit will be zero or $P_{D\_i\_max}$. Then, the other FESUs will be distributed as the residual discharging power.

4.4. Implementation Scheme of Discharging Control for the FAESS

After the calculation of distributed power value corresponding to the minimum loss based on the EIP, the directive discharging order will be sent to each bottom controller to realize the discharging control coordination of each FESU, as shown in Figure 13. The main steps are as follows:

(1) The directive value of the total charging power $P^{*}{}_{array}$ and the rotating speed of each flywheel ${\omega }_{m\_i}$ are sampled, and systematic parameters are invoked.

(2) The upper limit of discharging power calculated to prevent overcharging and overcurrent of FESU, $P_{OCh\_i\_max}$ and $P_{OCur\_i\_max}$.

(3) The upper discharging limit of each FESU, $P_{D\_i\_max}$, is determined.

(4) The coefficients, $\alpha {\prime }_i$, $\beta {\prime }_i$, and $\gamma {\prime }_i$ are calculated using the quadratic equation of the discharging loss of FESU $P_{loss\_i}(P_{D\_i})$.

(5) The minimum and maximum values of the incremental ratio of the discharging consumption of each FESU ${\lambda }_i$ are calculated using

$$\{\begin{array}{l}{\lambda }_{i\_min}={\beta }_i^{\prime }\\ {\lambda }_{i\_max}==2{\alpha }_i^{\prime }{P}_{D\_i\_max}+{\beta }_i^{\prime }\end{array}(i=1,2,\dots ,N).$$

(76)

(6) The minimum and maximum values of the incremental ratio of public consumption of FAESS are determined:

$$\{\begin{array}{l}{\lambda }_{min}=\min ({\lambda }_{1\_min},{\lambda }_{2\_min},\cdots ,{\lambda }_{N\_min})\\ {\lambda }_{max}=\max ({\lambda }_{1\_max},{\lambda }_{2\_max},\cdots ,{\lambda }_{N\_max})\end{array}.$$

(77)

(7) The incremental ratio of public consumption is set as ${\lambda }^{\left(0\right)}={\lambda }_{max}$.

(8) The charging power $P_{D\_i}^{(0)}$ corresponding to ${\lambda }^{\left(0\right)}$ is calculated using

$$P_{D\_i}^{(0)}=\frac{{\lambda }^{(0)}-{\beta }_i^{\prime }}{2{\alpha }_i^{\prime }}(i=1,2,\dots ,N).$$

(78)

If $P_{D\_i}^{(0)}$ is over the limit, then the limit value is used.

(9) The deviation between sum of discharging power and the total discharging power $P_{D\_array}$ corresponding to ${\lambda }^{\left(0\right)}$ is calculated using

$$\Delta P_D^{(0)}={\displaystyle \sum _{i=1}^N{P}_{D\_i}^{(0)}}-P_{D\_array}.$$

(79)

(10) Verification of whether the power deviation will satisfy the permitted range is carried out using

$$\Delta P_D^{(0)}\le \epsilon .$$

(80)

If (80) is satisfied, then $P_{D\_i}^{(0)}$ is the value of the distributed discharging power with the minimum loss, and $P_i^{\ast }$ is given the value of $P_{D\_i}^{(0)}$ as the directive value of the charging power of FESU.

(11) If (80) is not satisfied, and ${\lambda }^{\left(0\right)}>{\lambda }_{min}$, then

$${\lambda }^{\left(1\right)}>{\lambda }^{\left(0\right)}-\Delta \lambda .$$

(81)

The iterative computations will start again from step (8) until the power deviation satisfies the permitted range or the incremental ratio of public consumption is lower than or equal to ${\lambda }_{min}$.

The discharging coordination control program with EIP is shown in Figure 14.

5. Simulation and Experimental Verification

Examples are analyzed to verify the effectiveness of the charging–discharging control strategy for FAESS based on EIP. The IGBT parameters used in the converters of the standardized FESU are shown in

Table 1. Insulated gate bipolar transistor (IGBT) parameter of the converter.

IGBT On Energy, Eon_test	51 mJ	IGBT Off Energy, Eoff_test	45.5 mJ
Diode on Energy, Erec_test	32.5 mJ	Switching Frequency, fsw	6000 Hz
Test Voltage, Vtest	900 V	Test Current, Itest	225 A
VT0	1.14 V	VD0	1.1925 V
RCE	0.0036 ohm	RD	0.0027 ohm

. The parameters of the permanent magnet synchronous motor and flywheel are shown in

Table 2. Parameters of the Flywheel of the Permanent Magnet Synchronous Motor.

Number of Pole-pairs, p	2	Kc	0.11 ohm/rpm
Stator Resistance, Rs	0.097 Ω	Excitation flux linkage, Ψf	0.1286 Wb
D-axis inductance, Ld	1.435 mH	q-axis inductance, Lq	2.085 mH
Viscous friction coefficient, B	0.0035 Nm/(rad/s)	Flywheel rotating inertia, J	2.063 kg m2

. The effective value of the grid voltage is 270 V, the frequency is 50 Hz, and the DC bus voltage of converter $U_{dc}$ is 500 V.

According to Table 1 and Table 2, the loss constants of the standardized FESU can be calculated as shown in

Table 3. Total Loss Constant of the FESU.

b	5.8733	c	0.004725
d	0.0178	f	4.321 × 10−8
g	0.1455	h	0.3858
k1	2.4829 × 10−5	k2	−3.79 × 10−6
k3	0.094457	l	−0.9822

.

The rated power of the standardized FESU is $P_{rated}=40\mathrm{kW}$; the maximum q-axis current is $i_{qmax}=99\mathrm{A}$; the flywheel rotating speed range is [5000 rpm, 10,000 rpm], the upper speed limit is $n_{max}=10,000\mathrm{rpm}$, and the lower speed limit is $n_{min}=5000\mathrm{rpm}$. The time interval of the power distribution is $\Delta T=1\mathrm{s}$.

5.1. Charging Example Analysis of FAESS

The FAESS that contains three standardized FESUs is used as an example. A simulation analysis is conducted on three strategies: Equal distribution, distribution by chargeable energy, and distribution by the EIP. The initial rotating speeds of each FESU are $n_1=5000\mathrm{rpm}$, $n_1=7000\mathrm{rpm}$, and $n_1=8000\mathrm{rpm}$, respectively. The speed limits are all $n_{max}=10,000\mathrm{rpm}$. The total charging power is $P_{C\_array}=60\mathrm{kW}$.

(1) The results of charging power distributed by equal distribution

With equal distribution, the directive values of the charging power, the actual charging power, the rotating speed of the flywheel rotor, and the q-axis current distribution are as shown in Figure 15a–d.

According to Figure 15a,b, the directive values of charging power of the three units are all equal to −20 kW. The actual charging powers of the three units are also around −20 kW, which corresponds to the charging directive value. According to Figure 15d, the absolute values of the q-axis currents of the three motors are gradually decreasing, because the charging power is constant and the rotating speed of each motor is increasing. Therefore, the torque of each motor is gradually decreasing, and so is the absolute value of the q-axis current. In addition, as shown in Figure 15c, the rotating speeds of the three flywheels increase from 5000, 7000, and 8000 rpm to 7539, 8937, and 9699 rpm, respectively. After calculation, the total kinetic energy increment is $1049.3\mathrm{kJ}$, and the total energy loss is $150.7\mathrm{kJ}$.

(2) The results of power charging by chargeable energy distribution

For the chargeable energy distribution, the directive value of the charging power, the actual charging power, the rotating speed of the flywheel rotor, and the q-axis current distribution are as shown in Figure 16a–d.

As shown in Figure 16a,b, the lower the rotating speed is, the larger the chargeable energy of the FESU is and the greater the absolute value of the assigned directive value of charging power is. The charging power of the three units can follow the directive value. According to Figure 16d, when the q-axis current of the first flywheel motor exceeds its maximum value (99 A), overcurrent occurs. In addition, as shown in Figure 16c, the rotating speeds of the three flywheels increase from 5000, 7000, and 8000 rpm to 8250, 8805, and 9119 rpm, respectively. Therefore, the total kinetic energy increment is $1026.3\mathrm{kJ}$, and the total energy loss is $173.7\mathrm{kJ}$.

(3) The results of power charging distributed by the EIP

With optimal distribution by the EIP, the directive value of the charging power, the actual charging power, the rotating speed of the flywheel rotor, and the q-axis current distribution are as shown in Figure 17a–d.

As shown in Figure 17a, the directive charging value undergoes a larger change when the EIP is used. The absolute value of the directive charging value of the third unit begins to descend and approaches zero, and the absolute values of the directive charging values of the other two start ascending. As shown in Figure 17c, at 18 s, the rotating speed of the third flywheel reaches the upper limit. In terms of the constraints of the EIP, the charging power of the third flywheel is automatically limited. Since the total charging power is fixed, it will only be distributed between the first and second units. Thus, the absolute values of the directive charging values of the two FESUs begin to increase.

According to Figure 17d, the actual charging power of the three units can follow the directive value of the power charging. As seen in Figure 17c,d, there is no over-charging or over-discharging in any of the three FESUs.

In addition, as shown in Figure 17c, the rotating speeds of all three units increase from 5000, 7000, and 8000 rpm to 6797, 9214, and 10,000 rpm, respectively. The total kinetic energy increment is $1052.9\mathrm{kJ}$, and the total energy loss is $147.1\mathrm{kJ}$.

The calculation above shows that when the charging power is distributed by equal distribution, chargeable energy, and EIP, the total energy loss of the FAESS is $150.7\mathrm{kJ}$, $173.7\mathrm{kJ}$ and $147.1\mathrm{kJ}$, respectively. A comparison of these three distribution methods is shown in Figure 18.

According to Figure 18, when the charging power is distributed by EIP, the total energy loss of the FAESS is at the minimum; the loss can be lowered by 15.3% and 2.4% compared with distribution by chargeable energy and equal distribution, respectively. It can be seen from the above simulation analysis results that power charging distribution using the EIP can not only decrease loss, but also avoid overcharging and overcurrent. The simulation results verify the effectiveness of this charging control strategy.

5.2. FAESS Discharging Example

Four strategies—equal distribution, distribution by rotating speed, distribution by residual energy, and distribution by EIP—are simulated. The FAESS that contains three standardized FESUs is again used as an example. The initial rotating speeds of each FESU are $n_1=10,000\mathrm{rpm}$, $n_2=8000\mathrm{rpm}$, and $n_3=7000\mathrm{rpm}$, respectively. The lower limit of the rotating speed is $n_{min}=5000\mathrm{rpm}$ and the total discharging power is $P_{D\_array}=60\mathrm{kW}$.

(1) The results of discharging power distributed by equal distribution

With equal distribution, the directive value of discharging power, the actual discharging power, the rotating speed of the flywheel rotor, and the absolute value of the q-axis current distribution are as shown in Figure 19a–d.

According to Figure 19a,b, the directive values of the discharging powers of the three units are equal to 20 kW. The actual discharging powers of the three units are also around 20 kW, which corresponds to the discharging directive value. According to Figure 19d, the absolute values of q-axis currents of the three motors gradually increase, because the discharging power is constant and the rotating speed of each motor is decreasing. Therefore, the torque of each motor gradually increases, as does the absolute value of the q-axis current. When the q-axis currents of the second and third motors exceed the maximum value of 99 A, overcurrent occurs.

In addition, as shown in Figure 19c, the rotating speeds of the three flywheels are lowered from $10,000\mathrm{rpm}$, $8000\mathrm{rpm}$, and $7000\mathrm{rpm}$ to $7666\mathrm{rpm}$, $4931\mathrm{rpm}$, and $3079\mathrm{rpm}$, respectively. When the speeds of the second and third flywheels are lower than the lower limit, overcharging occurs. The total kinetic energy decrement is $1362.1\mathrm{kJ}$, and the total energy loss is $162.1\mathrm{kJ}$.

(2) The results of discharging power distributed by the rotating speed

When the discharging power is distributed by the rotating speed, the directive value of the discharging power, the actual discharging power, the rotating speed of the flywheel rotor, and the absolute value of the q-axis current distribution are as shown in Figure 20a–d.

According to Figure 20a,b, the higher the rotating speed, the larger the directive value of discharging power. The actual powers of the three units can follow the directive value of discharging power. As seen in Figure 20d, the q-axis currents of three motors are basically the same.

In addition, as shown in Figure 20c, the rotating speeds of the three flywheels are lowered from $10,000\mathrm{rpm}$, $8000\mathrm{rpm}$, and $7000\mathrm{rpm}$ to $7072\mathrm{rpm}$, $5114\mathrm{rpm}$, and $4123\mathrm{rpm}$, respectively. When the speed of the third flywheel is lower than the lower limit of 5000 rpm, over-discharging occurs. It is calculated that the total kinetic energy decrement is $1355.3\mathrm{kJ}$, and the total energy loss is $155.3\mathrm{kJ}$.

(3) The results of discharging power distributed by residual energy

According to Figure 21a,b, the FESU with more residual energy will be distributed with more discharging power. The actual discharging powers of the three units follow the directive value. According to Figure 21d, when the absolute value of the q-axis current of the first motor exceeds the maximum value of 99 A, overcurrent occurs.

In addition, as shown in Figure 21c, the rotating speeds of the three flywheels are lowered from $10,000\mathrm{rpm}$, $8000\mathrm{rpm}$, and $7000\mathrm{rpm}$ to $5960\mathrm{rpm}$, $5530\mathrm{rpm}$, and $5329\mathrm{rpm}$, respectively. The total kinetic energy decrement is $1340.2\mathrm{kJ}$, and the total energy loss is $140.2\mathrm{kJ}$.

(4) The results of discharging power distributed by the EIP

When the discharging power is distributed by the EIP, the directive value of the discharging power, the actual discharging power, the rotating speed of the flywheel rotor, and the absolute value of the q-axis current distribution are as shown in Figure 22a–d.

According to Figure 22, at 18 s, the value of the discharging power distributed by the EIP undergoes a large change, and the directive value of the third unit becomes zero. According to Figure 22c, at 18 s, the rotating speed of the third unit reaches the lowest limit. With the constraint for preventing over-discharging in the distribution strategy based on the EIP, the discharging power of the third unit is automatically limited to zero.

According to Figure 22b, the actual discharging powers of the three FESUs can follow the directive value of the discharging power. As seen in Figure 22c,d, there is no over-discharging and overcurrent in any of the three FESUs.

As shown in Figure 22c, the rotating speeds of the three flywheels are lowered from $10,000\mathrm{rpm}$, $8000\mathrm{rpm}$, and $7000\mathrm{rpm}$ to $6620\mathrm{rpm}$, $5279\mathrm{rpm}$, and $5004\mathrm{rpm}$, respectively. The total kinetic energy decrement is $1315.0\mathrm{kJ}$, and the total energy loss is $138.7\mathrm{kJ}$.

The above simulation analysis shows that when the discharging power is distributed by equal distribution, rotating speed, residual energy, and the EIP, the total energy loss of the FAESS is $162.1\mathrm{kJ}$, $155.3\mathrm{kJ}$, $140.2\mathrm{kJ}$, and $138.7\mathrm{kJ}$, respectively. A comparison of these four distribution methods is shown in Figure 23.

According to Figure 23, when the discharging power is distributed by the EIP, the total discharging energy loss of the system is at the minimum. Compared with the loss of the average distribution, the rotating speed distribution, and the residual energy distribution, the energy loss of can be reduced by 14.4%, 10.7%, and 1.1%, respectively, when the EIP is used. It can be seen from the above simulation analysis that the discharging power distribution based on the EIP can not only decrease the loss, but also avoid over-discharging and overcurrent.

The simulation results verify the effectiveness of the charging–discharging control strategy of FAESS based on the EIP.

6. Conclusions

In this paper, the charging–discharging control strategy for the FAESS is studied in-depth. Firstly, the power loss problem of the FESU during the charging and discharging operation is studied. The results show that the power loss is directly related to power charging–discharging. Based on a relevant theoretical analysis, this paper establishes an objective function of the FAESS with the minimum total power loss as the optimization target. Based on the EIP, it also proposes the constraints of overcharging, over-discharging, and overcurrent. Secondly, a specific implementation scheme for the charging–discharging control strategy for the FAESS based on EIP is proposed. Finally, the charging power distribution strategies (equal distribution, distribution by chargeable energy, distribution by the EIP) are analyzed and compared with examples. In addition, the discharging power distribution strategies (equal distribution, distribution by rotating speed, distribution by residual energy, distribution by the EIP) are analyzed and compared with examples. The results show that the control effect of the charging–discharging control strategy of the FAESS based on the EIP is better than the other control strategies.