Grid-Forming Control of Voltage Source Converters Based on the Virtual-Flux Orientation

The increasing penetration of renewable generation in power systems is causing growing concern about system stability. Grid-forming (GFM) control has been appointed as the technology required for achieving a high penetration of renewables in the grid, as it successfully contributes to the power system stability. This paper proposes a novel grid-forming control scheme for voltage source converters (VSC) achieving the operational characteristics of a virtual synchronous machine (VSM). The proposed scheme is based on the orientation of a defined virtual flux to a reference axis obtained from the emulation of the synchronous generator swing equation. The virtual flux is obtained by the integration of the VSC internal voltage and it is oriented to the reference axis by means of a flux controller that also controls the flux magnitude. In this way, the VSC synchronism is maintained, operating as a voltage source connected to the grid. Moreover, the flux orientation in turn allows to control the VSC active power, while the flux magnitude control allows to regulate the VSC reactive power or terminal voltage in isolated mode. The paper also demonstrates that the stability of GFM converters is negatively affected when emulating high inertia values. So, a stabilizer is also proposed for damping these modes. Finally, due to fact that GFM converters behave as voltage sources, their low voltage ride through capability becomes an issue, because they naturally respond to faults with high short-circuit currents, that cannot be withstand by power electronics converters. Therefore, a method for limiting fault currents is proposed as well in the paper. The proposed grid-forming control scheme has been implemented in a real-time control platform and validated using comprehensive simulation models in a real-time simulator for assessing its grid-forming capability.


I. INTRODUCTION
It is a well-known fact that power generation has changed significantly due to the massive integration of renewable energy sources (RES) in power systems worldwide. Renewable power plants are mainly based on inverter-based resources (IBR) for their grid connection, unlike conventional synchronous generators (SG) that use an electromechanical device for such a purpose.
According to the International Renewable Energy Agency (IRENA) [1] more renewable power capacity was added to the grid annually than all fossil fuels and nuclear combined The associate editor coordinating the review of this manuscript and approving it for publication was Md. Rabiul Islam .
between 2013 a 2020. Renewable power technologies now dominate the global market for new electricity generation capacity. So, it can be concluded that electricity generation growth based on RES in detrimental of non-renewables power plants is a fact nowadays. It seems clear that electrical systems are moving from classical centralized power system where large synchronous generators (SG) units ensure power stability and robustness, towards a decentralized power system where the generation shifts into the distribution level. This type of generation is based on IBR that can be grouped into two categories: a) grid-following converters (GFLs) and b) grid-forming converters (GFCs).
GFL inverter controls are used in most grid-connected IBRs today. These converters typically use a phase-locked loop (PLL) and a fast-current controller. A PLL is used to synchronize the IBR ''following'' the grid voltage by extracting its phase information. This angular position of the grid voltage vector is needed for the internal current control loops. For that reason, GFL-VSC converters are considered as controlled current sources. Due to its control structure, GFL control cannot instantaneously respond to load changes on the grid as it is the case of SGs. There is always an inherent delay in the GFL-VSC response caused by the PLL, so that IBR's active and reactive currents are injected when the disturbance has passed. In addition, it has been shown that GFLs can become unstable under certain low-strength systems [2], [3]. As a consequence of SG declining due to the high penetration of GFL converters in power systems, GFM converters arise as a promising alternative to solve the shortcomings of GFL converters.
GFM controls have the main purpose of maintaining internal voltage phasor in the sub-transient and transient time frame [4]. This allows the IBR to immediately respond to changes in the external network maintaining stability during severe disturbances. Besides, the internal voltage vector must be controlled to maintain synchronism with other generators by regulating active and reactive power appropriately to support the grid. GFM-VSC are recommended to provide some services such as, operation in low system-strength conditions, grid frequency and voltage stabilization, power oscillations damping, re-synchronization capability, system restoration and black-start [5].
On the other hand, as indicated in the literature, [58], [59] there are basically two methods of limiting the current in a GFM-VSC during a voltage dip or severe disturbances in the grid: switching the control mode to a GFL mode during the fault, or limiting the current by using virtual impedances. In the first case, IBR loses all functionalities of the grid-forming control during the fault periods and uses a back-up PLL for the synchronization [60], [61], [62], which has its own stability issues in weak-grids. Current limiting can be achieved directly by saturating the regulators of the inner current loops (if GFM-VSC uses a current regulator), however this causes a reduction of the synchronous stability margin [63] and a wind-up effect in the outer power loop which can lead to instability [64]. To avoid this, some researchers propose the use of virtual resistors, either linear [65] or nonlinear [66]. Also, the influence of the virtual impedance on the grid-forming control for current limitation is analyzed in [67], where it is shown that the current limitation is largely depending on the fault location and the selected virtual impedance. It has also been shown that the use of virtual impedances can lead to stability problems in parallel operation [66]. However, improved virtual impedance methods have been proposed, such as the one reported in [67], which uses a new current limiting method based on the threshold virtual impedance (TVI) during balanced and unbalanced faults.
Therefore, GFM-VSC act as real voltage sources where the internal voltage is kept under control. The module of this internal voltage is determined by an external reactive power control loop, or a voltage regulator, and the voltage angle is controlled by a synchronization loop using the active power. On the one hand, this mode of operation is similar to that of a SG, and has certain advantages such as simplicity. However, it has some disadvantages in regards to its dynamics. In this configuration, no system state variable is directly controlled, which means a lack of robustness. In order to solve this problem, a scheme of several cascaded controllers has been proposed, with a current inner loop and a voltage outer loop [28]. Although it improves the dynamic robustness, the setting of several cascaded controllers has the drawback of slowing down the dynamic response of the system [68].
This paper presents a novel solution using a scheme without any internal current control loops and regulating a state variable, called virtual-flux, that is obtained as a function of the output current and the integral of the voltage measured at converter's terminals. This virtual-flux is proportional to the internal voltage, so that the control of the converter is achieved in a simple and robust way. In the past, virtual-flux was proposed for inverter droop control in microgrids [69], also using direct control in three-phase rectifiers [70], and applying predictive control [71]. Unlike these proposals, the present virtual-flux control is used to perform a GFM control. In addition, in this work a novel active and reactive current limitations are proposed, as well as a power stabilizer in order to avoid power oscillations produced when designing the GFM-VSC with a high inertia constant and a low damping factor. This paper is organized as follows. The VSC system description is presented in Section II. The dynamic electric equations of the power converter are also presented by using virtual-fluxes as state-variables. Section III is aimed to the GFM-VSC control by orientation of the virtual-flux vector to a rotating reference frame. This section also includes the small-signal stability analysis and the flux regulators tuning. Section IV shows a method for obtaining the system control angle using an active power synchronization loop including a power stabilizer (PSS) and an active current limiter. Likewise, in Section V a reactive-power/voltage controller with a reactive current limiter are presented. In Section VI the synchronization process is addressed. Subsequently, in Section VII results of a comprehensive real-time simulation are shown and discussed. The complete GFM-VSC control is executed using a hardware-in-the-loop (HiL) experimental set-up connected to a real-time digital simulator (RTDS) where the response of the VSC is analyzed under different grid disturbances. Finally, in Section VIII the final conclusions of this work are presented.  Fig. 1 shows the schematic of a three-phase VSC connected to a grid using an LC filter, with parameters R f , L f and C f . The internal voltage of the VSC is indicated as e and the grid voltage as v. The filter current is denoted as i. A constant voltage source V dc is assumed at the converter DC bus.

II. SYSTEM DESCRIPTION
A circuit breaker is connected between the VSC terminals and the grid, and its state determines the operation mode of the VSC (on: grid connected, off: islanded). At VSC terminals a local load is also connected. As it will be shown later, the virtual-flux control allows the VSC to be controlled in a robust way regardless of the output filter. Obviously, the type of filter has some effect on dynamics and power quality of the converter, but unlike other GFM controls which use internal current loops, the use of capacitors is not mandatory.
GFM control block shown in Fig. 1 uses the three-phase voltages measured at the VSC terminals (v a , v b , v c ), the output currents in each phase (i a , i b , i c ), as well as the DC voltage V dc to generate the switching pattern S 1...6 that determines the on and off states of the six switches of a twolevel VSC. The VSC can be controlled as a PQ node, setting the active P * and reactive Q * power setpoints at the converter output, or as a PV node. In this case, the reactive power reference is replaced by a voltage setpoint v * .

A. DYNAMIC EQUATIONS
The electrical equations corresponding to each of the three phases of the output filter of the VSC are as follows where the index k corresponds to the phases a, b, c.
The voltages v k can be defined as the derivative of a virtual flux ψ k at the VSC terminals as Substituting (2) in (1), grouping the derivative terms and considering L f as constant, (1) is equal to where the derivative term is defined as the converter virtualflux in each phase ψ kv which is expressed as a function of the output current and the terminal voltage v k as that is, in each phase the converter virtual-flux, ψ kv , is obtained as the sum of the instantaneous current, i k , multiplied by L f and the integral of the voltage, v k , measured at the output terminals of the VSC. The direct measurement of the voltage and current, indicated in Fig. 1, allows the virtual-flux ψ kv to be calculated directly according to (5). The precision in the measurement of ψ kv is not significantly affected by the variation of parameters, since the filter inductance L f is a very stable magnitude. Calculating in (4) the current as a function of the fluxes, it is obtained that By substituting (6) in (3) and taking (2) into account, the system of differential equations that define the dynamics of the VSC in Fig. 1 is obtained as follows where T f is the filter time constant defined as L f /R f . As shown in (7), the state-variables are the virtual-fluxes, ψ kv and ψ k , being e k and v k the control variables and the input variables, respectively. Equation (7) can be also expressed in vector form by transforming the three-phase variables, indicated with subscript k, into space vectors referred to a rotating dq-axes. These axes rotate at ω rad/s, being θ the d-axis angular position respect to a stationary reference frame. A generic spatial vector ⃗ x, in this rotating reference frame, is obtained from the instantaneous three-phase components x a , x b ,x c using the well-known expression where a = e j2π 3 and a 2 = e − j2π 3 . By applying the space-vector definition of (8) in (7), the following VSC dynamics equations are obtained Note that in (9) there are two terms proportional to the angular velocity ω, as a result of the dq-axes rotation. Separating into real and imaginary components the first equation of (9) can be expressed as And the second one as As shown in (10), e d controls the virtual-flux component ψ dv . The dynamic relationship between these two variables responds to a first-order function neglecting the crosscoupling term ωψ qv , and the resistive term ψ d R f /L f , since 1/T f = R f /L f . This term is usually neglected as it is much smaller than the rest of terms. A similar reasoning can be made for the dynamics of ψ qv and its control variable e q in the second equation of (10). Fig. 2 shows the block diagram corresponding to equations (10) and (11) in the Laplace s-domain.

B. STEADY-STATE OPERATION
In this subsection, the operation of the VSC in steady state is analyzed. Considering in (9) that the derivatives are zero and the filter resistance R f ∼ 0 (R f ≪ ωL f ) the internal voltage ⃗ e can be expressed in terms of the virtual-flux as and likewise, for the grid voltage ⃗ v By transforming (4) into a vectorial expression using (8), and multiplying by jω it is obtained that Substituting (12) and (13) in (14), and given that the reactance X f = ωL f , this equation is expressed in terms of voltages as In other words, the VSC has an internal voltage, ⃗ e, equal to sum of the voltage drop in the reactance, jX f ⃗ i, and the terminal voltage, ⃗ v. Fig. 3 shows the equivalent circuit corresponding to (15), which is analogous to the equivalent circuit of a SG represented by a voltage source behind a reactance. While the synchronous reactance of a SG is around 1 p.u, X f in a VSC is much smaller (typically between 0.10-0.20 p.u.) The active and reactive power injected by the VSC into the grid, can be expressed in terms of vectors ⃗ v and ⃗ i as follows where the operator ( †) indicates complex conjugate. Calculating the current vector in (15) as a function of ⃗ v and ⃗ e, and separating into real and imaginary part, the expression for active power, P, is obtained as where δ is the power angle between the voltage vectors, ⃗ e and ⃗ v. Likewise, the reactive power expression is equal to Equations (17) and (18) are the classic expressions of the active and reactive power transmitted by a SG to the grid. According to (15) the internal voltage ⃗ e is proportional to the virtual-flux ⃗ ψ v , and therefore, its control allows to regulate the active and reactive power exchanged by the VSC.   4 shows the VSC vector diagram in steady-state operation. The voltage vector ⃗ vangle is taken as the angular reference. The current vector lags by an angle ϕ, which means that the VSC is delivering active and reactive power to the grid. The internal voltage ⃗ e is calculated according to (15), its angular position corresponds to a positive power angle (δ > 0) and its module, e, is higher than v, (e > v). A similar triangle to that representing the voltage balance in the VSC can be represented, lagging 90 • , using fluxes. In this case, the virtual-flux vector is obtained as ⃗ ψ ψ v is leading an angle δ respect to ⃗ ψ and has a higher module (ψ v > ψ). Both conditions, as in the triangle of voltages, determine that the VSC generates positive active and reactive power, according to (17) and (18).
In Fig. 4 it has been considered that all vectors rotate synchronously with the dq-axes at ω rad/s. The d-axis is aligned to the virtual-flux ⃗ ψ v , and the q-axis to the internal voltage vector ⃗ e.

III. VIRTUAL-FLUX ORIENTATION CONTROL
The VSC control scheme for the virtual-flux orientation is shown in Fig. 5. It is composed of four blocks.
In this section, the virtual-flux measurement (VFM) and its orientation control (VFOC) are presented below. The active power control block (APC) and the reactive power control block (RPC) are explained in sections IV and V, respectively.

A. VIRTUAL-FLUX MEASUREMENT
As indicated in (5) the virtual-flux in each phase, ψ kv , is calculated as the product of the current i k and the filter inductance L f plus the integral of the voltage v k . By applying a Clarke transformation abc/αβ to these three-phase virtualfluxes the αβ-components (ψ αv and ψ βv ) are obtained.
where i α , i β and v α , v β are the currents and voltages αβcomponents. To avoid using derivatives and the offset from the unknown initial conditions, the integral is approached by a first-order filter, whose cutoff frequency is close to a value of 1 Hz [73]. This gives rise to an error in the low-frequency range that is very small when the virtual flux is estimated at the fundamental frequency (50 or 60 Hz). Fig. 6 shows the Bode-plot corresponding to the response of a pure integrator 1/s and the first order function 1/(s+2π). As shown at 50 Hz (314.16 rad/s) the frequency response of both functions is almost the same Fig. 7 shows the practical way to calculate the virtualfluxes ψ αv and ψ βv . Both components are obtained by passing v α and v β through two filters, such as 1/(s + 2π), and adding to each output the signals L f i α and L f i β , respectively.  From the αβ-components of voltages and currents, the active and reactive power can be calculated as and the voltage module as

B. VIRTUAL-FLUX CONTROL
The virtual-flux orientation control scheme is shown in Fig 8. Initially, the stationary components ψ αv and ψ βv are calculated in dq-axes by applying a αβ/dq transformation, where the rotation angle, θ, is obtained from the APC block through a synchronization loop.
Once the virtual-flux components in dq-axes are obtained, a control is performed in order to keep the virtual-flux vector, ⃗ ψ v , oriented to the d-axis of the rotating dq-axes. Fig. 9 shows the position of vectors ⃗ ψ v and ⃗ ψ, which are separated by and angle δ, and the dq-axes rotating to an angular speed ω. The angular position of the d-axis is θ respect to the α-axis. Only when the virtual-flux vector, ⃗ ψ v , is aligned with the d-axis it can be said that the VSC is operating synchronously In order to achieve this alignment, a control of the virtualflux through the internal VSC voltage must be carried out.
As explained in section I.A (Fig. 2), the relationship between the dq-components of both magnitudes responds to a first-order function provided that the cross-coupling terms are compensated.
In Fig. 8 the control loops to regulate the virtual-flux components are shown. On the one hand, in the q-axis the reference ψ * qv = 0, in order to align ⃗ ψ v to the d-axis. This reference is compared to ψ qv in order to obtain e q as the output of the PI regulator and an additional feedforward signal, ωψ dv , which is added to compensate the crosscoupling terms. On the other hand, e d is obtained in a similar way by passing through a PI regulator the error between the virtual-flux module reference ψ * v and ψ dv and compensating the cross-coupling term −ωψ qv . By applying an inverse Park transformation to e d and e q using the angle θ, the voltage references, e a , e b , e c , are obtained. From these three-phase voltages the switching pattern of the VSC, S 1...6 , is calculated by using a well-known pulse with modulation (PWM) algorithm.

C. REGULATOR TUNING
The process for tuning the PI regulators is the same regardless of the axis (d or q) considered. After compensating crosscoupling terms, the transfer function G p (s) between the virtual-flux and the internal voltage in each dq components is equal to Fig. 10 shows the block diagram of the close-loop control system where the reference signal is any of the virtualflux components, generically represented as ψ * jv , being the disturbance signal, as indicated in (10) The transfer function corresponding to the PI regulator is expressed as follows where k p is the PI-gain and T c its time constant. The openloop transfer function is represented by the product of G p (s) and G r (s) so that If the time constant of the regulator is chosen equal to the time constant of the filter T c = T f , the open-loop transfer function of (24) becomes in an integral function being then, the close-loop transfer function G ′ (s) So that, the PI-gain is k p = 1/τ i , being τ i the system bandwidth. k p is expressed in unit of s −1 , and in p.u. as When T c = T f a k p = 1 p.u. means a regulator bandwidth of ω 0 , that for a 50 Hz is equal to 314.16 rad/s.

D. SMALL-SIGNAL STABILITY
In (10) and (11) the dynamic equations of the VSC were presented, using as state-variables the dq-components of vectors ⃗ ψ v and ⃗ ψ. When these four equations are expressed in a matrix form, its state-matrix presents 2 complex-conjugated eigenvalues as shown in Fig. 11 where ω 0 = 314.16 rad/s and the filter time constant, according to the parameters indicated in Appendix A, is T f = 159.2 ms. As shown in Fig. 11, there are two eigenvalues over the imaginary axis corresponding to state-variables ψ d and ψ q whose natural frequency ω n is equal to ω 0 = 314.16 rad/s. Likewise, there are two other eigenvalues corresponding to ψ dv and ψ qv with a negative real part equal to 6.28 rad/s, which is the inverse value of T f , and an imaginary value equal to ±ω 0 . These eigenvalues present a small damping factor ξ = 2%. As will be seen below, the virtual-flux regulator will allow increasing the damping factor of these eigenvalues. When the virtual-flux control is implemented, the VSC output voltage vector ⃗ e is calculated as a function of the virtual-flux ⃗ ψ v , and its reference ⃗ ψ * v , according to the following expression where k p is the PI-gain and k i its integral constant that can be expressed in term of T c as k i = k p T c . Separating (28) into real and imaginary parts, e d and e q are obtained as follows Note that in (29) So that (29) is expressed in terms of x d , x q as By substituting (31) in (10) the equations corresponding to dynamics of ψ dv and ψ qv are now expressed as In (32) the terms corresponding to ψ dv in the first equation and ψ qv in the second one, depend on k p affecting directly to the system damping, as will be indicated below.
Now the VSC dynamics is represented by six statevariables, the virtual-flux components ψ dv , ψ qv , the grid flux components ψ d , ψ q , as well as x d and x q . Fig. 12 shows the eigenvalues loci corresponding to the VSC dynamics applying a virtual-flux orientation control with T c = T f and k p varying from 0 to 2 p.u, being the parameters of the system the same than before (ω 0 = 314.16 rad/s and T f = 159.2 ms). The eigenvalues corresponding to ψ d and ψ q remain constant in the imaginary axis. However, the eigenvalues ψ dv , ψ qv are shifted towards the negative real axis when k p increases. As shown in Fig. 12, the imaginary component of theses eigenvalues is constant and equal to ω 0 . When k p = 1 p.u. the system bandwidth matches the fundamental frequency ω 0 = 314.16 rad/s and the damping factor is ξ = 70.7% (same real and imaginary components). The eigenvalues x d and x q remain practically constant in the negative half-plane of Fig. 12.

IV. ACTIVE POWER SYNCHRONIZATION LOOP
The active power controller (APC) is shown in Fig. 13. This controller establishes the angular position, θ, to which the virtual-flux vector ⃗ ψ v must be oriented. The inputs to APC block are: the active power reference, P * , and the frequency ω 0 , as well as the active power measurement, P, and the grid voltage module, v. For the calculation of θ, the APC try to reproduce the swing equation of a SG by calculating the internal frequency ω ′ from the difference between the active power reference P * and the active power P as follows where J is the virtual inertia moment, that is defined as twice the inertia constant H , (J = 2H ) when both parameters are expressed in seconds. Likewise, D is the damping factor which, in p.u., is inverse to the droop constant R (D = 1 R ). The rated frequency is represented by ω 0 rad/ s. As shown in Fig. 13 the control angle θ is obtained by integration of the frequency ω, which is calculated as the internal frequency ω ′ minus ω 1 plus ω 2 . These frequencies are the outputs of an active current limiting (ACL) block and a power system stabilizer (PSS), respectively. So that, the control angle θ is calculated as follows A. ACTIVE CURRENT LIMITER The ACL block prevents the active current I act from exceeding a certain limit I max act . The active power is calculated according to (20), but it can also be expressed as the product of the voltage v and I act as follows I act is proportional to P over v, so that the active current can be limited by limiting the active power. According to (17) the active power is proportional to sinδ, being δ the power angle that is defined as So, if δ is limited, P is also limited, and the way to do this is by acting on θ. Fig. 14 shows the active current limiting block. Firstly, I act is obtained according to (35), as a value proportional to P divided by the module of the voltage v, whose value is previously saturated at its lower limit to a value k + , greater than zero, in order to avoid indeterminacy when v = 0. When the active current is out of limits, e.g. I act > I max act the PI regulator generates a ω 2 < 0 that is saturated to −k, so that the control angle θ is limited, and therefore, the active current as well. The ACL stops working when I act is inside of limits. Similarly, the ACL acts generating a ω 2 > 0 when I act is below the lower limit −I max act .

B. APC LINEARIZED MODEL
The block diagram presented in Fig. 13 denotes the implementation of an APC, but it does not explain the synchronization process principle. For such purposes, the system is linearized around an equilibrium point. In this case equation (33) can be expressed in the s-domain as where overlined magnitudes denote that they are expressed in p.u. Likewise, the angular increment δ can be calculated as a function of ω ′ and ω g as follows where ω′ and ω g are the increment of the internal frequency and the grid frequency, respectively, both expressed in p.u, and ω 0 the rated frequency in rad/s. The dynamic equations (37) and (38) are represented in Fig. 15 as a block diagram, where it has been considered that the active power increment P is function of δ.
Equation (17) denotes a non-linear relationship between the active power P and the angleδ. Considering increments in this equation, the following linear relationship for the active power and angle is obtained where K s is the synchronization constant whose value is calculated by applying the partial derivative of P with respect to δ (15) at the point of equilibrium (''0'') The (3/2) constant does not appear in this equation, since the active power P is expressed in p.u. Considering the VSC parameters in Appendix A (X f = 0.15 p.u.) and assuming that, at the point of equilibrium, the internal voltage e and the grid voltage v are close to 1, (e 0 = v 0 = 1 p.u.) and the angle δ 0 is small (cosδ 0 ∼ 1), the synchronization constant is equal to K s = 6.67 p.u. This value is significantly higher than the value for conventional SGs, which is close to 1 [74]. This causes that the power transmission is produced at lower active power angles.
When observing the block diagram of Fig. 15 and considering the synchronization constant K s , the transfer function between δ and the reference torque P, which is equivalent to the oscillation equation of a SG.
Reordering (41) yields When equating the terms of (42) with regards to the normalized expression of a second-order function s 2 + 2ξ ω n s + ω 2 n , it is obtained that where ω n is the undamped natural frequency and ξ is the damping coefficient. In Fig. 15, the dynamics of the VSC depend on the changes in the active power reference, P * , and the grid frequency, ω g . The choice of the parameters D and J is not free. The damping constant D indicates the active power variation that the VSC has to produce to maintain the system synchronized when the frequency changes. For a common value of R= 0.05 p.u., the damping constant is D= 20 p.u., which denotes that, if there is a frequency variation of 0.05 p.u. (2.5 Hz in a 50 Hz grid), the VSC would increase its active power by 1 p.u. If, in (42), D is considered constant and J = 0, the response of P under changes of P * , corresponds to a first-order function with a negative pole located on the real axis The smaller D is, the further the pole will be from the origin and the higher the natural frequency ω n will be. Therefore, the system's response will be faster and also more sensitive to frequency changes. Since K s ω 0 is a constant, the coefficient D is the first choice when determining the dynamics of the active power synchronization loop. If the inertia constant J is gradually increased from 0 to greater values, the loci of the system poles is as presented in Fig. 16. Initially, the system is overdamped (ξ > 1), and the poles are negative and real until they meet at −2σ where ξ = 1. From this point, when J increases, the poles will become conjugated complexes following the path of a circumference, as shown in Fig. 16. It can be observed that this circumference has a radius of σ and the center is at (−σ , 0). When the inertia constant increases and the value of D is held constant, the system becomes oscillatory and slows down. Similarly, the transfer function between δ and ω g when the power reference is P * = 0 can also be achieved. In this case the power increment P is equal to and considering that the internal frequency increment can be expressed in term of P as Substituting (46) in (45) it is obtained that and reordering, (47) can be expressed as follows P = − ω 0 K s Js The first term corresponds to the inertial response where the product s ω g is defined as the rate of change of frequency (RoCoF), so that the power injected in steady-state (s → 0) due to the inertial response is equal to Unlike this, the second term corresponds to the droop response that in steady-state is equal to In an application where the VSC presents only inertial response (D = 0), the second term of (48) is equal to zero and considering (39), the dynamic relationship between δ and ω g is expressed as follows which presents two complex-conjugated poles on the imaginary axis, being its response critically stable and oscillatory. In this case the use of a power stabilizer is required in order to avoid a poorly damped oscillatory response.  Fig. 17 shows the PSS block integrated on the active power controller. As explained above, a PSS is required when the swing equation of the APC is designed with a high inertia constant and a low damping gain, as when D=0. In these cases, any disturbance in the grid frequency, or any change in the active power reference, produce a poorly damped active power response that must be compensated by a PSS. If (51) is expressed as a differential equation it is obtained that

C. POWER SYSTEM STABILIZER
The PSS transfer function correspond to a washout filter with gain K w , a zero in the origin and a pole in −1/T w , so that Considering the relationship between the frequency and the angle, and according to (39), equation (53) can be expressed in terms of δ as By adding (52) and (54) the following equation is obtained where the new parameter J ′ and D ′ after integrating the PSS are now and, (39) presents a new damping term According to (51) J ′ depends on T w , but its value is similar to J since usually J is higher than T w . J is in the order of seconds when T w < 0.5 s. So, when T w is chosen, J ′ is defined and the value of D ′ depends directly on the PSS-gain K w since the rests of parameters are constants. Note that when no PSS is connected (T w = 0) equation (57) is the same as (51).

FIGURE 18.
Poles loci of the transfer function δ/ ω g under K w variation. Fig. 18 shows the pole placement loci of (57) considering J = 16 s, K s = 6.67 p.u, T w = 0.1s as a function of K w . As shown in this figure when K w = 0 two complexconjugated poles are located on the imaginary axis and are shifted to the negative half-plane following a circumference when K w increases up to reach a value where the poles are on the negative real axis.
As stated above, it seems essential to implement a power stabilizer to ensure an adequate response of the VSC during disturbances in the grid frequency, or when the active power reference changes. Fig. 19 shows the VSC active power response under a frequency change from 50 Hz to 47.5 Hz with a rate of change of frequency (RoCof) of 1 Hz/s. The VSC has been designed with J = 16 s and D = 20 p.u. (droop constant R = 0.05). As shown in that figure, for this droop value, when the frequency is reduced a 5% (2.5 Hz/50 Hz) the active power increases from 0 up to 1 p.u.
The dynamic response is very different when a PSS is used, without PSS the active power presents a high oscillatory response and exceed its rated value. On the contrary, the use of a PSS allows a smoother response of the active power but with a certain delay. The PSS has been tuning with K w = 1 p.u. and T w = 100 ms   20 shows the VSC active power response for different RoCof when using the PSS's parameters indicated above. The lower the RoCof is, the slower the response. The active power response for a RoCof = 1 Hz/s is the same as in Fig. 19, the other two graphs correspond to the response for a higher and lower RoCof than the previous one. According to (50) a negative frequency variation of 5% and a droop D=20 p.u. gives rise to a positive power increase of 1 p.u The VSC can be designed as a static synchronous compensator (SSC), so that it presents inertial response, providing only active power during frequency variations. In this case D = 0 and according to (49) the active power increment is proportional to minus the product of J and the RoCof. Fig. 21 shows the inertial response of a VSC designed as a SSC for different J values during a frequency change with RoCof = 0.5 Hz/s. This frequency change is equivalent to a s ω g = −0.01 p.u./s, so that according to (49) the inertial active power increment is equal to 0.3 p.u. when J = 30 s. As shown in Fig. 21 the VSC only exchange active power when the frequency varies, otherwise is zero. During the period of frequency variation, the active power reaches the steady state without oscillations due to the use of the PSS in the active power synchronization loop.  Fig. 22 shows the reactive power controller that produces the virtual-flux module reference ψ * v . In this case, the inputs to RPC block are the reactive power reference, Q * , or the voltage reference v * , if the VSC is operated as a PV node. The reactive power Q or the voltage module, v, measurements are also required in order to fulfill the control. According to (18) the reactive power exchanged by the VSC mainly depends on the difference between the module of the internal voltage e and the grid voltage v. As indicated in (12) e is proportional to ψ * v , so the reactive power control is achieved by regulating the virtual-flux module.

V. REACTIVE POWER AND VOLTAGE CONTROLLER
Indeed, the RPC can become in a voltage controller in the case the VSC operates as a PV node. Anyway, the control principle is basically the same, it is based on the well-known droop control where the virtual-flux module reference ψ * v is calculated as the product of droop gain n q and the difference between the reactive power reference, Q * , and its actual value, Q, plus an initial virtual-flux ψ v0 This equation can be interpreted as follows. When the reactive power is higher than its reference (Q < Q * )ψ * v should be reduced in order to decrease the reactive power generated. On the contrary, ψ * v should be increase when (Q > Q * ). As shown in (58), the last term ψ * v is a control signal that remains equal to zero when the reactive current I react is inside The reactive current limiter (RCL) block shown in Fig. 22 calculates ψ * v as a function of Q and v. Fig. 23 shows the logic of the RCL. Initially, I react is obtained using (59), where v is previously limited to a value greater than zero in order to avoid indeterminacy when v = 0. Once I react is calculated it is compared to I max react , so that when I react > I max react a PI regulator calculates a negative ψ * v that reduce the virtual-flux module reference ψ * v , limiting as well I react . If the reactive current is inside limits, e.g. I react < I max react the PI-regulator output is ψ * v = 0. A symmetrical regulator is used for negative values of I react

VI. GRID SYNCHRONIZATION CONTROL
As in SGs, GFM-VSCs currently uses a synchronization method when connected to a grid in order to avoid overcurrent or an unstable operation in the connection instant. During synchronization the main breaker shown in Fig. 1 is open, the VSC tries to generate an instantaneous three-phase voltage system identical to the grid, so that when the breaker closes the connection current is practically zero. Fig. 24 shows the synchronization method based on the virtual-flux orientation control (VFOC). In this case, the synchronization is achieved when the virtual-flux vector ⃗ ψ v is oriented to the grid flux vector ⃗ ψ g generated by the three-phase voltages v ag , v bg , v cg , and besides the module of both vectors are the same. Firstly, the three-phase voltages of the grid are passed through a the first-order filter 1/(s + 2π) in order to obtain the gridfluxes in each phase. From these fluxes, and applying a Park transformation, ψ dg and ψ qg components are obtained. Then, in order to orientate the virtual flux to the grid flux a PI regulator is used to maintain the grid flux q-component equal to zero (ψ qg = 0). When this condition is satisfied, the control angle θ indicates the angular position of vector ⃗ ψ g , being its module is equal to ψ dg . On the other hand, the VFOC block shown in Fig. 24 tries to generate a virtual-flux vector whose references are: the module ψ * v that is equal to ψ dg and the angular position θ, so that when the control objectives are achieved ⃗ ψ g and ⃗ ψ v are equal. This is the synchronization condition that ensures that the three-phase voltage systems, before and after the circuit breaker, are the same. So that, the breaker can be closed without overcurrent.

VII. CASE STUDIES AND RESULTS
This section presents the results of the study cases, which include grid synchronization, response to a grid frequency disturbance and phase jump, response to a grid fault and response to a transition from grid-tied to islanded mode. These study cases have been chosen to assess the grid-forming capabilities of the proposed control system, according to the grid code draft of National Grid ESO, ''Minimum Specification Required for Provision of GB Grid Forming Capability''. Finally, the proposed control system has been tested in a 39-node network and compared with conventional grid-forming algorithm. To carry out these studies, the plant illustrated in Fig. 1 has been modelled in RSCAD, for real-time simulation in the RTDS (Real Time Digital Simulator) platform, shown in Fig. 25. Tests, such as phase jumps or frequency changes, proposed by National Grid ESO, require the use of an ideal AC source in the test bench shown in Fig 1. Moreover, the control system, depicted in Fig. 1, has been implemented in a real-time dSpace platform, also shown in Fig. 25, to perform a Hardware in the Loop (HIL) real-time simulation. The parameters used in the real-time simulation of both, the plant and the control system, are given in the Table 1 and Table 2 of Appendix.

A. SYNCHRONIZATION
Before, performing any test for the assessing of the gridforming capability, it is needed to connect the inverter to the grid. In this subsection the grid synchronization control presented above is used for such purpose. Fig 26 and 27 show the synchronization results. The difference between both figures is that in the first one the inverter uses a LC filter while in the second one, only a L filter is used. However, in both figures the state of the breaker is included. So, when the synchronization control has established a virtual flux equal to the grid flux, the breaker can be closed. Fig. 26 and 27 show the instantaneous voltages at both sides of the breaker in one phase of the system. Also, the voltage across the breaker is shown. It is demonstrated that the virtual flux control allows the synchronization even without using a LC filter, because despite the voltage error due to the PWM, the virtual flux components are nearly constant as they are obtained by integration of such voltage.  This can be observed in the results presented in Fig. 28 and Fig. 29, where the direct and quadrature components of the virtual flux and the grid flux, respectively, are presented during the synchronization. Recall, that the qcomponents are zero, which means that both components have the same phase angle, oriented to the reference axis, while the d-components are one, which means they have the same magnitude, equal to the nominal value.
Finally, Fig. 30 shows the inverter currents during the synchronization

B. FREQUENCY AND PHASE-JUMP RESPONSE
Once the inverter is synchronized, the response to a frequency disturbance and a phase jump is obtained. The top graph of Fig. 31 shows the frequency disturbance, which consist of a frequency rise of 0.5 Hz followed by a frequency drop of 2.5 Hz, both over the nominal frequency of 50 Hz and with a ROCOF of 1 Hz/s. Fig. 31 shows the inverter response, in terms of the active and reactive power (in the middle graph) and the active and reactive current (in the bottom graph). The  responses show that the inverter reduces its active power and current following the frequency rise and increases its active power and current following the frequency drop. Moreover, the figure shows that the reactive power is hardly affected by the frequency disturbance, demonstrating the decoupled control of active and reactive power shown in Fig. 5. Furthermore, according to Fig. 13, the power increment following the frequency disturbance is proportional to the frequency deviation multiplied by the constant D, which was set to 50 in this test. So, the frequency rise of 1%, produces a power drop of 50%, while the frequency drop of 5% produced would produce a power increment of 250%. Nevertheless, the latter is not possible, because the inverter power (or current) cannot be increased overrated. Therefore, in order to avoid this overcurrent, the active current limiter, presented in Section IV-A has to limit current to rated. So, as shown in Fig. 31, when current reach its nominal value, it doesn't increase beyond this value, regardless frequency continues dropping.
On the other hand, Fig. 32 shows a 10 • phase jump in the grid voltage. Fig. 33 shows the inverter active and reactive power and current responses to such event. This test is employed to demonstrate the inertial response capability of the proposed control system, because a phase jump means a frequency pulse. First of all, the inverter maintains the synchronism following the event, responding with a  peak of power reduction to the phase jump. Besides, the reactive power is also affected by the event, but in a lower extend. Moreover, Fig. 34 shows a 60 • phase jump in the grid voltage, while Fig. 35 shows the inverter response in this case. The response is quite similar to the previous case, but in this case the active power peak is higher, as a consequence of the higher disturbance. However, both cases show the capability of the proposed grid-forming control to provide inertial response to the grid phase-jump disturbance. Fig. 34 and Fig. 35 also show the instantaneous current response.

C. VOLTAGE RIDE-THROUGH
In this section, the response to a voltage dip is obtained, first when the voltage dip is provoked by the starting of a highpower motor and second when the voltage dip is provoked by a three-phase fault. In the first case, Fig. 36 shows the inverter terminal voltage and active and reactive current during the motor start-up at not load, followed by the loading of the motor at t = 2 s. It can be observed that the voltage drop starts at t = 1 s and lasts about 600 ms due to the high reactive current demanded by the motor during start-up and the low grid SCR, equal to 1. However, the reactive current supplied by the inverter is limited to a maximum value of 1.15 p.u. using the reactive current limiter presented in Fig. 23. Fig. 37 shows the instantaneous currents supplied by the inverter and grid, while Fig. 38 shows the inverter and grid active and reactive power. Finally, Fig. 39 shows the motor torque and rotational speed during the start-up and also the mechanical torque applied to the motor.
In the second case, a 20% voltage dip produced by a three-phase fault is applied to the inverter. Fig. 40 shows the voltage dip and the inverter response in terms of active and reactive current and instantaneous current as well. This response demonstrates the low voltage ride through capability of the proposed grid-forming control scheme. Also, it can be   observed that reactive current is limited to the maximum set value during the fault, as in the previous case. It has to be noted that if the limiter was not implemented, an overcurrent would have been obtained due to the grid-forming capability of the inverter. However, this overcurrent is not allowed in the inverter and the overcurrent protection would have disconnected the inverter, preventing its LVRT capability.

D. HOT SWAP FROM GRID-CONNECTED TO ISLANDED
The transition from grid-tied to islanded mode is performed in this test. The system is first supplying a dynamic load   demanding 1 MW and 300 kVAr. At t = 3 s, the breaker of Fig. 1 is open, separating the inverter and the load from the grid. Fig. 41 shows the inverter voltage and frequency during the transition. A disturbance can be observed in both magnitudes, but nominal voltage and frequency is recovered after a short time. Fig. 42 shows the inverter instantaneous voltages. Moreover, Fig. 43 shows the inverter and grid active and reactive power response. When the inverter is in islanded mode, it supplies the active and reactive power demanded by the load. Therefore, its active and reactive power increases until the load demand is met. Finally, Fig. 44 shows the inverter and grid instantaneous currents, while Fig. 45 shows the motor electromagnetic torque presents a disturbance due to the voltage disturbance during the transition.

E. FREQUENCY AND FAULT RESPONSES IN THE IEEE 39-NODE SYSTEM
The performance of the proposed control system has been tested in the IEEE 39-node network as shown in Fig. 46, composed by 9 SGs and one GFC connected in bus 31. In order to check the frequency response of the GFC, a 140 MW load increase has been applied to the system on bus 15. Likewise, a fault on bus 12 has been applied in order to evaluate also the GFC fault response. The network parameters have been included in Table 3 of the Appendix.  Regarding the frequency response analysis all SGs have their primary regulation system disabled except G3, which will be the generator that assume the load increase. The GFC has also the primary frequency regulation disabled providing active power only during frequency variations contributing exclusively to the inertial response. Fig. 47 shows the active power response of G3 and GFC during a load step of 140 MW at t = 0.1 s. G3 has a rated power of 1000 MVA and a 2% droop constant so that increases its output from 650 MW to 790 MW with a frequency drop from 50 Hz to 49.9 Hz. During the first instants after the load change GFC increases its output power providing inertial response, supporting the system ROCOF. As G3, the initial GFC power is 650 MW and oscillates in the opposite direction to the frequency in order to reduce the ROCOF.
In order to analyze the GFC fault response during a balanced three-phase short-circuit, a 100 ms fault has been applied in bus 12. Figure 48 shows the active power exchanged by G3 and the GFC during the fault, as well as the voltage magnitude at their terminals buses 32 and 31, respectively. The fault response of the SG and the GFC differs because the during the voltage dip, the GFC current is limited to rated, in order to protect the converter. Therefore, in the GFC, active power is proportional to the voltage during the fault. After the fault is cleared, both, the SG and GFC, reach the pre-fault state in about 150 ms, proving the low voltage ride through capability of the GFC.

F. COMPARISON WITH CONVENTIONAL GRID-FORMING ALGORITHM
Finally, a comparison has been made between the performance of the proposed grid-forming control system and a conventional one using inner current loops for current limitation, as shown in Fig. 49. For comparison, the second test of section C was reproduced employing the conventional grid-forming control system, obtaining the response to a 20% voltage dip provoked by a balanced three-phase fault. In addition, the frequency disturbance shown in Fig. 31 of section B has also been reproduced for the conventional control system. According to the literature, current-limited grid-forming converters can lose system stability under a voltage dip occurs. In addition, depending on the depth and the duration of the dip, the converter may not be able to re-synchronize to the grid once the fault is cleared and the voltage has been recovered. Fig. 50, shows how the same voltage dip as in section C provokes the current saturation and as result the regulator operates in open loop, losing the stability which causes that the converter is not able to re-synchronize to the grid once the fault has been cleared. This is explained by the existence of an internal inertia in the system, which causes that the control angle obtained for regulating the active power is not related to power under current limitation. In opposition, the proposed GFC control scheme is able to re-synchronize to the grid once the fault has been cleared without losing the converter stability.  Fig. 51 shows the response of the conventional control system to the same frequency disturbance test of Fig. 31. During the frequency rise, the system response shows that the inverter absorbs active power increasing its current without reaching the limit. However, during the frequency drop, the inverter injects active power until the current limitation is reached, which causes the regulator to operate in open loop losing again the stability.

VIII. CONCLUSION
This paper has proposed a novel grid-forming control strategy for VSCs based on the virtual flux orientation. The virtual flux is defined as the integral of the VSC internal VOLUME 11, 2023 voltage, obtaining in this was a state variable for the VSC control. Then, the virtual flux is aligned to an axis obtained from the emulation of the synchronous generator swing equation. Through this alignment, two objectives are met, first, controlling the synchronism with the grid and, second, controlling the VSC active power. Moreover, the control of the virtual flux magnitude allows controlling the VSC reactive power.
The dynamic model of the VSC using the virtual flux was first introduced and based on this dynamic model a control scheme is derived, demonstrating its analogy to the synchronous generator. A small-signal stability analysis has also been performed to demonstrate the system stability.
The operation as a virtual synchronous machine presents two drawbacks, inherited from the SG. First, the oscillatory nature of the system with high inertia. Hence, a power system stabilizer is proposed for damping these low frequency oscillations. Second, the voltage source nature of the system leads to high short-circuit currents as a response to grid faults. Here, a current limiter is proposed as the VSC cannot withstand currents overrated.
Finally, the assessment of these capabilities has been performed using real-time hardware in the loop simulation employing comprehensive simulation models. Different tests have been performed for assessing the grid-forming capability of the proposed control scheme, including the response to frequency disturbances and voltage phase jump. Also, the low voltage ride through capability has been demonstrated, where the proposed current limiter is key. In addition, the capability to swap from grid-tied to islanded mode has also been assessed. The proposed control scheme has also been tested in a 39-node system to assess its capabilities in an electrical grid comprising SGs. Finally, a comparison has been made between the proposed control system and a traditional GFM control system.
As conclusion, the grid-forming capability of the VSC employing the proposed control scheme have been demonstrated, proving also its excellent dynamic response.