Improved State-space Modelling for Microgrids Without Virtual Resistances

Paranagamage S. A. Peiris; Graduate; Shaahin Filizadeh; Dharshana Muthumuni

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Improved State-space Modelling for Microgrids Without Virtual Resistances PDF

- ORCID：
Paranagamage S. A. Peiris ¹
✉
- ORCID：
Graduate (Student Member, IEEE)
- ORCID：
Shaahin Filizadeh ¹ (Senior Member, IEEE)
✉
- ORCID：
Dharshana Muthumuni ² (Member, IEEE)
✉

1. the Department of Electrical and Computer Engineering, University of Manitoba, Winnipeg MB, R3T 5V6, Canada； 2. the Manitoba HVDC Research Center, 211 Commerce Drive, Winnipeg, MB ,R3P 1A3, Canada

Updated：2024-03-26

DOI：10.35833/MPCE.2023.000085

OUTLINE

Abstract

Power converters and their interfacing networks are often treated as modular state-space blocks for small-signal stability studies in microgrids; they are interconnected by matching the input and output states of the network and converters. Virtual resistors have been widely used in existing models to generate a voltage for state-space models of the network that require voltage inputs. This paper accurately quantifies the adverse impacts of adding the virtual resistance and proposes an alternative method for network modelling that eliminates the requirement of the virtual resistor when interfacing converters with microgrids. The proposed nonlinear method allows initialization, time-domain simulations of the nonlinear model, and linearization and eigenvalue generation. A numerically linearized small-signal model is used to generate eigenvalues and is compared with the eigenvalues generated using the existing modelling method with virtual resistances. Deficiencies of the existing method and improvements offered by the proposed modelling method are clearly quantified. Electromagnetic transient (EMT) simulations using detailed switching models are used for validation of the proposed modelling method.

Keywords

Low-voltage converter; state-space modelling; dynamic phasor; time-domain simulation; eigenvalue analysis

I. Introduction

RAPID expansion of renewable energy generation systems has introduced new requirements for power system modelling and simulation methods. Conventional power sytems with centralized synchronous generation and distributed loads have been transformed with decentralized generation with a mix of converter-based generation, nonlinear loads, and microgrids. Conventional synchronous machine models have been rigorously studied, and their generic controls and interactions are well understood. Conversely, converter-tied renewable energy resources and loads in microgrids consist of advanced control loops and different dynamic characteristics that cannot be easily generalized. These fundamental changes call for detailed modelling methods for careful control design of converter-intensive microgrids that consist of dynamic devices with a wide range of operational time constants [

1].

The need for improved modelling techniques for modern power systems to overcome the limitations of constant-frequency phasor models is highlighted in [

1]-[4]. Dynamic phasor-based modelling has been introduced as a method to capture the dynamics of high-frequency elements [5], [6]. The majority of these modelling methods are geared towards electromagnetic transient (EMT) simulations, and therefore, do not offer analytical capabilities such as linearization, eigenvalue analysis, and other state variable-based tools. Modelling methods based upon state-variable methods [1], [4] that employ differential-algebraic equations (DAE) in a decoupled domain provide better analytical capabilities, e.g., eigenvalue analysis, which are paramount for controller tuning tasks.

Oscillatory stability is an important concern in microgrids due to the variety of dynamic characteristics present in loads and converters. Their controls must be tuned to meet the expected fast response while maintaining stability. Oscillatory stability and mitigatory control modifications are best tackled through small-signal stability studies involving eigenvalue analysis. For such studies, full-order state-space modelling was initially proposed in [

4], clearly pointing out the requirement of small-signal modelling to study the interactions among converters in microgrids.

The converters and loads in microgrids are often connected using short line segments. Converters are normally interfaced to the grid using inductors or inductor-capacitor-inductor (LCL) filters at their outputs depending on harmonic attenuation requirements. When modelled in the state-space domain, both the line and the converter require a voltage input variable, making the interconnection challenging. A solution is achieved by adding a large virtual resistance at the interconnection point and taking the currents from the lines and converters as inputs to generate a voltage. Initially introduced in [

4], this method has been widely utilized in a number of studies that involve optimization [7]-[9], model order reduction [10], [11], dynamic load studies [12], and small-signal stability studies [13]-[15].

While the virtual resistor connection changes the network structure with minimal impact on modelling accuracy, its impacts on stability have not been widely and quantitatively studied. It has been shown only qualitatively [

14] that the value of the virtual resistance affects the fundamental frequency oscillatory modes in small-signal stability studies. A method has been proposed in [16] to eliminate the virtual resistance by ignoring the LCL filter dynamics and using an impedance-based model. This eliminates any dynamics associated with the LCL filters and is detrimental to small-signal stability studies as it ignores important dynamics of the microgrid. Therefore, critical gaps in the modelling and analysis of converter-intensive microgrid exist. The contributions of this paper are two-fold: ① quantitative, detailed assessment of the impact of virtual resistance on small-signal stability studies of a microgrid using eigenvalue analysis; ② mathematical derivation and validation of a microgrid modelling method without using virtual resistances.

The proposed modelling method develops a state-space representation of the overall microgrid with only voltage inputs being the connection points for converters. The developed model is capable of generating computationally efficient time-domain responses while also providing accurate, numerically-derived small-signal models for eigenvalue analysis and controller tuning studies.

The rest of the paper is organized as follows. In Section II, a microgrid is modeled using virtual resistances at nodes where the voltage is required as an input to the state-space equations. The eigenvalues of this system are then analyzed with their impact on the stiffness of the system model. Section III proposes a modelling method without virtual resistances. The improvements of the proposed modelling method are then fully compared with the virtual resistance method. Section IV presents comparison and model validation through time-domain EMT simulations and eigenvalue analysis for validating the proposed modelling method. Finally, Section V concludes this paper.

II. Modelling of Microgrid Using Virtual Resistances

The conventional modelling method of microgrids with virtual resistances is exemplified using a test system consisting of four converters and two loads, as shown in Fig. 1, where $i_{k}$ is the current from k^th inverter. More model parameters used in the simulations in this paper are given in the following text and in Appendix A.

Fig. 1 Schematic diagram of a microgrid with virtual resistances.

A. Dynamic Phasor Formulation

The short line segments connecting the buses in Fig. 1 are modelled in the time domain as given in (1) for a line connecting buses $m$ and $n$ .

{\dot{i}}_{m n} (t) = \frac{1}{L_{m n}} (v_{b u s m} (t) - v_{b u s n} (t) - i_{m n} (t) R_{m n})

(1)

where $v_{b u s m} (t)$ , $v_{b u s n} (t)$ , and $i_{m n} (t)$ are the three-phase bus voltages and line currents in the form of $A (t) s i n (α (t) + ϕ)$ , respectively. $\dot{α} (t) = ω (t)$ is the frequency and may vary with time, $ϕ$ is the initial phase of the waveform, $A (t)$ is the time-varying amplitude of the waveform; and $L_{m n}$ and $R_{m n}$ are the line inductance and resistance between the nodes m and n, respectively. The dot above a variable denotes the time derivative of the variable.

The dynamic phasor form of (1) can be written as (2) and (3) shows the general form of a dynamic phasor quantity. The time domain quantities are denoted in the dynamic phasor form as upper case variables with ~.

{\dot{\tilde{I}}}_{m n} (t) = \frac{1}{L_{m n}} ({\tilde{V}}_{b u s m} (t) - {\tilde{V}}_{b u s n} (t) - {\tilde{I}}_{m n} (t) R_{m n} - j ω (t) L_{m n} {\tilde{I}}_{m n} (t))

(2)

\tilde{X} (t) = X^{q} (t) + j X^{d} (t)

(3)

where superscripts d and q denote the d and q components of the original dynamic phasor, respectively. The formulation of (2) is a resultant of the dynamic phasor differentiation property [

17] of the phasor equivalent of the time-varying signal

x (t)

given by (4a):

P (\frac{d}{d t} x (t)) = \frac{d}{d t} \tilde{X} (t) + j ω (t) \tilde{X} (t)

(4a)

P (x (t)) = \tilde{X} (t)

(4b)

where $P (\cdot)$ is defined as the phasor operator [

17], which transforms the time-varying signal

x (t)

into a time-varying phasor (i.e., dynamic phasor)

\tilde{X} (t)

, as shown by (4b).

Expanding (2) results in two time-domain equations as:

\{\begin{array}{l} {\dot{I}}_{m n}^{d} (t) = \frac{1}{L_{m n}} (V_{b u s m}^{d} (t) - V_{b u s n}^{d} (t) - I_{m n}^{d} (t) R_{m n} - I_{m n}^{q} (t) ω (t) L_{m n}) \\ {\dot{I}}_{m n}^{q} (t) = \frac{1}{L_{m n}} (V_{b u s m}^{q} (t) - V_{b u s n}^{q} (t) - I_{m n}^{q} (t) R_{m n} + I_{m n}^{d} (t) ω (t) L_{m n}) \end{array}

(5)

The model for a series RL load is identical to this line model with only one voltage (i.e., terminal voltage) term. To construct the state-space model of the system for eigenvalue analysis, linearized state-space models of the converters and lines need to be combined. The linearized model for a converter may be expressed in the form of (6), which will be presented in the next subsection.

B. State-space Formulation

The converter models are developed as voltage-dependent current injection components. The inputs of the state-space model are voltages of bus $n$ given by ${Δ V_{b u s n}^{q}, Δ V_{b u s n}^{d}}$ and references given to the converter controller; the output currents, which are injected to the bus, are state variables given by ${Δ I_{n}^{q}, Δ I_{n}^{d}}$ . Since the current injections from the converter are already states of the converter model, the coefficients of the output matrix C are 1, and as a result, ${Δ I_{n}^{q}, Δ I_{n}^{d}}$ , which are the decoupled components of the output current dynamic phasor, become the injections to the bus. The rest of the state variables that include the converter control states and filter states are indicated by ${Δ x_{1}, Δ x_{2}, . . ., Δ x_{n}}$ . $Δ P_{r e f n}$ , $Δ Q_{r e f n}$ , $Δ V_{r e f n}$ , and $Δ ω_{r e f n}$ denote the external setpoints for the converter. The matrices $A_{c o n v n}$ , $B_{c o n v n}^{1}$ , and $B_{c o n v n}^{2}$ contain the elements of the linearized state-space equations given in Appendix B, and are not shown here due to space limitations.

\{\begin{array}{l} [\begin{matrix} Δ {\dot{I}}_{n}^{q} \\ Δ {\dot{I}}_{n}^{d} \\ \begin{array}{l} Δ {\dot{x}}_{1} \\ Δ {\dot{x}}_{2} \end{array} \\ \begin{matrix} ⋮ \\ Δ {\dot{x}}_{n} \end{matrix} \end{matrix}] = A_{c o n v n} [\begin{matrix} Δ I_{n}^{q} \\ Δ I_{n}^{d} \\ \begin{array}{l} Δ x_{1} \\ Δ x_{2} \end{array} \\ \begin{matrix} ⋮ \\ Δ x_{n} \end{matrix} \end{matrix}] + B_{c o n v n}^{1} [\begin{matrix} Δ V_{b u s n}^{q} \\ Δ V_{b u s n}^{d} \end{matrix}] + B_{c o n v n}^{2} [\begin{matrix} Δ P_{r e f n} \\ Δ Q_{r e f n} \\ Δ V_{r e f n} \\ Δ ω_{r e f n} \end{matrix}] \\ [\begin{matrix} Δ I_{n}^{q} \\ Δ I_{n}^{d} \end{matrix}] = C [\begin{matrix} Δ I_{n}^{q} \\ Δ I_{n}^{d} \\ \begin{array}{l} Δ x_{1} \\ Δ x_{2} \end{array} \\ \begin{matrix} ⋮ \\ Δ x_{n} \end{matrix} \end{matrix}] \end{array}

(6)

The linearized state-space model of a line described in (5) is given by:

[\begin{matrix} Δ {\dot{I}}_{m n}^{q} \\ Δ {\dot{I}}_{m n}^{d} \end{matrix}] = A_{l i n e} [\begin{matrix} Δ I_{m n}^{q} \\ Δ I_{m n}^{d} \end{matrix}] + B_{l i n e}^{1} [\begin{matrix} Δ V_{b u s m}^{q} \\ Δ V_{b u s m}^{d} \\ Δ V_{b u s n}^{q} \\ Δ V_{b u s n}^{d} \end{matrix}] + B_{l i n e}^{2} Δ ω

(7)

where $A_{l i n e}$ , $B_{l i n e}^{1}$ , and $B_{l i n e}^{2}$ contain the matrices that have linearized elements from (5); and $Δ ω$ is the linearized term of the frequency term $ω (t)$ from (5). It is noticed that the line and converter models contain the d and q components of the bus voltage as inputs. Since these are not inputs of the actual system, the bus voltages need to be either defined as states or eliminated through substitution. This poses a significant challenge when combining the state-space models for eigenvalue analysis and dynamic simulation studies. A solution is provided by adding a sufficiently large fictitious resistance $R_{N}$ at each bus. The resultant model can be described by:

v_{b u s n} (t) = R_{N} I_{s u m}^{}

(8a)

where $I_{s u m}^{}$ is the sum of currents excluding $R_{N}$ entering or leaving bus n.

Note that (8a) produces different outcome depending on what elements are connected to bus $n$ . For example, if bus $n$ only connects two lines from buses $n + 1$ and $n - 1$ , (8a) becomes:

v_{b u s n} (t) = R_{N} (i_{n, n - 1} (t) - i_{n, n + 1} (t))

(8b)

If bus $n$ also contains a load and a converter, (8a) becomes:

v_{b u s n} (t) = R_{N} (i_{n, n - 1} (t) + i_{k} (t) - i_{l o a d n} (t) - i_{n, n + 1} (t))

(8c)

where $i_{n, n - 1} (t)$ and $i_{n, n + 1} (t)$ are the currents entering and leaving the bus n, respectively; and $i_{l o a d n} (t)$ and $i_{k} (t)$ are the currents to the n^thload and and current from the $k^{t h}$ converter, respectively.

The value of $R_{N}$ is often selected to be in the range of 0.1-10 k $Ω$ [

4], [14] to minimize its impact on the current flow in the network and hence on the accuracy of the model. This model takes in current injections from the connected sources and lines as inputs and produces a voltage that is used as an input for the source and line models. This method has been widely used in the development of microgrid models [7]-[15]. Finally, the d and q components of the small-signal form of (8b) are obtained as:

[\begin{matrix} Δ V_{b u s n}^{q} \\ Δ V_{b u s n}^{d} \end{matrix}] = R_{N} [\begin{matrix} Δ I_{n, n - 1}^{q} \\ Δ I_{n, n - 1}^{d} \end{matrix}] - R_{N} [\begin{matrix} Δ I_{n, n + 1}^{q} \\ Δ I_{n, n + 1}^{d} \end{matrix}]

(8d)

where R_N is an diagonal matrix with the values of virtual resistances along its diagonal.

And the small-signal form of (8c), where a load and converters are connected, is given as:

[\begin{matrix} Δ V_{b u s n}^{q} \\ Δ V_{b u s n}^{d} \end{matrix}] = R_{N} [\begin{matrix} Δ I_{n, n - 1}^{q} \\ Δ I_{n, n - 1}^{d} \end{matrix}] + R_{N} [\begin{matrix} Δ I_{k}^{q} \\ Δ I_{k}^{d} \end{matrix}] - R_{N} [\begin{matrix} Δ I_{l o a d n}^{q} \\ Δ I_{l o a d n}^{d} \end{matrix}] - R_{N} [\begin{matrix} Δ I_{n, n + 1}^{q} \\ Δ I_{n, n + 1}^{d} \end{matrix}]

(8e)

Equations (8d) and (8e) can be substituted in (6) or (7) to eliminate the bus voltage components.

C. Reference Frame Transformation

The term $ω (t)$ in (5) denotes the frequency at which the system equations are decoupled. Since there are multiple frequency generation sources in the system, a single value cannot be used in the development of various component models. As such, the abc-domain quantities of the medium-voltage system (i.e., the line and load currents) are decoupled in the reference frame of converter 1 with a rotational speed of $ω_{1} (t)$ , which is taken as the common reference frame; other converters, including their converter-interfaced transformers, are decoupled in their own reference frames rotating at their respective speeds. Decoupled voltage and current quantities from one reference frame are transformed to and from other reference frames using (9) and its inverse.

[\begin{matrix} X_{g r i d}^{q} (t) \\ X_{g r i d}^{d} (t) \end{matrix}] = [\begin{matrix} c o s δ_{n} (t) & - s i n δ_{n} (t) \\ s i n δ_{n} (t) & c o s δ_{n} (t) \end{matrix}] [\begin{matrix} X_{c o n v}^{q} (t) \\ X_{c o n v}^{d} (t) \end{matrix}]

(9)

{\dot{δ}}_{n} (t) = ω_{1} (t) - ω_{n} (t)

(10)

The state variables in the grid-side reference frame (i.e., converter 1) are denoted by $X_{g r i d}^{d} (t)$ and $X_{g r i d}^{q} (t)$ , and those in the reference frame of converter $n$ are by $X_{c o n v}^{d} (t)$ and $X_{c o n v}^{q} (t)$ . These variables are dq components of either the voltage of the connection bus or the output current of converter, which are based on the transformation used. The transformation angle $δ_{n} (t)$ is generated using the frequency of reference frame of converter $n$ $ω_{n} (t)$ and that of the converter 1 $ω_{1} (t)$ . Since converter 1 is already in the grid-side reference frame, and does not have a transformation angle, for a system with $N$ converters, the number of transformation angles will always be $N - 1$ . This is identical to the reference frame transformation method used in [

4] and all subsequent publications.

D. Converter Model

A block diagram of the converter model used is shown in Fig. 2 and its state-space model equations are provided in Appendix B based on the models given in [

4] and [18]. In Fig. 2, v_c is the three-phase (abc) form of the terminal voltage, and v_c_d and v_c_q are its d and q components, respectively, which are generated from the current controller.

Fig. 2 Block diagram of converter model.

Each line and series RL load element consist of two state equations representing the derivatives of the dq components of their currents, amounting to a total of 12 state equations for the considered microgrid. Each converter model consists of 17 state equations. Since the reference frame transformation takes place for three converters (converter 1 is the common reference frame), the reference frame transformation in (10) generates three state equations, giving a total of 83 state equations for the entire considered microgrid.

E. Initialization and Eigenvalue Analysis

The linearized model of the overall system can be obtained either analytically by linearizing the explicit state equations or numerically by perturbing the systems around an operating point. For this system, the linearized model is obtained numerically using (11):

J_{f} = {[\frac{δ f_{i}}{δ x_{j}}]}_{n \times n}

(11)

where J_f is the system matrix; x is the vector of state variables; and $f$ is the nonlinear system function that contains all the state equations of the system; and subscripts i and j are used to denote the i^th function and the j^th state variable used to calculate the partial derivative term of a particular i^th function with respect to a j^th variable for a system with n differential equations and state variables.

In order to obtain the eigenvalues of the system, the model has to be initialized and brought to a stable operating point. In conventional stability studies, this is achieved using a power flow solver. However, for a state-space model in the dq domain, this method is not adequate due to the large number of dynamic network states that need to be initialized [

16]. Two distinct methods can be used to initialize the system. One method is to use a nonlinear optimization solver, which takes advantage of optimization to minimize the derivatives of the system. This method does not guarantee a solution and can provide an undesired operating point. The other method is to utilise an integration method such as Runge Kutta’s 4^th-order (RK4) method, to run a time-domain simulation till a stable operating point is reached. Using such a method allows the time-domain responses to be validated with EMT simulation tools.

For the given system, the time-domain model can be initialized by setting the values of the state variables to a value close to zero (e.g., $1 \times 10^{- 10}$ ) to make sure that the states remain in the domain of attraction at initialization. For much larger systems that contain many converters and synchronous generators, the sources need to be equipped with voltage sources and solved using a power flow program to obtain the angles and voltage magnitudes of the voltage sources for a desired operating point. Once this is done, the proposed network model shall also be initialized with voltage sources and the obtained bus voltages and angles. This is also similar to the method used in many EMT simulators.

Additionally, once a stable operating point is reached using the time-domain simulation, initialization using an optimization solver can be done for small changes in parameters for stability studies related to eigenvalue analysis. Therefore, in this paper, the second method is used to initialize the system, whilst the optimization solver $f s o l v e$ in Python’s SciPy [

19] optimization library is used to generate subsequent stable operating points. The solution procedure used in this paper is given in Fig. 3.

Fig. 3 Solution procedure used in this paper.

To assess the impact of the virtual resistances, the system in Fig. 1 was initialized for different virtual resistance values ranging from 0.1-10 k $Ω$ . The model was then linearized around the operating point, and the eigenvalue plots for different resistance values are shown in Fig. 4.

Fig. 4 Eigenvalue plots for different virtual resistance values.

The zoomed-in region in Fig. 4 contains all the eigenvalues related to the converters. The eigenvalues closest to the origin are associated with the time constants of the control loops of the converters, whilst the eigenvalues with real values around 500-2500 have higher participation in the filter states. The eigenvalues with real values of around 1000 participate in the power measurement filter states.

It is important to notice that the value of the virtual resistance has minimal impact on the eigenvalues associated with control states, whilst it slightly affects the filter states further away from the origin. However, similar to the observations in [

14], a number of fundamental frequency oscillatory modes can be seen with real values in the range of

10^{7}

10^{8}

, which participate in the states associated with the lines of the microgrid. This is fundamentally inaccurate as the system dynamics do not contain such large coefficients in real world. To fully understand the impact of the virtual resistance on these eigenvalues, the real parts of the eigenvalues participating in system states are listed along with the respective virtual resistances in Table I.

Table I Real Parts of Eigenvalues Participating in System States

Virtual resistance $R_{N}$ (kΩ)	Real part of eigenvalues
0.1	-2278082.04, -1677812.14, -896914.26, -296868.81, -48844.00
1.0	-22777530.58, -16774886.50, -8965904.09, -2965962.24, -491239.16
10.0	-227772014.66, -167745626.10, -89655794.09, -29656781.17, -4911614.37

It can be readily noticed that the real part of these eigenvalues is significantly impacted by the virtual resistance value. As the magnitude of the virtual resistance is increased by a factor of 10, the real part of every eigenvalue is also increased 10 times.

F. Analysis of Impact of Virtual Resistance

The impact of the virtual resistances may be explained further using the simple RL circuit with a large virtual resistance shown in Fig. 5.

Fig. 5 Simple RL circuit with a large virtual resistance.

The time-domain equations for this circuit are given in (12). The notation for a function of time is omitted for simplicity.

\{\begin{array}{l} {\dot{i}}_{1} = \frac{1}{L_{1}} (v_{1} - v_{n} - i_{1} R_{1}) \\ {\dot{i}}_{2} = \frac{1}{L_{2}} (v_{2} - v_{n} - i_{2} R_{2}) \end{array}

(12)

where the variables shown in (12) correspond to Fig. 5. The voltage across the virtual resistance $R_{N}$ is given by:

v_{n} = i_{n} R_{N} = (i_{1} + i_{2}) R_{N}

(13)

Substituting (13) in (12), the state-space model for the simple circuit can be formulated as:

\dot{x} = A x + B u

(14)

\{\begin{array}{l} A = [\begin{matrix} \frac{- (R_{N} + R_{1})}{L_{1}} & \frac{R_{N}}{L_{1}} \\ \frac{R_{N}}{L_{2}} & \frac{- (R_{N} + R_{2})}{L_{2}} \end{matrix}] \\ B = [\begin{matrix} \frac{1}{L_{1}} \\ \frac{1}{L_{2}} \end{matrix}] \end{array}

(15)

\{\begin{array}{l} x = [\begin{matrix} i_{1} \\ i_{2} \end{matrix}] \\ u = [\begin{matrix} v_{1} \\ v_{2} \end{matrix}] \end{array}

(16)

Since (12) is linear, the poles of the system can be obtained by solving $d e t (λ I - A) = 0$ , as given in (17).

\begin{array}{l} λ_{12} = - \underset{C_{1}}{\underset{︸}{(\frac{R_{1} + R_{N}}{L_{1}} + \frac{R_{2} + R_{N}}{L_{2}})}} \pm \\ \sqrt[]{\underset{C_{2}}{\underset{︸}{{(\frac{R_{1} + R_{N}}{L_{1}} - \frac{R_{2} + R_{N}}{L_{2}})}^{2}}} + \underset{C_{3}}{\underset{︸}{\frac{4 R_{N}^{2}}{L_{1} L_{2}}}}} \end{array}

(17)

For the minimal impact on current flows as a result of including the virtual resistance, one must ensure that:

R_{N} > > R_{1}, R_{2}, L_{1}, L_{2}

(18)

Assuming that the line lengths are equal, i.e., $R_{1} = R_{2} = R$ and $L_{1} = L_{2} = L$ , $C_{1}$ , $C_{2}$ , and $C_{3}$ can be written as:

\{\begin{array}{l} C_{1} = \frac{2 (R + R_{N})}{L} \\ C_{2} = 0 \\ C_{3} = \frac{4 R_{N}^{2}}{L^{2}} \end{array}

(19)

Based on (19), the roots of (17) may be approximated as:

\{\begin{array}{l} λ_{1} = - \frac{2 R + 4 R_{N}}{L} \\ λ_{2} = - \frac{2 R}{L} \end{array}

(20)

where $λ_{1}$ is a significantly large number; and $λ_{2}$ is a relatively small number. For example, for the values listed in Appendix A (i.e., $R_{1} = R_{2} = 0.502$ p.u. and $L_{1} = L_{2} = 0.512$ p.u., and with $R_{N} = 84033$ p.u. (=10 $k Ω$ ), the poles will be:

λ_{12} = {- 324782.23, - 0.980}

(21)

If the expressions in (12) are converted to the dq domain with a 60 Hz (377 rad/s) rotating frame similar to (5), the roots shown in (21) will result in complex conjugate pairs of eigenvalues, as shown in (22).

λ_{1 - 4} = {- 324782.23 \pm j 377.67, - 0.980 \pm j 377.67}

(22)

The magnitude of these eigangalues can be related to that of the real components of eigenvalues shown in Table I. The real part of these eigenvalues denotes time constants that are significantly larger than those associated with the system dynamics. For different values of $X / R$ and unequal line lengths, the values of $λ_{1}$ and $λ_{2}$ can be observed to have a similar distinction in magnitude.

Also of note is an important measure to assess the solvability of a system of differential algebraic equations (DAEs) known as stiffness. Stiffness plays an important role when determining the initial values of a set of differential equations. One measure to determine the stiffness of a system is using the stiffness ratio [

20]

R_{s t i}

shown in (23), where

n

is the size of the system and

R e (λ_{i})

is the real part of the

i

^th eigenvalue

λ_{i}

; hence, (23) provides the ratio between the largest and smallest real components of the system eigenvalue. Since the eigenvalues are also associated with the time constants of the system, the stiffness ratio also provides insight into the fastest and slowest dynamics observed in the system.

R_{s t i} = \frac{\underset{1 \leq i \leq n}{m a x} |R e (λ_{i})|}{\underset{1 \leq i \leq n}{m i n} |R e (λ_{i})|}

(23)

Table II provides the stiffness ratios and time-step sizes of the system for three different values of virtual resistance. Along with the maximum time-step, the RK4 solver can be used to initialize the system at a stable operating point.

Table II Stiffness Ratios and Time-step Sizes

$R_{N}$ (kΩ)	Largest eigenvalue	Smallest eigenvalue	Stiffness ratio	Simulation time-step size with RK4 (s)
0.1	$- 2278082 \pm j 377.67$	$0.63 \pm j 18.131$	3587531.00	$10^{- 7}$
1.0	$\begin{array}{l} - 22777530.6 \pm \\ j 377.67 \end{array}$	$0.62 \pm j 18.131$	35870126.90	$10^{- 8}$
10.0	$\begin{array}{l} - 227772014.7 \pm \\ j 377.67 \end{array}$	$0.619 \pm j 18.131$	358696086.08	$10^{- 9}$

It is observed that the simulation time-step size needed to solve the problem is in the order of $10^{- 7}$ - $10^{- 9}$ s. As a result, it is concluded that the virtual resistance has a direct impact on the stiffness of the system, and as its value is increased to obtain more accurate results, the system tends to become increasingly ill-defined, making it difficult to solve for initial conditions and dynamic simulations. These are the core problems that the commonly-used virtual resistance method faces, and the proposed modelling method in the next section will be able to solve them.

III. Proposed Modelling Method Without Virtual Resistances

In order to model the system without virtual resistances, consider a generic system shown in Fig. 6, where a microgrid is connected to a larger high-voltage power system at the topmost node. This terminal can be substituted with a converter (Conv.), load or any other voltage source in order to convert the overall system to a standalone microgrid, as shown in the example provided in Section IV. A procedure to derive a state-space model for this network is provided in this section, which can be extended to a microgrid with $n$ branches and levels.

Fig. 6 Low-voltage network diagram without virtual resistances.

Whilst the introduction of the fictitious virtual resistance provides a voltage input for the currents in the RL elements of the network, the actual power system does not contain such larger resistances at their respective buses. As a result, the requirement for a bus voltage in the microgrid system can be eliminated by substitution and using KCL. Since the model for virtual resistance given in (8) takes in the states of currents in the connected elements as inputs to generate the bus voltage, eliminating the requirement for the bus voltage avoids modelling the line currents that are used as inputs to (8). The procedure to formulate the systems in this form is provided in this section.

To model the system without a virtual resistance, the overall state-space model of the network needs to be formulated to contain only the states of the currents entering the network (in Fig. 6, these are the currents with three numeral subscripts). The voltages at the converter terminals and $v_{g}$ become the inputs to the network. Since voltages are inputs, the voltage across the converter filter capacitance $v_{f}$ is used as the input, whilst all inductances and resistances associated with the LCL filter, the line after the capacitance, and the step-up transformer are included into the microgrid model. Each inducance or resistance between two nodes is described with the same numerical subscripts used for the current travelling across it. The subscript of each voltage describes the number of the paricular node.

Similar to (1), the time derivative of current entering the topmost node can be written as:

{\dot{i}}_{g} (t) = \frac{1}{L_{g}} (v_{0} (t) - v_{g} (t) - i_{g} (t) R_{g})

(24)

where the currents and voltage in (24) correspond to Fig. 6 and (1); and $R_{g}$ and $L_{g}$ are the resistance and inductance between the nodes, respectively. Equation (24) can be written in functional notation as:

{\dot{i}}_{g} = f_{1} (v_{0}, v_{g}, i_{g})

(25)

Given that (24) is linear, it can be rearranged as in (26) and can be written in functional notation as (27).

v_{0} (t) = L_{g} {\dot{i}}_{g} (t) + v_{g} (t) + i_{g} (t) R_{g}

(26)

v_{0} = f_{2} ({\dot{i}}_{g}, v_{g}, i_{g})

(27)

Similar to (24) and (26), an equation can be formulated for node 1 voltage $v_{1}$ as (28) and $v_{0}$ substituted from (26), resulting in (29) with its functional notation form given by (30).

v_{1} (t) = L_{1} {\dot{i}}_{1} (t) + v_{0} (t) + i_{1} (t) R_{1}

(28)

v_{1} (t) = L_{1} {\dot{i}}_{1} (t) + L_{g} {\dot{i}}_{g} (t) + v_{g} (t) + i_{g} (t) R_{g} + i_{1} (t) R_{1}

(29)

\begin{matrix} v_{1} = f_{3} ({\dot{i}}_{1}, i_{1}, \overset{v_{0}}{\overset{⏞}{{\dot{i}}_{g}, v_{g}, i_{g}}}) \end{matrix}

(30)

Similarly, formulating $v_{11}$ and substituting for $v_{1}$ gives (31) and its functional notation by (32):

\begin{array}{l} v_{11} (t) = L_{11} {\dot{i}}_{11} (t) + L_{1} {\dot{i}}_{1} (t) + L_{g} {\dot{i}}_{g} (t) + v_{g} (t) + i_{g} (t) R_{g} + \\ i_{1} (t) R_{1} + i_{11} (t) R_{11} \end{array}

(31)

\begin{matrix} v_{11} = f_{4} ({\dot{i}}_{11}, i_{11}, \overset{v_{1}}{\overset{⏞}{{\dot{i}}_{1}, i_{1}, {\dot{i}}_{g}, i_{g}, v_{g}}}) \end{matrix}

(32)

Since the analytical formulation becomes increasingly complex with the coefficients, the derivations shall be represented in functional notations in the following sections. The notation for function of time is omitted and bus voltages are indicated by their numerical subscripts for simplicity. Substituting (31) in the state equation of the output current of converter 1 $i_{111}$ gives the state equation in (33):

\begin{matrix} {\dot{i}}_{111} = f_{5} (v_{111}, i_{111}, \overset{v_{11}}{\overset{⏞}{{\dot{i}}_{11}, {\dot{i}}_{1}, {\dot{i}}_{g}, i_{11}, i_{1}, i_{g}, v_{g}}}) \end{matrix}

(33)

Similarly, for the output current of converter 2, we can write:

\begin{matrix} {\dot{i}}_{112} = f_{6} (v_{112}, i_{112}, \overset{v_{11}}{\overset{⏞}{{\dot{i}}_{11}, {\dot{i}}_{1}, {\dot{i}}_{g}, i_{11}, i_{1}, i_{g}, v_{g}}}) \end{matrix}

(34)

Now using KCL for bus at 11, the current $i_{11}$ can be written as:

i_{11} = i_{111} + i_{112}

(35)

Taking the derivative of $i_{11}$ gives:

\begin{matrix} {\dot{i}}_{11} = {\dot{i}}_{111} + {\dot{i}}_{112} \end{matrix}

(36)

Substituting (33) and (34) into (36) and rearranging for ${\dot{i}}_{11}$ gives:

{\dot{i}}_{11} = f_{7} (v_{111}, v_{112}, v_{g}, i_{111}, i_{112}, i_{g}, i_{1}, i_{11}, {\dot{i}}_{g}, {\dot{i}}_{1})

(37)

Now substituting for $i_{11}$ from (35) eliminates the $i_{11}$ terms from (38), which results in:

{\dot{i}}_{11} = f_{8} (v_{111}, v_{112}, v_{g}, i_{111}, i_{112}, i_{g}, i_{1}, {\dot{i}}_{g}, {\dot{i}}_{1})

(38)

Similarly, using KCL at bus 12, and formulating ${\dot{i}}_{121}$ and ${\dot{i}}_{122}$ similar to (33) and (34), ${\dot{i}}_{12}$ can be formulated as (39). It is noted that the only remaining derivative terms in the equations are ${\dot{i}}_{1}$ and ${\dot{i}}_{g}$ .

{\dot{i}}_{12} = f_{9} (v_{121}, v_{122}, v_{g}, i_{121}, i_{122}, i_{g}, i_{1}, {\dot{i}}_{g}, {\dot{i}}_{1})

(39)

Using KCL at bus 1 and obtaining derivatives similar to (35) and (36) give:

\begin{matrix} i_{1} = i_{11} + i_{12} \end{matrix}

(40)

\begin{matrix} {\dot{i}}_{1} = {\dot{i}}_{11} + {\dot{i}}_{12} \end{matrix}

(41)

Substituting (38) and (39) to (41) and rearranging the equation for ${\dot{i}}_{1}$ give:

{\dot{i}}_{1} = f_{10} (v_{111}, v_{112}, v_{121}, v_{122}, v_{g}, i_{111}, i_{112}, i_{121}, i_{122}, i_{g}, i_{1}, {\dot{i}}_{g})

(42)

Equation (40) can be expanded by substituting for the $i_{11}$ and $i_{12}$ terms using KCL at buses 11 and 12, resulting in:

\begin{matrix} i_{1} = i_{11} + i_{12} = i_{111} + i_{112} + i_{121} + i_{122} \end{matrix}

(43)

Substituting to (42) from (43) eliminates the remaining $i_{1}$ terms, resulting in:

{\dot{i}}_{1} = f_{11} (v_{111}, v_{112}, v_{121}, v_{122}, v_{g}, i_{111}, i_{112}, i_{121}, i_{122}, i_{g}, {\dot{i}}_{g})

(44)

Using the same procedure from (28) onward for the bus 2 voltage, the current ${\dot{i}}_{2}$ can be formulated as:

{\dot{i}}_{2} = f_{12} (v_{211}, v_{212}, v_{221}, v_{222}, v_{g}, i_{211}, i_{212}, i_{221}, i_{222}, i_{g}, {\dot{i}}_{g})

(45)

Using KCL for bus 0 and obtaining equations similar to (36) and (41), the current ${\dot{i}}_{g}$ can now be written as (46) by substituting for ${\dot{i}}_{1}$ and ${\dot{i}}_{2}$ .

\begin{array}{l} {\dot{i}}_{g} = f_{13} (v_{111}, v_{112}, v_{121}, v_{122}, i_{111}, i_{112}, i_{121}, i_{122}, \\ v_{211}, v_{212}, v_{221}, v_{222}, i_{211}, i_{212}, i_{221}, i_{222}, i_{g}, v_{g}) \end{array}

(46)

The remaining $i_{g}$ can be eliminated by substituting from (47).

\begin{array}{l} i_{g} = i_{1} + i_{2} = i_{11} + i_{12} + i_{21} + i_{22} = \\ i_{111} + i_{112} + i_{121} + i_{122} + i_{211} + i_{212} + i_{221} + i_{222} \end{array}

(47)

It results in a final equation for ${\dot{i}}_{g}$ as:

\begin{array}{l} {\dot{i}}_{g} = f_{14} (v_{111}, v_{112}, v_{121}, v_{122}, i_{111}, i_{112}, i_{121}, i_{122}, \\ v_{211}, v_{212}, v_{221}, v_{222}, i_{211}, i_{212}, i_{221}, i_{222}, v_{g}) \end{array}

(48)

The expression in (48) provides a state equation that only consists of the desired inputs and states, with all the intermediate branch currents and voltages eliminated. This can now be substituted into equations of ${\dot{i}}_{1}$ and ${\dot{i}}_{2}$ as the only derivative term present in (44) and (45) is ${\dot{i}}_{g}$ , thus formulating them in the desired form given in (48). The obtained equations can be then substituted into equations of ${\dot{i}}_{11}$ , ${\dot{i}}_{12}$ , ${\dot{i}}_{21}$ and ${\dot{i}}_{22}$ and finally to the equations of converter output currents ${\dot{i}}_{111}, {\dot{i}}_{112}, \dots, {\dot{i}}_{222}$ to obtain them in the desired form of (48).

As a result of this modelling method, every state equation present in the microgrid model becomes a combination of the states of the currents entering or leaving the microgrid, with the terminal voltages becoming the inputs. The coefficients of these terms are combinations of the resistances and inductances of the lines.

Formulation of this modelling method manually becomes time-consuming and impractical as the network size and complexity grow. However, this procedure can be easily computerized using computational symbolic calculation methods by assigning numerical values to resistance and inductance elements and using symbolic variables for currents and voltages. The code for the generation of the microgrid network for this example is given in [

21], and is generated using Python SymPy [22] computational library. The time-domain equations can then be transformed to the decoupled phasor domain, as shown in (5). The proposed method can be extended to any network that consists of only RL network elements. Any RL loads can be represented with the same line model. Systems can also be partitioned into sections where a parallel capacitive element such as a

π

-section, capacitor bank, or a parallel resistive load is present. A simple two-converter system is shown in Appendix C, which shows that even though it is possible to derive explicit formulations for simple cases, the task becomes rapidly too cumbersome and needs computerized methods for derivation and implementaion.

IV. Comparison and Model Validation

In order to assess the accuracy improvements afforded by the proposed modelling method, the microgrid system shown in Fig. 1 is formulated as per Section III, and the resultant state-space model is given in Appendix D. This model is written in the dynamic phasor domain similar to (2) and (3). Each state equation is then expanded to obtain its d and q components as in (5).

Given that the converter output currents (currents through the transformer) are included in the microgrid, each converter model now has 15 states, with a total of 60 states for the four converters. The microgrid consists of 10 states including 2 states for the load current at bus 3, and 10 states for each of the converter output currents. The load current 1 is substituted to all other currents similar to (48), and hence, is eliminated from the model. The total system model has 73 states including 3 states from the transformation angles. Compared with the system with virtual resistances, the order of the system model is reduced by 10 states. The reduction in order will be more pronounced as the network size grows since the state equations of the interconnecting branches of the microgrid are substituted and eliminated when the system is modelled using the proposed method in this paper.

The developed system model is initialized to the same operating point as the system with virtual resistances was initialized to. The eigenvalue plot of the proposed modelling method and that with virtual resistance value of 10 k $Ω$ is presented in Fig. 7.

Fig. 7 Eigenvalue plot of proposed modelling method and that with virtual resistance of 10 kΩ.

The eigenvalue plot shows an identical match between the eigenvalues present in the two systems (note that the proposed method does not include the 10 erroneous eigenvalues produced by the virtual resistance method), which suggests that the represented dynamics of the system remains intact. Stiffness ratios and the largest simulation time-steps are shown in Table VIII.

Table VIII Stiffness Ratios and Largest Simulation Time-steps

Method	Largest eigenvalue	Smallest eigenvalue	Stiffness ratio	Largest simulation time-step with RK4 (s)	CPU time to simulate 1 s (s)
Proposed	$- 2337.93 \pm j 6409.89$	$0.619 \pm j 18.131$	3812.19	$3 \times 10^{- 4}$	0.7
With $R_{N} = 100$ kΩ	$- 227772014.71 \pm j 377.67$	$0.619 \pm j 18.131$	358696086.08	$1 \times 10^{- 9}$	1642.0

Since the proposed method does not contain the higher-order eigenvalues as in that with virtual resistances, it can be simulated at much larger time-steps (up to $10^{6}$ times higher), allowing it to efficiently run time-domain simulations faster with significantly reduced computational intensity.

To compare the time-domain responses predicted by the proposed method, its dynamic response is compared with an EMT solver. The initialization of the system and instability due to a parameter change are presented to validate the capabilities of the developed model. An detailed model of the microgrid and the converters is developed in PSCAD/EMTDC simulator with the parameters given in Appendix A.

A. System Initialization

The system of four converters is initialized with the loads and power reference values given in Appendix A. For accurate representation, the converters of the PSCAD/EMTDC model are equipped with both average-value voltage source models and detailed switching models with a switching frequency of 3 kHz and an ideal DC voltage source of 1.85 kV. The active power response of converter 1 is shown in Fig. 8 during the system startup.

Fig. 8 Active power response of converter 1 during system startup.

The responses for active power of converter 1 shown in Fig. 8 indicate that the system initialized using RK4 is identical to the average-value voltage source model of the EMT simulation. Slight variations observed in the detailed switching model are due to the switching dynamics of the pulse-width modulator. This shows that the system can be initialized accurately and brought to a stable operating point for linearization without using a power flow program or an optimization solver. Additionally, since the system is modelled without using virtual resistances, its model does not contain the eigenvalues with extremely large magnitudes. Since the eigenvalues in the system represent the actual time constants, the model can be solved for dynamic response studies using larger time-steps that are restricted only by the dynamics of the physical system and not its mathematical model.

B. Instability due to Parameter Variations

In order to investigate small-signal instabilities that can occur due to parameter variations, the gain of the reactive power controller $K_{q}$ is gradually increased from 0.2 to 0.35. The change in eigenvalues for changes in reactive power controller gain is shown in Fig. 9. For simplicity of presentation, only the eigenvalues affected by the parameter are shown.

Fig. 9 Change in eigenvalues for changes in reactive power controller gain.

The eigenvalues that cross the imaginary axis participate in the output currents of converter 3. The increased reactive power droop coefficient attempts to aggressively control the voltage across the filter based on the error in reactive power; as a result, the converters interact with each other, leading to instability. A similar movement of eigenvalues associated with the reactive power controller is shown in [

4]. To validate the prediction of the proposed method with an EMT simulator, the EMT model is initialized with the reactive power gain of

K_{q} = 0.35

The resulting responses are shown in Fig 10. The active power responses indicate a close match between the PSCAD-EMTDC’s average value model and the proposed model; growing oscillations at a frequency of 400 rad/s ( $\approx$ 63 Hz) are predicted as per Fig. 9. The same oscillatory behaviour with the same frequency is observed in the detailed switching model with slight deviations due to the inclusion of switching dynamics.

Fig. 10 Active power responses of converter 1 for $K_{q} = 0.35$ .

C. Modelling for a Meshed Network

The microgrid in Fig. 1 is radial; however, the method can be extended to meshed networks as well. Consider the network in Fig. 11, which is a modified version of the radial network in which the addition of lines between buses 1 and 5 and between buses 2 and 4 has resulted in a meshed topology.

Fig. 11 Schematic diagram of a meshed network without virtual resistances.

The procedure to develop the state-space model for this meshed network is listed as follows.

1) Similar to the procedure for the generalized case in Fig. 6, the equations describing the derivatives of the currents ( $i_{1}$ , $i_{2}$ , $i_{3}$ , $i_{4}$ , and $i_{l o a d 2}$ ) and bus voltages ( $v_{b u s 1}$ , $v_{b u s 2}$ , $v_{b u s 3}$ , $v_{b u s 4}$ , and $v_{b u s 5}$ ) are first formulated as functions of the derivatives of states of $i_{l o a d 1}$ , $i_{42}$ , and $i_{51}$ and the rest of the states.

2) Writing an equation for the derivative of the current $i_{51}$ , and substituting $v_{b u s 5}$ and $v_{b u s 1}$ , the derivative of $i_{51}$ is written in terms of the derivatives of the states of $i_{l o a d 2}$ , $i_{42}$ and the rest of the states.

3) The derivative of $i_{51}$ is now substituted into the expressions of the derivatives of the currents ( $i_{1}$ , $i_{2}$ , $i_{3}$ , $i_{4}$ , and $i_{l o a d 2}$ ) and bus voltages ( $v_{b u s 1}$ , $v_{b u s 2}$ , $v_{b u s 3}$ , $v_{b u s 4}$ , and $v_{b u s 5}$ ).

4) Similar to 2), writing an equation for the derivative of $i_{42}$ , and substituting $v_{b u s 2}$ and $v_{b u s 4}$ , the derivative of $i_{42}$ is written in terms of the derivatives of the states of $i_{l o a d 1}$ and the rest of the states.

5) The derivative of $i_{42}$ is now re-substituted to the derivatives of the currents ( $i_{1}$ , $i_{2}$ , $i_{3}$ , $i_{4}$ , and $i_{l o a d 2}$ ) and bus voltages ( $v_{b u s 1}$ , $v_{b u s 2}$ , $v_{b u s 3}$ , $v_{b u s 4}$ , and $v_{b u s 5}$ ).

6) Using Kirchhoff’s current law (KCL), the derivative of $i_{l o a d 1}$ is formulated to contain only the states of the converter current injections ( $i_{1}$ , $i_{2}$ , $i_{3}$ , and $i_{4}$ ), load current $i_{l o a d 2}$ , and the meshed branch currents $i_{42}$ and $i_{51}$ .

7) The derivative of $i_{l o a d 1}$ is now finally substituted to the state derivatives of the current injections from the converters ( $i_{1}$ , $i_{2}$ , $i_{3}$ , and $i_{4}$ ), load current $i_{l o a d 2}$ , and the branch currents $i_{42}$ and $i_{51}$ to formulate the state-space system.

The state-space model now has state equations for $i_{1}$ , $i_{2}$ , $i_{3}$ , $i_{4}$ , $i_{l o a d 2}$ , $i_{42}$ , and $i_{51}$ . The variables for currents in the intermediate branches can be eliminated using KCL, as described in the derivation in Section III. Compared with the model of the radial network, the model of this meshed network has two extra states for the branch currents of the meshed network. In the case of a radial network, the derivatives of all the intermediary currents can be eliminated using substitution, as each branch current can be represented using the currents entering or leaving the microgrid. If any branch current cannot be solely derived using the currents entering or leaving the microgrid model, it needs to be represented as a state equation due to their dependency on other currents and vice versa, which is the case in a meshed network.

Full formulation of the meshed network is also given in [

21], detailing the steps for obtaining the system of equations to construct the state-space model along with the PSCAD/EMTDC case for validation. The parameters used for existing line 23 and line 34 are used for the new line 42 and line 51, as given in Appendix A.

The time-domain equations describing the meshed microgrid are then transformed into the dynamic phasor form, resulting in 14 state equations for the network model, with the same number of input voltages as in the radial system. Similar to the example in Section IV-A, the meshed network is initialized and validated against an EMT simulator. The active power responses of converter 1 during the startup in the meshed network are shown in Fig. 12.

Fig. 12 Active power responses of converter 1 during startup in meshed network.

V. Conclusion

This paper proposed a novel method to analytically model a microgrid with RL line elements in the state-space domain. Due to the short line segments in microgrids, the RL line elements require a voltage as the input when modelled in the state-space domain. Existing modelling methods use a large virtual resistance to generate this voltage, resulting in an ill-defined system of equations that requires very small time-steps to initialize for linearization and dynamic simulation studies. The new method proposed in this paper models the microgrid state equations by only taking in the converter output voltages as inputs, and its states only consist of inputs and output currents of the network as well as the controller states. The proposed modelling method is validated against EMT simulations by comparing time-domain responses during initialization. Additionally, the numerically linearised model of the system accurately predicts the eigenvalue-based small-signal stability, which is also validated using EMT simulations. Significant time savings in the simulation of the microgrid model have also been observed.

Appendix

Appendix A

Table AI Converter Parameters

Component	Parameter	Value
Converter	Capacity (MVA)	2
	Terminal voltage (kV)	0.69
	Impedance (p.u.)	0.1
	Transformer copper loss (p.u.)	0.021
Filter	$L_{f}$ (H)	0.000335
	$C_{d}$ (F)	0.0007
	$L_{d}$ (H)	0.000621
	$R_{d}$ ( $Ω$ )	1.332
Power	$H$	2
	$K_{q}$	0.2
	$ω_{c}$	1000
	$K_{p}$	1
	$K_{d}$	20
	$ω_{0}$	377
Voltage controller	$K_{p v}$	1
Voltage controller	$K_{i v}$	10
Current controller	$K_{p i}$	5
Current controller	$K_{i i}$	5
Input	$P_{r e f}$ (p.u.)	0.52
	$ω_{r e f}$ (p.u.)	1
	$Q_{r e f}$ (p.u.)	0.3
	$V_{r e f}$ (p.u.)	1

Table AII Microgrid Parameters (p.u. Values are in 34.5 kV, 100 MVA Base)

Parameter	Value (p.u.)	Parameter	Value (p.u.)
Line resistance $R_{m n}$	0.502	Line inductance $L_{m n}$	0.512
Load resistance $R_{l o a d n}$	44.528	Load inductance $L_{l o a d n}$	26.550

Appendix B

The control diagrams related to the equations in this section are based on [

18], and are not shown here due to space limitations. The definitions of variables in this Appendix could also be found in [18] and are not given here. The converter model states are shown as below.

A. Filter States

Given the filter components are represented as RLC elements, they are presented here as single dynamic phasor equations similar to (1). It should be noted that each of these equations generates two state-space equations for the $d$ and $q$ components similar to (2).

\{\begin{array}{l} {\dot{\tilde{v}}}_{f} = \frac{1}{C_{f}} ({\tilde{i}}_{c} - {\tilde{i}}_{S} - {\tilde{i}}_{1 f}) \\ {\dot{\tilde{i}}}_{S} = \frac{1}{L_{g}} ({\tilde{v}}_{f} - {\tilde{v}}_{g} - R_{g} {\tilde{i}}_{S}) \\ {\dot{\tilde{i}}}_{1 f} = \frac{1}{R_{d} C_{f}} ({\tilde{i}}_{c} - {\tilde{i}}_{S} - {\tilde{i}}_{1 f}) - \frac{{\tilde{i}}_{1 f}}{R_{d} C_{d}} + \frac{1}{L_{d}} ({\tilde{v}}_{f} - {\tilde{v}}_{2 c}) \\ {\dot{\tilde{v}}}_{2 c} = \frac{1}{C_{d}} {\tilde{i}}_{1 f} \end{array}

(B1)

B. Power Controller States

\{\begin{array}{l} {\dot{ω}}_{i n t} = \frac{- K_{D}}{2 H} ω_{i n t} + \frac{P_{r e f}}{2 H} - \frac{P_{m e a s}}{2 H} \\ {\dot{P}}_{m e a s} = - ω_{c} P_{m e a s} + ω_{c} (v_{f d} i_{c d} + v_{f q} i_{c q}) \\ {\dot{Q}}_{m e a s} = - ω_{c} Q_{m e a s} + ω_{c} (v_{f d} i_{c q} - v_{f q} i_{c d}) \end{array}

(B2)

The converter frequency is given by:

ω_{n} = K_{p} ω_{i n t} + ω_{r e f}

(B3)

C. Voltage Controller States

\{\begin{array}{l} {\dot{ϕ}}_{q} = v_{r e f} + K_{q} Q_{r e f} - K_{q} Q_{m e a s} - v_{f q} \\ {\dot{ϕ}}_{d} = - v_{f d} \end{array}

(B4)

D. Current Controller States

\{\begin{array}{l} \begin{array}{l} {\dot{γ}}_{q} = (v_{r e f} + K_{q} Q_{r e f} - K_{q} Q_{m e a s} - v_{f q}) K_{p v} + \\ ϕ_{q} K_{i v} + i_{s q} - v_{f d} ω_{0} C_{f} - i_{c q} \end{array} \\ {\dot{γ}}_{d} = - v_{f d} K_{p v} + ϕ_{d} K_{i v} + i_{s d} + v_{f q} ω_{0} C_{f} - i_{c d} \end{array}

(B5)

E. Converter Output Current States

\{\begin{array}{l} {\dot{i}}_{c d} = \frac{ω_{b a s e}}{L_{f}} (- v_{f d} K_{p v} + ϕ_{d} K_{i v} + i_{s d} + v_{f q} ω_{n} C_{f} - i_{c d}) K_{p i} + γ_{d} K_{i i} + \\ v_{f d} + i_{c q} ω_{0} L_{f} - \frac{v_{f d}}{L_{f}} - ω_{i n t} K_{p} i_{c q} - ω_{r e f} i_{c q} \\ {\dot{i}}_{c q} = \frac{ω_{b a s e}}{L_{f}} [(v_{r e f} + K_{q} Q_{r e f} - K_{q} Q_{m e a s} - v_{f q}) K_{p v} + (ϕ_{q} K_{i v} + i_{s q} - \\ v_{f d} ω_{n} C_{f} - i_{c q}) K_{p i}] + γ_{q} K_{i i} + v_{f q} - i_{c d} ω_{0} L_{f} - \frac{v_{f q}}{L_{f}} + \\ ω_{i n t} K_{p} i_{c d} + ω_{r e f} i_{c d} \end{array}

(B6)

The state vector of the converter model is given by:

x = [ω_{i n t}, P, Q, ϕ_{q}, ϕ_{d}, γ_{q}, γ_{d}, i_{c q}, i_{c d}, v_{f q}, v_{f d}, i_{s q}, i_{s d}, i_{1 f q}, i_{1 f d}, v_{2 c q}, v_{2 c d}]

(B7)

Appendix C

The analytical formulation of the model is provided for two converters shown in Fig. C1.

Fig. C1 Analytical formulation of model with two converters.

The inductance and resistance of converter transformer windings are used in the system equations as $R_{2}$ , $L_{2}$ and $R_{1}$ , $L_{1}$ . The current $i_{g}$ through the branch connected to the external grid can be given as:

{\dot{i}}_{g} = \frac{1}{L_{g}} (v_{1} - v_{g} - R_{g} i_{g})

(C1)

Formulating in terms of $v_{1}$ yields:

v_{1} = L_{g} {\dot{i}}_{g} + v_{g} + R_{g} i_{g}

(C2)

Similarly, the derivative of the output current of the converter 1 may be written as:

{\dot{i}}_{1} = \frac{1}{L_{1}} (v_{c 1} - v_{1} - R_{1} i_{1})

(C3)

Substituting for $v_{1}$ from (C2) yields:

{\dot{i}}_{1} = \frac{1}{L_{1}} [v_{c 1} - (L_{g} {\dot{i}}_{g} + v_{g} + R_{g} i_{g}) - R_{1} i_{1}]

(C4)

The voltage $v_{2}$ at bus 2 can be written as:

v_{2} = L_{21} {\dot{i}}_{21} + v_{1} + R_{21} i_{21}

(C5)

Using KCL at bus 1 yields:

\{\begin{array}{l} i_{21} = i_{g} - i_{1} \\ {\dot{i}}_{21} = {\dot{i}}_{g} - {\dot{i}}_{1} \end{array}

(C6)

Substituting for $i_{21}$ , ${\dot{i}}_{21}$ and $v_{1}$ in (C5) yields:

\begin{array}{l} v_{2} = L_{21} \{{\dot{i}}_{g} - \frac{1}{L_{1}} (v_{c 1} - [L_{g} {\dot{i}}_{g} + v_{g} + R_{g} i_{g}) - R_{1} i_{1}]\} + \\ L_{g} {\dot{i}}_{g} + v_{g} + R_{g} i_{g} + R_{21} (i_{g} - i_{1}) \end{array}

(C7)

The derivative for the output current for converter 2 can be written as:

{\dot{i}}_{2} = \frac{1}{L_{2}} (v_{c 2} - v_{2} - R_{2} i_{2})

(C8)

Substituting for $v_{2}$ from (C7) yields:

{\dot{i}}_{2} = \frac{1}{L_{2}} \{v_{c 2} - L_{21} {\dot{i}}_{g} - \frac{L_{21}}{L_{1}} [v_{c 1} - (L_{g} {\dot{i}}_{g} + v_{g} + R_{g} i_{g}) - R_{1} i_{1}] + (L_{g} {\dot{i}}_{g} + v_{g} + R_{g} i_{g}) + R_{21} (i_{g} - i_{1}) - R_{2} i_{2}\}

(C9)

Using KCL at bus 2 and bus 1 yields:

\{\begin{array}{l} i_{g} = i_{1} + i_{2} \\ {\dot{i}}_{g} = {\dot{i}}_{1} + {\dot{i}}_{2} \end{array}

(C10)

{\dot{i}}_{g} = \frac{1}{L_{1}} [v_{c 1} - (L_{g} {\dot{i}}_{g} + v_{g} + R_{g} i_{g}) - R_{1} i_{1}] + \frac{1}{L_{2}} \{v_{c 2} - L_{21} {\dot{i}}_{g} - \frac{L_{21}}{L_{1}} [v_{c 1} - (L_{g} {\dot{i}}_{g} + v_{g} + R_{g} i_{g}) - R_{1} i_{1}] + (L_{g} {\dot{i}}_{g} + v_{g} + R_{g} i_{g}) + R_{21} (i_{g} - i_{1}) - R_{2} i_{2}\}

(C11)

The above equation can now be solved for the derivative of $i_{g}$ , resulting in an equation with converter output currents and terminal voltages as the rest of the variables. This can then be substituted to the state equations of derivatives of $i_{1}$ (C4) and $i_{2}$ (C8) to obtain the state-space model for the two converters. This method can be extended for multiple converters; however, due to the nature of the equations being extensively complicated, a symbolic solver can be easily utilized with the values substituted for the inductances and resistances.

Appendix D

The proposed state-space model for the microgrid shown in Fig. 1 without virtual resistances is derived and formulated in the form shown in (D1).

\dot{x} = A x + B u

(D1)

A = [\begin{matrix} - 239.8 & - 152.8 & - 24.2 & - 133.0 & - 128.6 \\ - 137.8 & - 215.7 & - 3.2 & - 123.9 & - 118.9 \\ - 25.2 & - 26.7 & - 628.3 & - 26.9 & - 25.6 \\ - 100.7 & - 106.5 & 15.0 & - 200.3 & - 123.5 \\ - 92.3 & - 97.7 & 13.8 & - 119.6 & - 208.0 \end{matrix}]

(D2)

B = [\begin{matrix} 46.4 & - 16.3 & - 11.9 & - 10.9 & - 4.3 \\ - 16.3 & 48.4 & - 13.2 & - 12.1 & - 3.5 \\ - 3.0 & - 3.3 & - 3.4 & - 3.1 & - 0.6 \\ - 11.9 & - 13.2 & 48.0 & - 17.0 & - 2.5 \\ - 10.9 & - 12.1 & - 17.0 & 45.4 & - 2.3 \end{matrix}]

(D3)

The state vector x and input vector u are given by:

\{\begin{array}{l} x = [i_{1}, i_{2}, i_{l o a d 2}, i_{3}, i_{4}] \\ u = [v_{c o n v 1}, v_{c o n v 2}, v_{c o n v 3}, v_{c o n v 4}, v_{g}] \end{array}

(D4)

References

U. Markovic, O. Stanojev, P. Aristidou et al., “Understanding small-signal stability of low-inertia systems,” IEEE Transactions on Power Systems, vol. 36, no. 5, pp. 3997-4017, Sept. 2021. [Baidu Scholar]

F. Milano, F. Dorfler, G. Hug et al., “Foundations and challenges of low-inertia systems (invited paper),” in Proceedings of 20th Power Systems Computation Conference, Dublin, Ireland, Jun. 2018, pp. 1-25. [Baidu Scholar]

D. Ramasubramanian, Z. Yu, R. Ayyanar et al., “Converter model for representing converter interfaced generation in large scale grid simulations,” IEEE Transactions on Power Systems, vol. 32, pp. 765-773, Jan. 2017. [Baidu Scholar]

N. Pogaku, M. Prodanović, and T. C. Green, “Modeling, analysis and testing of autonomous operation of an inverter-based microgrid,” IEEE Transactions on Power Electronics, vol. 22, no. 2, pp. 613-625, Mar. 2007. [Baidu Scholar]

K. Mudunkotuwa, S. Filizadeh, and U. Annakkage, “Development of a hybrid simulator by interfacing dynamic phasors with electromagnetic Transactionsient simulation,” IET Generation, Transmission & Distribution, vol. 11, pp. 2991-3001, Aug. 2017. [Baidu Scholar]

J. Rupasinghe, S. Filizadeh, A. M. Gole et al., “Multi-rate co-simulation of power system transients using dynamic phasor and EMT solvers,” The Journal of Engineering, vol. 2020, pp. 854-862, Oct. 2020. [Baidu Scholar]

J. He, X. Wu, X. Wu et al., “Small-signal stability analysis and optimal parameters design of microgrid clusters,” IEEE Access, vol. 7, pp. 36896-36909, Feb. 2019. [Baidu Scholar]

V. N. Kumar and S. K. Parida, “Parameter optimization of universal droop and internal model controller for multi inverter-fed DGs based on accurate small-signal model,” IEEE Access, vol. 7, pp. 101928-101940, Jul. 2019. [Baidu Scholar]

K. Yu, Q. Ai, S. Wang et al., “Analysis and optimization of droop controller for microgrid system based on small-signal dynamic model,” IEEE Transactions on Smart Grid, vol. 7, no. 2, pp. 695-705, Mar. 2016. [Baidu Scholar]

W. Hu, Z. Wu, and V. Dinavahi, “Dynamic analysis and model order reduction of virtual synchronous machine based microgrid,” IEEE Access, vol. 8, pp. 106585-106600, Jun. 2020. [Baidu Scholar]

M. Rasheduzzaman, J. A. Mueller, and J. W. Kimball, “Reduced-order small-signal model of microgrid systems,” IEEE Transactions on Sustainable Energy, vol. 6, no. 4, pp. 1292-1305, Oct. 2015. [Baidu Scholar]

A. Mahdavian, A. A. Ghadimi, and M. Bayat, “Microgrid small-signal stability analysis considering dynamic load model,” IET Renewable Power Generation, vol. 15, no. 13, pp. 2799-2813, May 2021. [Baidu Scholar]

S. Leitner, M. Yazdanian, A. Mehrizi-Sani et al., “Small-signal stability analysis of an inverter-based microgrid with internal model-based controllers,” IEEE Transactions on Smart Grid, vol. 9, no. 5, pp. 5393-5402, Mar. 2018. [Baidu Scholar]

Y. Wang, X. Wang, Z. Chen et al., “Small-signal stability analysis of inverter-fed power systems using component connection method,” IEEE Transactions on Smart Grid, vol. 9, no. 5, pp. 5301-5310, Sept. 2018. [Baidu Scholar]

M. Rasheduzzaman, J. A. Mueller, and J. W. Kimball, “An accurate small-signal model of inverter-dominated islanded microgrids using (dq) reference frame,” IEEE Journal of Emerging and Selected Topics in Power Electronics, vol. 2, no. 4, pp. 1070-1080, Jul. 2014. [Baidu Scholar]

F. Cecati, R. Zhu, M. Liserre et al., “Nonlinear modular state-space modeling of power-electronics-based power systems,” IEEE Transactions on Power Electronics, vol. 37, no. 5, pp. 6102-6115, May 2022. [Baidu Scholar]

V. Venkatasubramanian, H. Schattler, and J. Zaborszky, “Fast time-varying phasor analysis in the balanced three-phase large electric power system,” IEEE Transactions on Automatic Control, vol. 40, no. 11, pp. 1975-1982, Nov. 1995. [Baidu Scholar]

P. S. A. Peiris and S. Filizadeh, “An open-source platform for modelling, simulation, and performance analysis of multi-converter, mixed-generation power systems,” in Proceedings of 18th International Conference on AC and DC Power Transmission (ACDC 2022), virtual event, Jul. 2022, pp. 206-213. [Baidu Scholar]

P. Virtanen, R. Gommers, T. E. Oliphant et al., “SciPy 1.0: fundamental algorithms for scientific computing in Python,” Nature Methods, vol. 17, pp. 261-272, Mar. 2020. [Baidu Scholar]

B. Stott, “Power system dynamic response calculations,” Proceedings of the IEEE, vol. 67, no. 2, pp. 219-241, Feb. 1979. [Baidu Scholar]

S. Peiris. (2023, Feb.). Microgrid model. [Online]. Available: https://github.com/shiroshpeiris/PYSIM.git [Baidu Scholar]

A. Meurer, C. P. Smith, M. Paprocki et al., “SymPy: symbolic computing in python,” PeerJ Computer Science, vol. 3, p. e103, Jan. 2017. [Baidu Scholar]

Address:No.19 Chengxin Avenue, Jiangning District, Nanjing 211106, China

E-mail: mpce@alljournals.cn

Tel:86-25-81093060

Fax:86-25-81093040

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher