A Power Calculation Algorithm for Single-Phase Droop-Operated-Inverters Considering Linear and Nonlinear Loads HIL-Assessed

Jorge El Mariachet 1,* , Jose Matas 1 , Helena Martín 1 , Mingshen Li 2, Yajuan Guan 2 and Josep M. Guerrero 2 1 Electric Engineering Department, Polytechnic University of Catalonia (EEBE-UPC), 08019 Barcelona, Spain; jose.matas@upc.edu (J.M.); m.helena.martin@upc.edu (H.M.) 2 Energy Teknik Department, Aalborg University (ET-AAU), 9220 Aalborg, Denmark; msh@et.aau.dk (M.L.); ygu@et.aau.dk (Y.G.); joz@et.aau.dk (J.M.G.) * Correspondence: jorge.el.mariachet@upc.edu


Introduction
The parallelization of single-phase inverters without communications between them has usually been performed using the droop method, which drives the sinusoidal references of the inverters for sharing the common loads [1][2][3]. The method introduces proportional droops in the inverter frequency ω* and voltage amplitude V* references, respectively, according to the P and Q load consumed powers. These powers are usually obtained as the product of the measured output current i o (t) with the measured output voltage v o (t) and its quadrature version, v o⊥ (t), respectively. The droop method works well for the sharing of linear loads but not for nonlinear ones, due to the nonlinear currents drawn by these loads, which highly distort P and Q, and the provided droop references. Therefore, LPFs in P and in Q with a very low cut-off frequency are used to deal with this problem. The LPFs provide the average powers consumed by the load Pav and Qav and removes the double frequency components

Materials and Methods
The scheme of a single-phase voltage source inverter (VSI) operated with the droop method when sharing a nonlinear load Z NL with another inverter is depicted in Figure 1. From the scheme, the main parts considered in this work can be seen: the power stage with inner controller, the pulse width modulation block (PWM), the LC output filter, the vo and io sensing, the Pav and Qav calculation block, and the droop generator producing the voltage reference νref [24,25]. In Figure 1, the scheme of the second inverter is not depicted. Rather, it is only outlined with discontinuous lines and indicated as #Inv. 2, due to the fact that this proposal is only concerned with the dynamic behavior and accuracy of the P-Q calculation block. For this reason, the simulations and experimental results shown from now on will only correspond to a single VSI.
In the simulations and experiments, a diode bridge rectifier (DBR) supplying a RC load is used as a nonlinear load Z NL (see Figure 2). At steady state, Z NL draws a distorted and symmetrical nonlinear current with peak levels reaching ±2.48A for the Z NL specified in Table 1, which induces a high distortion in Pav and Qav, as well as in v re f . Likewise, a switch S1 is inserted in the Z NL scheme (see Figure 2), allowing for step-perturbations for testing the dynamic behavior of Pav and Qav and for assessment and comparison purposes. In the simulations and experiments, a diode bridge rectifier (DBR) supplying a RC load is used as a nonlinear load ZNL (see Figure 2). At steady state, ZNL draws a distorted and symmetrical nonlinear current with peak levels reaching ±2.48A for the ZNL specified in Table 1, which induces a high distortion in Pav and Qav, as well as in . Likewise, a switch S1 is inserted in the ZNL scheme (see Figure 2), allowing for step-perturbations for testing the dynamic behavior of Pav and Qav and for assessment and comparison purposes. Simulations of the proposed system for the design and comparisons with [6] and [17] were performed with Matlab/Simulink/Simscape Power Systems software. Table 1 shows the main parameters of the inverter and ZNL. 960 Ω All the power calculation schemes were tuned to achieve the same power-ripple in order to obtain a fair comparison between the methods and to accurately measure the settling time of their transient responses when ZNL changes suddenly. The simulation results were obtained using Matlab/Simulink/Simscape Power Systems software, and were contrasted with HIL results at the In the simulations and experiments, a diode bridge rectifier (DBR) supplying a RC load is used as a nonlinear load ZNL (see Figure 2). At steady state, ZNL draws a distorted and symmetrical nonlinear current with peak levels reaching ±2.48A for the ZNL specified in Table 1, which induces a high distortion in Pav and Qav, as well as in . Likewise, a switch S1 is inserted in the ZNL scheme (see Figure 2), allowing for step-perturbations for testing the dynamic behavior of Pav and Qav and for assessment and comparison purposes. Simulations of the proposed system for the design and comparisons with [6] and [17] were performed with Matlab/Simulink/Simscape Power Systems software. Table 1 shows the main parameters of the inverter and ZNL. 960 Ω All the power calculation schemes were tuned to achieve the same power-ripple in order to obtain a fair comparison between the methods and to accurately measure the settling time of their transient responses when ZNL changes suddenly. The simulation results were obtained using Matlab/Simulink/Simscape Power Systems software, and were contrasted with HIL results at the  Simulations of the proposed system for the design and comparisons with [6,17] were performed with Matlab/Simulink/Simscape Power Systems software. Table 1 shows the main parameters of the inverter and Z NL .
All the power calculation schemes were tuned to achieve the same power-ripple in order to obtain a fair comparison between the methods and to accurately measure the settling time of their transient responses when Z NL changes suddenly. The simulation results were obtained using Matlab/Simulink/Simscape Power Systems software, and were contrasted with HIL results at the inverter-based intelligent Microgrid Laboratory (iML) of the department of Energy Technology at the Aalborg University (iML-AAU) in Denmark ( Figure 3). inverter-based intelligent Microgrid Laboratory (iML) of the department of Energy Technology at the Aalborg University (iML-AAU) in Denmark.
(a) (b) A second similar test was performed with another diode bridge rectifier nonlinear load drawing an asymmetrical nonlinear current in order to further test the calculation power block. This current reached a peak of 4.2 A in the positive half-cycle of the VSI and a negative one of −2.2 A in the VSI negative half-cycle after an abrupt load change. This asymmetry in the current introduced an extra distortion to the calculated powers. Experimental results for this load with a VSI inverter of the iML-AAU were also obtained. The inverter used in the iMG is a Danfoss© FC302, 2.2 kW rated, interfaced to a real-time dSPACE 1006. The algorithms for operating the VSI are developed in Matlab/Simulink software and compiled into the dSPACE. Figure 4 depicts the scheme of this experimental setup.

Droop Principle-Based Control Scheme.
The power calculation block of Figure 1 provides the frequency and voltage amplitude references * and * , respectively, as can be seen in the following equations: A second similar test was performed with another diode bridge rectifier nonlinear load drawing an asymmetrical nonlinear current in order to further test the calculation power block. This current reached a peak of 4.2 A in the positive half-cycle of the VSI and a negative one of −2.2 A in the VSI negative half-cycle after an abrupt load change. This asymmetry in the current introduced an extra distortion to the calculated powers. Experimental results for this load with a VSI inverter of the iML-AAU were also obtained. The inverter used in the iMG is a Danfoss© FC302, 2.2 kW rated, interfaced to a real-time dSPACE 1006. The algorithms for operating the VSI are developed in Matlab/Simulink software and compiled into the dSPACE. Figure 4 depicts the scheme of this experimental setup. inverter-based intelligent Microgrid Laboratory (iML) of the department of Energy Technology at the Aalborg University (iML-AAU) in Denmark.
(a) (b) A second similar test was performed with another diode bridge rectifier nonlinear load drawing an asymmetrical nonlinear current in order to further test the calculation power block. This current reached a peak of 4.2 A in the positive half-cycle of the VSI and a negative one of −2.2 A in the VSI negative half-cycle after an abrupt load change. This asymmetry in the current introduced an extra distortion to the calculated powers. Experimental results for this load with a VSI inverter of the iML-AAU were also obtained. The inverter used in the iMG is a Danfoss© FC302, 2.2 kW rated, interfaced to a real-time dSPACE 1006. The algorithms for operating the VSI are developed in Matlab/Simulink software and compiled into the dSPACE. Figure 4 depicts the scheme of this experimental setup.

Droop Principle-Based Control Scheme.
The power calculation block of Figure 1 provides the frequency and voltage amplitude references * and * , respectively, as can be seen in the following equations:

Droop Principle-Based Control Scheme
The power calculation block of Figure 1 provides the frequency and voltage amplitude references ω * and V * , respectively, as can be seen in the following equations: where m and n are the droop coefficients and ω n and V n are the nominal frequency and voltage amplitude, respectively. These references are used to generate the sinusoidal voltage reference v re f for the inverter inner control loops: v re f = V * ·sin(ω * ·t) Assuming that the output voltage and current of the inverter are [15] v o (t) = V·sin(ω o ·t) where V and I are the voltage and current amplitudes, ω o is the fundamental frequency and ϕ o is the phase angle between v o (t) and i o (t), the quadrature voltage, with a −π/2 delay, is defined as Therefore, the instantaneous active and reactive powers can be formulated as where p and q are the oscillating components that pulsate at twice the fundamental frequency ω o . These equations reveal that for a linear load that draws a sinusoidal current, the instantaneous powers oscillate around the averaged powers P av and Q av . However, if the load is nonlinear, the instantaneous powers are going to be highly distorted by harmonics, for which the nonlinear current can be expressed as: leading to where the subscript h represents the harmonic order, N is the maximum considered value for h and I h , h·ω o and ϕ h are the amplitude, the frequency and the phase-shift of the current harmonic components, respectively. As can be seen in Equations (10) and (11), p i (t) and q i (t) contain higher harmonic order components, in addition to the DC P av and Q av and the double frequency components p and q that were already present in the linear case. Then, for a nonlinear Z NL , Equations (1) and (2) can be expressed as follows when they are calculated employing the instantaneous powers in Equations (10) and (11): and the voltage reference is v re f (t) = V * (t)·sin(ω * (t)·t) (14) Figure 5 shows the conventional P-Q calculation scheme for obtaining P av and Q av , in which the instantaneous powers p i (t) and q i (t) are obtained as the products of i o (t) with v o (t) and v o⊥ (t), respectively [6]. The quadrature voltage v o⊥ (t) is obtained by delaying v o (t) by π/2, and P av and Q av by using LPFs with a low cut-off frequency value fc, to filter the multiple frequency components in the instantaneous powers [18]. When sharing linear loads, the value of fc is usually set to one or two orders of magnitude lower than the inverter fundamental operating frequency [22,23], which determines the transient dynamic performance of the system. However, in the case of nonlinear loads, the value of fc should be further reduced (usually to less than 1 Hz), to avoid strong distortions in the inverter output current and in the instantaneous powers [24]. Conversely, the distortion in the current induces excessive ripple in the averaged powers, which in turn is translated to the droop references ω * (t) and V * (t), and then to v re f (t), causing bad operation of the system. Nevertheless, this bandwidth reduction slows down the transient dynamic behavior of the system. The transfer functions of the averaged P av and Q av are shown in Appendix A. and by using LPFs with a low cut-off frequency value fc, to filter the multiple frequency components in the instantaneous powers [18]. When sharing linear loads, the value of fc is usually set to one or two orders of magnitude lower than the inverter fundamental operating frequency [22,23], which determines the transient dynamic performance of the system. However, in the case of nonlinear loads, the value of fc should be further reduced (usually to less than 1 Hz), to avoid strong distortions in the inverter output current and in the instantaneous powers [24]. Conversely, the distortion in the current induces excessive ripple in the averaged powers, which in turn is translated to the droop references * ( ) and * ( ), and then to ( ), causing bad operation of the system.

The SOGI and DSOGI Approach
A SOGI is a special linear filter with one input, vin (t), and two outputs, vd (t) and vq(t), one in phase and the other delayed π/2 with respect to the input, respectively, see Figure 6. The outputs of the SOGI filter, vd(t) and vq(t), have the following band-pass filter (BPF) and LPF transfer function relationships regarding to the input: where ξi is the filter damping factor and ωi is the tuning center frequency. In addition, a DSOGI is a four-order filter that consist in the cascade connection of two SOGI filters, as seen in Figure 7 [26]:

The SOGI and DSOGI Approach
A SOGI is a special linear filter with one input, v in (t), and two outputs, v d (t) and v q (t), one in phase and the other delayed π/2 with respect to the input, respectively, see Figure 6. and by using LPFs with a low cut-off frequency value fc, to filter the multiple frequency components in the instantaneous powers [18]. When sharing linear loads, the value of fc is usually set to one or two orders of magnitude lower than the inverter fundamental operating frequency [22,23], which determines the transient dynamic performance of the system. However, in the case of nonlinear loads, the value of fc should be further reduced (usually to less than 1 Hz), to avoid strong distortions in the inverter output current and in the instantaneous powers [24]. Conversely, the distortion in the current induces excessive ripple in the averaged powers, which in turn is translated to the droop references * ( ) and * ( ), and then to ( ), causing bad operation of the system.

The SOGI and DSOGI Approach
A SOGI is a special linear filter with one input, vin (t), and two outputs, vd (t) and vq(t), one in phase and the other delayed π/2 with respect to the input, respectively, see Figure 6. The outputs of the SOGI filter, vd(t) and vq(t), have the following band-pass filter (BPF) and LPF transfer function relationships regarding to the input: where ξi is the filter damping factor and ωi is the tuning center frequency. In addition, a DSOGI is a four-order filter that consist in the cascade connection of two SOGI filters, as seen in Figure 7 [26]: The outputs of the SOGI filter, v d (t) and v q (t), have the following band-pass filter (BPF) and LPF transfer function relationships regarding to the input: where ξ i is the filter damping factor and ω i is the tuning center frequency. In addition, a DSOGI is a four-order filter that consist in the cascade connection of two SOGI filters, as seen in Figure 7 [26]:  The transfer functions of the DSOGI with respect to the input [26]: When ωi is tuned to match , and for harmonic components higher than the fundamental (h ≫ 1), the gain of the BPF characteristic in Equation (17) can be simplified to: On the other hand, the frequency and damping factor are the parameters that determine the settling time (2% criterion) of the transient response of the BPF in Equation (15) for a step input: The SOGI filter has an inherent trade-off relationship between bandwidth (rejection capability to harmonics) and settling time response, Equation (20). This trade-off cannot be overcome, and relies on the damping factor parameter . However, in [26], it was shown that the DSOGI has a better trade-off than the SOGI and can achieve a faster transient response when it is designed to have the same bandwidth behavior. In this paper, this characteristic is used to achieve a faster response when extracting the fundamental component of the nonlinear load current and thereby to improve the droop transient response to load changes.  (10) and (11). This figure does not show the method for generating the π/2 delay, since it is not mentioned in [17]. Therefore, another SOGI, SOGI0, tuned at and = = 0.707, is used for generating this delay, as shown in Figure 9, for avoiding the delay issues reported in [15,16]. The transfer functions of the DSOGI with respect to the input [26]:

Advanced P-Q Calculation Scheme.
When ω i is tuned to match ω o , and for harmonic components higher than the fundamental (h 1), the gain of the BPF characteristic in Equation (17) can be simplified to: On the other hand, the frequency and damping factor are the parameters that determine the settling time t s (2% criterion) of the transient response of the BPF in Equation (15) for a step input: The SOGI filter has an inherent trade-off relationship between bandwidth (rejection capability to harmonics) and settling time response, Equation (20). This trade-off cannot be overcome, and relies on the damping factor parameter ξ i . However, in [26], it was shown that the DSOGI has a better trade-off than the SOGI and can achieve a faster transient response when it is designed to have the same bandwidth behavior. In this paper, this characteristic is used to achieve a faster response when extracting the fundamental component of the nonlinear load current and thereby to improve the droop transient response to load changes. Figure 8 shows a P-Q calculation method based on the proposed scheme in [17] for accelerating the computation of P av and Q av . In this figure, the SOGI1 and SOGI2 blocks are used for extracting the double frequency pulsating power components p and q, which are then removed from p i and q i . These SOGIs are both tuned at ω i = 2ω o and ξ 1 = ξ 2 = 1. The LPFs are used for improved filtering and for providing the averaged powers P av and Q av , by attenuating the higher harmonics components reported in Equations (10) and (11). This figure does not show the method for generating the π/2 delay, since it is not mentioned in [17]. Therefore, another SOGI, SOGI0, tuned at ω o and ξ 0 = ξ v = 0.707, is used for generating this delay, as shown in Figure 9, for avoiding the delay issues reported in [15,16].  The transfer functions for the calculated averaged active and reactive power are included in Appendix A, with the numbers (24) and (25). Please note that the current is not filtered. Figure 10 shows the simulation results using the P-Q algorithms of the schemes in Figures 5 and  9 when sharing a linear load that produces a current perturbation from peak 2A to peak 4A at a time of 3 s.  Figure 8. Scheme for the calculation of P av and Q av proposed in [20].  The transfer functions for the calculated averaged active and reactive power are included in Appendix A, with the numbers (24) and (25). Please note that the current is not filtered. Figure 10 shows the simulation results using the P-Q algorithms of the schemes in Figures 5 and  9 when sharing a linear load that produces a current perturbation from peak 2A to peak 4A at a time of 3 s. The transfer functions for the calculated averaged active and reactive power are included in Appendix A, with the numbers (24) and (25). Please note that the current is not filtered. Figure 10 shows the simulation results using the P-Q algorithms of the schemes in Figures 5 and 9 when sharing a linear load that produces a current perturbation from peak 2A to peak 4A at a time of 3 s.

Advanced P-Q Calculation Scheme
For the sake of simplicity, Figure 10 only shows the dynamic behavior of P av , obviating the representation of Q av . The dynamics of P av obtained with the advanced scheme in Figure 9 (hereinafter called P adv ) are compared with those obtained for P av by the conventional droop method depicted in Figure 5 (hereinafter called P droop ). The cut-off frequency of the LPFs in Figure 9 was designed to be fc = 3.7 Hz, whereas it was set to fc = 0.37Hz for the scheme in Figure 5. As shown in Figure 10, the double frequency component p was removed from P adv . Likewise, the higher cut-off frequency of its LPFs results in much faster dynamics than those of P droop . Also, the ripple corresponding to the double frequency component p, which is not completely filtered by the LPFs in Figure 5, can be observed in P droop . These results are compatible with those reported in [21]. However, the good dynamic behavior depicted in Figure 10 for P adv vanishes when a nonlinear load is used, as is shown in Figure 11. In this case, the nonlinear load is a DBR that draws a peak current of ±2.48 A and suffers a perturbation that pushes the peak to ±4.15 A. The simulation parameters are listed in Table 2.
The transfer functions for the calculated averaged active and reactive power are included in Appendix A, with the numbers (24) and (25). Please note that the current is not filtered. Figure 10 shows the simulation results using the P-Q algorithms of the schemes in Figures 5 and  9 when sharing a linear load that produces a current perturbation from peak 2A to peak 4A at a time of 3 s.  Table 2. Simulation parameters for Figure 10.

Parameter Value
As shown in Figure 11, the dynamics of the method proposed in [20] worsen using a nonlinear load, similarly to [26]. Thus, in the presence of a nonlinear load, the method has excessive steady state ripple that corrupts the calculated powers, opposite to what is stated in [20]. To reduce the ripple, the filtering capabilities of the LPFs in Figure 9 can be improved by reducing fc. By decreasing fc from 3.7 Hz to 1.1 Hz, the ripple of P adv matches that of P droop , as can be seen in Figure 12. Although the advanced method of Figure 9 calculates faster P av than the conventional droop controller, it presents less effectivity than initially argued. It can be clearly seen in Figures 11 and 12 that there is a trade-off between the filtering capability of the P av calculation scheme and the transient speed of its dynamical response. Note also that there is a positive offset in the calculated P av at steady state, since the mean value of P adv is slightly higher than that of P droop (see Figure 11). Comparing the transfer functions for the P-Q calculation by the conventional, Equations (22) and (23), or by the advanced method, Equations (24) and (25), it can be seen the filtering capabilities of each algorithm. Please note that the magnitude of the ripple of P adv forces the system's dynamic response to slow down by reducing more the cut-off frequency.

Parameter
Value 311V 2π50 rad/s DBR resistor load for t < 3s 950 Ω DBR resistor load for t > 3s 471.  As shown in Figure 11, the dynamics of the method proposed in [20] worsen using a nonlinear load, similarly to [26]. Thus, in the presence of a nonlinear load, the method has excessive steady state ripple that corrupts the calculated powers, opposite to what is stated in [20]. To reduce the ripple, the filtering capabilities of the LPFs in Figure 9 can be improved by reducing fc. By decreasing fc from 3.7 Hz to 1.1 Hz, the ripple of Padv matches that of Pdroop, as can be seen in Figure 12. Although the advanced method of Figure 9 calculates faster Pav than the conventional droop controller, it presents less effectivity than initially argued. It can be clearly seen in Figures 11 and 12 that there is a trade-off between the filtering capability of the Pav calculation scheme and the transient speed of its dynamical response. Note also that there is a positive offset in the calculated Pav at steady state, since the mean value of Padv is slightly higher than that of Pdroop (see Figure 11). Comparing the transfer functions for the P-Q calculation by the conventional, Equations (22) and (23), or by the advanced method, Equations (24) and (25), it can be seen the filtering capabilities of each algorithm. Please note that the magnitude of the ripple of Padv forces the system's dynamic response to slow down by reducing more the cut-off frequency.

Proposed P-Q Calculation Scheme.
The scheme of the proposed P-Q calculation is shown in Figure 13. A DSOGI approach is applied to ( ) in order to extract the fundamental component , see Equation (28) in Appendix A. There, the DSOGI composed of the SOGI3 and SOGI4 blocks filters the distorting high-order harmonics of ( ) by means of its higher BPF capability and avoids coping with a highly distorted current signal. Consequently, the instantaneous powers ( ) and ( ) are obtained as the product of the in-phase ( ) or the quadrature ⊥ ( ) voltages with the fundamental output current iOF, respectively. This produces a result with only double frequency components and without third or higher order harmonics.
Later, the SOGI1 and SOGI2 blocks were used as in Figure 9 for removing only the double frequency components with the help of the subtracting blocks. Therefore, the LPFs can now be removed from the scheme, since they are no longer necessary. The transfer function of the proposed scheme in Figure 13 is shown in Appendix A as Equations (26) and (27). This overcomes the main limitation of previous schemes and further accelerates the dynamic response of Pav and Qav. To distinguish Pav and Qav obtained with this proposed scheme, in the following they will be referred to as PDSOGI and QDSOGI.

Proposed P-Q Calculation Scheme
The scheme of the proposed P-Q calculation is shown in Figure 13. A DSOGI approach is applied to i o (t) in order to extract the fundamental component i OF , see Equation (28) in Appendix A. There, the DSOGI composed of the SOGI3 and SOGI4 blocks filters the distorting high-order harmonics of i o (t) by means of its higher BPF capability and avoids coping with a highly distorted current signal. Consequently, the instantaneous powers p i (t) and q i (t) are obtained as the product of the in-phase v o (t) or the quadrature v o⊥ (t) voltages with the fundamental output current i OF , respectively. This produces a result with only double frequency components and without third or higher order harmonics. In this case, because a DSOGI is used for filtering the current in Equation (28), and considering that the center frequencies provided by the droop method ω* vary in a small range around the nominal frequency ωn, the transient response speed is determined mainly by the DSOGI transient response, which is related to  ; see [26]. The damping factors of the DSOGI are here tuned to produce a power ripple identical in amplitude to that of the conventional droop controller, which is achieved for ξc = ξ3= ξ4 = 0.129. Later, the SOGI1 and SOGI2 blocks were used as in Figure 9 for removing only the double frequency components with the help of the subtracting blocks. Therefore, the LPFs can now be removed from the scheme, since they are no longer necessary. The transfer function of the proposed scheme in Figure 13 is shown in Appendix A as Equations (26) and (27). This overcomes the main limitation of previous schemes and further accelerates the dynamic response of P av and Q av . To distinguish P av and Q av obtained with this proposed scheme, in the following they will be referred to as P DSOGI and Q DSOGI .
In this case, because a DSOGI is used for filtering the current in Equation (28), and considering that the center frequencies provided by the droop method ω* vary in a small range around the nominal frequency ω n , the transient response speed is determined mainly by the DSOGI transient response, which is related to ξ c ; see [26]. The damping factors of the DSOGI are here tuned to produce a power ripple identical in amplitude to that of the conventional droop controller, which is achieved for ξ c = ξ 3 = ξ 4 = 0.129. Figure 14 shows the simulation results of the same nonlinear DBR load drawing a peak current of ±2.48 A and suffering a perturbation that pushes the peak to ±4.15 A. It can be observed that the proposed method is faster for calculating P av than the previously considered methods. In this case, because a DSOGI is used for filtering the current in Equation (28), and considering that the center frequencies provided by the droop method ω* vary in a small range around the nominal frequency ωn, the transient response speed is determined mainly by the DSOGI transient response, which is related to ξ ; see [26]. The damping factors of the DSOGI are here tuned to produce a power ripple identical in amplitude to that of the conventional droop controller, which is achieved for ξc = ξ3= ξ4 = 0.129. Figure 14 shows the simulation results of the same nonlinear DBR load drawing a peak current of ±2.48 A and suffering a perturbation that pushes the peak to ±4.15 A. It can be observed that the proposed method is faster for calculating Pav than the previously considered methods.  Figure 15 depicts the simulation results using a Z NL drawing an asymmetrical current, reaching a positive peak of 5.21 A and a negative peak of −3.07 A after the perturbation. In this case, the asymmetry in the nonlinear current (see Figure 15a), introduces further distortion into the calculated P av , showing higher ripple at steady state (see Figure 15c).  Figure 15 depicts the simulation results using a ZNL drawing an asymmetrical current, reaching a positive peak of 5.21 A and a negative peak of −3.07 A after the perturbation. In this case, the asymmetry in the nonlinear current (see Figure 15a), introduces further distortion into the calculated Pav, showing higher ripple at steady state (see Figure 15c). Please note that unlike Padv, PDSOGI does not exhibit a positive offset error in steady state. This means that the proposed method is also more accurate than that in [20]. Table 3 contains the measured settling time for the Pav transient responses depicted in Figure 15, which shows that the PDSOGI calculated with the proposed method implies a 60.00% and a 79.69% reduction in the response time regarding Padv and Pdroop, respectively.  Please note that unlike P adv , P DSOGI does not exhibit a positive offset error in steady state. This means that the proposed method is also more accurate than that in [20]. Table 3 contains the measured settling time for the P av transient responses depicted in Figure 15, which shows that the P DSOGI calculated with the proposed method implies a 60.00% and a 79.69% reduction in the response time regarding P adv and P droop , respectively.

Experimental Results
Hardware in Loop, as well as experimental results with an VSI inverter of the iML-AAU, were obtained. The HIL setup consisted of a dSpace 1006© digital Real-Time Interface platform that carried out the processor-based simulations with parameters configured by the dSpace Configuration Desk © tool. The setup supports models based on the physical modelling libraries of electrical plants designed by Matlab/Simulink/SimPowerSystems©. These models are integrated into the dSpace tool chain and are used for building the electrical system under test along with the Electronic Central Unit, ECU, of the dSpace. The control algorithms described in this paper were loaded and executed in real time in the ECU. In this case, the H-Bridge Inverter, the LCL filter and the Nonlinear Load depicted in Figure 4 are the electrical system implemented in the ECU under the HIL test. The experimental setup consisted of an inverter Danfoss© FC302, 2.2 kW rated, interfaced to the real-time dSPACE 1006 digital platform. The algorithms for operating the VSI were developed in Matlab/Simulink software and compiled in the dSPACE multiprocessor core, with its parameters configured and controlled through the Configuration-Desk software. The sampling frequency for the system was 10 kHz, which was the switching frequency for the VSI. A third-order method for discretizing the SOGI and DSOGI algorithms was employed, whereby the integrator 1 s was approximated as: with Ts being the sample time of 100 µs, which is consistent with the 10 kHz frequency. The experimental setup parameters for the LCL filter and for the nonlinear load are listed in Table 4. In this section, HIL results are shown first for a Z NL drawing a symmetrical current, and then the results for a VSI with a Z NL drawing an asymmetrical current. Figure 16 depicts the HIL P av results for a Z NL drawing a symmetrical current with peak values transitioning from ±2.48 A to ±4.15 A after a step perturbation at 1.28 s. The dynamics of P av obtained with the different considered methods are represented in green for P droop , in blue for P adv and in red for P DSOGI . Figure 17 shows a detail of the P av ripple waveforms at steady state, evidencing the different nonlinearities resulting from each power calculation method. There are some DC errors in the calculated powers that can mainly be attributed to the discretization method in Equation (21) and to the sampling frequency. This phenomenon was not evidenced in the simulations shown in Figure 15c, but is reflected in the HIL results.   Figure 17 shows a detail of the Pav ripple waveforms at steady state, evidencing the different nonlinearities resulting from each power calculation method. There are some DC errors in the calculated powers that can mainly be attributed to the discretization method in Equation (21) and to the sampling frequency. This phenomenon was not evidenced in the simulations shown in Figure  15c, but is reflected in the HIL results. Figure 18 shows the output voltage vo(t) waveform of the VSI (Figure 18a), the current waveform drawn for the DBR load and used for obtaining the VSI results ( Figure 18b) and the fundamental current component iOF extracted by the DSOGI (Figure 18c). Note the current asymmetry, with positive peaks that reach 4.2 A and negative ones that reach −2.8 A. Note how a fundamental current component t of 1.5 A amplitude has been achieved. Figure 16. HIL P av results for a Z NL drawing a symmetrical current with peak values from ±2.48 A to ±4.15 A after a step perturbation at 1.28 s (P droop in green, P adv in blue, P DSOGI in red).   Figure 17 shows a detail of the Pav ripple waveforms at steady state, evidencing the different nonlinearities resulting from each power calculation method. There are some DC errors in the calculated powers that can mainly be attributed to the discretization method in Equation (21) and to the sampling frequency. This phenomenon was not evidenced in the simulations shown in Figure  15c, but is reflected in the HIL results. Figure 18 shows the output voltage vo(t) waveform of the VSI (Figure 18a), the current waveform drawn for the DBR load and used for obtaining the VSI results ( Figure 18b) and the fundamental current component iOF extracted by the DSOGI (Figure 18c). Note the current asymmetry, with positive peaks that reach 4.2 A and negative ones that reach −2.8 A. Note how a fundamental current component t of 1.5 A amplitude has been achieved. Figure 17. Detail of HIL P av ripples at steady-state from Figure 14 (P droop in green, P adv in blue, P DSOGI in red). Figure 18 shows the output voltage v o (t) waveform of the VSI (Figure 18a), the current waveform drawn for the DBR load and used for obtaining the VSI results ( Figure 18b) and the fundamental current component i OF extracted by the DSOGI (Figure 18c). Note the current asymmetry, with positive peaks that reach 4.2 A and negative ones that reach −2.8 A. Note how a fundamental current component t of 1.5 A amplitude has been achieved.
Finally, Figure 19 depicts the P av experimental VSI results with the DBR load when a perturbation is applied at 8.9 s. Note in this figure that the results are coherent with those obtained by HIL, as shown in Figure 16, and by the simulations in Figure 15. Table 5 presents the settling time of P av obtained with all the considered P-Q calculation methods, as measured in Figure 15 for Matlab/Simulink ® simulations results, in Figure 16 for the HIL results, and in Figure 19 for the experimental VSI results. In this table, the achieved percentage of settling time reduction with respect to the conventional droop method is also indicated. Note how the proposed method achieves a faster response. The transfer functions for the active and reactive averaged powers of the droop, advanced and DSOGI methods are shown in Appendix A. Finally, Figure 19 depicts the Pav experimental VSI results with the DBR load when a perturbation is applied at 8.9 s. Note in this figure that the results are coherent with those obtained by HIL, as shown in Figure 16, and by the simulations in Figure 15.  Table 5 presents the settling time of Pav obtained with all the considered P-Q calculation methods, as measured in Figure 15 for Matlab/Simulink ® simulations results, in Figure 16 for the HIL results, and in Figure 19 for the experimental VSI results. In this table, the achieved percentage of settling time reduction with respect to the conventional droop method is also indicated. Note how the proposed method achieves a faster response. The transfer functions for the active and reactive averaged powers of the droop, advanced and DSOGI methods are shown in Appendix A.  Finally, Figure 19 depicts the Pav experimental VSI results with the DBR load when a perturbation is applied at 8.9 s. Note in this figure that the results are coherent with those obtained by HIL, as shown in Figure 16, and by the simulations in Figure 15.  Table 5 presents the settling time of Pav obtained with all the considered P-Q calculation methods, as measured in Figure 15 for Matlab/Simulink ® simulations results, in Figure 16 for the HIL results, and in Figure 19 for the experimental VSI results. In this table, the achieved percentage of settling time reduction with respect to the conventional droop method is also indicated. Note how the proposed method achieves a faster response. The transfer functions for the active and reactive averaged powers of the droop, advanced and DSOGI methods are shown in Appendix A. Figure 19. VSI experimental result of the P droop transient response for the nonlinear DBR load current perturbation in Figure 18 (P droop in green, P adv in blue, P DSOGI in red).

Discussion
In this work, a P-Q calculation method was proposed for single-phase inverters with the purpose of improving the speed and accuracy of the power calculation when they are sharing linear and nonlinear loads. The dynamic response of the power calculation used in the conventional droop method and in another advanced method was first analyzed to show their limitations in speed and accuracy when sharing a DBR-type nonlinear load. For this reason, a novel calculation method was proposed, bearing in mind the dynamic transient response velocity of the system to sudden perturbations in the shared load. The method was studied and compared with previous ones by obtaining results from Matlab ® -based simulations, the HIL platform, and a VSI experimental setup. These results show how the proposed method achieves, under the same distorting conditions, a faster calculation settling time with a measured time reduction over the conventional droop method that arrives around 80%, which is higher than the achieved by the advanced method. This improvement supposes an enhancement in the droop speed operation under linear and nonlinear load conditions that leads to a better dynamic performance of the system when parallelizing inverters or microgrid operation in islanded mode.
The aim of this work was to identify what the main limitation of the power calculation methods was for sharing linear and nonlinear loads, and to propose a new approach based on the DSOGI scheme. This issue can be further investigated by using the SOGI approach and other ones for obtaining the fundamental component of the distorted current signals with less effort, faster, and best quality. This would help in the operation and stability of the inverters when sharing linear and nonlinear loads. Acknowledgments: Thanks to Alexander Micallef, University of Malta, for his crucial support during the experimental campaign. Thanks also to the E3PACS laboratory staff from UPC, that provided technical support during the redaction of this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The transfer functions of the active and reactive averaged powers regarding the instantaneous powers of the methods described here are as follows: where ω c = 2π f c is the cut-off frequency in rad/s of the LPFs used in the droop-and advanced-based methods. The instantaneous powers, p i (t) and q i (t), are derived from the product between i o (t) and v o (t) and v o⊥ (t), respectively. The transfer functions for v o (t) and v o⊥ (t) correspond to Equations (15)