HIL-Assessed Fast and Accurate Single-Phase Power Calculation Algorithm for Voltage Source Inverters Supplying to High Total Demand Distortion Nonlinear Loads

Abstract

The dynamic performance of the local control of single-phase voltage source inverters (VSIs) can be degraded when supplying to nonlinear loads (NLLs) in microgrids. When this control is based on the droop principles, a proper calculation of the active and reactive averaged powers (P–Q) is essential for a proficient dynamic response against abrupt NLL changes. In this work, a VSI supplying to an NLL was studied, focusing the attention on the P–Q calculation stage. This stage first generated the direct and in-quadrature signals from the measured load current through a second-order generalized integrator (SOGI). Then, the instantaneous power quantities were obtained by multiplying each filtered current by the output voltage, and filtered later by utilizing a SOGI to acquire the averaged P–Q parameters. The proposed algorithm was compared with previous proposals, while keeping the active power steady-state ripple constant, which resulted in a faster calculation of the averaged active power. In this case, the steady-state averaged reactive power presented less ripple than the best proposal to which it was compared. When reducing the velocity of the proposed algorithm for the active power, it also showed a reduction in its steady-state ripple. Simulations, hardware-in-the-loop, and experimental tests were carried out to verify the effectiveness of the proposal.

1. Introduction

The Smart Grid, as a concept, can be defined in terms of its outcomes as an electrical system that operates in an efficient manner, that provides a reliable energy supply and a power quality for the needs of a digital economy, that demonstrates a resilient performance against uncertainties or grid faults, and that integrates a large variety of distributed energy resources (DERs) in the conventional electrical grid, especially renewable energy sources (RESs), according to the Department of Energy of the United States [1] and the European Commission [2]. In addition, following the United Nations Sustainable Development Goals [3], the energy supply must be affordable, reliable, sustainable, and accessible for all users, especially when involving DERs based on RESs, as can be seen in References [4,5]. These DERs can be composed of a mixture of energy production units, energy storage systems (ESSs), and loads that operate jointly in clusters that are connected or not to the main electrical grid infrastructure [6]. This set of interconnected distributed energy resources and loads operating as a single entity that is connected or not to the grid is known as a microgrid (MG), according to References [7,8]. Other similar definitions and analyses of MGs’ signal remarkable aspects, such as their electrical boundaries, the energy flow between DERs when sharing loads, and the islanded and grid-connected operation [9,10,11,12]. A specific type of MGs is the hybrid renewable MG (HRMG), which can be based on RESs and conventional generators to provide supply reliability against the intermittent nature of some RESs, such as photovoltaic, wind, and tidal energy, as seen in References [13,14,15,16]. These HRMGs also contribute to the desired resiliency of the grid in front of strong perturbations. The energy sources present in the MGs are usually interfaced with an AC bus through parallelized voltage source inverters (VSIs). When operating in islanded or connected mode, the MGs’ control is generally achieved through a hierarchical structure based on three levels of control [17,18,19]. At the tertiary level of control, the primary purpose is to manage the power injection to the grid according to economic and efficiency optimization parameters or the regulatory framework. The secondary level of control coordinates all the MG elements through the AC bus voltage amplitude and frequency. Finally, a primary level of control is achieved for the self-regulation of the VSI inside of an MG in islanded or grid-connected mode. The variables employed here for the regulation are the amplitude and frequency of the AC bus voltage and the phase-shift between this voltage and the injected current. This control is usually performed without intercommunication, and the time of response can reach the order of milliseconds. In contrast, the secondary and tertiary levels need a reliable intercommunication layer, and response times are substantially higher. In Reference [20], the impact of the lack of intercommunication between the distributed generators at all levels of control was analyzed by comparing diverse control methods, i.e., consensus-based, agent-based, decomposition-based, finite-time convergence control, aperiodic sampling, sliding mode control, and droop-based control. From these techniques, the droop-based methods presented some remarkable advantages: They do not require intercommunication, they offer a low cost of implementation, and they are suitable for grid-connected and islanded modes of operation [6,18,19,20]. These qualities facilitate the use of droop-based methods in the local control for both hierarchically and non-hierarchically controlled DC, three-phase MGs [17,18,21], or single-phase MGs [22,23,24]. The VSIs droop-operated contribute to the MG resilience, since they can be employed to stabilize the AC bus voltage and frequency with proper control. However, in the case of a VSI supplying to nonlinear loads in islanded mode, slow responses and degradation of power quality have been reported [20,25]. Therefore, the objective of this work was to contribute to the enhancement of the dynamic performance of droop-based single-phase VSI when supplying to nonlinear loads by doing the following:

• Improving the velocity of the transient response of the averaged power calculation.
• Improving the steady-state accuracy of the averaged power calculation.

This paper is organized as follows: In Section 2, the droop-based control techniques suitable for linear or nonlinear loads are described. Then, a single-phase VSI system supplying to a characterized nonlinear load (NLL) is presented in Section 3. Section 4 presents a set of averaged active and reactive power (P–Q) calculation algorithms, and Section 5 describes the novel proposed method. The results of Matlab/Simulink simulations are presented in Section 6, to compare the dynamic performances of the different calculation algorithms. Section 7 gives the results of hardware-in-loop (HIL) tests to assess the simulation results. The results of the experimental tests for the evaluation of the proposal are provided in Section 8.

2. Droop-Based Local Control Techniques against Nonlinear Loads

The droop-based local control VSIs, which are parallelized for sharing loads, need to calculate P–Q to generate the sinusoidal references to follow, as can be seen in References [17,26,27,28,29,30]. Moreover, considering a mainly inductive line transmission is desirable for the P–Q calculations and has yielded beneficial results in low-voltage systems [18,30,31,32,33].

The calculation of the P–Q parameters is fundamental for these droop-operated VSIs. When functioning in islanded mode, single-phase MGs are weaker in terms of stability, as compared with three-phase VSI-based systems, especially when sharing NLLs. In these conditions, accurate and fast calculations of P–Q become crucial [34,35]. Different solutions for the calculation of single-phase P–Q can be found in the literature. It is usual to see that, for the P–Q calculations, a voltage signal delayed 90° to the measured one in the point of common coupling (PCC) is employed. After that, the instantaneous active and reactive powers are calculated and later conditioned to obtain the P–Q parameters. For the in-quadrature voltage, one method consists of applying a transport delay by digital means without filtering the signal in magnitude [36,37]. Other approaches are based on the p-q theory, which implies a dq-Synchronous Reference Frame technique [38,39]. In this case, it is also acquired the in-quadrature signal of the PCC current. Moreover, second-order generalized integrators (SOGIs) were used thanks to their low-pass filter capability (SOGI-LPF), as in References [40,41]. Coming up next, the calculations of the instantaneous power quantities signals are conditioned for the extraction of the averaged active power (P), and the averaged reactive power (Q). One method consists of the application of SOGI-LPF [42]. Thus, a better strategy for improving the transient response against abrupt load changes is performed by adding a final stage for removing the double frequency power components, either with a final low-pass filter (LPF) stage in Reference [41] or without final LPF [35]. In the same conditions, in Reference [43] a different approach based on the discrete Fourier Transformation was utilized for directly obtaining P and Q. All of these strategies have in common that were only designed for droop-operated single-phase VSI sharing linear loads, introducing a significant delay that constraints the velocity of the response in front of abrupt load changes. The common goal of all those techniques is to obtain accurate and fast P–Q calculations to enhance the robustness of the parallelization of VSIs when sharing loads. Then, in Reference [35], it was proposed a method that introduced a pre-filtering of the measured current before the calculation of the instantaneous powers that resulted in being faster in front of NLL. This pre-filtering was achieved by the band-pass filter (BPF) capability of a double SOGI approach (DSOGI), followed by the remotion of the double frequency power components like Reference [34], but without the LPF final stage. The fastest solution from those benchmarked was [35], but at the cost of increased complexity in the calculation scheme and a worsening of the reactive averaged power calculation. The algorithm was designed based on the high distortion of the NLL, increasing the signal conditioning blocks with respect to previous proposals. Thus, the NLL was presented as a highly polluted in harmonics current, but was not quantitatively characterized.

The algorithm proposed in this work is designed by firstly characterizing the NLL accordingly to the well-known standard IEEE std 519-2014 [44]. Consequently, the topology of the calculation structure is modified, and the conditioning signal blocks are reduced with respect to Reference [35]. The objective is to obtain a faster and more accurate calculation of the P–Q parameters.

3. Description of the System under Test: VSI Supplying to a Nonlinear Load

Figure 1 represents a basic scheme of a single-phase MG topology containing two VSI sharing an NLL. The block “P–Q Power Calculation” is the research object of this paper. The voltage v_o(t) and the current i_o(t) are measured at the PCC when the switches S₀ and S₁ are open. The averaged active and reactive power, P and Q, are respectively calculated in the “P–Q Power Calculation” block from these measurements. Then a voltage reference v_ref(t) is generated. Later, the inner control loops based on [45] employ this reference and finally generate a pulse-width modulation (PWM) for the switching of the H-Bridge.

The NLL is a load which voltage–current characteristic is not linear. Different types of NLLs are found in the literature, which may be classified as attending different criteria. They are either characterized only considering the V–I characteristic [46], or based on power quality parameters associated with the active, reactive, and distorted power [47] or on the measured current [48]. In References [49,50], the classification is focused on the harmonic components only present in the current through the NLL. Similarly, Reference [51] studies an NLL accordingly to its current harmonic components and its total harmonic distortion in the current. The NLL study in Reference [52] is based on the frequency spectrum and the total harmonic distortion in voltage. The particular NLL load of this work (Z_NL) consists of an unbalanced diode-bridge rectifier (DBR) that supplies power to a resistive–capacitive (R–C) load, shown in Figure 2. The R–C load parameters are listed in

Table 1. Parameters for Z_NL.

Parameter	Value
RonD1-D3/RonD2-D4	0.01 Ω/1 Ω
LL/C1 = C2	84 µH/470 µF
RC1 = RC2/RL1 = RL2	37 kΩ/1560 Ω

and are characterized accordingly to power quality parameters listed in

Table 2. Characterization of Z_NL considering IEEE std 519-2014 [44].

Parameter	Value
Individual THDv; Total THD	<5%; <8%
Individual THDi	>4% for 3 < odd harmonic < 11 >1% for even harmonics
Total TDD in current	>5%

, based on IEEE std 519-2014 [44]:

After reaching a steady-state, the switch S2 closes at t = t_on, triggering an abrupt change in the value of P and Q.

This load is the same as the one employed in Reference [35], and the voltage measured in the PCC follows Expression (1). The current in the PCC can be described by (2), and draws a highly distorted waveform:

$$v_o\left(t\right)=V_o·sin(w_o·t)$$

(1)

$$i_o\left(t\right)=I_{DC}+I_0·sin\left({\omega }_{o}t-{\phi }_o\right)+{\displaystyle \sum }_{h=2}^N{I}_h·sin\left(h·{\omega }_o·t-{\phi }_h\right)$$

(2)

where V_o and I_o are the voltage and current amplitudes, respectively; ${\omega }_o$ is the fundamental frequency of the system (100π rad/s); ${\phi }_o$ is the phase-shift between the fundamental components of the voltage and the current; h is the harmonic index; I_h is the amplitude of the harmonic components of the intensity; and the phase-shift ${\phi }_h$ corresponds to each current harmonic component. The term $I_{DC}$ corresponds to a DC offset present in the load current.

Figure 3a detailed the local control structure for a single-phase VSI supplying the NLL of Figure 2. This model is based on References [45,53], and this work only focuses on the block named “POWER CALCULATION BLOCK,” that generates the voltage reference for the voltage and current inner loops. Those consist of a proportional-integral (PI) block for the voltage and a proportional-resonant (PR) loop for the current. A final stage for the PWM for the switching of an H-bridge is shown.

Figure 3b displays the harmonic distribution of the measured current. Regarding the recommended limitations in [44], the individual harmonic distortion in the current (THDi) and the total distortion demand (TDD) are out of limits. Moreover, a DC component in (2) and pictured in Figure 3b is not allowed by Reference [44]. Finally, the individual harmonic distortion in voltage, THDv, and its total harmonic distortion, THD, are within limits. Table 2 summarizes this characterization of Z_NL:

4. P–Q Calculation Algorithms

The droop equations for a mainly inductive system are as follows:

$${\omega }^{\ast }={\omega }_n-m·P$$

(3)

$$V^{\ast }=V_n-n·Q$$

(4)

where ${\omega }^{\ast }$ is the calculated angular frequency of the system, $V^{\ast }$ the calculated voltage amplitude, $V_n$ is the rated value for the voltage, ${\omega }_n$ is the frequency rated value, and m and n are the droop coefficients. With these parameters, it is generated the sinusoidal voltage reference $v_{ref}\left(t\right)$ necessary for the inverter inner control loops of Figure 3a:

$$v_{ref}\left(\mathrm{t}\right)=V^{\ast }·sin\left({\omega }^{\ast }·t\right)$$

(5)

The calculation of P–Q is done as follows, in the time domain:

$$p\left(\mathrm{t}\right)=v_o\left(\mathrm{t}\right)·i_o\left(\mathrm{t}\right)=\mathrm{P}+\tilde{p}\left(\mathrm{t}\right)+v_o\left(t\right)·{\displaystyle \sum }_{h\ne 1}^N{I}_h·sin\left(h·{\omega }_o·t-{\phi }_h\right)$$

(6)

$$\mathrm{P}=\frac{1}{2}{V}_o·I_o·cos\left({\phi }_o\right)$$

(7)

$$\tilde{p}\left(t\right)=\frac{1}{2}{V}_o·I_o·cos\left(2{\omega }_0·t-{\phi }_o\right)$$

(8)

where (6) is the calculated active instantaneous power [38], P is the averaged active power in (7), and a double frequency pulsating component, $\tilde{p}\left(t\right)$, is represented by (8). For the instantaneous reactive power calculation, an in-quadrature signal for the voltage is employed, represented by (9). Then an instantaneous reactive power quantity is calculated:

$$v_{o\mathrm{q}}\left(t\right)=V_o·sin\left({\omega }_o·t-\frac{\pi }{2}\right)$$

(9)

$$q\left(\mathrm{t}\right)=v_{oq}\left(\mathrm{t}\right)·i_o\left(\mathrm{t}\right)=\mathrm{Q}+\tilde{q}\left(\mathrm{t}\right)+v_{oq}\left(t\right)·{\displaystyle \sum }_{h\ne 1}^N{I}_h·sin\left(h·{\omega }_o·t+{\phi }_h\right)$$

(10)

$$\mathrm{Q}=\frac{1}{2}{V}_o·I_o·sin\left({\phi }_o\right)$$

(11)

$$\tilde{q}\left(t\right)=\frac{1}{2}{V}_o·I_o·sin\left(2{\omega }_0·t-{\phi }_o\right)$$

(12)

where q(t) is the instantaneous reactive power, (10), Q is the averaged reactive power in (11), and $\tilde{q}\left(t\right)$ is a double frequency component. Only in Reference [35], the highly distorted load current was pre-filtered prior to the instantaneous power calculations in (6) and (10), using a DSOGI approach and its BPF capability. Thus, the LPF of a SOGI was employed as a quadrature signal generator (QSG) delaying in $\frac{\pi }{2}$ rad the voltage signal. The BPF and LPF transfer functions for a SOGI are described in (13) and (14), respectively:

$$H_d\left(s\right)=\frac{2\xi \omega s}{s^2+2\xi \omega s+{\omega }^2 }$$

(13)

$$H_q\left(s\right)=\frac{2\xi {\omega }^2}{s^2+2\xi \omega s+{\omega }^2}$$

(14)

where $\xi$ is the damping factor, and $\omega$is the center frequency of the system. Figure 4 shows the structure of a SOGI, with its BPF and LPF magnitude and phase Bode plots.

In Figure 4b, the SOGI-BPF is more selective while reducing the damping factor. The same occurs with its LPF capability. The attenuation presents a rate of −30 dB/decade for the BPF and −60 dB/decade for the LPF. However, the LPF transfer function presents a drawback in attenuation for frequencies below the fundamental.

The following Figure 5 shows three averaged active and reactive power calculation:

The obtained averaged active powers are named as PavC, PavT, and PavM, for Figure 5a–c, respectively. In the same manner, the averaged reactive powers are named as QavC, QacT, and QavM.

In Figure 5a, $v_{oqC}\left(t\right)$ is the $\frac{\pi }{2}$ delayed voltage signal, as described in (9), achieved by using a time-delay block [47]. In Figure 5b,c, a SOGI-LPF was employed as QSG for the obtention of $v_{oqT}\left(t\right)$ and $v_{oqM}\left(t\right)$, as well as BPF for the direct component of the voltage.

The next common stage of the three calculation algorithms consists of obtaining the instantaneous active and reactive powers, as described in (6) and (10). Note that the measured current is directly applied, except in Figure 5c, where it utilizes a BPF for the current [35].

In Figure 5a, Pavc and Qavc are obtained after an LPF stage. On Figure 5b, the oscillatory double frequency components are extracted by a SOGI-BPF tuned at $2{\omega }_0$ and then removed from the calculated instantaneous powers. Lastly, an LPF similar to Figure 5a is employed to reduce the steady-state ripple. The latest algorithm in Figure 5c also removes the double frequency components but eludes the final LPF stage.

Besides, Figure 5 shows the increasing complexity of the algorithms, especially Figure 5c, when supplying to the NLL in Reference [35]. Therefore, the next section proposes a novel algorithm that calculates in a faster manner the averaged active power and, in a more accurate manner, the averaged reactive power, while reducing the complexity of the calculation scheme.

5. Proposed P–Q Calculation Algorithm

The proposed algorithm for the calculation of P–Q is presented in Figure 6. This new algorithm is intended to ease the dynamic performance of the system by reducing the settling-time during abrupt load changes.

Figure 6a shows that the measured current, $i_o\left(t\right)$, is conditioned through SOGI-0, obtaining a BPF filtered component, $i_{od}\left(t\right)$, and an LPF filtered component, $i_{oq}\left(t\right)$. Then, each current component is directly multiplied by the voltage, computing the instantaneous active power, $p_i^{\prime }\left(t\right)$, and the instantaneous reactive power, $q_i^{\prime }\left(t\right)$. Those quantities are expected to contain the averaged active and reactive powers Pav and Qav, respectively, plus a certain amount of undesired harmonic components, similarly to Reference [35] and following Equations (6) and (10). Therefore, a final LPF stage is applied to each instantaneous power to acquire active and reactive quantities with the lesser possible steady-state harmonics. Those last active and reactive power quantities are, respectively, $p_F\left(t\right)$ and $q_F\left(t\right)$, and contain the desired Pav and Qav. Comparing this scheme to those pictured in Figure 5, note that the voltage signal is not conditioned due to the specific NLL characterized according to Table 2. The proposed scheme also presents fewer signal conditioners and control parameters than in Reference [35], in Figure 5c, showing a simplified calculation structure.

Figure 6b is the frequency-domain analytical representation of the calculated quantities and the transfer functions for each SOGI employed in Figure 6a. The signals reported in this scheme are the Laplace Transform of those indicated in Figure 6a: $I_0\left(s\right)=\mathfrak{L}\left\{i_0\left(t\right)\right\}$, $I_{0d}\left(s\right)=\mathfrak{L}\left\{i_{0d}\left(t\right)\right\}$, $I_{0q}\left(s\right)=\mathfrak{L}\left\{i_{0q}\left(t\right)\right\}$, $V_0\left(s\right)=\mathfrak{L}\left\{v_0\left(t\right)\right\}, P_i^{\prime }\left(s\right)=\mathfrak{L}\left\{p_i^{\prime }\left(t\right)\right\}$, $Q_i^{\prime }\left(s\right)=\mathfrak{L}\left\{q_i^{\prime }\left(t\right)\right\}$, $P_F\left(s\right)=\mathfrak{L}\left\{p_F\left(t\right)\right\}$, and $Q_F\left(s\right)=\mathfrak{L}\left\{p_F\left(t\right)\right\}$. Analytically, the LPF and BPF transfer functions of SOGI0 are represented, respectively, by $H_{d0}^{\prime }\left(s\right)$ and $H_{q0}^{\prime }\left(s\right)$:

$$H_{q0}^{\prime }\left(s\right)=\frac{2{\xi }_i{\omega }_0^2}{s^2+2{\xi }_i{\omega }_{0}s+{\omega }_0^2 }$$

(15)

$$H_{d0}^{\prime }\left(s\right)=\frac{2{\xi }_i{\omega }_{0}s}{s^2+2{\xi }_i{\omega }_{0}s+{\omega }_0^2 }$$

(16)

Those are tuned at ${\omega }_0$, and their selectivities are controlled through the damping factor, ${\xi }_i$. Hence, the direct and the in-quadrature filtered currents, $i_{od}\left(s\right)$ and $i_{oq}\left(s\right)$, respectively, are as follows:

$$i_{od}\left(s\right)=\left(\frac{2{\xi }_i{\omega }_o·s}{s^2+2{\xi }_i{\omega }_o·s+{\omega }_0^2}\right)i_o\left(s\right)$$

(17)

$$i_{oq}\left(s\right)=\left(\frac{2{\xi }_i{\omega }_o^2}{s^2+2{\xi }_i{\omega }_o·s+{\omega }_0^2}\right)i_o\left(s\right)$$

(18)

and the time-domain expressions are:

$$i_{od}\left(t\right)=I_0·sin\left({\omega }_{o}t-{\phi }_o\right)+{\displaystyle \sum }_{h\ne 1}^N I_{dh}·sin\left(h·{\omega }_o·t-{\phi }_{dh}^{}\right)$$

(19)

$$i_{oq}\left(t\right)=I_0·sin\left({\omega }_{o}t-{\phi }_o-\frac{\pi }{2}\right)+{\displaystyle \sum }_{h\ne 1}^N{I}_{qh}·sin\left(h·{\omega }_o·t-{\phi }_{qh}^{}\right)$$

(20)

where $h\ne 1$ is a harmonic index. For (19), $I_{dh}$ is the harmonic amplitude, and ${\phi }_{dh}$ is its phase-shift. In (20), $I_{qh}$ is the harmonic amplitude and ${\phi }_{qh}$, is its phase-shift. Note that $I_{dh}\ne I_{qh}$ and ${\phi }_{dh}\ne {\phi }_{qh}$.

Later, following the scheme in Figure 6a, the instantaneous active and reactive powers, $p_i^{\prime }\left(t\right)$ and $q_i^{\prime }\left(t\right)$, are as follows:

$$p{\prime }_i\left(t\right)=P_{av}+v_o\left(t\right)·{\displaystyle \sum }_{h\ne 1}^N{I}_{dh}·sin\left(h·{\omega }_o·t-{\phi }_{dh}\right)$$

(21)

$$q{\prime }_i\left(t\right)=-Q_{av}+v_o\left(t\right)·{\displaystyle \sum }_{h\ne 1}^N{I}_{qh}·sin\left(h·{\omega }_o·t+{\phi }_{qh}\right)$$

(22)

$$P_i^{\prime }\left(s\right)=\mathfrak{L}\left\{p_i^{\prime }\left(t\right)\right\}$$

(23)

$$Q_i^{\prime }\left(s\right)=\mathfrak{L}\left\{q_i^{\prime }\left(t\right)\right\}$$

(24)

Pav and Qav are the averaged active and reactive power outputs, respectively. $P_i^{\prime }\left(s\right)$ and $Q_i^{\prime }$ are the domain frequency expression of (21) and (22), related in Figure 6b.

Therefore, for SOGI1 and SOGI2, they utilize their LPF capability as follows:

$$H_P\left(s\right)=\left(\frac{2{\xi }_p{h}_1^2{\omega }_o^2}{s^2+2{\xi }_p{h}_1{\omega }_{o}s+h_1^2{\omega }_o^2}\right)$$

(25)

$$H_Q\left(s\right)=\left(\frac{2{\xi }_p{h}_2^2{\omega }_o^2}{s^2+2{\xi }_p{h}_2{\omega }_{o}s+h_2^2{\omega }_o^2}\right)$$

(26)

where $H_P\left(s\right)$ and $H_Q\left(s\right)$ are the transfer functions for the LPF capability of SOGI1 and SOGI2, respectively. Note that both transfer functions are essentially the same, only differentiated by the h₁ and h₂ coefficients. Those coefficients $\left(h_1,h_2\right)\complement \left[0.05,0.5\right]$ are employed for the attenuation of subharmonics reducing the SOGI1 and SOGI2 LPF cutoff frequencies. Next, Figure 7 shows their magnitude and phase plots. For the sake of simplicity, $h_1=h_2=h_i$:

Then, $P_F\left(s\right)$ and $Q_F\left(s\right)$ are the result of the filtering of (23) and (24) by their respective LPF-SOGI:

$$P_F\left(s\right)=P_i^{\prime }\left(s\right)·H_P^{\prime }\left(s\right)$$

(27)

$$Q_F\left(s\right)=Q_i^{\prime }\left(s\right)·H_Q^{\prime }\left(s\right)$$

(28)

Finally, back to the time domain, the following can be found:

$$p_F\left(t\right)={\mathfrak{L}}^{-1}\left\{P_F\left(s\right)\right\}=P_{av}+{\displaystyle \sum }_{k=1}^X{p}_k\left(t\right)$$

(29)

$$q_F\left(t\right)={\mathfrak{L}}^{-1}\left\{Q_F\left(s\right)\right\}=-Q_{av}+{\displaystyle \sum }_{k=1}^X{q}_k\left(t\right)$$

(30)

The expressions $p_F\left(t\right)$and $q_F\left(t\right)$ contain Pav and Qav. The parameter k is a harmonic index similar to h in (21), being $p_k\left(t\right)$ and $q_k\left(t\right)$ the undesired oscillatory components. For the attenuation of those, SOGI1 is tuned at a frequency $\omega =h_1{\omega }_0$ and SOGI2 at $\omega =h_2{\omega }_0$. Low values for h₁ and h₂ will lead to a substantial reduction of these components. The next section includes simulations to study the values of ${\xi }_P$, ${\xi }_i$, h₁, and h₂ for a more accurate and faster calculation of the P–Q parameters.

6. Simulation Results

The proposed algorithm is simulated to compare its dynamic performance against the structures shown in Figure 5, with an abrupt load change after closing S₁. A similar steady-state ripple for the active power calculation is set as a reference for the analysis. The parameters of the simulations are listed in

Table 3. Simulation parameters.

Parameter	Value
$V_o$	311 V
$w_o$	100 π rad/s
TDD for io(t)	116.41%
${\xi }_i/{\xi }_p$	0.1 to 0.7075/0.7075
$h_1$	0.10 to 0.30 (in steps of 0.05)
$h_2$	0.10
${\omega }_C$; ${\omega }_{CT}$	0.74 π rad/s; 2.20 π rad/s
${\xi }_{v0}/{\xi }_1={\xi }_2$	0.7/1
${\xi }_{MV}/{\xi }_{MI}/{\xi }_A={\xi }_B$	0.7/0.14/1
C1 = C2/ RC1 = RC2	470 µF/37 kΩ
RL1 = RL2 = RL3 = RL4	1.8 mH; 0.01 Ω
RC BRANCH	25 µF; 1 Ω
SWITCHING FREQUENCY, fs	10 kHz

.

In Reference [35], it was demonstrated that the more suitable algorithm in the presence of an NLL was the calculation of PavM. The following figure shows a family of Pav plots after varying its control parameters and comparing it against PavM.

Figure 8a shows a family of P_av plots varying, $0.1\le {\xi }_i\le 0.7075$, while keeping constant $h_1=0.25$ and ${\xi }_p=0.7075$. It can be seen that, when increasing the damping factor, the transient response is faster. However, there is an undesired overshoot when ${\xi }_i$> 0.2. Therefore, the fastest configuration avoiding the overshoot is the one with ${\xi }_i$= 0.2.

Figure 8b shows the comparison between P_avM and a family of P_av calculations, varying $h_1$from 0.1 to 0.30 and ${\xi }_i$= 0.2. It can be observed how the rapidity of the proposed calculation is reduced when h₁ is increased. For $h_1=\left[0.1,0.15\right]$ the Pav transient results slower concerning PavM, and the steady-state ripple results smoother. When $h_1\ge 0.2$, the velocity of the transient increases, as well the ripple in steady-state. Only for $h_1=0.25$, it results in being faster, without introducing a remarkable overshoot. At higher values of $h_1$, it presents an undesired ripple in steady-state and an overshoot that is better to avoid.

Figure 9a shows the calculation of PavC, PavT, and PavM, compared with the proposed P_av with ${\xi }_{i }$ = 0.2. and $h_1=0.25$. Figure 9b is the detail of the steady-state ripple, where it can be seen that the fastest algorithm, Pav, shows a similar ripple to the other algorithms. In these conditions, the calculation of Pav results faster than PavM. Therefore, these parameters were chosen for the study of the employed for this last simulation.

In Figure 10, the spectrum and the steady-state THDs of Pav and PavM are compared. From the observation of Figure 9 and Figure 10, it is deduced that, for $0.15\le h1<0.25,$ there is a family of Pav calculations that result to be faster and more accurate than PavM.

The study of the reactive power quantities is done with h₂ = 0.1, ${\xi }_p$ = 0.7075 and with ${\xi }_i$ = 0.2. Figure 11 shows the calculation of Qav, QavM, QavT, and QavC when the abrupt load change occurs.

Figure 11a shows that QavM is the worst option in terms of steady-state ripple when supplying the NLL. Moreover, Qav presents a similar time response than QavM and the lowest steady-state ripple (see Figure 11b). Thus, similarly to the active power calculation, if the transient time needs to be reduced, that will succeed at the cost of a higher ripple. Nevertheless, although the reactive power calculation is necessary for the droop control, a variation of less than 10 VAr barely influences the droop reference generation due to the nature of the NLL. Considering the calculated mean value of Pav and Qav, it can be extracted a power factor (PF) equals to 0.9976. If the value of Qav increases up to 35 VAr, then PF = 0.9928. Although the load has not been characterized according to the PF, it indicates that the reactive power variations in mean value barely influence its value.

Then, it is deduced that the proposed calculation method for Pav can be faster than PavM for a range of values of ${\xi }_i$ and h₁, keeping ${\xi }_p$ = 0.7075. Moreover, when increasing the SOGI1-LPF capability, it maintains a better or similar settling-time while reducing the steady-state ripple. Thus, when Pav reaches a similar settling-time of that of PavM by reducing the h₁ parameter, the THD falls from 1.32% to 0.59% (Figure 10). The commented results are shown in

Table 4. THD, settling-time, and time-delay for the simulation of Pav.

	Calculated THD	Settling-Time (ms)/% Reduction with Respect to PAVM	Time-Delay (ms)/% Increasing with Respect to PAVM
$P_{avM}$	1.13%	120	38
$P_{av}\to$ ${\xi }_i=0.2$ $h_1=0.25$	1.32%	75/32%	40/5%
$P_{av}\to$ ${\xi }_i=0.2$ $h_1=0.15$	0.59%	90/18%	50/20%

:

As it can be seen in Table 4, Pav settling-time is 37.5% shorter than PavM when ${\xi }_i=0.2 and h_1=0.25$ with a similar ripple (Figure 11c). The time-delay for both calculations is almost the same, 38 ms for P_avM and 40 ms for P_av. Then, when $h_1=0.15$, the THD falls drastically down to 0.59% while keeping an 18% shorter settling-time. However, in this last case, the time-delay is higher in 20% for the Pav calculation concerning PavM. For this final reason, the chosen set of parameters for comparing Pav against PavM is are ${\xi }_i=0.2, {\xi }_p=0.7075, h_1=0.25$.

On the other hand, the simulation results for the reactive power calculation algorithms are compared in

Table 5. Settling-time and THD for the simulation of reactive power algorithms.

Reactive Power Calculation Algorithm	Calculated THD	Settling-Time (ms)	Reduction of THD with Respect to QavM
$Q_{avM}$	7.85%	140	Not applicable
$Q_{av}:{\xi }_i=0.2 h_{\mathrm{i}}=0.25;h_2=0.1$	2.46%	150	68.66%
$Q_{avT}$	3.49%	250	55.64%
$Q_{avC}$	1.87%	780	76.18%

, in terms of settling-time and comparing its steady-state ripple through a THD, with respect to the DC component analysis.

Table 5 shows how the conventional droop method, QavC, is the best option for reducing steady-state ripple in reactive power. However, its settling-time is the worst with a value of 780 ms. On the other hand, the best time is achieved by the QavM algorithm but at the cost of a higher THD = 7.85%. The conclusion is that the calculation of Q is more accurate through the proposed algorithm, with similar transient velocity. Therefore, the chosen set of parameters for the proposed algorithm is as follows:

$${\xi }_i=0.2 , {\xi }_p=0.7075, h_1=0.25 \mathrm{and} h_2=0.1$$

The next section pretends to assess the obtained results from the carried out simulations.

7. Hardware in the Loop Assessment

For the assessment of the proposed algorithm, HIL tests are carried out. Those tests compare the calculation of PavM and QavM against Pav and Qav with the chosen parameters from the simulations.

For this purpose, a real-time interface platform based on dSPACE 1006© (dSPACE Inc.50131 Pontiac Trail Wixom, MI, USA 48393-2020) digital platform is operated. The control structure presented in Figure 2 is firstly discretized in Matlab/Simulink/SimPowerSystems© (The MathWorks Inc., Natick, MA, USA) and then compiled in C code for its download in the dSPACE. Moreover, this RTI platform supports the model libraries of physical/electrical plants from Matlab/Simulink/SimPowerSystems© (The MathWorks Inc., Natick, MA, USA). Those libraries correspond to the modeled H-bridge, the LCL filter, and the NLL under test. The electronic central unit (ECU) of the dSPACE compiles the control algorithms on its multiprocessor core. The control-desk software permits the configuration and control of the tests, acting as a human–machine interface. The switching frequency of 10 kHz for the VSI is emulated, setting the sample time at Ts = 100 µs. Note that the discretization of the integrators employed in the SOGI has been achieved through a third-order method:

$$\frac{T_s}{12}·\frac{5z^{-3}-16z^{-2}+23z^{-1}}{1-z^{-1}}$$

(31)

A first HIL Test-1 is then carried out to compare the proposed algorithm against PavM, PavT, and PavC. The parameters for Pav are listed in

Table 6. P_av and Q_av parameters for the hardware-in-loop (HIL) tests.

Test	${\xi }_i$	${\xi }_p$	$h_1$; $h_2$
HIL TEST-1	0.2	0.7075	0.25; 0.10
HIL TEST-2	0.2	0.7075	0.15; 0.10

.

Figure 12 shows the load current plots; the active and reactive power after an abrupt load change is done manually.

Figure 12a corresponds to the current in the asymmetric NLL, with peaks of +2.85A and −1.8A after the load step.

Then, in Figure 12b, Pav is compared against PavM, PavT, and PavC. The shortest transient response corresponds to Pav, while the steady-state ripple is kept constant.

Figure 12c shows the transient response for the calculated reactive powers. There, it can be appreciated that the response rapidity is similar between Qav and QavM. Moreover, the lowest steady-state ripple corresponds to Qav. The relevant comparison here is between the proposed algorithm and that based on Reference [35]. The results are exposed in

Table 7. Settling-time and time-delay for active and reactive power.

HIL Calculation Algorithm	Settling-Time (ms)	Time-Delay (ms)
$P_{av}$	90	38
$P_{avM}$	140	40
$Q_{av}$	130	50
$Q_{avM}$	140	45

.

The settling-time for Pav is shorter in 35.7% with respect to PavM, similar to the simulation results with the same control parameters. The time-delay was found to be similar for Pav and PavM. However, in the reactive power, this time-delay is larger for Qav than for QavM, even when the Qav settling-time is a 7.1% minor than QavM.

A second HIL Test-2 is carried out, reducing the h₁ to 0.15 for Pav, to assess the simulation tests achieved in this sense. Figure 13 shows the active power responses during an abrupt load change.

Figure 13 shows that a low h1 coefficient in Pav allows a smaller steady-state ripple while maintaining a similar settling-time both for active and reactive powers, Figure 13 and Figure 12c, respectively. However, the time-delay for Pav is 50 ms, superior to the reported 40 ms in Table 7 for PavM.

8. Experimental Results

An experimental test was carried out, employing the load described in Table 2. The experimental setup is prepared to evaluate the model simulated and HIL-assessed for the active power calculation in the presence of a measured TDD = 124.9% in current at the PCC. The experimental setup is shown in Figure 14. It is composed of a VSI Danfoss© FC302, 2.2 kW rated, interfaced to a real-time dSPACE 1006 platform, for the switching signals for the H-bridge and the measured parameters. The current of the NLL was monitored by using a Fluke 435-II Power Quality and Energy Analyzer. The power calculation algorithms tested correspond to PavM and Pav. The results obtained are displayed in Figure 15 and

Table 8. Simulation, HIL, and experimental results for active power.

Algorithm	Settling-Time Matlab (ms)/(%Reduction)	Settling-Time HIL (ms)/(%Reduction)	Settling-Time Experimental (ms)/(%Reduction)
Pav	75/37.5%	90/(35.7%)	125/(30%)
PavM	120	140	180

.

First of all, the load draws an asymmetrical current, as shown in Figure 15a, with +2.19A/−1.46A peak values. The difference between the measured and the simulated and HIL tested current is due to power losses in the whole system.

Figure 15b shows the harmonic spectrum measured by employing a Fluke 435-II Power Quality and Energy Analyzer, which yields a TDD = 124.9%, compatible with the simulated in Figure 3b.

Figure 15c presents the PavM calculation, with the detail of the steady-state ripple. In the same manner, Figure 15d shows the proposed Pav calculation, with the detail of its steady-state ripple. Note that the transient for Pav results faster than for PavM, as expected. The measured settling-times are exposed in Table 8, along with the simulation and HIL results for active power.

The Simulation, the HIL test and experimental results for the settling-time in active power calculation are resumed in Table 8. The results compared here correspond to the Pav and PavM keeping the same steady-state ripple. The proposed Pav results to be 37.5% faster than PavM in the Simulation, 35.7% in the HIL tests, and 30% in the experimental test. The settling-time measured in the experimental test is a 40% higher than the obtained in the simulation, for Pav. In contrast, the measured settling-time for the HIL test is 16.67% higher than the simulated for Pav. Those differences are attributed to the 3rd order integrator employed both in HIL and the experimental setup. For the experimental setup, it may be considered latencies due to the internal communications and data acquisition boards, as well as nominal values biases of components and energy losses.

9. Conclusions

The proposed method enhances the dynamical performance in terms of rapidity and accuracy of the droop-based local control, which degrades in the presence of NLLs like the employed in this work, which was characterized considering IEEE std 519–2014. Only PavM demonstrated its suitability in the presence of NLL in Reference [35], focusing the calculation effort on the obtention of the fundamental component of the current and avoiding a final LPF stage. However, this previous work did not differentiate types of NLL. Oppositely, the proposed method characterizes the NLL, and then the algorithm architecture is decided. Thus, the implemented algorithm results less complex than those compared with, when supplying an NLL.

The proposed algorithm was compared with the previously studied algorithms and assessed through Matlab/Simulink simulation and a HIL test. Finally, an experimental test for the active power Pav and PavM was carried out to evaluate the proposed model. The main conclusions after analyzing the simulation, HIL, and experimental results are summarized as follows, in terms of transient response velocity and steady-state accuracy:

Velocity:

• Reduction of settling-time in 30% for the calculation of Pav with respect to Pav_M [35] while keeping a similar steady-state ripple (Table 8).

Accuracy:

• Active Power: Reduction in 47.78% of the steady-state calculated THD with respect to DC in the simulations for Pav, when the settling-time is similar (Table 4).
• Reactive Power: Reduction in 68.66% of the steady-state calculated THD with respect to DC in the simulations for Qav, when the settling-time is similar (Table 5).

As expected, the settling-time for the Pav calculation during an abrupt load change resulted in being smaller than the other compared methods. Moreover, from Table 8, it can be seen how the relative reduction of settling-time was preserved in all the scenarios, i.e., between 30% and 37.5%. That leads to a faster operation in the droop controlled VSI in the presence of high TDD NLLs, which points to an increase in single-phase MG stability when sharing NLLs.

Concerning the enhancement of the accuracy, it is noteworthy that the calculation of Qav results to be more accurate than the other methods, see Table 5. Regarding the active power calculation, the steady-state ripple can be smoother by reducing the settling-time for Pav, see Figure 8b and Figure 13.

Future investigations are intended to study the same issues when other NLL types are present, considering a well-known standard as IEEE std 519–2014. Further studies are also planned for the parallelization of single-phase VSI against different types of NLL to study its dynamic performance and control stability.