Robust Tuning of Multiresonant Current Controllers for Grid-Tied Converters and Erroneous Use of the Naslin Polynomial Method

The majority of the paper is devoted to debunking flawed methods proposed in the literature within the context of multiresonant control systems for grid-connected converters and then a novel approach is proposed. Several flawed approaches to grid tied power electronic converters (rectifiers and shunt active power filters) are discussed to initiate a critical discussion regarding some tuning methods reported in the topical literature. In this paper we analytically show that some of them are even erroneous from the control theory point of view or are based on mathematically erroneous derivations. That is especially prevalent for papers in which the Naslin polynomial method is employed to tune proportional-multiresonant controllers. Some authors recognize that the resulting control systems are not robust, i.e. not practical, and use obtained gains only as a starting point for further tuning using the trial and error method or evolutionary global optimization. Robustness of such systems is often not guaranteed and if it occurs at all, it is by chance or thanks to a combination of luck and expert knowledge of the engineer (not by deliberate design). Therefore, the absence of such a guarantee motivates the search for better ways to tune such controllers. As a result, a practical robust controller tuning method for multiresonant grid current controllers is proposed and verified experimentally. We are convinced that one of the solutions can be based on the disk margin stability analysis and a global search algorithm. The required robustness is expressed as stability margins. The design procedure is very user-friendly and the disk size is the key design parameter to be selected by an engineer. The practicality of the tuning method is demonstrated empirically in a 10 kVA physical grid-tied converter.


SSE
sum of squared errors. VOC voltage oriented control. VSC voltage source converter.

I. INTRODUCTION
Multiresonant controllers have become an industrial standard in power electronic converters interfacing grid-tied systems nowadays. Such controllers are used to comply with harmonic distortion limits specified for electrical appliances [1], [2]. Multiresonant controllers are adopted to enable sinusoidal currents regardless of the distorted grid voltage harmonic content. Selected frequencies are rejected by incorporating relevant oscillators into the controller. The number of the resulting controller gains to be tuned grows linearly with the number of frequencies to be dealt with. The problem is fairly trivial for a proportional-resonant (PR) controller with two gains to be determined, whereas it becomes quite challenging for a proportional-multiresonant (P-MR) controller. The problem is also equally difficult for repetitive controllers [3]. The proposed tuning procedures often require good expert knowledge to get robust, thus practical, controller gains. That is due to the presence of many parameters that have to be selected, usually by guessing and checking, for the tuning algorithm itself to make it produce robust gains for a given plant at hand. Those optimisation methods tend to shift the guessing and checking activity from direct guessing controller gains to guessing successful settings for the tuning procedure. This motivates some (if not many or even most) practitioners to go back to direct guessing controller gains for the multioscillatory part. It should not be concluded though, that among the reported analytical tuning methods there are no potential candidates for becoming industrial standards.
Regarding grid current control of grid-tied converters, the most commonly discussed methods of solving the problem of the multi-resonant (multi-oscillatory) controller tuning are probably pole placement and linear quadratic regulator [4], [5], [6]. These methods require good expert knowledge and they do not necessarily ensure the expected robustness of the system. Therefore, some designers try to use LQG and LMI to address this issue [7], [8]. It is also possible to use global optimization tools such as PSO to support a non-intuitive selection of weighting matrices in a multi-oscillatory current controller designed using the LQR method [9]. Alternatively, online optimisation may be applied as in model predictive control (MPC) [10], [11], still taking into account that such methods may also be parameter dependant and often are very computation intense. The latter has been addressed in numerous publications and several improvements reducing the computational burden have been proposed, e.g. [12].
With all of the above in mind, it is important to remember that our goal should be not just to do things differently, but first and foremost innovatively, and by innovative solutions we mean not just novel ones but first and foremost practical and improving ways we do things. At the same time it is worth mentioning that innovation often comes from applying already well-established tools in a slightly different arrangement to solve a given problem more effectively. And this is the case for the proposed method in which the disk margin analysis and the particle swarm optimizer are employed to determine robust settings for the multioscillatory controller. Yet, this is discussed in the second part of the paper. The first part of the paper is devoted to different reported attempts to tune multioscillatory controllers using the Naslin polynomial method. That is to serve two purposes. First of all, several of these methods are flawed (if not entirely erroneous) and they do not seem to have been assessed critically in the literature. As a consequence, consecutive generations of engineers try to employ them in their systems to little or no practical effect, simply wasting their time. What is also worrying, new papers that claim to use these approaches are still being published. Our second purpose in pointing out these flawed methods is to illustrate that the problem is not trivial and there is a need for more thorough discussion regarding tuning multiresonant controllers in the context of power electronic converters. We would like to initiate that discussion by presenting our most successful (up to now) attempt to produce a multioscillatory control system for a grid-tied converter with guaranteed robustness measured in stability margins.

II. OVERVIEW OF THE CONTROL SYSTEM
This section presents a model of the system and the control system design process. The physical plant consists of a voltage source converter (VSC) connected to the grid using an L-type filter (Fig. 1). Selected parameters of the system (for the stage of simulation and physical experiments) are given in Table 1. The three-level neutral point clamped (NPC) topology is used; therefore, the DC-link voltage (v dc1 + v dc2 ) is divided with a use of two capacitors of equal capacity. The DC-link capacitor voltage balancing strategy is realized by a proper pulse width modulation (PWM) algorithm [13]. The analyzed control system of the converter is based on voltage oriented control (VOC). The applied cascade control system involves two control loops. The inner grid current loop is equipped with a proportional-integral-multiresonant (PI-MR) controller, while the outer DC-link voltage loop includes a proportional-integral (PI) controller. The grid current feedback signals are transformed to the dq frame using angle θ received from a phase lock loop (PLL) based on the delayed signal cancellation method (DSC-PLL) [14].
The main objective of using the PI-MR current controller is to obtain sinusoidal symmetrical grid currents under distorted grid voltage conditions. The 5 th , 7 th , 11 th and 13 th harmonics of the grid voltage usually have the largest amplitudes and the L filter has the weakest attenuation for these harmonics. In order to reduce the sinusoidal component of the error in current tracking, caused by this type of voltage distortion, the resonant controllers are incorporated. The compensation for the unbalance and the 5 th , 7 th , 11 th and 13 th harmonics in the natural reference frame takes place in the dq frame through the compensation of the 2 nd , 6 th and 12 th harmonics of the grid currents. A gain selection process for the PI-MR current controller is the key step to which we devote most of our effort in this paper.
The outer DC-link voltage control loop includes a proportional-integral (PI) controller. In this loop, there are the following filters: a low-pass (LPF) one and two bandstop (NF 2ω , NF 6ω ) ones with the resonant pulsatances of 2ω and 6ω. The DC-link voltage controller is tuned with a use of the symmetrical optimum method as shown in detail in [15].

III. NASLIN POLYNOMIAL METHOD
The Naslin polynomial method uses as the reference transfer function for a closed loop system [16], where ω n = a n a n+1 (2) and are the characteristic pulsatances and the characteristic ratios. More on the characteristic ratios and the characteristic pulsatances can be found in [17], [18], and [19]. The basic original method then gives all the characteristic ratios a common value α, which possesses the value of a genuine damping factor (α = 4ζ 2 ) in the second order system. The resulting coefficients for n = 1 . . . N are then Such systems have very similar step responses for a given α, regardless of their order N , which is demonstrated in Fig. 2. In this way we obtain a family of standard systems with characteristic polynomials with adjustable damping. These polynomials are then used to design controller gains by solving a system of equations resulting from equating coefficients of the characteristic polynomial of the closed-loop system to the coefficients of the reference characteristic polynomial with adjustable damping ζ . The key remark here is that the method is successful only if the system of equations has a solution, i.e. is consistent.

IV. PROPORTIONAL-RESONANT CURRENT CONTROLLER
For the grid-tied converter with an L filter shown in Fig. 1 a simplified transfer function of the plant to be controlled in terms of the current control loop is where the pulse width modulator is configured to introduce a unity gain along with the converter. The simplification neglects the delay introduced by the modulator and the digital controller. Applying a proportional-resonant controller we get the open loop transfer function The characteristic polynomial of the closed loop system is then The relevant reference Naslin polynomial has coefficients given by (4) and is as follows Polynomials (9) and (10) are equivalent iff (if and only if) where L > 0, R > 0, ω 1 > 0 and α > 0 (often set to 2) are specified by the designer (known parameters), whereas a 0 , ω 0 > 0, k P and k 1 are unknowns to solve for. Note that the system of equations is exactly determined, i.e. having the same number of equations as unknowns. Substituting k P + R = a 0 ω 2 1 from the last equation we get Substituting a 0 = α 3 ω 3 0 L we get the PR controller gains   where ω 0 = ω 1 α −0.5 . This result can also be found in [20] (pages 140 and 141) along with a numerical example. The conclusion is that the Naslin polynomial method can be indeed deployed for tuning the PR current controller in PWM rectifiers. The same page 141 of [20] states that the presented design procedure can be extended for two or more different resonant frequencies, resulting in a multiresonant controller. Unfortunately, that statement is incorrect. We will demonstrate in the next section that the method already fails in the case of the proportional-multiresonant controler with only two resonant terms.

V. PROPORTIONAL-MULTIRESONANT CURRENT CONTROLLER
Applying a proportional-resonant controller to control the plant described by the transfer function (5) The characteristic polynomial of the closed loop system is then D closed (s) = k P (s 2 + ω 2 1 )(s 2 + ω 2 2 ) +sk 1 (s 2 + ω 2 2 ) + sk 2 (s 2 + ω 2 1 ) Arranging by descending order gives D closed (s) = Ls 5 + (k P + R)s 4 The matching 5 th order Naslin polynomial has coefficients given by (4) and is as follows 88214 VOLUME 10, 2022 Polynomials (17) and (18) where L > 0, R > 0, ω 1 > 0, ω 2 > 0 and α > 0 (often set to 2) are specified by the designer (known parameters), whereas a 0 , ω 0 > 0, k P , k 1 and k 2 are unknowns to solve for. From the algebraic point of view, the case is already substantially different in comparison with the previous case of the PR controller. An this is because the system of equations is not exactly determined any more -it is overdetermined, i.e. having more equations (to be specific, six of them) than unknowns (to be specific, five of them). Overdetermined systems may still be consistent. Nevertheless, especially in real life engineering problems, they should be regarded as strong candidates for inconsistent systems of equations, thus having no set of values for the unknowns that satisfies all of the equations. Unfortunately, it is the case here, i.e. (19) is inconsistent and will not produce a feasible set of gains for the controller. As we are going to show the inconsistency of the parametrized system, we need to point out that inconsistency may be parameter dependent, i.e. there may exist specific sets of parameter values that turn the system into a consistent one. However, in our problem the designer has little to no room to adjust the parameters: L is designed for current ripples, R is parasitic, ω 1 > 0 and ω 2 > 0 are determined solely by the expected harmonic content of the grid voltage and the reference frame chosen for the control system, whereas α > 0 should be freely adjustable to get the desired damping, thus all these cannot be used to potentially turn the system of equations into a consistent one. We will demonstrate the inconsistency of (19) by showing two subsets of its equations that produce a disjoint set of solutions. Let these subsets be and Solving (20) for ω 0 > 0 we get Step response of the reference system with the correct Naslin polynominal compared with the step response of the system with incorrectly performed coefficient fitting.
whereas solving (21) for Pulsatances (22) and (23) would be the same pulsatance only for α = ω 2 1 +ω 2 2 ω 1 ω 2 , but as mentioned above, in practice we do not have the flexibility of choosing α as a function of ω 1 and ω 2 , because α is an independent design parameter that sets the damping for the closed loop system. Therefore, the conclusion is that the Naslin polynomial method cannot be used to tune the proportional-resonant controller with two oscillatory terms of the form (14). The severity of an error made when such equation omitting step is carried out to conceal inconsistency of the system of equations can be illustrated with the following example. Let us suppose L = 2 mH, R = 0.2 , ω 1 = 2π · 50 · 6 rad s −1 and ω 2 = 2π · 50 · 12 rad s −1 . If one chooses to calculate ω 0 from (22), it would result in a coefficient mismatch yielding the response shaped impractically far from the expected one as shown in Fig. 3.
Obviously, one cannot infer from this single example that adding more oscillatory terms to the controller will always result in inconsistent equation systems. Neverteless, we verified algebraically that P-MR with three resonant terms also introduces a contradiction. Moreover, the more resonant terms, the more overdetermined the resulting equation system becomes. This makes us believe that checking beyond three resonant terms would not be successful.
Unfortunately, the overoptimistic statement in a very popular handbook [20] has kind of a ripple effect. For example, in [21] the Naslin polynomial technique was used to tune P-MR controller with three resonant terms for an APF. The presented final formulas may suggest that the overdetermined and inconsistent equation system was incorrectly turned into VOLUME 10, 2022 an exactly determined and consistent one by omitting some equations. The error was repeated in [22] and potentially also in [23], where the Naslin polynomial tuning method was reported for the P-MR controller. Also other research groups report the Naslin polynomial method to be successful for P-MR controllers by means of some reiteration [24] without providing any details on inner workings of the procedure. Another example comes from [25], where a controller with two resonant terms was tuned reportedly using the technique described in [26], which is confusing because the latter presents formulas only for a proportional-integral-resonant controller with a single resonant therm. A step by step derivation for the PIR controller is presented in [27] and will not be repeated here.
Another peculiar attempt to use the Naslin polynomial method that goes against mathematical reasoning is presented in [28]. The PI-MR with four resonant term is being pretuned for a shunt APF with the help of the tenth order Naslin polynomial. The resulting system of equations is overdetermined: 10 equations and just 8 unknowns, including 6 controller gains for a general case and just 6 unknowns, including 4 controller gains for the presented case with PI gains tuned a priori using a different method. Despite that, the system is not checked for potential inconsistencies. In fact the inconsistencies were found but instead of being interpreted as such they are presented as viable solutions -all the six obtained controller gain quadruplets. The paper states: ''With a 0 , six values of ω 0 can be identified, since there are six equations related to ω 0 as a function of a 0 .'' An equation system was incorrectly interpreted as a logical disjunction of its equations, instead of being solved as a set of equations linked by a logical conjunction. The existence of six equations that produce at least two different values of ω 0 means that the system of equations is inconsistent and that for a given controller there are no settings able to match the coefficients given by (4), thus the Naslin polynomial method cannot be applied to tune the controller. This is similar to the case discussed in Section V. The existence of two different (inconsistent) formulas (22) and (23) nullifies the solution set. Nevertheless, it should be acknowledged that the overall tuning procedure presented in [28] may be successful. The Naslin method was used incorrectly at the pretuning stage. The main tuning stage involves differential evolution optimization that could be run even on randomized initial controller gains, thus it is also able to work with the ones produced by the erroneous derivation. Two stage tuning procedures are not uncommon among practitioners. If the Naslin polynomial method is to be deployed at the first stage, a PR (not P-MR) controller is to be used. An example of such a workflow is mentioned in [29], in which selected controllers are tuned employing the trial and error method with the starting point obtained from the Naslin polynomial approach. The second stage is there to increase robustness in the presence of less significant plant poles neglected during the first stage (the analytical solution assumes a first order plant (5), thus only the dominant pole is taken into account).

VI. SUPERPOSITION PRINCIPLE AND SUPERPOSITION OF CONTROLLERS
Another category of tuning procedures that go against the control theory and also often fail to deliver robust enough multiresonant controller settings involves superposition of individually tuned PR controllers to form a P-MR controller. These workflows are often flawed because they lack a critical assessment of amplitude-frequency characteristics of the individually tuned PR controllers. Supposably, the root cause for these errors is a misunderstood concept of superposition of controllers -presumably wrongly envisaged as something similar to the superposition principle. The superposition principle states that the net response of an LTI system caused by two or more stimuli is the sum of the responses that would have been caused by each stimulus individually [30]. For example, consider a linear time-invariant dynamic system G(s) fed by the signal then On the other hand, superposition of controllers does not usually produce a response of the closed-loop system that would be equal to the sum of responses for the individual controllers. For example, consider an LTI plant G p and LTI controllers G c = G c1 +G c2 connected to form three feedback systems all fed by the same reference signal r. If then in general To get exactly y = y 1 + y 2 one very impractical condition would have to be met, i.e. G c1 G c2 = 0 for all frequencies. This can be obtained from: which is held under G c1 G c2 = 0. This implies that the separate tuning of controllers aimed at their parallel operation in the target control system should be accompanied with an assessment of their selectiveness in the frequency domain.
In practice no controllers give G c1 G c2 = 0, thus the condition is replaced with G c1 G c2 ≈ 0 and the overlapping of an amplitude-frequency characteristic of the controllers should be intentionally kept negligible to render the method viable. In summary, the concept of superposition of separately tunable individual controllers requires their selectiveness in order to avoid overlapping of the neighbouring controllers. This requirement is well addressed in e.g. [31], but totally missed in e.g. [21], [22]. The latter proposes a set of separately tuned PR controllers to be superposed into a P-MR controller. The proportional part of the P-MR controller is calculated as the sum of proportional parts of the individual PR controllers. This does not make much sense either from the theoretical or practical standpoints and should not be regarded as a viable heuristics for such systems. This is because the Naslin polynomial method was not designed with the separation requirement in mind and the proportional gain is not kept negligible -on the contrary, please note that the proportional gain would increase proportionally with the increasing resonant pulsatance of consecutive resonant terms in the P-MR controller (13). Therefore, adding consecutive resonant terms would decrease stability margins till the moment when the critical gain is crossed, destabilizing the physical system. One should not fall here into a trap of a simplified simulation model with (5) as the plant transfer function. Such a simplified system has an infinite critical gain when controlled with a proportional controller. The closed loop system is stable for any positive k P . This is no longer true for higher order systems such as digitally controlled PWM converters with delays introduced for example by the modulator and the digital control system. Let us demonstrate that with the help of where L = 2 mH, R = 0.2 , τ controller = 0.1 ms and τ modulator = 1 2 τ controller . Designing individual PR controllers using (13) to reject disturbance at ω 1 = 2π · 50 · 6 rad s −1 , ω 2 = 2π · 50 · 12 rad s −1 and ω 3 = 2π · 50 · 18 rad s −1 would give k P 1 = 10.46, k P 2 = 21.13 and k P 3 = 31.79. The critical gain of the system is k critical = 60.90 as calculated using the margin() function [32]. A superposition by addition of proportional gains would destabilize the system, because k P 1 + k P 2 + k P 3 > k critical . Our conclusion is that the superposition of PR controllers proposed in [21] and [22] is by no means practical and should not be presented as a viable heuristics.
A correct heuristic for P-MR controller could involve tuning the PR controller using the Naslin technique and superposing it with resonant controllers (not proportional-resonant ones) tuned individually with a use of other techniques, for example by the trial and error method. This is because the required selectivity can be obtained for a PR controller and consecutive resonant controllers if they are configured to target different resonant frequencies, which is the case in any P-MR controller. To the best of authors' knowledge such an incremental design procedure is very popular among practitioners. Nevertheless, superposing individually designed resonant controllers makes it hard (if not impossible) to get a specified robustness. This is because even two ideal (not damped) oscillatory terms set for two different resonant frequencies still have overlapping amplitude spectra. This motivated us to propose a novel optimization procedure for a P-MR controller that produces systems with a guaranteed level of robustness assessed by means of the disk margin analysis.

VII. OPTIMIZATION WITH GUARANTEED ROBUSTNESS BY MEANS OF DISK MARGIN ANALYSIS
Nowadays, engineers of all fields have a variety of global optimization tools in their toolboxes. They help us to tackle problems that we do not know how to solve analytically, i.e. by means of manipulation of symbolic algebraic expressions. Such solutions should not be regarded as less elegant -their elegance is expressed in their practicality. Control engineering, especially controller tuning, is one of such fields where analytical solutions are known only for a relatively narrow set of problems. We are able to analytically derive gains for relatively narrow combinations of plant types, controller types and cost functions (i.e. performance indices). That is why in real-life situations we often encounter the trial and error method. It is worth noticing that the trial and error method does not necessarily mean that the trial part is totally randomized. It is an iterative process and the more guesses we make, the more educated next iterations of guessing can be. In fact, this is the main idea behind global optimization tools -they automatize this educated guessing procedure. An engineer can focus solely on problem definition, which in our case involves building a mathematical model of the physical system (here a grid connected converter), selecting its parameters to be optimized (here controller gains) and choosing a performance index to be minimized or maximized (here the objective is to draw sinusoidal currents). The type of optimisation tool is of secondary importance as long as it provides satisfactory exploration and exploitation capabilities for the problem at hand. As numerous tools can be used to solve a given optimization problem this choice often reflects preferences of a particular engineer. For example, we prefer particle swarm optimization (PSO) over other search algorithms, e.g. genetic algorithms or ant colony optimization. The PSO is a well-documented tool which is why we refrain from describing it. Comprehensive documentation can be found e.g. in MathWorks Help Center [33], whereas practically countless implementations can be downloaded at no charge from MathWorks File Exchange, from the very basic code [34] to the code run using parallel computing [35] and a wide variety of tutorials, including video ones, e.g. [36].

A. THE CURRENT CONTROL SYSTEM
The optimizer uses a mathematical model of the current control system to assess solutions. A schematic diagram of the modelled control system is shown in Fig. 4.
The main features of the simplified model used during the optimization are: • a linear discrete time model of the plant (the zero-order hold (ZOH) method); • six discrete time resonant controllers, three per d and q axis, (Tustin with pre-warping method) of PI-MR type with the following resonant pulsatances: ω 1 = 2ω, ω 2 = 6ω, ω 3 = 12ω; • v dc = const and ω = const; • ideal synchronisation with the grid voltage vector (no PLL); • a grid modelled using ideal voltage source; • no pulse width modulator, lack of current ripples. The goal of the simplified model is to accelerate the optimisation in terms of wall-clock time. Not modelling some aspects, e.g. PWMs, significantly speeds up tests without compromising the reliability of controller gain optimization. This is because even in the target system with PWMs the controllers work on signals averaged over the sampling period.

B. THE PLANT DESCRIPTION
For v dc = const and ω = const the plant state-space model in the dq rotating reference frame is as follows (in the continuous-time domain): where There are two state variables i d and i q collected in the state vector x, two control variables u d and u q collected in the control vector u and two disturbance variables v d and v q collected in the disturbance vector v. There are also three matrices: A -the state matrix, B -the control matrix and E -the disturbance matrix. The continuous-time model (31) is discretized by using a ZOH method implemented in the Matlab environment to produce its equivalent discrete representation taking into account the fact that the physical system performs ZOH action within the modulator. It should be stressed that the ZOH discretisation method is the only logical choice for this part of the system and this choice has nothing to do with the discretisation method chosen for the controller -it is perfectly fine to choose a different one for the controller itself to meet some specific design requirements such as resonant frequency matching by prewarping. The ZOH method, as well as the ZOH element in the modulator, introduce a half sample delay [6]. The obtained discrete-time model is where F, G, Z are the discrete-time system matrices.

C. THE PI-MR CURRENT CONTROLLER DESCRIPTION
The grid currents i d and i q are controlled in the dq rotating reference frame using two PI-MR controllers which can be expressed in the Laplace domain: where k P , k I , k 1 , k 2 and k 3 are the decision variables (gains to be tuned). In order to convert the controller model from the continuous to the discrete time domain, Tustin approximation with pre-warping is used. Therefore, as it has been shown in Fig. 4, H PI (z), H Rω 1 (z), H Rω 2 (z), H Rω 3 (z) are discrete time equivalents of G PI (s), G Rω 1 (s), G Rω 2 (s), G Rω 3 (s), respectively. The other key parameters are compatible with the ones given in Table 1.

D. THE COST FUNCTION
Taking into account that optimality does not imply robustness, the main effort that differentiates successful robust solutions from just optimal ones goes into designing a proper cost functional. The trade-off between improving reference tracking and disturbance rejection and keeping a required robustness to plant uncertainty is fundamental for any control system [37]. One of the common approaches to combining conflicting objectives into a single objective function is the weighted sum method with constant weights that have to be selected by the designer prior to the optimization. For the current controllers in the grid-tied converter the cost function could be where e d and e q are current control errors, u d and u q represent control signals (reference voltages for the converter), T s is the sampling period and t stop denotes the end of the interval over which the performance is assessed. The penalty factor β is to enable the designer to control for the robustness -the response slows down with an increasing penalty factor. Choosing a sufficiently high penalty for control signal dynamics would produce practical controller gains, which is demonstrated e.g. in [38]. However, the robustness of the resulting system is not tested until the optimisation is completed -stability margins are not explicitly included in the cost functional (34). The originality of the present work is that the stability margins are explicitly stipulated in the cost functional. First of all, defining a way to measure stability margins for a multiple-input multiple-output (MIMO) system, which is the case here, is not a trivial task. Limitations of classical margins, as well as the proposed solution to overcome those limitations by using disk margins, are presented in [39]. The stability analysis using disk margins was recently included in Robust Control Toolbox (Matlab) [40], significantly simplifying its deployment in engineering practice. In the disk margin analysis the size of the disk is the robustness measure -it represents how far from becoming unstable the system is. A disk-based uncertainty model can be applied at inputs and/or outputs of the plant. It was tested that in order to get practical tuning results for a grid-tied converter it is necessary to compute the stability margins when considering variations at both the plant inputs and plant outputs. This is because in a physical application, there is uncertainty at both the inputs and outputs of the plant, related to the control delay, parameter identification errors and measurement noise. In general, calculating a disk margin for an input and output independent concurrent uncertainty gives a more conservative result than taking into account just the one at the input. Also, the loopat-a-time disk margin (employed for a grid-tied converter e.g. in [41]) is less conservative than the multiloop disk margin. The latter is able to capture maximum tolerable uncertainties across all the feedback loops, which is a major advantage over the classical margin analysis. In the discussed converter the uncertainties may occur in all the loops simultaneously, therefore, the multiloop analysis is selected as a more practical one. In this study unbiased margins are used (the skew is set to zero for diskmargin() [42]). The designer specifies the required robustness by selecting the circular exclusion region centred near the critical point, which in turn guarantees that the open-loop response stays at a specified safe distance from the critical point at all frequencies. In this study the disk size of 0.35 is selected as a compromise between robustness on one hand and reference tracking and disturbance rejection capabilities on the other hand.
The goal for the optimizer is to find controller gains that minimize  Fig. 4.
The Greek letter α is selected to denote the disk size in order to keep it consistent with the topical literature -not to be confused with α present in the Naslin polynomial.
To get a solution with the guaranteed stability disk margin, the cost function for the swarm should strongly penalize solutions that cross α threshold . The penalty should be high enough to render cost function values for all the solution that VOLUME 10, 2022 cross α threshold greater than for any solution that obeys the constriction. This makes clear distinction between feasible (α candidate ≥ α threshold ) and infeasible (α candidate < α threshold ) solutions. The swarm itself is not limited to testing only feasible solutions, but the region of infeasible solutions is naturally left by the swarm as the search transitions from the exploration stage to the exploitation stage. A function of the form where α candidate is the disk size of the candidate solution and b has to be set to a number beyond reach of J SSE , creates a funnel-shaped area with all the feasible solutions at the bottom of the funnel. The funnel guides particles towards feasible solutions mainly during the exploration phase. At the exploitation stage particles operate predominantly at the bottom of the funnel and the resulting solution is found near the wall created by b. One might conclude that this cost function introduces one new parameter to be selected. However, the solution is practically not sensitive to it as far as it conforms to the following rule: b is greater than the biggest possible J SSE for any feasible solution. The condition can be met by assessing J SSE for any stable solution found by guessing and checking and then setting b several orders of magnitude greater to be on the safe side. A schematic diagram of the control system connected to the optimizer is depicted in Fig. 4.

VIII. RESULTS
PSO is run with the cost function (36) and the mathematical model of the converter system (shown in Fig. 4). It should be remembered that PSO is a stochastic search algorithm and it may happen that some runs of even well-posed and well-conditioned optimization rounds end up in a dead end. One should run any optimizer of that kind several times to assess the reproducibility of the final result. For well-posed and well-conditioned problems unsuccessful runs are outnumbered by the successful ones and this has been found to be the case for the stated problem.

A. PARTICLE SWARM OPTIMIZATION
Solutions proposed by particles in each iteration are rated in the control system stimulated by reference signals and disturbances shown in Fig. 5. These signals include: step changes in reference currents i ref d and i ref q ; at 0.1 s grid voltage becomes distorted (asymmetry of 3% and presence of 5 th , 7 th , 11 th and 13 th harmonics); from 0.4 s to 0.55 s symmetrical voltage sag of 10% during active power load; from 0.9 s to 1.05 s symmetrical voltage sag of 10% during reactive power load; at 1.3 s grid voltage is restored to sinusoidal. The optimization problem is well-conditioned which is demonstrated in Fig. 6 showing four consecutive consistent optimization runs.  Evolution of the swarm is depicted in Fig. 7. The corresponding evolution of the grid current shape is shown in Fig. 8.

B. SIMULATION VERIFICATION USING FULL MODEL
The main features of the full model used to verify the final solution before implementing it in a physical converter: • a full switching model (NPC three-level converter); • the control system implemented in C-code including six discrete time resonant controllers • a PI type v dc controller with anti-windup • PLL based on Delayed Signal Cancellation • a grid modelled using voltage source with internal impedance • PWM with zero-sequence signal. The controller settings found by the swarm using the simplified model are then fully verified in the numerical model that includes PWM, DSC-PLL, filters and dead-times. The steady state under distorted grid voltages is shown in Fig. 9.
The response to the step load change is demonstrated in Fig. 10. Loading the converter with reactive power is shown in Fig. 11. Reaction to voltage sags of 10% and 20% are presented in Figs 12 and 13, respectively. The converter behaves as expected.

C. EXPERIMENTAL VERIFICATION
This paper is primarily theoretical, i.e. its goal is to show fundamental flows in some attempts to use the Naslin polynomial method reported in the literature in regard to multiresonant controller tuning -the attempts wrongly reported as successful ones -and then to demonstrate the developed method that provides required robustness of the solution and does not abuse algebra. In order to demonstrate practicality of the proposed method an experimental verification is also  provided. A physical grid-tied converter photographed in Fig. 14 of the parameters assumed for the optimization stage (see Table 1) is used. It should be stressed that no additional fine tuning takes place after transferring all the gains into the physical system. This is one of the features of a robust tuning method -possible discrepancies between the mathematical model and the physical system do not render the dynamic response of the two systems significantly different. For comparison purposes the system is always first run with a PI type current controllers for each reference and disturbance test scenario. This is to justify the need of resonant controllers for a given converter if sinusoidal symmetrical currents are expected. The performance verification under nonsinusoidal grid voltage conditions takes place at THD of grid voltage equal to 11.6 %. The response to a 10 kW load step change is shown in Fig. 15. The resulting grid current at the steady state is visibly nonsinusoidal and has THD of 23.7 %. This is expected as the structure of the controller does not follow the internal model principle -it does not implement the  Then the structure of the controller is changed to the proportional-integral-multiresonant one and all its gains are set exactly as obtained during the optimization procedure. The distortion of the grid voltage in the experimental setup is introduced with the help of a weak grid and additional nonlinear loads. As a result, THD of the voltage depends also on load conditions and for the PI-MR controllers amounts to 15.4 %. In comparison with the case in which PI controlled converter is deployed, the voltage distortion increases, which is to be expected. This in fact results from better conditions introduced by the PI-MR controlled converter. The voltage   is now higher for the nonlinear loads due to lack of voltage drops from higher harmonics effectively suppressed on the controlled rectifier side. Additional nonlinear loads draw now  higher currents, distorting the voltage even deeper. But even for a higher distortion of the grid voltage the current distortion drops significantly in comparison to the PI controlled system and is equal to 9.1 %. The relevant oscillographs are depicted in Fig. 16. This steady-state behaviour is expected and rather obvious in view of the internal model principle. What is crucial here is that transient states are also of high quality and that the controller is robust to inevitable identification errors -there was no need to tinker with the gains after moving into the physical power electronic converter.
Similarly, satisfactory results are obtained both for transient and steady-state behaviour under unbalanced grid voltages. The asymmetry of voltages is 11.2 % and 16.7 %, respectively. The root cause of that difference is exactly the same as for the previous case scenario. The current THD drops from 17.3 % to 9.0 %. The relevant oscillographs are shown in Figs. 17 and 18. Again, this is not something surprising as we implemented an oscillatory term specifically aimed at asymmetry of the disturbance. What is significant is that the transients are of high performance system (see Fig. 18).

IX. DISCUSSION
Robustness of a control system is a fundamental requirement to make the system practical. Controller tuning methods should explicitly address this issue to be regarded as good candidates for actual real-life applications by the practitioners. Our results indicate that the disk-based stability margin multiloop analysis is an absolutely complete theoretical tool that is able to produce fully applicable controller settings. Moreover, the design is straightforward and the required expert knowledge to successfully deploy it is rather moderate. A single value representing the disk size has to be chosen by the user and the system with the requested robustness is found by means of particle swarm optimisation. The effectiveness of the method confirmed in the case of multiresonant controllers for a grid-tied converters should encourage future research within the context of other repetitive control systems in power electronics and drives, including true/pure sine wave inverters and repetitive motion control. The research could be potentially expanded to cover not only multiresonant controllers but also iterative learning controllers [38].

X. CONCLUSION
It has been demonstrated that several Naslin polynomial based tuning methods reported in the literature for multiresonant controllers are in fact flawed and should not be regarded as viable approaches. A novel tuning procedure has been proposed. The resulting control system for a grid-tied converter provides an explicitly stipulated stability margin. At no stage of the design procedure this margin is compromised and this is thanks to the cost function definition that rates all solutions with smaller margins as infeasible. Guaranteed robustness of the system as well as observed robustness of the optimization process itself manifested in reproducibility of found minima make the method a good candidate for industrial grade products. VOLUME 10, 2022