Research on Improved VSG Control Algorithm Based on Capacity-Limited Energy Storage System

Abstract: A large scale of renewable energy employing grid connected electronic inverters fail to contribute inertia or damping to power systems, and, therefore, may bring negative effects to the stability of power system. As a solution, an advanced Virtual Synchronous Generator (VSG) control technology based on Hamilton approach is introduced in this paper firstly to support the frequency and enhance the suitability and robustness of the system. The charge and discharge process of power storage devices forms the virtual inertia and damping of VSG, and, therefore, limits on storage capacity may change the coefficients of VSG. To provide a method in keeping system output in an acceptable level with the capacity restriction in a transient period, an energy control algorithm is designed for VSG adaptive control. Finally, simulations are conducted in DIgSILENT to demonstrate the correctness of the algorithm. The demonstration shows: (1) the proposed control model aims at better system robustness and stability; and (2) the model performs in the environment closer to practical engineering by fitting the operation state of storage system.


Introduction
In recent years, penetration of distributed generation (DG) in power systems is increasing rapidly, accompanied with the reduction of total system inertia.To play their potential fully and reduce their negative effects at the same time, specific functionalities are designed for the converter interfacing DG and power grid [1][2][3][4].Virtual synchronous generator (VSG) technology [5,6], which imitates the output characteristic of synchronous generator (SG) to enhance the stability of power grid, is one of the effective methodologies and attracting more and more attention.
In the project of Virtual Synchronous Machines (VSYNC) in Poland, an idea named VSG was put forward the first time, aiming at making the converter imitate the mechanical equation for a more stable grid frequency [7].A further study in excitation control was introduced later, completing and developing the VSG system [8].However, since the immaturity of the early VSG technology, the dynamic response of the systems are commonly slow and accompanied with slight oscillations.
Small-signal modeling is one of the mainstream methods in linear system analysis and it is widely used in VSG for system stabilization.Comparisons are illustrated in [9] between droop control and VSG control.Superiority of VSG has been proved by comparing dynamic responses via small-signal modeling and state equations.In [10], every specific elements of the inverter control modeling are listed and analyzed.The corresponding small-signal model is then established for obtaining controller sensitivity and revealing the effects of different parameters towards VSG system.
One of the merits of VSG compared to the SG is that the former one enjoys flexible parameters selection.Analysis in [11] shows that the virtual inertia enjoys self-adaption depending on disturbance.In [12,13], an adaptive alternating inertial control algorithm is introduced.They verify the validity that the VSG can work at the state of negative inertia, being more flexible and adjustable than the SG.With a further study in [14,15], virtual inertia varies to meet the demands according to advanced control algorithm rather than a fixed coefficient in [12,13].Based on the studies above, VSG is commonly regarded as the second-order system for simplification by neglecting parts of SG, resulting in a less precise modeling.In this paper, a more accurate model is established and performs better in mimicking SG characteristics.The generalized Hamilton function method enjoying its superiority in nonlinear system research for the merits in depicting the dynamic energy flow [16][17][18].The technology is attracting more and more attention gradually.In [19], a fifth-order model of the SG is established with an energy-based Lyapunov function.In [20,21], a further study is conducted that the modeling of generator integrated AVR (Automatic Voltage Regulator) and PSS (Power System Stabilization) shows the superiority of the Hamilton approach in improving system robustness and stabilization.Considering the superiority, an advanced control algorithm based on the Hamilton approach is applied to the VSG for better suitability and robustness.Energy-based control laws are conducted via applying the Hamiltonian approach to valving control and excitation system in this paper.
Electronic inverters fail to contribute inertia or damping to power grid, leading to a necessity that the energy storage devices should be equipped to compensate the power vacancy [22].In current study, DG and storage devices are commonly simplified as DC power sources.In [23,24], conclusions can be drawn that it is the charge and discharge progress of the storage devices that form the virtual inertia.In [25], an optimal control strategy of the storage systems by setting the VSG's output active power depending on the grid frequency deviation is put forward.In [26], storage devices act as power buffers to compensate power shortage in short-term impact for transient stability.Existing studies focus on the significance of storage devices and binary control strategies between VSG and storage devices, neglecting the connections between the storage devices and the VSG coefficients.Furthermore, storage capacity is generally considered as unlimited in theoretical analysis, being considerably different from the situation in practical engineering.Given the insufficiencies, corresponding analyses and solutions are proposed.
Contributions of this paper are listed as follows: 1.
A fourth-order model of VSG, being more detailed and accurate in mimicking SG, is established.
Based on the second-order model, the steam-valving system and excitation system are added, taking the coordination between the two systems into consideration while avoiding the complicated transition period.

2.
A nonlinear control method, Hamilton approach is applied to VSG for better system suitability and stability.Control laws are deduced and the system enhancements are illustrated.

3.
Small-signal model is established for the obtaining boundaries of VSG coefficients.

4.
Relations between the capacity limits of storage devices and the VSG coefficients, such as the virtual inertia and damping, are discussed.An adaptive VSG control is designed via the energy control algorithm for accelerating the recovery from the transient period.
This paper is organized as follows: In Section 2, VSG modeling aspects are illustrated.Energy controller unit is then conducted in Section 3. Simulations are given in Section 4. Conclusions are drawn in Section 5.

VSG Control System Overview
Detailed classification methods in VSG modeling are illustrated in [27] based on the order of SG modeling.In this paper, a fourth-order model contains the steam-valving system and steam turbine generator excitation system is built, performing better in revealing the coordination between the valve Energies 2018, 11, 677 3 of 17 and the exciter.Therefore, a control algorithm based on energy function is deduced by making use of Hamilton control system with virtual synchronous generator excitation (VSG) system and turbine valve control, in which the physical significance is clearer than before.Figure 1  valve and the exciter.Therefore, a control algorithm based on energy function is deduced by making use of Hamilton control system with virtual synchronous generator excitation (VSG) system and turbine valve control, in which the physical significance is clearer than before.Figure 1 depicts the overall block diagram of VSG control strategy.As indicated in Figure 1, Udc, which is assumed that the dc-link of the VSG is connected to an energy storage unit or to a source with sufficient available buffer capacity, acts as the virtual prime mover.The inductance of the filter inductor is Lf and filter capacitor is Cf.The output voltage for the converter is denoted by e, representing the transient voltage of synchronous generator while the current reference, tracked by the controller, is given by if.Meanwhile, u0 is the voltage at the filter capacitors and i0 is the current flowing into the grid or load equivalent, which are equivalent to the terminal voltage and current of stator winding of a synchronous generator (SG), making the dynamic parameters at grid connection point, such as u0, i0, f and Pe, mimic the characteristic of a SG.

VSG Modeling Based on Dissipative Hamilton Theory
VSG shares the same mathematical model as the SG, and its mechanical equations are given as Equation ( 1) where D is damping, M is virtual inertia, Pe is the system output power, Pm the mechanical power, ω the angular velocity, ω0 the angular velocity references.
The mechanical power Pm considering the steam-valving is given as Equation (2).
where Ts represents the time constant, Pm0 is the mechanical power reference and the valving control signal.
The output power Pe of the VSG is given as Equation (3).
where δ is the power angle of VSG, Eq' the transient electrodynamic force (E.M.F) on q-axis, Us the voltage of infinite bus-system and x'dΣ the total stator-winding reactance of d-axis.
The excitation control equation is given as Equation ( 4).As indicated in Figure 1, U dc , which is assumed that the dc-link of the VSG is connected to an energy storage unit or to a source with sufficient available buffer capacity, acts as the virtual prime mover.The inductance of the filter inductor is L f and filter capacitor is C f .The output voltage for the converter is denoted by e, representing the transient voltage of synchronous generator while the current reference, tracked by the controller, is given by i f .Meanwhile, u 0 is the voltage at the filter capacitors and i 0 is the current flowing into the grid or load equivalent, which are equivalent to the terminal voltage and current of stator winding of a synchronous generator (SG), making the dynamic parameters at grid connection point, such as u 0 , i 0 , f and P e , mimic the characteristic of a SG.

VSG Modeling Based on Dissipative Hamilton Theory
VSG shares the same mathematical model as the SG, and its mechanical equations are given as Equation (1) .
where D is damping, M is virtual inertia, P e is the system output power, P m the mechanical power, ω the angular velocity, ω 0 the angular velocity references.The mechanical power P m considering the steam-valving is given as Equation (2). .
where T s represents the time constant, P m0 is the mechanical power reference and µ g the valving control signal.The output power P e of the VSG is given as Equation (3).
where δ is the power angle of VSG, E q ' the transient electrodynamic force (E.M.F) on q-axis, U s the voltage of infinite bus-system and x' dΣ the total stator-winding reactance of d-axis.
The excitation control equation is given as Equation (4).
where T d0 is the time constant in excitation winding, x d the self-inductance resistance, x dΣ the total stator-winding reactance of d axis and U f the excitation voltage.i d in Equation ( 4) are noted as Equation ( 5).
According to Equations ( 1)-( 5), system state equations are given in Equation ( 6): where x 1 = δ, x 2 = ω − ω 0 , x 3 = E q ', x 4 = P m and Equation ( 6) can be transformed into Equation ( 7) In Equation (7), p 1 and p 2 are parameters with perturbation, u 1 and u 2 are system input and w 1 , w 2 and w 3 disturbances of the system.After obtaining the equilibrium points, Equation (7) can be transformed to Equation (8).Details are shown in Equation (9). . where Equation ( 8) is a dissipative Hamilton approach where g 1 (x), g 2 (x) and v are all column vectors of m dimension and H(x) is regarded as the Hamilton function.Matrix J is antisymmetric and matrix R is positive semidefinite.
Energies 2018, 11, 677 5 of 17 where v 1 and v 2 are the inputs after the introduction of state feedback, whose relations with u 1 and u 2 , the original input of the system, can be given by Equation (10).
In Equation (10), U f is the voltage of the excitation system, defined by the equation U f = U f0 + u f , in which u f is the regulation of the excitation voltage.
Equation ( 9) fits in the form of the dissipative Hamilton system, which can be regarded as one of the Hamilton approaches of VSG.Then, Hamilton function can be deduced after solving ∇H(x) in Equation ( 9) out and noted as Equation ( 12) [28,29].
A penalty function z, an interference suppression level γ and a weight matrix r(X) are firstly defined and given in Equation (13).
If it satisfies the condition in Equation ( 14) in any circumstances, Then, the control law of VSG based on Hamilton approach is given in Equation ( 15): where θ is an intermediate variable.
Details are shown in Equation ( 15) when parameters are substituted to Equation (16).

Improved VSG Control Strategy Considering Storage Capacity Limits
In existing studies, energy storage capacity is generally assumed as sufficient in suppressing oscillations.However, in practical engineering, storage capacity is commonly limited, resulting in fluctuations in system outputs.To keep fluctuations of output active power P e stay within the limits (±20% variation from the steady state value of P e ), a self-adaptive control strategy is designed and coefficients of VSG will change to meet the output demands and survive the transient period.
Connections between the VSG parameters and charge and discharge power/energy will be discussed in the first step and the solutions will be illustrated in the following parts.

Linearized Small-Signal Model
To simplify the analysis, it is assumed that the renewable energy sources have a stable output power.
The output active power of VSG around the steady operating points of the states can be described as P e = P m = P m0 .The small-signal state-space equations of the VSG can be deduced from Equations ( 8) and ( 9), presented as Equation (17).Then, according to Equation ( 16), small-signal model of H∞ control algorithms around the steady operating points of the state are given in Equation (18).
where E' q0 and δ 0 are transient E.M.F and power angle respectively at steady-state point.They are presented in Equation ( 20) deduced from Equation (19).
Regulation of the excitation voltage at steady-state point u f0 is represented as Equation ( 21) based on Equation (16).
With the Laplace transform of Equations ( 17) and ( 18), the transfer function of small-signal modeling is given in Equation (23), revealing the relationship between the small deviation of output active power ∆P e and the small deviation of the system frequency ∆ω g .
Energies 2018, 11, 677 7 of 17 G 1 (s) and G 2 (s) in Equation ( 23) are elaborated in Equation (24).Related elements and equations of the small-signal state-space equation of Equation (23) are presented in the Appendix A.

Boundaries of VSG Parameters
The characteristic equation, given as Equation (25), can be drawn from the denominator of Equation (23).
Five roots of Equation ( 25) are given in Equation ( 26), defined as x 1-5 , worked out via Cardan and Tartaglia formulae.Elements (a 1 , b 1 ), (a 2 , b 2 ) and k are presented in Appendix A, Appendix B and Table A2.
As shown in Equation ( 26), Appendix B and Table A2, x 1-4 are connected to M and D.x 5 is only connected with M and located far away from x 1-4 .In Figure 2a, M = 7.6 while D is variable.When D < 0, the real parts of x 1 , x 2 are positive, resulting in an unstable system.When D increases, x 1 is approaching approximately to zero while x 4 is moving away from zero.When D > 18, x 3 and x 4 become negative real numbers and system stabilizes gradually.
In Figure 2b, D = 3 while M is variable.When M increases, x 3 , x 4 and x 5 are approaching approximately to zero.When M increases, the real parts of x 3 , x 4 and x 5 become negative and the damping ratio of the system decreases.When M ≥ 13, the imaginary parts of x 1 , x 2 , x 3 and x 4 are not 0 anymore.The system is destabilizing for turning to the mode of oscillating decay.As shown in Equation ( 26), Appendix B and Table A2, x1-4 are connected to M and D.x5 is only connected with M and located far away from x1-4.In Figure 2a, M = 7.6 while D is variable.When D < 0, the real parts of x1, x2 are positive, resulting in an unstable system.When D increases, x1 is approaching approximately to zero while x4 is moving away from zero.When D > 18, x3 and x4 become negative real numbers and system stabilizes gradually.
In Figure 2b, D = 3 while M is variable.When M increases, x3, x4 and x5 are approaching approximately to zero.When M increases, the real parts of x3, x4 and x5 become negative and the damping ratio of the system decreases.When M ≥ 13, the imaginary parts of x1, x2, x3 and x4 are not 0 anymore.The system is destabilizing for turning to the mode of oscillating decay.Based on Figure 2, VSG fails to contribute enough inertia to the grid if M is too small while the overshoots occur if M is too big.Proper D should be selected for achieving the overdamping condition without slowing down the response time too much.As a deduction, boundaries of D and M are given in Equation (27).0< 10 It is not very easy to obtain the accurate boundaries of virtual inertia and virtual damping, and in Equation ( 27), margins exist for ensuring the system to keep stability when coefficients meet the certain formulations.

Design of the Energy Controller with the Restriction of Storage Capacity
As mentioned at the beginning of Section 3, objective formulations are given as Equation ( 28).Based on Figure 2, VSG fails to contribute enough inertia to the grid if M is too small while the overshoots occur if M is too big.Proper D should be selected for achieving the overdamping condition without slowing down the response time too much.As a deduction, boundaries of D and M are given in Equation (27).
Energies 2018, 11, 677 8 of 17 It is not very easy to obtain the accurate boundaries of virtual inertia and virtual damping, and in Equation ( 27), margins exist for ensuring the system to keep stability when coefficients meet the certain formulations.

Design of the Energy Controller with the Restriction of Storage Capacity
As mentioned at the beginning of Section 3, objective formulations are given as Equation (28).
where P k.min = 0.8 p.u., P k.max = 1.2 p.u. (±20% variation from the steady state P e , 1.0 p.u.) and ∆P e is the output active power.E k.min is the minimum, E k.max the maximum and ∆E the output energy.E k.min and E k.max will not be assigned and the reasons are illustrated in the following analysis.

Relations between Storage Capacity Limits and VSG Coefficients
The step response curve of Equation ( 23) is figured in the condition of underdamping and overdamping, where S 1 , S 2 , . . ., S 5 and S 1 ' are the areas enclosed by the curves and axes, representing the energy in Figure 3.In Figure 3a, S1 is the largest among all the areas.Either of the response curves reaches their first peak, ∆ e reaches its maximum; the first time response curves reach zero, the energy ∆ (S1 and S1') meet the maximum.Through the analysis, the time power and energy of the storage devices reach their limits can be calculated.Meanwhile, the minimum capacity of the storage devices is convenient to access.
Transfer function in Equation ( 23) can be given in the form containing poles and zeros as Equation (29).
[ ] For the case of a perturbing with ωg, an inverse transformation of Laplace combined with residue analysis is applied to Equation ( 29) to obtain Equation (30).Related elements of ki (i = 1, 2, …, 6) are presented in the Appendix C.
The smallest root of the positive roots in Equation ( 31) is given as t1, the time when ∆ e reaches its maximum.( ) The maximum of ∆ can be calculated in Equation (32).In Figure 3a, S 1 is the largest among all the areas.Either of the response curves reaches their first peak, ∆P e reaches its maximum; the first time response curves reach zero, the energy ∆E (S 1 and S 1 ') meet the maximum.Through the analysis, the time power and energy of the storage devices reach their limits can be calculated.Meanwhile, the minimum capacity of the storage devices is convenient to access.
Transfer function in Equation ( 23) can be given in the form containing poles and zeros as Equation (29).
For the case of a perturbing with ω g , an inverse transformation of Laplace combined with residue analysis is applied to Equation ( 29) to obtain Equation (30).Related elements of k i (i = 1, 2, . . ., 6) are presented in the Appendix C.
After taking derivatives of Equation (30), The smallest root of the positive roots in Equation ( 31) is given as t 1 , the time when ∆P e reaches its maximum.
∆P emax = ∆P e (t 1 ), The maximum of ∆E can be calculated in Equation (32).
where t 2 represents the time when ∆P e becomes 0.00 the first time in underdamping condition, or the time when ∆P e reaches ±0.05 in overdamping condition (∆P e will infinitely approach zero in overdamping condition and the precision of energy deviation is set as ±0.05 for approximate calculation).It is noticeable that t 2 is generally regarded as a very short period when faults occur and, therefore, ∆E is small and formulations in (b) of Equation ( 28) can be generally satisfied, explaining the reasons that E k.min and E k.max without assignment.
Combined with the Appendixs B and C, both ∆P emax and ∆E max are related to M and D and therefore, detailed studies on their relationship are conducted.In Figure 4, the curves show the influences that M/D exert on ∆E max and ∆P emax while, in Figure 5, the curves show the effect that ∆E max /∆P emax exert on M and D.
Energies 2018, 11, x FOR PEER REVIEW 10 of 18 Combined with the Appendix B and C, both ∆ emax and ∆ max are related to M and D and therefore, detailed studies on their relationship are conducted.In Figure 4, the curves show the influences that M/D exert on ∆ max and ∆ emax while, in Figure 5, the curves show the effect that ∆ max/∆ emax exert on M and D. In Figure 4a, relation curves of ∆ max and M, ∆ emax and M are fitted in MATLAB (MATLAB 2012a, MathWorks, USA) when D is constant.When M increases, it is a monotone increasing relationship between ∆ max (or ∆ emax) and M approximately.The bigger the virtual inertia becomes, the larger capacity is needed.In Figure 4b, relation curves of ∆ max and D, ∆ emax and D are fitted when M is constant.When D increases, ∆ max and ∆ emax drop rapidly in a short period, then turn to a mild decrease.Influences that D exert on ∆ max and ∆ emax is smaller than M on ∆ max and ∆ emax.
In Figure 5, it can be concluded that the larger ∆ max and ∆ emax become, the wider the M and D cover.Meanwhile, in any working conditions, relationship between M and D is monotonically when either ∆ max or ∆ emax is constant.Meanwhile, when ∆ max is constant, the variation of M is smaller than those when ∆ emax is constant.In other words, if ∆ e meets the demands in (a) of Equation ( 28), ∆ follows in (b) of Equation ( 28) and the illustrations in defining the Ek.min and Ek.max are proved.In Figure 4a, relation curves of ∆E max and M, ∆P emax and M are fitted in MATLAB (MATLAB 2012a, MathWorks, USA) when D is constant.When M increases, it is a monotone increasing relationship between ∆E max (or ∆P emax ) and M approximately.The bigger the virtual inertia becomes, the larger capacity is needed.In Figure 4b, relation curves of ∆E max and D, ∆P emax and D are fitted when M is constant.When D increases, ∆E max and ∆P emax drop rapidly in a short period, then turn to a mild decrease.Influences that D exert on ∆E max and ∆P emax is smaller than M on ∆E max and ∆P emax .
In Figure 5, it can be concluded that the larger ∆E max and ∆P emax become, the wider the M and D cover.Meanwhile, in any working conditions, relationship between M and D is monotonically when either ∆E max or ∆P emax is constant.Meanwhile, when ∆E max is constant, the variation of M is smaller than those when ∆P emax is constant.In other words, if ∆P e meets the demands in (a) of Equation ( 28), ∆E follows in (b) of Equation ( 28) and the illustrations in defining the E k.min and E k.max are proved.a mild decrease.Influences that D exert on ∆ max and ∆ emax is smaller than M on ∆ max and ∆ emax.
In Figure 5, it can be concluded that the larger ∆ max and ∆ emax become, the wider the M and D cover.Meanwhile, in any working conditions, relationship between M and D is monotonically when either ∆ max or ∆ emax is constant.Meanwhile, when ∆ max is constant, the variation of M is smaller than those when ∆ emax is constant.In other words, if ∆ e meets the demands in (a) of Equation ( 28), ∆ follows in (b) of Equation ( 28) and the illustrations in defining the Ek.min and Ek.max are proved.

Adaptive Control via Energy-Based Algorithms
The energy-based algorithm should work within the restrictions in Equations ( 27) and ( 28).The diagram with an additional energy-based control algorithm is shown in Figure 6.

Adaptive Control via Energy-Based Algorithms
The energy-based algorithm should work within the restrictions in Equations ( 27) and ( 28).The diagram with an additional energy-based control algorithm is shown in Figure 6.In Figure 6, ME is the virtual inertia when D = DE = 20 and MP the virtual inertia when D = DP = 40.
According to Equation ( 27), a detailed algorithm is given in Figure 7.It is noticeable that the value of D is fixed at 20 or 40, the minimum and maximum of its boundaries.The simplification is reasonable for the influences that D exert on ∆ max and ∆ emax is smaller than that M on ∆ max and ∆ emax, proved in Section 3.3.1.When (ME, DE) and (MP, DP) are calculated, the smaller one of ME and MP will be selected and, therefore, ∆ max and ∆ emax, energy and power needed in suppressing oscillations, is smaller.The effeteness of the algorithm will be proved in Section 4.3.

Simulations
Simulations based on VSG of single-machine infinite system were performed on the model in DIgSILENT (DIgSILENT/Power Factory 15.1.7,Germany) of Figure 8  In Figure 6, M E is the virtual inertia when D = D E = 20 and M P the virtual inertia when D = D P = 40.
According to Equation ( 27), a detailed algorithm is given in Figure 7.In Figure 6, ME is the virtual inertia when D = DE = 20 and MP the virtual inertia when D = DP = 40.
According to Equation ( 27), a detailed algorithm is given in Figure 7.It is noticeable that the value of D is fixed at 20 or 40, the minimum and maximum of its boundaries.The simplification is reasonable for the influences that D exert on ∆ max and ∆ emax is smaller than that M on ∆ max and ∆ emax, proved in Section 3.3.1.When (ME, DE) and (MP, DP) are calculated, the smaller one of ME and MP will be selected and, therefore, ∆ max and ∆ emax, energy and power needed in suppressing oscillations, is smaller.The effeteness of the algorithm will be proved in Section 4.3.It is noticeable that the value of D is fixed at 20 or 40, the minimum and maximum of its boundaries.The simplification is reasonable for the influences that D exert on ∆E max and ∆P emax is smaller than that M on ∆E max and ∆P emax , proved in Section 3.3.1.When (M E , D E ) and (M P , D P ) are calculated, the smaller one of M E and M P will be selected and, therefore, ∆E max and ∆P emax , energy and power needed in suppressing oscillations, is smaller.The effeteness of the algorithm will be proved in Section 4.3.

Simulations
Simulations based on VSG of single-machine infinite system were performed on the model in DIgSILENT (DIgSILENT/Power Factory 15.1.7,Germany) of Figure 8 with the parameters of f base = 50 Hz, P emax = 0.6 MW, E = 0.3 MW•s, both of which are abstracted from the photovoltaic power plant equipped with storage devices in Qinghai Province, and the renewable energy power was replaced by the DC source.In the simulation, a single lumped model, the replacement of the photovoltaic power plant composed of 34 sets of solar subsystem of the Tucson Power Company, USA is adopted for simplification.Related parameters of the VSG control system and the photovoltaic power station are presented in the Hz, Pemax = 0.6 MW, E = 0.3 MW•s, both of which are abstracted from the photovoltaic power plant equipped with storage devices in Qinghai Province, and the renewable energy power was replaced by the DC source.In the simulation, a single lumped model, the replacement of the photovoltaic power plant composed of 34 sets of solar subsystem of the Tucson Power Company, USA is adopted for simplification.Related parameters of the VSG control and the photovoltaic power station are presented in the Table A1 of Appendix D.

Comparisons between VSG Based on Droop Control and Hamilton Approach
To verify the superiority of Hamiltonian approach, comparisons are made between different VSG control algorithms-the droop control [10] and the Hamiltonian approach.A three-phase metallic short-circuit fault is set at the middle of the transmission line 3-4 at t = 1 s and the fault is cleared at t = 1.05 s.In Figure 9, responding curves of the output active power Pe, grid frequency f, transient E.M.F Eq' and output voltage U0 during fault period are shown in Figure 9a-d  In these figures, solid curves represent the output of VSG with Hamilton approach while the dotted curves represents the output of VSG with droop control.It is obvious that the advanced algorithm enjoys the enhancement in suppressing oscillations, contributing inertia and shortening

Comparisons between VSG Based on Droop Control and Hamilton Approach
To verify the superiority of Hamiltonian approach, comparisons are made between different VSG control algorithms-the droop control [10] and the Hamiltonian approach.A three-phase metallic short-circuit fault is set at the middle of the transmission line 3-4 at t = 1 s and the fault is cleared at t = 1.05 s.In Figure 9, responding curves of the output active power P e , grid frequency f, transient E.M.F E q ' and output voltage U 0 during fault period are shown in Figure 9a-d.Hz, Pemax = 0.6 MW, E = 0.3 MW•s, both of which are abstracted from the photovoltaic power plant equipped with storage devices in Qinghai Province, and the renewable energy power was replaced by the DC source.In the simulation, a single lumped model, the replacement of the photovoltaic power plant composed of 34 sets of solar subsystem of the Tucson Power Company, USA is adopted for simplification.Related parameters of the VSG control system and the photovoltaic power station are presented in the Table A1 of Appendix D.

Comparisons between VSG Based on Droop Control and Hamilton Approach
To verify the superiority of Hamiltonian approach, comparisons are made between different VSG control algorithms-the droop control [10] and the Hamiltonian approach.A three-phase metallic short-circuit fault is set at the middle of the transmission line 3-4 at t = 1 s and the fault is cleared at t = 1.05 s.In Figure 9, responding curves of the output active power Pe, grid frequency f, transient E.M.F Eq' and output voltage U0 during fault period are shown in Figure 9a-d  In these figures, solid curves represent the output of VSG with Hamilton approach while the dotted curves represents the output of VSG with droop control.It is obvious that the advanced algorithm enjoys the enhancement in suppressing oscillations, contributing inertia and shortening In these figures, solid curves represent the output of VSG with Hamilton approach while the dotted curves represents the output of VSG with droop control.It is obvious that the advanced algorithm enjoys the enhancement in suppressing oscillations, contributing inertia and shortening the fluctuation period compared to the VSG based on droop control.Validity can be proved that system stability and robustness are strengthened via the VSG based on a Hamilton approach.

Effects of Different Controller Coefficients on the System Output Responses
To verify the correctness of the controller coefficient boundaries in Section 3.2, simulations are conducted at the same short-circuit condition above.Output responding curves are presented under different coefficients including M, the virtual inertia and D, the system damping.the fluctuation period compared to the VSG based on droop control.Validity can be proved that system stability and robustness are strengthened via the VSG based on a Hamilton approach.

Effects of Different Controller Coefficients on the System Output Responses
To verify the correctness of the controller coefficient boundaries in Section 3.2, simulations are conducted at the same short-circuit condition above.Output responding curves are presented under different coefficients including M, the virtual inertia and D, the system damping.As we can see in Figure 10, it is remarkable in suppress oscillations with a stronger D. When D = 1, it is too weak to suppress system oscillations compared to D = 10.When D = 30, it can be found that an excessive damping constant results in a slower response speed compared to D = 10.When M = 1, it is too small for the inertia to suppress the oscillations.When M = 30, the system tends to be in underdamping condition with a long oscillating period for a small damping ratio.As we can see in Figure 10, it is remarkable in suppress oscillations with a stronger D. When D = 1, it is too weak to suppress system oscillations compared to D = 10.When D = 30, it can be found that an excessive damping constant results in a slower response speed compared to D = 10.

Output Responses with Different M (Virtual Inertia)
In this simulation, D is fixed as 20 while M varies from 1 to 20.Curves are shown in Figure 11a,b.
Energies 2018, 11, x FOR PEER REVIEW 13 of 18 the fluctuation period compared to the VSG based on droop control.Validity can be proved that system stability and robustness are strengthened via the VSG based on a Hamilton approach.

Effects of Different Controller Coefficients on the System Output Responses
To verify the correctness of the controller coefficient boundaries in Section 3.2, simulations are conducted at the same short-circuit condition above.Output responding curves are presented under different coefficients including M, the virtual inertia and D, the system damping.As we can see in Figure 10, it is remarkable in suppress oscillations with a stronger D. When D = 1, it is too weak to suppress system oscillations compared to D = 10.When D = 30, it can be found that an excessive damping constant results in a slower response speed compared to D = 10.When M = 1, it is too small for the inertia to suppress the oscillations.When M = 30, the system tends to be in underdamping condition with a long oscillating period for a small damping ratio.When M = 1, it is too small for the inertia to suppress the oscillations.When M = 30, the system tends to be in underdamping condition with a long oscillating period for a small damping ratio.
Without taking the limits of the storage capacity into consideration, M varies from 0 to 10 while D from 20 to 40 according to Section 3. The results in Figure 10 verify the boundaries of M and D discussed above, while laying a foundation for a further study.

Comparison between VSG Equipped with Capacity-Limited and Limit-Free Storage Devices
Under the same simulation conditions as above, the variation of M is shown in Figure 12a.
Energies 2018, 11, x FOR PEER REVIEW 14 of 18 Without taking the limits of the storage capacity into consideration, M varies from 0 to 10 while D from 20 to 40 according to Section 3. The results in Figure 10 verify the boundaries of M and D discussed above, while laying a foundation for a further study.

Comparison between VSG Equipped with Capacity-Limited and Limit-Free Storage Devices
Under the same simulation conditions as above, the variation of M is shown in Figure 12a.Meanwhile, response curves of output active power Pe and grid frequency f under the storage capacity limits displayed in Figure 12b,c.In those figures, solid curves represent the output of VSG with capacity-limit storage devices while the dotted curves represents the output of VSG with limitfree storage devices.As can be seen in Figure 12a, during the fault period, the adaptive control of VSG takes effect for keeping the system output to stay within the limits.M drops to 4.4 during the fault and returns to the initial value (M = 4.6) when the fault is cleared.It is noticeable that the variation of D is too small to be recorded.Hence, the variation of D will not be given in figures.
In Figure 12b,c, results are expected that VSG with limit-free storage devices performs better in suppressing oscillations and shortening response time than the one with capacity-limit devices.However, with the advanced algorithm, the virtual inertia of the VSG proposed in the paper is adaptive, falls to 4.4 for the limitation in storage capacity during the transient period.As a result, either ∆ max or ∆ emax becomes smaller and less power is needed to contribute from the storage devices.Meanwhile, although the inertia becomes a smaller value, the VSG system can still keep stabilized and the differences from the counterpart is acceptable.In a word, the results prove the validity and effeteness of the adaptive VSG control algorithm and the energy-based control algorithm.Meanwhile, response curves of output active power P e and grid frequency f under the storage capacity limits displayed in Figure 12b,c.In those figures, solid curves represent the output of VSG with capacity-limit storage devices while the dotted curves represents the output of VSG with limit-free storage devices.
As can be seen in Figure 12a, during the fault period, the adaptive control of VSG takes effect for keeping the system output to stay within the limits.M drops to 4.4 during the fault and returns to the initial value (M = 4.6) when the fault is cleared.It is noticeable that the variation of D is too small to be recorded.Hence, the variation of D will not be given in figures.
In Figure 12b,c, results are expected that VSG with limit-free storage devices performs better in suppressing oscillations and shortening response time than the one with capacity-limit devices.However, with the advanced algorithm, the virtual inertia of the VSG proposed in the paper is adaptive, falls to 4.4 for the limitation in storage capacity during the transient period.As a result, either ∆E max or ∆P emax becomes smaller and less power is needed to contribute from the storage devices.Meanwhile, although the inertia becomes a smaller value, the VSG system can still keep Energies 2018, 11, 677 14 of 17 stabilized and the differences from the counterpart is acceptable.In a word, the results prove the validity and effeteness of the adaptive VSG control algorithm and the energy-based control algorithm.

Conclusions
In this paper, an improved VSG model was put forward to enhance system suitability and robustness.A corresponding strategy based on Hamilton theory was then applied to meet the complexity of the model.Compared to the VSG based on droop control via simulations, the advanced model performed better in suppressing oscillations, contributing inertia and shortening the fluctuation period.Boundaries of virtual inertia and damping corresponding to the system stable region were deducted by the study on the small-signal model of VSG.Simulations on the effects that coefficients variation exert on system output responses proved the correctness of the boundaries.To develop the control algorithm with the limits in storage capacity, relations between VSG parameters and charge and discharge capacity were analyzed.The dynamic response of the virtual inertia showed that the proposed VSG made an adaption to satisfy the capacity restriction.With the acceptable cost of losing part of the ability in suppressing oscillations during fault period, the simulations based on improved model worked in the condition closer to the one in practical engineering.

Figure 1 .
Figure 1.Block diagram of the VSG control strategy.

Figure 1 .
Figure 1.Block diagram of the VSG control strategy.

Figure 2 .
Figure 2. Roots of the poles with the variations of the parameters.(a) Pole Loci when M a constant value; (b) Pole Loci when D a constant value.

Figure 2 .
Figure 2. Roots of the poles with the variations of the parameters.(a) Pole Loci when M a constant value; (b) Pole Loci when D a constant value.

18 Figure 3 .
Figure 3. Step response with different damping ratio.(a) Active power response in under-damping condition (b) Active power response in over-damping condition.

Figure 3 .
Figure 3. Step response with different damping ratio.(a) Active power response in under-damping condition (b) Active power response in over-damping condition.

Figure 4 .
Figure 4. Relationship between the Capacity limits of the energy storage equipment and the parameters.(a) Storage capacity limits with M; (b) Storage capacity limits with D.

Figure 4 .
Figure 4. Relationship between the Capacity limits of the energy storage equipment and the parameters.(a) Storage capacity limits with M; (b) Storage capacity limits with D.

Figure 5 .
Figure 5. Relation of M and D under different storage limits.(a) Relation of M and D when ∆ emax = Const; (b) Relation of M and D when ∆ max = Const.

Figure 5 .
Figure 5. Relation of M and D under different storage limits.(a) Relation of M and D when ∆P emax = Const; (b) Relation of M and D when ∆E max = Const.

4. 2 . 1 .
Output Responses with Different D (Virtual Damping) In this simulation, M is fixed as 7 while D varies from 1 to 30.Curves are shown in Figure 10a,b.Energies 2018, 11, x FOR PEER REVIEW 13 of 18

4. 2 . 2 .
Output Responses with Different M (Virtual Inertia) In this simulation, D is fixed as 20 while M varies from 1 to 20.Curves are shown in Figure 11a,b.

Figure 11 .
Figure 11.Output responses with different virtual inertia under short-circuit fault condition.(a) Comparisons of Pe/p.u.; (b) Comparisons of f/p.u.

Figure 10 .
Figure 10.Output responses with different damping under short-circuit fault condition.(a) Comparisons of P e /p.u.; (b) Comparisons of f /p.u.

Figure 11 .
Figure 11.Output responses with different virtual inertia under short-circuit fault condition.(a) Comparisons of P e /p.u.; (b) Comparisons of f /p.u.

Figure 12 .
Figure 12.Figures of system output with limited/limit-free energy storage capacity.(a) Variation of the virtual inertia; (b) Comparison of Pe/p.u;(c) Comparison of f/p.u.

Figure 12 .
Figure 12.Figures of system output with limited/limit-free energy storage capacity.(a) Variation of the virtual inertia; (b) Comparison of P e /p.u; (c) Comparison of f /p.u.
depicts the overall block diagram of VSG control strategy.