Hardware-in-the-Loop Methods for Stability Analysis of Multiple Parallel Inverters in Three-Phase AC Systems

—Modern electric distribution systems typically con- tain several feedback-controlled parallel inverters that together form a complex power distribution system. Consequently, a num- ber of issues related to stability arise due to interactions among multiple inverter subsystems. Recent studies have presented methods where the stability and other dynamic characteristics of a paralleled inverter system can be effectively analyzed using impedance measurements. This paper presents implementation techniques for a comprehensive online stability analysis of grid- connected paralleled inverters using power hardware-in-the-loop measurements based on an OPAL-RT real-time simulator. The analysis is based on simultaneous online measurements of current control loop gains of the inverters and the grid impedance, and aggregated terminal admittance measurements of the inverters. The analysis includes the measurement of the inverters’ ag- gregated output impedance, inverters’ loop gains, global minor loop gain, and grid impedance. The presented methods make it possible to rapidly evaluate the system on both global and local level in real time, thereby providing means for online stability monitoring or adaptive control of such systems. Experimental measurements are shown from a high-power energy distribution system recently developed at DNV GL, Arnhem, Netherlands.


I. INTRODUCTION
A N increasing number of applications of AC-distributed power systems have been made possible due to advances in semiconductor technologies and inverter topologies. These applications can be found in several new and emerging fields, including renewable energy generation [1], hybrid and electric vehicles [2], smart grids [3], electric aircrafts [4], and electric ships [5].
AC-distributed power systems are most often dependent on the operation of paralleled inverters. In some applications, such as large-scale photovoltaic or wind power plants, the system may include hundreds or even thousands of inverters operating in parallel for scaling up the total power generation capacity. Fig. 1 shows a typical configuration of multiparalleled inverters connected to a power grid through a H. Alenius (PCC). Usually, the inverters are individually designed based on the stability requirements of the inverter standalone operation. However, the use of multiple inverters in parallel results in dynamic interaction causing system performance degradation or even instability [6], as reported from distribution systems that include photovoltaic plants [7], wind power systems [8], [9], and data centers [10]. Stability issues involved in multi-paralleled grid-connected inverters have been studied extensively in recent years. In [11], a small-signal stability assessment method was presented, where the stability is assessed based on the outermost voltage and frequency droops while neglecting the inner controller dynamics. The work in [12] presented a method in which multiple paralleled inverters were modeled as a multivariable system, where an equivalent inverter described the totality of multiple inverters. The work in [12] was based on the hypothesis that all of the current references of each inverter are the same. However, considering that the reference of each inverter is independently controlled, such as in a PV plant where current references vary along with maximum-powerpoint-tracking algorithm, the presented model lacks some essential physical significance and cannot comprehensively describe the characteristics of a PV plant. The authors in [13] applied a similar approach and presented models based on interactive current and common current to describe the interaction among multiple paralleled inverters.
Impedance-based small-signal analysis provides another method to evaluate the system stability [14], [15]. In this method, the impedances of source and load subsystems are measured and Nyquist stability criterion is applied to assess the system stability. While this method is characterized by simplicity, it does not provide any information on the system internal stability poles [16]. The system's internal stability was considered by the authors in [17], which provided an interpretation of the stability analysis presented in [12] using the source and load impedances. An alternative method for the impedance-based analysis was presented in [18], where the stability of a grid-connected inverter was characterized by the measured load-affected current controller loop.
In many modern distribution systems, the high share of renewable production combined with multiple parallel inverters introduce changes to the power system [19]. Especially in micro grids, the system operation status and configuration may fluctuate along with changes in loading state, power generation, or grid topology. As a consequence, it may become challenging to predict the detailed system characteristics required in most stability analysis methods. In highly fluctuating conditions, a stability analysis based on typical operation status may be insufficient. To tackle this issue, real-time analysis methods based on continuous online measurements are highly desirable. The real-time stability indications can be applied in advanced system monitoring, fault prevention, or adaptive control. The online identification and analysis can be integrated into the inverter controllers applying digital analyzer techniques [20], [21].
The present paper considers real-time stability analysis of paralleled grid-connected inverters by utilizing an OPAL-RT power hardware-in-the-loop (PHIL) setup that was developed recently at DNV GL, Arnhem, Netherlands. The OPAL-RT is a multi-purpose real-time simulator that is widely used in real-time analysis and control of various power-electronics applications including wind-turbine emulation [22], fuel-cell modeling [23], and analysis of smart-grid performance [24]. A PHIL setup consists of a hardware-in-the-loop simulator, such as OPAL-RT, and power hardware, such as physical power cables and devices [25]. Recently, PHIL setups have been applied in versatile experimenting of various power-electric systems [26]- [28].
The primary goal of this work is to provide detailed steps towards implementing a comprehensive real-time stability analysis for a system that has multiple paralleled inverters connected to a grid. The second goal of the work is to present the proof-of-concept for online stability analysis performed with parallel approaches simultaneously, and to compare the obtained stability indications. The method includes simultaneous use of loop gain analysis and impedance-based analysis applying an online multivariable measurement. The concept is applied to an experimental system at 50 kW power level.
The remainder of the paper is organized as follows. Section II reviews the theory of the frequency-response identification using orthogonal wideband techniques applied in the paper. Section III introduces the experimental high-power PHIL setup recently developed at DNV GL and presents the implementation of the online frequency-response measurement method. Section IV presents experimental results based on paralleled high-power inverters connected to a high-impedance grid and operated with different controller tunings. Finally, Section V draws conclusions.

A. Frequency-response measurement in the dq domain
Direct-quadrature (dq) transformation rotates the reference frame of three-phase systems in order to simplify the analysis of three-phase circuits. In the case of balanced three-phase circuits, the transformation reduces three AC quantities to two DC quantities, which makes it possible to utilize controllers with simpler structures and lower dynamic orders [29]. The dq-frame representation allows straightforward small-signal analysis as the non-linear characteristics can be linearized around the steady-state operation point [30].   2 shows a typical measurement setup where the system under test is to be identified in the dq domain [31]. Examples of such identification include the output impedance of a threephase inverter, inverter loop gain(s), and grid impedance [29]. In the setup, the system is perturbed by d-and q-component injections x d (n) and x q (n), which yield the corresponding output responses y d (n) and y q (n). The measured input and output signals, [x d (n),x q (n)] and [ŷ d (n),ŷ q (n)], are corrupted by input noise and output noise, respectively. The noise signals are assumed to resemble white noise and are uncorrelated with the other system variables. All of the signals are assumed to be zero mean sequences. The measured signals are buffered and segmented, and the signals are transformed to the frequencydomain by applying discrete Fourier transformation (DFT), given asX whereX(jω) is the signal transformed to the frequency domain. From the Fourier-transformed input and output signals, the frequency response (impedance or loop gain) can be obtained for each input/output pair by applying whereX(jω) is the Fourier-transformed input signal and Y (jω) is the Fourier-transformed output signal. In impedance measurements, for example, the input is a current signal (containing the nominal current and injected current perturbation) and the output is the resulting voltage signal. In noisy environments, the logarithmic averaging procedure [31] is often applied to compute the frequency response between desired variables as where P denotes the number of injected excitation periods. The method tends to cancel out the effect of uncorrelated noise from both input and output sides, so that the frequency response is obtained more accurately compared to conventional cross-correlation techniques [31]. In the dq-domain analysis the method is applied to each input-output couple, resulting in a 2×2 frequency-response matrix that represents a frequency response from each input (input d and q component) to each output (output d and q component). That is, the matrix includes the direct components as well as the cross-coupling components.

B. Maximum-length binary sequence
The perturbation design plays a very important role in obtaining the desired frequency response (or any other parametric or nonparametric system model) through the experiment described in Fig. 2. An optimal design leads to maximally informative experiments that is, experiments that extract the maximum amount of information and reduces operational costs associated with the identification procedure. For a linearsystem identification of sensitive systems, a binary signal such as maximum-length pseudo-random binary sequence (MLBS) most often offers the best possible choice in terms of maximizing signal power within time-domain-amplitude constraints [32]. The MLBS is a periodic broadband signal that has a largely controllable spectral energy distribution, and consequently, the measurements can be averaged over multiple periods to increase the signal-to-noise ratio [32]. The averaging procedure enables accurate online measurements even with very small injection amplitude, which may be a requirement in sensitive systems where large injection amplitude, such as in impulse injection, may disrupt the system operation. Another major advantage of the MLBS over the other types of signals such as sinusoids is that the sequence can be implemented with a low-cost system whose output can only generate a small number of signal levels. Due to the several favorable characteristics, the MLBS has become a popular perturbation signal in stability analysis of both AC and DC power distribution systems [33]- [37].

C. Orthogonal perturbations
Considering the multiple-input-multiple-output (MIMO) system in Fig. 2, the inputs and outputs are most often coupled. For such a system, the conventional approach is to apply the superposition theorem to the frequency-response measurements, where the excitation signal is separately injected to each system input one at a time, and all the output responses are measured for each input excitation [31]. Then, (3) is applied to each combination of input and output. However, the superposition approach requires multiple consecutive measurements, which increases the measurement time and may make the results prone to variations in the system under measurement. In this work, the MIMO system is measured applying orthogonal binary sequences, where multiple orthogonal injections are used to simultaneously excite all the system inputs. As the injections are orthogonal (that is, they have energy at different frequencies), the injections do not disturb each other and the MIMO system can be identified simultaneously within one measurement cycle [38]. Consequently, the technique has considerable advantages over the method using superposition theorem and sequential perturbations, as all the frequency responses are measured simultaneously under the same system operation conditions. The synthesis of orthogonal (periodic) binary sequences has been well documented [32]. One of the most popular techniques is to apply a modulation with rows of a Hadamard matrix [38]- [40]. In this method, a periodic binary sequence, such as the maximum-length binary sequence (MLBS), is used as a base signal. The second signal is obtained by adding, modulo 2, the sequence [0 1 0 1 ...] to the MLBS. The third signal is obtained by adding the sequence [0 0 1 1 0 0 1 1 ...] to the MLBS, and so on. Note, that the sequence length of the i th orthogonal sequence is doubled compared to the length of the (i − 1) th sequence. Fig. 3 shows an example of two orthogonal binary sequences generated by the method. The length of the MLBS is 63 bits, and generated at 1 kHz. The second sequence is known as an inverse-repeat binary sequence (IRS) because the modulation inverts every other digit of the MLBS. As Fig. 3 shows, the energy of every other harmonic has zero value in the IRS, which means that the MLBS and IRS have energy at different frequencies.

III. SETUP IMPLEMENTATION
The power hardware-in-the-loop (PHIL) setup shown in active front-end. The power amplifiers can be freely configured to act as sources or loads depending on what kind of power system architecture is studied.
A. PHIL Implementation Fig. 4 presents a schematic diagram of the proposed realtime method for system stability analysis, which is implemented in OPAL-RT. The diagram consists of three sections: preprosessing, actual hardware, and postprocessing. The preprocessing section is responsible for the injection design and synthesis, where the orthogonal wideband perturbation signals are generated. The MLBS is implemented through shift registers and exclusive-or feedback, where the unit-delay blocks have a delay that corresponds to the generation frequency [41]. The orthogonal inverse-repeated sequence (IRS) is similarly generated in shift registers by processing the MLBS through Hadamard modulation [32]. Then, the perturbations are simultaneously injected to the physical system to d and q components, and the input and output signals are continuously measured and buffered. The injection point is the currentor voltage-references of the device used in measurements; for example, the loop and grid impedance measurements are performed by superimposing the perturbation to the current reference of an inverter. In the postprocessing, the data sequences are Fourier-transformed and averaged by using (3) yielding the frequency-response data. The obtained frequencyresponse matrix is then (simultaneously) used to complete the Bode plot(s), Nyquist countour(s), system loop gains, and system poles and zeros. The refresh rate of the frequencyresponse matrix is P N/f g , where P is the number of injection periods, N is the injection(s) length, and f g is the injection(s) generation frequency. In this implementation, the postprocessing takes place in the OPAL-RT. In real implementations, however, the postprocessing can be integrated on the inverter controller structure, which also enables easy system scaling for higher number of inverters, as the computational effort for an independent inverter is unchanged. Online techniques   with similar computational demands have been successfully implemented on the converter controller hardware [20], [21], [42]. The designed perturbation signals are injected to the system through an adjustable gain (K). The optimal excitation amplitude depends on the system under measurement, for example noise level and nonlinear behavior may deteriorate the measurement accuracy. In challenging measurement conditions, the excitation can be tuned to have more energy to ensure sufficient signal-to-noise ratio. It is also possible to automatically or adaptively adjust the injection amplitude; for example, one can continuously measure the system total harmonic distortion (THD), and select the injection amplitude according to the THD limitations.

IV. EXPERIMENTS AND ANALYSIS
A. Experimental setup A paralleled inverter system connected to a power grid was implemented using the PHIL setup (Fig. 5). Fig. 6 shows a schematic diagram of the system that consists of two inverters connected to a point of common coupling (PCC). Each inverter consists of an amplifier unit, output filter, internal current controller and phased-locked loop (PLL). Fig. 7 shows a detailed dq-domain control system of Inverter 1, where the injection and measured variables are highlighted with red. The filters were physically implemented by using a 0.5 mH threephase inductor in Inverter 1 and 3.2 mH inductor in Inverter 2. The third amplifier group acts as a grid emulator, which provides the stiff grid voltages for the system and sinks the power from the inverters. An inductor of 1.0 mH was applied for the grid impedance emulation. The inverters operate at 50 kW total power. The detailed system parameters are shown in Table I in Appendix.

B. Experimental measurements
Two orthogonal binary perturbations were designed and applied in all experiments: the MLBS for measuring the desired d current/voltage component (together with the cross-coupling component from d to q) and the IRS for measuring the desired q current/voltage component (together with the cross-coupling component from q to d). The MLBS was synthesized using a 11-bit-length shift register. Thus, the MLBS was 2047-bit long, so it yielded a frequency resolution (i.e. the spacing between adjacent measured frequencies) of approximately 2.4 Hz. The length of the IRS was, by definition, twice the length of the MLBS. In each experiment, the MLBS and IRS were simultaneously injected into the system at 5 kHz. The injection amplitudes were adjusted such that output voltages and currents did not deviate from their nominal values by more than 5 %. Considering the inherent harmonics present in most grids, the injection can be designed so that the decrease in power quality is negligible.
Three sets of experiments were performed to analyze the system stability in different conditions. In the experiments, the inverters' PLL bandwidth was tuned differently: in the first experiment the bandwidth was set to 20 Hz, in the second to 60 Hz, and finally in the last experiment to 85 Hz. As previously reported [43], increasing the PLL bandwidth in weak grid generally decreases the system robustness, and altering the PLL tuning was chosen so that the proposed approach can be tested in versatile stability conditions. In all experiments, the full 2x2 frequency-response matrix including the direct components (d and q) and cross-coupling components (dq and qd) was measured.
The experiments were performed by first utilizing the inverters in measurements (injection points 1 and 2 in Fig.  6), and then performing the grid measurements (injection point 3). The injection point in the inverter measurements allows measuring both the current controller loop gains and grid impedance simultaneously with the same injection. As a consequence, this process eliminates the need for separate measurements for loop gains and grid impedance. In the loop gain measurements, the input vector is the current controller reference i in = [i in-d , i in-d ] T and the output vector is the current response i o = [i d , i q ] T (as shown in Fig. 7). Similarly, in the grid impedance measurements, the output current i o is the excitation and the output vector is the voltage response v pcc = [v d , v q ] T . The measurements are sampled at 10 kHz and the impedances and loop gains are extracted from the measured signals by applying equations (1)-(3) to each input/output pair. In order to avoid spectral leakage, the measured signals are segmented to consist of 4094 samples for the MLBSperturbed signals and 8188 samples for the IRS-perturbed signals (as the sampling frequency is twice as high as the generation frequency of the perturbations). At first, Inverter 2 was disconnected and the loop and impedance measurements were performed with Inverter 1. After this, Inverter 2 was reconnected and the injection was applied to Injection Point 2 to obtain the current control loop of the second inverter. It would have been possible to simultaneously measure the current loop gains from both inverters but this approach would have required the use of four orthogonal sequences. Finally, the injections were applied to Injection Point 3 for obtaining the aggregated output impedance of the paralleled inverters. In all measurements, (1) was applied over 100 injection periods for the MLBS and over 50 injection periods for the IRS (as the length of the IRS is doubled compared to the MLBS). The refreshment rate, which is equal to the duration of averaged measurement cycle, is 41 seconds during which the computation is performed. With modern inverter controllers, transfer function calculation and stability analysis can be easily performed within this time frame enabling real-time stability assessment. When very fast (millisecond range) measurement cycles are applied, the computational requirements become stricter and the hardware performance must be considered in the method tuning. the system with 20 Hz PLL bandwidths, the phase margin is approximately 57 degrees at the intersection frequency around 190 Hz. Increasing the bandwidth to 60 Hz reduces the phase margin to 13 degrees, and increase to 85 Hz decreases phase margin to 4 degrees indicating almost marginal stability. However, this single-input-single-output consideration ignores the inherent multivariable nature of the system, and reliable analysis requires the use of complete 2x2 matrices [6], which is presented in the following section. Similarly, the 2x2 matrices of the load-affected internal current control loops are shown in Fig. 9 for Inverter 1 and in Fig. 10 for Inverter 2.

C. Stability analysis
The stability analysis can be performed on the measured PCC impedances or internal loop gains, where similar methods can be applied. In this work, the stability analysis is performed by generalized Nyquist criterion (GNC) for the impedances, and pole-zero fitting for the internal loop gains. These methods can be applied simultaneously, as the measurements can be obtained with a single injection cycle as discussed previously. In the GNC, the distance from the eigenvalue contour to the critical point shows the system robustness. Fig. 11 presents the eigenvalue contours for the system for different experiments. A drastic change towards instability is seen when the PLL bandwidth is increased as the eigenvalue 1 contour shifts towards the critical point. With 85 Hz PLL bandwidth, the system is almost marginally stable.
In addition to the impedance-based GNC analysis, the system stability was analyzed through system poles and zeros. The poles and zeros were continuously calculated from the same impedance data, and additionally, from the internal   Fig. 11: Impedance-based Nyquist contours with different PLL tunings. current controller loops for each inverter. Consequently, the method enabled the use of two parallel approaches. For evaluating the values of zeros and poles, parametric models were continuously estimated based on the measured loop gains [44]. Fig. 12 shows the estimated poles and zeros based on the minor loop (impedance ratio) and loop gains from inverters 1 and 2. The critical poles shown in bottom right subfigure shift towards the imaginary axis and approximately the same stability margins can be obtained through all the methods.
The use of parallel approaches increases the flexibility and usability of the stability analysis method, as both approaches have their advantages. The impedance-based approach excels in system-level analysis and monitoring, but becomes more time consuming in very complex systems with many inverters. On the other hand, the loop-gain-based analysis is readily available even in such systems, and additionally provides detailed information on individual-inverter level. This local data can be utilized in, for example, fault prevention by identifying the low-stability devices or in independent adaptive control of the inverters.

D. Time-domain verification
The presented methods assess the stability in the frequency domain, as the frequency-domain results are significantly more straightforward to derive. However, ultimately the interest on system performance often relies in the time domain operation, and the frequency domain is applied as an auxiliary domain. In order to validate the obtained stability margins and to verify the proposed methodology, a comparison between transient performance of the actual system and the predicted stability margins is performed in the time domain. In the following analysis, the minor loop (impedance ratio) poles and zeros were used to construct the corresponding dynamics in the time domain by predicting the time-domain step response from the obtained frequency-domain margins. The validation was performed offline by first constructing transfer functions from the pole zero data (as shown in (4) in the footnote for the 85 Hz PLL system). Then, the time-domain prediction was obtained by applying MATLAB step-function on the pole-zero data, where the step response of a transfer function is calculated numerically. The calculated response is compared to the actual response of the experimental system which shown in Fig. 6. Fig. 13 presents the predicted system behavior and the actual response to a step change in the current reference of the Inverter 1. The q channel current reference is changed from -10 G(s) = 2.03e05 * (s + 3.57e03)(s + 1.99e03) (s + 1.31e06)(s + 3.89e01 + j1.15e03)(s + 3.89e01 − j1.15e03) = 2.03e05s 2 + 1.13e09s + 1.45e12 s 3 + 1.31e06s 2 + 1.03e08s + 1.74e12 (4) to 10 A, and the response is measured for all the test parameter sets. The figure shows that the actual system response is very accurately predicted with the proposed method, as both the oscillatory frequency and the damping ratio closely match the experimental waveform. A steady state error that persists after the oscillations have settled results from the minor deviations in the fitted locations of the poles and zeros, which accumulate to the DC-gain of the system. From the stability perspective, the significant characteristics that quantify the system stability are the overshoot, oscillatory frequency, and the damping ratio (that is, the rate of decay in the oscillations). With the proposed methods, these were accurately predicted even for the system that was close to marginal stability.

V. CONCLUSIONS
Paralleled inverters play an important role in the operation of most grid-connected system. A number of stability issues arise from the interactions among multiple inverter subsystems and the power grid. Recent studies have presented techniques, such as impedance-based analysis, to assess the stability and dynamic characteristics of multi-inverter systems. Frequencyresponse measurements are often required to extract information from these systems, as the internal dynamics of some subsystems may be unknown. Especially online methods for measurements are analysis are highly desirable, in order to provide real-time data on the system status, which can be utilized in advanced system monitoring, fault protection, and adaptive control. This paper has presented a practical realtime approach to perform a comprehensive stability analysis of a three-phase grid-connected system containing paralleled inverters by applying a power hardware-in-the-loop platform. The method applies orthogonal injections for multivariable system identification and utilizes parallel approaches for stability analysis by applying impedance-based analysis and loop-gain assessment simultaneously. Consequently, a detailed information can be obtained from both local and global point of view for holistic stability assessment. Experimental measurements were shown from a high-power energy distribution system recently developed at DNV GL, Arnhem, Netherlands. The presented methods can be used to modify various system characteristics, such as impedance behavior and control dynamics, in real time, thereby providing means for several stability and adaptive control design tools for grid-connected systems containing paralleled inverters. His research interests include grid-connected threephase power converters for renewable energy applications and microgrids, dynamic modeling, control design, impedance-based interactions in three-phase systems, and impedance design of grid-connected inverters.