Small-signal stability analysis of Energy Internet through differential inclusion theory

: In this study, small-signal stability issue associated with uncertainty excitation is investigated by using differential inclusion theory. Specifically, a polytopic linear differential inclusion (PLDI) model for Energy Internet is developed, in which a modified non-sequence Monte Carlo method is introduced to identify a series of time-variant operation states. Additionally, a simplified small-signal model of renewable energy sources (RES) with virtual synchronous generator control is proposed, and the outputs of RES are modelled as time-varying elements in PLDI model to reflect the inner stochastic excitation in the linearised matrix. The stability criterion for the stochastic time-varying system is mathematically deduced based on convex hull Lyapunov function. Simulation demonstrates the benefits of the proposed model in describing system stochastic characteristics and reducing computational burden.


Introduction
Unpredictable events and disturbances caused by high penetration of renewable energy sources (RES) have brought new challenges to small-signal stability analysis of Energy Internet [1]. To address the challenges, the following two research efforts have been elaborated by the engineering community to identify uncertainties in stability analysis: (i) stochastic time-varying equilibrium points [2] caused by time-varying power flows with RES output fluctuations, and (ii) stochastic state matrix coefficients also caused by fluctuations of RES but may not affect steady-state operating.
The emerging issue, usually based on several critical operation states, is widely studied. Based on a certain probability distribution function (PDF), a large number of scenarios can be generated by Monte Carlo (MC) [3] or latin hypercube [4] method for deterministic simulation. Gram-Charlier expansion and cumulants methods are utilised in [5] to obtain scenarios and identify the distribution density of the critical eigenvalues so as to derive the possibilities of system stability. Those methods, however, are still based on the deterministic model of the Riemann differential equation that ignores the uncertainty during time series. Therefore, stochastic differential equation (SDE) theory [6] is introduced to better describe the dynamic behaviours of the intrinsic nature of the stochastic components and hence demonstrates the stability criterion for systems with external excitation. However, the method ignores the inclusion of perturbations. For this concern, the authors in [7,8] developed a 'region-wise' approach to probe into the impact of uncertainty on large-scale system stability. Nevertheless, the consideration of varying equilibrium point is ignored in these approaches. Recently, enabled by massive historical data, random matrix theory is introduced to dig out the evolution law of system operating trajectory and predict the tendency of operation state [9,10]. Need not to consider varying equilibrium point issue, this promising method provides a possibility to grasp the behaviour of a power system from the data-driven perspective without accurately establishing detailed system models. However, it still needs some works on establishing mathematic stability analysis models to figure out and characterise intrinsic regularity of a stochastic system.
To address the issues mentioned above, this paper proposes an advanced small-signal model for Energy Internet based on linear differential inclusion (LDI) theory. The LDI only requires a motion-changed region of the uncontrolled resource in the time period [11]. This advantage provides a possibility to apply a stochastic forecast model to deal with the stability analysis with rational prediction errors. Meanwhile, the uncertainty system can also be described by a set-valued mapping of operating states. In [12], polytopic LDI (PLDI) is adopted to achieve a robust control design of flexible AC transmission systems (FACTS) for the power system. Through a linearising system under a number of typical conditions, a satisfactory global optimum operation condition can be obtained. In [13], the definition and property of composite quadratic hull Lyapunov functions (CHLFs) are investigated. This method is subsequently applied to construct a stabilising feedback law for saturated systems, such as AC-DC converters [14] and battery boost converter [15]. Inherited from this innovative idea, this paper tackles the stability analysis problem due to the following salient features: i. Introducing an operation-point identification method based on a modified non-sequence MC method, which enables a stochastic time-varying PLDI model for system analysis. ii. Proposing a simplified small-signal stability model for RES, considering virtual synchronous generator (VSG) control strategy. iii. Demonstrating a stability criterion for Energy Internet based on CHLF, and a resultant stochastic-stability constraint is established for robustness control.
This paper is organised as follows. Sections 2 and 3 theoretically present the principle of LDI as well as the simplified small-signal model of RES and introduce the operation-point identification method to establish the time-varying PLDI model. The stability criterion based on CHLF is demonstrated in Section 4. Section 5 is devoted to the simulation analysis of a two-machine system. Section 6 presents conclusions and future works.

Simplified small-signal model of RES
For simplicity and without loss of generality, this paper mainly focuses on solar and wind resources. The linearised mechanical powers of wind power P wg and solar power P pv are characterised as follows: where ρ is the air density coefficient kg/m 2 , R is the length of wind turbine blade, V w is the wind velocity. For PV, A and G r are the respective array of photovoltaic solar panels m 2 and solar radiation W/m 2 , respectively. The conversion efficiency parameters of wind generator and PV under normal operation state are set as C p λ, β = 47%, η w = 0.97, η pv = 15%.
Since the time-scale of converter controllers is inconsistent with the mechanical adjustment of synchronous machines (SMs), Jaime et al. [16] proved that the integration of RES will not involve in the original electromagnetic oscillation mode (EOM) of system. This conclusion indicates that the corresponding SMs replaced by the RES in EOM do not need detailed model converter control strategies and provide a possibility to represent RES as simple as a time-varying controlled current source. However, when the phaselocked loop (PLL) controller is considered, the RES becomes coupled with dynamic characteristics especially in a weak network [17]. Therefore, the simplified small-signal stability model for RES should include the PLL: where K p_PLL and K i_PLL are the proportion integration (PI) parameters of the PLL controller. Assume the controller is ideally accurate, then U q ≃ U, which is the terminal voltage of integrated node and meets the steady equation if the network appears inductive: where P sys is the system electrical output governed by generator.
P ref and E are the respective output and terminal voltage of RES, and X out is the equivalent impedance between the two terminals. Then, according to (3) and assuming P sys = constant, ΔU q ≃ ΔU should satisfy In (4), ΔP ref satisfies (1) and δ 0 is the initial rotor angle of system. Then, (2) can be expressed as As traditional current/voltage source inverters lack inertia, VSG control algorithm shown in Fig. 1 is introduced to mimic SM characteristics. It is able to provide frequency and voltage support during power fluctuations or failures [18]. As the VSG control strategy results from the movement of SM, the VSG-based converters participate in EOM of system and hence influence the small-signal stability of power systems [19]. Therefore, this type of controller should be considered in the small-signal model of RES. Combined with ω − P droop strategy and assuming P ref = P wg or P ref = P pv , the linearisation model for a VSG yields In (3), ω and ω n are operational and rated rotor angular velocity, respectively. J is the virtual inertia, D is the damping coefficient, and K ω is the active droop coefficient. Assume the power system appears inductive according to (3), the linearised form of the real powers is expressed as In (5) and (6), ΔP ref meets (1) and can be expressed as where ΔP refivar and ΔP reficos are the varying and constant parts of electromagnetic power provided by RES, respectively. x i represents the external excitation that is either wind velocity V w or solar radiation in this case. x imin represents the minimum stochastic external excitation. Equation (8) indicates that the external excitation can be divided into the constant part and the variable part. Hence, (8) can be represented as a matrix with a bounded time-varying region. Note that the controller often uses a PLL to synchronise its output voltage with the grid voltage before actually connecting to the grid. In current control approaches, the PLL remains in the control loop all the time. Therefore, especially under weak grid conditions, it can cause oscillations and stability issues for the system. This issue does not exist in VSG-type controls where the PLL (even if it exists initially, it) does not remain in the loop during normal operation. Once the switching is done, the VSG control is properly initialised and activated to ensure a smooth power output and the PLL can be removed [20]. Therefore, the PLL is only for one-time use and is not considered in the smallsignal model of RES. Hence, the small-signal model of renewable generation is revised as According to (8), the detailed matrix identified by (9) can be divided into certainty and uncertainty parts. Then, the constant state variables and control variables take wind farm, for instance, with respect to (9) The related matrices of RES based on PLL and VSG are (see (10)) x t : 0, T → R n is an unknown function satisfying x 0 = x 0 , then a system in the following form is called LDI: If the set-valued mapping F in R n → R n does not explicitly contain T, we have a time-invariant system LDI as follows: Γ ∈ R N represents a closed convex set in the first quadrant, and then the LDI system with the set co A is called PLDI. For a typical time-varying system containing uncertainty element, it meets the following form: where x ∈ R n and u ∈ R m are the respective state and input variables of the system, respectively. A i and B i are the constant matrices of the system. K i σ and H i σ represent the uncertainty of system, which satisfies the following: In (17) and K ir ∈ R n × n and H ir ∈ R n × m are the respective constant-varying matrices at each state according to uncertainty functions σ i , i ∈ 1, N . Theorem 1: Suppose a set-valued mapping F: R n → R m is a bounded closed convex, solutions in both (10) and (11) exist. Denote solutions set R x 0 and the LDI system is exponentially stable, only if ∀x t ∈ R x 0 , ∃β > 0, c > 0 satisfies x t ≤ cexp −βt . For a typical PLDI, the system has a globally uniform stability when it satisfies exponentially stable Theorem 1 [13]. This criterion indicates that, for a certain closed compact convex set, the global stability of this PLDI system can be estimated by analysing the boundary property of the convex set. This stable characteristic provides a possibility to apply PDLI theory to qualitatively judge the stability of a stochastic time-variant system, and simplifies the analysis process of uncertainty disturbances, further improving the modelling accuracy. For instance, the essence of the kernel function of multi-equilibrium linearisation means in [21] can be theoretically elaborated by PLDI.

Modelling of Energy Internet based on PLDI
For a given perturbation, e.g. wind velocity V imin , V imax , power system exists a series of representative equilibrium points. This characteristic indicates that the range of excitation, such as light intensity, can be divided into several intervals with respect to its equilibrium points, and also indicates that the operation state of interconnect system is constant under the external excitation of each interval of perturbation.
Therefore, the stochastic time-variant power system can be represented as a series of typical states with their internal excitations by utilising (9) and (16). The system is described as a convex hull of multiple operation states, and the perturbation excitation at each interval is presented as bounded excitation in each state. Then, to obtain a series of system operation points under large fluctuations caused by RES, a method that judges the trigger of SM additional adjustment should be developed.
Instead of calculating the PDF of uncertainties at each bus, a modified non-sequence MC method is introduced to identify various operation states of the uncertainty system. According to frequency constraints, the deviation threshold of additional adjustment is set to ± 0.2 Hz, which triggers the secondary frequency modulation [22]. The adjustment of SM and additional device will change system equilibrium point and hence the state of the system. As the procedure is similar to the bulk power system reliability evaluation, this proposed algorithm is based on heuristic adjustment strategy, which is summarised as follows: (a) Generate stochastic excitation and the mechanical power input of RES using MC; if it is NOT the first iteration, go to step (g). Using the above-mentioned method, a number of typical operating states of power system and the corresponding interval of excitations can be identified. Then, form (16) can be introduced to represent these random variables and formulate vertex systems of a PLDI model. As mentioned earlier, RES and synchronous system are mainly coupled with each other through power flow. However, for each equilibrium point identified by non-sequence MC method, the excitation of RES under this scenario does not trigger the movement, which indicates the independence of the power system and RES. Similar to [6], the polytopic linear differential matrices A i and B i in (16) yield Equation (18) includes linearised interconnect system matrices A sysi and B sysi under operation state i. Based on the Heffron-Philips model, the classic third-order small-signal model of SM is introduced and A sysi and B sysi are described as follows [22]: A RES cos i and B RES cos i are the constant parts of the small-signal model with stochastic excitation in (18). Then, the convex combination of linear differential matrices A i and B i identified by the modified non-sequence MC method will be a closed set.
The time-varying elements in PLDI can be obtained by using (9), (16), and (17). In this case, the varying parts in the RES model (9) only constitute uncertainty matrix H i σ . Then, for each external variable, let σ i j ∈ x i min , x i max , the matrices H ir ∈ R n × m can be rewritten as (20) In (20), x imin and x imax represent the maximum and minimum of the stochastic external excitation, i.e. wind velocity and light intensity. Then, with the above-presented model, the small-signal model based on PLDI can be expressed as (16).

Definition and general properties of CHLF
For a set of positive-definite symmetric matrices P i > 0 ∈ R n × n , , if a vector γ ∈ Γ and satisfies (15), then define [13] It is easy to see that both Q γ > 0 and P γ > 0 are continuous surjective linear mappings in γ ∈ Γ. CHLF is defined as [13] V c (x) = min Obviously, V c x is quadratic and satisfies V c αx = α 2 V c x . According to the Schur complement rule, we have Then, we give a simple preliminary that will be used throughout this section [13]. Lemma 1: Let x ∈ bdL c . For simplicity and without loss of generality, assume γ k * > 0 for k ∈ 1, N 0 and γ k * = 0 for k ∈ N 0 + 1, N . Combine (21) and denote Q γ x k = Q k Q γ * −1 x, k ∈ 1, N 0 , then we have [13] γ * x = arg min

Stability criterion for PLDI-based CHLF
Consider x ∈ bdL V c and F k ∈ R m × n both satisfy Lemma 2, for simplicity and without loss of generality, assume F γ * = ∑ k = 1 N 0 γ k * F k and construct u = F γ * Q −1 γ * x [13]. Then, substituting (25) and (26) into (22) yields Hence, as According to the properties of the symmetric quadratic matrix, (28) can be represented as Hence, for the linear control system, the problem V˙x < 0 is equivalent to Multiply (32) from the left and from the right with Q k , the PDLI system will be stable only if it satisfies

Tests and results
To validate the proposed approach, a two-machine bus system test case is utilised. The case is shown in Fig. 2. In this case, a wind farm based on a permanent magnet synchronous generator (PMSG) of 10% equivalent MVA rating of the SM is added in bus 1. The specification and setting of PMSG are presented in Table 1. The wind velocity range is set as V m ∈ 8 − 6 m/s and V m ∈ 8 − 2 m/s, respectively. As previously mentioned, the stability analysis procedure can be summarised as follows: Step 1: Calculate stochastic system equilibrium points and identify respective fluctuation range of excitation resource for each operation state.
Step 2: Generalise a series of linearised matrices according to system operation states. Describe the uncertain system using the PLDI model with time-varying elements.
Step 3: Estimate the stability of the stochastic interconnect system based on CHLF.
According to power flow analysis results, there are two equilibrium points that correspond to two fluctuation intervals, i.e. V wind1 ∈ 8 − 10.47 m/s, V wind2 ∈ 10.47 − 12 m/s. The bounds of voltage are shown in Fig. 3, and the active power of SM is limited to 0.72, 0.78 p . u .
The calculation results with or without the VSG control strategy are shown in Tables 2 and 3, respectively. According to Theorem 2, as exist and only exist positive-definite matrices Q 1v , F 1v , Q 2v , F 2v , under VSG-based RES satisfying the above inequalities, the stochastic PLDI system is exponentially stable. However, for PLL control, there only exist Q 1p , F 1p , Q 2p , F 2p , during V wind1 ∈ 8 − 10.47 m/s. This result indicates that, without auxiliary control, the PLDI system is not stable under PLL control during V m ∈ 8 − 12 m/s.
To verify the consistency of theoretical analysis and numerical simulation results, an electromagnetic transient model is established using MATLAB/SIMULINK. A step change from 8 to 16 m/s in wind velocity happens at 2 s. To adequately demonstrate the validity of the proposed method, the output of the PMSG after step change is modelled as a persistent disturbance and meets the four component-model [4] shown in Fig. 4. Fig. 5 shows the dynamic response of respective state variables under the VSG control strategy using simulation. For traditional PLL control, the system is not stable after the applied disturbance shown in Fig. 6. According to the simulation results, the VSG-based system is stable under fluctuation. In contrast, the PLL -based system is unable to reach a desired steady-state operating point. This finding is in accord with the theoretical analysis.  To highlight the benefits of this advanced model in stability assessment and simulation speedup, two test cases including case 1 is introduced to compare with MC and SDE models. The additional case, i.e. case 2, is based on IEEE 145-bus system, and the detailed system parameters can be found in [6]. The magnitudes of small disturbance caused by wind power are set to 2.0 for these three test systems. The simulation is carried out on a personal laptop (Intel   Core CPU i7-5500U 2.40 GHz, 16 GB RAM). The results are presented in Table 4. According to Table 4, both MC and SDE methods yield inconsistent results with the PLDI method in case 1. Meanwhile, the previous simulation has proved the instability of the PLL-based system in case 1. As those types of methods are based on constant equilibrium-point assumption, they cannot accurately reflect all the dynamic responses caused by RES. Besides, as shown in Table 4, the proposed method has an advantage of reducing computational burden in MC simulation. In contrast, though the SDE method has a faster convergence rate, it requires engineering experiences to quantify stochastic excitation, and this process may be prone to introduce man-made error. The results of this case study illustrate the effectiveness and correctness of the proposed algorithm in both small-scale and medium-sized systems.
To be mentioned, typical power system models are often composed of several hundreds of state variables. One disadvantage of this proposed method is that an iterative BMI solver with high dimension may be time consuming. Of course, the computational cost of optimisation has to be further reduced to be applicable for practical use. Therefore, the technique presented in [7] can be utilised to generate the reduced-order models of real-sized power systems. Hence, the previously mentioned model in this paper can be successfully applied to a large-scale power system.

Conclusion
This paper has introduced differential inclusion theory to investigate the static angle stability. Based on theoretical analysis and simulation study, the major conclusions include the following: (i) The authors delicately derive the linearised model of a stochastic time-varying power system containing RES in accordance with VSG-based and PLL-based strategies by utilising the PLDI model, and the corresponding asymptotic stability criterion is proposed. (ii) Compared to MC and SDE models, the proposed model can accurately reflect the impact of the inner stochastic excitation and time-varying operation states of system caused by RES. Meanwhile, analysis efficiency is significantly improved.
Admittedly, concerning the proposed small-signal the statespace model, a few aspects should be enriched, such as the design and optimise robust controller by utilising asymptotic stability criterion, and the process optimisation of BMI systems calculation. These research tasks will be carried out in future.