Jump-Diffusion Modeling of a Markov Jump Linear System with Applications in Microgrids

—This paper proposes a jump-diffusion model for the analysis of a microgrid operating in grid-tied and standalone modes. The model framework combines a linearized power system dynamic model with a continuous-time Markov chain (CTMC). The power system has a stochastic input modeled with a one-dimensional Wiener process and multiplicative diffusion coefﬁcient. The CTMC gives rise to a compound Poisson pro- cess that represents jumps between different operating modes. Starting from the resulting stochastic differential equation (SDE), the proposed approach uses Itˆo calculus and Dynkin’s formula to derive a system of ordinary differential equations (ODEs) that describe the evolution of the moments of the system. The approach is validated with the IEEE 37-bus system modiﬁed to form a microgrid. The results match a Monte Carlo simulation but with far lower computational complexity. In addition, the ODEs are amenable to further analysis, such as stability analysis and determination of operational bounds.


I. INTRODUCTION
T HE stochastic jump-diffusion model is well suited to dynamic systems subject to random disturbances of a limited amplitude, as well as to random and abrupt perturbations of a relatively higher magnitude [1]. Disturbances of limited (small) amplitude correspond to the diffusion component of the model and are described using the Wiener process. Higher (large) magnitude disturbances apply to the jump component of the model and are represented by the Poisson process [1], [2], [3]. The objective of this study is to develop such jump-diffusion model for a microgrid that oscillates between standalone (islanded) and grid-tied operating modes. Preceding works focused on deterministic systems controlled by stochastic inputs of the jump type [4], [5] and, in other cases, on stochastic systems subject to switching behaviors of a pure-jump type [6], [7]. This work extends previous models: 1) the entire system is considered switching between distinct stochastic modes, 2) this switching behavior involves all dynamic and algebraic states, and not just the system inputs as in [4], [5], 3) the modeling of the stochastic process includes a Wiener process as well Poisson jump processes in the derivation of the ordinary differential equations that govern the evolution of the system statistics.
The application of stochastic models to the analysis of the dynamics of power systems is relatively new [4], [8]. This representation is particularly relevant for a class of stochastic processes called Stochastic Hybrid Systems (SHSs). In this framework, the linearized power system dynamic model is combined with discrete transitions of the system variables triggered by stochastic events [9], [10], [11], [4], [5], [7]. The linearized model is represented by a system of differential algebraic equations (DAEs) that includes the evolution of the dynamic states, and of the output variables expressed as a linear combination of the dynamic states and the control variables [4], [7], [5]. In a particular class of SHSs called Markov Jump Linear Systems (MJLSs), the stochastic discrete transitions of the system are explicitly modeled as continuous time Markov chains. The jump-diffusion model developed in this study is based on the Markovian property of the jumps, appropriately supplied by a Poisson process, while the diffusion component is supplied by a Wiener process.
For a microgrid, a jump in a dynamic state is defined as a variation of significant magnitude that cannot be represented by a linear drift or a noise perturbation. For instance, when an abrupt change occurs in the loading condition or in the characteristics of distributed power sources, and because of the finite and limited inertia of the system, the bus voltages and the system frequency can experience drastic variations that, in the absence of adequate control, could result in instabilities and even in the collapse of the entire system [12], [7]. In the context of this paper, discrete jumps result in the switching of the microgrid between one grid-tied operating mode and two distinct standalone (islanded) operating modes. The difference between the two standalone modes results from a difference in load conditions and corresponds to two distinct sets of dynamic states and output variables [7], [13]. Figure 1 shows bus voltages and load currents at an arbitrary bus (labeled 26) of the IEEE 37-bus microgrid. This plot was obtained from a Simulink R /PLECS R simulation of the IEEE 37-bus microgrid and from results from an experimental testbed [14].
The derivation of a small-signal model of an inverter-based microgrid system was carried out in [14]. For the IEEE 37bus microgrid, the resulting state space model is of the 225 th order and represent a two-time scale system where states with slow dynamics are mixed with those with fast dynamics [15], [7]. A reduction method based on the singular perturbation technique was presented in [16], [15], [17]. In the case of the 225 th order system (IEEE 37-bus system), the technique allowed to reduce the system to the 56 th order by eliminating the fast states. Accuracy of this method was evaluated in [15], where the system's dynamic response was evaluated in gridtied and islanded modes, as well as during transitions between the two modes. Results obtained using the reduced model were found to match experimental results from a hardware testbed [15].
A SHS model for a power system was developed in [4] to represent the dynamic behavior of a system subject to stochastic inputs. The resulting ODEs representing the evolution of the conditional moments were derived using stochastic inputs and did not include random variations of the Wiener type. This model was used in [7] for a system switching between different equilibrium points. It was shown in [5], [18], [7] that a matrix representation of the system can be used along with analytical tools to compute the system statistics. In this study, the model is adapted for a slow-switching system, and includes a Wiener process component. The derivation of the model is carried out directly from the jump-diffusion stochastic differential equation (SDE), rather than from the SHS framework as in [4], [6], [7]. This paper is organized as follows. In section II, the development of a jump-diffusion model of a power system is presented. In section III, the Dynkin's formula is applied to the Itô formula of the Stochastic Differential Equation to derive the ODEs that represent the evolution of the system statistics. In section IV, a numerical application of the jump-diffusion model to the IEEE-37 bus microgrid system is discussed. A conclusion and avenues for future work are presented in Section V.

II. JUMP-DIFFUSION MODEL OF A POWER SYSTEM
Let X = {X(t), t ≥ 0} be a one-dimensional continuous stochastic process. In the most general sense, a jump-diffusion process can be represented by a stochastic differential equation with jumps dX(t) = f t, X(t) dt + g t, X(t) dW (t) for t ≥ 0 with initial value X(0) = X 0 . Here, W = {W (t), t ≥ 0} represents a standard Wiener process and N = {N (t), t ≥ 0} a Poisson process with intensity λ. The processes f , g, and h are called respectively drift, diffusion, and jump coefficient functions. These coefficient functions could be deterministic, time dependent and non linear, in general. In accordance with power system models, this study considers deterministic, state dependent and linear coefficient functions. Furthermore, the existence and uniqueness of the solution to the SDE (1) are guaranteed only if the drift, diffusion and jump coefficient functions satisfy the Lipschitz conditions, as well as the linear growth conditions, the definition of which can be found in [1], [19], [20]. When given in the differential form, the SDE (1) is considered an informal or short hand notation for its integral form. The integral form yields the solution to (1) and is called an Itô process with jumps or a jump-diffusion process The first integral is an ordinary Riemann integral. The second integral is an Itô integral with respect to the Wiener process W (t). The third term represents a compound Poisson process where, as noted above, {N (t), t > 0} is a Poisson process with intensity λ.
The notation X τ k − corresponds to the almost sure left-hand limit of X(t) at time τ k , and if a jump occurred at time τ k , the jump size is defined as The generalization of (1) and (2) to a multidimensional system is straightforward. The evolution of the stochastic vector X = {X(t) ∈ R n , t ≥ 0}, is described by an ndimensional SDE with jumps [1] and d independent compound Poisson processes, N , with parameters λ , = 1, ..., d, and jump sizes h (t, X(t−)). The integral form is derived from (4) as To develop a jump-diffusion model for a microgrid system, the conventional dynamic model needs to be cast into the stochastic differential equation with jumps defined by (4). The deterministic dynamic model is represented by the differential algebraic equation (DAE) [5], [4] where x(t) ∈ R n is referred to as the dynamic states; y(t) ∈ R v denotes the algebraic states or outputs of the system; and u(t) ∈ R w represents the inputs of the system. The linearization of (6) around a given stable equilibrium point results in a linear affine model [5], [4] x = Ax + Bu + C Following [7], [14], [15], the equilibrium points correspond to the steady-state operating points of a microgrid operating in grid-tied mode (X 1 ) or in two distinct standalone modes (X 2 , X 3 ). The difference between the two standalone modes is based on different load conditions on the system. The affine term, C, is defined in such a way that small variations are modelled around the equilibrium points X 1 , X 2 , X 3 In the output equation, F is set to zero: the outputs depend on the dynamic states and input vector. The first step in developing a stochastic model for a power system is to convert the deterministic variables in equation (7) into stochastic variables To establish a correspondence between (4) and (9), the state equation of the dynamic model is cast as the drift term of the SDE with jumps, i.e., The diffusion term is modeled as a multiplicative diffusion component, with a noise generated ripple proportional to the dynamic state. The diffusion coefficient is a constant adequately chosen to limit the ripple to about 10% of the steady-state value of the dynamic state, where β = 0.1 . The jump size is defined as where X(t) is set to the steady-state operating point corresponding to the destination mode (i.e. mode the system has jumped to), and X(t−) is the value of the dynamic state just prior to the jump event With all the terms of the SDE (4) defined, the jumpdiffusion model of a power system based on the deterministic model (7) can now be represented by the SDE with jumps where W (t) is a one-dimensional Wiener process 1 and N , = 1, ..., d, represents the th Poisson process with intensity λ , and i ∈ {1, 2, 3} represents the system operating mode. 1 A one-dimensional Wiener process is used here to simplify the equations, without loss of generality.

III. DERIVATION OF THE CONDITIONAL MOMENTS
The procedure to derive the conditional moments include the following steps: • Derivation of the Itô formula for the SDE with jumps • Application of the Dynkin's formula to the Itô formula • Derivation of the conditional moments dynamics • Derivation of a matrix representation A. Itô formula for the SDE with jumps The aim of this sub-section is to derive the expression of the stochastic differential equation for a process that is function of the solution of equation (4). In a more formal way, the problem can be stated as follows: Given a stochastic differential equation with jumps of the type (4), and a process ψ(t) which is a function of X(t) where the function ψ(t, X(t)) is continuously differentiable in t and twice continuously differentiable in X, determine the stochastic differential equation for the process ψ(t) Equation (15) is called the Itô formula (or Itô rule) and the definitions of the coefficient functionsf ,g,h are based on the rules of the stochastic calculus. The Itô formula is the equivalent to the chain rule in classical calculus. It is used to quantify the changes in ψ(t, X(t)) that are caused by changes in X(t).

B. Application of Dynkin's formula
In stochastic analysis, Dynkin's formula is a theorem that relates the expectation of a function of a jump-diffusion process and a functional of the backward jump-diffusion operator [3]. It can be interpreted as a stochastic generalization of the fundamental theorem of calculus. For the jump-diffusion SDE of the type (4), the Dynkin's formula consists in taking the expectation of (16) Notice the term with dW (t) has disappeared from (17). In fact, the expected value of the Wiener process E[dW (t)] = 0. In addition, the differential term dN (t) has been replaced by its expected value λ dt where λ represents the parameter of a Poisson process N (t).

C. Derivation of the conditional moments dynamics
The Dynkin's formula, combined with the law of total expectation allows for a derivation of the ODEs that describe the evolution of the conditional moments.
The law of total expectation can simply be defined as follows: given a set of stochastic events A i , i = 1, ...N , the expectation of a random variable X equals the sum of the expectations of X given A i where P r(A i ) denotes the probability of event A i . The conditional moment of a process function ψ(t) given a stochastic event A i is expressed as It follows from (18) that the total expectation of ψ(t, X(t)) equals the sum of all conditional moments given the events and the evolution of a conditional moment: The right hand side of equation (21) represents the Dynkin's formula (17), expressed here with respect to a stochastic mode i, where i = 1, ..., Ṅ For a power system jump-diffusion model (13), jump sizes correspond to jumps in dynamic and algebraic states when the system switches between one equilibrium point to another one. They can subsequently be represented as follows ψ(t, X(t)) − ψ(t, X(t−)) = ψ j (t, X(t)) − ψ i (t, X(t)) (23) where "i" represents the origin mode and "j" the destination mode. Therefore, these poissonian jumps between different modes of the stochastic system can be characterised by a set of jump parameters (or transition rates) λ ij , i, j = 1, ..., N . The full transition rate matrix for a stochastic system with N stochastic modes is defined as where, in accordance with the Markov chain theory, the selftransitions 2 The jump term in (22) can now be written as The function process ψ i (t, X(t)) of a stochastic variable X(t) has not been defined yet. It can be expressed in many different ways. For instance, for ψ i (t, X(t)) = X(t)   i , equation (22) represents the evolution of the conditional moment of the first order (mean) and will be denoted byμ (1) i (t), and for ψ i (t, X(t)) = X(t) 2 i , equation (22) represents the evolution of the conditional moment of the second order (also called uncentered second moment) and will be denoted bẏ µ (2) i (t). In general, conditional moments of order m will be denoted by µ It results, for the 0 th order moment For the k th element of X(t), X k , the first moment is defined as The uncentered second moment can be defined with respect to one element (k th ) of X(t) or two elements, (k th ) and ( th ) It follows from (28), (29), (30), (31) and from the jumpdiffusion model of a power system in (13) A more detailed derivation of the conditional moments of a stochastic hybrid system is described in [4]. However, the resulting model in [4] is different from the one developed in this paper, for the reason that the power system in [4] is subjected to stochastic inputs whereas in this study, the entire system is switching stochastically between distinct equilibrium points.
In (32), ψ j (t, X(t) is set to the equilibrium operating point of the destination mode (slow switching system) after a jump event has occurred. In addition, the second term of the right hand side of (32) is equal to zero for m < 2.
The resulting system of ODEs that describes the evolution of the conditional moment of the jump-diffusion model of a power system iṡ ∀ i = 1, ..., N , where v p represents the p th element of the vector Bu + C and e p , e r are two unit vectors with a 1 at the p th , r th entry, respectively.
1) The 0 th moments represent the occupational probabilities, which are the probabilities for the system to be in a particular operating mode. The matrix representation of the ODEs (33) is described in equation (34). This corresponds to the expression that was obtained in [5], [7], [18] where Λ T is the transpose of the full transition rate matrix (24). Equation (34) corresponds to the well-known Chapman-Kolmogorov equation.
2) The 1 st conditional moments are the statistical means of the stochastic distribution. The general form of the 1 st moment ODEs derived in [5], [7] was as followṡ where the evolution of the first moments is function of the moments of the 1 st and 0 th orders. Using Matlab notation, the coefficient functions G (1) and H (1) in [5], [7] were then defined as Because of the way the jump term is defined in the jumpdiffusion model developed here, the matrices G (1) and H (1) are modified as follows where • represents the Hadamard product (element-wise multiplication, Matlab '.*'), Λ trans = Λ − diag(Λ), and 1 n×1 is an n × 1 matrix of all ones (Matlab ones(n, 1)).
3) The 2 nd conditional moments represent the uncentered second moments. From [5], [7] where the evolution of the second moment is function of the moments of the 2 nd and 1 st orders. The structure of the coefficient functions G (2) and H (2) is pretty complex and was described in length in [5]. However, for the jump-diffusion model, because of the way the jump term is defined and of the inclusion of the Wiener term, the equation of the ODEs representing the evolution of the second moment is as followṡ Matrix G (2) contains the contributions from the state matrices A i , i = 1, 2, 3, from the diagonal elements of the transition rate matrix (λ ii , self-transitions), and from the square of the Wiener coefficient 1 2 β 2 . Each block entry of G (2) is defined as + diag(Λ ⊗ eye(n + n z )) + β 2 * eye(n * (n + 1)/2) where n z is the binomial coefficient of n and 2, T x , T y (and T z below) are reordering matrices. These matrices are not unique, they ensure the elements of G ( 2) are ordered in a desired fashion.
Matrix H (2) is constructed with the affine term in the drift component of the SDE (13), −A i X i + B i u, i = 1, 2, 3. Each bloc entry is Matrix J (2) is a Kronecker product of the transition rates (non-diagonal elements of Λ) with the steady-state vectors. The block elements of J (2) are of the form where i, j = 1, 2, 3 and i = j.

IV. NUMERICAL APPLICATION
In this section, a numerical application of the jump-diffusion model to the IEEE 37-bus microgrid is presented and discussed. The ODEs representing the conditional moments are solved using the matrix representation of the jump-diffusion model. Computation of the conditional moments is limited to the lower order moments (zeroth, first and second orders). These lower order moments correspond respectively to the occupational probabilities, the statistic means and the uncentered second moments, and are in general sufficient to characterise a probability density function of a distribution. Conceivably, moments of higher order can be calculated, but the analytical expression of the system of ODEs becomes too complex and the solutions too computational expensive, they are not addressed in this paper.
The results presented in this study were limited to 20,000 runs and the computer execution time was 55 hr 20 min 23 s, on a PC with a 3.2 GHz Intel R Core TM i7-800 CPU processor with 32GB memory in the MATLAB's environment. The computer execution time increases significantly as the number of runs is increased. In this regard, the jump-diffusion method developed here is vastly superior to the averaged Monte Carlo method. The execution time to compute the zeroth, first and second order conditional moments was 2 hr 18 min 54 s, a reduction of 95.8% relative to the Monte Carlo approach.
1) Zeroth moment results For the zeroth moments, the jump-diffusion model corresponds to the Chapman-Kolmogorov equations, the differential expression of which is provided in (34). The solutions for the zeroth order moments are identical to those obtained in [7], they are represented here again, in Fig. 2 The first moments represented in Fig. (3) correspond to the dynamic states associated with Inverter 2 of the modified IEEE 37-bus power system [21], [18], [14], [15]. The plots in red represent the averaged Monte Carlo simulations. They closely match the jump-diffusion plots in dashed-black. The error observed in Fig. 3 (b) and (h) can be considerably reduced by increasing the number of Monte Carlo runs, but at the expense of the execution time.

3) Second moment results
To obtain the solutions for the uncentered second moments, the moment order is set as m = 2 in matrix ODEs (41). It results from (41), G (2) is 4788×4788, H (2) is 4788×168, and J (2) is a 4788 × 3. Following the argument developed for the first order moments, the uncentered second moments for the dynamic states is a 1596 × 1, and the second moments for the algebraic states is a 1540 × 1 vector. The results are shown in Fig. (4) for the dynamic states (inverter 2). An examination of Fig. (4) indicates the averaged Monte Carlo simulation results converge to the jump-diffusion results. The comments related to the comparison between the first moments obtained from the Jump-Diffusion model and form the MJLS model from [7] are also valid for the second moments and the variances.

V. CONCLUSION AND FUTURE WORK
This study presented a jump-diffusion model to analyze the dynamics of a microgrid operating in grid-tied and standalone modes. For a switching system, a statistical approach based on the computation of the conditional moments is suitable for the analysis of the system dynamics. The SDE representing the system included the Wiener process as well as Poisson processes to represent the stochastic jumps between different operating modes of the power system. The solution of the jump-diffusion SDE led to a system of ODEs that represent the evolution of the conditional moments of the system. A matrix representation of the system of ODEs, and the numerical solutions were shown to converge to the results of the averaged Monte Carlo simulation. Future work will include using the Jump-Diffusion model developed in this study to analyse the stability of a microgrid system subjected to random and abrupt transitions between different operating modes.