Zero-Phase FIR Filter Design Algorithm for Repetitive Controllers

Abstract

Repetitive controllers (RCs) are linear control structures based on the internal model principle. This control strategy is known for its ability to control periodic reference signals, even if these signals have many harmonic components. Despite being a solution that results in a good performance, several parameters of the repetitive controller need to be correctly tuned to guarantee its stability. Among these parameters, one that has high impact on the system performance and stability is the finite impulse response (FIR) filter, which is usually used to increase the stability domain of RC-based controllers. In this context, this paper presents a complete tutorial for designing the zero-phase FIR filter, which is often used to stabilize control systems that use RC-based controllers. In addition, this paper presents a Matlab^® application developed for performing the stability analysis of RC systems and designing its FIR filter. Simulation and experimental results of a shunt active power filter are used to validate the algorithm and the Matlab^® application.

1. Introduction

Repetitive controllers (RCs) are becoming more and more popular due to their theoretical ability to ensure zero steady-state error for periodical reference signals and periodical disturbances. During the last decade, several researchers have been focused on proposing fast, accurate, and stable RCs, applicable in different areas, such as robotics and power electronics [1,2,3,4,5,6,7,8].

The mathematical basis for evaluating the stability of a RC was first presented in [9], where the concept of stability domain was used. The RC structure proposed in [9] contains a periodic signal generator in parallel to a second direct path, whose gain (a) enables enlarging the stability domain of the RC system. For being stable, the RC system stability domain must contain a Nyquist trajectory determined by the controlled plant [9]. However, regardless of the value of a, the stability domain does not contain the point $(0,0)$ of the complex plane, which is part of the Nyquist path of most practical systems. For this reason, a gain attenuation block $Q\left(z\right)$ must be introduced in the RC’s structure to further enlarge its stability domain to include the point $(0,0)$. The most common solutions used in this scenario are a constant attenuation (such as $Q\left(z\right)=0.98$[10] or $Q\left(z\right)=0.95$ [11]) or a low-pass filter.

Many types of low-pass filters can be used to improve the stability of RC systems. As these filters are inserted in the periodic signal generator of the repetitive controller, any phase lag caused by them will negatively affect the performance of the controller [12], possibly detuning the frequencies to be controlled. Therefore, the main solutions of filters that can be used in RCs are either zero-phase or linear-phase low-pass filters. However, when using a linear-phase filter, a phase shift compensation must be carried out.

Since RC-based controllers are based on periodic signal generators, they are often implemented as digital controllers. In this scenario, the low-pass filter to be included in its structure must also be a digital filter. As indicated in [13], this digital filter can be either a finite impulse response (FIR) or an infinite impulse response (IIR) filter. When using a FIR filter, one can still choose between a non-causal FIR filter [13] (which presents a zero-phase shift characteristic) and a causal FIR filter [14] (which presents a linear-phase shift characteristic). On the other hand, when using an IIR filter, several different solutions can be found in the literature, such as Butterworth [15] and Chebyshev [12] low-pass IIR filters; however, the phase shifts caused by most of them are non-linear [14]. Alternative approaches using quasi-linear-phase IIR filters are also viable solutions [14], but they are not yet as widely used as FIR filters. Depending on the chosen structure, the frequency spectrum of the controller may change, affecting the main advantages of RCs in regulating specific harmonic frequency components [16], or the stability may be compromised, such as when using IIR filters [17]. Against this background, zero-phase FIR filters are considered the best choices for RC-based controllers, since they are BIBO stable and do not apply any phase shift to the controller [17].

The design of the FIR filter to include the origin of the complex plane in the stability domain of RC systems is essential. It not only impacts the system stability but also affects its control performance. The filter passband and order need to be carefully selected to effectively attenuate components above the maximum harmonic frequency to be regulated. Additionally, depending on the complexity of the RC structure, it should be broken down into primitive repetitive cells (PRCs) in parallel to determine where the FIR filter should be allocated [18]. To the best of the authors’ knowledge, no literature presents a detailed step-by-step process for designing FIR filters for RC systems or a tool that automates this project.

It is not the aim of this paper to propose a new RC scheme. However, given the different structures of repetitive controllers available and the impact of the FIR filter on the stability of the system and on the control performance, we present a tutorial and a complete algorithm for designing zero-phase FIR filters to be used in any RC scheme, ensuring the necessary stability margin with the best control performance achievable. Section 2 presents the fundamentals of RCs, the concept of the stability domain, and the decomposition of non-trivial RC structures into PRCs in parallel. In Section 3, the impacts of the FIR filter on the stability and performance of the RC are analyzed. Although the subjects of Section 2 and Section 3 are objects of several previous papers, they are summarized here for clarity. The paper’s main contribution, i.e., the proposed algorithm for the zero-phase FIR filter design, is described in Section 4. As a second contribution, a Matlab^® application for the automatic implementation of the proposed design procedure was developed, which is described in Section 5. In Section 6, the results of the experiments (which were performed to validate the design tool) are shown. Conclusions are presented in Section 7.

2. Fundamentals of Repetitive Controllers

In 1981, Inoue et al. [19] proposed a linear control strategy based on the internal model principle [20] that became known as the repetitive controller. The main characteristic of this control strategy is the ability to control periodic reference signals, regardless of their harmonic content. Although the repetitive controller was initially developed for robotics applications, it proved to be effective in power electronics applications, especially when controlling electrical variables that have high harmonic content, such as in active power filters [16,21,22].

As shown in Figure 1, the repetitive controller has a periodic signal generator in its structure. This feature enables the control system to obtain zero steady-state error for periodic reference signals. Seeking to expand the use of repetitive control, Hara et al. [9] added a second direct path (with gain a) in the controller, which is also shown in Figure 1. Gain a enables the control system designer to obtain better stability configurations for the repetitive controller. The block $z^{-N}$ implements a delay of N samples, where N is the number of samples per fundamental period of the reference signal.

Taking into account the repetitive controller presented in Figure 1 (with repetitive gain $K_{rc}$) and using it to control a plant $G_p\left(z\right)$, then, for an unit negative feedback system, the following conditions must be true in order to obtain an input–output stable system (${\ell }_2$-stability analysis) [9]:

• $G_m\left(z\right){[1+a·G_m\left(z\right)]}^{-1}$ must be a proper stable rational transfer function, in which $G_m\left(z\right)=K_{rc}·G_p\left(z\right)$; and
• ${∥\frac{1+(a-1)·G_m\left(z\right)}{1+a·G_m\left(z\right)}∥}_{{\mathcal{H}}_{\infty }}<1$ must hold.

Hara et al. [9] also show, from developing Condition 2, that the Nyquist contour of $G_m\left(z\right)$ must be contained in an area of the complex plane known as the stability domain of the repetitive control system. These stability domains depend on parameter a, which may be different for each application. As parameter a increases, the stability domain is enlarged, as presented in Figure 2.

To show the importance of evaluating the stability domain of a RC-based system, one can consider a unit negative feedback system with a hypothetical plant given by

$$G_{p1}\left(z\right)=\frac{z-0.94}{z-0.975},$$

(1)

which is controlled by a repetitive controller with $K_{rc}=1$. The Nyquist contour of $G_{m1}\left(z\right)=K_{rc}·G_{p1}\left(z\right)$ and the stability domain of the RC-based system for $a=0$ are presented in Figure 3a (Case 1), and for $a=0.5$ in Figure 4a (Case 2). This hypothetical plant was chosen because its Nyquist contour is entirely contained in the right half-plane of the complex plane, allowing the RC-based system to be stable without having to add a gain attenuation block (Q(z)) in the periodic signal generator of the repetitive controller.

As can be seen in Figure 3a, the Nyquist contour of $G_m\left(z\right)$ is not contained in the stability domain of the RC-based system for Case 1, leading to an unstable system, as shown in Figure 3b for a sinusoidal input of frequency $f_1=50\mathrm{Hz}$. On the other hand, the same does not happen for Case 2 (Figure 4b) since the Nyquist contour of $G_m\left(z\right)$ is entirely contained in the new stability domain of the RC-based system (Figure 4a).

However, the stability domains presented in Figure 2 do not contemplate the point $(0,0)$ of the complex plane. As a consequence, the conventional repetitive controller cannot be used for controlling strictly proper plants, which have more poles than zeros and a Nyquist contour converging to the point $(0,0)$. In this sense, it is important to evaluate alternatives that allow enlarging the stability domain of RC-based systems.

2.1. Enlarging the Stability Domain of a Repetitive Control System

As mentioned before, most practical control systems have plants with nonzero relative degrees, i.e., plants with strictly proper transfer functions. For these systems, the Nyquist contour of the plant $G_p\left(s\right)$ (and, hence, $G_m\left(z\right)$) tends toward point $(0,0)$ of the complex plane as the frequency increases. Thus, it is necessary to enlarge the stability domain of the RC system to meet the stability requirements.

In order to do so, a block $Q\left(z\right)$ is usually added in series with the delay of the repetitive controller, in the periodic signal generator. Figure 5 shows a discrete-time version of the repetitive controller proposed in [9], with $Q\left(z\right)$. As presented in [23], when adding the block $Q\left(z\right)$, the condition used to obtain the stability domain becomes

$${∥Q\left(z\right)·\frac{1+(a-1)·G_m\left(z\right)}{1+a·G_m\left(z\right)}∥}_{{\mathcal{H}}_{\infty }}<1.$$

(2)

Equation (2) can be developed to obtain a new stability domain equation, which will now also depend on $Q\left(z\right)$, as demonstrated in [23]. Therefore, (2) can be evaluated as:

$$\underset{\left|\omega \right|<\frac{{\omega }_s}{2}}{sup}\left|Q\left(e^{j\omega T_s}\right)·\frac{1+a·G_m\left(e^{j\omega T_s}\right)-G_m\left(e^{j\omega T_s}\right)}{1+a·G_m\left(e^{j\omega T_s}\right)}\right|<1.$$

(3)

Then, (3) can be rewritten as:

$$\underset{\left|\omega \right|<\frac{{\omega }_s}{2}}{sup}\sqrt{\frac{{\left[A·(1+aX-X)-B·(aY-Y)\right]}^2+{\left[A·(aY-Y)+B·(1+aX-X)\right]}^2}{{(1+aX)}^2+{\left(aY\right)}^2}}<1.$$

(4)

Note that, for simplification, the notations $X=Re\left\{G_m\left(e^{j\omega T_s}\right)\right\}$, $Y=Im\left\{G_m\left(e^{j\omega T_s}\right)\right\}$, $A=Re\left\{Q\left(e^{j\omega T_s}\right)\right\}$ and $B=Im\left\{Q\left(e^{j\omega T_s}\right)\right\}$ are adopted.

As (4) represents the stability boundary condition, the stability domain can be obtained considering the entire area defined by this boundary curve. In order to do so, one can suppress the supremum of (4). Consequently, the stability domain is determined through the inequality:

$${\left[A·(1+aX-X)-B·(aY-Y)\right]}^2+{\left[A·(aY-Y)+B·(1+aX-X)\right]}^2<{(1+aX)}^2+{\left(aY\right)}^2.$$

(5)

The inequality shown in (5) can still be simplified in a region that is delimited by a circle equation. By doing so, the stability domain of a RC system with constant $a\in \mathbb{R}$ can be calculated by (for $|Q\left(e^{j\omega T_s}\right){|}^2=A^2+B^2$):

$$X^2·f_1(a,|Q\left(e^{j\omega T_s}\right)\left|\right)+Y^2·f_1(a,|Q\left(e^{j\omega T_s}\right)\left|\right)+X·f_2(a,|Q\left(e^{j\omega T_s}\right)\left|\right)<\left[1-|Q\left(e^{j\omega T_s}\right){|}^2\right],$$

(6)

in which

$$\begin{array}{cc}\hfill & f_1(a,|Q\left(e^{j\omega T_s}\right)\left|\right)=a^2|Q\left(e^{j\omega T_s}\right){|}^2-2a|Q\left(e^{j\omega T_s}\right){|}^2-a^2+{\left|Q\left(e^{j\omega T_s}\right)\right|}^2,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & f_2(a,|Q\left(e^{j\omega T_s}\right)\left|\right)=2a|Q\left(e^{j\omega T_s}\right){|}^2-2{\left|Q\left(e^{j\omega T_s}\right)\right|}^2-2a.\hfill \end{array}$$

(8)

This means that (6) must hold for all $(X,Y)=(Re\left\{G_m\left(e^{j\omega T_s}\right)\right\},Im\left\{G_m\left(e^{j\omega T_s}\right)\right\})$, in which $\omega \in [-\frac{{\omega }_s}{2},\frac{{\omega }_s}{2}]$ (${\omega }_s$ is the sampling frequency).

Although the block $Q\left(z\right)$ can be used to increase the stability domain of the RC system, making it contain the point $(0,0)$ of the complex plane, the block $Q\left(z\right)$ may also cause attenuation and phase shifts in the signal synthesized by the periodic signal generator, making it no longer periodic. As a consequence, this phase shift can disrupt the proper functioning of the repetitive controller. Thus, $Q\left(z\right)$ is generally recommended to be a constant attenuation (i.e., a constant gain lower than 1) or a low-pass zero-phase filter with a finite impulse response (FIR) [9]. The main characteristics of these solutions are described below.

• $Q\left(z\right)$ as a constant attenuation: This solution promotes a constant reduction in the amplitude of the signal that is being fed back by the periodic signal generator of the repetitive controller, compromising its operation in relation to the internal model principle [20]. As a consequence, despite improving the system stability, it does not have a zero steady-state error. In order for the influence of $Q\left(z\right)$ to be small on the stationary error, it is usually chosen to be as close to one as possible. According to (6), the stability domain of the RC system is increased for $Q\left(z\right)<1$ (Figure 6). Examples of $Q\left(z\right)$ as a constant attenuation can be found in [10] (where $Q\left(z\right)=0.98$) and [11] (where $Q\left(z\right)=0.95$).
• $Q\left(z\right)$ as a FIR low-pass filter: This solution promotes the use of a low-pass filter in the periodic signal generator of the repetitive controller. Thus, it is expected that the dynamics of the control system will not be significantly altered for low frequencies, still making it possible to obtain a zero steady-state error (or as close to zero as possible). On the other hand, the repetitive controller will have a pass band defined by $Q\left(z\right)$, thus, it will no longer control exogenous signals whose harmonic frequencies are beyond the cutoff frequency of the low-pass filter $Q\left(z\right)$. The most common solution applied in the literature is $Q\left(z\right)=0.25·z+0.5+0.25·z^{-1}$ (some examples are [24,25,26]).

Considering $Q\left(z\right)$ as a zero-phase FIR filter, its transfer function will be given by:

$$Q\left(z\right)=a_0·z^{\frac{L}{2}}+a_1·z^{\frac{L}{2}-1}+\cdots +a_{\frac{L}{2}}+\cdots +a_L·z^{-\frac{L}{2}},$$

(9)

in which $a_0=a_L$, $a_1=a_{L-1}$ and so on. However, the FIR filter presented in (9) is not implementable. In order to solve this problem, one can see that (9) can be rewritten as

$$Q\left(z\right)=z^{\frac{L}{2}}·\left(a_0+a_1·z^{-1}+\cdots +a_L·z^{-L}\right),$$

(10)

and since $z^{-N}$ is cascaded in series with $Q\left(z\right)$, one can implement them as

$$z^{-N}·Q\left(z\right)=\underset{➀}{\underbrace{z^{-N+\frac{L}{2}}}}·\underset{➁}{\underbrace{\left(a_0+a_1·z^{-1}+\cdots +a_L·z^{-L}\right)}},$$

(11)

where ➀ and ➁ are implementable for $L/2<N$. Note that L must be even.

To apply $Q\left(z\right)$ to more elaborate repetitive control structures (here referred to as non-trivial RC-based controllers), one can use the unified approach to evaluate the repetitive controllers proposed in [18]. A brief summary of this analysis is presented in the following subsection.

2.2. Unified Approach to Evaluation of Repetitive Controllers

Non-trivial RC-based controllers, such as $nk\pm m$ RCs [4,27,28], may also contain FIR filters to improve system stability. A simple way to understand where to place these filters is achieved by decomposing the RC-based structure into an association of primitive repetitive cells (PRCs) in parallel and, then, placing a FIR filter on each PRC [18].

The PRC is a complex control structure (i.e., it is a block with complex signals as input and output) that allows the system to control only a selected set of harmonic components, here identified as set $H=\{nk+m$|$k\in \mathbb{Z}\}$. It is important to note that, for $n\ne 1$, the PRC does not control all harmonics, as the conventional repetitive controller does (Figure 7). However, multiple PRCs can be parallelized to behave as the conventional RC [29].

In this paper, the PRC presented by Neto et al. [18] is used to evaluate where the control system designer must place the block $Q\left(z\right)$ to improve the stability characteristics of non-trivial RC-based systems. This PRC is illustrated in Figure 8, which is composed by the following blocks:

• A generic delay of $N/n$ samples ($z^{-N/n}$): The PRC has a generic delay of $N/n$ samples in its structure, in which N is the number of samples per fundamental period. The parameter n can be used to select the periodicity of the set of harmonic components that the PRC applies high gain (Figure 7). This means that, disregarding the other blocks presented below, the adequate selection of n allows the RC-based system to control periodic signals whose harmonics are in the set $H=\{nk$|$k\in \mathbb{Z}\}$ instead of all harmonic components;
• A complex gain $e^{j2\pi \frac{m}{n}}$: If a complex gain $e^{j2\pi \frac{m}{n}}$ is cascaded with the generic delay $z^{-N/n}$ presented above, then the frequency response changes suffer frequency shifts. Therefore, the adequate selections of m and n allow for the RC-based system to control periodic signals whose harmonics are in the set $H=\{nk+m$|$k\in \mathbb{Z}\}$ instead of $H=\{nk$|$k\in \mathbb{Z}\}$, where $m\in \mathbb{N}$ is restricted to $0\le m\le (n-1)$.
• A constant gain a: As done for the conventional repetitive controller [9], the PRC has a second direct path with the constant gain a. Gain a establishes a constant proportion between the repetitive action and a proportional action, and, as a consequence, it can be used to enlarge the stability domain of the control system [23]. Changes in this parameter do not make the stability domain contain the point ($0,0$); thus, changes on this are not enough to make the repetitive controller applicable to strictly proper plants.
• A low-pass FIR filter $Q\left(z\right)$: As described in Section 2.1, this block is usually a zero-phase FIR filter used to enlarge the stability domain of the control system.

Since the target application defines the exogenous signals of the control system (and, as a consequence, the harmonic components they have), it also defines the parameters m and n. For instance, when evaluating the harmonic content of the currents demanded by a three-phase diode rectifier, one can see that they have all harmonics of the set $H=\{1,5,7,11,13,\cdots \}$. Thus, for this application, a PRC that controls the harmonics that belong to the set $H_1=\{6k+1$|$k\in \mathbb{Z}\}$ can be used in parallel with a PRC designed to control the harmonic components of $H_2=\{6k-1$|$k\in \mathbb{Z}\}$ as unique controllers, as conducted in [4].

As presented in [23], the ${\ell }_2$-stability analysis of a control system with PRC can be done using the same equations presented in [9] for the conventional repetitive controller. Thus, Equation (6) is for the evaluation of the stability domain.

3. Impact of the FIR Filter on the Stability and Performance of RC-Based Control Systems

As presented in [23], the bandwidth of a PRC can be evaluated from the ${\mathcal{L}}_2$-stability analysis (or the ${\ell }_2$-stability analysis in digital implementation) of the control system. This means that, from this stability analysis, the control system designer can estimate up to which harmonic a PRC can compensate, which depends on the cutoff frequency of the FIR filter used to enlarge the stability domain.

On the other hand, as shown by Neto et al. [18], RC-based controllers can be decomposed in multiple PRCs in parallel (Figure 9). Thus, considering a linear control system, if one of the loops becomes unstable due to its PRC, the entire control system also becomes unstable. Therefore, instead of evaluating the entire RC-based controller, which may be mathematically complex, the analysis can be simplified by evaluating multiple closed-loop systems that use a single PRC as a controller (as conducted in [23]).

As presented in Section 2, the stability domain of a PRC must contain the Nyquist contour of $G_m\left(z\right)=K_{rc}·G_p\left(z\right)$, otherwise, the system becomes ${\ell }_2$-unstable [23]. However, most plants are modeled by transfer functions that have a relative degree greater than zero. Thus, for these plants, the Nyquist contour of $G_m\left(z\right)$ contains the point ($0,0$), making it necessary to include $Q\left(z\right)$ in the structure of the repetitive controller to enlarge its stability domain.

To explain the impact of the FIR filter on the stability and performance of RC-based control systems, a hypothetical system with discrete PRC is considered below (initially, without $Q\left(z\right)$, i.e., $Q\left(z\right)=1$). Figure 10a shows the stability domain of a PRC with parameters $K_{rc}=2$, $a=0.5$, $m=0$, and $n=1$ (which leads to the conventional RC [9]) and the Nyquist contour of $G_{m2}\left(z\right)$ for (ZOH discretization with $T_s=50$$\mathsf{\mu }$s):

$$\begin{array}{c}\hfill G_{p2}\left(s\right)=0.8·\frac{1.21\times {10}^7}{s^2+3000s+1.21\times {10}^7}\to G_{p2}\left(z\right)=\frac{0.01149z+0.01093}{z^2-1.833z+0.8607}.\end{array}$$

(12)

As a consequence of the Nyquist contour of $G_{m2}\left(z\right)$ not being contained in the stability domain, the system becomes unstable (Figure 10b). Therefore, a FIR filter $Q\left(z\right)$ should be used to enlarge the stability domain.

When considering a system with a PRC that has a FIR low-pass filter $Q\left(z\right)$ in its structure, the cutoff frequency of an ideal filter to this application ($f_c$) is limited to the highest frequency that makes the Nyquist contour of $G_m\left(z\right)$ remain in the stability domain. For instance, it can be seen in Figure 11a that the Nyquist contour of $G_{m2}\left(z\right)$ reaches the stability domain boundary at 530 Hz ($|G_{m2}\left(e^{j2\pi 530}\right)|$). This means that 530 Hz is the maximum frequency at which the system can operate, ${\ell }_2$-stable, without changes in the gain of the periodic signal generator.

The stability domain of the repetitive control system must be increased in order to cover the rest of the Nyquist contour of $G_{m2}\left(z\right)$. For this hypothetical system, the typical solution applied in the literature

$$Q\left(z\right)=0.25·z+0.5+0.25·z^{-1}$$

(13)

is not enough to achieve system stability (Figure 11b). This happens because the typical solution is not optimized for all situations. After using the algorithm proposed in this paper, it is possible to achieve a custom $Q\left(z\right)$ that provides system stability and the zero steady-state error (Figure 11c).

If the repetitive gain $K_{rc}$ is increased, then the magnitude of $G_m\left(z\right)$ also increases. As a result, the Nyquist contour of $G_m\left(z\right)$ reaches the stability domain boundary at a lower frequency, leading to a lower $f_c$ for the custom FIR low-pass filter that should be designed to make the system stable (assuming that other parameters of the repetitive controller, such as a, are not changed). Thus, when evaluating the control system, one can observe that a reduction in $f_c$ impacts the system performance in two ways:

• The frequency response of the controller will show a magnitude reduction for the frequency components beyond $f_c$), which means that the controller will reduce its ability to compensate for these frequency components (i.e., its bandwidth will be reduced);
• Since FIR low-pass filters are not ideal, a small reduction in magnitude is expected for frequency components that are below but near the cutoff frequency $f_c$. This reduction in magnitude means that the internal model principle [20] is not fully met, resulting in the controller not leading to zero steady-state error (Figure 11d). However, it is expected that the higher the $f_c$, the lower the steady-state error should be.

Therefore, after selecting the parameters of the repetitive controller, the cutoff frequency of the FIR low-pass filter should increase up to the limit allowed by the system stability domain to optimize the system performance.

4. Zero-Phase FIR Filter Design Algorithm

Based on the impact of the FIR filter $Q\left(z\right)$ on the stability and performance of RC-based control systems (Section 3), and knowing how to obtain the stability domain equation of these RC systems (Section 2), this section presents a proposed algorithm that obtains the best FIR filter $Q\left(z\right)$ to be applied in a RC system (Section 4.1). Regarding the position of the FIR filter, Neto et al. (2021) [18] show where to position the FIR filter in non-trivial RC-based controllers. This methodology for positioning the FIR filter is reviewed in Section 4.2.

4.1. Obtaining the Specifications of a FIR Filter for a Repetitive Controller Using Its Stability Domain

In order to obtain the specification of an optimized FIR filter $Q\left(z\right)$ for a RC controller, after selecting a RC structure, the control system designer should follow the steps below:

(a)

STEP 1: Firstly, the control system designer must evaluate the exogenous signals of the control system and determine what are the harmonic components that must be controlled. Using this information, one can determine the passband of the repetitive controller. In this paper, it is considered that the RCs must work on exogenous signals with harmonic content up to $f_{max}$.

(b)

STEP 2: Then, the control system designer must obtain the transfer function of the plant ($G_p\left(s\right)$ or $G_p\left(z\right)$). According to [23], both continuous and discrete approaches leads to inequalities that result in the same stability domains, thus, any of them can be used.

(c)

STEP 3: As shown in Section 2, non-trivial RC-based controllers can be decomposed into PRCs in parallel (Figure 9). Thus, one must decompose the selected RC structure into PRCs in parallel, which can be done following the guidelines presented in [18]. Then, based on this decomposition, one must take note of parameter a obtained for these PRCs. This parameter remains constant during the execution of the algorithm and it will be used later to plot the stability domain of the RC system.

(d)

STEP 4:$K_{rc}=1$ and $Q\left(z\right)=1$ are initially set. Then, the Nyquist contour of $G_m\left(z\right)=K_{rc}·G_p\left(z\right)$ (or $G_m\left(s\right)=K_{rc}·G_p\left(s\right)$) together with the stability domain (which should be calculated for parameter a obtained in the previous step) must be plotted. Examples of this step are presented in Figure 12a,b.

(e)

STEP 5: As shown in Section 3, the bandwidth of the repetitive controller can be obtained by evaluating at which frequency the Nyquist contour of $G_m\left(z\right)=K_{rc}·G_p\left(z\right)$ (or $G_m\left(s\right)=K_{rc}·G_p\left(s\right)$) reaches the stability domain boundary. Thus, the repetitive gain $K_{rc}$ can be tuned to select the highest frequency (here called as $f_c$) that, for $Q\left(z\right)=1$, will be contained in the stability domain. Phase-margin and gain-margin can also be used to select the repetitive gain. Note that $f_c$ must be greater than $f_{max}$. An example of this step is presented in Figure 12c;

(f)

STEP 6: As the magnitude of $Q\left(z\right)$ decreases, the maximum frequency of the Nyquist contour that is kept inside the stability domain increases (e.g., Figure 13). Based on this characteristic, the control system designer must gradually decrease the magnitude of $Q\left(z\right)$ in regular steps (here referred to as $\Delta q$), making it work as a constant attenuation ($Q\left(z\right)<1$), while evaluates the highest frequency on the Nyquist contour of $G_m\left(z\right)=K_{rc}·G_p\left(z\right)$ (or $G_m\left(s\right)=K_{rc}·G_p\left(s\right)$) that still is inside the stability domain. By doing this, it becomes possible to determine the superior magnitude limit of $Q\left(z\right)$ for the entire frequency spectrum. An example of this step is shown in Figure 12d. This step is further detailed in the flowchart presented in Figure 14.

(g)

STEP 7: With the curve plotted in STEP 6 (Figure 12d), one can estimate the order of the FIR filter (M) from the decay observed after the frequency $f_c$. Firstly, one should calculate the line that tangents the decay (line in red in Figure 12e). From this line, one can obtain the parameters $A_{dB}$ and $\Delta f$ (Figure 12e). Then, adapting from the Harris method [30], the order of the FIR filter can be estimated as

$$M\approx ceil\left[\frac{f_s}{\Delta f}\left(\frac{A_{dB}}{22}\right)\right]+2,$$

(14)

in case $ceil\left[\frac{f_s}{\Delta f}\left(\frac{A_{dB}}{22}\right)\right]$ is even, or

$$M\approx ceil\left[\frac{f_s}{\Delta f}\left(\frac{A_{dB}}{22}\right)+1\right]+2,$$

(15)

in case $ceil\left[\frac{f_s}{\Delta f}\left(\frac{A_{dB}}{22}\right)\right]$ results in an odd number. In (14) and (15), $A_{dB}$ is the desired stop-band attenuation in dB; $\Delta f$ is the transition bandwidth desired to the filter; and $f_s$ is the sampling frequency. It must be noted that M must be even to maintain the FIR filter symmetry. As can be seen in (14) and (15), the filter order is increased by “+2” so that the decay slope is steeper around the filter cutoff frequency, making the system more robust and stable. This step is further detailed in the flowchart presented in Figure 15.

(h)

STEP 8: The specifications of the FIR filter $Q\left(z\right)$ are its order M and its cutoff frequency. In order to improve system performance, from this point on, the FIR filter cutoff frequency is considered as the frequency at which the superior magnitude limit curve of $Q\left(z\right)$ crosses $-3$ dB ($f_{3dB}$). The frequency $f_{3dB}$ is shown in Figure 12e. In case the superior magnitude limit curve does not cross $-3$ dB, the frequency $f_c$ can be used as the cutoff frequency for the FIR filter. These data can be used to design the FIR filter using numeric computing software, such as Matlab. One can compare the magnitude response of the designed FIR filter with the curve plotted in STEP 6. An example of this comparison is presented in Figure 12f.

STEP 6 of the algorithm can be implemented as presented in the flowchart of Figure 14. The proposed flowchart for STEP 6 scans the Nyquist contour of $G_m\left(s\right)$ (or $G_m\left(z\right)$) checking up at which frequency the inequality that defines the stability domain is no longer satisfied. When finding a frequency that does not satisfy the inequality (which occurs when the Nyquist contour is outside the stability domain), the magnitude limit for $Q\left(z\right)$ is decreased in constant steps of $\Delta q$ (as conducted in the blue curve in Figure 12d in the frequency $f_c$). The inequality of the stability domain (presented in (6)) is updated with the new $Q\left(z\right)$, which increases the stability domain. Then, the algorithm continues scanning the Nyquist contour of $G_m\left(s\right)$ (or $G_m\left(z\right)$) and performing the same process each time the inequality that defines the stability domain is not satisfied.

STEP 7 of the algorithm can be implemented as presented in Figure 15. This flowchart has the vectors $f[·]$ and $q[·]$ (each one with k elements) as inputs and, based on them, it calculates the order of the FIR filter. This is done as follows:

(i)

Firstly, $A_{dB}$ is calculated from the difference between $q\left[k\right]$ and $q\left[1\right]$. One must convert the subtraction result to a magnitude in dB.

(ii)

Then, one must obtain the equation of the line that passes through the points ($f_c,1$) and ($f\left[k\right],q\left[k\right]$). This equation is referred to here as $l\left(f\right[·\left]\right)$.

(iii)

To obtain the line tangent to the decay of the magnitude limit curve of $Q\left(z\right)$ (line in red in Figure 12e), the flowchart calculates the error $e\left[j\right]=q\left[j\right]-l\left(f\right[j\left]\right)$, with $j\in \{k_1,\cdots ,k\}$ (where $f\left[k_1\right]=f_c$). If the vector $e[·]$ has any negative element, the line passes above the superior magnitude limit curve and, as a consequence, it is not yet the desired tangent line. Thus, the slope of the “tangent line” must be reduced, which can be done by changing the point ($f\left[k\right],q\left[k\right]$) to ($f[k-1],q[k-1]$) in the item $\left(ii\right)$ and recalculating the line equation. This procedure is repeated until the vector $e[·]$ has no negative elements.

(iv)

The parameter $\Delta f$ is calculated from the difference between $f_c$ and $f_{stopband}$, where $l\left(f_{stopband}\right)$ is equal to $q\left[k\right]$.

(v)

Finally, the parameter M can be calculated from $A_{3dB}$ and $\Delta f$ using (14).

The Matlab function “fir1()” uses a Hamming window to make a linear-phase causal filter. This function can be used in STEP 8, however, it becomes necessary to compensate for its linear-phase shift. This compensation can be done by reducing the order of the delay present in the periodic signal generator by $M/2$, as shown in (11).

From the presented methodology, the control system designer obtains an operational value for the parameter $K_{rc}$, in addition to obtaining the specifications of a custom FIR filter $Q\left(z\right)$ optimized for the evaluated RC system. However, relative stability can still be used to fine-tune the gain $K_{rc}$ [31].

4.2. Positioning the Zero-Phase FIR Filters in Non-Trivial RC-Based Controllers

As discussed in Section 3, non-trivial RC-based controllers, such as $nk\pm m$ RC [4,27,28], may contain FIR filters to improve system stability characteristics. A simple way to know where to place these filters is to decompose them into multiple PRCs in parallel, and, then, place a FIR filter on each PRC in the structure [18].

To demonstrate this methodology, the $nk\pm m$ RC proposed by Lu et al. [27] is considered below (Figure 16a). As presented in [18], this non-trivial repetitive controller can be decomposed into two PRCs with $a=0.5$ in parallel (Figure 16b). Then, using the block diagram presented in Figure 8 as reference, one must place the FIR filter $Q\left(z\right)$ on each PRC (Figure 17a). Furthermore, to implement the FIR filter, it is necessary to decrease the number of delayed samples of the periodic signal generator of the PRC in $M/2$, where M is the order of the FIR filter (as shown in (11)). Finally, one can recompose the original non-trivial RC-based structure (Figure 17b).

5. Development of a Matlab app for Automatic Implementation of the Proposed Algorithm

In this section, we present a Matlab application designed to support the design of control systems that utilize repetitive controllers. This Matlab app has three main screens: (i) the first screen, referred to as the “start screen”, allows the user to load the mathematical model of the plant into the app (only after that it is possible to proceed to the other two screens); (ii) the second screen, here referred to as the “stability acreen”, allows the user to validate the stability domain equations of a RC system, enabling visualization of the influence of the controller’s parameters on system stability; and (iii) the third screen, here referred to as the “FIR Filter Screen”, presents a tool for designing a FIR filter support tool (parameter $Q\left(z\right)$ of the repetitive controller). More details of each screen are given below. This Matlab app, developed by the authors of this paper, is available for download at https://www.mathworks.com/matlabcentral/fileexchange/74759-repetitive-controller-designer (accessed on 2 January 2023). In order to understand how to install the developed Matlab^® application, one can see the video and documentation in [32] and [33], respectively.

To make this Matlab app as general and useful as possible, we used the generic PRC equations to validate the stability domain. These equations can be manipulated to obtain the same equation shown in (6), (7), and (8), as described in [23]. Therefore, when using non-trivial RC-based controllers, the control system designer must first decompose them into PRCs in parallel (as shown in [18]) before utilizing the Matlab app presented in this section.

5.1. Start Screen

The start screen of the Matlab app is shown in Figure 18. In Figure 18a, the block diagram of a generic PRC is illustrated. The first step for filling this screen is to define the repetitive gain $K_{rc}$ of the control system in Figure 18b and the control system plant. For the plant, two approaches are possible: manually type the parameters of the plant in Figure 18c; or choose a variable already defined in the Matlab workspace, in Figure 18d. After filling these fields, the user must press the OK button (Figure 18e). Then, the two buttons present in Figure 18f are enabled, allowing the user to choose between evaluating the stability domain of the RC system (“system stability” button) or designing the FIR filter for the RC system (“Filter Design” button). Pressing the button on Figure 18g opens up a popup screen that contains general information about the Matlab app and a small summary of its operation.

To make the Matlab app more robust, we implemented a validation data routine to avoid possible typos and errors when filling in the fields on this screen. In cases where the user manually types the parameters of the plant in Figure 18c, the following tests are performed:

• It is verified if the repetitive gain is a real number of dimension 1 × 1;
• It is verified if the numerator and denominator fields are filled with row vectors; and
• When working with a discrete plant, it is verified if the sampling period is a real number of dimension 1 × 1.

When the user chooses a variable from the Matlab workspace, the following tests are performed:

• It is verified if the repetitive gain is a real number of dimension 1 × 1; and
• It is verified if the chosen variable exists and if it is one of the following types: tf (transfer function), zpk (zero-pole-gain), or ss (state space).

The user is only allowed to continue with the application if all conditions are satisfied. Otherwise, an error message appears that makes it impossible to continue until the error is corrected.

5.2. Stability Screen

After defining the plant, it becomes possible to evaluate the system stability as a function of the controller parameters. To do so, the user must press the “system stability” button in the start screen, leading to the stability screen. The stability screen of the Matlab app is shown in Figure 19.

As the reader already know, the PRC parameters are $K_{rc}$, m, n, N, a, and $Q\left(z\right)$ (Figure 8). At this point, the first parameter, $K_{rc}$, was already computed in the Matlab app (it was conducted in the start screen). The following three parameters are constants used to select the harmonic components to be controlled. These parameters are freely chosen by the user and they do not influence the system stability. Therefore, Parameters a and $Q\left(z\right)$ remain, which directly interfere with the stability domain of the RC system, as shown in (6), (7), and (8).

As shown in Figure 19, the user is provided with two options to observe the influence of Parameters a and $Q\left(z\right)$ (Figure 19a). However, $Q\left(z\right)$ is evaluated here as a constant attenuation. The user can either vary Parameters a and $Q\left(z\right)$ using sliders or freely define real values. The impact of these variations can be immediately observed in the graph shown in Figure 19c. In addition to the stability domain in green, the Nyquist contour of $G_m\left(s\right)$ (or $G_m\left(z\right)$) is also displayed in gray and red. The small red circles in the graph indicate the threshold frequency up to which the Nyquist contour belongs to the stability domain of the system (part of the Nyquist contour in red). This frequency is shown in Figure 19b, along with an indication of whether the system is stable or not with such a configuration.

In case any configuration is useful to the user, it is possible to capture such a graph and choose the location where one wants to save the generated image (Figure 19d). Some features of Matlab itself are also available to analyze the generated graph, such as zoom in, zoom out, data tips, and pan options (Figure 19e). Still, by pressing the “configurations” button in the tool bar (Figure 19f), it is possible to adjust some parameters of the sliders (minimum step, maximum value, and minimum value), some parameters of the text boxes are available for freely defining the values (resolution for rounding) and configuring the minimum and maximum limits of the graph axes.

To ensure no unexpected behavior, some error handling has also been implemented, as conducted for the start screen.

5.3. FIR Filter Screen

After defining the plant, it is also possible to use the Matlab app to support the design of the FIR filter $Q\left(z\right)$; both algorithms are presented in Figure 14 and Figure 15. The FIR filter screen of the Matlab app is shown in Figure 20.

By using this functionality of the created Matlab app, the control system designer can obtain the magnitude limit of $Q\left(z\right)$ for the entire frequency spectrum (such as in Figure 12d). However, to execute the algorithms, the user must fill in the available fields with the following information (Figure 20a):

• The magnitude of $Q\left(z\right)$ for low frequencies, called the “max initial value” and the desired decrement for parameter $Q\left(z\right)$. The authors recommend keeping the “max initial value” as one to make the Matlab app evaluate all possible values of $Q\left(z\right)$);
• The initial and final values for the frequency evaluation and the number of frequencies to be evaluated between these extremes;
• Parameter a of the evaluated PRC;
• With respect to the graph that will be plotted, one must choose the Y-axis display mode, which can be “absolute” or “decibels” modes; and
• If the user defines the plant as continuous in the start screen, the sampling frequency of the discrete controller will also be required.

The curve obtained from the algorithm proposed in this paper is presented in the FIR filter screen (in Figure 20b) and indicates the magnitude limit for designing $Q\left(z\right)$. In Figure 20c, the order and the cut-off frequency of an implementable FIR filter that meets the requirements is shown. Finally, by clicking on the button “export data to workspace...” (Figure 20d), it is possible to extract the vectors that plot the magnitude`s limit graph, and the values of order (M), the $-3$ dB cutoff frequency ($f_{3dB}$), and coefficients of a linear-phase FIR filter that meet the requirements of the application (obtained using the algorithms proposed in Section 4.1 of this paper).

6. Experimental Validation

To validate the Matlab app and the proposed algorithms for designing FIR filters for RCs, the experimental results of a real system are presented below. The chosen system is a three-wire three-phase shunt active power filter (SAPF), whose block diagram is shown in Figure 21.

Table 1. Parameters of the repetitive controller and experimental setup.

Parameters of the Repetitive Controller			Parameters of the Experimental Setup (Figure 21)
n	m	$N=f_s/f_g$	$V_{g\left(line\right)}$	$L_g$	$R_g$	$L_l$	$L_f$	$R_f$	$R_{load}$	$V_{dc}$	$f_s$ *	$f_g$
6	1	288	380 V${}_{rms}$	$186.17$$\mathsf{\mu }$H	$31.7$ m$\Omega$	$1.483$ mH	$2.563$ mH	$307.5$ m$\Omega$	$48.4$$\Omega$	600 V	$17.28$ kHz	60 Hz

* Sampling and switching frequency.

shows the parameters of the SAPF and the parameters of the repetitive controller presented in [23], which were used to validate the proposed algorithm.

6.1. Description of the Experimental Setup

A SAPF is a power electronic device used to mitigate the harmonic contamination caused by non-linear loads. In this paper, a three-wire three-phase inverter is used as SAPF. As can be seen in the block diagram of the complete control system (Figure 21), a three-phase diode rectifier with resistive load (no filter on the DC side of the rectifier) is used as the non-linear load. In this configuration, an inductive filter is used on the AC side of the three-phase rectifier ($L_l$). Even though industrial three-phase rectifiers are typically connected to predominantly inductive DC loads, several works in the literature use experimental platforms to evaluate three-phase APFs, considering loads comprised of three-phase rectifiers with small filter inductances at the AC side and a DC resistive load, as conducted in this paper. Examples of works that use this assemble are [3,16,34,35].

In this application, the output currents of the SAPF are the controlled variables of the evaluated control system. Several papers have shown that the current control strategy applied in this system has a high impact on the system’s performance [22,34,36,37]; therefore, in order to design an adequate controller for this application, it is necessary to obtain its mathematical model. As can be seen in [16], when controlling the three-phase output currents of a SAPF, the converter can be modeled by its output filter. Thus, in an $\alpha \beta$ reference frame, the plant transfer function is given by:

$$G\left(s\right)=\frac{{\overrightarrow{I}}_f\left(s\right)}{\overrightarrow{D}\left(s\right)}=\frac{V_{dc}}{L_{f}s+R_f},$$

(16)

in which ${\overrightarrow{i}}_f={\mathcal{L}}^{-1}\left\{{\overrightarrow{I}}_f\left(s\right)\right\}$ is the space-vector obtained from the SAPF output currents and $\overrightarrow{d}={\mathcal{L}}^{-1}\left\{\overrightarrow{D}\left(s\right)\right\}$ is the space-vector obtained from the duty cycles.

The reference generator and the current controller were implemented in a dSPACE platform (shaded area of Figure 21). The block $H_l\left(z\right)$ is a lead compensator used to attenuate the effects of the computational delay [31] at the desired 0 dB crossover frequency (1.73 kHz—around 10 times lower than the sampling frequency) of the open loop transfer function. The transfer function of $H_l\left(z\right)$ is given by:

$$H_l\left(z\right)=\frac{0.6526z-0.4301}{z-0.08271},$$

(17)

which was designed to cause a phase lead of $36.0^{\circ }$ around 1.73 kHz. The magnitude attenuation caused by the lead compensator for low frequencies is naturally compensated when tuning the repetitive gain $K_{rc}$. Figure 22 presents the phase shift caused by the computational delay ($z^{-1}$) and the phase shift caused by both the computation delay and the lead compensator together ($z^{-1}·H_l\left(z\right)$), indicating that the $H_l\left(z\right)$ used in the experiment can successfully attenuate for the phase shift caused by the computational delay for frequencies up to the desired 0 dB crossover frequency. The DC bus voltage control loop is not represented in Figure 21) for simplicity purposes. However, in practice, there is an external voltage control loop, which is significantly slower than the internal current control loop.

The dSPACE platform is a modular hardware system that allows the user to program a control system application using Simulink/Matlab. This platform is composed of I/O and processor boards. The dSPACE platform used in the experimental setup has the DS1005 PPC board featuring a PowerPC 750GX processor running at 1 GHz. Figure 23 shows the complete experimental setup.

6.2. Description of the Evaluated Repetitive Control System

Since the control system must synthesize the non-fundamental frequency currents demanded by the non-linear load, the harmonic spectra of these currents must be evaluated in order to choose which controller can be used in this application. However, as it is a three-phase current, one can evaluate the space-vector obtained from the load currents and, thus, represent all three phases in a single plot [22]. Figure 24 presents the harmonic spectrum of the load’s current space-vector.

When evaluating the harmonic spectrum of a space-vector signal, one must know that positive frequencies represent positive-sequence harmonic components while negative frequencies represent negative-sequence components [22]. Therefore, based on Figure 24, one can use a $6k$ + 1 RC, which is complex repetitive controller that applies high gain in the harmonic components of the set $H=\{\cdots ,-17,-11,-5,1,7,13,19,\cdots \}$. In this paper, the $nk+m$ RC proposed in [23] is used with $n=6$ and $m=1$, as shown in Table 1. It is important to know that, according to [18], this $nk+m$ RC is equivalent to the PRC described in Section 2; therefore, parameter a is an additional degree of freedom for the control system designer.

Based on the delimitation presented in the previous paragraph, and considering the repetitive controller presented in [23], the current controller used to validate the impact of $Q\left(z\right)$ (indicated in Figure 21) is given by:

$$C\left(z\right)=K_{rc}\left[a+\frac{e^{j2\pi \frac{1}{6}}{z}^{-\frac{N}{6}}Q\left(z\right)}{1-e^{j2\pi \frac{1}{6}}{z}^{-\frac{N}{6}}Q\left(z\right)}\right],$$

(18)

for $Q\left(z\right)$ being a zero-phase low-pass FIR filter. Since $N=f_s/f_g$ is the number of samples per period of the fundamental frequency, which is given in Table 1, then (18) becomes

$$C\left(z\right)=K_{rc}\left[a+\frac{e^{j\frac{\pi }{3}}{z}^{-48}Q\left(z\right)}{1-e^{j\frac{\pi }{3}}{z}^{-48}Q\left(z\right)}\right].$$

(19)

Parameters a, $K_{rc}$, and $Q\left(z\right)$ are tuned using the created Matlab app to validate the proposed algorithm, as indicated in the following subsections.

The created Matlab app requires the user to define which repetitive gain $K_{rc}$ is going to be used in the control system. This is done at the start screen of the Matlab app, before any analysis of the system stability and performance. In the evaluation presented in this paper, $K_{rc}=0.06$ is used.

However, one can interactively use the Matlab app to tune the best $K_{rc}$ for the control system. This interactive use would be conducted as follows:

(i)

Firstly, the user would define a $K_{rc}$ gain for the initial analysis in the start screen of the Matlab app;

(ii)

Then, the user would evaluate the stability of the system using parameter a and a constant attenuation—using the stability screen—or a FIR filter—using the FIR filter screen;

(iii)

Using the obtained parameters, the system must be simulated to verify whether the performance requirements are met or not;

(iv)

If the requirements are not met, the control system designer must increase or decrease the $K_{rc}$ gain to make the system faster/slower and redo the process from step (ii).

The effects of the computational delay and PWM delay can make the system unstable for most repetitive gain values. In order to solve this problem, in step (i), the control system designer can add a lead compensator tuned to work around the desired 0 dB crossover frequency (as conducted in [31]). By doing so, the effects of the delays are mitigated, allowing the use of higher $K_{rc}$ gains, resulting in faster dynamic responses [31]. When using the lead compensator $H_l\left(z\right)$, its transfer function must be included in $G_m\left(s\right)$, i.e., for a discrete control system,

$$G_m\left(z\right)=K_{rc}{H}_l\left(z\right){\mathcal{Z}}_{zoh}\left\{G\left(s\right)\right\}{z}^{-1}.$$

(20)

In Equation (20), $z^{-1}$ represents the computation delay, and ${\mathcal{Z}}_{zoh}G\left(s\right)$ represents the discretized version of $G\left(s\right)$ using the ZOH (zero-order hold) method. According to [38], the ZOH method introduces some dynamic characteristics into the system that are similar to those obtained when using a PWM (pulse-width modulation) modulator. Therefore, it can be used to represent the PWM dynamics in a simple way.

6.3. Validation of the Matlab app and the Proposed Algorithms

Firstly, it is important to obtain the $G_m\left(z\right)$, which can be used to evaluate the repetitive control system. Based on equations (16) and (20), and in the data presented in Table 1, $G_m\left(z\right)$ is given by

$$G_m\left(z\right)=K_{rc}{H}_l\left(z\right){\mathcal{Z}}_{zoh}\left\{\frac{V_{dc}}{L_{f}s+R_f}\right\}{z}^{-1},$$

(21)

where

$$\left\{\begin{array}{c}{\mathcal{Z}}_{zoh}\left\{\frac{V_{dc}}{L_{f}s+R_f}\right\}\approx \frac{13.5}{z-0.9931}\hfill \\ H_l\left(z\right)=\frac{0.6526-0.4301z^{-1}}{1-0.08271z^{-1}}\hfill \\ K_{rc}=0.06\hfill \end{array}\right..$$

(22)

These data are entered in the Matlab app using the start screen. To do so, one must enter $K_{rc}$ in the field shown in Figure 18b and $G_m\left(z\right)/K_{rc}$ in Figure 18c. Parameter $T_s$ can be calculated as $1/f_s$, which is given in Table 1. Based on these inputs, the validation of the Matlab app is organized into four steps, which are described below.

6.3.1. System Stability Analysis Varying Only parameter a (With $Q\left(z\right)=1$ Being a Constant Attenuation)

The stability domain of the repetitive control system was first observed by considering $Q\left(z\right)=1$ and $K_{rc}=0.06$. For this analysis, three distinct values for parameter a were considered: $a=0.4$; $a=0.5$; and $a=0.8$. The analyses were performed using the stability screen of the Matlab app presented in this paper.

As can be seen in Figure 25 (graphs exported from the Matlab app), despite the fact that the value of parameter a directly impacts the geometry of the stability domain, none of the obtained domains fully contemplate the Nyquist contour of $G_m\left(z\right)$. In fact, the Nyquist contour of $G_m\left(z\right)$ extends to the region on the left of the complex plane’s origin, thus, it is not possible to make the system stable for constant $Q\left(z\right)=1$. Since this system is unstable, only the simulation results are presented in Figure 25.

From the above analysis, one can see that the Matlab app can be used to understand and evaluate the effects of parameter a on the system stability.

6.3.2. System Stability Analysis Varying Parameters a and $Q\left(z\right)$ (With $Q\left(z\right)$ Being a Constant Attenuation)

As presented in Section 2, an easy solution to make the system stable is achieved when using a constant attenuation $Q\left(z\right)<1$ in the repetitive controller. The Matlab app proposed in this paper can be used to conduct this evaluation. Once parameter a is chosen, the control system designer can vary the constant value of $Q\left(z\right)$ in the stability screen. By doing so, one can see that the stability domain of the repetitive control system is altered. Based on Section 2, it can be seen that as $\left|Q\right(z\left)\right|$ decreases, the system becomes more stable, however, it loses its capacity to achieve the zero steady-state error.

For the purpose of comparison, a value of $a=0.5$ is used for the following analysis. To evaluate the impact of $Q\left(z\right)<1$, two distinct values of $Q\left(z\right)$ are considered: $Q\left(z\right)=0.9$ and $Q\left(z\right)=0.6$. From Figure 26a,c, it is expected that the system only remains stable in the second case (for $Q\left(z\right)=0.6$), which is confirmed when observing the corresponding phase-a grid current (Figure 26d). In the first case ($Q\left(z\right)=0.9$), despite the visibly better behavior (compared to when $Q\left(z\right)=1$), the current does not exhibit stable behavior. In the second case, although the grid current is still not an ideal sinusoidal signal, it remains stable.

As explained in Section 2, the error signal does not diverge or converge to zero when using a constant value for the parameter $Q\left(z\right)$, which can be observed in Figure 27. In fact, this solution can reduce the total harmonic distortion (THD) of the grid currents, as shown in Figure 26d; however, the resulting THD does not meet grid requirements (it must be lower than $5\%$) [39]. Therefore, using $Q\left(z\right)$ as a zero-phase FIR filter is an interesting choice to repetitive control systems.

From the analysis presented in this subsection, one can see that the Matlab app can be used to understand and evaluate the effects of parameter $Q\left(z\right)$ (as a constant attenuation) on the system stability.

6.3.3. Designing a FIR Filter for the RC System

In this section, the proposed algorithms that define the magnitude limit of $Q\left(z\right)$ will be used. The parameters entered in the FIR filter screen and the results are displayed in

Table 2. Parameters and results obtained from executing the proposed algorithms using the created Matlab app.

Parameters
Parameter a	1
Maximum initial value	1
Decrement step ($\Delta q$)	$0.005$
Sampling rate	$17.280$ kHz
Initial frequency	100 Hz
Final frequency	10 kHz
Number of evaluated frequencies	1000
Results
Filter Order (M)	6
$-3$ dB Cutoff Frequency ($f_{3dB}$)	$2.744$ kHz

.

The best filter to be used in this application is the one whose magnitude response is as close as possible to the magnitude limit curve. However, sometimes this ideal filter can be too costly to be implemented. For this reason, the Matlab app allows the exportation of the generated curve, in order to leave the user free to design a less costly filter (lower order and cutoff frequency). Based on the magnitude limit curve of $Q\left(z\right)$, after some tests using the MATLAB function fir1(), suitable for generating FIR filters, it was observed that a FIR filter of order 6 and cutoff frequency equaling $1.8$ kHz meets the requirements defined by the curve generated by the Matlab app (Figure 28). The estimation of the order of the FIR filter by the proposed algorithm is not applied in this example (STEP 7 is not being evaluated for now).

Based on the strategy presented in the previous paragraph, the following symmetrical FIR filter was obtained:

$$Q\left(z\right)=0.01269z^3+0.07715z^2+0.2415z+0.3372+0.2415z^{-1}+0.07715z^{-2}+0.01269z^{-3}.$$

(23)

This filter was implemented in a RC to control the output currents of the SAPF. However, as presented in Section 2.1, in order to make this filter implementable, the number of delays in the periodic signal generator of the RC must be reduced to $N-M/2$, i.e., it must be reduced in 3. This means that (23) can be implemented as:

$$Q\left(z\right)=z^3\left(0.01269+0.07715z^{-1}+0.2415z^{-2}+0.3372z^{-3}+0.2415z^{-4}+0.07715z^{-5}+0.01269z^{-6}\right).$$

(24)

The experimental results are shown in Figure 29, where $i_{ga}$ is the phase-a grid current, $i_{la}$ is the phase-a load current, and $i_{fa}$ is the phase-a output SAPF current.

It is still possible to observe the impacts of varying the parameter a on system stability. For that, the values of all other parameters are maintained constant while testing the following configurations: $a=0.8$; $a=0.6$; and $a=0.4$. Figure 30 shows the results of these tests. Figure 30a,c,e show the magnitude limit for parameter $Q\left(z\right)$ (in green) and the magnitude response curve of the filter used in the experiment (in red). Figure 30b,d,f show the phase-a grid (in red), load (in blue), and SAPF (in orange) currents. As seen, for the first two scenarios, the FIR filter is still within the magnitude limits, resulting in stable and well-controlled currents, unlike when $a=0.4$, in which the grid current shows increasing oscillation, indicating that the system tends towards instability.

From the analysis above, one can see that the Matlab app can be used to generate the magnitude limit of $Q\left(z\right)$, which is necessary to design the FIR filter of a RC. Experimental validation is used to show the effectiveness of the proposed algorithm and the created Matlab app.

6.3.4. Error Analysis for Distinct Choices of Filter Characteristics

In the previous subsection, the FIR filter used in the control system was chosen from the magnitude limit curve of $Q\left(z\right)$, generated by the Matlab app. Although the designed FIR filter meets the application requirements, it was not the solution presented by the Matlab app, whose parameters are presented in Table 2. Thus, in this subsection, the performance of the repetitive control system is compared for two distinct FIR filters, which are (i) the solution presented in (23); and (ii) a FIR filter designed based on the “result data” in Table 2.

According to the Matlab app (Table 2), the best FIR filter for this application is the one with an order of 6 and a cutoff frequency of 2744 Hz (Figure 31). The algorithm implemented in the Matlab app to obtain these values is presented in Figure 15.

Using this algorithm, the Matlab app can find a FIR filter that satisfies the requirements imposed by the magnitude limit curve while keeping the filter gain as close to unity as possible for the widest possible range of frequencies. It is worth noting that, as mentioned in [23], in frequencies where the filter has a non-unity gain, the system tends to have a non-zero steady-state error. This behavior is discussed in Section 3 of this paper. Therefore, the best filter is the one with the widest possible unity-gain passband.

Figure 32 shows the magnitude limit curve of $Q\left(z\right)$ (in blue), the magnitude response of the FIR filter suggested by the Matlab app (in red), and the magnitude response of the FIR filter implemented by (23). From these curves, one can see that both evaluated FIR filters meet the magnitude limit graph demanded by the RC system. The magnitude response of the FIR filter suggested by the Matlab app has a larger passband, leading to a wider range of frequencies that the filter gain is close to unity. As a consequence, this solution is expected to result in a smaller steady-state error, which is proved to be true when evaluating the THD of the two repetitive control systems (Figure 33). Simulation results are used to ensure similar electrical grid characteristics and the same window of data to calculate the THD for all evaluated cases.

From the analysis above, one can see that the Matlab app can suggest a FIR filter that keeps the repetitive control system fully functional. Due to the effectiveness of the proposed algorithm, the obtained FIR filter leads to an optimized steady-state behavior. A performance comparison between all results obtained in this section is presented below.

6.4. Performance Comparison

The THD results of the phase-a grid currents obtained in this section are presented in

Table 3. Performance comparison of the evaluated SAPF considering the RC-based controller proposed in [23]. THD of the phase-a grid current obtained for all $Q\left(z\right)$ evaluated in this paper. IAE and ITAE performance indices obtained comparing phase-a SAPF and reference currents for all $Q\left(z\right)$ evaluated in this paper. Results are obtained from the simulation.

Cases	Evaluated Scenarios	THD	ISE ($\times {10}^3$)	ITAE ($\times {10}^3$)
$\#1$	Grid currents without SAPF	$27.10\%$	–	–
$\#2$	SAPF controlled by a RC with $Q\left(z\right)$ as a constant attenuation ($Q\left(z\right)=0.6$).	$9.57\%$	$87.46$	$7.01$
$\#3$	SAPF controlled by a RC with FIR filter $Q\left(z\right)$ presented in Section 6.3.3 ($M=6$ and $f_{3dB}=1.800$ kHz)	$4.77\%$	$15.27$	$2.88$
$\#4$	SAPF controlled by a RC with FIR filter $Q\left(z\right)$ suggested by the Matlab app ($M=6$ and $f_{3dB}=2.744$ kHz)	$4.60\%$	$14.27$	$2.80$

. All THD data were obtained from simulation results in Matlab/Simulink software. This was done to ensure similar electrical grid characteristics and the same window of data to calculate the THD for all evaluated cases. From Table 3, one can see that the use of the proposed Matlab app (cases $\#3$ and $\#4$), especially when using the suggested FIR filter, results in a significant reduction in the THD of grid currents.

Although cases $\#3$ and $\#4$ present similar THDs, it is worth mentioning that in case $\#3$, the Matlab app was only used to generate the superior magnitude limit curve of $Q\left(z\right)$ (the FIR filter was obtained via manual tuning). Therefore, the duration of the design process was significantly higher than the duration of the design of $Q\left(z\right)$ for case $\#4$, which was the filter suggested by the Matlab app (based on STEP 7 and STEP 8 of the algorithm proposed in this article). Furthermore, as expected, when evaluating Figure 32, the $Q\left(z\right)$ used in case $\#4$ results in a lower THD than the $Q\left(z\right)$ used in case $\#3$, given its larger bandwidth. Finally, it is also worth noting that the designer can fine-tune the filter suggested by the Matlab app (increasing its cutoff frequency and/or order) to further improve the performance of the control system.

Two other performance indices are also compared in Table 3, they are: (i) the integral of the square of the error (ISE); and (ii) the integral of time multiplied by absolute error (ITAE). The ISE index is calculated using the squared error of the control system, having as its main contribution the transient error. This happens because the steady-state error tends to zero for controllers based on the internal model principle (such as the repetitive controller), resulting in a significantly lower squared error for steady-state response when compared to transient response. This means that the ISE index can be used to compare the transient performance of control systems. On the other hand, the ITAE index weighs the absolute error using time, as a result, having the stationary error as the main contribution of this index. In the analysis presented in Table 3, an interval of 0.2 seconds (after enabling the SAPF operation) was used for computing the ISE and ITAE indices (the lower the better).

From Table 3, it is observed that case $\#4$ results in the lowest ISE and ITAE indices among all. This means that when using the proposed algorithm in this paper, it was possible to achieve the lowest transient and stationary errors among evaluated cases.

7. Conclusions

This paper presents a detailed, step-by-step algorithm for designing zero-phase FIR filters to be used in RC-based control schemes. Additionally, this paper presents a Matlab app developed to perform the stability analysis of RC systems and the design of their FIR filters. The proposed algorithm enables the user to obtain functional parameters for the FIR filter of an RC-based controller ($Q\left(z\right)$), thus avoiding gross errors that could lead to system instability. A shunt active power filter is used to validate the proposed algorithm. The simulation and experimental results confirm that:

• The greater the passband of the FIR filter $Q\left(z\right)$, the smaller the steady-state and transient errors of the RC system;
• Despite the characteristics mentioned in the previous item, the cutoff frequency of the FIR filter $Q\left(z\right)$ must not be high to the point of violating its superior magnitude limit. This limit can be obtained from the system stability analysis using the proposed Matlab app;
• The algorithm presented in this article is functional and it results in a FIR filter $Q\left(z\right)$ that does not violate its magnitude limits. Furthermore, the algorithm chooses a cutoff frequency for the filter that enables better steady-state results than most filters.