A Unified Model for Voltage-Source and Current-Source Virtual Synchronous Generators

1. Introduction

Virtual synchronous generators (VSGs) are widely adopted in renewable generation systems to provide virtual inertia for supporting grid voltage and frequency [1]. Two primary control strategies are derived from the PQ dynamics and synchronous generator dynamic equations, respectively [1,2,3,4,5]. The former employs second-order PQ dynamics to regulate active and reactive power for grid frequency and voltage stabilization; however, it cannot fully replicate the dynamic characteristics of synchronous generators. In contrast, the latter serves as a prominent inertia emulation approach that comprehensively captures flux linkage dynamics, the swing equation, and additional electromechanical transients.

VSGs are classified into voltage-controlled and current-controlled types, which function as voltage sources and current sources, respectively [6]. Voltage-controlled VSGs are better suited to weak grids, whereas current-controlled VSGs perform better in strong grids. Conventional voltage-controlled VSGs are based on the mathematical model of synchronous generators, incorporating stator voltage dynamics, flux linkage equations, and mechanical torque relationships [7,8,9,10].

To further simplify the control strategy, a permanent-magnet synchronous generator (PMSG) model can be adopted for VSG implementation, requiring fewer equations while preserving essential dynamic characteristics [11]. However, under strong grid conditions, voltage-source inverter (VSI)-based VSGs may induce instability. Although current-controlled VSGs have been proposed to mitigate this issue, they remain inherently reliant on VSIs, whose fundamental voltage-output nature imposes limitations on current regulation fidelity [12,13,14]. Moreover, the phase-locked loop (PLL), an essential component of conventional voltage-controlled VSGs, can itself become a source of instability [15,16].

Given the duality between voltage-source inverters (VSIs) and current-source inverters (CSIs) illustrated in Figure 1 [17,18], a CSI-based VSG can be derived as the dual counterpart to the virtual permanent-magnet synchronous generator (VPMSG). Therefore, leveraging the three-phase CSI topology shown in Figure 1(b), this paper proposes the virtual permanent-capacitor synchronous generator (VPCSG) in Section 2. The duality between the VPMSG and the VPCSG is further elaborated. Section 3 establishes a unified mathematical model that encompasses both the VPMSG and the VPCSG. Section 4 derives the stability condition for the unified model in standalone mode via a Lyapunov function candidate, and applies eigenvalue analysis to assess its grid-connected stability. The frequency-locking condition is also derived in this section. Subsequently, Section 5 validates the dynamic stability of both the VPMSG and the VPCSG through time-domain simulations.

2. Derivation of the Virtual Permanent-Capacitor Synchronous Generator Model

Following the space vector formulation of PMSG stator voltage equations in [19,20,21], the VPMSG counterpart is derived analogously as:

$${\overrightarrow{v}}_{\mathrm{s}}=-R_{\mathrm{s}}{\overrightarrow{i}}_{\mathrm{s}}+{\displaystyle \frac{\mathrm{d}{\overrightarrow{\phi }}_{\mathrm{s}}}{\mathrm{d}t}}$$

(1)

Here, the virtual stator voltage, current and flux vectors are defined as

$${\overrightarrow{v}}_{\mathrm{s}}={\displaystyle \frac{2}{3}}(v_a+v_b{\mathrm{e}}^{j{\displaystyle \frac{2\pi }{3}}}+v_{\mathrm{c}}{\mathrm{e}}^{j{\displaystyle \frac{4\pi }{3}}})$$

(2)

$${\overrightarrow{i}}_{\mathrm{s}}={\displaystyle \frac{2}{3}}(i_a+i_b{\mathrm{e}}^{j{\displaystyle \frac{2\pi }{3}}}+i_{\mathrm{c}}{\mathrm{e}}^{j{\displaystyle \frac{4\pi }{3}}})$$

(3)

$${\overrightarrow{\phi }}_{\mathrm{s}}=-L_{\mathrm{s}}{\overrightarrow{i}}_{\mathrm{s}}+{\psi }_{\mathrm{r}}{\mathrm{e}}^{j{\theta }_{\mathrm{e}}}$$

(4)

where v_a, v_b and v_c are the virtual phase voltages; i_a, i_b and i_c are the virtual phase currents; L_s, ψ_r and θ_e are the virtual stator inductance, the virtual flux amplitude and the virtual electrical angle of the rotor, respectively. θ_e=ω_et where ω_e is the virtual electrical angular frequency of the rotor. ω_e=n_pω_m where ω_m and n_p is the virtual mechanical angular frequency and the virtual number of pole pairs, respectively. By the duality principle [22], dual circuits share identical mathematical structures; thus, the dual of Equation (1) is expressed as:

$${\overrightarrow{\widehat{i}}}_{\mathrm{s}}=-{\widehat{G}}_{\mathrm{s}}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}+{\displaystyle \frac{\mathrm{d}{\overrightarrow{\widehat{q}}}_{\mathrm{s}}}{\mathrm{d}t}}$$

(5)

The duality correspondence between Equations (1) and (5) is summarized in Table 1.

Here, the dual counterparts of the virtual stator voltage, current, and flux vectors—namely, the current, voltage, and charge vectors—are defined as follows:

$${\overrightarrow{\widehat{i}}}_{\mathrm{s}}={\displaystyle \frac{2}{3}}({\widehat{i}}_a+{\widehat{i}}_b{\mathrm{e}}^{j{\displaystyle \frac{2\pi }{3}}}+{\widehat{i}}_{\mathrm{c}}{\mathrm{e}}^{j{\displaystyle \frac{4\pi }{3}}})$$

(6)

$${\overrightarrow{\widehat{v}}}_{\mathrm{s}}={\displaystyle \frac{2}{3}}({\widehat{v}}_a+{\widehat{v}}_b{\mathrm{e}}^{j{\displaystyle \frac{2\pi }{3}}}+{\widehat{v}}_{\mathrm{c}}{\mathrm{e}}^{j{\displaystyle \frac{4\pi }{3}}})$$

(7)

$${\overrightarrow{\widehat{q}}}_{\mathrm{s}}=-{\widehat{C}}_{\mathrm{s}}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}+{\widehat{q}}_{\mathrm{r}}{\mathrm{e}}^{j{\widehat{\theta }}_{\mathrm{e}}}$$

(8)

Henceforth, the hat symbol (ˆ) denotes dual quantities. ${\widehat{i}}_a$, ${\widehat{i}}_b$ and ${\widehat{i}}_{\mathrm{c}}$are the dual phase currents; ${\widehat{v}}_a$, ${\widehat{v}}_b$and ${\widehat{v}}_{\mathrm{c}}$are the dual phase voltages; ${\widehat{C}}_{\mathrm{s}}$, ${\widehat{q}}_{\mathrm{r}}$and ${\widehat{\theta }}_{\mathrm{e}}$ are the dual capacitance, the dual charge amplitude and the dual electrical angle, respectively. ${\widehat{\theta }}_{\mathrm{e}}={\widehat{\omega }}_{\mathrm{e}}t$ where ${\widehat{\omega }}_{\mathrm{e}}$ denotes the dual electrical angular frequency. Projecting the abc-frame space vectors onto the dq frame yields the dq-axis current and voltage components, namely,

$${\overrightarrow{\widehat{i}}}_{\mathrm{s}}{}^{\text{'}}={\overrightarrow{\widehat{i}}}_{\mathrm{s}}{\mathrm{e}}^{-j{\widehat{\theta }}_{\mathrm{e}}}={\widehat{i}}_{\mathrm{d}}+j{\widehat{i}}_{\mathrm{q}}$$

(9)

$${\overrightarrow{\widehat{v}}}_{\mathrm{s}}{}^{\text{'}}={\overrightarrow{\widehat{v}}}_{\mathrm{s}}{\mathrm{e}}^{-j{\widehat{\theta }}_{\mathrm{e}}}={\widehat{v}}_{\mathrm{d}}+j{\widehat{v}}_{\mathrm{q}}$$

(10)

Here, ${\widehat{i}}_{\mathrm{d}}$ and ${\widehat{i}}_{\mathrm{q}}$ are the dq-axis components of the dual current vector; ${\widehat{v}}_{\mathrm{d}}$ and ${\widehat{v}}_{\mathrm{q}}$ are the dq-axis components of the dual voltage vector.

Substituting Equation (8) into Equation (5) yields

$${\overrightarrow{\widehat{i}}}_{\mathrm{s}}=-{\widehat{G}}_{\mathrm{s}}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}-{\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\overrightarrow{\widehat{v}}}_s}{\mathrm{d}t}}+j{\widehat{\omega }}_{\mathrm{e}}{\widehat{q}}_{\mathrm{r}}{\mathrm{e}}^{j{\widehat{\theta }}_{\mathrm{e}}}$$

(11)

Multiplying Equation (11) by ${\mathrm{e}}^{-j{\widehat{\theta }}_{\mathrm{e}}}$ yields

$${\overrightarrow{\widehat{i}}}_{\mathrm{s}}{\mathrm{e}}^{-j{\widehat{\theta }}_{\mathrm{e}}}=-{\widehat{G}}_{\mathrm{s}}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}{\mathrm{e}}^{-j{\widehat{\theta }}_{\mathrm{e}}}-{\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\overrightarrow{\widehat{v}}}_s}{\mathrm{d}t}}{\mathrm{e}}^{-j{\widehat{\theta }}_{\mathrm{e}}}+j{\widehat{\omega }}_{\mathrm{e}}{\widehat{q}}_{\mathrm{r}}$$

(12)

Differentiating Equation (10) with respect to t yields

$${\displaystyle \frac{\mathrm{d}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}}{\mathrm{d}t}}{\mathrm{e}}^{-j{\widehat{\theta }}_{\mathrm{e}}}={\displaystyle \frac{\mathrm{d}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}^{\text{'}}}{\mathrm{d}t}}+j{\widehat{\omega }}_{\mathrm{e}}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}^{\text{'}}$$

(13)

Substituting Equation (13) into Equation (12) yields

$${\overrightarrow{\widehat{i}}}_{\mathrm{s}}{}^{\text{'}}=-{\widehat{G}}_{\mathrm{s}}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}{}^{\text{'}}-{\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}^{\text{'}}}{\mathrm{d}t}}-{\widehat{C}}_{\mathrm{s}}j{\widehat{\omega }}_{\mathrm{e}}{\overrightarrow{\widehat{v}}}_{\mathrm{s}}{}^{\text{'}}+j{\widehat{\omega }}_{\mathrm{e}}{\widehat{q}}_{\mathrm{r}}$$

(14)

Equations (9) and (10) are substituted into Equation (14), giving

$$\left\{\begin{array}{l}{\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{d}}}{\mathrm{d}t}}=-{\widehat{i}}_{\mathrm{d}}-{\widehat{G}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{d}}+{\widehat{\omega }}_{\mathrm{e}}{\widehat{C}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{q}}\\ {\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{q}}}{\mathrm{d}t}}=-{\widehat{i}}_{\mathrm{q}}-{\widehat{G}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{q}}-{\widehat{\omega }}_{\mathrm{e}}{\widehat{C}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{d}}+{\widehat{q}}_{\mathrm{r}}{\widehat{\omega }}_{\mathrm{e}}\end{array}\right.$$

(15)

The dual torque equation for the VPMSG is given by

$$J{\displaystyle \frac{\mathrm{d}{\widehat{\omega }}_{\mathrm{e}}}{\mathrm{d}t}}=n_{\mathrm{p}}{T}_{\mathrm{m}}-n_{\mathrm{p}}{\widehat{T}}_{\mathrm{e}}-B_{\mathrm{f}}{\widehat{\omega }}_{\mathrm{e}}$$

(16)

Here, J, ${\widehat{T}}_{\mathrm{e}}$, T_m, and B_f are the virtual total inertia, the virtual electromagnetic torque, the virtual mechanical torque and the virtual viscous friction coefficient, respectively; ${\widehat{\omega }}_{\mathrm{e}}=n_{\mathrm{p}}{\widehat{\omega }}_{\mathrm{m}}$ where ${\widehat{\omega }}_{\mathrm{m}}$is the dual mechanical angular frequency.

Multiplying Equation (8) by ${\mathrm{e}}^{-j{\widehat{\theta }}_{\mathrm{e}}}$ yields

$${\overrightarrow{\widehat{q}}}_{\mathrm{s}}{}^{\text{'}}={\overrightarrow{\widehat{q}}}_{\mathrm{s}}{\mathrm{e}}^{-j{\widehat{\theta }}_{\mathrm{e}}}=-{\widehat{C}}_{\mathrm{s}}({\widehat{v}}_{\mathrm{d}}+j{\widehat{v}}_{\mathrm{q}})+{\widehat{q}}_{\mathrm{r}}={\widehat{q}}_{\mathrm{d}}+j{\widehat{q}}_{\mathrm{q}}$$

(17)

Thus,

$$\left\{\begin{array}{l}{\widehat{q}}_{\mathrm{d}}=-{\widehat{C}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{d}}+{\widehat{q}}_{\mathrm{r}}\\ {\widehat{q}}_{\mathrm{q}}=-{\widehat{C}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{q}}\end{array}\right.$$

(18)

The output electrical power is given by

$$\begin{array}{l}{\widehat{p}}_{\mathrm{s}}={\displaystyle \frac{3}{2}}({\widehat{i}}_{\mathrm{d}}{\widehat{v}}_{\mathrm{d}}+{\widehat{i}}_{\mathrm{q}}{\widehat{v}}_{\mathrm{q}})\\ ={\displaystyle \frac{3}{2}}{\widehat{q}}_{\mathrm{r}}{\widehat{\omega }}_{\mathrm{e}}{\widehat{v}}_{\mathrm{q}}-{\displaystyle \frac{3}{2}}({\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{d}}}{\mathrm{d}t}}{\widehat{v}}_{\mathrm{d}}+{\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{q}}}{\mathrm{d}t}}{\widehat{v}}_{\mathrm{q}})-{\displaystyle \frac{3}{2}}({\widehat{G}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{d}}{}^2+{\widehat{G}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{q}}{}^2)\\ ={\widehat{p}}_{\mathrm{em}}-{\displaystyle \frac{3}{2}}({\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{d}}}{\mathrm{d}t}}{\widehat{v}}_{\mathrm{d}}+{\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{q}}}{\mathrm{d}t}}{\widehat{v}}_{\mathrm{q}})-{\widehat{p}}_{\mathrm{loss}}\end{array}$$

(19)

In Equation (19), the first term represents the virtual input mechanical power, the second term denotes the capacitive power stored in the electric field, and the third term accounts for the power losses. The virtual input mechanical power is defined as

$${\widehat{p}}_{\mathrm{em}}={\widehat{T}}_{\mathrm{e}}{\widehat{\omega }}_{\mathrm{m}}$$

(20)

Thus, the virtual electromagnetic torque is defined as

$${\widehat{T}}_{\mathrm{e}}={\displaystyle \frac{3}{2}}{n}_{\mathrm{p}}({\widehat{q}}_{\mathrm{d}}{\widehat{v}}_{\mathrm{q}}-{\widehat{q}}_{\mathrm{q}}{\widehat{v}}_{\mathrm{d}})={\displaystyle \frac{3}{2}}{n}_{\mathrm{p}}{\widehat{q}}_{\mathrm{r}}{\widehat{v}}_{\mathrm{q}}$$

(21)

Equations (15), (16) and (21) form the VPCSG model, which has a mathematical structure similar to that of the VPMSG. Thus, the VPCSG replicates dominant dynamic characteristics of the VPMSG when the dual operating parameters are set to identical values. Finally, the duality mappings between the VPMSG and VPCSG are summarized in Table 2. Depending on the circuit topology, the swing equations listed in Table 2 are either emulated through virtual inertia control or physically realized via the actual electromechanical dynamics of the circuit. However, implementing the virtual swing equation increases control complexity due to additional state variables and tuning parameters, whereas the physical electromechanical swing equation inherently simplifies controller design by leveraging natural system dynamics. For the topologies shown in Figure 1, the dc-link capacitor equation of the VSI and the dc-link inductor equation of the CSI serve as the respective swing equations. Therefore, explicit implementation of the swing equation in the control strategy becomes unnecessary.

Thus, the duality between VPMSG and VPCSG for the topologies in Figure 1 is summarized in Table 3. Here, i_dc and v_dc denote the input current and capacitor voltage of the VSI, respectively; ${\widehat{i}}_{\mathrm{L}}$ and ${\widehat{v}}_{\mathrm{dc}}$ denote the inductor current and input voltage of the CSI, respectively.

To derive control strategies for the VPMSG and the VPCSG, the dynamic equations of the VSI and CSI in Figure 1 are formulated as follows.

VSI: $\left\{\begin{array}{l}{L}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{i}_{\mathrm{d}}}{\mathrm{d}t}}=-v_{\mathrm{d}}-R_{\mathrm{s}}{i}_{\mathrm{d}}+{\omega }_{\mathrm{e}}{L}_{\mathrm{s}}{i}_{\mathrm{q}}+u_{\mathrm{d}}\\ L_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{i}_{\mathrm{q}}}{\mathrm{d}t}}=-v_{\mathrm{q}}-R_{\mathrm{s}}{i}_{\mathrm{q}}-{\omega }_{\mathrm{e}}{L}_{\mathrm{s}}{i}_{\mathrm{d}}+u_{\mathrm{q}}\\ C{\displaystyle \frac{\mathrm{d}{v}_{\mathrm{dc}}}{\mathrm{d}t}}=I_{\mathrm{s}}-i_{\mathrm{dc}}-G{v}_{\mathrm{dc}}\end{array}\right.$                       (22)

CSI: $\left\{\begin{array}{l}{\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{d}}}{\mathrm{d}t}}=-{\widehat{i}}_{\mathrm{d}}-{\widehat{G}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{d}}+{\widehat{\omega }}_{\mathrm{e}}{\widehat{C}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{q}}+{\widehat{i}}_{\mathrm{id}}\\ {\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{q}}}{\mathrm{d}t}}=-{\widehat{i}}_{\mathrm{q}}-{\widehat{G}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{q}}-{\widehat{\omega }}_{\mathrm{e}}{\widehat{C}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{d}}+{\widehat{i}}_{\mathrm{iq}}\\ \widehat{L}{\displaystyle \frac{\mathrm{d}{\widehat{i}}_{\mathrm{L}}}{\mathrm{d}t}}={\widehat{V}}_{\mathrm{s}}-{\widehat{v}}_{\mathrm{dc}}-\widehat{R}{\widehat{i}}_{\mathrm{L}}\end{array}\right.$                      (23)

Here, u_d and u_q are the dq-axis components of the output voltages of the VSI in Figure 1(a); ${\widehat{i}}_{\mathrm{id}}$ and ${\widehat{i}}_{\mathrm{iq}}$ are the dq-axis output current components of the CSI in Figure 1(b).

Comparison of Table 2 and Table 3 with Equations (22) and (23) shows that the VSI and CSI are equivalent to the VPMSG and VPCSG under the following control strategies:

VPMSG: u_d^*=0 and u_q^*= n_pψ_rv_dc (24)

VPCSG: ${\widehat{i}}_{\mathrm{id}}^{*}=0$ and ${\widehat{i}}_{\mathrm{iq}}^{*}=n_{\mathrm{p}}{\widehat{q}}_{\mathrm{r}}{\widehat{i}}_{\mathrm{L}}$                     (25)

Here, u_d^* and u_q^* represent the dq-axis reference voltages of the VSI; ${\widehat{i}}_{\mathrm{id}}^{*}$ and ${\widehat{i}}_{\mathrm{iq}}^{*}$ denote the dq-axis reference currents of the CSI; n_p, ψ_r and ${\widehat{q}}_{\mathrm{r}}$ denote the control parameters. The control strategies for the VPMSG and VPCSG, given by Equations (24) and (25), are implemented via the SPWM modulation applied to the VSI and CSI in Figure 2. The corresponding dq-frame circuit equivalent models are derived as shown in Figure 3, confirming consistency with the models in Table 2 and Table 3. Here, ${\theta }_{\mathrm{e}}={\displaystyle \int n_{\mathrm{p}}{v}_{\mathrm{dc}}dt}$, ${\widehat{\theta }}_{\mathrm{e}}={\displaystyle \int n_{\mathrm{p}}{\widehat{i}}_{\mathrm{L}}dt},$

3. A Unified Modeling and Analysis Framework for Virtual Synchronous Generators

Based on Table 2 and Table 3, the VPMSG and VPCSG models are reformulated as follows:

VPMSG: $\left\{\begin{array}{l}{L}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{i}_{\mathrm{d}}}{\mathrm{d}t}}=-v_{\mathrm{d}}-R_{\mathrm{s}}{i}_{\mathrm{d}}+{\omega }_{\mathrm{e}}{L}_{\mathrm{s}}{i}_{\mathrm{q}}\\ L_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{i}_{\mathrm{q}}}{\mathrm{d}t}}=-v_{\mathrm{q}}-R_{\mathrm{s}}{i}_{\mathrm{q}}-{\omega }_{\mathrm{e}}{L}_{\mathrm{s}}{i}_{\mathrm{d}}+{\psi }_{\mathrm{r}}{\omega }_{\mathrm{e}}\\ C{\displaystyle \frac{\mathrm{d}{\omega }_{\mathrm{e}}}{\mathrm{d}t}}=n_{\mathrm{p}}{I}_{\mathrm{s}}-{\displaystyle \frac{3}{2}}{n}_{\mathrm{p}}{}^2\psi _{\mathrm{r}}{i}_{\mathrm{q}}-G{\omega }_{\mathrm{e}}\end{array}\right.$

VPCSG: $\left\{\begin{array}{l}{\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{d}}}{\mathrm{d}t}}=-{\widehat{i}}_{\mathrm{d}}-{\widehat{G}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{d}}+{\widehat{\omega }}_{\mathrm{e}}{\widehat{C}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{q}}\\ {\widehat{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{\widehat{v}}_{\mathrm{q}}}{\mathrm{d}t}}=-{\widehat{i}}_{\mathrm{q}}-{\widehat{G}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{q}}-{\widehat{\omega }}_{\mathrm{e}}{\widehat{C}}_{\mathrm{s}}{\widehat{v}}_{\mathrm{d}}+{\widehat{q}}_{\mathrm{r}}{\widehat{\omega }}_{\mathrm{e}}\\ \widehat{L}{\displaystyle \frac{\mathrm{d}{\widehat{\omega }}_{\mathrm{e}}}{\mathrm{d}t}}=n_{\mathrm{p}}{\widehat{V}}_{\mathrm{s}}-{\displaystyle \frac{3}{2}}{n}_{\mathrm{p}}{}^2\widehat{q}_{\mathrm{r}}{\widehat{v}}_{\mathrm{q}}-\widehat{R}{\widehat{\omega }}_{\mathrm{e}}\end{array}\right.$

Through the following transformations:

VPMSG: $t={\tau }_1\tilde{t}$ and $x={\xi }_1\tilde{x}$                       (26)

VPCSG: $t={\tau }_2\tilde{t}$ and $y={\xi }_2\tilde{x}$                       (27)

the unified model is derived as follows:

$$\left\{\begin{array}{l}{\displaystyle \frac{d{z}_1}{d\tilde{t}}}=-z_1+z_3{z}_2-{\Lambda }_{\mathrm{d}}\\ {\displaystyle \frac{d{z}_2}{d\tilde{t}}}=-z_2-z_3{z}_1+\gamma z_3-{\Lambda }_{\mathrm{q}}\\ {\displaystyle \frac{d{z}_3}{d\tilde{t}}}={\Gamma }_{\mathrm{m}}-\sigma (z_2+z_3)\end{array}\right.$$

(28)

where τ₁=L_s/R_s, $x={[i_{\mathrm{d}}\hspace{0.33em}\hspace{0.33em}{i}_{\mathrm{q}}\hspace{0.33em}\hspace{0.33em}{\omega }_{\mathrm{e}}]}^{\mathrm{T}}$, $\tilde{x}={[z_1\hspace{0.33em}\hspace{0.33em}{z}_{2\hspace{0.33em}}\hspace{0.33em}{z}_3]}^{\mathrm{T}}$, $k={\displaystyle \frac{2}{3}}{\displaystyle \frac{G}{n_{\mathrm{p}}{}^2\tau _1{\psi }_{\mathrm{r}}}}$, ${\xi }_1=\left[\begin{array}{l}k\hspace{1em}0\hspace{1em}0\\ 0\hspace{1em}\hspace{0.33em}\hspace{0.33em}k\hspace{1em}0\\ 0\hspace{1em}\hspace{0.33em}\hspace{0.33em}0\hspace{1em}1/{\tau }_1\end{array}\right]$, $\gamma ={\displaystyle \frac{{\psi }_{\mathrm{r}}}{k{L}_{\mathrm{s}}}}$, ${\Lambda }_{\mathrm{d}}={\displaystyle \frac{1}{R_{\mathrm{s}}k}}{v}_{\mathrm{d}}$, ${\Lambda }_{\mathrm{q}}={\displaystyle \frac{1}{R_{\mathrm{s}}k}}{v}_{\mathrm{q}}$, ${\Gamma }_{\mathrm{m}}={\displaystyle \frac{n_{\mathrm{p}}{\tau }_1^2}{C}}{I}_{\mathrm{s}}$ and $\sigma ={\displaystyle \frac{{\tau }_{1}G}{C}}$ for VPMSG; ${\tau }_2={\widehat{C}}_{\mathrm{s}}/{\widehat{G}}_{\mathrm{s}}$, $y={[{\widehat{v}}_{\mathrm{d}}\hspace{0.33em}\hspace{0.33em}{\widehat{v}}_{\mathrm{q}}\hspace{0.33em}\hspace{0.33em}{\widehat{\omega }}_{\mathrm{e}}]}^{\mathrm{T}}$, $\widehat{k}={\displaystyle \frac{2}{3}}{\displaystyle \frac{\widehat{R}}{n_{\mathrm{p}}{}^2\tau _2{\widehat{q}}_{\mathrm{r}}}}$,${\xi }_2=\left[\begin{array}{l}\widehat{k}\hspace{1em}0\hspace{1em}0\\ 0\hspace{1em}\hspace{0.33em}\hspace{0.33em}\widehat{k}\hspace{1em}0\\ 0\hspace{1em}\hspace{0.33em}\hspace{0.33em}0\hspace{1em}1/{\tau }_2\end{array}\right]$, $\gamma ={\displaystyle \frac{{\widehat{q}}_{\mathrm{r}}}{\widehat{k}{\widehat{C}}_{\mathrm{s}}}}$, ${\Lambda }_{\mathrm{d}}={\displaystyle \frac{1}{{\widehat{G}}_{\mathrm{s}}\widehat{k}}}{\widehat{i}}_{\mathrm{d}}$, ${\Lambda }_{\mathrm{q}}={\displaystyle \frac{1}{{\widehat{G}}_{\mathrm{s}}\widehat{k}}}{\widehat{i}}_{\mathrm{q}}$, ${\Gamma }_{\mathrm{m}}={\displaystyle \frac{n_{\mathrm{p}}{\tau }_2^2}{\widehat{L}}}{\widehat{V}}_{\mathrm{s}}$ and $\sigma ={\displaystyle \frac{{\tau }_2\widehat{R}}{\widehat{L}}}$ for VPCSG.

Furthermore, setting the time derivative in Equation (28) to zero yields the equilibrium point:

$$\left(\begin{array}{l}{z}_1^{*}\\ z_2^{*}\\ z_3^{*}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{{\Gamma }_{\mathrm{m}}{\xi }^{*}}{\sigma }}-{\xi }^{*}{}^2-{\Lambda }_{\mathrm{d}}\\ {\displaystyle \frac{{\Gamma }_{\mathrm{m}}}{\sigma }}-{\xi }^{*}\\ {\xi }^{*}\end{array}\right)$$

(29)

where ξ^* is the solution of the following equation:

$${\xi }^3+a_1{\xi }^2+a_2\xi +a_3=0$$

(30)

According to the cubic formula [23], the solutions of Equation (30) can be analytically derived as follows:

$$\left\{\begin{array}{l}{\xi }_1{}^{*}=T_{+}+T_{-}-{\displaystyle \frac{1}{3}}{a}_1\\ {\xi }_{2,3}{}^{*}=-{\displaystyle \frac{1}{2}}(T_{+}+T_{-})-{\displaystyle \frac{1}{3}}{a}_1\pm {\displaystyle \frac{1}{2}}\mathrm{i}\sqrt{3}(T_{+}-T_{-})\end{array}\right.$$

(31)

where $a_1=-{\displaystyle \frac{{\Gamma }_{\mathrm{m}}}{\sigma }}$, $a_2=1+\gamma +{\Lambda }_{\mathrm{d}}$, $a_3=-{\displaystyle \frac{{\Gamma }_{\mathrm{m}}+{\Lambda }_{\mathrm{q}}\sigma }{\sigma }}$, $Q={\displaystyle \frac{3a_2-a_1^2}{9}}$ and $R={\displaystyle \frac{9a_1{a}_2-27a_3-2a_1^3}{54}}$, $T_{+}=\sqrt[3]{R+\sqrt{\Delta }}$, $T_{-}=-Q/T_{+}$ and $\Delta =Q^3+R^2$. Therefore, three equilibria E₁^*, E₂^* and E₃^* correspond to ξ₁^*, ξ₂^* and ξ₃^*, respectively. When $\Delta \le 0$, the unified model (28) admits at least two distinct real equilibria; otherwise, it has a unique equilibrium.

In grid-connected mode, the grid parameters are specified as follows:

Grid-Side voltages of the VPMSG: v_a=U_msin(ω_gt+θ_v), v_b=U_msin(ω_gt+θ_v−2/3π) and v_c=U_msin (ω_gt+θ_v+2/3π), where θ_v, ω_g and U_m are the initial angle, electrical angular frequency and amplitude of the grid voltages. Here, ω_g=2πf_g, f_g denotes the grid frequency;

Grid-side currents of the VPCSG: ${\widehat{i}}_{\mathrm{a}}=I_m\sin ({\omega }_{\mathrm{g}}t+{\theta }_i)$, ${\widehat{i}}_{\mathrm{b}}=I_m\sin ({\omega }_{\mathrm{g}}t-{\displaystyle \frac{2}{3}}\pi +{\theta }_i)$ and ${\widehat{i}}_{\mathrm{c}}=I_m\sin ({\omega }_{\mathrm{g}}t+{\displaystyle \frac{2}{3}}\pi +{\theta }_i)$, where θ_i and I_m are the initial angle and amplitude of the grid currents, respectively.

Applying the Park transformation yields:

$\left\{\begin{array}{l}{v}_{\mathrm{d}}=U_{\mathrm{m}}\sin \left[\left({\omega }_{\mathrm{g}}-{\omega }_{\mathrm{e}}\right)t+{\theta }_{\mathrm{v}}\right]\\ v_{\mathrm{q}}=-U_{\mathrm{m}}\cos \left[\left({\omega }_{\mathrm{g}}-{\omega }_{\mathrm{e}}\right)t+{\theta }_{\mathrm{v}}\right]\end{array}\right.$ for the VPMSG and

$\left\{\begin{array}{l}{\widehat{i}}_{\mathrm{d}}=I_{\mathrm{m}}\sin \left[\left({\omega }_{\mathrm{g}}-{\widehat{\omega }}_{\mathrm{e}}\right)t+{\theta }_{\mathrm{i}}\right]\\ {\widehat{i}}_{\mathrm{q}}=-I_{\mathrm{m}}\cos \left[\left({\omega }_{\mathrm{g}}-{\widehat{\omega }}_{\mathrm{e}}\right)t+{\theta }_{\mathrm{i}}\right]\end{array}\right.$ for the VPCSG.

Whether the dq-frame grid voltages and currents are constant or time-varying depends on the relationship between ω_g and ω_e (or ω_g and ${\widehat{\omega }}_{\mathrm{e}}$). Standalone mode is a special case of grid-connected operation where ω_g=ω_e for the VPMSG and ω_g =${\widehat{\omega }}_{\mathrm{e}}$ for the VPCSG; consequently, v_d, v_q, ${\widehat{i}}_{\mathrm{d}}$ and ${\widehat{i}}_{\mathrm{q}}$ remain constant. Therefore, the stability of both the VPMSG and the VPCSG is analyzed under standalone and grid-connected modes.

3.1. Standalone Mode

3.1.1. Λ_d=Λ_q=Γ_m=0

By removing all excitation terms, Equation (28) simplifies to

$$\left\{\begin{array}{l}{\displaystyle \frac{d{z}_1}{d\tilde{t}}}=-z_1+z_3{z}_2\\ {\displaystyle \frac{d{z}_2}{d\tilde{t}}}=-z_2-z_3{z}_1+\gamma z_3\\ {\displaystyle \frac{d{z}_3}{d\tilde{t}}}=-\sigma (z_2+z_3)\end{array}\right.$$

(32)

Setting the time derivative in Equation (32) to zero yields the equilibrium point:

$\left(\begin{array}{l}{z}_1^{*}\\ z_2^{*}\\ z_3^{*}\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right)$

The Jacobian matrix is given by

$J_{\mathrm{ac}}=\left(\begin{array}{ccc}-1& z_3^{*}& z_2^{*}\\ -z_3^{*}& -1& \gamma -z_1^{*}\\ 0& -\sigma & -\sigma \end{array}\right)$

Solving the characteristic equation det(λI−J_ac)=0 yields the eigenvalues:

$\left\{\begin{array}{l}{\lambda }_1=-1\\ {\lambda }_{2,3}={\displaystyle \frac{-(1+\sigma )\pm \sqrt{{(\sigma -1)}^2-4\hspace{0.17em}\gamma \sigma }}{2}}\end{array}\right.$

The signs of the eigenvalues are analyzed as follows:

CASE 1: (σ−1)²−4γσ≥0

λ₂ and λ₃ are real eigenvalues. Since λ₂×λ₃=σ(1+γ)>0 and λ₂+λ₃=−(σ+1)<0, both λ₂ and λ₃ are negative real numbers.

CASE 2: (σ−1)²−4γσ<0

λ₂ and λ₃ form a complex conjugate pair with negative real parts since –(1+σ)/2<0.

Therefore, all eigenvalues have negative real parts, implying local asymptotic stability of the equilibrium point. The representative phase trajectory of system (32) is plotted in Figure 4.

3.1.2. Λ_d≠0, Λ_q≠0 and Γ_m≠0

In standalone mode, it can be concluded that ω_g=ω_e for the VPMSG and ω_g=${\widehat{\omega }}_{\mathrm{e}}$ for the VPCSG, so we can have Λ_d=constant and Λ_q=constant. For convenience of analysis, by applying the following transformations:

${\tilde{z}}_1=z_1-z_1^{*}$, ${\tilde{z}}_2=z_2-z_2^{*}$, ${\tilde{z}}_3=z_3-z_3{}^{*},$

Equation (28) is recast into the following dynamical system:

$$\dot{\tilde{z}}=\left(\begin{array}{c}{\dot{\tilde{z}}}_1\\ {\dot{\tilde{z}}}_2\\ {\dot{\tilde{z}}}_3\end{array}\right)=A\tilde{z}+N=\left(\begin{array}{ccc}-1& z_3^{*}& z_2^{*}\\ -z_3^{*}& -1& \gamma -z_1^{*}\\ 0& -\sigma & -\sigma \end{array}\right)\left(\begin{array}{c}{\tilde{z}}_1\\ {\tilde{z}}_2\\ {\tilde{z}}_3\end{array}\right)+\left(\begin{array}{c}{\tilde{z}}_2{\tilde{z}}_3\\ -{\tilde{z}}_1{\tilde{z}}_3\\ 0\end{array}\right)$$

(33)

Choose the Lyapunov function candidate as:

$V(\tilde{z})={\tilde{z}}^{\mathrm{T}}P\tilde{z}$

Here, $\tilde{z}$=(${\tilde{z}}_1$,${\tilde{z}}_2$,${ \mathrm{z} ˜}_3$)^T, p>0, q>0, r=q/p and

$P=\left(\begin{array}{ccc}p& 0& 0\\ 0& p& 0\\ 0& 0& q\end{array}\right).$

So, the time derivative of $V(\tilde{z})$ can be derived as:

$\begin{array}{l}\dot{V}(\tilde{z})={\dot{\tilde{z}}}^{\mathrm{T}}P\tilde{z}+{\tilde{z}}^{\mathrm{T}}P\dot{\tilde{z}}\\ =(A\tilde{z}+N{)}^{\mathrm{T}}P\tilde{z}+{\tilde{z}}^{\mathrm{T}}P(A\tilde{z}+N)\\ \hspace{0.33em} ={\tilde{z}}^{\mathrm{T}}(A^{\mathrm{T}}P+PA)\tilde{z}\\ =-{\tilde{z}}^{\mathrm{T}}Q\tilde{z}\end{array}$

Here,

$Q=-(A^{\mathrm{T}}P+PA)=\left(\begin{array}{ccc}2\hspace{0.17em}p& 0& -p\hspace{0.17em}{z}_2^{*}\\ 0& 2\hspace{0.17em}p& q\sigma -p\hspace{0.17em}\left(\gamma -z_1^{*}\right)\\ -p\hspace{0.17em}{z}_2^{*}& q\sigma -p\hspace{0.17em}\left(\gamma -z_1^{*}\right)& 2\hspace{0.17em}q\hspace{0.17em}\sigma \end{array}\right).$

The equilibrium point ${\tilde{\mathbf{z}}}^{*}$ is globally asymptotically stable if $\dot{V}$ is negative definite, which requires all leading principal minors of Q to be positive. The first leading principal minor Q₁=2p>0, The second leading principal minor is Q₂=4p²>0, and the third leading principal minor is

$$Q_3=\det (Q)=2p^3[4r\sigma -{(\gamma -z_1^{*}-r\sigma )}^2-{(z_2^{*})}^2]>0.$$

(34)

Consequently, the stability condition is established as follows:

$\Xi =4r\sigma -{(\gamma -z_1^{*}-r\sigma )}^2-{(z_2^{*})}^2>0$

It is known that $r^{*}=\arg \underset{r>0}{\max }\Xi (r)={\displaystyle \frac{2+(\gamma -z_1^{*})}{\sigma }}$, so we have

$$\Xi (r^{*})=4+4(\gamma -z_1^{*})-{(z_2^{*})}^2$$

(35)

Therefore, when Ξ>0, the unified system is locally asymptotically stable at the equilibrium point ${\tilde{\mathbf{z}}}^{*}$. The representative phase trajectory of the system (28) is plotted in Figure 5.

3.2. Grid-Connected Mode

In grid-connected mode, Equation (28) can be rewritten as:

$$\left\{\begin{array}{l}{\displaystyle \frac{d{z}_1}{d\tilde{t}}}=-z_1+z_3{z}_2-{\Theta }_{\mathrm{m}}\sin \left[\left({\tilde{\omega }}_{\mathrm{g}}-z_3\right)\tilde{t}+{\theta }_0\right]\\ {\displaystyle \frac{d{z}_2}{d\tilde{t}}}=-z_2-z_3{z}_1+\gamma z_3+{\Theta }_{\mathrm{m}}\cos \left[\left({\tilde{\omega }}_{\mathrm{g}}-z_3\right)\tilde{t}+{\theta }_0\right]\\ {\displaystyle \frac{d{z}_3}{d\tilde{t}}}={\Gamma }_{\mathrm{m}}-\sigma (z_2+z_3)\end{array}\right.$$

(36)

where ${\Theta }_{\mathrm{m}}={\displaystyle \frac{1}{R_{\mathrm{s}}k}}{U}_{\mathrm{m}}$, θ₀=θ_v and ${\tilde{\omega }}_{\mathrm{g}}={\omega }_{\mathrm{g}}{\tau }_1$ for the VPMSG; ${\Theta }_{\mathrm{m}}={\displaystyle \frac{1}{{\widehat{G}}_{\mathrm{s}}\widehat{k}}}{I}_{\mathrm{m}}$, θ₀=θ_i and ${\tilde{\omega }}_{\mathrm{g}}={\omega }_{\mathrm{g}}{\tau }_2$ for the VPCSG.

In steady state, the angular frequency in (36) is synchronized with the grid angular frequency; thus, ${\tilde{\omega }}_{\mathrm{g}}=z_3$. Therefore, the following conditions are satisfied:

$\varphi (\tilde{t})=\left({\tilde{\omega }}_{\mathrm{g}}-z_3\right)\tilde{t}+{\theta }_0$, ${\dot{z}}_3=0$ and $\dot{\varphi }(\tilde{t})=0$

Substitution of ${\dot{z}}_3=0$ into $0={\Gamma }_{\mathrm{m}}-\sigma (z_2+{\tilde{\omega }}_{\mathrm{g}})$ yields

$${\tilde{z}}_2{}^{*}={\displaystyle \frac{{\Gamma }_{\mathrm{m}}}{\sigma }}-{\tilde{\omega }}_{\mathrm{g}}$$

(37)

Furthermore, by setting ${\displaystyle \frac{d{z}_1}{d\tilde{t}}}=0$ and ${\displaystyle \frac{d{z}_2}{d\tilde{t}}}=0$, it is derived that

$\left\{\begin{array}{l}0=-z_1+{\tilde{\omega }}_{\mathrm{g}}{\tilde{z}}_2{}^{*}-{\Theta }_{\mathrm{m}}\sin {\theta }_0\\ 0=-{\tilde{z}}_2{}^{*}-{\tilde{\omega }}_{\mathrm{g}}{z}_1+\gamma {\tilde{\omega }}_{\mathrm{g}}+{\Theta }_{\mathrm{m}}\cos {\theta }_0\end{array}\right.$

So,

$$\left\{\begin{array}{l}{z}_1={\tilde{\omega }}_{\mathrm{g}}{\tilde{z}}_2{}^{*}-{\Theta }_{\mathrm{m}}\sin {\theta }_0\\ {\tilde{z}}_2{}^{*}+{\tilde{\omega }}_{\mathrm{g}}{z}_1-\gamma {\tilde{\omega }}_{\mathrm{g}}={\Theta }_{\mathrm{m}}\cos {\theta }_0\end{array}\right.$$

(38)

Substituting the first expression in (38) into the second yields

$${\displaystyle \frac{{\tilde{z}}_2{}^{*}(1+{\tilde{\omega }}_{\mathrm{g}}{}^2)-\gamma {\tilde{\omega }}_{\mathrm{g}}}{{\Theta }_{\mathrm{m}}\sqrt{1+{\tilde{\omega }}_{\mathrm{g}}^2}}}=\sin ({\theta }_0+\arctan ({\displaystyle \frac{1}{{\tilde{\omega }}_{\mathrm{g}}}}))$$

(39)

From Equation (39), the initial angle θ₀ admits two analytical solutions:

$\left\{\begin{array}{l}{\theta }_0{}^{(1)}=-\arctan ({\displaystyle \frac{1}{{\tilde{\omega }}_{\mathrm{g}}}})+\arcsin {\displaystyle \frac{{\tilde{z}}_2{}^{*}(1+{\tilde{\omega }}_{\mathrm{g}}{}^2)-\gamma {\tilde{\omega }}_{\mathrm{g}}}{{\Theta }_{\mathrm{m}}\sqrt{1+{\tilde{\omega }}_{\mathrm{g}}^2}}}\\ {\theta }_0{}^{(2)}=-\arctan ({\displaystyle \frac{1}{{\tilde{\omega }}_{\mathrm{g}}}})+\pi -\arcsin {\displaystyle \frac{{\tilde{z}}_2{}^{*}(1+{\tilde{\omega }}_{\mathrm{g}}{}^2)-\gamma {\tilde{\omega }}_{\mathrm{g}}}{{\Theta }_{\mathrm{m}}\sqrt{1+{\tilde{\omega }}_{\mathrm{g}}^2}}}\end{array}\right.$

Substituting Equation (37) into Equation (39) yields

$$\Pi =\left|{\displaystyle \frac{(\frac{{\Gamma }_{\mathrm{m}}}{\sigma }-{\tilde{\omega }}_{\mathrm{g}})(1+{\tilde{\omega }}_{\mathrm{g}}{}^2)-\gamma {\tilde{\omega }}_{\mathrm{g}}}{{\Theta }_{\mathrm{m}}\sqrt{1+{\tilde{\omega }}_{\mathrm{g}}^2}}}\right|\le 1$$

(40)

Equation (40) gives the necessary condition for the angular frequency of the unified model to be locked to the grid angular frequency, as illustrated for the VPMSG and the VPCSG in Figure 6. It can be seen that the regions corresponding to ${\tilde{\omega }}_{\mathrm{g}}=z_3$ should appear among the regions with Π≤1. However, when Π≤1, the stability of the unified model requires further discussion, as follows.

In the frequency-locked state, the equilibrium point of the grid-connected unified system (36) can be rewritten as

$$\left\{\begin{array}{l}{\tilde{z}}_1{}^{*}={\tilde{\omega }}_{\mathrm{g}}{\tilde{z}}_2{}^{*}-{\Theta }_{\mathrm{m}}\sin {\theta }_0\\ {\tilde{z}}_2{}^{*}={\displaystyle \frac{{\Gamma }_{\mathrm{m}}}{\sigma }}-{\tilde{\omega }}_{\mathrm{g}}\\ {\tilde{z}}_3{}^{*}={\tilde{\omega }}_{\mathrm{g}}\end{array}\right.$$

(41)

Here, there are two equilibrium points E₁^* and E₂^* corresponding to θ₀⁽¹⁾ and θ₀⁽²⁾.

Therefore, the Jacobian matrix can be obtained as

$${\tilde{J}}_{\mathrm{ac}}=\left(\begin{array}{ccc}-1& {\tilde{\omega }}_{\mathrm{g}}& {\tilde{z}}_2^{*}\\ -{\tilde{\omega }}_{\mathrm{g}}& -1& \gamma -{\tilde{z}}_1^{*}\\ 0& -\sigma & -\sigma \end{array}\right)$$

(42)

The characteristic equation $\det (\lambda I-{\tilde{J}}_{\mathrm{ac}})=0$ is derived as

$${\lambda }^3+b_1{\lambda }^2+b_2\lambda +b_3=0$$

(43)

Here, $b_1=\sigma +2$, $b_2={\tilde{\omega }}_{\mathrm{g}}{}^2+2\hspace{0.17em}\sigma +\sigma (\gamma \hspace{0.17em}-{\tilde{z}}_1{}^{*})+1$, $b_3=\sigma \hspace{0.17em}({\tilde{\omega }}_{\mathrm{g}}{}^2-\hspace{0.17em}{\tilde{z}}_2{}^{*}\tilde{\omega }_{\mathrm{g}}+1+\gamma \hspace{0.17em}-{\tilde{z}}_1{}^{*})$.

Applying the cubic formula, the eigenvalues are analytically obtained as follows:

$$\left\{\begin{array}{l}{\lambda }_1={\widehat{T}}_{+}+{\widehat{T}}_{-}-{\displaystyle \frac{1}{3}}{b}_1\\ {\lambda }_{2,3}=-{\displaystyle \frac{1}{2}}({\widehat{T}}_{+}+{\widehat{T}}_{-})-{\displaystyle \frac{1}{3}}{b}_1\pm {\displaystyle \frac{1}{2}}\mathrm{i}\sqrt{3}({\widehat{T}}_{+}-{\widehat{T}}_{-})\end{array}\right.$$

(44)

where$\widehat{Q}={\displaystyle \frac{3b_2-b_1^2}{9}}$ and $\widehat{R}={\displaystyle \frac{9b_1{b}_2-27b_3-2b_1^3}{54}}$, ${\widehat{T}}_{+}=\sqrt[3]{\widehat{R}+\sqrt{\widehat{\Delta }}}$, ${\widehat{T}}_{-}=-\widehat{Q}/{\widehat{T}}_{+}$ and $\widehat{\Delta }={\widehat{Q}}^3+{\widehat{R}}^2$. Providing that all eigenvalues possess negative real parts, the system is asymptotically stable.

In the stable frequency-locked state, the typical phase trajectory of the unified model (36) is shown in Figure 7. The system is stabilized at the equilibrium point E₂^*. The state variable z₃ converges to the grid angular frequency ${\tilde{\omega }}_{\mathrm{g}}$ as shown in Figure 8.

4. Numerical Simulations

Without loss of generality, the operating parameters of the unified model are chosen as follows: L_s=2mH, R_s=200mΩ, C=470μF and G=2.1S for the VPMSG; ${\widehat{C}}_{\mathrm{s}}=200\mathsf{\mu }\mathrm{F}$, ${\widehat{G}}_{\mathrm{s}}=20\mathrm{mS}$, $\widehat{L}=8\mathrm{mH}$ and $\widehat{R}=50\mathrm{m}\mathsf{\Omega }$ for the VPCSG. To verify the stability of the unified model, the system is analyzed in standalone mode and grid-connected mode, respectively.

3.1. Standalone Mode

Based on Equation (35), Ξ is calculated as 223.5>0 for the VPMSG and 20.6>0 for the VPCSG, so they are stable. The typical waveforms are shown in Figure 9. As observed, the output three-phase currents and voltages are symmetrical, and the output angular frequency is stabilized at f_g=50Hz for both generators.

3.2. Grid-Connected Mode

Under the control strategies (24) and (25), the VSI and the CSI function as the VPMSG and the VPCSG, respectively. When the unified model (36) falls within the frequency-locked regions shown in Figure 6, its stability must be further verified through eigenvalue analysis to ensure that the VSG outputs stable three-phase currents or voltages. According to Equation (44), the real part of the rightmost eigenvalue λ_max is numerically computed as a function of U_m and I_m, as illustrated in Figure 10. Therefore, when Π<1 and Re(λ_max)<0, the system is stabilized at the grid frequency f_g=50Hz. The typical waveforms are illustrated in Figure 11. Obviously, the electrical angular frequency ω_e can maintain synchronization with the grid frequency across a broad operating range of U_m and I_m.

5. Conclusions

The VSI and CSI are dual topologies that share a similar mathematical model. Under specific control conditions, the VSI is mathematically equivalent to a VPMSG; likewise, the CSI is mathematically equivalent to a VPCSG. Therefore, the VPMSG and VPCSG form a pair of dual VSGs that share an analogous mathematical model. A unified dimensionless model can be established to characterize their dynamic behaviors. By selecting an appropriate Lyapunov function candidate, the stability conditions for the unified system in standalone mode are analytically derived. Through eigenvalue analysis, the stability conditions for the unified model in grid-connected mode are also derived. Nevertheless, grid-connected operation necessitates that the VSGs remain locked to the grid frequency. Both frequency synchronization and stability conditions must be simultaneously satisfied to enable grid connection of the VSGs.