Black-Box Impedance Prediction of Grid-Tied VSCs Under Variable Operating Conditions

Impedance/admittance models (IM/AMs) have been widely used to analyze the small-signal stability of grid-tied power electronic devices, such as the voltage source converter (VSC). They can be either derived from theoretical analysis under white-box conditions, where all parameters and control structures are fully known, or measured based on experiments under black-box conditions, where the topology and parameters of the controllers are completely unknown. As the IM/AMs depend on specific operating conditions, it is highly desirable to develop fast algorithms for IM/AM prediction (or estimation) under the black-box and variable-operating-point conditions. This article extends the nearly-decoupled AM method for sequence AMs proposed recently by Liu et al to fit any unknown control structure, including not only grid-following VSC, but also grid-forming VSC. It is, therefore, referred to as the fully-decoupled IM (FDIM) method. Furthermore, a curve fitting method for the transfer function is proposed to expedite the algorithm, based on the information of a few disturbance frequencies only. Finally, the algorithm is verified by wide simulations and experiments under different situations, including the direct-drive wind turbine generator. The whole approach is expected to be broadly applicable to the stability analysis of power-electronic-based power systems under variable operating conditions.

The intermediate variables in (9). Z The impedance matrix, including Z dd , Z dq , Z qd , Z qq .

gv
The intermediate variables in (12). a kn , b kn The all coefficients in (13). x The operating point vector in (15). a k , b k The system parameters in (15).

R,R r
The intermediate variables in (33) and its real part. C The DC capacitance of VSC. r The number of operaing points. Q Output reactive power. ω 0 Fundamental frequency of the system. k p,ip Proportional and integral coefficients of PI controller in the power flow control. K p,q Gain coefficients of the basic droop control. k p,iv Proportional and integral coefficient of PI controller in the basic droop control. D p,q , J Gain coefficient of the virtual synchronous generator control.

N f
The number of disturbance frequencies in Equ. (36).

Z FSIM
The Impedance calculated by frequency-scan impedanc model in (36).

E j
The intermediate variables in (44).

I. INTRODUCTION
Due to increasing pressure of environmental protection and energy resource, a great quantity of solar and wind powers with the form of distributed sources have been connected to grid recently [1], and meanwhile the traditional power system has being gradually transformed into a semiconducting power system (also called power-electronic-based power system) [2]- [5]. As one of most common devices, the voltage source converter (VSC) has been widely used in photovoltaic inverters, wind farms, and high voltage dc transmission system so on. Recently we have been faced with various types of serious stability and oscillation problems in system level caused by power electronic devices, and hence accurate mathematical modeling and analysis of VSC is highly desirable [6]- [8].
Until now there are two major methods for the small-signal stability analysis of power-electronic-based power systems, including the state-space model and the impedance/admittance model (IM/AM) [9]- [11]. The state-space model is a time-domain analytic method based on the system state equation. The state equation is usually established according to the circuit configuration and then linearized on the steady-state operating points. Further the stability of the system can be studied by analyzing the associated Jacobian matrix. But the state-space model is workable only under white-box conditions, in which the system structure and parameters should be fully known. In contrast, generally the IM/AM can be divided into two types. The first is the IM/AM based on theoretical modeling. The transfer function is established according to the system topology to obtain the theoretical IM/AM. Unfortunately, it is also only applicable under white-box conditions. The second is the IM/AM based on measurement. There the IM/AM can be obtained through experiments by frequency sweeping measurement via either series voltage or shunt current injection, under black/graybox conditions [12]- [15], where usually we cannot get all the structure and parameters of devices for the sake of commercial secrecy. Based on these IM/AM results from either theoretical modeling or measurement, the small-signal stability analysis can be further studied based on the generalized Nyquist stability criterion [16]- [19]. As the black-box or gray-box is more realistic, the IM/AM measurement has been widely used under various practical conditions and has become a dominant method in electric power engineering. In addition, the method of amplitude-phase motion equations has been proposed [20], and the equivalence and relation of all these important methods have been investigated very recently [21], [22].
It is well known that the small-signal stability relies on a linearization around the sufficiently small neighboring of a certain operating point, and the IM/AM of devices should depend on a specific operating point. Even if there is no any change within the power electronic device, any external change would also change the operating point and further the IM/AM result of the device. This makes the IM/AM measurement quite time-consuming under variable operating conditions, especially when we want to exhaustively study the system's stability boundary or fault-set searching for various working conditions. Therefore, it becomes a great challenging to develop proper, fast algorithms for IM/AM prediction (or IM/AM estimation) under both black-box and variableoperating-point conditions.
So far, there is not any work on the IM/AM prediction under fully black-box conditions, to the best knowledge of the authors. Relevant to this target, recently Amin and Molinas proposed a gray-box method to estimate controller parameters based on the already-known control structure [23]. Zhang et al. and Gong et al. realized admittance estimation based on neural network [24] and multidimensional interpolation [25], respectively. Liu et al. proposed a novel ''nearly-decoupled AM'' (NDAM) model, in which the sequence-domain AM of VSC by using the phase-locking loop (PLL) synchronization can be obtained under variable operating conditions [26]. The major limitation of the latter method is that it is still a gray-box method, on the assumption that the control structure of VSC should include PLL. Recently, Liu et al. [27] further proposed a sequence-domain frequency-coupled impedance model (FCIM) with the corresponding stability criterion, and well studied the problem of subsynchronous oscillation caused by interaction between power electronic equipment and weak AC grids.
This article will extend the NDAM in the sequence domain [26] to the dq domain and make some great improvements by applying it to a VSC under an arbitrary control structure and completely black-box conditions. The whole prediction algorithm is, therefore, referred to as fully-decoupled IM (FDIM). It is further improved by a transfer-function curve fitting algorithm. All these are expected to be helpful for system stability analysis and applicable for various practical problems, in particular, when a large number of repetitive IM measurements under variable operating conditions are needed.
The following contents are organized as follows. Section II introduces the detailed algorithm of the fully-decoupled IM. In Section III, a curve fitting method of the transfer function is proposed. In Section IV, the accuracy of our IM prediction under different control strategies is verified by simulations in MATLAB/SIMULINK. Section V is dedicated to validating the FDIM method through hardware-in-the-loop (HIL) experiments of direct-drive wind turbine generators under different operating conditions. Finally, Section VI is devoted to conclusions.

II. FDIM FOR GRID-TIED VSC
A. TYPICAL STRUCTURE OF GRID-TIED VSC Figure 1(a) schematically shows a grid-tied VSC with an unknown control system under black-box conditions. The symbols L f and R represent the interfacing inductance and resistance of the VSC, respectively, C denotes the capacitor at the DC side, P in is the input active power at the DC side, e abc is the internal potential of the VSC, and V tabc and i abc are the three-phase voltages and currents at the point of common coupling (PCC). We use V t and e t to represent the voltage amplitude of V tabc and e abc , respectively. The PWM is for the pulse width modulation technique. The key objective of the IM method is to study the terminal characteristics of devices, i.e., the linearized relation between the terminal voltages and currents at the PCC, which can be further used for small-signal stability analysis.
The control methods denoted by ''control system'' in Fig. 1(a) can be diverse. In particular, a typical control structure of the VSC is shown in Fig. 1(b); it includes the vector control based on the Park transformation of the phaselocking loop, which is a typical synchronization technique in renewable energy integration. The outer control loops include the terminal voltage control (TVC) and the DC-link voltage  control (DVC). The inner control loop is the alternating current control (ACC). θ PLL is the phase output of the PLL. V tref and V dcref are the voltage references, and k p1−4 and k i1−4 represent the corresponding proportional and integral coefficients of the four PI controllers in the DVC, TVC, ACC, and PLL, respectively.

B. IM OF VSC IN THE DQ DOMAIN
Transformation of the three-phase variables into the dq reference frame, namely, from AC signals to DC ones, is an useful technique for the control design of converters. No matter what kind of control system the VSC contains in Fig. 1(a), the inverter system usually has two dq frames: one is the system dq frame (denoted by d s -q s ), and the other is the controller dq frame (denoted by d c -q c ); their relation is illustrated in Fig. 2, where the angle mismatch between the two dq frames in dynamical process is denoted by θ. In the steady state, the two coordinates are identical, i.e., θ = 0. However, if an external disturbance is applied, the two coordinate systems will no longer coincide. VOLUME 10, 2022 The conversion of the voltage and current vectors in the two coordinate systems is determined by the matrix T θ , and are the steady-state values of (V d , V q , I d , I q ) of the two coordinate systems, respectively. Under the steady-state situation, Considering a small perturbation θ ( θ 1), we have and, similarly where the small perturbations of voltage and current at the PCC within the system and the controller dq frames are denoted by ( V s d , V s q , I s d , I s q ) and ( V c d , V c q , I c d , I c q ), respectively. Therefore, no matter what kind of control system is used in the VSC, the conversion relation between the small perturbations of the control dq frame and the system dq frame is determined by (5) and (6). The steady-state variables V d , V q , I d and I q at the PCC can be obtained directly by measurements. Since the dq impedance model in this article is established in the dq coordinate system with the angle of the terminal voltage as the reference angle of the Park transformation, in the steady state, V d = V t and V q = 0.
Before developing our black-box impedance model, let us take a closer look at some typical control schemes of the VSC, which have been widely proposed in the literature. Generally we can divide the control structure into two parts, namely the synchronization control part and the other auxiliary control part. Among them, the role of the synchronization control is to generate a reference angle θ in the control system, and the other part is to design a suitable control loop. Generally there are two major divisions for synchronization, namely the grid-following control and the grid-forming control [15], [28]. The grid-following control is to follow the terminal voltage phase angle through the phase-locking loop, as we have seen in Fig. 1(b). In contrast, the grid-forming control provides the phase information mainly based on the active  [29]. We will see that this universality in their control structures could help us to perform IM prediction conveniently.
After performing small-signal derivations on the PLL in Fig. 1(b) and the four different controllers in Fig. 3, we find that θ can always be written in a unified form, These eight variables including f din , f did , f qin , f qid , f dvn , f dvd , f qvn , and f qvd are intermediate variables generated by linearization. They are all implicit functions of the operating points, that is, the implicit function formula does not change with the change of the operating points. For the detailed derivations, see Appendix A. As usually V d = V t and V q = 0, only three variables of the operating points: V t , I d , and I q are needed.
We emphasize that only based on these general relations, can we get deeply into the internal structure of impedance for any black-box controlled VSC. Although we cannot exhaustively study all existing different control schemes in the literature, the above relation is expected to be generally applicable. It is also worth noting that although under the VSG control in Fig. 3(d), the output angle θ is the internal voltage angle of the VSC (e t ) instead of the terminal voltage angle at the PCC (V t ), θ still represents the angle mismatch between the system and the controller dq frames in the dynamical process, and does not affect our impedance modeling.
For the PLL synchronization control, the theoretical impedance expression of VSC in the dq domain has already been deduced in Refs. [16] and [18]. After linearizing the controllers, all the intermediate variables except θ are a linear function of the operating points. Namely, f dθ , and f qθ have a linear relation with the operating points (V t , I d , and I q ). Again this relation in (8) is believed as workable for any general black-box VSC. Substituting (7) into (8), we obtain Based on these equations, one can find that all g i1−4 and g v1−4 have a polynomial, nonlinear relation with the operating points.
Hence according to the impedance definition of the VSC, we can explicitly express each element of the impedance matrix Z as where gv = g v1 g v4 − g v2 g v3 .

C. DECOMPOSITION OF IM IN THE FDIM
Based on the above theoretical analysis, one can find that the four elements of the IM in (12) are complicated nonlinear functions of system parameters, operating points, and disturbance frequency ω (as s = jω). In contrast, all f 's in (7) and (8) are linear functions of the operating points, and all g i1−4 and g v1−4 in (10) are their polynomial functions. Based on these observations, as the first step, we like to reshape the IM in (12) by decoupling the operating points.
As g i1−4 and g v1−4 can be written as polynomials of the operating points (I d , I q , and V t ), including one constant term, three first-order terms (I d , I q , and V t ), six second-order terms and I q V t ), and many third-order terms . . ), etc, and importantly in the per-unit system, I d ≤ 1, I q ≤ 1, and V t ≤ 1, which indicates that the lower-order terms should have more contribution. Then g i1−4 and g v1−4 can be written as based on their different orders, where all coefficients a kn and b kn for k = 1, 2, 3, 4 and n = 1, 2, 3, 4, . . . depend on the system parameters and disturbance frequency only, and they are irrelevant with the operating-point information.
Equations (13) can be rewritten in the compact vector forms: and where a k and b k denote the induced system parameters. They contain the combined information of the system parameters and disturbance frequency, and x represents the operating point vector separately. Clearly now the operating-point information has been completely decoupled. Further substituting (14) and (15) into (12) gives with the matrices A 0 and A ij (for i, j = 1, 2) expressed by Hence the four IM elements (Z dd , Z dq , Z qd , and Z qq ) in (16) have been explicitly divided into the operating point vector VOLUME 10, 2022 x, and the matrices A 0 and A ij (for i, j = 1, 2), which depend on the induced device parameters a 1−4 and b 1−4 only. Now the operating point in the IM prediction has been fully decoupled, and this point is believed as a great improvement, in comparison to the NDAM method in [26]. Based on (16), we can get IM for any operating condition immediately, under the condition that we have obtained all the induced system parameters of the device, a k and b k , which should be further achieved by parameter identification from multiple IM measurements.

D. IDENTIFICATION METHOD IN THE FDIM
For the aim of parameter identification of all induced system parameters (a k and b k ), first we should obtain their possible restriction conditions, based on the already-known information of the operating point vector x and the IM from several IM measurements. Inputting (14) into (9), we have the following homogeneous equations According to the definition of the dq-domain impedance in (11), we have Substituting (19) into (18) gives Since both I s d and I s q represent external disturbances in the small-signal perturbation tests, they should be independent. Thus the relations in (20) hold true only if their coefficients 1−4 are all zero, i.e., 1−4 = 0. Then inputting (14) into (21) gives Now we have obtained the restricted relations between the impedances (Z dd , Z dq , Z qd , and Z qq ), the operating point vectors x, and the system parameters (a k and b k ). These equations in (22) would become the basis for the parameter identification of a k and b k .
Before doing this, we like to express Eqs. (22) in a more compact form. Taking the first and third ones in Eqs. (22) into one group, and the second and forth ones into the other group, we have and Observing Eqs. (23) and (25), one can see that ρ 1 and ρ 2 are the two solutions of the same linear equation m T 1 ρ = 0. After solving (23), we can obtain a 1 , a 3 , and b 1−4 . Similarly, we can further solve (24) and obtain a 2 and a 4 . Therefore, in principle we can obtain all these necessary parameters a 1−4 and b 1−4 , and accomplish the IM prediction in (16).

E. DETAILED ALGORITHM OF THE FDIM
The detailed algorithm of the FDIM is listed as follows.
i) Since the system control topology and system parameters are all unknown, and the real combination of the operating points (I d , I q , V t ) in x is undecided under the black-box conditions, we first need to fix the form of x, which usually starts from a constant term, first-order terms (I d , I q , and V t ), and second-order terms (I 2 d , I 2 q , V 2 t , I d I q , I d V t , and I q V t ), based on the fact that lower-order terms usually have a stronger impact, i.e., clearly the size of x is 10 × 1, the same as that of a 1−4 and b 1−4 . ii) Suppose that the dq-domain impedance elements (Z dd , Z dq , Z qd , and Z qq ), the terminal currents (I d and I q ), and the terminal voltage (V t ) of different operating points have been already obtained by N perturbation tests in the impedance measurements. Substituting them into (23) and (24), we get the following linear algebraic equations: and  (27), representing either ρ 1 or ρ 2 , and the size of the vector ρ is 30 × 1.
iii) Next we solve the homogeneous linear equations (27). Obviously, ρ 1 and ρ 2 are linearly independent, which means that the dimension of the vector space of solutions of (27) is at least 2. In view of this, the rank of M 1 must be no more than 28 (30 − 2 = 28). Furthermore, the vector space of solutions of (27), i.e., ρ 1 and ρ 2 can be obtained only if the rank of M 1 reaches the maximum, and the total number of perturbation experiments, N , should be large enough. For convenience, we choose N = 28 different groups of operating points accompanying with their corresponding impedance matrices from the IM measurement. Hence we obtain b 1−4 , a 1 , and a 3 from (27), and then substitute them into (28) to get a 2 and a 4 . iv) Finally, based on these identified induced system parameters of a 1−4 and b 1−4 , we can easily calculate the IM in (16), for any varying operating point.
It is worth noting that since the assumed composition of x in (26) may be too sufficient, the addition of more terms would add more new free variables to the solution vectors ρ. Then the number of the basic solution set of (27) may be more than two. In fact, under certain preconditions, any two linearly independent solutions of (27) are equivalent to the model parameter vectors ρ 1 and ρ 2 in calculating the FDIM. For the detailed proof and precondition, see Appendix B. If the precondition is not met, it means that the selection of the operating points are inappropriate, and the known operating points and their corresponding impedances need to be increased. If the basic solution set of (27) is less than two, it means that the assumed element of x is insufficient. Then we need to increase the combination of operating points in x in (26), such as the third-order terms (including I 3 d , I 3 q , V 3 t , I 2 d I q , . . .). However, in all the studied cases in the paper, this does not happen.
In addition, in order to ensure the accuracy of the FDIM, any two sets of operating points used for the parameter identification should be significantly different from each other, and the operating points near operating boundary should be selected as possible as we can.

III. FURTHER IMPROVEMENT: CURVE-FITTING BASED FDIM
Basically, we have already accomplished the IM prediction by applying the above algorithm on any fixed disturbance frequency and step by step on any frequency region which we are interested in; this manner is similar to the impedance measurement by frequency sweeping. Actually, we can also realize this in a more efficient manner by using a transfer function curve fitting technique. Namely we can estimate the IM within the whole frequency region, once and for all, based on several already-known IM data under certain frequencies only, if we can predict the specific transfer function curve form. For a comparison, we call the former method as pointby-point FDIM, whereas the latter one as curve-fitting based FDIM.
Below let us introduce the curve-fitting based FDIM in detail. First, we assume that the highest power of s in the numerator and denominator of the transfer function is n, and write the transfer function as f (s) = p 0 + p 1 s + · · · + p n s n q 0 + q 1 s + · · · + q n s n (29) where both p n = 0 and q n = 0 are not allowed. Express (29) in the form of the homogeneous linear equation: where Then suppose that we have S group data of amplitude frequency response under different disturbance frequency ω's, namely, s (s = jω) and f (s) in (30) are different for different ω's. Since there are 2n + 2 unknown parameters: p 0−n , q 0−n in (30), we may choose S = 2n + 2 to have the square matrix forms of K 1 , K 2 , H 1 , and H 2 , and obtain where the specific expressions of K 1 , K 2 , H 1 , and H 2 are presented in detail in Appendix C.
Since K 1 is a Vandermonde matrix and the selected s's in K 1 are independent, the determinant of K 1 is not equal to 0 and hence K 1 is invertible. Similarly, K 2 , H 1 , and H 2 are all invertible.
Rearranging (32), we have where Here R is a complex matrix. Denote R r as the real part and R i the imaginary part of R. As q is a real vector, we obtain Therefore, q can be obtained by calculating the eigenvector corresponding to R r with an eigenvalue of 1. Furthermore, p can be solved in any equation in (32), and hence the parameters p and q in the transfer function in (29) can be fully obtained. It should be noted that in the calculation, as we do not know the specific form of the transfer function, we usually VOLUME 10, 2022 start from n = 1, increase n gradually, and stop the calculation until the eigenvector corresponding to the eigenvalue 1 of R r in (35) is unique.
Wide numerical simulation results have shown that the curve-fitting based FDIM can indeed greatly improve the calculation efficiency, as it deals with the estimation of the transfer function in a global manner. In contrast, the point-bypoint FDIM costs longer calculation time. Compared to the classical frequency-domain transfer function curve fitting, which is usually approximated by a rational function of factorizing numerator and denominator; see, e.g., [30] and [31], our curve-fitting based FDIM method is superb, when the order of the transfer function is low.

IV. SIMULATION AND VERIFICATION A. SYSTEM CONFIGURATION
Without losing generality, the system parameters in Fig. 1 To test the FDIM algorithm, we have studied several major control forms, such as the PLL combined with ACC, DVC, and TVC, the PLL combined with ACC and power flow control, the droop control, the virtual synchronous generators control, etc., as shown in Figs. 1 and 3. For the black-box study, we do not know any information about the system parameter and the control forms in advance. The whole system is established and simulated in MATLAB/SIMULINK. Totally 28 sets of operating points are necessary for the parameter identification and the last operating point (set r = 29 in the bottom row) is used for verification; they are all listed in Table 1. It is worth noting that in the simulations, the information of any operating point of V t , P (P = V t I d ), and Q (Q = −V t I q ) is recorded instead of that of I d , I q , and V t , based on the convention of electric power engineering. As V d = V t and V q = 0, only three variables are independent and both (V t , P, Q) and (I d , I q , V t ) can be used.
First, the IM of the VSC under 28 sets of different operating points are obtained by the frequency-scan perturbation tests, the same as the usual IM measurement. These data serve as our basis for the parameter identification of a 1−4 and b 1−4 . The perturbation frequencies f set are obtained by randomly selecting from 1 Hz to 1000 Hz. Second, the induced parameters, a 1−4 and b 1−4 , at each disturbance frequency are identified, as described in Section II. Then by using the point-by-point FDIM, the impedance of each disturbance frequency under the operating point of the 29th set (r = 29 in TABLE 1. Operating points used to identification (r = 1 − 28) and verification (r = 29). Table 1) is obtained. Finally, by using a faster curve-fitting based FDIM, the IM under the operating point of the 29th set can be predicted, as described in Section III. To clearly exhibit the whole calculation process, a detailed flow chart of the curve-fitting based FDIM prediction method is given in Appendix D.
In simulations, we found that the running time of the FDIM depends on the time consumed by impedance measurement. If the test operating point is selected properly, for a single frequency point, the running time of the FDIM is the time of 28 groups of impedance measurements.
As one typical example, the results of the frequency-scan IM (FSIM) (solid line), point-by-point FDIM (cross points), and curve-fitting based FDIM (dashed line) are compared in Fig. 4, where clearly both estimations from the pointby-point based FDIM and the curve-fitting based FDIM fit with the FSIM from measurement almost perfectly. This figure well verifies the accuracy of the FDIM. Note that due to space constraints, only the result of Z dq is presented here and all other elements show similar excellent results.

C. PREDICTION OF IM UNDER THE ACC, PFC, AND PLL
The control strategy is now changed to the ACC, PFC, and PLL, as shown in Fig. 5. Without losing generality, the system parameters are k pp = 0.0028, k ip = 7, k p3 = 0.3, k i3 = 160, k p4 = 50, and k i4 = 2000 [18].
The procedure to compute the FDIM is the same, and the amplitude and phase response results by using the frequency-scan IM, point-by-point FDIM, and curve-fitting based FDIM are presented in Fig. 6. Again we see that the IM prediction results match with the IM measurement result well.

D. PREDICTION OF IM UNDER THE BASIC DROOP CONTROL
In this section, the basic droop control is adopted, as shown in Fig. 7. Since the small-signal output angle θ's of the PSC,  , here we only take the basic droop control as an example. Without losing generality, we choose the following system parameters: k p = 0.02, k q = 0.1, k pv = 0.6, k iv = 700, k p3 = 0.4, k i3 = 15, and ω 0 = 2π × 50 [32], [33]. The amplitude and phase response comparison results are presented in Fig. 8. Clearly although there is a jump at around 500 Hz in its phase response curve, their matching is still good.

E. PREDICTION OF IM UNDER THE VSG CONTROL
We have also tested our FDIM prediction algorithm by using the virtual synchronous generator control [34], shown in Fig. 9. VSG is a control scheme applied to inverters of distributed power generation units by imitating the behavior of a synchronous generator. Without losing generality, the following system parameters are chosen: D p = 10, D q = 1/15, J = 0.5, and ω 0 = 2π × 50 [35], [36]. The final comparison results are shown in Fig. 10. Again they match well.

A. TEST SYSTEM CONFIGURATION
In order to further verify the accuracy of the FDIM method, we apply the algorithm to a larger power system consisting VOLUME 10, 2022  of direct-drive wind turbine generators (WTGs). Different with the previous tests in the article, here we consider the case of non-constant active power on the DC-link capacitor, which is more reasonable in practical scenarios. We carry out hardware-in-the-loop (HIL) experiments based on the RTLAB. The digital signal processing (DSP) used is a WTG simulation produced by GoldWind company. In the experiments, the RTLAB simulates the dynamics of power electronics and other electromagnetic components, and outputs analog signal to the DSP control system to realize the control algorithm and modulation, and then the DSP feedbacks the PWM pulse signal to the RTLAB. In addition, because the controller is an integrated dual-WTG controller, the operating point of the second WTG (WTG2) should be consistent with that of the first WTG (WTG1) and the test is only on the WTG1 to verify the accuracy of the fully decoupled impedance model.
The whole experimental platform and the schematic diagram are illustrated in Fig. 11. Without losing generality, the system parameters are all listed in Table 2.

B. PREDICTION RESULTS OF THE FDIM
The procedure to compute the FDIM is completely consistent with the above simulation. Due to the voltage amplitude and power limitation in the experimental setup, again 28 sets of operating points for identification were re-selected in the   Table 3. In order to verify the accuracy of FDIM, three sets of operating points, shown in Table 4, are arbitrarily chosen and used for the test. Since we do not know the internal structure of the controller, we treat it as a black box. Therefore, we like to verify whether the predicted impedance at the given operating points is consistent with the measured impedance. The normalized error ε is defined as    Table 4.

experiment, as shown in
where N f is the number of disturbance frequencies and the final ε values are given in Table 4. We find that the error between the predicted FDIM and the measured FSIM is small; ε < 5%. The detailed comparison results of amplitude and phase response of these three sets of operating points are presented  Table 4. in Figs. 12, 13, and 14, respectively. In these plots, the red solid line represents the impedance value from the measurement, and the blue cross point represents the point-by-point based FDIM. Although the prediction results for the pointby-point FDIM are not as perfect as the simulation results above, they still exhibit a sufficient accuracy. Now due to the high order of the IM, the curve-fitting based FDIM does not work.

VI. CONCLUSION
As is well-known, the impedance of a white-box device can be theoretically derived, but that of a black-box device can only be obtained from measurements. Until now the IM measurement is generally believed as a dominant workable method in electric power engineering. Since in the IM measurement, the impedance is highly correlated with operating points, it becomes difficult to quickly obtain impedances for a large number of operating points. Especially, as any impedance measurement is possible only for stable working conditions and the generalized Nyquist criterion can only be used when the impedance amplitude-frequency response is known, it is difficult to predict its small-signal stability before the VSC's connection to the grid, which might be either stable or unstable. All these might restrict the IM measurement application in many practical problems.
The main purpose of this paper is to solve this challenging problem for a fast impedance prediction under both the black-box and variable operating conditions, by relying on a completely decoupled parameter identification algorithm and several already-known impedance measurement data only. In particular, based on this IM prediction algorithm, the impedance becomes apparent immediately, on the basis of only 28 independent IM measurements. The fully-decoupled impedance prediction algorithm has been found accurate and efficient, and it has been widely verified by the simulations VOLUME 10, 2022 and experiments. When the active power into the DC capacitor side is not constant in a more realistic environment, such as in the case of direct-drive wind turbines, it is found that the algorithm also works well. In addition, as the transfer function curve fitting algorithm, working as an auxiliary tool, is useful only when the order of the system is low, the original pointby-point based FDIM should be used in experiments.
Finally, to make the relation of our algorithm: FDIM, with all existing major IM methods in the literature clearer, a comparison is presented in Appendix E, with an accompanying Table 5. Therefore, the developed impedance prediction FDIM is expected to be capable of working as a supplemental technique of the IM measurement. It is useful for studies on relevant problems in various practical engineering environments, such as IM prediction for unknown operating conditions, stability boundary determination, and small-signal fault discrimination under a large amount of working conditions.

APPENDIX A SMALL SIGNAL MODELING OF SYNCHRONIZATION CONTROL PARTS IN A UNIFIED FORM A. PHASE LOCKING LOOP (PLL)
For the PLL synchronization in Fig. 1(b), recall that the PLL output angle in the small signal is Substituting (37) into (5), we obtain and rewrite it in a unified form, where G PLLn = k p4 + k i4 /s and G PLLd = s + V d (k p4 + k i4 /s).

B. POWER-SYNCHRONIZATION CONTROL OR BASIC DROOP CONTROL
For the power-synchronization control in Fig. 3(a) and the basic droop control in Fig. 3(b), their θ's are the same after the small-signal derivations. Namely,

C. DROOP CONTROL WITH LOW-PASS FILTERS
Under the droop control with the low-pass filters in Fig. 3 Table 4.
Therefore, for all these typical synchronization controls, we can have a general relation between the synchronization phase mismatch θ and the small-signal variables of terminal voltage and current, as shown in (7).

APPENDIX B PROOF AND PRECONDITION OF EQUIVALENCE RELATION BETWEEN THE MODEL PARAMETER VECTORS AND THE SOLUTIONS OF EQUATIONS
Define vectors α and β They can be expressed as α = k 1 η 1 + k 2 η 2 +· · · + k n η n and β = k 1 η 1 + k 2 η 2 +· · · + k n η n , where k i and k i , i = 1, 2, . . . , n are constants and η 1 , η 2 , . . . , η n are the fundamental solutions of (27) and (28), they are linearly independent obviously. Let Then we obtain Substituting (43) and (44) into (17) gives Flow chart for the point-by-point based FDIM prediction algorithm, mainly consisting of the impedance prediction part (left) and the curve fitting part (right).
Take the calculation of Z dd as an example. After substituting (45) into (16), we have To make Z dd independent of the coefficients k 1 , k 2 , . . . , k n and k 1 , k 2 , . . . , k n , only when Then, the precondition is that for η 1 , η 2 , . . . , η n , any two of the solutions as α and β are substituted into the Z dd calculated by (22) and (21) are equal.

APPENDIX C
The expressions of K 1 , K 2 , H 1 , and H 2 in (32) are as follows:

APPENDIX D FLOW CHART FOR THE FDIM
A flow chart for illustrating the major steps of the point-bypoint based FDIM is shown in Fig. 15.

APPENDIX E COMPARISON OF DIFFERENT IMPEDANCE ACQUISITION METHODS
In order to make the properties of FDIM clearer, the advantages and disadvantages with five other major impedance identification methods are presented in Tab. 5, including the FSIM, artificial neural network-based impedance estimation (ANN-based IE) [24], multidimensional interpolation based impedance estimation method (MI-based IE) [25] gray-boxbased controller parameter estimation (GCPE) [23], and the nearly decoupled impedance model (NDIM) [26]. The FSIM is one of experiment-driven and black-box impedance measurement methods without systematic error. It is simple and accurate. What is more, it has wide and mature applications in engineering practice. However, it cannot be used to predict impedance values. The ANN-based IE is one of data-driven and black-box impedance prediction methods with systematic error, but it lacks physical meanings and is often greatly influenced by the training data. The MI-based IE is one of data-driven and black-box impedance prediction methods with systematic error as ANN-based IE. It also lacks physical meanings and the predicted values outside the interpolation interval are inaccurate. The GCPE is one of model-driven and gray-box impedance identification methods without systematic error, but the control structures must be known in advance. The NDIM is one of model-driven and black-box impedance prediction methods without systematic error. However, it is valid under the PLL synchronization control structure, with the dynamic behaviors of the DC side completely neglected. In contrast, the proposed FDIM in this paper is a model-driven and black-box impedance method without systematic error. As we have seen, besides the VSC and VSG, it is also valid for wind turbines and other control structures. Hence it is unique for electric power engineering applications. His research interests include the ways to analyze the stability of power system with renewable energy generation, and in particular impedance modeling and analysis of the voltage-source converter, and stability and control of power systems with renewable energy generation, such as impedance model, harmonic linearization, and small-signal stability analysis. He was a Postdoctoral Researcher with the National University of Singapore, Singapore, and the University of Toronto, Canada, from 2001 to 2006. He worked as a Full Professor with the Wuhan Institute of Physics and Mathematics, and the Chinese Academy of Sciences, from 2006 to 2015. In 2015, he joined the State Key Laboratory of Advanced Electromagnetic Engineering and Technology, School of Electrical and Electronic Engineering, Huazhong University of Science and Technology. He has long been engaged in the study of nonlinear dynamics theory of complex systems in multiple directions, such as coupled nonlinear systems, chaos synchronization and control, pattern formation, and complex network dynamics. He has published more than 100 SCI articles in internationally peer-reviewed journals. His current research interests include renewable energy integration stability, power-electronic-based power system dynamics, and nonlinear analysis of power systems.