Markov Jump Linear System Analysis of a Microgrid Operating in Islanded and Grid-Tied Modes

—The analysis of power system dynamics is usually conducted using traditional models based on the standard nonlinear differential algebraic equations (DAEs). In general, solutions to these equations can be obtained using numerical methods such as the Monte Carlo simulations. The use of methods based on the Stochastic Hybrid System (SHS) framework for power systems subject to stochastic behavior is relatively new. These methods have been successfully applied to power systems subjected to stochastic inputs. This study discusses a class of SHSs referred to as Markov Jump Linear Systems (MJLSs), in which the entire dynamic system is jumping between distinct operating points, with different local small-signal dynamics. The numerical application is based on the analysis of the IEEE 37-bus power system switching between grid-tied and standalone operating modes. The Ordinary Differential Equations (ODEs) representing the evolution of the conditional moments are derived and a matrix representation of the system is developed. Results are compared to the averaged Monte Carlo simulation. The MJLS approach was found to have a key advantage of being far less computational expensive.

Abstract-The analysis of power system dynamics is usually conducted using traditional models based on the standard nonlinear differential algebraic equations (DAEs). In general, solutions to these equations can be obtained using numerical methods such as the Monte Carlo simulations. The use of methods based on the Stochastic Hybrid System (SHS) framework for power systems subject to stochastic behavior is relatively new. These methods have been successfully applied to power systems subjected to stochastic inputs. This study discusses a class of SHSs referred to as Markov Jump Linear Systems (MJLSs), in which the entire dynamic system is jumping between distinct operating points, with different local small-signal dynamics. The numerical application is based on the analysis of the IEEE 37-bus power system switching between grid-tied and standalone operating modes. The Ordinary Differential Equations (ODEs) representing the evolution of the conditional moments are derived and a matrix representation of the system is developed. Results are compared to the averaged Monte Carlo simulation. The MJLS approach was found to have a key advantage of being far less computational expensive.
Index Terms-Markov Jump Linear Systems, Monte Carlo simulations, Stochastic Hybrid Systems, Power System statistics.

I. INTRODUCTION
T HE objective of this paper is to analyze a microgrid using advanced techniques that apply to Markov Jump Linear Systems. The system under study is a modified IEEE 37-bus power system where seven inverters were added at selected buses [1] to represent Distributed Energy Resources (DERs). The application of stochastic methods to the analysis of microgrids has gained increased attention with a lot of research papers focused on the control and stability analysis of power systems subject to random disturbances [2], [3], [4]. However, previous works did not consider large signal disturbances such as those resulting from a microgrid switching between different operating equilibrium points.
In this paper, a microgrid is defined as a small-scale power grid that can operate independently or collaboratively with other small power grids [5]. This architecture based on the use of microgrids in a power system is also known as distributed, dispersed, decentralized, district or embedded power generation. Any small-scale and localized power system that includes localized generation resources and, but not always, storage capability, with well defined boundaries can be identified as a microgrid. In the case a microgrid is allowed to integrate with the main power grid, it is referred to as a grid-connected or hybrid microgrid. Typically, power production in microgrids is supported by generators or renewable energy systems such as wind and solar resources, and when integrated with the main grid they provide backup power or supplemental power during periods of heavy demand. This strategy provides redundancy for critical and essential services and makes the main grid more immune to temporary disasters [5].
In real-time operation, the control of distributed energy resources in a microgrid and the preservation of an adequate system inertia pose significant challenges to the stability of microgrids' operation. This is particularly true in the absence of a stiff grid, when a microgrid is operating in standalone or islanded mode. Because of the combination of factors such as low microgrid inertia, limited power storage capability and tight coupling among various elements in the system, fluctuations in the output power of distributed energy resources or changes in local load may lead to power quality or frequency/voltage stability concerns [6]. For example, if an important load is switched on and off, frequency and voltage may change dramatically, and could even result in a collapse of the entire system. Additionally, with an increasing penetration of microgrids, the main power grid could face a similar challenge when subject to active/reactive power injections resulting from a random coupling and decoupling of multiple DERs [2].
To characterize the randomness in system operation as described above, methods based on Stochastic Hybrid Systems were proposed [2], [7], [8]. These methods combine the conventional power system dynamic model and a stochastic representation. In this study, focus is on one particular class of SHSs referred to as Markov Jump Linear Systems. An MJLS is composed of two coupled sub-systems [9].
The first sub-system is based on a linearized power system dynamic model. Each mode of the stochastic model is mapped to a linear system with continuous state variables. These linear systems correspond to continuous-time small-signal models of the microgrid linearized around well defined operating points.
The second sub-system is a continuous-time Markov process, which has a defined number of discrete states, N , also described as stochastic modes of the Markov process. These discrete states correspond to the steady-state operating points of the microgrid. Hence, the system is described as "jumping" between different modes (i, j) at a transition rate λ ij .
The two coupled models are represented in Fig. 1. The general SHS model framework for a power system is discussed in [2], where a power system is subject to stochastic power injections resulting from the coupling and decoupling of microgrids. The power system is operating around an identified stable operating point, only the control variables are randomly switching between a finite set of discrete values. The method results in a set of ODEs that describe the evolution of the conditional moments of the power system. The integration of these ODEs yield the statistics of the dynamic and algebraic states. This method is also used in [7] to predict the influence of random load behavior on the dynamics of dc microgrids and distribution systems.
Following a method described in [2], a set of ODEs representing the evolution of the conditional moments is derived to represent the dynamics of a microgrid system oscillating between distinct operating equilibrium points, with different local small-signal dynamics.
For the system under study, the IEEE 37-bus microgrid, the analysis begins with a description of the power system dynamics by standard nonlinear differential algebraic equations. Based on the derivation in [1], [10], a total of 225 dynamic states and a collection of non-linear equations are identified to represent the microgrid dynamics. These non-linear equations are then linearized around equilibrium points, and the result is a 225 th order state space representation. The system is considered switching between three equilibrium points (one grid-tied mode and two standalone modes). Therefore, there are three state space representations of the 225 th order corresponding to each operating mode. However, the linearized system is a two-time scale system that combines fast and slow dynamic states. By applying the singular perturbation method [1], [11], slow states can be separated from the fast states. The procedure results in a system order reduction to 56 dynamic states. The stochastic modelling consists in characterising the random switching of the system from one operating mode to another. This switching behavior is appropriately represented as a Continuous Time Markov Chain [2], [8], [9]. In practice, this corresponds to including a jump term in the state equations to represent jump events at Poisson distributed jump times [8].
To solve this combined system (linear dynamics and jumps), there are two possible methods: the numerical integration of the system of ODEs representing the dynamics of the conditional moments, or the averaged solution of repeated Monte Carlo simulations. The solutions obtained using either method represent the conditional moments of the stochastic system. This study computes these moments using the two methods mentioned above and demonstrates that the MJLS results converge to the averaged Monte Carlo simulations. The steps discussed above are represented in Fig. 2. This paper is structured as follows. Section II describes the system under study: the modified IEEE 37-bus power system. In Section III, the power system model is cast as an MJLS and, using the equations developed in [2], [12], a matrix representation of the system is derived for the computation of the conditional moments. In Section IV, the solutions to the ODEs of the system dynamics are computed and discussed using the two methods described above (MJLS and averaged Monte Carlo simulations). Conclusion and avenues for future work are summarized in Section V.

II. SYSTEM UNDER CONSIDERATION
The method developed in this study is applied to the IEEE 37-bus power system. The standard IEEE 37-bus system was modified with the addition of seven inverters at selected nodes to represent inverter-based DERs. The one-line diagram is represented in Fig. 3, where the larger dots at designated bus locations represent the inverters. Static loads are spot loads at various nodes and represent consumers of constant power (P Q), constant current or constant impedance. Details of the modified system including the load and line parameters are given in [1], [11].
For the purpose of this study, the system is considered switching between the grid-tied and two distinct islanded operating modes. The connection to the main power grid is made at the point of common coupling (PCC). In the gridtied mode, all the bus voltages and the system frequency are maintained by the stiff main grid. However, in the islanded mode, voltage and frequency are controlled by the DERs using the droop control strategy [13]. The numerical application considers one grid-tied operating mode, and two distinct islanded (standalone) operating modes. The two standalone operating modes correspond to two distinct equilibrium points characterized by two distinct sets of bus loads and inverter settings. A sample of numerical values for the dynamic states of the IEEE 37-bus microgrid operating in grid-tied and standalone modes are presented in Table 1. The definition of all parameters can be found in [1], [10], [11]. The small-signal model of the overall system was developed in [1], [11] using the dq reference frame and resulted in a 225 th order system. The derivation of the linearized dynamic model is also presented in details in the same papers. Each inverter system contains 15 states. Load and line models contain two states: the load currents (d, q) and the line currents (d, q). These dynamic states are presented below for the islanded operating mode For the grid-tied system [10], [11], the inverter states φ P and φ Q in the voltage controller of the islanded system are replaced with φ d and φ q in the power controller, respectively. It can be shown that in the two cases (grid-tied and islanded), the system matrices are compatible and therefore a transition islanded mode to grid-tied mode (and vice versa) can be envisaged. In the standalone mode, φ d represents the integral of the error between the angular frequency ω P LL generated by the PLL and the reference frequency ω * . Similarly, φ q is equivalent to the integral of the error between the reference voltage and the measured q-axis voltage v oq . In the grid-tied mode, φ P is the integral of the error between the reference active power P * and the measured active power P , in the power controller, and φ Q is the integral of the error between the reference reactive power Q * and the measured reactive power Q.
Another difference between the two operating modes is in the inputs to the system. In the grid-tied mode, the inputs are the active and reactive power references (P * , Q * ). However, in the islanded mode, the inverters are droop controlled and the inputs can be defined as the bus voltages. As explained in [11], to accurately predict the effects of load perturbations when bus voltages are used as inputs, a method based on a virtual resistor is used. The virtual resistor with high resistance is connected at the inverter bus and has a negligible effect on the system dynamics. With this method, the terms related to the bus voltages are included in the system matrix.
An analysis of the dynamic model reveals the existence of a two-time-scale behavior [14], which requires small time steps for the fast dynamics and a long simulation time for the slow dynamics. Using the singular perturbation method [1], [11], fast transients can be eliminated and, as a result, the system order is reduced. In the case of the modified IEEE 37-bus microgrid, the system order is reduced from 225 to 56 dynamic states [1], [11]. The eliminated fast states can always be recovered through algebraic equations as linear combinations of the slow dynamics states.
Following the procedure described above: non-linear model identification, linearization around steady-state operating points, system order reduction, and with help of the Matlab symbolic toolbox, the small-signal mathematical model of the grid-tied system isẋ where the control vectorũ = [P * Q * ] T . In the case of the islanded operating modes, the small-signal dynamic model iṡ where, as indicated above and in [10], [11], the input vector of bus voltages has been expressed in terms of the dynamic states.

III. IEEE 37-BUS POWER SYSTEM AS AN MJLS
In this section, the IEEE-37 bus power system is framed as an MJLS. The system is switching between three operating modes in a random manner. The conjoint modeling of the linear dynamic system and the stochastic system was developed in [2], [7], [8], [15]. The resulting system of ODEs that represent the evolution of the conditional moments were derived in [2]. The ODE for the m th order conditional moment with respect to the stochastic mode i is represented aṡ denotes the time derivative of the conditional moment of the m th order associated with mode i, m p is the p th element of the moment index vector m = (m 1 , ..., m N ), a pr are elements of the state matrix A i , v i p represents the elements of the control vector augmented with an affine term, λ ji is the transition rate from j to i, Q − i is the set of incoming transitions, and Q + i is the set of outgoing transitions, e p , e r are unity vectors with a 1 at the p th (r th , respectively) position and 0 elsewhere. The solutions to the MJLS model can be computed using two approaches. The first method is a discrete time approximation using the averaged Monte Carlo simulation, and the second method is based on a matrix representation of equation (6). The two approaches are explained below.
1) Averaged Monte Carlo simulation: The method consists in solving the dynamic equations in (4) or (5) numerically along the paths resulting from the CTMC. Therefore, the procedure begins with generating the CTMC integration paths based on a given transition rate matrix. The procedure to generate these integration paths is summarized in Algorithm 1, where j = max(j) and the operation: mode = k + 1 is defined such as 1 + 1 = 2, 2 + 1 = 3, and 3 + 1 = 1, for k = 1, 2, 3, respectively. The method described in Algorithm 1 is an adaptation of the Gillepsie's Stochastic Simulation Algorithm (SSA) [16]. The next step is to solve the dynamic states using the Euler method along each path generated by Algorithm 1. This procedure is repeated during successive runs, the number of which will determine the accuracy of the final result. The final solutions for the dynamic and algebraic states are calculated by averaging the results from all runs. This yields the statistics of the system: zeroth, first and second moments, which are usually sufficient to characterize a distribution. The procedure for the repetitive MC simulations is described in Algorithm 2

Algorithm 2: Repetitive Monte Carlo simulations
Require: T s , T f inal , M ode sequence f rom Algorithm 1 Initialize: max iter, X init , x iter = 0 for iter = 1 : max iter do for k = 1 : length(t) do 2) Matrix representation of the MJLS: For small systems, with a very small number of dynamic states (for instance the SMIB system described in [2]), an analytical solution of the equation above is pretty straightforward. However, for larger systems with a higher number of dynamic states, an analytical solution becomes quickly very complex. It was discussed in [17], [18] that, for a larger system, the solutions are more conveniently computed when equation (6) is put in matrix form. This matrix representation is derived below, following the methods discussed in [17], [7]. Fortunately, the system of ODEs (6) representing the evolution of the conditional moments are finite-dimensional. Each term on the left hand side of (6) depends only on moments of equal order or lower on the right hand side. Therefore, moment-closure methods are not necessary. It can be seen from (6), thatμ      (t), associated with modes i and j are of the same dimension m. Regrouping all terms associated with moments of the same order [7], the system of ODEs represented in (6) can now be expressed as followṡ This expression is equivalent to (6) and matrices G (m) and H (m) can be derived using the method detailed in [7]. They are essentially block diagonal combinations of sub-matrices related to each mode of the stochastic system. The general structure of G (m) and H (m) are provided in (8) and (9). G (m) is a sum of two terms. The first term is a block diagonal matrix which elements G  H (m) is also a block-diagonal matrix which elements correspond to the input term of equation (4) augmented with an affine term: The method to derive these matrices is described in length in [7], [17] and succinctly presented here for the lower order moments: zeroth, first and second.
The 0 th order moments are the occupational probabilities of the stochastic model, i.e., they represent the probability for the system to be in a particular mode i. To obtain the 0 th order moments, we replace m = 0 in equations (6) and (7). It follows It is worth noting that µ (0−1) is not permitted as it would result in a negative order moment.
The first order moments are the statistical means of the system. It can be shown that The second order moments are the uncentered second moments of the system. G (2) and H (2) elements result from the Kronecker product of I(N ) and A i (respectively, v i ), subsequently multiplied by W m , which is a transformation matrix that describes the structure and ordering of the second moment [7] The next step in the calculation of the second moment is to eliminate the redundant second order moments. The procedure results in the reduction of the order and size of G (2) and H (2) . After the redundant second order moments are eliminated, the equations (8) and (9) can now be rewritten as follows: It follows that the second order matriceŝ In (16), (18), and (19), N u (N ) is the binomial coefficient Finally, the expressions for the low order moments are summarized as follows: Distributions are usually described with means and variances. For a stochastic variable, X, the variance is derived from the uncentered second moment as follows An example is provided below to illustrate the derivation of the low order moments (0 th , 1 st , and 2 nd ) for a small system (2 states, 2 modes).
Example: The 2-state, 2-mode system were discussed at length in [2]. The matrix representation for this small system is shown here.
For the 0 th order moment (occupational probabilities), the matrix representation of the ODE (6) is The computation of equation (23) requires the knowledge of an initial condition. Assuming the system is in mode 2 at time t = 0, the initial condition is [0 1] T . The solution to (23) yields a 2 × 1 vector The evolution of the 1 st order moments is represented in the matrix equation (25) The first order moment is a n × 1 vector that represents the statistic means of the dynamic states, and is computed as The matrix equation of the 2 nd order moments (uncentered second moments) is too large to be include here. The second order moment is an n(n + 1)/2×1 vector whose elements are the uncentered second moments of the dynamic states, and is computed as The computation of the ODEs for the first order and second order moments requires the knowledge of the transition matrix, which contains the transition rates between different stochastic modes. These rates are either calculated or evaluated based on empirical observations of the system. In this study, the transition rates are considered given. For a 3-mode system, the incoming and outgoing transition rates are represented by a 3 × 3 matrix Notice that the diagonal elements are equal to zero. To obtain the full transition matrix, the self-transition elements λ ii are added in such a way that the sum of row elements is equal to zero. This matrix Λ should not be confused with the Transition Probability Matrix (TPM) in a discrete-time Markov process where the sum of row elements is equal to 1. The full transition matrix is The numerical values of the state vectors corresponding to the three modes of the stochastic model and the transition rate matrix are presented in [18]. The values of the transition rate matrix are representative of a slow switching system, with the states spending most of the time at the equilibrium operating points of the different stochastic modes and barely any time in-between.     The solutions for the first moments of the dynamic and algebraic states are represented in Fig. (5) and (6). The Monte Carlo results converge pretty well to the MJLS solutions. However, small spikes are noticeable, which will be greatly exacerbated in the second moments plots, as it can be seen in Fig. (7) and (8). This phenomena is due to the way stochastic jumps are characterized in the model (6). In the Monte Carlo simulation, this translates to unbounded transient overshoots during switching between two stochastic modes. The jump term in the model represented by (6) needs to be refined to accurately represent transient events when the power system switches from one equilibrium point to another one.
In term of computational performance, the computer execution time was 25 hr 44 min 22 s for 20,000 Monte Carlo runs compared to 1 hr 51 min 41 s for the MJLS model on a PC with a 3.2 GHz Intel R Core TM i7-800 CPU processor with 32GB memory in the MATLAB's environment.

V. CONCLUSION AND FUTURE WORK
A method based on the MJLS was applied to a microgrid operating in grid-tied and standalone operating modes. The microgrid is represented by the modified IEEE-37 bus power system, that includes seven inverters at selected bus locations. The ODEs representing the time evolution of the conditional moments were solved to derive the statistics of the system (zeroth, first and second moments). These solutions were compared to the results of the averaged Monte Carlo simulations. The study shows that the model (6) applies well to a switching system and that the results remain valid for a larger system (56 dynamic states). It was observed that the averaged Monte Carlo simulations of the lower order conditional moments converge pretty well to the solutions of the MJLS. Avenues for future works include improving the model to better characterize the jumps between distinct operating modes and by doing so, eliminate the spikes that occur during the Monte Carlo simulation. Another study will focus on the stability of the system when subject to stochastic switching between distinct operating equilibrium points.