Figure 1 shows a HVDC system connected with two WFs. This paper mainly focuses on the power sending terminal containing the rectifier of the HVDC and the GSCs of the WFs. In such a system, P and Q represent active power and reactive power. E ̇ and I ̇ represent the voltage and current vectors of each equipment. U t ̇ represents the voltage vector at the point of common coupling. X represents the reactance between each equipment and the grid connection point. Subscripts1, 2, 3, g represent GSC 1, GSC 2, a rectifier, and infinite bus, respectively.

According to (Lu et al., 2020), Figure 2A shows the typical control strategy of the rectifier of a HVDC system. DC current control (DCC) is applied to control the DC current and generate an appropriate order triggering angle α _ord for the thyristors. At the same time, a phase-locked loop (PLL) is used to capture the phase of the AC terminal voltage. The value of α _ord is the sum of actual triggering angle and PLL’s angle.

Figure 2B depicts the typical control strategy of a GSC, including DC voltage control (DVC), reactive power control (RPC), inner current control (ICC) and PLL. The active branch controls the DC voltage and the AC current in d-axis. The reactive branch controls the reactive power and the AC current in q-axis. Similar with the rectifier of HVDC, all the above control are controlled on the PLL synchronization reference frame.

Utilizing the above control strategies, the small-signal models of a rectifier and a GSC in DVC timescale (about 10 Hz) are proposed in (Yuan et al., 2017; Lu et al., 2020), respectively. Figure 3A depicts the original small-signal model of a rectifier. Figure 3B depicts the original small-signal model of a GSC. In these models, the active and reactive power are the input signals and the phase and magnitude of the internal voltage or current are the output signals. They describe the dynamic characteristics of the rectifier and the GSC with clear mechanism understanding based on the motion equation concept.

The basic theory of ν-gap metric is proposed (Zhou, 2010) to analyze the robust stability of a feedback system. Consider a nominal feedback system depicted in Figure 4A. P (s) is the forward channel and C (s) is the feedback channel. The robust stability margin to quantify the relative stability is denoted as b [ P (s), C (s)]. If the system is stable, b [ P (s), C (s)] is in the interval (0,1]. If the system is unstable, b [ P (s), C (s)] = 0.

For the corresponding feedback system with uncertainties as shown Figure 4B. P ̃ ( s ) and C ̃ ( s ) are sets of uncertainties around the given system P (s) and C (s). Denote the distance between P (s) and P ̃ ( s ) as ν P ( s ) , P ̃ ( s ) . Denote the distance between C (s) and C ̃ ( s ) as ν C ( s ) , C ̃ ( s ) . If ν P s , P ̃ s + ν C s , C ̃ s < b P s , C s (1) the uncertain system is stable. Otherwise, the uncertain system is unstable.

It is worth noticing that the system can be a SISO system or a MIMO system. Furthermore, the stability criterion is a necessary and sufficient condition without any conservatism.

In this subsection, we build the forward channel of the MIMO system model reflecting the dynamic characteristics of the rectifier and GSCs.

On the basis of the rectifier model shown in Figure 3A the GSC model shown in Figure 3B, we establish the small-signal model of the HVDC system connected with two GSCs in Figure 5. In this model, ΔP _in represents the active power input of the DC inductance and DC capacitors of each device. ΔP _out and ΔQ _out represent the active and reactive power output of the DC inductance and DC capacitors of each device. Δθ and ΔE represent the phase and magnitude of the internal voltage. Denote Δθ and ΔI as phase and magnitude of the internal current of the rectifier. Transfer functions M(s) and D(s) represent the equivalent inertia and damping, respectively. G _EQ (s) is the transfer function related to the reactive power control. G _Pr3(s) is a partial transfer function between ΔP _out and Δθ of the rectifier. Denote Δi ₃ as AC current of the rectifier. G _Pθr (s) is a partial transfer function between Δi ₃ and Δθ of the rectifier. K ₃ is a proportionality transfer function between DC current dynamics and ΔI of the rectifier. J _{6 × 6} is an admittance matrix of the power flow in AC network. Specific expressions of the matrix and the transfer functions are provided in Supplementary Appendix SA according to (Yuan et al., 2017; Lu et al., 2020).

For the GSCs in the model, the transfer function between Δθ of the internal voltage and ΔP _out of GSC 1 can be expressed as G P θ 1 s = s D 1 s + 1 s 2 M 1 s . (2) The transfer function between Δθ of the internal voltage and ΔP _out of GSC 2 can be described as G P θ 2 s = s D 2 s + 1 s 2 M 2 s . (3) For the rectifier in the model, denote the transfer function G _Pθ3(s) between Δθ of the internal current and ΔP _out of the rectifier as G P θ 3 s = − s D 3 s + 1 s 2 M 3 s G P θ r s + G P r 3 s . (4) The transfer function between ΔI and ΔP _out of the rectifier can be described as G P I 3 s = K 3 s M 3 s . (5) Specific expressions of G _EQ1, G _EQ2, G _Qθ are given in (Lu et al., 2020). For compact expressions, denote Δ V = Δ θ 1 Δ θ 2 Δ θ 3 Δ E 1 Δ E 2 Δ I 3 T , Δ S = Δ P o u t 1 Δ P o u t 2 Δ P o u t 3 Δ Q 1 Δ Q 2 Δ Q 3 T . (6) Then, the forward channel of the MIMO system can be described as Δ V = H s Δ S . (7) where H ( s ) = G P θ 1 ( s ) 0 0 0 0 0 0 G P θ 2 ( s ) 0 0 0 0 0 0 G P θ 3 ( s ) 0 0 G Q θ ( s ) 0 0 0 G E Q 1 ( s ) 0 0 0 0 0 0 G E Q 2 ( s ) 0 0 0 G P I 3 ( s ) 0 0 0 .

Next, let’s establish the model of the AC network, i.e., the feedback channel of the MIMO system.

As seen from Figure 5, the feedback channel reflects the power flow of the AC network. Since the input signals of the AC network contains the phase and magnitude of both the internal voltage and internal current, we need to derive the Jacobian matrix of the mixed system.

First, we calculate the branch circuit across each device. For the system shown in Figure 1, the current of each branch can be expressed as I ̇ 1 = E ̇ 1 − U ̇ t j X 1 , I ̇ 2 = E ̇ 2 − U ̇ t j X 2 , I ̇ 3 = I ̇ 3 , I ̇ g = U ̇ g − U ̇ t j X g . (8) According to the Kirchhoff’s Current Law, the current flowing through the rectifier is I ̇ 1 + I ̇ 2 + I ̇ g = I ̇ 3 . (9) Based on (8) and (9), the AC voltage on the bus can be calculated as U ̇ t = − j X 1 X 2 X g  I ̇ 3 + E ̇ 1 X 2 X g  + E ̇ 2 X 1 X g  + U ̇ g  X 1 X 2 X 2 X g + X 1 X g + X 1 X 2 . (10) By substituting (10) into (8), the current of the two GSCs can be obtained as I 1 ̇ = Y 11 E 1 ̇ + Y 12 E 2 ̇ + G 13 I 3 ̇ + Y 1 g U g ̇ , I 2 ̇ = Y 21 E 1 ̇ + Y 22 E 2 ̇ + G 23 I 3 ̇ + Y 2 g U g ̇ , (11) where Y ₁₁, Y ₁₂, G ₁₃, Y _1g , Y ₂₁, Y ₂₂, G ₂₃, Y _2g are the elements of the admittance matric as shown in Supplementary Appendix SA. The internal voltage of the rectifier can be described as E ̇ 3 = U ̇ t − j I ̇ 3 X 3 . (12) By substituting (10) into (12), we have E ̇ 3 = Y 31 E 1 ̇ + Y 32 E 2 ̇ + G 33 I 3 ̇ + Y 3 g U g ̇ , (13) where the elements of the admittance matric Y ₃₁, Y ₃₂, G ₃₃, Y _3g are as shown in Supplementary Appendix SA.

Then, the apparent power of the GSCs and the rectifier can be calculated as. S 1 = E 1 ̇ I 1 ̇ ̄ = P o u t 1 + j Q 1 , (14) S 2 = E 2 ̇ I 2 ̇ ̄ = P o u t 2 + j Q 2 , (15) S 3 = E 3 ̇ I 3 ̇ ̄ = P o u t 3 + j Q 3 , (16) where I 1 ̇ ̄ , I 2 ̇ ̄ , I 3 ̇ ̄ are conjugate of I 1 ̇ , I 2 ̇ , I 3 ̇ .

According to (11) and (13) into (14), the active power and reactive power of each branch can be expressed as. P o u t 1 = G 13 E 1 I 3 ⁡ cos θ 1 − θ 3 − j Y 12 E 1 E 2 ⁡ sin θ 1 − θ 2 − j Y 1 g E 1 U g ⁡ sin θ 1 , (17) Q 1 = j Y 11 E 1 2 + G 13 E 1 I 3 ⁡ sin θ 1 − θ 3 + j Y 12 E 1 E 2 × cos θ 1 − θ 2 + j Y 1 g E 1 U g ⁡ cos θ 1 , (18) P o u t 2 = G 23 E 2 I 3 ⁡ cos θ 2 − θ 3 − j Y 21 E 1 E 2 ⁡ sin θ 2 − θ 1 − j Y 1 g E 2 U g ⁡ sin θ 2 , (19) Q 2 = j Y 22 E 2 2 + G 23 E 2 I 3 ⁡ sin θ 2 − θ 3 + j Y 21 E 1 E 2 × cos θ 2 − θ 1 + j Y 2 g E 2 U g ⁡ cos θ 2 , (20) P o u t 3 = Y 31 E 1 I 3 ⁡ cos θ 1 − θ 3 + Y 32 E 2 I 3 × cos θ 2 − θ 3 + Y 3 g I 3 U g ⁡ cos θ 3 , (21) Q 3 = Y 31 E 1 I 3 ⁡ sin θ 1 − θ 3 + Y 32 E 2 I 3 ⁡ sin θ 2 − θ 3 − Y 3 g U g I 3 ⁡ sin θ 3 − j G 33 I 3 2 . (22) Linearizing (17) at the steady-state operating point yields. Δ P o u t 1 = K P θ 11 Δ θ 1 + K P θ 12 Δ θ 2 + K P θ 13 Δ θ 3 + K P E 11 Δ E 1 + K P E 12 Δ E 2 + K P I 13 Δ I 3 , (23) Δ Q 1 = K Q θ 11 Δ θ 1 + K Q θ 12 Δ θ 2 + K Q θ 13 Δ θ 3 + K Q E 11 Δ E 1 + K Q E 12 Δ E 2 + K Q I 13 Δ I 3 , (24) Δ P o u t 2 = K P θ 21 Δ θ 1 + K P θ 22 Δ θ 2 + K P θ 23 Δ θ 3 + K P E 21 Δ E 1 + K P E 22 Δ E 2 + K P I 23 Δ I 3 , (25) Δ Q 2 = K Q θ 21 Δ θ 1 + K Q θ 22 Δ θ 2 + K Q θ 23 Δ θ 3 + K Q E 21 Δ E 1 + K Q E 22 Δ E 2 + K Q I 23 Δ I 3 , (26) Δ P o u t 3 = K P θ 31 Δ θ 1 + K P θ 32 Δ θ 2 + K P θ 33 Δ θ 3 + K P E 31 Δ E 1 + K P E 32 Δ E 2 + K P I 33 Δ I 3 , (27) Δ Q 3 = K Q θ 31 Δ θ 1 + K Q θ 32 Δ θ 2 + K Q θ 33 Δ θ 3 + K Q E 31 Δ E 1 + K Q E 32 Δ E 2 + K Q I 33 Δ I 3 . (28) In (23), K _Pθ11 ∼ K _Pθ33, K _PE11 ∼ K _PI33, K _Qθ11 ∼ K _Qθ33, K _QE33 ∼ K _QI33 are shown in Supplementary Appendix SA. Denote K P θ s = − K P θ 11 − K P θ 12 − K P θ 13 − K P θ 21 − K P θ 22 − K P θ 23 K P θ 31 K P θ 32 K P θ 23 , K P E s = − K P E 11 − K P E 12 − K P l 13 − K P E 21 − K P E 22 − K P l 23 K P E 31 K P E 32 K P l 33 , K Q θ s = K Q θ 1 K Q θ 12 K Q θ 13 K Q θ 21 K Q θ 22 K Q θ 23 K Q θ 31 K Q θ 32 K Q θ 33 , K Q E s = K Q E 11 K Q E 12 K Q l 13 K Q E 21 K Q E 22 K Q l 23 K Q E 31 K Q E 32 K Q l 33 . (29) Then, the feedback channel (23) can be written as Δ S = L s Δ V . (30) where L ( s ) = K P θ ( s ) K P E ( s ) K Q θ ( s ) K Q E ( s ) .

Then, the entire model of the system represented by transfer function matrix can be illustrated by Figure 6. In this model, Δ S represent the set of the input signals and Δ V represent the set of the output signals. H (s) and L (s) represent the forward and feedback channels, respectively.

In this subsection, we propose a new index, i.e., the robust stability margin based on the ν-gap metric to evaluate the relative stability of the MIMO system.

For a standard feedback system as depicted in Figure 6A, if the system is stable, the robust stability margin b _H,L of Figure 6A are defined as b H , L = max ω ∈ R ∪ ∞ σ ̄ H j ω I I − L j ω H j ω − 1 − L j ω I − 1 . (31) where σ ̄ refers to the maximum singular value. For a MIMO system represented by the transfer function matrix G(s), the maximum singular value equals to the maximum eigenvalue of the product of the transfer function matrix G (jω) and its conjugate transpose G *(jω). In addition, H _∞ norm denotes the maximum singular value of a system represented by transfer function matrix. Therefore, for the MIMO system with H (s) as the forward channel and L (s) as the feedback channel as shown in Figure 6A, if the system is stable, the robust stability margin can be expressed as b H , L = H s I I − L s H s − 1 − L s I ∞ − 1 . (32)

The proposed robust stability margin can be used to assess the relative stability of the system. If the system is stable, b _H,L is in the interval (0,1]. The larger b _H,L is, the more stable the system is. When the system is unstable, b _H,L equals to 0. Compared with other stability margin, i.e., phase margin and amplitude margin, the proposed robust stability margin can be directly computed in a MIMO systems.

In this subsection, we will clarify the methodology of calculating the stable region of parameters. From the point view of control, the stable region for the control systems can be obtained by calculating the ranges of uncertain parameters preserving the system stability. In this paper, we mainly focus on assessing the stable region of the control parameters in DVC, DCC and PLL of the GSCs and the rectifier. The structure and the power flow of the network are considered to be fixed. It means that the uncertainties exist only in H (s).

Consider the feedback system with uncertainties as shown in Figure 6B. Denote H ̃ ( s ) as the set of uncertainties around the given system H (s). Denote the distance between H (s) and H ̃ ( s ) as ν H ( s ) , H ̃ ( s ) . According to the theory of ν-gap (Zhou, 2010), the ν-gap between the uncertain system in Figure 6B and the nominal system in Figure 6A can be calculated as ν H s , H ̃ s = 1 − H s I I + H − s T H s − 1 H − s T I ∞ − 2 1 2 . (33)

Now, we present the stability criterion of the MIMO system with uncertainties in Figure 6B. Suppose the nominal MIMO system in Figure 6A is stable with the stability margin ν H ( s ) , H ̃ ( s ) . The uncertain system in Figure 3B is stable if and only if a r c s i n ν H s , H ̃ s < b H , L . (34) Otherwise, the uncertain system is unstable.

We can use the proposed stability criterion in (34) to assess the stable region of parameters. Figure 7 shows the algorithm for calculating the stable region of parameters. Specifically, the flowchart is divided into three steps as follows.

Step 1 Check the stability of the system. If the system is unstable, set b _H,L = 0. Otherwise, compute b _H,L of the nominal MIMO system in Figure 6A and continue.

Step 2 Assume the parameter k = k ₀ in the nominal H (s) and the corresponding uncertain parameter k = k ₁. Calculate the ν-gap distance ν H ( s ) , H ̃ ( s ) .

Step 3 If a r c s i n ν H ( s ) , H ̃ ( s ) < b H , L , the uncertain system with k = k ₁ is stable. Otherwise, the system with k = k ₁ is unstable.

Step 4 Using the incremental search, find the smallest k ₁ and the largest k ₁ so that the system is stable. Then the stable region for parameter k is (k _1smallest , k _1largest ).

Recall the HVDC system connected with GSCs as shown in Figure 1. Common operation points and control parameters are provided in Supplementary Appendix SA. With five different sets of control parameters, we calculate the robust stability margin of the systems. The results are shown in Table 1. We can see that the second system has the largest robust stability margin, which means it has the best stability performance. The fifth system has the smallest robust stability margin, which indicates the worst stability performance.

Next, we compare the results with eigenvalue analysis. Table 2 shows the results of the eigenvalues and participation factors with two given sets of parameters. We can see that both of the two systems are stable since the eigenvalues are all in the left half plane. In addition, the real parts of the dominant poles of the two systems are close to each other. It seems that the relative stability of the systems are close. However, the relative stability of the systems are different according to the robust stability margin calculated by our proposed method. It can be verified by time-domain simulation based on detailed nonlinear model. According to the participation factors, k _Ppllr and k _p1a have the greatest influence on the dominant poles of the two systems, respectively. Figure 8 shows the time domain response of θ ₃ of the rectifier with different k _Ppllr . When k _Ppllr changes from 2 to 0.05, the system remains stable. Figure 9 shows the time domain responses of θ ₁ of GSC 1 with different k _p1a . When k _p1a changes from 2.3 to 0.23, the system changes from stable to unstable. The results show that the system with k _Ppllr = 2 and k _p1a = 50 can tolerate a larger range of parameter variation than the system with k _Ppllr = 20 and k _p1a = 2.3, which means the former is more robust.

In this subsection, we show how to calculate the stable region of the parameters of k _p1a and X _g in the MIMO system.

First, we focuse on the proportional parameter k _p1a of PLL control in GSC 1. Consider the MIMO system with parameters given in Supplementary Appendix SA. According to the expressions of robust stability margin (32), the robust stability margin of the system is 0.19. Based on the algorithm proposed in Section 2.3.2, when k _p1a varies, we can compute the ν-gap ν [ H ( s ) , H ̃ ( s ) ] between H (s) and H ̃ ( s ) . The results are shown in Table 3. It can be seen that ν [ H ( s ) , H ̃ ( s ) ] gradually increases when k _p1a decreases. ν [ H ( s ) , H ̃ ( s ) ] equals to 0.1791, which is close to 0.19 when k _p1a = 0.6. ν [ H ( s ) , H ̃ ( s ) ] is slightly over the stability margin 0.19 when k _p1a = 0.1. Time-domain simulation based on the nonlinear model in Figure 10 verifies the results. The system is in a critically stable state when k _p1a = 0.6. The system oscillates when k _p1a = 0.1. Therefore, to maintain the stability of the system, the value of the proportional parameter k _p1a of PLL control of GSC 1 should be greater than 0.6.

Then, we compute the stable region of the reactance between the infinite bus and the point of interconnection X _g . As we know, X _g affects the short-circuit ratio (SCR) of the system. With the parameters shown in Supplementary Appendix SA (where X _g = 0.35), the robust stability margin of the system is 0.19. Table 4 shows the ν-gap ν-gap ν [ H ( s ) , H ̃ ( s ) ] between H (s) and H ̃ ( s ) when X _g changes. We can see that ν [ H ( s ) , H ̃ ( s ) ] increases when X _g increases. ν [ H ( s ) , H ̃ ( s ) ] equals to 0.1846, which is close to 0.19 when X _g = 0.43. It indicates that the system is critically stable with the corresponding SCR = 3.49. Then, when X _g = 0.44, ν [ H ( s ) , H ̃ ( s ) ] is larger than 0.19, which means that the system is unstable. Time domain simulation shown in Figure 11 based on the nonlinear model verifies the results.

The above two cases show the effectiveness of the proposed method to compute the stable regions of parameters. It is worth noticing that the results do not have conservatism since the stability criterion is a necessary and sufficient condition.