DC-Side Impedance Estimation of a Modular Multilevel Converter Through System Identification of a Partially Black-Boxed Control System

The stability of a power electronicssystem can be assessed by means of the impedance-based stability criterion. Impedance modeling is a useful tool to analyze the effect of different circuit parameters and control schemes on the behavior of a converter. Modeling the input impedance of a power electronics converter is often successful when having full knowledge of the converter topology, the circuit parameters, and the parameters and implementation of the control system. However, due to the proprietary nature of voltage source converter-based high voltage direct current systems, their exact control structure is often concealed. This complicates the calculation of the impedance of a modular multilevel converter, known for its complex internal dynamics. This paper proposes a method to estimate the impedance of a modular multilevel converter with partially black-boxed converter control. A discussion on partitioning the control system into open and closed parts is made, and the results are verified with simulations in time and frequency domains.

V OLTAGE source converter (VSC)-based high voltage direct current (HVdc) transmission allows low-loss transfer of bulk power over long distances and interconnection of asynchronous ac grids [1]. To date, the converter stations of nearly all the existing VSC-HVdc systems in operation are provided by single vendors. These systems are anticipated to interconnect in the future and form a supergrid [2]. It is fair to assume that in such an event, converter stations designed by different vendors would need to be interconnected. The behavior of these converters is mainly impacted by their circuit design and control systems, the specifics of which are almost always proprietary and protected by intellectual property (IP) rights [3].
HVdc systems are prone to control interactions requiring careful investigation into causes and measures of mitigation [4], [5]. In multi-vendor HVdc systems, the legal difficulties for vendors to exchange information and expertise complicates the analysis into the root causes behind the control interoperability issues [6]. This is troublesome for the transmission system operator as the owner of the HVdc converters and the responsible party for power delivery and system maintenance. So far, the practice has been to seek assistance from an independent party that may be given access to the proprietary information but is not allowed to disclose any to either vendor. The independent party then has to collaborate closely with the vendors to find and alleviate the issue(s) [7]. This procedure is cumbersome and time-consuming as close collaboration with an individual vendor while coordinating with another runs the risk of unintentional IP disclosure.
Harmonic stability of a power electronics-based system, including VSC-HVdc, can be studied using several tools, one of which is the impedance-based stability criterion [8]. To this end, the system is partitioned into a source and a load subsystem modeled by their input and output impedances, and the system stability is assessed via the ratios of these impedances. For an accurate stability assessment, either accurate measurement of the input impedance of the converter at its terminals or an accurate model of the converter input impedance is required. Impedance modeling is a useful tool that although initially time-consuming, allows computing the impedance under various circuit and control parameters instead of a fixed set of parameters at a single operating point.
Research in impedance modeling of VSC-based HVdc systems and modular multilevel converters (MMCs), which are the state-of-the-art in such technology [9], is abundant [10]- [17]. Nevertheless, without full knowledge of the control system implementation and the main circuit parameters, the impedance cannot be modeled accurately. A few attempts have been made This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ to derive the impedance model without prior knowledge of the control system. In [18], a gray-box method is presented to estimate the assumed unknown control parameters of the control system of a wind energy conversion system, and mitigate the arising instability issues. While the parameter estimation is beneficial, an impedance transfer function estimation does not add value as the stability can also be assessed using the measured impedance values without the need for a transfer function. In [19], a novel transfer learning method based on artificial neural networks is presented with an aim to adapt the impedance of a known VSC system to a different but related system. An accurate model is obtained with an insignificant error, which is capable of estimating the impedance at various operating points. Reference [20] presents a method for parametric identification of the impedance from time-domain data, which mitigates the coupling effect of the grid impedance introduced by impedance modeling in the dq-frame. In all of the above, it is assumed that the control system implementation of the VSC is fully known, but the controller settings are sometimes unspecified.
In industrial HVdc projects, the transmission system operator has full knowledge of the converter topology and parameters and limited or no knowledge about the converter control system implementation and its parameters. In particular, in multi-vendor installations, the suppliers are reluctant to disclose their entire control system structures as this information essentially gives them the competitive edge over other involved suppliers. Nevertheless, it is conceivable that some (high-level) outer control loops of the converter control systems can be designed or treated as open blocks [21].
A proposed approach to study control interoperability is to divide the converter control system into open and closed parts [22], [23]. Thereafter, the proprietary information, which includes vendor IP, is placed within the closed segment, and more generic controllers are released as open blocks. This would provide a degree of freedom for the system operator to run control interoperability studies while reducing the risk of vendor IP disclosure. In general, the innermost control blocks (modulation techniques, current controllers, etc.) contain the highest degree of innovation of which the vendors are most protective. These blocks can remain within the closed-source (black-box) part of the control system, whereas the outermost high-level control loops (dc voltage, power control, etc.) intended mainly as system-related control functions may possibly be disclosed, designed, or tuned independently.
Research on system identification of three-phase VSCs for control and stability assessment purposes is scarce. This paper attempts to fill this gap. Compared to the state-of-the-art in the literature, the novelty and original contributions of this work are r applied subspace-based system identification on a partially black-box converter control system, r dc-side impedance model estimation for a modular multilevel converter with black-box control components, r investigation into the effect of the outer loop (open blocks) control parameters on the dc-side impedance without the need for additional frequency sweeps, r and verification of the proposed method by means of simulations.
A method is proposed to estimate an impedance model for a modular multilevel converter in which a part of the control system (structure as well as parameters) is black-boxed. A combination of input/output measurement and system identification is used to estimate the transfer functions of these control stages and incorporating them within the impedance model of the known system. In doing so, information about the black-box control blocks is extracted and used in tuning the open control system with the aim of tackling interoperability issues. This has the benefit of calculating the converter impedance with a high degree of accuracy while keeping the intellectual property and the exact implementation of the black-box control blocks intact. Finally, in a case study, the usefulness of the proposed approach in handling harmonic stability issues in a back-to-back MMC-based VSC-HVdc system is demonstrated.
The paper is organized as follows: Section II describes the converter dynamics and discusses conventional control design and how the converter control is partitioned. Section III explains the methods of system identification for the black-box control system considered in this work. Section IV describes the derivation of the estimated dc-side impedance model and discusses the effect of outer loop design on the impedance of a modular multilevel converter with black-box control components. Section V presents and discusses the results of the study, and finally, Section VI draws conclusions based on the results.

II. SYSTEM MODEL
The MMC for ac/dc conversion comprises three phase legs with two arms within each leg. The converter arms consist of a (large) number of cascaded submodules and an arm inductance designed to prevent high transient currents and reduce harmonic distortion. A half-bridge submodule consists of two semiconductor switches with anti-parallel freewheeling diodes and a capacitor C SM . Fig. 1 shows the topology of an MMC in which R and R g denote the arm resistance and grid resistance, respectively, L and L g denote the arm inductance and grid inductance, respectively, and v ul,j , j ∈ {a, b, c}, the upper-and lower-arm voltages of the three phases. Some of the system parameters, if not given, can be estimated. For instance, the arm inductance dominates the high frequency input impedance [17] and thus, its value can be approximated by injecting high frequency perturbations. The arm capacitance can be approximated using either the second-order harmonic ripple of the capacitor voltages [24], or by data driven methods [25]. In this paper, full knowledge of the converter circuit structure and parameters has been assumed, and only parts of the control system have been treated as a black-box. Throughout the text, a per-phase modeling approach is taken where the subscript indicating specific phases is dropped.

A. Converter Dynamics
A continuous-variable dynamic model of the MMC is adopted, which neglects the switching operations and assumes balanced submodule capacitor voltages within an arm [26]. Thus, the dynamics of the MMC arm current can be obtained using Kirchhoff's voltage law as where i u (i l ) is the upper (lower) arm current, i s is the ac-side current, v du (v dl ) is the upper (lower) dc-side voltage, e is the point of common coupling (PCC), and the grid voltage is defined as and f 1 (ω 1 ) as the (angular) fundamental frequency. Assuming v du = v dl = v dc /2, adding and subtracting (1a) and (1b) yields with where v s is the voltage driving ac-side current i s and v c the voltage driving the circulating current i c . The arm voltages can be described as where n u,l , the (upper and lower arm) insertion indices, are the outputs of the control system which are computed by closed-loop control, and v Σ Cu,l are the (upper and lower arm) sum capacitor voltages. Assuming balanced capacitor voltages, the expression describing the dynamics of the sum capacitor voltages can be obtained using the relationship between the instantaneous power and the stored energy in the arm [26] as where N is the number of submodules and v Σ C0 is the initial sum capacitor voltage.

B. Converter Control
MMC control is often implemented in a cascaded manner, with outer control loops setting the reference for the inner loops, as shown in Fig. 2. For instance, it is common practice that inner current controllers are nested within the outer voltage or power control loops [27]. The outputs of the inner control loops are then modulated to generate the switching signals applied to the converter. Due to their effective high-frequency nature in HVdc systems, the modulation and balancing algorithms have little impact on the average behavior of the converters [28] and, consequently, the control interactions which can emerge in multi-vendor converter installations. Therefore, these stages have been disregarded in the modeling in this work. In addition, the MMC topology requires an additional current controller to suppress the second-order harmonic circulating current, which has a strong impact on the MMC impedance. A phase-locked loop (PLL) estimates the angle of the voltage at the PCC, used in Park and inverse Park transformations in the control system.
The effect of different control schemes on the small-signal stability of the modular multilevel converter has been a topic of extensive research [29]- [31]. The dynamics of a few typical outer control loops are described below.

1) Dc-Bus Voltage Control:
The dc-bus voltage is controlled via a proportional-integral (PI) controller. The controller sets the current reference i sd to where k p,v and k i,v denote the proportional and integral gains of the controller, respectively.

2) Ac-Bus Voltage Control:
The magnitude of the ac-bus voltage can be controlled by reactive power compensation. The controller sets the q-component current reference to

3) Active Power Control:
The active power controller similarly sets the d-component of the ac-side current reference. Accordingly the dynamics are defined as where k p,a and k i,a denote the proportional and integral gains of the controller, respectively, and P is calculated via

4) Reactive Power Control:
The reactive power controller sets the current reference i sq to where Q is calculated as and k p,r and k i,r denote the proportional and integral gains of the controller.

C. Control System Partitioning
As discussed above, each converter station includes multiple control loops, which are strongly coupled, making the HVdc system susceptible to control interactions [32], [33]. Moreover, multi-vendor control interoperability is also of concern given the fact that many existing single-vendor HVdc installations may need to be extended to form multi-vendor systems in the future [6]. Methodologies to identify and mitigate control interactions have been developed and investigated in the literature [18], [30], [34]. Analyzing the stability of internal converter control loops in the control design stage serves a useful purpose as it attempts to ensure that harmonic instabilities do not occur due to incorrect tuning of the inner control loops [35], [36]. Nevertheless, harmonic stability can mainly be evaluated analytically if full knowledge of the control system is rendered. In the context of stability assessment, impedance modeling is more efficient than impedance measurement as a model provides the opportunity to directly assess the impact of different system and control parameters on the impedance without the need for further frequency scans. Therefore, system identification of the black-box controllers proves relevant and useful for impedance modeling in cases where parts of the control system are blackboxed.
Analysis of the effects of different control loops on converter stability and impedance is a prerequisite to partition the control system. Several papers have analyzed the impact of different control schemes on the impedance of the MMC [14], [17], [30], [35], [37]. In [37], it is shown that the impact of inner current control loops on the dc-side impedance of an MMC is minimal compared to that of the outer loops. To demonstrate the impedance estimation method in this study, the ac-side current control, the circulating current control loop, and the insertion index calculation block are treated as a multiple-input multipleoutput (MIMO) black-box system, and an attempt has been made to estimate the transfer function of this MIMO system by means of input/output measurement and system identification. Subsequently, the estimated transfer functions of the now assumed proprietary part of the control system have been incorporated in the impedance model of the converter seen from the dc-side. The method can as well be applied to other black-box control systems than those considered in this work.

D. Assumptions and Approach
In impedance measurement of VSCs, it is standard practice to superimpose a perturbation frequency on the (dc or ac) terminal voltage or current and measure the response of the applied perturbation under frequency sweeps of the perturbation source [38]. The MMC is a nonlinear system with complex dynamics. A perturbation frequency superimposed on any converter variable induces oscillations at the same or different frequencies on other converter and control variables. E.g., applying a perturbation frequency f p on the dc-bus of an MMC induces harmonics at f = {mf p ± nf 1 }, m ∈ N, n ∈ Z on converter and control variables.
If a stiff three-phase ac grid is assumed, the dynamics of the ac grid and the phase-locked loop do not affect the dc-side converter impedance [17]. However, in case of a weak grid, the dynamics of the PLL become relevant. The design of PLL is also often proprietary. Here, it is assumed that the PLL is black-boxed and only the angle estimate of the PCC voltage is available. With these premises, the inputs of the assumed black-box control system are the ac-side and circulating currents in the synchronous reference frame, i s,dq and i c,dq , and their references, and the outputs are the computed insertion indices n u,l , respectively, see Fig. 2. In impedance modeling in the frequency domain, it is sufficient to consider the variables in a single arm, e.g., the upper-arm of phase a, because of the symmetry of the MMC topology. Thus, the output of the MIMO black-box system is defined as n ua .
As an example, Fig. 3 shows the ac-side currents i s,dq , the circulating currents i c,dq , the upper arm phase a insertion index n ua , and their harmonic spectra when a perturbation frequency of 106 Hz is superimposed on the dc-side voltage. As seen in the figure, the major perturbation frequency components of the currents i s,dq and i c,dq are at f p , and the most significant perturbation frequency components of n ua are at Picking a handful of these frequency components provides an impedance model with good accuracy [17], [30]. Thus, in this work, the black-box system inputs are considered and assessed at f p and the outputs at f p ± f 1 and f p ± 2f 1 . Considering the  above, the control system can be partitioned in the frequency domain as shown in Fig. 4.
To identify the transfer function of the MIMO black-box control system, its inputs and outputs are measured while a perturbation frequency is injected on the dc-side of the converter. The arm currents are normally monitored in HVdc systems and the outputs of the assumed black-box system can be obtained directly from (5) if the arm and sum capacitor voltages are measured, or indirectly via (1a) with an error linked to the accuracy of the values of arm inductance and resistance. Thereafter, the measured inputs (currents) are transformed into the synchronous reference frame using Park's transformation which requires the PCC angle estimate. A discrete Fourier transform is performed on the currents in the dq-frame and the insertion index n ua to obtain their frequency responses during impedance measurement at relevant frequencies.
The open control system can be designed vendorindependently. The outputs of the open control system i s,dq and i c,dq , as shown in Fig. 2, are assumed to be accessible and measurable. Thus, as seen in Fig. 4, the frequency components of ΔI s,dq (jω p ) and ΔI c,dq (jω p ), the true inputs to the black-box control system can be calculated. In case the outputs of the open control system are not accessible, the converter can be operated in constant power mode, i.e., constant P and Q , for which I s,dq (jω p ) = I c,dq (jω p ) = 0 for ω p > 0, and thus ΔI s,dq (jω p ) and ΔI c,dq (jω p ) can again be calculated. Having measured the inputs and outputs of the black-box control system in the frequency domain, system identification algorithms can be used to estimate the MIMO transfer function. Section III discusses a method of system identification for the black-box part of the control system.

III. BLACK-BOX SYSTEM IDENTIFICATION
A black-box system provides no physical insight into the relationship between its inputs and outputs. Thus, system identification of black-boxes is inevitably a process of data fitting rather than precise modeling. Assuming that input-output data of a black-box system are available, classical prediction error-based or subspace-based system identification methods can be applied for system estimation [39]. Subspace-based system identification does not require explicit parametrization, does not involve nonlinear optimization, and is non-iterative.
System identification of MIMO systems using samples of frequency response data has been studied in detail since the 90's [40]. Non-iterative subspace-based system identification algorithms can estimate linear state-space models from frequency spectra of inputs and outputs [41], [42]. Since the MMC system is excited by periodic perturbations during the impedance measurement, Fourier transforms of the inputs and the outputs of the black-box system can be obtained with high quality.
The estimation of a rational transfer function G(jω s,k ) ∈ C p×m involves minimizing where U k ∈ C m are the measured inputs, Y k ∈ C p are the measured outputs, and ω s,k ∈ R + are the angular frequencies at which data is measured.
The subspace-based system identification can be applied to minimize (13) [41]. The identification algorithm estimates statespace matrices {A, B, C, D} of the state-space representation of the linear black-box system from which the transfer function of the black-box system can be obtained by where I is the identity matrix of appropriate dimensions, n u , n y and n x denote the number of input, output, and state variables, respectively, and s is the Laplace variable. The subspace-based identification algorithm determines the state-space matrices noniteratively. To this end, a desired integer model order n and an auxiliary integer model order q is assumed for the system such that q > n. Let us define auxiliary input U q ∈ C q×z and output Y q ∈ C q×z matrices as where z denotes the length of the sampled data. An orthogonal projection Π ⊥ q is utilized such that U q Π ⊥ q = 0. The singular value decomposition of the Hankel matrix Y q Π ⊥ q can be defined as where S and V are unitary matrices of appropriate dimensions and Σ is a diagonal matrix. Subsequently, the system matricesÂ andB can be obtained asÂ and I i and 0 i×j denote the i × i identity matrix and i × j zero matrix, respectively. Having estimatedÂ andĈ, matrices B and D can be estimated by solving a least squares problem of the form where F denotes the Frobenius norm. The estimation can be improved by iterating the least squares problem [43]. Having obtained the estimated state-space matrices {A, B, C, D}, the transfer function of the black-box control system can be found by simply using (15).

IV. IMPEDANCE MODELING
To model the dc-side impedance of the converter, the harmonic responses of the converter and the control variables are calculated assuming that a small-signal perturbation voltage with an amplitude of e p and a frequency of f p is superimposed on the dc-bus voltage, see  (5) and (6) leads to a combination of frequency components through addition and subtraction in the frequency domain [16]. The most significant of these harmonic responses appear at a combination of the perturbation frequency f p and multiples of the fundamental frequency mf 1 , i.e., f = f p + mf 1 , m ∈ Z. Thus, to obtain an accurate frequency domain model of the system, converter and control responses need to be evaluated not only at f p , but also the whole range of the affected harmonic spectra. To keep the complexity of the resulting impedance model at bay, m < 3 is assumed in this work. The capitalized letters indicate the Fourier coefficient of the (converter and control) variables.

A. Response of the Phase-Locked Loop
Considering a weak grid on the ac-side of the converter makes the dynamics of the PLL relevant. Thus, its response to the dcside perturbation needs to be calculated. The PCC voltage can be obtained via Kirchhoff's voltage law as The Fourier coefficient of the ac-side current i s is equal to that of the difference of the upper-and lower-arm currents, i.e., In [30], it is shown that differential-mode components of upper-and lower-arm currents will flow to the ac-side. Arm currents have differential-mode components at ω p ± ω 1 and common-mode components at ω p and ω p ± 2ω 1 . Thus, the perturbation frequency component of the ac-side current can be obtained as No perturbation frequency component appears at the (stiff) ac grid voltage v g . Therefore, the Fourier coefficient V g is equal to zero at all perturbation frequencies. Accordingly, perturbation frequency components will appear at the PCC voltage as The PLL is assumed to be a black-box with only its output, the estimated angleθ, available, see Fig. 6. We define a = cos(θ) and b = sin(θ) .
From Park transformation, the dq-components of the PCC voltage can be obtained as Frequency-domain analysis of Clarke transformation indicates the existence of differential-mode components Thus, from (26), the Fourier coefficients E d and E q can be obtained as Perturbation frequency components A(jω 1 ), A(jω p ± jω 1 ), B(jω 1 ), B(jω p ± jω 1 ) and steady-state components of the PCC voltage in the alpha-beta frame can be computed following the perturb-and-measure of quantities a and b, and post-processing them to obtain their Fourier coefficients at the relevant frequencies.

B. Response of the Converter Variables
For the converter variables, i u , v u , and v Σ Cu , the harmonic responses to the applied perturbation are calculated at jω p . The rest of the expressions for other frequencies (including the steady-state frequencies) are given in [17].
1) Converter Currents: The response of the upper-arm phase a current to the applied perturbation can be calculated from (1a). As seen in Fig. 5, v du (the upper dc-side voltage) has a Fourier coefficient V du (jω p ) = e p /2 2 at the perturbation frequency. No perturbation frequency appears at the ac grid voltage, thus the Fourier coefficient of v g is equal to zero at all frequencies.
Considering all the above, the Fourier coefficient of i u at ω p can be calculated as and rewritten as The impact of the grid impedance is thus captured in the differential-mode perturbation components of arm currents as and the Fourier coefficients of i u at jω p ± j2ω 1 can be calculated as Similar to (28), the responses of ac-side current in the dqframe to the applied perturbation can be obtained as I sd (jω p ) = I sα (jω p −jω 1 )A(jω 1 ) + I sβ (jω p −jω 1 )B(jω 1 ) +A(jω p −jω 1 )I sα (jω 1 )+B(jω p −jω 1 )I sβ (jω 1 ) +I sα (jω p +jω 1 )A(jω 1 )+I sβ (jω p +jω 1 )B(jω 1 ) +A(jω p +jω 1 )I sα (jω 1 )+B(jω p +jω 1 )I sβ (jω 1 ) (33a) where I sα (jω p ± jω 1 ) = I s (jω p ± jω 1 ) = 2I u (jω p ± jω 1 ) (34a) The circulating currents have a dominant negative sequence double line frequency component. We define g = cos(−2θ) and h = sin(−2θ) , (35) used in Park transformation of the circulating currents. The responses of the circulating currents in the synchronous reference frame to the applied dc-side perturbation is also calculated similarly as 2) Arm Voltage: From (5), the perturbation frequency component of the arm voltage at jω p can be derived as where the multiplication in the time domain has translated to the convolution of V Σ Cu and N u in the frequency domain. The bar denotes the complex conjugate operation. The steady-state frequency component V Σ Cu (j2ω 1 ) is calculated as while the rest of the steady-state frequency components are given in [17].

3) Sum-Capacitor Voltage:
In a similar fashion, the perturbation frequency component of V Σ Cu at jω p is calculated using (6) as

C. Response of the Control Variables
The response of the outputs of the open control system comprising the dc-bus voltage control, the active power control, and the reactive power control to the applied perturbation is analyzed in the frequency domain.
1) Dc-Bus Voltage Control: From (7), the response of i sd to the applied perturbation can be calculated as 2) Ac-Bus Voltage Control: The perturbation frequency response of i sq can be calculated from (8) as where E d (0) and E q (0), the steady-state values of dqcomponents of the PCC voltage can be measured.

3) Active Power Control:
Similarly, from (9), the response to the applied perturbation can be calculated as In case of an ac stiff grid, the PCC voltage is not affected by the perturbation frequency, therefore, the response of i sd is simplified to

4) Reactive Power Control:
The response of i sq to the applied perturbation can be calculated from (11) by For a stiff ac grid, the response is simplified to

5) Black-Box Control System:
For the black-box system, the structure and parameters are unknown. The responses to the applied perturbation are measured, and the transfer function of the black-box system is estimated using the methods outlined in Section III. The estimated transfer function matrix is then evaluated at the frequencies at which the impedance model needs to be derived.

D. Impedance Derivation
Based on the expressions describing the frequency responses of the converter and control variables to the applied perturbation, a linear system of equations with 28 equations and 28 unknowns of the form is formulated where U p comprises the coefficients of the linear system, V p consists of the constant terms in the equations, and x p contains the system variables The elements of U p and V p for I u , V u , and V Σ Cu at f ∈ {f p ± f 1 , f p ± 2f 1 } are given in [17], while those of the rest of the variables are obtained in (33)-(36b) and (48). Solving the linear system of (50) yields I u (jω p ) from which the dc-side impedance of the converter can be obtained by .
(51) Fig. 7. Flowchart of impedance estimation using subspace system identification.  Fig. 7 shows the flowchart of the process of impedance modeling using this method.

V. RESULTS AND DISCUSSIONS
The estimated impedance is compared to the impedance determined via time-domain simulations of the same system in MATLAB/SIMULINK. The known system and control parameters are given in Tables I and II, respectively, for a single converter. The impedance is determined in the time domain simulation by superimposing a perturbation frequency f p on the dc-bus voltage, see Fig. 5. During the frequency sweep for the impedance measurement, system variables that are inputs and outputs to the black-box control system, i.e., i s,dq , i c,dq and n ua , are recorded and post-processed to obtain their harmonic spectra. Thereafter, subspace-based system identification described in Section III is utilized to estimate the transfer function matrix of the black-box system G. A uniformly spaced set of frequencies between 1 Hz and 1000 Hz with an interval of 1 Hz is considered at which the variables of interest are acquired and recorded. Fig. 8. Bode diagram of the harmonic spectra of the outputs to the black-box system. The solid lines demonstrate the measured values and the circles illustrate the estimated outputs via subspace-based system identification. Fig. 9. Bode diagram of the harmonic spectra of the measured inputs to the black-box system. Fig. 9 shows the Bode diagram of the harmonic spectra of the frequency components of the inputs to the black-box system at f p measured during a frequency sweep. The outputs of the blackbox system are measured in a similar fashion and thereafter, the data is fed to the system identification algorithms for a transfer function matrix to be estimated.

A. Subspace-Based Identification
Subspace system identification estimates a transfer function which maps the outputs shown in Fig. 8 to the inputs shown in Fig. 9. The algorithm is "correct", i.e., given a finite number of samples of data, an exact transfer function can be estimated, albeit possibly with an error in the presence of noise in the system. The estimated transfer function may have a higher order than the real system. Thus, inevitably, it may have some poles and zeros that are not present in the real system. Therefore, incorporating the estimated transfer function in the impedance Fig. 10. The error between the estimated and the measured outputs for various model orders. The estimated output error is below 2% for orders 13 and above. model in the Laplace domain will lead to an impedance which may have unstable poles. To remedy this, instead of incorporating the estimated transfer function in the Laplace domain in the impedance model, the estimated transfer function is evaluated numerically at the frequencies at which an impedance model needs to be derived. The end result will be the impedance computed numerically at arbitrary frequencies within the frequency range of estimation. In case an impedance transfer function in the Laplace domain is needed, a system identification technique based on vector fitting can be applied to the discrete frequency response of the final estimated impedance.
To determine the best model order, the percentage of the error between the estimated and the measured outputs is computed for various model orders. Fig. 10 shows the error in estimated outputs versus model order as As seen in Fig. 10, the estimated outputs can well approximate the real outputs when the black-box system order is assumed n = 13. This implies that an estimated transfer function of order 13 maps the inputs to the outputs of the black-box system with high accuracy. Fig. 8 demonstrates the estimated outputs of the blackbox system (based on the estimated transfer function) versus the actual (measured) outputs of the system with the aforementioned model order.

B. Dc-Side Impedance Model Estimation and Shaping
Having identified the black-box control system, the dc-side impedance of the MMC can be modeled and shaped by the open control system. To this end, the effects of the grid strength and the controller settings of the dc-bus voltage control, the active power control loop, and the reactive power control loop on the dc-side impedance are investigated. Fig. 11 shows the effect of grid weakness, characterized by short-circuit ratio (SCR), on the dc-side impedance. The largest impact is seen around the resonant peaks, but impedance at sub-fundamental frequencies is also affected. In the following, a strong grid is assumed. Fig. 12 shows the impact of different parameters of the dc-bus voltage controller on the impedance. As seen in the figure, the estimated impedance closely follows that of the measured impedance, implying that the that subspace-based system identification has estimated a sufficiently accurate transfer function of the black-box system. The impedance is measured in simulations Fig. 11. Effect of grid strength on the dc-side impedance of an MMC. The solid lines show the measured impedance, and the circles demonstrate the estimated impedance using subspace-based system identification. Fig. 12. Effect of the proportional gain increase of the dc-bus voltage controller on the dc-side impedance of an MMC. The solid lines show the measured impedance, and the circles demonstrate the estimated impedance using subspace-based system identification. of the same system taking into account all perturbation frequency components compared to a few in the estimated analytical model. It can be observed that the gains of the dc-bus voltage controller affect the dc-side impedance mainly in the low-frequency range.
Finally, Fig. 13 shows the effect various power factors on the dc-side impedance under active and reactive power control. As seen in the figure, active and reactive power control increase the dc-side impedance in the low-frequency range and reduce phase compared to dc-bus voltage control. The impedance shape of the active power controlled converter has lower magnitudes around the valleys compared to when dc-bus voltage control was active. The estimated impedance closely tracks the measured impedance, both at low frequencies and around the peaks. The effect of other outer control loops such as

C. Case Study
The system shown in Fig. 14 is considered as case study to investigate harmonic stability. The back-to-back MMC-based VSC-HVdc system is often used to couple two high voltage asynchronous ac networks with the same or different fundamental frequencies. In such a setup, a master converter controls the dc-bus voltage at the PCC, and the slave converter controls the flow of active power, see Fig. 2. Reactive power is controlled separately at each converter station. Table III summarizes the parameters of the system. System identification has been applied to the black-box control system of each converter station. Subsequently, the dc-side impedances of the converters have been estimated.
Thus, the small-signal stability of the system can be evaluated by means of applying the Nyquist stability criterion to the loop gain Z vs (jω p )/Z cs (jω p ). Fig. 16 shows the dc-side impedances of the master and slave MMCs. The impedance curves intersect at three points which are the possible interaction points. With correct tuning of the gains of the open control system, the system is stable as the phase difference at these points are all below 180 • . Incorrect tuning of the open control system, represented here by increasing the proportional gain of the dc-bus voltage controller, can lead to instability. As seen in Fig. 17, the intersection at 70 Hz has now a phase difference of just above 180 • . This is verified via the Nyquist plot in Fig. 18(a) where the loop gain now encircles the critical point (−1 + j0), i.e., the system is unstable. Timedomain simulation of the same system also indicates that the system becomes unstable after the gain change, see Fig. 18(b). The harmonic spectrum of the dc-bus voltage indicates that the frequency of the unstable oscillations corresponds to that predicted by the Bode plot, see Fig. 18(c). This demonstrates the efficacy of the proposed method in estimating the impedance model of an MMC and handling interactions caused by open control system elements. Fig. 17. Bode diagrams of the estimated dc-side impedances of the back-toback converters using subspace-based system identification. The proportional gain of the dc-bus voltage controller has been increased.

VI. CONCLUSION
A method is proposed to obtain the impedance model of an MMC when parts of the control system are black-boxed, i.e., no knowledge about the design specification of those control blocks is available. The application of system identification techniques on a MIMO system for the purpose of impedance modeling is investigated. An impedance model is derived for which the effect of different control schemes and parameters of the open control system, which can be designed vendor independently, is analyzed. The validity of the method is verified by means of simulations. The developed impedance model, which includes both white-box and black-box systems, can be used to handle control interoperability issues and instabilities in HVdc systems.