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

: The widely used ﬂywheel 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 ﬂywheel array energy storage system (FAESS), which includes the multiple standardized ﬂywheel energy storage unit (FESU), is an e ﬀ ective 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 speciﬁc implementation scheme of the optimal coordinated control strategy. Lastly, the charging–discharging coordinated control strategy is veriﬁed by examples. The results show that the coordinated control strategy can e ﬀ ectively 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 e ﬀ ect than the existing charging–discharging control strategies.


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

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)

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.

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)

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): P sw_G = 6 × (P sw_IGBT + P sw_D ) 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 o f f _test represents the switching-off loss of the IGBT; E rec_test represents where _ s w G P represents the switching loss 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 2 0 where I G B T P represents the conduction loss of the IGBT in the one-phase leg of the grid-side converter; 0 T V represents the IGBT threshold voltage; m I represents the output phase current in the grid-side converter; m represents the modulation ratio; ϕ represents the power factor angle; and C E R represents the IGBT conduction resistance. By calculating the diode conduction loss, we can acquire 2 0 where D P represents the conduction loss of the diode in the one-phase leg of the grid-side converter; 0 D V represents the diode threshold voltage; m I represents the output phase current in the grid-side converter; m represents the modulation ratio; ϕ represents the power factor angle; and D R represents the diode conduction resistance. The total conduction loss of grid-side converter c o n d P is The absolute value of mcosϕ is less than 1, and 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 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; ϕ represents the power factor angle; and R CE represents the IGBT conduction resistance. By calculating the diode conduction loss, we can acquire 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; ϕ represents the power factor angle; and R D represents the diode conduction resistance. The total conduction loss of grid-side converter P cond is The absolute value of m cos ϕ is less than 1, and Therefore, the total conduction loss can be simplified to Energies 2019, 12, 2844

of 27
According to (1) and (7), the total conduction loss in the grid-side P G_tot is Since the relation between the amplitude of output phase current I m and the active power directive value P * is the total loss in the grid-side converter P G_tot is 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 i q of the q axis current i q . Similarly, the total loss in the motor-side converter P M_tot is

Motor Loss of FESU
As shown in Figures 2 and 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 where 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. where s R 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 c R . The equivalent steady state circuit of the permanent magnet synchronous motor that includes the iron loss resistor c R is shown in Figure 4.  From Figure 4, we can acquire (13): According to (13), the iron loss of motor Fe P is From 0 ( ) we can acquire Bringing (16) into (14), we can acquire From Figure 4, we can acquire (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 ; ω e represents the electric angular speed of the motor; ψ 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 From we can acquire Bringing (16) into (14), we can acquire where ω m is the mechanical angular speed of the flywheel rotor and The friction loss of the flywheel P f r can be represented by where B represents the viscous friction coefficient and ω m represents the mechanical angular speed of the flywheel rotor.  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 where  During charging in the grid, the total loss lossin P of FESU is  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 where Energies 2019, 12, x FOR PEER REVIEW 8 of 28 electromagnetic power of the flywheel motor; and F P is the actual output power of the flywheel.
During charging of the grid, the total loss lossout P of FESU is 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.  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.

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: 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

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

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 (2) Inequity constraints of overcharging prevention The current rotating speed of FESU i is set as ω m_i , and the kinetic energy is set as E i . After operating cycle ∆T, the rotating speed is changed to ω 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 ω m_i_next should not exceed the maximum permitted rotating speed ω max , as shown in Figure 7.
is the current input power of i of the flywheel rotor. Equation (25) can be represented by According to (26), It can be seen from (27) and (25) that if (28) is established, (25) is also established: 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).
where P F_i is the current input power of i of the flywheel rotor. Equation (25) can be represented by According to (26), It can be seen from (27) and (25) that if (28) is established, (25) is also established: From (28), (29) is acquired: Then, According to (32), 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 ∆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 ∆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: According to (34), and then, where  According to (24), (33), and (36), the upper limit of charging power should be the minimum of these three values, which is Comprehensively, when the FAESS is charging, the charging power of each unit should meet the following inequality constraints:

Research of the Charging Control Method for FAESS
The charging power distribution of FAESS can be presented with the equality constraint of Let the objective function: where λ represents the LaGrange multiplier.
The necessary condition for LaGrange function L to take the extreme value is Comprehensively, when the FAESS is charging, the charging power of each unit should meet the following inequality constraints:

Research of the Charging Control Method for FAESS
The charging power distribution of FAESS can be presented with the equality constraint of Let the objective function: be the minimum, where p loss_i (P C_i ) = α i P 2 C_i + β i P C_i + γ 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): where λ represents the LaGrange multiplier. The necessary condition for LaGrange function L to take the extreme value is 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, where λ i represents the incremental rate of consumption of each FESU. Equation (42) can be rewritten as 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.

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 charging loss loss_i _ ( ) C i P P of each unit is only the function that relates to its charging power where i λ represents the incremental rate of consumption of each FESU.

Equation (42) can be rewritten as
There are 1 N + equations and 1 N + variables in (44), which means the value of _ 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 _ _ max C i P . Then, the residual charging power will be distributed to other FESUs.

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 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 ω 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, α i , β i , and γ 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, λ i , are determined: (6) The minimum and maximum values of the incremental ratio of public consumption of the FAESS are determined: The incremental ratio of public consumption is set as C_i corresponding to λ (0) is calculated using If P (0) C_i 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 λ (0) is calculated: (10) Verification of whether the power deviation will satisfy the permitted range is carried out using If (49) is satisfied, then P C_i is the value of distributed charging power with minimum loss, then P * i is given a negative value of P C_i as the directive value of the charging power of the FESU. (11) If it is not satisfied, and λ (0) > λ min , then 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 λ min .
The charging coordination control program based on EIP is shown in Figure 10.

Discharging Objective Function for the FAESS
When the FAESS is discharging, the loss _ loss i P of FESU i is related to its active discharging In the case that the total discharging power is fixed, the distributed power of FESU i _ D i P is different, and so is its loss loss_i P and the total loss of the system. With a fixed total discharging power, in order to minimize the total loss loss_tot P , the objective function is defined as follows: where

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: where

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

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 (2) Inequity constraint of preventing over-discharging The current rotating speed of FESU i is set as ω m_i and the kinetic energy is set as E i . After operating cycle ∆T, the rotating speed is changed to ω 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 ω m_i_next should not be lower than the minimum permitted rotating speed ω min , as shown in Figure 11.

Discharging Constraints of the FAESS
(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 _ D a r r a y P , which is (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 (2) Inequity constraint of preventing over-discharging The current rotating speed of FESU i is set as m_i ω and the kinetic energy is set as i E . After operating cycle T Δ , the rotating speed is changed to _ _ m i next ω and the kinetic energy is changed to i_next E . In order to prevent over-discharging, the rotating speed of the next moment _ _ m i next ω should not be lower than the minimum permitted rotating speed min ω , as shown in Figure 11. The kinetic energy of FESU in the next moment i _ n e x t E should be higher than its minimum limit min E , and Equation (55) should be guaranteed: In (55), _ F i P represents the current output power of flywheel rotor i : According to (55) and (56), which is Figure 11. Prevention of over-discharging.
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: In (55), P F_i represents the current output power of flywheel rotor i: According to (55) and (56), which is 1 2 and (60) can be acquired from (59): where 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 ∆T, the q axis current in the next moment will be i q_i_next . Since it is within operating cycle ∆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 i q_i_next 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.
and (60) can be acquired from (59): 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 q_i i is set. After operating cycle T Δ , the q axis current in the next moment will be _ _ i q i next . Since it is within operating cycle 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 i q _ i _ n e x t of the q axis current in the next moment cannot exceed the maximum permitted q axis current i q _ m a x , as shown in Figure 12. During discharging, there is and and Figure 12. Prevention of overcurrent during discharging.
From (61) and (62) we can acquire From (66), we can acquire 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 ).
Generally, when the FAESS is discharged, the discharging power of each unit should meet the following inequity constraints: N). (69)

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: and the objective function is the minimum. P loss_i (P D_i ) = a i P 2 D_i + β i P D_i + γ 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).
where λ is the LaGrange multiplier. The necessary conditions for LaGrange function L at the extremum are 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, where λ i represents the incremental rate of consumption of each FESU. Equation (73) can be rewritten as 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.

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 ω 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, α i , β i , and γ 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 λ i are calculated using (6) The minimum and maximum values of the incremental ratio of public consumption of FAESS are determined: The incremental ratio of public consumption is set as λ (0) = λ max . (8) The charging power P D_i corresponding to λ (0) is calculated using If P (0) D_i 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 λ (0) is calculated using (10) Verification of whether the power deviation will satisfy the permitted range is carried out using If (80) is satisfied, then P D_i is the value of the distributed discharging power with the minimum loss, and P * i is given the value of P 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 λ min .
is over the limit, then the limit value is used. (9) The deviation between sum of discharging power and the total discharging power 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 m in λ .  The discharging coordination control program with EIP is shown in Figure 14. The discharging coordination control program with EIP is shown in Figure 14. Calculate the power deviation ( )

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. The parameters of the permanent magnet synchronous motor and flywheel are shown in Table 2. The effective value of the grid voltage is 270 V, the frequency is 50 Hz, and the DC bus voltage of converter d c U is 500 V.

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. The parameters of the permanent magnet synchronous motor and flywheel are shown in Table 2. 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 Tables 1 and 2, the loss constants of the standardized FESU can be calculated as shown in Table 3. The rated power of the standardized FESU is P rated = 40 kW; the maximum q-axis current is i qmax = 99 A; the flywheel rotating speed range is [5000 rpm, 10,000 rpm], the upper speed limit is n max = 10, 000 rpm, and the lower speed limit is n min = 5000 rpm. The time interval of the power distribution is ∆T = 1 s.

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 rpm, n 1 = 7000 rpm, and n 1 = 8000 rpm, respectively. The speed limits are all n max = 10, 000 rpm. The total charging power is P C_array = 60 kW. According to Tables 1 and 2, the loss constants of the standardized FESU can be calculated as shown in Table 3. The rated power of the standardized FESU is

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 1    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 kJ, and the total energy loss is 150.7 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.
Energies 2019, 12, x FOR PEER REVIEW 22 of 28 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 1 0 4 9 .3 k J , and the total energy loss is 150.7 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 1 0 2 6 .3 k J , and the total energy loss is 173.7 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 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 kJ, and the total energy loss is 173.7 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 1 0 5 2 .9 k J , and the total energy loss is 1 4 7 .1 k J .
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 kJ , 173.7 kJ and 1 4 7 .1 k J , respectively. A comparison of these three distribution methods is shown in Figure 18.  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 kJ, and the total energy loss is 147.1 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 kJ, 173.7 kJ and 147.1 kJ, respectively. A comparison of these three distribution methods is shown in Figure 18. 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 1 0 5 2 .9 k J , and the total energy loss is 1 4 7 .1 k J .
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 kJ , 173.7 kJ and 1 4 7 .1 k J , 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.

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 rpm, n 2 = 8000 rpm, and n 3 = 7000 rpm, respectively. The lower limit of the rotating speed is n min = 5000 rpm and the total discharging power is P D_array = 60 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 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.

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 = 10 000 rpm (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 rpm ， , 8000 rpm , and 7000 rpm to 7666 rpm , 4931 rpm , and 3079 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 1 3 6 2 .1 k J , and the total energy loss is 1 6 2 .1 k J .
(2) The results of discharging power distributed by the rotating speed 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 rpm, 8000 rpm, and 7000 rpm to 7666 rpm, 4931 rpm, and 3079 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 kJ, and the total energy loss is 162.1 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- 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 rpm ， , 8000 rpm , and 7000 rpm to 7072 rpm , 5114 rpm , and 4123 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 1 3 5 5 .3 k J , and the total energy loss is 155.3 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. 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 rpm, 8000 rpm, and 7000 rpm to 7072 rpm, 5114 rpm, and 4123 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 kJ, and the total energy loss is 155.3 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. 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 rpm ， , 8000 rpm , and 7000 rpm to 7072 rpm , 5114 rpm , and 4123 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 1 3 5 5 .3 k J , and the total energy loss is 155.3 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 rpm , 8000 rpm , and 7000 rpm to 5960 rpm , 5530 rpm , and 5329 rpm , respectively. The total kinetic energy decrement is 1 3 4 0 .2 k J , and the total energy loss is 140.2 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 rpm ， , 8000 rpm , and 7000 rpm to 6620 rpm , 5279 rpm , and 5004 rpm , respectively. The total kinetic energy decrement is 1 3 1 5 .0 k J , and the total energy loss is 138.7 kJ . In addition, as shown in Figure 21c, the rotating speeds of the three flywheels are lowered from 10, 000 rpm, 8000 rpm, and 7000 rpm to 5960 rpm, 5530 rpm, and 5329 rpm, respectively. The total kinetic energy decrement is 1340.2 kJ, and the total energy loss is 140.2 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. In addition, as shown in Figure 21c, the rotating speeds of the three flywheels are lowered from 10, 000 rpm , 8000 rpm , and 7000 rpm to 5960 rpm , 5530 rpm , and 5329 rpm , respectively. The total kinetic energy decrement is 1 3 4 0 .2 k J , and the total energy loss is 140.2 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 rpm ， , 8000 rpm , and 7000 rpm to 6620 rpm , 5279 rpm , and 5004 rpm , respectively. The total kinetic energy decrement is 1 3 1 5 .0 k J , and the total energy loss is 138.7 kJ . 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 rpm, 8000 rpm, and 7000 rpm to 6620 rpm, 5279 rpm, and 5004 rpm, respectively. The total kinetic energy decrement is 1315.0 kJ, and the total energy loss is 138.7 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 kJ, 155.3 kJ, 140.2 kJ, and 138.7 kJ, respectively. A comparison of these four distribution methods is shown in Figure 23. 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 1 6 2 .1 k J , 155.3 kJ , 140.2 kJ , and 138.7 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.

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.