Characterization of limit cycle oscillations induced by Fixed Threshold Samplers

In this work, a generalized study of the conditions for the appearance of limit cycle oscillations induced by any kind of sampler with multilevel fixed thresholds is presented. These kinds of samplers, which will be referred to as Fixed Threshold Samplers (FTS), are characterized by a series of parameters, which, when selected properly, allow obtaining some of the most used forms of quantization in Event-Based Control (EBC). Because of some sampler characteristics, the obtained limit cycle oscillations can present a bias, therefore, to characterize them the Dual Input Describing Function (DIDF) method is used. The obtained DIDF is analyzed revealing some interesting properties allowing to simplify the robustness analysis. The analysis takes into account the effect of the disturbance and reference signal influence on the system, generally overlooked in DF analysis. Guidelines about how to perform the robustness analysis are given, showing their application through some study cases.


I. INTRODUCTION
In recent days, Event-Based Control (EBC) is becoming a more and more popular control alternative. This change of paradigm is motivated by the apparition of a new industrial framework named Industry 4.0. This new communication benchmark requires efficient data flow through the industrial networks, fostering the implementation of EBC algorithms.
Several authors have already applied the EBC to well known control algorithms, notably to the PID algorithm, which constitutes one of the most simple and reliable algorithms, and therefore, it is placed among the most used in the industry [1]. In addition, its role within the context of Industry 4.0 has been brought forward in [2].
One of the first works referring to Event-Based PID control algorithm was Årzén's contribution [3], in which it was shown how this algorithm could be used as a tool to reduce the CPU usage to perform the control of a process without affecting significantly its performance. In that work some key aspects about EBC were identified, specially, the error involving the calculation of integral and derivative terms when the elapsed time between samples increases.
Subsequently, some works have been addressed to resolve these problems, mainly regarding the integral time calculation. In that sense, the contributions made by Durand [4], [5] and Vasyutynskyy [6], [7] must be highlighted. Other works like [8] show the validity of this kind of algorithms when implemented with the standard IEC 61499 [9], which supports the implementation of distributed control systems.
As important as the control algorithm is the sampling strategy used to regulate them, which is in charge of generating events whenever significant changes are detected in the state of the controlled system. Among the most used strategies the Fixed Threshold Samplers (FTS) can be found. Some examples of this kind of sampling strategies are the Regular Quantization (RQ), the Symmetric-Send-On-Delta SSOD [10] or the Regular Quantization with Hysteresis (RQH) [11] which represents an intermediate case.
A lot of works have been addressed considering these strategies into the loop. There can be found works in which they have been used for identification purposes [12]- [14], but more importantly, for the control of processes [15]- [17] and, for a correct performance, several tuning rules and procedures have been proposed [18]- [21].
Most of the aforementioned works deal with a drawback of these kinds of FTS, which is that they can induce limit cycle oscillations. Under the perspective of EBC these oscillations degrade the loop performance, accelerate the wear out of actuators and overload the communication networks with unnecessary events. Therefore, the study of the conditions for the existence of limit cycles as well as the characterization of the oscillation by predicting the amplitude and frequency is a key point in the analysis and design of event-based control systems.
The presence of limit cycle oscillations when a FTS is included in the control loop has been study following different approaches and assumptions. The characterization of limit cycles for different kind of systems, such as integrator processes plus time delay (IPTD), first order processes plus time delay (FOPTD), and second order processes plus time delay (SOPTD), when using a SSOD sampling strategy has been presented in [22]. In [19], [23] and [24], tuning methods for PI controllers with SSOD sampler have been developed based in new robustness margins for limit cycles, that were obtained by applying the Describing Function (DF) technique and entail with the classical concepts of phase and gain margins. The same approach was used in [25], where a unified design of SSOD-PID control architecture for selfregulating and integral processes was investigated. In [26] and [20], the authors proposed a new robustness measure to avoid limit cycle on SSOD based PI controllers. The proposal is based on the Tsypkin's method [27], which has been widely used to study relay control systems.
In this work, which extends the results in [28], the limit cycle oscillations induced by any kind of FTS are characterized. To that end, a variant of the Describing Function (DF) known as Dual Input Describing Function (DIDF) [29] is used as analysis tool. To perform the study the FTS is parametrized and the general DIDF is obtained. The results presented here are relevant in twofold: first, it allows avoiding the assumption that input signals to the control loop, that is, reference and disturbance, are equal to zero, which is a generalized assumption when applying the Describing Function technique in the previous works about FTS. However, this is an unrealistic premise for most of control systems. We prove that these inputs strongly determine the existence of limit cycles depending on the model structure of plant and controller. Consequently, this result shed light about the proper selection of the controller structure for a given plant in order to effectively avoid limit cycles. Secondly, the characterization of limit cycle oscillations presented in this paper could improve the use of FTS for system identification. As aforementioned, some papers have shown the effectiveness of φ x(φ) x(φ) FIGURE 1: Sinusoidal oscillation sampled with a given FTS (parameters: ∆δ/δ=1/2, ∆ / =1/4, /δ=2, h/δ=1/2, explained below). using the oscillation induced by FTS for estimating process models. In those works, symmetric and unbiased oscillation are assumed. It is known, however, that this kind of oscillation are hardly ever obtained in actual plants, furthermore, previous works has demonstrated that the asymmetric behavior can improve the identification results [30]. In that sense, the study presented in this paper will allow to extend the previous works on FTS based identification for the case of asymmetric and biased oscillation.
The paper is organized as follows. Section II presents the generalized FTS, defining its characteristic parameters, and shows how other samplers can be obtained by choosing the correct parameters. Section III presents the DIDF approach and some preliminary properties necessary to develop further analysis. Section IV studies the behavior of different type of process under FTS loops and how to analyze them. An experimental case validating the analysis methodology is presented in Section V. In Section VI it is discussed the effectiveness of the presented DIDF analysis for the design of control systems. Finally, the conclusions about this work are drawn in Section VII.

II. PROBLEM STATEMENT
The choice of a proper sampler, or the task to adequate it to the given necessities, can be challenging, specially considering that the mere fact of quantifying can induce limit cycle oscillations.
In addition, as the complexity of the sampler increases, the resulting oscillations can become difficult to predict. For example consider the oscillation presented in Fig. 1 produced by a given FTS. In this figure, it can be observed a huge difference between the input signal x(φ) and the output of the FTSx(φ). Among the main differences it can be seen that the quantization step is higher than the increase in the signal, the switching thresholds do not start at the origin, the output of the FTS is not symmetric and the number of levels crossed differs between upwards (four) and downwards (three) direction. All this characteristics can be parametrized by regarding the relationship between input (x) and output (x) of a general FTS, which is summarized graphically in Fig. 2. In this figure, it can be seen all the parameters that characterize a FTS, namely, the detection width δ, the detection hysteresis h, the detection asymmetry ∆δ, the output quantification and the output asymmetry ∆ . The asymmetries have been measured with regard to the center of the rectangle δ × and the hysteresis h is vertically centered on that rectangle. It is considered without loss of generalization that ∆ ∈ [− /2, /2] and ∆δ ∈ [−δ/2, δ/2].
These parameters already prognosticate the behavior of the sampler and the system. For example, the sampler will contain a dead zone only if ∆ = ± /2, it will present a symmetric output only if ∆δ = ±δ/2 or it will present a certain amplification or attenuation depending on the relationship between and δ. Additionally, it is worth noticing that by choosing a certain combination of parameters the most common FTS samplers traditionally studied can also be found. A summary with some of them can be found in Table 1. Obviously, this table only presents some cases of multileveled samplers, there exist many other options that can be obtained with the appropriate configuration.

III. DIDF APPROACH
As it has been previously commented, one of the most important phenomena when using a FTS is that they can induce limit cycle oscillations. To study this kind of nonlinear behavior there exist several tools, among them, the Describing Function approach excels due to its ease of use and accuracy in the obtained results [23], [31].
The generic FTS under study includes, depending on the chosen parameters, a wide variety of commonly used samplers, SSOD [10], RQH [11] or an RQ sampling among them, as it has been shown in the precedent section. Therefore, by addressing an analysis considering the most generic case a wide variety of samplers can be characterized. Considering the parameters previously described, the relationship Multilevel relay with hysteresis ∆ = 0, =δ=h between the input signal to the sampler x(t) and its output x(t) is defined in (1). To apply the DF technique consider a generic control loop like the one presented in Fig. 3. Consider that the controller C(jω) and the process G(jω) can be grouped in a single block G ol (jω)=C(jω)G(jω) representing the open-loop transfer function which gathers the linear behavior on the loop. Because of the inherent generalization induced by the FTS, there may exist some possible bias in the output of the sampler which may not be attenuated through the loop and, therefore, it must be considered at the input of the non-linearity. Thus, in this type of oscillatory steady-state the input to the sampler corresponds to the form x(t) = A sin(ωt) + B. Therefore, to include in the study the cases where a bias appears, a variant of the DF is used called Dual Input Describing Function (DIDF), which studies the sustainability of VOLUME 4, 2016 Applying the DIDF principles, it is known that the condition for the maintainability of an unattenuated oscillation, called limit cycle, in this kind of systems is given by: where N A is the describing function of the non-linear element related to the dynamic behavior. This condition can be interpreted graphically as: if it exists intersection between G ol (jω) and the inverse negative of N A in the Nyquist diagram, a limit cycle can occur. In addition to the condition for the dynamic part presented in (2), the sustainability of the bias through the loop should also be studied. The conditions involving the bias maintenance were presented in [32] and, as in the oscillation maintenance, also depend on the linear part of the process under study. However, those conditions did not take into account the effect of the disturbance and reference input, which were considered to be r(t) = 0 and p(t) = 0, therefore, they have to be adjusted to take into account these signals.
Consider the loop structure presented in Fig. 3 and assume reference and disturbance inputs to be step-like signals. Assume also that a biased oscillation appears, so that the bias on x(t) is B. Then, considering only the non-oscillatory parts of the elements in the loop: where r(∞) and p(∞) are the magnitude of the changes in reference and disturbance signals, G(0) and C(0) the steady-state gains of the process and the open-loop transfer function, B the bias and N B the describing function for the maintenance of the bias. Arranging the previous expression, the condition for the maintenance of the bias is obtained: where G ol (0) is the steady-state gain of the open-loop transfer function. This expression depends on the form of G(s) and C(s) since it depends on the steady-state gain of both.
For several combinations between the form of G(s) and C(s), expression (3) has been further developed in Table 2 of Appendix B.
The DIDF of the general FTS involving the bias maintenance and the dynamic oscillations are given respectively by expressions (4)-(5) (calculated in Appendix A), where i 0 is the initial or base level, around which the oscillations are produced. The oscillations can have a number of m sup levels crossed in upwards direction and m inf in downwards direction from i 0 . These parameters are calculated as follows: and: The DF for the bias and for the sinusoidal part can be expressed with a series of meaningful dimensionless ratios, namely, /A, /B, δ/A, δ/B, ∆ / , ∆δ/δ and h/δ. The ratios that describe the sampler are bounded: ∆ / ∈ [−1/2, 1/2], ∆δ/δ ∈ [−1/2, 1/2] and it will only be considered h/δ ∈ [0, 1] even though it can take a wider range. The other ratios involve the input signal and, therefore, they cannot be bounded.
By substituting in the previous expressions the parameter configurations presented in Table 1, the dynamic describing function (N A ) considering B = 0 can be found for the RQ [31], RQH [11] or SSOD [23] for example.

A. CHARACTERISTICS OF THE DYNAMIC PART OF THE DIDF
Due to the generality of the expression presented above and the number of parameters involved, the analysis study will result in a complex research. Therefore, to bound this study, N A is evaluated analyzing its behavior under certain conditions. This evaluation leads to a better comprehension of the behavior of N A and to a significant reduction of the analysis study effort.
The first behavior has been already presented in Fig. 1, where the oscillations present different number of levels crossed upwards and downwards from i 0 , i.e. m sup = m inf .
Nevertheless, the same sampler can provide oscillations with the same number of levels crossed upwards and downwards, m sup =m inf , depending on the amplitude A of the oscillation. This behavior is reflected in the dynamic part of the DF by presenting two alternating "series" of bands, one where m sup = m inf and another with m sup =m inf alternatively. That can also be seen evaluating the expressions in (6) for a range of values of A.
On the top image of Fig. 4, the inverse negative of the dynamic DF N A (see (5)) of a sampler with ∆ / =∆δ/δ=1/3 and =δ=h and without bias has been represented. The presented traces are grouped in two series that have been plotted in different colors. In that image, two squares, one in each series, have been marked, which correspond to the oscillation presented on the images below. As it can be seen, the orange oscillation only has one upper level and no levels are crossed downwards from the initial level, whereas the violet oscillation has the same number of levels crossed in both directions.
Despite the ratio /δ not being a ratio used to define the DF, as it can be seen from (4) and (5), it has a great influence on the placement of the DF traces. This ratio describes the relationship between the detection and the output thresholds of the sampler, therefore, it has the same behavior as a gain. Thus, the DF traces present a radial shrink when this ratio increases and a radial expansion when it decreases. Therefore, this parameter is of great influence to the robustness of the loop since it restricts the frequency response of the linear elements. This radial variation of −1/N A has been represented in Fig. 5 where the inverse negative of a sampler with ∆ / =1/4, ∆δ/δ=1/3 and h/δ=1/2 has been presented using three different ratios /δ.
Another parameter which plays a paramount role is the relative hysteresis, presented by the ratio h/δ, which prevents the generation of unnecessary events due to the presence of noise in the sampled signal [11]. In addition, regarding this parameter from the non-linear analysis point of view, it makes the sampler take into consideration from where does the sampled signal come, i.e. it adds a memory effect. As it was stated in [32], non-linearities with memory present an imaginary part in the inverse negative of their DF and those without memory only have a real part. This fact is  Fig. 6, in which a sampler with ∆ / =1/5, ∆δ/δ=1/2 and /δ=0.75 has been tested with three different levels of relative hysteresis, showing that the traces fold gradually towards the real axis as the hysteresis is reduced.
In addition to the influence of the parameters of the sampler, in those cases where a bias is present, i.e. the input to the non-linearity is of the form x(t)=A sin(ωt)+B with B = 0, it plays an important role on the shape of the DF under study.
The presence of this bias can change completely the be-VOLUME 4, 2016  havior expected of a given sampler. For example, imagine any FTS with a deadband (SSOD, RQH, etc.). If the loop contains any element that accomplishes a bias B such as ∆δ = B, then the behavior of the FTS will not correspond to a deadband, but it will be more similar, for example, to a multilevel relay, i.e. a sampler with no deadband. This fact is logical because the value of the bias B compensates the asymmetry ∆δ, and therefore the same behavior is expected.
Traduced to the analysis of N A , this fact implies that for a given sampler in which ∆δ = B a straight line that reaches the imaginary axis will appear, as it appears for other nonlinearities without deadband, for example, the ideal relay. This effect of the asymmetry compensation with the bias is illustrated in Fig. 7, where the parameters of the SSOD are used and it can be seen how an horizontal line that reaches the imaginary axes appears. The apparition of a bias can also propitiate a lag between the upwards and downwards crossed levels. This can be seen from expressions in (6): having a given sampler makes That is to say that the apparition of a bias is equivalent to an asymmetry.
For the analysis study it can be very concerning this big variance of the traces of −1/N A with the bias B. Nevertheless, it has been observed that the placement of these traces is periodically repeated with the change of B. Concretely, the placement of the traces will be the same for: where B = B+kδ and k ∈ Z. This result is logical because it only represents a shift, and it can be observable from the expressions of N A and of the levels crossed m sup and m inf where the ratio δ/B and the levels are presented together. Besides this periodicity, it has also been observed that −1/N A presents the same traces for given variations of B lesser than δ. Concretely, the obtained traces are the same if B changes with regard to some symmetry axes, being these symmetry axes the extreme and the center of the square δ × . Therefore, the obtained traces will be the same for B=∆δ − δ/2+α and for B =∆δ − δ/2 − α, α ∈ . And regarding the other symmetry axis, the same traces will be obtained considering B=∆δ+β and B =∆δ − β, β ∈ .
With regard to N B , no generic properties can be extracted from its study as for the case of N A . The application of the part of the DIDF related to the bias will be treated deeply in the next section.

IV. OSCILLATION CONDITIONS STUDY
The study will pursue different objectives depending on the studied sampler, its application and the process involved. For example, it may be desired that a given system presents an oscillatory response around a given set point with controlled oscillations or that a system presents a robust behavior against these limit cycle oscillations that the DIDF describes. This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. To perform the analysis it is necessary to know the bias and if it is not attenuated through the loop even though the exact bias is not known a priori. However, it can be estimated by considering that the changes produced on the loop, whether their source is the disturbance or the reference, are big enough with regard to the detection threshold δ. In that case, the effect of the sampler can be substituted by /δ, which is the slope of input-output relationship (see Fig. 2). Then, an estimation of the bias, which will be referred to as bias central value B c , can be obtained from expression (3): From this central value, a set of candidate biases with amplitude δ will be considered defined as The amplitude of the set is chosen to be δ due to the symmetry property and periodicity of the dynamic part explained above.
Then, using the maintenance condition of the bias, which will be one of the presented in Table 2 depending on the elements of the loop, for the values in B s , pairs {A, B} which are candidate to present limit cycle oscillations are obtained for the given process, controller and sampler. Nevertheless, the pairs {A, B} must be evaluated to determine whether they conduce to an unattenuated propagation of the sinusoidal part of the signal. This evaluation will be done with the expression presented in (2), which also depends on the bias B. If the pairs {A, B} that fulfill the corresponding expression in Table 2 are also those that intersect the open-loop transfer function in the Nyquist diagram, analytically: then a limit cycle oscillation of the type A sin(ω o t) + B will be obtained. This concept will be introduced through an example in the following subsection.

Consider a process with transfer function
and for illustrative purposes, consider that this process is placed in a loop with a sampler with ∆ / = 0, ∆δ/δ = 0, /δ = 1 and h/δ = 1/2, whose input-output relationship is shown in Fig. 8, with a controller C(s) = 1.
This case will correspond to the stability condition for the bias where both controller and process have neither poles nor zeros at the origin presented in Table 2. Adapting that equation to the current case, L 1 (0) = 1 and L 2 (0) = G(0), resulting in: Five scenarios will be considered: 1) Unitary reference input without disturbance. 2) Reference input of magnitude 0.95 without disturbance.
3) Unitary reference input and disturbance of magnitude 0.95. 4) Same as scenario 3 but increasing the gain of the process. 5) Same as scenario 3 but doubling the value of the sampler. Scenario 1: In this case, the equation for the maintenance of the bias is simplified, resulting in: The central value B c can be obtained from (7) as: Thus, the set of values to evaluate will be B s ∈ [B c − δ/2, B c + δ/2]. A width of the set equal to δ is chosen because of the symmetry properties of N A explained in the previous section. The set of values B s has been evaluated with (9). Regarding that expression, the left hand term varies with each term of the set B s evaluated, while the right hand term is constant. Therefore, the solutions A and B to that expression can be easily found by representing the left hand term, and the horizontal line representing the right hand term and evaluating the intersections between them. The graphical representation of the left hand term has been presented in Fig. 9a, each trace representing an item of the set. As it can be seen, it only exists one trace that equals to r(∞), the right hand term, which has been represented by a dashed black line. The trace that exactly matches the dashed black line corresponds to the trace obtained for B = 0.5. As all the values of that trace are equal to the searched solution, if the inverse negative of N A intersects at any point with the open-loop transfer function, then an oscillation will appear. The oscillation condition for the dynamic part has been presented in Fig. 9b where it can be seen that −1/N (A, B = 0.5) intersects with G(jω).
The temporal response of the process under the conditions presented in this scenario has been obtained, and the error and sampled error signals have been presented in Fig. 10.    The oscillation period has been measured in this figure and G(jω o ) has been represented in the Nyquist diagram in Fig.  9b with a red circle, proving the validity of the method. Scenario 2: To evaluate the maintenance of the bias, (9) remains valid because only the magnitude of the signal has changed. The procedure to obtain the set of values of bias B s is also the same, firstly the central value is obtained, in this case B c = 0.475, and then a margin of ±δ/2 for symmetry reasons is established.
The set of values B s has been evaluated with (9) and the graphical representation of the left hand term has been presented in Fig. 11a, each trace representing an item of the set. A dashed black line has been added representing the right hand term of (9). The intersections between right and left hand terms constitute a set of pairs {A, B} that will allow the maintenance of the bias in the oscillation. However, it must be evaluated if the dynamic part is also    N A (A, B), where the pairs {A, B} are solution of (9), have been encircled. The plots in gray are illustrative and not necessary for the stability analysis. As it can be seen, the open-loop transfer function does not intersect with any of the highlighted points, therefore, no oscillations will take place, and so it has been observed in the simulation presented in Fig. 12. Scenario 3: In this scenario, a disturbance signal is present on the loop, therefore, the expression to determine the central 8 VOLUME 4, 2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and  value B c obtained from (7) is: A set of values for the bias as in the previous cases have been evaluated, obtaining the representation of the left hand term of (8) presented in Fig. 13a. In this figure, the dashed black line presented corresponds to the right hand term of (8). A set of pairs {A, B} is obtained, and as in the previous case, it can be observed in Fig. 13b that no intersection between the inverse negative of the DIDF evaluated at those points and the open-loop transfer function is observed. Hence, no oscillation will take place. This result evinces that, in this type of loop architectures, the disturbance input can act both as a stabilizing or destabilizing agent. In Fig. 14, the previous statement can be seen. The first half of the experiment corresponds to the temporal response of the error signal to the reference input change, which also corresponds to Scenario 1, where it was proven the oscillatory behavior. However, by applying the disturbance change at t = 40 s the process manages to stabilize within the evaluated range of biases.
Scenario 4: In this scenario, there will be reference and disturbance changes and the process gain will be 4 times the original. For the sake of brevity the repetitive steps followed above are omitted.
The graphical representation of the oscillation conditions is presented in Fig. 15. In Fig. 15a, the graphical representation of the left hand term of (8) is presented with different colored lines and the right hand term of that equation is presented with a dashed black line. The intersections between both terms represent the solutions of the equations from  it can be observed that there exist an intersection between G(jω) and some of the highlighted points, revealing that a limit cycle oscillation can occur. The temporal response of this system has been obtained and it is represented in Fig. 16. In this simulation, the reference change has been applied at t = 0 s and the disturbance at t = 40 s. In this scenario, it can be seen that two different limit cycle oscillations are obtained, the first one with m sup + m inf = 1 resulting from the application of a step change, which has not been evaluated theoretically, but it could be done as in the Scenario 1. The second limit cycle is the expected result of the study in this scenario, with m sup =m inf =1 after the application of the disturbance change. The period of the resulting oscillation VOLUME 4, 2016 9 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and has been measured and G(jω o ) has been represented in the Nyquist diagram in Fig. 15b with a red square, validating the prediction of the DIDF method. The apparition of different type of oscillations modifying the disturbance or reference input signals can be interesting for the identification of processes, in which different types of oscillations with different oscillation periods can be used to identify multiple points of the frequency response of a process. However, for those cases the correct case study regarding the equations must be applied.
Scenario 5: In this final scenario, the ratio /δ, that so far in precedent scenarios was 1, has been doubled. Thus, this change must be introduced in the calculation of the central value: As in previous scenarios, a set of values B s is obtained to compute the stability conditions, obtaining the graphical representations presented in Fig. 17. From these figures it can be seen that the DIDF method predicts an oscillation. The system has been tested in simulation obtaining the temporal response presented in Fig. 18. The period of the obtained oscillation has been measured and G(jω o ) has been represented in the Nyquist diagram of Fig. 17b with a red square proving the validity of the prediction.
In this scenario, the opposite case to the presented in the Scenario 3 can be seen, starting from a stable temporal response as a consequence of the reference change, the application of the disturbance input to the system acts as a destabilizing agent, resulting in the apparition of a limit cycle oscillation.

V. EXPERIMENTAL CASE STUDY
In this section, experimental results are presented to validate the theoretical outcomes obtained in the previous sections. The experimental setup is composed of a DC motor and a PLC unit. The goal of the application is to control the speed of the motor, which is read from an encoder. To actuate on the motor a PWM signal is used, being able to regulate the functioning with its duty rate. The PLC rack contains the corresponding modules for reading encoder inputs and   producing PWM signals. Some electronics are required for adapting the power levels of motor and PLC. All the described elements are presented in Fig. 19. The loop configuration corresponds to the one presented in Fig. 3. The process G(s) describes the relationship between the DC motor speed and the PWM duty rate which is modeled by the following transfer function: , the controller C(s) is considered to be a proportional controller with K p = 768 and for the sampler a SSOD with δ = 3 is chosen. These last two elements are implemented within the PLC. The reference and disturbance signals are introduced also within the PLC scope. Three scenarios will be presented, which will be analyzed with the principles developed in previous sections: 10 VOLUME 4, 2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and  Since G(s) and C(s) do not have poles nor zeros at the origin, and regarding to Table 2, the sustainability of the bias is studied by the following equation: Scenario A: In this case, only a reference change is introduced, therefore, the theoretical study is conducted as in Scenario 1 and 2 from Section IV.
The DIDF study is presented graphically in Fig. 20, where the sustainability of the bias is presented in Fig. 20a and the sustained oscillations are studied in Fig. 20b.
The temporal response of the process under the conditions presented in this scenario is presented in Fig. 21, where the error signal x(t) and the sampled error signalx(t) have been plotted in black and teal respectively. From this oscillation, the pair (A, B) has been measured and, using this pair, the left hand term of equation (10) has been presented on Fig.  20a, and −1/N (A, B) in Fig. 20b, both with a red square.
As it can be seen, in Fig. 20a, the square that represents the measured limit cycle is placed near to the critical line, defined by the right hand member of equation (10). Besides, in   the observed frequency. Obviously, the accuracy on this estimation heavily relies on having an accurate process model, which lies beyond the scope of this contribution. Scenario B: The starting point for this scenario is the oscillatory state presented in the previous one, which has been already analyzed. In this case, a disturbance is introduced within the PLC as an increase of 500 units in the control action.
The theoretical analysis has been performed as in Scenarios 3-5 of the previous section, in which both reference and disturbance signals are present. The DIDF study is presented graphically in Fig. 22, and as in previous cases the dynamic part of DIDF is presented in Fig. 22b and the sustainability of the bias is studied in Fig. 22a.
As it can be seen in this case, even if the open-loop transfer function intersects some of the traces of the dynamic part of the DIDF, it is still far from the candidate points highlighted in Fig. 22b. These points have been obtained from the intersection in Fig. 22a between the left hand term VOLUME 4, 2016 11 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and   and the right hand term of equation (10), being the left hand term presented with different colors and the right hand term with a dashed black line. Therefore, from this study no limit cycle oscillation are expected under these conditions, which is corroborated experimentally from the temporal response of the system in Fig.  23. In this figure, it can be seen the error signal in black and the sampled error signal in teal, which start from the oscillatory state predicted in Scenario A and, once the disturbance is introduced, the limit cycle oscillation disappears.
Scenario C: This last scenario is similar to the previous one in the fact that reference and disturbance signals are involved in the analysis. Therefore, the theoretical analysis is performed as in Scenarios 3-5 of the previous section. The DIDF study is presented graphically in Fig. 24, in which the the dynamic part of DIDF is presented in Fig. 24b and the sustainability of the bias is studied in Fig. 24a.
The temporal responses of the error and sampled error signals under the conditions presented in this scenario are   presented in Fig. 25, where the error signal is presented in black and the sampled error signal in teal. As it can be seen, the system starts from a non-oscillatory state which is reached as a consequence of a reference change. However, once the disturbance is introduced, a limit cycle oscillation is induced on the temporal response. From the data obtained of these oscillations, the pair (A, B) has been measured and it has been used to calculate the left hand term of equation (10), which has been presented on Fig. 24a, and to compute −1 /N (A, B), which is shown in Fig. 24b, both points with a red square. As it can be seen, these squares are placed on a vicinity of the candidate oscillation points. With regard to the oscillation frequency, from the theoretical study presented in Fig. 24b, the open-loop transfer function crosses nearby the candidate points −1 /N (A, B) at a frequency ω o = 46.1 [rad/s]. This frequency is similar to the oscillation frequency observed experimentally, which has been measured to be of ω o = 47.242 [rad/s].
Regarding to the estimation results, it can be said that 12 VOLUME 4, 2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and the proposed DIDF method correctly predicts the apparition of limit cycle oscillations, since the real oscillation point is placed near the candidate points and the open-loop transfer function, and the forecast oscillation frequency is also similar to the experimental oscillation frequency.

A. COMPARISON WITH OTHER METHODS
Taking advantage of the theoretical and experimental results presented previously on this section, in this subsection, it would be presented a comparison between the proposed DIDF method and other methods that analyze the apparition of limit cycle oscillations when a SSOD sampler is used. In recent years, several contributions have published addressing the robustness issues arisen from SSOD sampling. The comparison will be conducted using as analysis method the DF method presented in [11], which defines robustness margins in terms of classical gain and phase margins, and [20], which analyzes the apparition of limit cycle oscillations using the Tsypkin method, which does not neglect the contribution of high order harmonics in the apparition of limit cycles as the DF method does.
The analysis has been conducted according to the principles stated in the cited references, obtaining the graphical representation of the robustness analysis presented in Fig. 26, in which the traces of the inverse negative of the describing function of the SSOD sampler have been presented in red and the critical Tsypkin branch B T (ω min ), to which the robustness margin M T is measured, in light green.
Graphically, it can be seen that no intersection occurs between the inverse negative of the describing function traces and the open-loop transfer function at any frequency, which indicates that a stable response is expected according to the DF method. Similarly, it can be seen that the critical Tsypkin branch B T (ω min ) does not intersect the open-loop transfer function point for which it has been obtained (G ol (jω min )) which has also been indicated.
Numerically, the classical DF method approach presents a phase and gain margins of Φ h/δ = 10.59 • and γ h/δ = 0.8271 [dB] respectively, and the Tsypkin approach presents a Tsypkin margin of M T = 0.08. Therefore, according to both methods, a non-oscillatory behavior is expected theoretically.
It could be stated that the obtained margins are insufficient to provide a robust response, justifying the limit cycle apparition observed experimentally. However, the margins that these methods provide evaluate the apparition of a type of oscillation that does not correspond with the type of oscillation observed experimentally. More specifically, these methods evaluate the apparition of limit cycle oscillations without bias, B = 0, and with the same number of levels crossed in each direction m inf = m sup (see Fig. 4), whereas the observed oscillation presents bias and different number of levels crossed in each direction m inf = m sup , thus, justifying the development of the analysis presented in this work.

VI. DIDF FOR EFFECTIVE DESIGN OF CONTROL SYSTEMS
Even thought the design of control systems with FTS is out of the scope of this paper, this section addresses how the results presented in the previous sections can be applied to avoid limit cycle in closed loop control systems or to regulate the limit cycle oscillations in systems controlled in a relay-like fashion.
Focusing on the avoidance of limit cycle oscillations, there exist several approaches that can be followed, including the modification of the input/output characteristic of the FTS or selecting and tuning a proper controller.
Regarding to the input-output characteristic of the FTS, it can be modified according to the dynamic DF behavior presented in Section III-A to avoid the intersection with the open-loop transfer function G ol (jω). For example, the ratio /δ can be reduced, which will push the inverse negative traces of the dynamic DF deeper into the third quadrant (see Fig. 5). Also, increasing the value of the hysteresis can lead to similar effects (see Fig. 6).
Regarding to the selection and tuning the controller, it must be taken into account that the model structure of G ol (jω) as well as the signals involved in the loop play an important role in the appearance of limit cycle, as it is summarized in Table  2, where three kinds of G ol (jω) can be differentiated: has not integrative nor derivative predominance, which correspond to cases 4, 7 and 12 in the table. One of these cases has been presented in Sections IV-A and V, which requires a full analysis with the DIDF as shown.
In these cases, the FTS can be modified as mentioned, by adjusting the ratio /δ and the hysteresis, to avoid limit cycle oscillations. • Cases where G ol (jω) is predominantly derivative, corresponding to cases 1, 2, 3, 6 and 11. In those cases, there exist only one candidate bias which is a function of reference and disturbance signals not depending on the describing function for the maintenance of the bias, N B . Therefore, only the inverse negative of the dynamic VOLUME 4, 2016 13 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. DF has to be analyzed, and the controller can be tuned to avoid intersection between −1/N A and G ol (jω). • Cases where G ol (jω) is predominantly integrative, which correspond to cases 5, 8, 9, 10 and 13. In cases 5, 9, 10 and 13 there is only one candidate bias, B = 0, and consequently the limit cycles can be sidestepped by tuning the controller to avoid intersection between −1/N A and G ol (jω). This result is the base of several tuning methods for PI under different FTS strategies [11], [19], [24]. In fact, cases 9, 10 and 13, where the controller has an integrator, corresponds to the structure of G ol (jω) assumed in those works, which include a PI controller. Only in case 8 the usage of N B is required, therefore, the analysis must be done considering both N B and N A . On the other hand, for systems controlled in a relay-like fashion, it may be desired to control a given limit cycle oscillation. Assuming that an oscillatory state has been reached, which is accomplished by fulfilling the DIDF conditions, it may be pursued a modification in the amplitude or frequency of the oscillations.
Regarding to the amplitude of oscillation, it can be modified by directly changing the output quantification and keeping the dimensionless ratios of the FTS detailed in Section III. This modification will result in a proportional change of the amplitude keeping the current oscillation frequency. This is sustained by the dimensionless properties of the DF, which will result in the same −1/N A traces without changing the point where they intersect with G ol (jω).
However, regarding to the oscillation frequency, its modification would entail a change in the intersection point between the dynamic DF traces and G ol (jω), implying a change also in the amplitude. Taking for example the Nyquist diagram presented in Fig. 9b, the intersection between −1/N A and G ol (jω) takes place at a given frequency. Changing that frequency would entail either to modify the DF traces to intersect another point, or modify the open-loop transfer function. Some of the strategies that can be applied are: change the hysteresis h or to introduce a proper compensator in the linear part of the system.
Regarding to the hysteresis, its effect on −1/N A has been explained in Section III-A, expanding or folding the traces. Alternatively, a compensator can be introduced into the loop in order to guarantee the intersection between G ol (jω) and −1/N A at a given frequency.
Both of the approaches to modify the oscillation frequency entail a modification of the amplitude, therefore, the described procedure to modify the oscillation amplitude must be applied afterwards.

VII. CONCLUSION
In this work, the tools to characterize the limit cycle oscillations induced by any kind of multilevel Fixed Threshold Sampler (FTS) have been provided. To that end, the structure of these kinds of samplers have been parametrized and it has been shown that by choosing some specific set of parameters some of the most used quantifiers can be obtained. Some notable specific cases are RQ quantization, Symmetric-Send-On-Delta or any kind of multi-level relays.
In order to characterize the oscillations induced by these samplers the Describing Function technique has been used. Since depending on the set of parameters chosen for the sampler oscillations with a bias can be obtained, a variant called the Dual Input Describing Function (DIDF) has been used, which already contemplates the apparition of this bias.
The DIDF has been obtained for the generalized parameters describing any FTS. In addition, the DIDF involving the dynamic part has been deeply analyzed, revealing some interesting effects of the parameters on its traces and some properties such as periodicity and symmetry depending on the bias.
The guidelines to perform the robustness analysis have been provided. In the analysis study, it has been shown the paramount importance of the loop configuration and the signals involved in the loop, because they have a direct effect on the limit cycle oscillation conditions. Some examples with specific cases have been included to show how the analysis methodology should be applied. From these examples some features like the appearance of different types of oscillations produced by the application of some signals have been revealed, which can be used in further works for identification purposes. .

APPENDIX A DIDF CALCULATION
The input-output relationship of a FTS whose main parameters are represented in Figure 2 is described in (1). Let the input to the FTS be described by x(φ) = A sin(φ) + B. Then, the output can be expressed as: where i 0 is the level around which the oscillation is centered. Under this assumption of the form of x(φ) and because of a possible output bias that propitiates this type of oscillation the Dual Input Describing Function (DIDF) must be used. This variant of the Describing Function has two components, one related to the sustainability of the bias, and another one related to the dynamic part of the signal, which can be respectively calculated with the following expressions: (φ)e −jφ dφ.
Developing, for example, for N A as can be seen in (11), where l = 2m sup + 2m inf being m sup and m inf the num- ber of levels crossed in upwards and downwards direction respectively. Arranging terms: The first integral equals 0, and for the second, it is known that: This leads to the following expression: The same procedure is followed for obtaining the expression of N B : The expressions for m sup and m inf can be calculated with: and the initial level i 0 with:

APPENDIX B CONDITIONS ON THE BIAS
Depending on the process and controller in the loop there exist different situations which define the necessary conditions to obtain a sustained bias. The disturbance and reference signals have been considered to be step-like signals of magnitude p(∞) and r(∞) respectively, resulting in the oscillation conditions to be a function of the process and controller involved. Then, the different combinations between these elements and their respective oscillation condition are presented in Table 2.  2: Expanded requirements for any combination of process and controller for the maintenance of the bias under step-like changes in reference and disturbance inputs. L 1 (s) and L 2 (s) have neither poles nor zeros at the origin and it is considered that n, m ≥ 1. From 1995 to 1998 taught in the department of automatics and computational systems of Universidad Central de las Villas and since 2004 he has been at Department of Industrial Systems Engineering and Design in University Jaume I (Spain) as a professor. His scientific activities include tuning and auto-tuning methods for industrial controllers, distributed control systems and event-based controllers.
16 VOLUME 4, 2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and JOSÉ SÁNCHEZ-MORENO received the computer sciences degree in 1994 from Madrid Polytechnic University and the PhD degree in sciences from UNED in 2001 with a thesis on the development of virtual and remote labs for teaching automatic control across the Internet. Since 1993, he has been at UNED Department of Computer Sciences and Automatic Control as an assistant professor. His main research interests for the time being: event-based control, networked control systems, remote and virtual laboratories in control engineering, and pattern recognition in nuclear fusion databases.
SEBASTIÁN DORMIDO received the BS degree in physics from Complutense University, Madrid, Spain, in 1968 and the PhD degree in the sciences from Basque Country University, Bilbao, Spain, in 1971. He received a Doctor Honorary degree from the Universidad de Huelva and Universidad de Almería. In 1981, he was appointed as a professor of control engineering at UNED, Madrid. His scientific activities include computer control of industrial processes, model-based predictive control, hybrid control, and web-based labs for distance education. He has authored or coauthored more than 250 technical papers in international journals and conferences and has supervised more than 35 PhD thesis. From 2001 to 2006, he was the president of the Spanish Association of Automatic Control, CEA-IFAC. He received the National Automatic Control prize from IFAC Spanish Automatic Control Committee. VOLUME 4, 2016 17 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3182794