Analysis and Design Guidelines for Current Control Loops of Grid-Connected Converters Based on Mathematical Models

: Having a method for analyzing and designing regulators of controls that contain many current loops such as active ﬁlters is not a trivial task. There can be many parameters of regulators and ﬁlters that must be carefully selected in order to fulﬁll certain desired requirements. For instance, these can be stability, dynamic response, robustness under uncertainty of parameters, and rejection capability to switching harmonics. Hence, this paper provides general analysis guidelines for designing current control loops by using mathematical models in an αβ reference frame. Then, by using the proposed modeling tool, a multi-objective tuning algorithm is proposed that helps obtain all the control loops’ regulator and ﬁlter parameters, meeting all the desired requirements. Thus, the proposed analysis and design methodology is illustrated by applying it to three di ﬀ erent controls conceived in a dq rotating reference frame with PI (Proportional Integral) regulators. The ﬁrst control presents two current loops (simple dq current control), the second control uses four current loops (dual vector control, for unbalanced loads), while the third control presents eight current loops (active ﬁlter controlling current harmonics). Several experimental and simulation results show the e ﬀ ectiveness and usefulness of the proposed method. Since the mathematical model employed is in the αβ reference frame, it can also be easily applied to controls conceived in a αβ reference frame using resonant regulators, providing also a common comparative framework.


Introduction
It can be affirmed that nowadays, either in industry or research, most of the current vector controls of three-phase grid-connected converters are made by using PI (proportional integral) controllers in synchronously rotating reference frame (dq) or with resonant controllers in a stationary reference frame (αβ) [1]. In addition, it can be observed that there is a clear and natural tendency to use mathematical models for analysis in a stationary reference frame (αβ), obviously in controls with resonant controllers [2]. While in controls with PI controllers, the mathematical models for analysis are more often developed in a synchronous rotating reference frame (dq), even if the final resulting or real currents exchanged between the grid and the power electronic converter are actually sinusoidal [3,4].
Recently, in [5], it has been shown that it is also possible and useful to develop mathematical models in a stationary reference frame (αβ) to analyze controls in a synchronous rotating reference frame (dq) with PI controllers. Although the modeling philosophy is not new [6,7], this modeling technique can be useful due to the following main reasons: it allows performing directly the analysis of the actual sinusoidal currents (stationary reference frame αβ) exchanged between the converter Figure 1. Grid-connected converter operating with current references in the dq reference frame; id*(t) and iq*(t) and simplified equivalent electric circuit in the αβ reference frame.
The control block diagram is depicted in Figure 2 [13]. The current control is effectively implemented in a rotating reference frame (dq) by using two PI controllers. As commented before, in most of the applications, the current references id*(t) and iq*(t) are created in a rotating reference frame (dq); therefore, in this case, the control is also implemented in a dq reference frame rotating at ω angular speed. The structure of the control is quite standard and often used in many applications [1,3,4]. As a not so common feature, it can be seen that optional low-pass filters have been included at both voltage and current measurements. Depending on the numerical scenario where this control is going to be used, if voltage and current ripples due to harmonics are significant at both measurements, it may be necessary to attenuate the ripple level so it does not interfere in the controls performance. Then, also in Figure 1, the αβ equivalent power circuit is depicted [3]. Note that it is composed by the grid voltage v gα (t) and v gβ (t) that in a first analysis of this paper is supposed to be purely sinusoidal. The equivalent impedances R g and L g represent the sum of the series impedances of the transformer (leakages impedances), the equivalent impedance of the PCC, and the Thevenin equivalent impedance of the grid. Then, the equivalent impedances R and L of the filter are also present together with the C c and damping resistance R c , while finally, the voltages created by the converter are also present. These voltages v convα (t) and v convβ (t) are automatically generated by the corresponding current control method employed. As justified before, this paper derives all the mathematical expressions in a αβ stationary reference frame. Note that the three components could be also used, abc; however, the results and conclusions would be equivalent.
The control block diagram is depicted in Figure 2 [13]. The current control is effectively implemented in a rotating reference frame (dq) by using two PI controllers. As commented before, in most of the applications, the current references i d *(t) and i q *(t) are created in a rotating reference frame (dq); therefore, in this case, the control is also implemented in a dq reference frame rotating at ω angular speed. The structure of the control is quite standard and often used in many applications [1,3,4]. As a not so common feature, it can be seen that optional low-pass filters have been included at both voltage and current measurements. Depending on the numerical scenario where this control is going * * 11 12

Mathematical Model of the Power Circuit and Control
The mathematical model of the power circuit and the control algorithm in the αβ reference frame is developed in Table 1. As it is seen, first of all, the model of the control is derived, starting from the measurements. Both voltage and currents are measured and passed through two low-pass filters F v (z) and F i (z); then, the phase shift alteration by the filters at the fundamental frequency is compensated. Note that both measured grid voltages and currents are needed to be applied as feed-forward terms [5]. Then, the different mathematical terms that constitute the control action are derived: C αα (s), C αβ (s), C βα (s), C ββ (s), C αd (s), C αq (s), C βd (s), C βq (s) in the αβ reference frame. For a deeper understanding of the mathematical development, it is possible to study the references [4,5,7]. Then, all these terms are discretized by the desired method (zero-order-hold ZOH in the present paper, with T s sample time). Once these terms are discretized, the v regαβ (z) voltages are first obtained and then the v controlαβ (z) terms. Finally, the v convαβ (z) output voltages provided by the control are obtained. In this case, the converter is modeled simply by an update delay (see Section 2.4 for more details) and its corresponding compensation at the fundamental frequency [20]. Once the converter voltages are derived, the filter and grid model equations are used to derive the currents i αβ (z) and voltages v RCαβ (z). To conclude, all the partial elements calculated must be connected in order to obtain the global model matrix equation G mm (z) that relates the output electric magnitudes, i αβ (z) and v RCαβ (z), in function of the inputs of the system, i dq *(z) and v gαβ (z). This last step is solved numerically by using the "connect" function from Matlab. The obtained closed loop global system by applying the numerical scenario of Appendix A by using the "connect" function from Matlab is represented in state-space form, in which "A" is a matrix of 16 × 16, "B" is a matrix of 16 × 4 and "C" is a matrix of 4 × 16. This order of matrices are obtained with the control structure of Figure 2 but not including the voltage and current filters and assuming also that there is not delay in the converter voltage update, i.e., F i (z) = 1, F v (z) = 1, Del(z) = 1. Table 1. Mathematical model of the control of Figure 2. It can be directly implemented sequentially in a Matlab program. c2d: continuous to discrete. The discretization method used in this paper is ZOH.

Summary of the Analysis Carried out with the Assistance of the Mathematical Model
In the subsequent sections, the following analyses are carried out. First of all, the dominant poles that determine the dynamic performance of the control are identified. Then, the most important parameters that determine the location of each pole (dominant and others) are found. After that, how the converter delay affects the location of the poles and therefore the dynamic behavior and stability of the system is studied. Then, we investigate how the current filter affects the poles' location and whether this filter is strictly necessary or not. Finally, we ascertain how the voltage filter affects the poles' location and when this filter is strictly necessary. At the end, the proposed mathematical model is validated by means of simulation and experiments.

Analysis of the Location of the Poles
Once the mathematical equations of the entire system are obtained, it is possible to analyze the locations of the poles to better understand the stability of the system and the dominant dynamics that defining the time domain responses. This section includes an example showing how the poles are located as a function of the numerical values of the system. The k p and T n gains of the PI regulators are tuned according the formulas described in [3] that consider the system in a dq reference frame and the closed loop current loops are ideally simplified to the following expressions: T n /L s 2 + 2ξ cl ω cl s + ω cl 2 (6) ω cl and ξ cl are chosen by the designer, as they are the desired natural frequency and the damping ratio of the closed loop time domain step response. Note that as these equations neglect several realities of the system (converter delay, filter measurements, etc.), the actually obtained closed loop behavior defers slightly from Equation (6).
However, it is a reasonably good and simple first approach to the tuning of the system. Then, once the k p and T n gains of the PI regulators are obtained in a first step, the αβ model presented in Table 1 is applied to the numerical scenario of Appendix A, obtaining thus the poles and zeros of the total closed loop system. This can be made for instance, by using the "pzmap" function from Matlab or any other equivalent software tool. Hence, Figure 3 depicts the poles and zeros of the  Table 1 applied to the numerical scenario of Appendix A. It is seen that there are five pole pairs. One pole pair is at the unit circle, so it is responsible for the steady-state current (i αβ currents). realities of the system (converter delay, filter measurements, etc.), the actually obtained closed loop behavior defers slightly from Equation (6).
However, it is a reasonably good and simple first approach to the tuning of the system. Then, once the kp and Tn gains of the PI regulators are obtained in a first step, the αβ model presented in Table 1 is applied to the numerical scenario of Appendix A, obtaining thus the poles and zeros of the total closed loop system. This can be made for instance, by using the "pzmap" function from Matlab or any other equivalent software tool. Hence, Figure 3 depicts the poles and zeros of the global closed loop system, corresponding to Figure 2, which are deduced from the mathematical model of Table 1 applied to the numerical scenario of Appendix A. It is seen that there are five pole pairs. One pole pair is at the unit circle, so it is responsible for the steady-state current (iαβ currents).

Figure 3.
Poles and zeros of i α (z)/i d *(z) and dependence on the parameters in the αβ reference frame for step input references, corresponding to the control of Figure 2, deduced from the mathematical model of Table 1, and applied to the numerical scenario of Appendix A. The voltage and current filters are not included, and it is assumed that there is no delay in the converter voltage update, i.e., F i (z) = 1, The poles for the rest of the controlled input-output relations are equal ( The remaining four pole pairs are responsible for the transient response, and as they are within the unit circle, they are damped. From these four pole pairs, it is seen that the dominant two pole pairs (poles with a bigger module, or with slower time response) depend mainly on the k p and T n gains of the PI regulators. The other remaining two pole pairs depend on mainly L, L g , C c , and R c . Sample time T s affects the four damped pole pairs. The dependence of each pole on each numerical parameter of the system can be deduced in several ways. In this article, it is simply done by varying the parameters one by one within a given limited range and plotting the pole-zero map as depicted in Figure 4. As can be seen, variation of some of the parameters makes some of the poles tend to be more stable and some others more instable (further to the unit circle). Note that the natural frequency of the non-dominant two pole pairs depends on the filter and grid parameters, together with the gains of the PI regulator. Therefore, the control itself and even a potentially changing grid impedance makes this frequency vary.

Study of the Converter Delay Effect (Current and Voltage Filters Not Included)
Once the basic and ideal system is analyzed in terms of its poles, this subsection studies the effect of the converter's delay effect. This update delay due to the modulation process is graphically represented in Figure 5. A classical asymmetrical regular sampled PWM is used [20]; thus, the reference is updated in the lower and upper peaks of the carrier waveform. In this way, a higher updated ratio and, therefore, lower delay is achieved compared to the symmetrical regularly sampled PWM (Pulse Width Modulation).
Energies 2020, 13, 5849 9 of 47 the parameters one by one within a given limited range and plotting the pole-zero map as depicted in Figure 4. As can be seen, variation of some of the parameters makes some of the poles tend to be more stable and some others more instable (further to the unit circle). Note that the natural frequency of the non-dominant two pole pairs depends on the filter and grid parameters, together with the gains of the PI regulator. Therefore, the control itself and even a potentially changing grid impedance makes this frequency vary. Figure 4. Variation of the poles and zeros of i α (z)/i d *(z) in the αβ reference frame for step input references, corresponding to the control of Figure 2, deduced from the mathematical model of Table 1, and applied to the numerical scenario of Appendix A. Each parameter's value is modified within a range. Voltage and current filters are not included, and it is assumed that there is no delay in the converter voltage update, i.e., Energies 2020, 13, x FOR PEER REVIEW 8 of 42 Figure 4. Variation of the poles and zeros of iα(z)/id*(z) in the αβ reference frame for step input references, corresponding to the control of Figure 2, deduced from the mathematical model of Table  1, and applied to the numerical scenario of Appendix A. Each parameter's value is modified within a range. Voltage and current filters are not included, and it is assumed that there is no delay in the converter voltage update, i.e., Fi(z) = 1, Fv(z) = 1, Del(z) = 1 and ϕi = ϕv = ϕd = 0. (a) kp variation; (b) Tn variation; (c) L variation; (d) Lg variation; (e) C variation; (f) Rc variation.

Study of the Converter Delay Effect (Current and Voltage Filters Not Included)
Once the basic and ideal system is analyzed in terms of its poles, this subsection studies the effect of the converter's delay effect. This update delay due to the modulation process is graphically represented in Figure 5. A classical asymmetrical regular sampled PWM is used [20]; thus, the reference is updated in the lower and upper peaks of the carrier waveform. In this way, a higher updated ratio and, therefore, lower delay is achieved compared to the symmetrical regularly sampled PWM (Pulse Width Modulation). Therefore, the present analysis is focused on the previously presented control block diagram of Figure 2, where filters of voltage and current measurements are not included (Fi(z) = 1, Fv(z) = 1), but the update delay of the converter and its compensation as proposed in [21] is modeled (widely adopted as a solution to mitigate the negative effect on the control performance of this delay). Consequently, the mathematical representation of the delay is included in the system model   Therefore, the present analysis is focused on the previously presented control block diagram of Figure 2, where filters of voltage and current measurements are not included (F i (z) = 1, F v (z) = 1), but the update delay of the converter and its compensation as proposed in [21] is modeled (widely adopted as a solution to mitigate the negative effect on the control performance of this delay). Consequently, the mathematical representation of the delay is included in the system model equations, as presented in Table 1. The resulting pole-zero maps are represented in Figure 6. Therefore, the present analysis is focused on the previously presented control block diagram of Figure 2, where filters of voltage and current measurements are not included (Fi(z) = 1, Fv(z) = 1), but the update delay of the converter and its compensation as proposed in [21] is modeled (widely adopted as a solution to mitigate the negative effect on the control performance of this delay). Consequently, the mathematical representation of the delay is included in the system model equations, as presented in Table 1. The resulting pole-zero maps are represented in Figure 6.  Table 1, and applied to the numerical scenario of Appendix A. The voltage and current filters are not included, i.e., Fi(z) = 1, Fv(z) = 1 and ϕi = ϕv = 0. Blue: without delay, Red: with delay, Yellow: with delay and compensation. The poles for the rest of the controlled input-output relations are equal (iα(z)/iq*(z), A l g o r i t h m k+2 Figure 6. Poles and zeros of i α (z)/i d *(z) in the αβ reference frame for step input references, corresponding to the control of Figure 2, deduced from the mathematical model of Table 1, and applied to the numerical scenario of Appendix A. The voltage and current filters are not included, i.e., F i (z) = 1, F v (z) = 1 and φ i = φ v = 0. Blue: without delay, Red: with delay, Yellow: with delay and compensation. The poles for the rest of the controlled input-output relations are equal ( It is seen that the delay produces a considerable movement of the location of the poles. Then, the delay compensation [21] produces a reduction of the module (becomes more stable) of the dominant pole, slightly improving the stability of the system. However, it is seen that the compensation does not return to the pole placement of the original system without delay.

Study of the Current Measurement Filter Effect (Voltage Filters Not Included)
In this section, the effect of the filter at the current measurement is studied as included in the block diagram presented previously in Figure 2. As it is commonly done also, the phase shift produced by the current filter at the 50 Hz frequency is compensated (steady-state phase shift compensation). Consequently, including a second-order filter as follows (ω LP = 9425 rd/s): The pole-zero map of the entire closed loops system yields is presented in Figure 7. The current filter effect, as occurred with the update delay of the converter, displaces in general all the poles of the system, in which some of them (including a dominant one) become more unstable or closer to the unit circle. It can also be remarked that although it is not shown for simplicity in the exposition, if the phase shift angle of the filter is not compensated, the pole locations become slightly worse from the stability point of view.
In addition, many other filters have been tested (first order and different cut-off frequencies) reaching very similar results, concluding that the nature of the filter is not so crucial and affects the location of the poles quite similarly.
On the other hand, the current measurement filter is normally included to avoid the current ripple due to the switching behavior of the converter being transmitted through the control loops and deteriorating the current quality amplifying the ripple, or even in extreme cases, making the system become unstable.
The pole-zero map of the entire closed loops system yields is presented in Figure 7. The current filter effect, as occurred with the update delay of the converter, displaces in general all the poles of the system, in which some of them (including a dominant one) become more unstable or closer to the unit circle. It can also be remarked that although it is not shown for simplicity in the exposition, if the phase shift angle of the filter is not compensated, the pole locations become slightly worse from the stability point of view.  Table 1, and applied to the numerical scenario of Appendix A. The voltage filter is not included, i.e., Fv(z) = 1 and ϕi = 0. Blue: without current filter, Red: with current filter and its phase shift compensation at 50 Hz.
In addition, many other filters have been tested (first order and different cut-off frequencies) reaching very similar results, concluding that the nature of the filter is not so crucial and affects the location of the poles quite similarly.
On the other hand, the current measurement filter is normally included to avoid the current ripple due to the switching behavior of the converter being transmitted through the control loops and deteriorating the current quality amplifying the ripple, or even in extreme cases, making the system become unstable.

Study of the Voltage Measurement Filter Effect (Current Filter and Delay Included)
As done with the current filter, this subsection examines how the poles affect the voltage filter, the stability of the entire system, and its dynamic. The studied control block diagram is exactly the one presented previously in Figure 2, while the filter choice for the voltage measurement is (ωLP = 9425 rd/s): Therefore, the pole-zero map of the entire closed loops system yields is presented in Figure 8. The voltage filter effect, as occurred in the previously studied two cases, displaces in general all the poles of the system, in which some of them become more instable or closer to the unit circle. It can  Table 1, and applied to the numerical scenario of Appendix A. The voltage filter is not included, i.e., F v (z) = 1 and φ i = 0. Blue: without current filter, Red: with current filter and its phase shift compensation at 50 Hz.

Study of the Voltage Measurement Filter Effect (Current Filter and Delay Included)
As done with the current filter, this subsection examines how the poles affect the voltage filter, the stability of the entire system, and its dynamic. The studied control block diagram is exactly the one presented previously in Figure 2, while the filter choice for the voltage measurement is (ω LP = 9425 rd/s): Therefore, the pole-zero map of the entire closed loops system yields is presented in Figure 8. The voltage filter effect, as occurred in the previously studied two cases, displaces in general all the poles of the system, in which some of them become more instable or closer to the unit circle. It can also be remarked that although it is not shown for simplicity in the exposition, if the phase shift angle of the filter is not compensated, the pole locations becomes slightly worse from the stability point of view. It has to be highlighted that compared with the previously studied two cases, the inclusion of this filter move more poles closer to the unit circle, compromising more severely the stability of the system. As done before, many other filters have been tested (first order and different cut-off frequencies), reaching very similar results.
Energies 2020, 13, x FOR PEER REVIEW 10 of 42 also be remarked that although it is not shown for simplicity in the exposition, if the phase shift angle of the filter is not compensated, the pole locations becomes slightly worse from the stability point of view. It has to be highlighted that compared with the previously studied two cases, the inclusion of this filter move more poles closer to the unit circle, compromising more severely the stability of the system. As done before, many other filters have been tested (first order and different cut-off frequencies), reaching very similar results.  Figure 2, which is deduced from the mathematical model of Table 1 and applied to the numerical scenario of Appendix A. Blue: without voltage filter, Red: with voltage filter and its phase shift compensation at 50 Hz. The poles for the rest of the controlled input-output relations are equal (iα(z)/iq*(z), iβ(z)/id*(z), iβ(z)/iq*(z)).
The voltage measurement filter is normally included to avoid the voltage ripple at Rc and Cc passing through the feed-forward terms and deteriorating the control performance, or even in extreme cases, making the system become unstable.

Ripple Effect on the Control Loops
It is important to analyze how the ripple present at the measured voltage and currents affects  Figure 2, which is deduced from the mathematical model of Table 1 and applied to the numerical scenario of Appendix A. Blue: without voltage filter, Red: with voltage filter and its phase shift compensation at 50 Hz. The poles for the rest of the controlled input-output relations are equal (i α (z)/i q *(z), i β (z)/i d *(z), i β (z)/i q *(z)). The voltage measurement filter is normally included to avoid the voltage ripple at R c and C c passing through the feed-forward terms and deteriorating the control performance, or even in extreme cases, making the system become unstable.

Ripple Effect on the Control Loops
It is important to analyze how the ripple present at the measured voltage and currents affects the control behavior. It is quite common to use current and voltage filters (F i (s), F v (s)) to avoid the problematic interactions that can cause these ripples and the control performance. Based on the mathematical model proposed in this article, it is possible to know and precisely quantify how this ripple affects the control performance. Thus, the current and voltage ripples affect the output voltage of the converter, according to the following expressions (extracted from Table 1): Therefore, on the one hand, the current ripple impacts the output converter voltage through terms C αα , C αβ , C βα , and C ββ . Meanwhile, the voltage ripple that comes directly from V RCαβ is directly converted into a converter voltage ripple and therefore into an output current ripple and V RCαβ voltage ripple again. This fact is numerically evaluated by the bode diagrams depicted in Figure 9 (without voltage and current filters). As it can be noticed, with the numerical values used in this example, C αα and C ββ present an attenuation of −15.2 dB to the frequencies where the ripples appear (multiples of f sw ). While as noticed also in Figure 9, the converter voltage that can be polluted by this ripple is again attenuated by the output LCL filter with −16.4 dB for the current and with −22 dB for the voltage (at f sw ). These consecutive attenuations are sufficient to avoid the necessity of an extra digital current and voltage filters, at least for this specific numerical example. The nature of the system itself and the selected gains k p and T n , autonomously dampen the effect of the ripple due to the switching behavior of the converter. Nevertheless, in a different scenario with different numeric values (bigger current ripples, lower switching frequencies, different L, L g , R c, C c values . . . ), the proposed model should be employed again in order to evaluate whether the current and voltage filters are necessary or not. Note that as shown in previous subsections, in general, the presence of the current filter and specially the voltage filter deteriorates the robustness of the control, making it more unstable (poles closer to the unit circle).

Analysis of the Poles Location with Update Delay of the Converter and Without Filters at Measurement Points
Finally, this subsection shows the pole location of the system for this specific numerical example, assuming that there is delay in the converter voltage update and the voltage and current filters are not included, i.e., Fi(z) = 1, Fv(z) = 1.
This fact is shown in Figure 10, as well as the sensitivity of the poles to the different parameters. Note that all the transient defining five pole pairs are strongly dependent on the L value, as well as on  Finally, this subsection shows the pole location of the system for this specific numerical example, assuming that there is delay in the converter voltage update and the voltage and current filters are not included, i.e., This fact is shown in Figure 10, as well as the sensitivity of the poles to the different parameters. Note that all the transient defining five pole pairs are strongly dependent on the L value, as well as on the k p control parameter and sample time T s as well. It must be emphasized that it is important to perform the tuning of the control loops in the 'z' domain rather than in the 's' domain, since the poles location strongly depends on the T s sample time and also influences quite significantly the discretization method (zero-order hold, first-order hold, or others). So, conclusions obtained in the 's' domain could significantly change if the sample time T s cannot be set sufficiently small.

Analysis of the Poles Location with Update Delay of the Converter and Without Filters at Measurement Points
Finally, this subsection shows the pole location of the system for this specific numerical example, assuming that there is delay in the converter voltage update and the voltage and current filters are not included, i.e., Fi(z) = 1, Fv(z) = 1.
This fact is shown in Figure 10, as well as the sensitivity of the poles to the different parameters. Note that all the transient defining five pole pairs are strongly dependent on the L value, as well as on the kp control parameter and sample time Ts as well. It must be emphasized that it is important to perform the tuning of the control loops in the 'z' domain rather than in the 's' domain, since the poles location strongly depends on the Ts sample time and also influences quite significantly the discretization method (zero-order hold, first-order hold, or others). So, conclusions obtained in the 's' domain could significantly change if the sample time Ts cannot be set sufficiently small.
A summary of the characteristics of each pole for this specific numeric example is provided in Table 2. The dominant pole pair in this case presents a module of |z| = 0.9726 and is quite poorly damped. Its natural oscillating frequency corresponds to 866 rd/s, and its location depends on basically all of the control parameters together with the L value. The rest of the poles present less amplitude, which means they are less dominant or their presence is less important at the transient response. It is seen that they present different damping and natural frequencies.  Table 1, and applied to the numerical scenario of Appendix A. It is assumed that there is a delay in the converter voltage update, and the voltage and current filters are not included, i.e., F i (z) = 1, F v (z) = 1. The poles for the rest of the controlled input-output relations are equal (i α (z)/i q *(z), i β (z)/i d *(z), i β (z)/i q *(z)).
A summary of the characteristics of each pole for this specific numeric example is provided in Table 2. The dominant pole pair in this case presents a module of |z| = 0.9726 and is quite poorly damped. Its natural oscillating frequency corresponds to 866 rd/s, and its location depends on basically all of the control parameters together with the L value. The rest of the poles present less amplitude, which means they are less dominant or their presence is less important at the transient response. It is seen that they present different damping and natural frequencies.

Simulation-Based Validation
In this section, the proposed mathematical model is validated by comparing several time-domain responses of the mathematical model with Simulink model blocks. Thus, on the one hand, the control block diagram of Figure 2 is built with Matlab-Simulink blocks by using ideal voltage sources as an ideal converter model. In other words, the harmonics created by the power electronic converter due to its switching nature are not considered. Then, on the other hand, the mathematical model corresponding to this control (Table 1) is implemented in state-space form and is excited in parallel by the same input variables. Thus, at the simulations shown, steps are applied to the input references at two time instances. One is at zero time, just at the beginning of the experiment when the converter is connected to the grid voltage as well. Then, the second is at 0.1, once the first transient has been damped and a second step input is performed. As can be noticed in Figure 11, where the superposition of both time-domain responses of both models is presented, there is an exact match in the most representative Energies 2020, 13, 5849 14 of 47 variables of both models. As it is seen in the zooms, both transient and steady-state responses are exactly the same, which means that the mathematical expressions developed in this article are correct.

Simulation-Based Validation
In this section, the proposed mathematical model is validated by comparing several timedomain responses of the mathematical model with Simulink model blocks. Thus, on the one hand, the control block diagram of Figure 2 is built with Matlab-Simulink blocks by using ideal voltage sources as an ideal converter model. In other words, the harmonics created by the power electronic converter due to its switching nature are not considered. Then, on the other hand, the mathematical model corresponding to this control (Table 1) is implemented in state-space form and is excited in parallel by the same input variables. Thus, at the simulations shown, steps are applied to the input references at two time instances. One is at zero time, just at the beginning of the experiment when the converter is connected to the grid voltage as well. Then, the second is at 0.1, once the first transient has been damped and a second step input is performed. As can be noticed in Figure 11, where the superposition of both time-domain responses of both models is presented, there is an exact match in the most representative variables of both models. As it is seen in the zooms, both transient and steadystate responses are exactly the same, which means that the mathematical expressions developed in this article are correct.  Table  1 and the control block diagram of Figure 2, built in Matlab-Simulink blocks by using ideal voltage sources as an ideal converter model. The same inputs and the same numerical parameters are used in both models (Appendix A).  Table 1 and the control block diagram of Figure 2, built in Matlab-Simulink blocks by using ideal voltage sources as an ideal converter model. The same inputs and the same numerical parameters are used in both models (Appendix A).
Then, once the mathematical model has been validated by using ideal voltage sources, it is confronted with a Simulink blocks-based model that uses a 3L-NPC power electronic converter switched at f sw . In this case, the switching nature of the voltage source will produce voltage and current harmonics as occurs in real application scenarios. Thus, exactly the same experiment as before is prepared, and the most important electric variables are registered in Figure 12.
In this case, the currents shown from the Simulink blocks are the ones sampled at the T s sample time (seen by the control blocks). It can be noticed that the proposed mathematical model follows the average value (at T s ) of the signals from Simulink blocks. So, both models produce the same average electric variables, but the mathematical model does not produce the harmonics due to the switching nature of the power electronic converter. This fact demonstrates that the mathematical model proposed reproduces the average behavior (at T s ) of the Simulink model. In this way, the model can be useful to deduce the stability issues of the real system as is shown in this paper.
Then, once the mathematical model has been validated by using ideal voltage sources, it is confronted with a Simulink blocks-based model that uses a 3L-NPC power electronic converter switched at fsw. In this case, the switching nature of the voltage source will produce voltage and current harmonics as occurs in real application scenarios. Thus, exactly the same experiment as before is prepared, and the most important electric variables are registered in Figure 12.  Table  1 and the control block diagram of Figure 2 (variables sampled at Ts) built in Matlab-Simulink blocks by using a 3L-NPC (neutral point clamped) power electronic converter switched at fsw. The same input variables are applied at both controls, and the same numerical parameters are used in both models as well (Appendix A).
In this case, the currents shown from the Simulink blocks are the ones sampled at the Ts sample time (seen by the control blocks). It can be noticed that the proposed mathematical model follows the average value (at Ts) of the signals from Simulink blocks. So, both models produce the same average electric variables, but the mathematical model does not produce the harmonics due to the switching nature of the power electronic converter. This fact demonstrates that the mathematical model proposed reproduces the average behavior (at Ts) of the Simulink model. In this way, the model can be useful to deduce the stability issues of the real system as is shown in this paper.

Experimental Validation of the Basic Equations of the Model
In this subsection, the experimental validation of the model equations is performed. For that purpose, the same validation idea followed in the simulations of the previous subsections is carried out. Figure 13 shows the experimental platform that has been used for the validation together with an example of the experimental results obtained. It is seen that the superposition of the time-domain step responses of the mathematical model proposed in Table 1 and the correspondent experimental results match very well. Note that only a step has been injected uniquely at iq current, since in this experimental setup, the id current is used for the control of the DC bus voltage (there is not a DC source at the DC side), so it is not an input variable that can be modified freely [3]. It has to be remarked that several other experiments have been carried out at different operating points and with different control parameters, obtaining always a good agreement between the experimental timedomain responses and the theoretical model. Hence, having validated the basic mathematical expressions developed and proposed in the paper, which are also used in the next controls studied in subsequent sections, in the following analyses and validations, only simulation-based validations will be presented for simplicity.  Table 1 and the control block diagram of Figure 2 (variables sampled at T s ) built in Matlab-Simulink blocks by using a 3L-NPC (neutral point clamped) power electronic converter switched at fsw. The same input variables are applied at both controls, and the same numerical parameters are used in both models as well (Appendix A).

Experimental Validation of the Basic Equations of the Model
In this subsection, the experimental validation of the model equations is performed. For that purpose, the same validation idea followed in the simulations of the previous subsections is carried out. Figure 13 shows the experimental platform that has been used for the validation together with an example of the experimental results obtained. It is seen that the superposition of the time-domain step responses of the mathematical model proposed in Table 1 and the correspondent experimental results match very well. Note that only a step has been injected uniquely at i q current, since in this experimental setup, the i d current is used for the control of the DC bus voltage (there is not a DC source at the DC side), so it is not an input variable that can be modified freely [3]. It has to be remarked that several other experiments have been carried out at different operating points and with different control parameters, obtaining always a good agreement between the experimental time-domain responses and the theoretical model. Hence, having validated the basic mathematical expressions developed and proposed in the paper, which are also used in the next controls studied in subsequent sections, in the following analyses and validations, only simulation-based validations will be presented for simplicity.

Summary of Conclusions of the Section
The most important conclusions of the analysis carried out in this section can be summarized as follows: there are two dominant pole pairs that determine the dynamic performance of the control. These dominant poles depend mainly on the following parameters of the system; k p , T n , L, L g , and T s . Thus, the tuning objective can be focused on choosing the appropriate k p and T n values to locate the dominant poles in a desired region. The other "non-dominant" poles of the system present different dependence on the parameters of the system. On the other hand, in the studied numerical scenario, it is seen that the converter's delay effect does not significantly affect the location of the dominant poles and therefore, the stability of the system. However, in a numerical scenario where for instance the switching frequency is much smaller and therefore the delay is bigger, the effect of this delay should be carefully quantified by means of the model proposed. Finally, again in the studied numerical scenario, it has been deduced that the current and voltage filters are better if they are not used, since they affect the stability of the system quite significantly. In addition, they are not required to attenuate the current and voltage harmonics that come through the measurements, since the control system itself can attenuate them naturally. To conclude, by comparison of time-domain responses, the proposed mathematical model has been validated with simulation and experimental results.

Summary of Conclusions of the Section
The most important conclusions of the analysis carried out in this section can be summarized as follows: there are two dominant pole pairs that determine the dynamic performance of the control. These dominant poles depend mainly on the following parameters of the system; kp, Tn, L, Lg, and Ts. Thus, the tuning objective can be focused on choosing the appropriate kp and Tn values to locate the dominant poles in a desired region. The other "non-dominant" poles of the system present different dependence on the parameters of the system. On the other hand, in the studied numerical scenario, it is seen that the converter's delay effect does not significantly affect the location of the dominant poles and therefore, the stability of the system. However, in a numerical scenario where for instance the switching frequency is much smaller and therefore the delay is bigger, the effect of this delay should be carefully quantified by means of the model proposed. Finally, again in the studied numerical scenario, it has been deduced that the current and voltage filters are better if they are not used, since they affect the stability of the system quite significantly. In addition, they are not required to attenuate the current and voltage harmonics that come through the measurements, since the control system itself can attenuate them naturally. To conclude, by comparison of time-domain responses, the proposed mathematical model has been validated with simulation and experimental results.

Control Scheme
In this section, a quite commonly used control strategy is modeled and studied. Compared with the control strategy studied in the previous section, this control is prepared to control the positive (idp*, iqp*) and negative (idn*, iqn*) sequence currents in a rotating reference frame. This type of control is very useful when for instance the converter is operating in an unbalanced voltage grid. It is useful as well to be able to provide the required positive and negative sequence currents as demanded by the application where this converter is being used. The schematic control block diagram is depicted in Figure 14 [1][2][3]. As it can be noticed, in essence, it is basically a "duplication" of the previously seen control ( Figure 2).

Control Scheme
In this section, a quite commonly used control strategy is modeled and studied. Compared with the control strategy studied in the previous section, this control is prepared to control the positive (i dp *, i qp *) and negative (i dn *, i qn *) sequence currents in a rotating reference frame. This type of control is very useful when for instance the converter is operating in an unbalanced voltage grid. It is useful as well to be able to provide the required positive and negative sequence currents as demanded by the application where this converter is being used. The schematic control block diagram is depicted in Figure 14 [1][2][3]. As it can be noticed, in essence, it is basically a "duplication" of the previously seen control ( Figure 2).
In this control strategy, some kind of filtering action is needed after transforming the currents at positive or negative sequences. Since at each sequence, the other sequence is coupled in an oscillation of 100 Hz, very commonly, a notch filter tuned at 100 Hz is used, as shown in the control block diagram of Figure 14 [1][2][3].
It is possible to find authors that use several other types of filters, such as Digital Signal Cancellation (DSC) methods or others [22] that are reasonably equivalent.

Mathematical Model of the Control
Hence, by following the same modeling methodology applied to the previous control, the model equations in the αβ reference frame of the block diagram of Figure 14 are summarized in Table 3. In this control strategy, some kind of filtering action is needed after transforming the currents a positive or negative sequences. Since at each sequence, the other sequence is coupled in an oscillation of 100 Hz, very commonly, a notch filter tuned at 100 Hz is used, as shown in the control block diagram of Figure 14 It is possible to find authors that use several other types of filters, such as Digital Signa Cancellation (DSC) methods or others [22] that are reasonably equivalent.

Mathematical Model of the Control
Hence, by following the same modeling methodology applied to the previous control, the mode equations in the αβ reference frame of the block diagram of Figure 14 are summarized in Table 3. Table 3. Mathematical model of the studied control. Control equations for loops dedicated to harmonics h = 1 (50 Hz) that can be implemented sequentially in a Matlab program for instance. c2d: continuous to discrete. The method used in this paper is ZOH.
Positive Sequence Negative Sequence Figure 14. Fundamental current control loops for positive and negative sequences using PI regulators in the dq reference frame: dual vector control.
Compared to the previous control strategy, the inclusion of the two new loops dq at the negative sequence, together with the usage of the notch filters at each current loop, provoke that the resulting global mathematical model is more complex. It is more complex in the sense that is composed by larger equations and of higher orders. However, in essence, the structure of the equations is very similar to the control strategy studied in the previous sections.

Summary of the Analysis Carried out with the Assistance of the Mathematical Model
In the subsequent sections, the following analyses are carried out. First of all, it is evaluated if the filter at the current measurement is necessary to avoid harmonics influence on the control. Then, which criterion should be used to choose the damping factor of the required notch filter is studied. After that, the poles that dominate the dynamic performance of the control are found and the most important parameters that determine the location of each pole are identified. Afterwards, a method for tuning the gains k p , T n , ξ is proposed, trying to optimize multiple objectives, maximize the stability and dynamic performance, reduce harmonic impact on the control, and other factors. To conclude, two illustrative analysis examples are provided, showing first that the algorithm effectively searches appropriate gain values in a standard numerical scenario. Secondly, it is studied whether it is possible to reduce the power losses of the system by reducing the damping R resistance of the filter without compromising the performance and stability of the system. Table 3. Mathematical model of the studied control. Control equations for loops dedicated to harmonics h = 1 (50 Hz) that can be implemented sequentially in a Matlab program for instance. c2d: continuous to discrete. The method used in this paper is ZOH.
Connection of all the elements including the grid, filter and converter models ( Table 1): Derivation of the 'Global model matrix', through the 'Connect' function from Matlab (the choice of the output variables is not unique):

Ripple Effect on the Control Loops
As in the previous section with the previous control, it is important to analyze how the ripple present at the measured currents affects the control behavior. In this section, the voltage ripple analysis is not repeated again, since it is already studied in the previous section and is only affected by the feed-forward terms, or in other words, it is not influenced by the feedback loops. Therefore, only how the current ripple affects the control has to be evaluated. The current ripple impacts the output converter voltage through terms C ααp , C αβp , C βαp , C ββp and C ααn , C αβn , C βαn , C ββn . This effect is numerically evaluated by the Bode diagrams depicted in Figure 15. As it can be noticed, with the numerical values used in this example, the attenuation at the frequencies where the ripples appear (multiples of f sw ) is −9.11 dB, for terms C ααp , C ββp , C ααn , C ββn while for the rest of the terms, the attenuation is even bigger. As in previous control analysis, these attenuations together with the voltage attenuations due to the LCL filter are sufficient to avoid the necessity of an extra digital current and voltage filter, at least for this specific numerical example. The nature of the system itself and the selected gains k p and T n autonomously dampen the effect of the ripple produced by the converter. Nevertheless, in a different scenario with different numeric values (bigger current ripples, lower switching frequencies, different L, L g , R c, C c values . . . ), the model should be employed again in order to evaluate whether the current and voltage filters are necessary or not.
As in the previous section with the previous control, it is important to analyze how the r nt at the measured currents affects the control behavior. In this section, the voltage r sis is not repeated again, since it is already studied in the previous section and is only aff e feed-forward terms, or in other words, it is not influenced by the feedback loops. Ther how the current ripple affects the control has to be evaluated. The current ripple impac ut converter voltage through terms Cααp, Cαβp, Cβαp, Cββp and Cααn, Cαβn, Cβαn, Cββn. This eff erically evaluated by the Bode diagrams depicted in Figure 15. As it can be noticed, wit erical values used in this example, the attenuation at the frequencies where the ripples ap tiples of fsw) is −9.11 dB, for terms Cααp, Cββp, Cααn, Cββn while for the rest of the terms uation is even bigger. As in previous control analysis, these attenuations together wit ge attenuations due to the LCL filter are sufficient to avoid the necessity of an extra d nt and voltage filter, at least for this specific numerical example. The nature of the system the selected gains kp and Tn autonomously dampen the effect of the ripple produced b erter. Nevertheless, in a different scenario with different numeric values (bigger current rip r switching frequencies, different L, Lg, Rc, Cc values…), the model should be employed ag r to evaluate whether the current and voltage filters are necessary or not.

Notch Filter Choice
In this subsection, some useful guidelines for choosing the damping value of the notch filter are provided. The notch filter is used mainly to eliminate the oscillations of 100 Hz produced by the sequences decompositions of the currents, and at the same time, it is desirable to have a minimal impact on the dynamics of the current loops. As graphically represented in Figure 16a, the sequence decomposition produces the superposition of two currents: a dc current plus a sinusoidal current at 100 Hz that is desired to be eliminated [3]. Thus, at steady state, the filter must try to eliminate the 100 Hz oscillations, but during transient of currents as well, the notch filter should produce as fast a transient as possible (fast damping) at both sinusoidal and step inputs. In addition, the overshoot (or peak) value during transients should be also as small as possible at both sinusoidal and step inputs, so the PI regulator that is coming after reacts as little as possible. All these requirements are evaluated for the studied case in Figure 16. It is seen that some of the requirements demand a big damping value, while some others demand a low damping value. as little as possible. All these requirements are evaluated for the studied case in Figure 16. It is seen that some of the requirements demand a big damping value, while some others demand a low damping value.
The information deduced from Figure 16 is summarized in Table 4. It can be inferred that a multi-objective choice of the damping ratio should be carried out, trying to find a balanced minimization of all the objectives. Then, note also that how the damping affects the overall control performance should be analyzed as well, but this is carried out in subsequent subsections.
Sinusoidal input response Bode Figure 16. (a) Graphical representation that shows how sequence decomposition produces the superposition of two currents; a dc current plus a sinusoidal current at 100 Hz, (b) step response of a notch filter tuned at 100 Hz with different damping ratios, (c) response to a sinusoidal input of a notch filter tuned at 100 Hz with different damping ratios, (d) Bode diagrams of a notch filter tuned at 100 Hz with different damping ratios.
The information deduced from Figure 16 is summarized in Table 4. It can be inferred that a multi-objective choice of the damping ratio should be carried out, trying to find a balanced minimization of all the objectives. Then, note also that how the damping affects the overall control performance should be analyzed as well, but this is carried out in subsequent subsections.

Analysis of the Poles Location
As seen in previous subsections, at least for this specific numeric scenario, it is possible to avoid the inclusion of voltage and current digital filters. In addition, having seen the criteria to choose the damping ratio of the notch filters, in this subsection, the pole location of the global systems is analyzed. Hence, the poles and zeros of the global closed loop system in the αβ reference frame are depicted in Figure 17. There are six pole pairs, of which five pole pairs are defining the transient response and one pole pair is responsible for the steady-state αβ currents. Figure 17 shows which parameters have a more important effect on the location of each pole pair. As done before, these parameters' dependency has been deduced with the parameter variations shown in Figure 18. In comparison with the control analyzed in the previous section, it can be said that the poles of the dual control studied here are closer to the unit circle in general. This means that the poles are closer to the instability, or if preferred, they are damped dynamically slower.
damping ratio of the notch filters, in this subsection, the pole location of the global systems is analyzed. Hence, the poles and zeros of the global closed loop system in the αβ reference frame are depicted in Figure 17. There are six pole pairs, of which five pole pairs are defining the transient response and one pole pair is responsible for the steady-state αβ currents. Figure 17 shows which parameters have a more important effect on the location of each pole pair. As done before, these parameters' dependency has been deduced with the parameter variations shown in Figure 18. In comparison with the control analyzed in the previous section, it can be said that the poles of the dual control studied here are closer to the unit circle in general. This means that the poles are closer to the instability, or if preferred, they are damped dynamically slower. Figure 17. Poles and zeros of the currents for the global closed loop system in the αβ reference frame for step input references, corresponding to the control of Figure 14, deduced from the mathematical model of Table 3, and applied to the numerical scenario of Appendix B. It is assumed that there is delay in the converter voltage update and the voltage and current filters are not included, i.e., Fi(z) = 1, Fv(z) = 1.
A summary of the characteristics of each pole for this specific numeric example is provided in Table 5. Note that the new ξ parameter introduced for the notch filter affects the location of practically all the poles. The dominant pole pair in this case presents a module of |z| = 0.99526 and is very poorly damped. Its natural oscillating frequency corresponds to 898 rd/s, and its location depends on basically all the control parameters together with the L value.    Figure 14, deduced from the mathematical model of Table 3, and applied to the numerical scenario of Appendix B. It is assumed that there is delay in the converter voltage update and the voltage and current filters are not included, i.e., F i (z) = 1, F v (z) = 1.
A summary of the characteristics of each pole for this specific numeric example is provided in Table 5. Note that the new ξ parameter introduced for the notch filter affects the location of practically all the poles. The dominant pole pair in this case presents a module of |z| = 0.99526 and is very poorly damped. Its natural oscillating frequency corresponds to 898 rd/s, and its location depends on basically all the control parameters together with the L value.

Multi-Objective Tuning
This section tries to give one of the various possible solutions to obtain a multi-objective tuning of the control loops. As it has been seen in previous subsections, a proper tuning of the loops should be able to achieve a good balance between the multiple objectives that must be fulfilled. As often occurs in a complex coupled system, the improvement of a few objectives implies the deterioration of some others, so a compromise needs to be found. Figure 19 graphically illustrates how the tuning algorithm could be conceived.
The inputs of the algorithm are the electric parameters of the system, i.e., the parameters of the LCL filter, the switching frequency, the electric parameters of the converter, and so on. Then, based on the mathematical expressions of the model summarized in Table 6, in this case, six different parameters or indicators are evaluated. These parameters or indicators are the ones that we want to optimize, or in other words, we can also say that these are the objectives we want to fulfill. In Figure 19, it has been represented as a Spider Chart, where at each axis (or objective), we want to obtain one specific numeric value. Then, the output values provided by the algorithm are the parameter values that fulfill the defined objectives. An explanation of the objective, the corresponding mathematical expression to be evaluated, an indicative numeric value for each objective, and an example of the code program in Matlab for implementation are provided in Table 6.
Energies 2020, 13, x FOR PEER REVIEW 20 of 42 Figure 18. Variation of the poles and zeros of the output currents of the global closed loop system for step input references, corresponding to the control of Figure 14, deduced from the mathematical model of Table 3, and applied to the numerical scenario of Appendix A. Each parameter is modified as specified. It is assumed that there is delay in the converter voltage update and the voltage and   Figure 14, deduced from the mathematical model of Table 3, and applied to the numerical scenario of Appendix A. Each parameter is modified as specified. It is assumed that there is delay in the converter voltage update and the voltage and current

Multi-Objective Tuning
This section tries to give one of the various possible solutions to obtain a multi-objective tuning of the control loops. As it has been seen in previous subsections, a proper tuning of the loops should be able to achieve a good balance between the multiple objectives that must be fulfilled. As often occurs in a complex coupled system, the improvement of a few objectives implies the deterioration of some others, so a compromise needs to be found. Figure 19 graphically illustrates how the tuning algorithm could be conceived. The inputs of the algorithm are the electric parameters of the system, i.e., the parameters of the LCL filter, the switching frequency, the electric parameters of the converter, and so on. Then, based on the mathematical expressions of the model summarized in Table 6, in this case, six different parameters or indicators are evaluated. These parameters or indicators are the ones that we want to optimize, or in other words, we can also say that these are the objectives we want to fulfill. In Figure  19, it has been represented as a Spider Chart, where at each axis (or objective), we want to obtain one specific numeric value. Then, the output values provided by the algorithm are the parameter values that fulfill the defined objectives. An explanation of the objective, the corresponding mathematical expression to be evaluated, an indicative numeric value for each objective, and an example of the code program in Matlab for implementation are provided in Table 6.
It has to be emphasized that probably the most determinant objective is the one that tries to find a module minimization of the pole with a bigger module (dominant pole) in the "z" domain. A graphical representation of this objective is depicted in Figure 20. Thus, in the "z" domain, the dominant poles are the ones that present bigger a module, not having a zero nearby that tends to cancel its effect. Consequently, the poles with a bigger module in the "z" domain would also present the biggest residue at the time domain response and also the biggest time constant. This means that they will be the poles whose effect needs more time to be extinguished during transients as well as the ones that are closer to the instability limit circle or have a radius equal to 1. Therefore, by setting the poles within a desired radio of a circle (less than 1), the settling time of the transient response is defined as well as the margin to reach the instability as well. Hence, the remaining five objectives are also useful to be supervised and set within certain values due to the reasons provided in previous subsections and because it is information that is not contained in the poles. Nevertheless, it is worth mentioning that these five objectives are not supposed to be minimized. We set a maximum limit that should not be exceeded in order to avoid problems in the performance of the control. However, their achieved final value could be anything below the set limit. The way in which the algorithm is implemented and finds the solution is not deeply explained here. Once the mathematical expressions and the objectives or minimizations  It has to be emphasized that probably the most determinant objective is the one that tries to find a module minimization of the pole with a bigger module (dominant pole) in the "z" domain. A graphical representation of this objective is depicted in Figure 20. Thus, in the "z" domain, the dominant poles are the ones that present bigger a module, not having a zero nearby that tends to cancel its effect. Consequently, the poles with a bigger module in the "z" domain would also present the biggest residue at the time domain response and also the biggest time constant. This means that they will be the poles whose effect needs more time to be extinguished during transients as well as the ones that are closer to the instability limit circle or have a radius equal to 1. Therefore, by setting the poles within a desired radio of a circle (less than 1), the settling time of the transient response is defined as well as the margin to reach the instability as well. Hence, the remaining five objectives are also useful to be supervised and set within certain values due to the reasons provided in previous subsections and because it is information that is not contained in the poles. Nevertheless, it is worth mentioning that these five objectives are not supposed to be minimized. We set a maximum limit that should not be exceeded in order to avoid problems in the performance of the control. However, their achieved final value could be anything below the set limit. The way in which the algorithm is implemented and finds the solution is not deeply explained here. Once the mathematical expressions and the objectives or minimizations to be achieved are known (Table 6), several search possibilities exist that can be applied. In this paper, simply three concatenated "for loops" have been used in an '.m' Matlab program, evaluating all the possible parameters to be analyzed (k p , T n and ξ) within the desired ranges and checking the actual k p , T n , and ξ values that optimally fulfill the specified objectives. Minimize the module of the pole with a bigger module, excluding the steady-state poles (in the αβ reference frame) Module of poles of: [G mm (z)] Find the minimum [poles,zeros] = pzmap(Gmm); Mod_pol_max = max(abs(poles)); Set a maximum limit to the settling time of notch filters to the step input (in the dq reference frame) s 2 +ω2n 2 s 2 +2ζω2ns+ω2n 2 <0.08 Set a maximum limit to the overshoot of the notch filters to the step input (in the dq reference frame) Set a maximum limit to the effect of the ripple at fsw in the control loops (in the αβ reference frame) C ααp , C αβp , C βαp , C ββp and C ααn , C αβn , C βαn , C ββn <0.4 [Magsw1,phase,wout]= bode(Cal_al_p, wsw); [Magsw2,phase,wout]= bode(Cbe_al_p, wsw); Set a maximum limit to the 100 Hz oscillations at output of the notch (in the dq reference frame) Finally, it has to be remarked that the reader can enrich this algorithm by increasing the number of objectives or increasing the restrictions not only of the dominant poles but also of some of the other poles in order to minimize the oscillatory behaviors, and many other options are available as well.  Table 7, with the obtained control parameters, we achieve a value for the pole with the bigger module of |z| = 0.99086. Compared with the initial module with the initial control parameters, which was 0.9952, it is noted that a considerable improvement in the dynamics and stability margin has been gained. Then, as it can be seen in Table 7, the rest of the objectives have been also achieved to be below the specification. Finally, it has to be remarked that the reader can enrich this algorithm by increasing the number of objectives or increasing the restrictions not only of the dominant poles but also of some of the other poles in order to minimize the oscillatory behaviors, and many other options are available as well.

Example 1: Search the Minimum Module Pole
This subsection shows as an example how the proposed multi-objective tuning method can search a set of control parameters according to the specifications. Thus, the initial set of parameters, are shown in Appendix B, are k p = 0.35, T n = 0.00457, and ξ = 0.08. By applying the proposed algorithm in the previous subsection and evaluating all the possible numerical options for k p , T n , and ξ (with a reasonable low precision for each parameter), the obtained values fulfilling the specifications are k p = 0.24, T n = 0.0065, and ξ = 0.096. As summarized in Table 7, with the obtained control parameters, we achieve a value for the pole with the bigger module of |z| = 0.99086. Compared with the initial module with the initial control parameters, which was 0.9952, it is noted that a considerable improvement in the dynamics and stability margin has been gained. Then, as it can be seen in Table 7, the rest of the objectives have been also achieved to be below the specification. Table 7. Obtained values after the multi-objective search: k p = 0.24, T n = 0.0065, and ξ = 0.096.

Objective Value for Specification Obtained Results
Minimize the module of the pole with the bigger module, excluding the steady-state poles (in the αβ reference frame) Hence, with the achieved optimum control parameters according to the defined objectives, the pole-zero map obtained is represented in Figure 21.
To conclude, the superposition of the time-domain step response with Matlab-Simulink blocks and the proposed mathematical model depicted in Figure 22a,b reveals an exact match of both. the output of the notch (in the dq reference frame) Set a maximum limit to any 100 Hz oscillations that transitorily can excite the PI (in the dq reference frame) <0.25 0.2483 (−12.1 dB) Hence, with the achieved optimum control parameters according to the defined objectives, the pole-zero map obtained is represented in Figure 21. In this subsection, one more example of analysis is provided. In this case, seeking to minimize the damping resistance of the LCL filter, the algorithm is applied with an R c = 250 × 10 −3 /50 = 50 × 10 −3 Ω, which is a value considerably lower than the previous cases. By applying the proposed searching algorithm and evaluating all the possible numerical options for k p , T n , and ξ (within a reasonable low precision for each parameter), the obtained values fulfilling the specifications are k p = 0.23, T n = 0.00487, and ξ = 0.12. As is summarized in Table 8, with the obtained control parameters, we achieve a value for the pole with a bigger module of |z| = 0.9878. Again, as occurred in the previous subsection, it is noted that a considerable improvement in the dynamics and stability margin has been obtained. Then, as can be seen in Table 8, the rest of the objectives have been also achieved to be below the specification.  Hence, with the achieved optimum control parameters according to the defined objectives, the pole-zero map obtained is represented in Figure 21.  Table 3 and the control block diagram of Figure 14 (variables sampled at Ts) built in Matlab-Simulink blocks. Same input variables are applied at both controls, and the same numerical parameters are used in both models as well (Appendix B). (a) Ideal voltage sources are used as an ideal power electronic converter, (b) A 3L-NPC power electronic converter switched at fsw is used as an power electronic converter.

Example 2: Search the Minimum Module Pole, Reducing the Damping Resistance of the LCL Filter to Rc = 50 m Ω
In this subsection, one more example of analysis is provided. In this case, seeking to minimize the damping resistance of the LCL filter, the algorithm is applied with an Rc = 250 × 10 −3 /50 = 50 × 10 −3 Ω, which is a value considerably lower than the previous cases. By applying the proposed searching algorithm and evaluating all the possible numerical options for kp, Tn, and ξ (within a reasonable low precision for each parameter), the obtained values fulfilling the specifications are kp = 0.23, Tn = 0.00487, and ξ = 0.12. As is summarized in Table 8, with the obtained control parameters, we achieve a value for the pole with a bigger module of |z| = 0.9878. Again, as occurred in the previous subsection, it is noted that a considerable improvement in the dynamics and stability margin has been obtained. Then, as can be seen in Table 8, the rest of the objectives have been also achieved to be below the specification. Table 8. Obtained values after the multi-objective search with Rc = 50 × 10 −3 Ω: kp = 0.23, Tn = 0.00487, and ξ = 0.12.

Objective Value for Specification Obtained Results
Minimize the module of the pole with bigger module, excluding the steady-state poles (in the αβ reference frame) Thus, with the achieved optimum control parameters, the pole-zero map obtained is represented in Figure 23.  Table 3 and the control block diagram of Figure 14 (variables sampled at T s ) built in Matlab-Simulink blocks. Same input variables are applied at both controls, and the same numerical parameters are used in both models as well (Appendix B). (a) Ideal voltage sources are used as an ideal power electronic converter, (b) A 3L-NPC power electronic converter switched at f sw is used as an power electronic converter. Table 8. Obtained values after the multi-objective search with R c = 50 × 10 −3 Ω: k p = 0.23, T n = 0.00487, and ξ = 0.12.

Objective Value for Specification Obtained Results
Minimize the module of the pole with bigger module, excluding the steady-state poles (in the αβ reference frame) Thus, with the achieved optimum control parameters, the pole-zero map obtained is represented in Figure 23.

Summary of Conclusions of the Section
The most important conclusions of the analysis carried out in this section can be summarized as follows: as evaluated in the previous control, it is concluded that the current filter is not necessary to avoid a current harmonics effect on the control loops. In this numeric scenario, the control itself naturally attenuates the harmonics of the system. Then, it is seen that there are three pole pairs that can be considered as dominant, since they present similar modules. These poles are located mainly by the

Summary of Conclusions of the Section
The most important conclusions of the analysis carried out in this section can be summarized as follows: as evaluated in the previous control, it is concluded that the current filter is not necessary to avoid a current harmonics effect on the control loops. In this numeric scenario, the control itself naturally attenuates the harmonics of the system. Then, it is seen that there are three pole pairs that can be considered as dominant, since they present similar modules. These poles are located mainly by the numerical values of k p , T n , ξ, L, L g , and T s . After that, it is shown that the proposed multi-objective tuning algorithm effectively finds the appropriate k p , T n , and ξ values, optimizing the desired objectives. Finally, it is also demonstrated that the damping R resistance of the filter can be significantly reduced, increasing the efficiency of the system, by choosing the appropriate k p , T n , and ξ values found by the proposed algorithm.

Control Scheme
This section models the last control strategy studied in the paper. It is oriented to control not only currents at the grid's main 50 Hz frequency but also to control current harmonics as well. It is often required in applications such as active filters, active impedances [23,24], etc., where some harmonic currents must be also controlled in order to mitigate associated problems due to their presence. The control block diagram is depicted in Figure 24 [4]. As occurred in the control studied in the previous section, in essence again, it can be understood as a replication of the previously seen control. It incorporates a duplication of the four current control loops dedicated at the 50 Hz frequency for another h 1 harmonic. Thus, four current control loops are used to control the 50 Hz currents and another four new loops are used to control the currents at harmonic h 1. Depending on the application scenario, it could be necessary to control more harmonics; h 1 , h 2, h 3 ... However, for simplicity in the exposition, only one harmonic is considered in this paper. Nevertheless, the analysis method presented and conclusions raised can be extended to any number of harmonics.

Mathematical Model of the Control
Hence, by following the same modeling methodology applied to the previously seen controls, the mathematical model equations in the αβ reference frame of the block diagram of Figure 24 are summarized in Tables 9 and 10. Compared to the previous control strategy, it is necessary to include a number of notch filters in series at each loop. This fact is studied in detail in the next subsection. Again, in order not to be repetitive, only the control model equations are covered in Tables 9 and 10. The filter and grid model equations are the same as presented in Table 1. The modeling process is equivalent to the one followed in the previous controls. Control terms (C ααp (s), C αβp (s), and C βαp (s), etc.) must be derived for each harmonic current control; then, we discretize all of them, and add all their effects in order to obtain the total output voltage references v controlαβ (z). Finally, the control model equations are combined with the grid and filter model equations (Table 1) and the global model matrix equation G mm (z) is achieved that relates the output electric magnitudes i αβ (z) and v RCαβ (z) as a function of the inputs of the system: i dqp *(z), i dqn *(z), i dqp1 *(z), i dqn1 *(z), and v gαβ (z).

Summary of the Analysis Carried out with the Assistance of the Mathematical Model
In the subsequent sections, the following analyses are carried out. First of all, it is evaluated how the different harmonics are coupled through the loops in order to know which notch filters are necessary at each control loop. Then, how the current ripple due to the switching affects the control loops is evaluated. It is detected that this control is more sensible to the switching harmonics than the previously studied controls. After that, the performance of the different combinations of the notch filters included at every loop is analyzed. Then, the poles location is studied, and the most important parameters that determine this location are identified. Afterwards, the multi-objective tuning method of the gains is proposed. Finally, several illustrative examples in different numerical conditions are provided, emphasizing a robust evaluation of the control.

Mathematical Model of the Control
Hence, by following the same modeling methodology applied to the previously seen controls, the mathematical model equations in the αβ reference frame of the block diagram of Figure 24 are summarized in Tables 9 and 10. Compared to the previous control strategy, it is necessary to include a number of notch filters in series at each loop. This fact is studied in detail in the next subsection. Again, in order not to be repetitive, only the control model equations are covered in Tables 9 and 10. The filter and grid model equations are the same as presented in Table 1. The modeling process is equivalent to the one followed in the previous controls. Control terms (Cααp(s), Cαβp(s), and Cβαp(s), etc.) must be derived for each harmonic current control; then, we discretize all of them, and add all their effects in order to obtain the total output voltage references vcontrolαβ(z). Finally, the control model equations are combined with the grid and filter model equations (Table 1) and the global model matrix equation Gmm(z) is achieved that relates the output electric magnitudes iαβ(z) and vRCαβ(z) as a function of the inputs of the system: idqp*(z), idqn*(z), idqp1*(z), idqn1*(z), and vgαβ(z).

Summary of the Analysis Carried out with the Assistance of the Mathematical Model
In the subsequent sections, the following analyses are carried out. First of all, it is evaluated how the different harmonics are coupled through the loops in order to know which notch filters are necessary at each control loop. Then, how the current ripple due to the switching affects the control loops is evaluated. It is detected that this control is more sensible to the switching harmonics than the previously studied controls. After that, the performance of the different combinations of the notch filters included at every loop is analyzed. Then, the poles location is studied, and the most important parameters that determine this location are identified. Afterwards, the multi-objective tuning method of the gains is proposed. Finally, several illustrative examples in different numerical conditions are provided, emphasizing a robust evaluation of the control.
Negative Sequence vregβn(t)+vregβp1(t)+vregβn1(t) Figure 24. Fundamental and harmonic current control loops for positive and negative sequences using PI regulators in the dq reference frame: active Filter.

Frequency of Currents Coupled among the Loops
The control strategy studied in this section incorporates control loops at different frequencies: two pairs of control loops for the h = 1 (50 Hz) frequency and two more pairs of control loops for the n 1 = 13 (650 Hz) frequency. This fact, together with the essential nature of the control loops, which presents positive and negative sequence decompositions, means that these frequencies are modified inside each loop. For this specific case, the harmonic frequencies that are present at each loop are specified in Table 11. Note that whenever there is a current reference value at each current control loop, a coupled current harmonic appears in another loop.
Thus, for instance, in the example covered in Table 11, at the positive sequence of harmonic h = 1 loops, its own current references appear as DC currents, but the other currents of the other loops are coupled as frequencies of h + h (100 Hz), h 1 -h (600 Hz), and h 1 + h (700 Hz). From this fact arises the necessity of incorporating different notch filters of different frequencies at each current loop. Table 9. Mathematical model of the studied control. Control equations for loops dedicated to harmonic h 1 = 1 (50 Hz) that can be implemented sequentially in a Matlab program (requires Table 10 to be completed).  Table 9 to be completed).
Energies 2020, 13, 5849 33 of 47 Connection of all the elements including the grid, filter and converter models ( Table 1): Derivation of the 'Global model matrix', through the 'Connect' function from Matlab (the choice of the output variables is not unique):

Ripple Effect on the Control Loops
As done in previous sections, it is important to analyze how the ripple present at the measured currents affects the control behavior. Once again, the voltage ripple analysis is not repeated here, since it is already studied in the first section and is only affected by the feed-forward terms, or in other words, it is not influenced by the feedback terms. Therefore, only how the current ripple affects the control has to be evaluated. This effect is numerically evaluated by the Bode diagrams depicted in Figure 25. Thus, the current ripple impacts the output converter voltage through terms C ααp , C αβp , C βαp , C ββp and C ααn , C αβn , C βαn , C ββn for harmonic n = 1. As can be noticed, the attenuation at the frequencies where the ripple appears (multiples of f sw ) is −9.11 dB for terms C ααp , C ββp , C ααn , and C ββn , which is something equal to what was obtained in the previous control case, since nothing has changed for these loops. However, for loops for harmonic h 1 = 13, it is seen that there are terms (C αβp1 , C βαp1 , C αβn1 , and C βαn1 ) that even amplify the current harmonics at 4.26 dB in this numerical example. It occurs mainly due to the terms that are incorporated to cancel the current coupling terms, which include a term such as ω n1 × L = 13·2 × π × 50 × L = 1.633, which produces a direct amplification of the current, at least for this studied numeric case. Therefore, to avoid the problematic control situations that can bring this fact, one possible solution could be to include a filter at the current measurement, but only at the loops where this problem arises, i.e., at the loops for h 1 = 13 harmonics.
Energies 2020, 13, x FOR PEER REVIEW 30 of 42 frequencies where the ripple appears (multiples of fsw) is −9.11 dB for terms Cααp, Cββp, Cααn, and Cββn, which is something equal to what was obtained in the previous control case, since nothing has changed for these loops. However, for loops for harmonic h1 = 13, it is seen that there are terms (Cαβp1, Cβαp1, Cαβn1, and Cβαn1) that even amplify the current harmonics at 4.26 dB in this numerical example. It occurs mainly due to the terms that are incorporated to cancel the current coupling terms, which include a term such as ωn1 × L = 13·2 × π × 50 × L = 1.633, which produces a direct amplification of the current, at least for this studied numeric case. Therefore, to avoid the problematic control situations that can bring this fact, one possible solution could be to include a filter at the current measurement, but only at the loops where this problem arises, i.e., at the loops for h1 = 13 harmonics.

Notch Filters Performance
In this control, each loop needs three filter notches to cancel oscillations at steady state at three different frequencies (Table 11). This means that the effect of the three notches together must be analyzed in detail instead of separately. Thus, Figure 26 shows the time-domain step response and the Bode diagrams of the three notches involved at the 50 Hz loops 2 . Regarding the settling time and peak of the step response, it is seen that the dominant filter is the 100 Hz filter or

Notch Filters Performance
In this control, each loop needs three filter notches to cancel oscillations at steady state at three different frequencies (Table 11). This means that the effect of the three notches together must be analyzed in detail instead of separately. Thus, Figure 26 shows the time-domain step response and the Bode diagrams of the three notches involved at the 50 Hz loops F 2h (s) · F h1+1 (s) · F h1−1 (s). Regarding the settling time and peak of the step response, it is seen that the dominant filter is the 100 Hz filter or F 2h (s). As occurs with only one filter, to reduce the peak value, a low damping is needed, while to reduce the settling time, a big damping is needed in all the notches. On the other hand, a big damping means better attenuation of the harmonics at steady state.

Notch Filters Performance
In this control, each loop needs three filter notches to cancel oscillations at steady state at three different frequencies (Table 11). This means that the effect of the three notches together must be analyzed in detail instead of separately. Thus, Figure 26 shows the time-domain step response and the Bode diagrams of the three notches involved at the 50 Hz loops 2 . Regarding the settling time and peak of the step response, it is seen that the dominant filter is the 100 Hz filter or 2 ( ) h F s . As occurs with only one filter, to reduce the peak value, a low damping is needed, while to reduce the settling time, a big damping is needed in all the notches. On the other hand, a big damping means better attenuation of the harmonics at steady state.
Step responses and Bode diagrams of the product of the notch filters; 2 for 50 Hz loops, with different damping ratios. (a) Step responses, (b) Bode diagrams.
Something equivalent occurs with the damping tendencies, with the product of three filters dedicated to the 650 Hz loops: F 2h (s) · F h1+1 (s) · F h1−1 (s), as is shown in Figure 27. Note that in this case, the time-domain step response is much faster than in the previous case, which would mean in principle that it has less of an effect on the dynamics of the corresponding current loops. Something equivalent occurs with the damping tendencies, with the product of three filters dedicated to the 650 Hz loops: 2 , as is shown in Figure 27. Note that in this case, the time-domain step response is much faster than in the previous case, which would mean in principle that it has less of an effect on the dynamics of the corresponding current loops. Finally, the choice of the damping ratio for each notch would be made by the multi-objective algorithm. In principle, for simplicity, all the notches will be tuned with the same damping ratio, and the actual value will be one degree of freedom to be searched by the algorithm. Step responses and Bode diagrams of product of notch filters: 2 1

Analysis of the Poles Location
In this subsection, the pole location of the global systems is analyzed in the αβ reference frame as done with the previous controls. In this specific control case, there is one fact that is quite important. This system present a pole pair (around 0.727 + 0.674i with ωn = 4.19 × 10 3 rd/s, its exact location depends on the exact values of Kp, Tn, Kp1, Tn1, etc.) very close to the unit circle, which means that this pole pair is very close to the instability. When a filter is included in the current measurements (Fi(s)), it has been seen that this pole pair tends to become instable very easily (module bigger than 1), being almost impossible to find a set of PI parameters that can stabilize the system. Having seen this and taking into Finally, the choice of the damping ratio for each notch would be made by the multi-objective algorithm. In principle, for simplicity, all the notches will be tuned with the same damping ratio, and the actual value will be one degree of freedom to be searched by the algorithm.

Analysis of the Poles Location
In this subsection, the pole location of the global systems is analyzed in the αβ reference frame as done with the previous controls. In this specific control case, there is one fact that is quite important. This system present a pole pair (around 0.727 + 0.674i with ω n = 4.19 × 10 3 rd/s, its exact location depends on the exact values of K p , T n , K p1 , T n1 , etc.) very close to the unit circle, which means that this pole pair is very close to the instability. When a filter is included in the current measurements (F i (s)), it has been seen that this pole pair tends to become instable very easily (module bigger than 1), being almost impossible to find a set of PI parameters that can stabilize the system. Having seen this and taking into consideration that this system needs a current measurement filter at least at the 650 Hz loops, to avoid problems with the current ripple, a current measurement filter has been included F LP (s) only at loops of 650 Hz and in the dq reference frame, as depicted in Figure 28 (see Tables 9 and 10). In this case, there are 18 pole pairs, which means that the system has become much more complex, since its behavior depends on many poles, and at the same time, their location depends on the coupling effect of all the control and physical parameters of the system. The dependence or sensitivity on main parameters is also shown in the figure (the filter's parameter dependence has not been included for simplicity). In comparison with the control analyzed in the previous section, it can be said that the poles of this control present a greater interdependence or coupling to all the parameters, and also in general, there are some poles closer to the unit circle. It is also curious but not so relevant that there are a few poles whose locations are not affected by the parameters evaluated. Thanks to this filter inclusion, it has been observed that this problematic pole pair becomes more stable and accepts control parameter variations without becoming unstable. In addition, as seen in Figure 28, the filtering capacity of the control loops at 650 Hz has been significantly improved. Therefore, a "robust stabilization" of the system has been achieved on the one hand, incorporating the necessity to attenuate the converter's harmonics on the other hand. Hence, the poles and zeros of the global closed loop system in the αβ reference frame are depicted in Figure 29, including the current measurement filter as commented.
In this case, there are 18 pole pairs, which means that the system has become much more complex, since its behavior depends on many poles, and at the same time, their location depends on the coupling effect of all the control and physical parameters of the system. The dependence or sensitivity on main parameters is also shown in the figure (the filter's parameter dependence has not been included for simplicity). In comparison with the control analyzed in the previous section, it can be said that the poles of this control present a greater interdependence or coupling to all the parameters, and also in general, there are some poles closer to the unit circle. It is also curious but not so relevant that there are a few poles whose locations are not affected by the parameters evaluated.
the coupling effect of all the control and physical parameters of the system. The dependence or sensitivity on main parameters is also shown in the figure (the filter's parameter dependence has not been included for simplicity). In comparison with the control analyzed in the previous section, it can be said that the poles of this control present a greater interdependence or coupling to all the parameters, and also in general, there are some poles closer to the unit circle. It is also curious but not so relevant that there are a few poles whose locations are not affected by the parameters evaluated. Figure 29. Poles and zeros for the global closed loop system in the αβ reference frame for step input references, corresponding to the control of Figure 24, deduced from the mathematical model of Tables  9 and 10, and applied to the numerical scenario of Appendix C. It is assumed that there is a delay in the converter voltage update and the voltage and current filters are not included, i.e., Fi(z) = 1, Fv(z) = 1. Dependence on parameters has been included as well (there are some poles whose location do not depend on the evaluated parameters).

Multi-Objective Tuning
This section presents one of the possible solutions of multi-objective tuning for the control loops. Figure 30 schematically illustrates how the tuning algorithm could be conceived. It follows the same idea of the previously presented algorithm for the previously studied control with the equivalent six objectives and the same inputs. Obviously, there are some other new parameters that have to be searched: Kp1, Tn1, ξLP and fLP. Figure 29. Poles and zeros for the global closed loop system in the αβ reference frame for step input references, corresponding to the control of Figure 24, deduced from the mathematical model of Tables 9 and 10, and applied to the numerical scenario of Appendix C. It is assumed that there is a delay in the converter voltage update and the voltage and current filters are not included, i.e., F i (z) = 1, F v (z) = 1. Dependence on parameters has been included as well (there are some poles whose location do not depend on the evaluated parameters).

Multi-Objective Tuning
This section presents one of the possible solutions of multi-objective tuning for the control loops. Figure 30 schematically illustrates how the tuning algorithm could be conceived. It follows the same idea of the previously presented algorithm for the previously studied control with the equivalent six objectives and the same inputs. Obviously, there are some other new parameters that have to be searched: K p1 , T n1 , ξ LP and f LP .
An explanation of the objective, the corresponding mathematical expression to be evaluated, an indicative numeric value for each objective, and an example of code program in Matlab for implementation are provided in Table 12. As commented in the previous section, the most determinant objective is the one that tries to find a module minimization of the pole with a bigger module (dominant pole) in the "z" domain. Note that a multi-objective optimization algorithm that searches the control parameters autonomously is a very suitable tool for a system so complex with such a large number of poles and so coupled or dependent on the parameters of the system. In other words, if a tool similar to this is not used in such coupled and complex systems, the tuning of the regulators "by hand" becomes very difficult. An explanation of the objective, the corresponding mathematical expression to be evaluated, an indicative numeric value for each objective, and an example of code program in Matlab for implementation are provided in Table 12. As commented in the previous section, the most determinant objective is the one that tries to find a module minimization of the pole with a bigger module (dominant pole) in the "z" domain. Note that a multi-objective optimization algorithm that searches the control parameters autonomously is a very suitable tool for a system so complex with such a large number of poles and so coupled or dependent on the parameters of the system. In other words, if a tool similar to this is not used in such coupled and complex systems, the tuning of the regulators "by hand" becomes very difficult.  Finally, it has to be remarked that having the mathematical model, the multi-objective optimization can be adapted to the needs or even style of the designer. Thus, for instance, a quite typical choice often is based on using several minimization objectives by defining a minimization cost function. Thus for instance, the minimization could incorporate the residue of the biggest pole also, making sure that is the dominant pole of the system, as its dominance is not reduced by a zero that could be near. Another possible choice, could be to minimize the module of the two (or three) more dominant poles, again to make sure that the system dynamic really depends on the dominance of the poles selected. Many other possible valid solutions could be selected as well.
To conclude, Figure 31 illustrates a conceptual block diagram that shows the procedure that can be followed to search the optimum set of parameters for the current control loops studied in this section. It can be implemented in a unique Matlab program, and it searches the optimum set of parameters (k p , T n , k p1 , T n1 , ξ, f LP . . . ) according to the defined multi-objectives.

Example 1: Search the Minimum Module Pole
In this subsection, an example that applies the proposed multi-objective optimization is shown. For simplicity in the exposition and in showing the results, first, a limited range and precision to the parameters is provided. The summary of the control parameters evaluated by the algorithm is depicted in Figure 32. The seven control parameters are degrees of freedom from the designer point of view and are varied or swept, as shown in the figure. Note that for simplicity, the same damping has been set in all the notch filters, easing the search because a unique damping has to be found. Then, at each control parameter set evaluated, the biggest module pole is collected as well. Checking the points evaluated, it is seen that every single point falls in a stable situation, with a module of the largest pole within the range [0.995, 0.990] approximately. The fastest (or further from instability) corresponds to parameters: K p = 0.266, T n = 0.0026055, K p1 = 0.55125, T n1 = 0.00096012, ξ Notch = 0.084, ξ LP = 1.8, f LP = 960 Hz, and |z| min = 0.9909. Then, it is also checked that this point meets the rest of the specifications of Table 12 (not included this information at the graph for an easier understanding and simplicity).
The pole-zero map resulting from these control parameters is shown in Figure 33, together with the time-domain response (Figure 34), validating its feasibility and good agreement between the theoretical model and the Simulink blocks-based model (the same simulations are carried out as in previous sections). Note that depending on the desired closed loop dynamics of the system, the designer can choose the most appropriate set of control parameters.
In relation to this, Figure 35 shows the spectrum of the currents exchanged with the grid at different f sw , i.e., the current after the grid-side filter that goes directly to the grid. It can be noticed that the quality of this current is improved (less undesired harmonics) with bigger switching frequencies.
The choice of the switching frequency should be made finding a balance in aspects such as the quality of currents exchanged with the grid (level of harmonics at the AC side), power losses of the power converter (determined by used semiconductors and the cooling system), level of harmonics at the DC side (amount of capacitance used at the DC side of the converter), stability/dynamic performance of the closed loop system (the pole locations are dependent on f sw or T s ), and many other aspects.
On the other hand, refining the range of parameters, with a narrower precision and with a higher range of variation, after a massive search evaluating many points, it is probably possible to find a set of parameters that reduces the module of the biggest pole. However, the search in this set of results is not done for simplicity and because the already obtained result could be good enough for many realistic application scenarios. This search by itself could be based on a parallel (or coordinated) multi-objective optimization that is out of the scope of this article.
Finally, as an illustrative example, Figure 36 shows a range of variation of control parameters, where stable and unstable points are found. Note that the proposed tool based on the mathematical model facilitates considerably the search and reduces the computational cost of searching the optimum set of parameters. By a simple "click", it is able to evaluate thousands of control parameters. <0.08 F 50 _inf.SettlingTime

F 650 _inf.SettlingTime
Set a maximum limit to the overshoot of the notch filters to the step input (in the dq reference frame)

F 650 _inf.Peak
Set a maximum limit to the effect of the ripple at fsw in the control loops (in the αβ reference frame) <0.001 Mag100_F 50 ,phase,w]= bode(F 50 , 100 × 2 × pi); To conclude, Figure 31 illustrates a conceptual block diagram that shows the procedure that can be followed to search the optimum set of parameters for the current control loops studied in this section. It can be implemented in a unique Matlab program, and it searches the optimum set of parameters (kp, Tn, kp1, Tn1, ξ, fLP…) according to the defined multi-objectives. Figure 31. Conceptual block diagram showing the procedure that can be followed to search the optimum set of parameters for the current control loops. For instance, it can be implemented together with the validation in a unique Matlab program. This procedure can be applied equivalently to the previously studied current controls.

Example 1: Search the Minimum Module Pole
In this subsection, an example that applies the proposed multi-objective optimization is shown. For simplicity in the exposition and in showing the results, first, a limited range and precision to the  Table 12 (not included this information at the graph for an easier understanding and simplicity). The pole-zero map resulting from these control parameters is shown in Figure 33, together with the time-domain response (Figure 34), validating its feasibility and good agreement between the theoretical model and the Simulink blocks-based model (the same simulations are carried out as in previous sections). Note that depending on the desired closed loop dynamics of the system, the designer can choose the most appropriate set of control parameters.  The pole-zero map resulting from these control parameters is shown in Figure 33, together with the time-domain response (Figure 34), validating its feasibility and good agreement between the theoretical model and the Simulink blocks-based model (the same simulations are carried out as in previous sections). Note that depending on the desired closed loop dynamics of the system, the designer can choose the most appropriate set of control parameters. In relation to this, Figure 35 shows the spectrum of the currents exchanged with the grid at different fsw, i.e., the current after the grid-side filter that goes directly to the grid. It can be noticed that the quality of this current is improved (less undesired harmonics) with bigger switching frequencies. Figure 33. Pole-zero map (i α (z)/i dp *(z)) corresponding to the control parameters with a minimum module pole resulting from the previous Figure 31: K p = 0.266, T n = 0.0026055, K p1 = 0.55125, T n1 = 0.00096012, ξ Notch = 0.084, ξ LP = 1.8, f LP = 960 Hz, |z| min = 0.9909. It is also checked that this point meets the rest of the specifications of Table 12 (this information is not included in the graph for simplicity).

Example 2: Robustness under Uncertainty of Parameters
In this example, how the system behaves under uncertainty of some of the parameters is studied. In this example, only two parameters are considered that are not accurately known. The first parameter is the L g value, which in real applications depends on the grid equivalent impedance. This impedance is not always known accurately. Despite the fact that the control does not use directly this parameter, it affects the stability of the entire system, as has been seen in previous sections. On the other hand, the other parameter that in this example is not precisely known is the filter inductance L. Note that this parameter is the only one that the control is using online through the cancelation of the coupling terms. Therefore, although it is not very common to have large uncertainty in this parameter, since it has such a direct impact on all the control loops, it is important to quantify how it affects the system.
Hence, the tuning algorithm is applied to the numeric scenario of Appendix C, where the parameters L g and L are varied in a range of ±5% from its actual value (uncertainty of ±5%). Both parameters are varied within this range, while the control uses the constant value of L = 400 µH. The result of the search is shown in Figure 37. As in the previous example, for simplicity in the exposition and in showing the graphical results, the steps in the variation of the parameters are not very small. In a real application, if much more accurate results are needed, it is recommended to evaluate bigger ranges of variations with bigger precision (smaller steps in the variation ranges). Note also that in this particular example, again for simplicity, the damping of the notch filters has been left constant to ξ = 0.084.
Analyzing the results, as the variation of the parameters is very similar to those presented in the previous example, it can be deduced that the L g and L variations do not significantly affect the poles with the bigger module. This means that the dynamic of the overall system as well as the stability is not affected significantly either, which means that is a quite robust system, at least within the evaluated uncertainty range. It has to be emphasized that the proposed modeling methodology permits evaluating and quantifying these facts whenever necessary. It can be noticed that the minimum module poles, as seen in the zoomed figure, are obtained with "wrong" values of L and L g .
The choice of the switching frequency should be made finding a balance in aspects such as the quality of currents exchanged with the grid (level of harmonics at the AC side), power losses of the power converter (determined by used semiconductors and the cooling system), level of harmonics at the DC side (amount of capacitance used at the DC side of the converter), stability/dynamic performance of the closed loop system (the pole locations are dependent on fsw or Ts), and many other aspects.   On the other hand, refining the range of parameters, with a narrower precision and with a higher range of variation, after a massive search evaluating many points, it is probably possible to find a set of parameters that reduces the module of the biggest pole. However, the search in this set of results is not done for simplicity and because the already obtained result could be good enough for many realistic application scenarios. This search by itself could be based on a parallel (or coordinated) multi-objective optimization that is out of the scope of this article.
Finally, as an illustrative example, Figure 36 shows a range of variation of control parameters, where stable and unstable points are found. Note that the proposed tool based on the mathematical model facilitates considerably the search and reduces the computational cost of searching the optimum set of parameters. By a simple "click", it is able to evaluate thousands of control parameters.

Example 2: Robustness under Uncertainty of Parameters
In this example, how the system behaves under uncertainty of some of the parameters is studied. In this example, only two parameters are considered that are not accurately known. The first parameter is the Lg value, which in real applications depends on the grid equivalent impedance. This impedance is not always known accurately. Despite the fact that the control does not use directly this parameter, it affects the stability of the entire system, as has been seen in previous sections. On the other hand, the other parameter that in this example is not precisely known is the filter inductance L. Note that this parameter is the only one that the control is using online through the cancelation of the coupling terms. Therefore, although it is not very common to have large uncertainty in this parameter, since it has such a direct impact on all the control loops, it is important to quantify how it affects the system. On the other hand, refining the range of parameters, with a narrower precision and with a higher range of variation, after a massive search evaluating many points, it is probably possible to find a set of parameters that reduces the module of the biggest pole. However, the search in this set of results is not done for simplicity and because the already obtained result could be good enough for many realistic application scenarios. This search by itself could be based on a parallel (or coordinated) multi-objective optimization that is out of the scope of this article.
Finally, as an illustrative example, Figure 36 shows a range of variation of control parameters, where stable and unstable points are found. Note that the proposed tool based on the mathematical model facilitates considerably the search and reduces the computational cost of searching the optimum set of parameters. By a simple "click", it is able to evaluate thousands of control parameters.

Example 2: Robustness under Uncertainty of Parameters
In this example, how the system behaves under uncertainty of some of the parameters is studied. In this example, only two parameters are considered that are not accurately known. The first parameter is the Lg value, which in real applications depends on the grid equivalent impedance. This impedance is not always known accurately. Despite the fact that the control does not use directly this parameter, it affects the stability of the entire system, as has been seen in previous sections. On the other hand, the other parameter that in this example is not precisely known is the filter inductance L. Note that this parameter is the only one that the control is using online through the cancelation of the coupling terms. Therefore, although it is not very common to have large uncertainty in this parameter, since it has such a direct impact on all the control loops, it is important to quantify how it affects the system. Hence, the tuning algorithm is applied to the numeric scenario of Appendix C, where the parameters Lg and L are varied in a range of ±5% from its actual value (uncertainty of ±5%). Both parameters are varied within this range, while the control uses the constant value of L = 400 μH. The result of the search is shown in Figure 37. As in the previous example, for simplicity in the exposition and in showing the graphical results, the steps in the variation of the parameters are not very small. In a real application, if much more accurate results are needed, it is recommended to evaluate bigger ranges of variations with bigger precision (smaller steps in the variation ranges). Note also that in this particular example, again for simplicity, the damping of the notch filters has been left constant to ξ = 0.084.
Analyzing the results, as the variation of the parameters is very similar to those presented in the previous example, it can be deduced that the Lg and L variations do not significantly affect the poles with the bigger module. This means that the dynamic of the overall system as well as the stability is not affected significantly either, which means that is a quite robust system, at least within the evaluated uncertainty range. It has to be emphasized that the proposed modeling methodology permits evaluating and quantifying these facts whenever necessary. It can be noticed that the minimum module poles, as seen in the zoomed figure, are obtained with "wrong" values of L and Lg.

Example 3: Robustness under Uncertainty of All Parameters
Finally, this example illustrates a quite typical and realistic situation. Once the algorithm has found the set of control parameters fulfilling the desired specifications, they are implemented at the control, and the system operates in its functionality. Then, by means of the proposed mathematical model, it is possible to evaluate how the poles of the system are moved when the electric parameters of the system are varied as well within a range (because their values are not accurately known or because they change their value due to phenomena such as the temperature, saturation, degradation, etc.). This analysis is often also known as evaluation of the robustness of the system under an uncertainty of parameters. Therefore, this evaluation is shown in Figure 38. The set of control parameters used are specified in the legend of the figure (search performed in Section 4.8), while the range of uncertainty of the electric parameters of the system has been set for all of them ±10% in first test, ±5% in second test, and ±2.5% in third test: L, Cc, Lg, and Rc. It is seen how the poles are moved in this circumstances, deducing also how the stability and

Example 3: Robustness under Uncertainty of All Parameters
Finally, this example illustrates a quite typical and realistic situation. Once the algorithm has found the set of control parameters fulfilling the desired specifications, they are implemented at the control, and the system operates in its functionality. Then, by means of the proposed mathematical model, it is possible to evaluate how the poles of the system are moved when the electric parameters of the system are varied as well within a range (because their values are not accurately known or because they change their value due to phenomena such as the temperature, saturation, degradation, etc.). This analysis is often also known as evaluation of the robustness of the system under an uncertainty of parameters.
Energies 2020, 13, 5849 44 of 47 Therefore, this evaluation is shown in Figure 38. The set of control parameters used are specified in the legend of the figure (search performed in Section 4.8), while the range of uncertainty of the electric parameters of the system has been set for all of them ±10% in first test, ±5% in second test, and ±2.5% in third test: L, C c , L g , and R c . It is seen how the poles are moved in this circumstances, deducing also how the stability and dynamic performance of the system is conditioned for the specific uncertainty that has been evaluated. Thus, with the control parameters that have been selected, for an uncertainty in the electric parameters of ±2.5%, the stability of the system is guaranteed. For an uncertainty in the electric parameters of ±10%, the stability is not ensured, since there are some combinations of parameter variations that locate some of the poles outside the unit circle. An uncertainty of ±5% is approximately the limit situation.

Summary of Conclusions of the Section
The main conclusions of the analysis performed with the proposed mathematical model can be summarized as follows. It is detected that the control loops dedicated to harmonic h1 are quite sensible to the switching harmonics, so it is necessary to include a low-pass filter on them. In relation to the location of the poles of this control system, it is seen that it is not easy to identify a few dominant poles. There are many pole pairs that can be considered as dominant (high module amplitude), and moreover, their dependence on the system's parameters is much more complex and coupled than in the previous controls studied. This fact strongly justifies the necessity of the proposed multi-objective tuning method, which autonomously searches the best control gains of the complex and interdependent control loops. Finally, the searched control gains are evaluated under several examples of uncertainty of parameters. In a given numerical example, we studied under which uncertainty level the control gains searched guarantee stability and under which uncertainty level the control falls on instability.

Conclusions
Designing and tuning controls that incorporate many current control loops requires selecting many associated gains and parameters. For instance, the most complex control that has been analyzed in this article (active filter compensating only one harmonic) presents eight loops, and at least six parameters must be tuned (even 10 parameters could be necessary to be tuned if the search wants to find different optimized damping factors of the filters). If the specific active filter application requires Figure 38. Robustness under uncertainty of all electric parameters. The control parameters are set constant (K p = 0.266, T n = 0.0026055, K p1 = 0.55125, T n1 = 0.00096012, ξ Notch = 0.084, ξ LP = 1.8, f LP = 960 Hz, |z| min = 0.9909) and all the electric parameters; L, C c , L g , R c , R, and R g are varied; (a) ±10% range (steps 5%), (b) ±5% range (steps 2.5%), (c) ±2.5% range (steps 1.25%).

Summary of Conclusions of the Section
The main conclusions of the analysis performed with the proposed mathematical model can be summarized as follows. It is detected that the control loops dedicated to harmonic h 1 are quite sensible to the switching harmonics, so it is necessary to include a low-pass filter on them. In relation to the location of the poles of this control system, it is seen that it is not easy to identify a few dominant poles. There are many pole pairs that can be considered as dominant (high module amplitude), and moreover, their dependence on the system's parameters is much more complex and coupled than in the previous controls studied. This fact strongly justifies the necessity of the proposed multi-objective tuning method, which autonomously searches the best control gains of the complex and interdependent control loops. Finally, the searched control gains are evaluated under several examples of uncertainty of parameters. In a given numerical example, we studied under which uncertainty level the control gains searched guarantee stability and under which uncertainty level the control falls on instability.

Conclusions
Designing and tuning controls that incorporate many current control loops requires selecting many associated gains and parameters. For instance, the most complex control that has been analyzed Energies 2020, 13, 5849 45 of 47 in this article (active filter compensating only one harmonic) presents eight loops, and at least six parameters must be tuned (even 10 parameters could be necessary to be tuned if the search wants to find different optimized damping factors of the filters). If the specific active filter application requires including a compensation of for instance three harmonics, which is quite common, the number of necessary loops would be 16, and the number of parameters that must be tuned would be at least 12 (more than 20 if the damping coefficients are included in the search). Hence, in order to have a clear idea of the performance of such complex controls with so many loops interacting one with each other, a manageable and at the same time powerful model is very useful.
It has been seen that with the proposed mathematical model in the αβ reference frame (the only common frame for all the loops), it is possible to accurately reproduce the behavior of the control and perform several types of analysis. The pole location analysis selected in this paper and the associated multi-objective tuning method have integrated the information related to useful design criteria such as the dynamic behavior, stability margin, robustness under uncertainty of parameters, and so on. In addition, it helps knowing which poles are more susceptible to becoming unstable, and with an associated sensibility analysis, it is possible to know which parameters have more impact on their location and could be redesigned for a better performance. Added to this, checking also the pole locations, it is possible to identify poorly damped poles that produce transitory oscillations that could be problematic and then effectively modify the necessary parameters to correct them. Thus, it can be concluded that the pole location analysis applied to this context is a powerful tool, since it provides a perspective of the global behavior of a complex control in a simple manner. In addition, by multi-parameter swept, it is possible to know all the possibilities or behaviors that the control can have, and with the multi-objective search, it is possible to find the most suitable and desired requirements.