Fitting of Generic Process Models by an Asymmetric Short Relay Feedback Experiment—The n -Shifting Method

: This paper presents the generalization of the shifting method for relay feedback identiﬁcation of dynamic systems of any order. The original shifting method enables the ﬁtting of a maximum of ﬁve parameters of a transfer function model from the information obtained from a short relay test and without prior knowledge of the process to identify. The generalization, known as n -shifting, allows the estimation of the parameters of transfer functions of any order by applying one short relay test to the process to identify. Without loss of generality, the n -shifting approach is applied to ﬁt an n -order plus time delay ( n -OPTD) model but the approach can be also developed to identify models with other structures (non-minimum phase, unstable, integrators). Some examples of estimations are presented.


Introduction
Experimental system identification of transfer function models using relays is a widely used technique in control engineering. Typskin originally described that technique in 1973 [1], although the most known publication applying the approach is [2], where the method is applied for the automatic tuning of simple PID (Proportional, Integral, Derivative) regulators.
Due to the high number of regular papers, conference contributions, technical reports and specific books on system identification by relay methods, it is not our aim to make a long introduction to the subject. The rationale of the method is that a linear system under an ideal symmetrical relay control (u A = −u B and e A = e B = 0 in Figure 1) oscillates, approximately, at its ultimate frequency, that is, ω osc ≈ ω u and the critical gain K u is derived from the describing function (DF) of the ideal relay. That is: where A osc is the oscillation amplitude and u A is the relay output. The reason of the oscillation is explained because it occurs when Equation (2): An easy way to get the parameters of a transfer function model from Equation (2) is to use the magnitude and argument of ( ), separating both components to solve them. That is: As an example of this approach, the following two equations are obtained for estimating the parameters of a first order process as : Supposing that the process gain is known, it is possible to calculate first the lag T using Equation (5) and next the time delay using Equation (6). If the process to identify is an underdamped second order plus time delay as ( + 1)( + 1) ⁄ , it is necessary to force the system to oscillate at another frequency to estimate the four parameters , , and of the model, by solving the four equations obtained from the magnitudes and arguments of: = . (8) There are two techniques for generating oscillations at other frequencies different from : (i) by adding a phase lag to the original critical point ( ) and (ii) by using a relay with hysteresis. The phase lag can be introduced in ( ) in two ways: by including an additional transport delay [3][4][5][6] or by an integrator [7][8][9]. Regarding the hysteresis technique, the DF corresponding to a relay with hysteresis is: where = = ≠ 0 and > . As −1/ ( , ) represents a horizontal line that can be moved down along the imaginary axis of the Nyquist plot changing . As the intersection of ( ) with −1/ ( , ) means an oscillation, the critical point where the intersection occurs can be modified by manipulating the width of the hysteresis. Replacing Equation (9) in Equation (7) and Equation (8), it is possible to decrease the phase lag of the critical point as needed. From Equation (9), the theoretical phase lag obtained by applying hysteresis to the relay is defined by: = arctan = arcsin . An easy way to get the parameters of a transfer function model from Equation (2) is to use the magnitude and argument of G(jω osc ), separating both components to solve them. That is: arg −1 N(A osc ) = argG(jω osc ). (4) As an example of this approach, the following two equations are obtained for estimating the parameters of a first order process as Ke −Ls Ts+1 : Supposing that the process gain K is known, it is possible to calculate first the lag T using Equation (5) and next the time delay L using Equation (6). If the process to identify is an underdamped second order plus time delay as Ke −Ls /[(T 1 s + 1)(T 2 s + 1)], it is necessary to force the system to oscillate at another frequency ω osc 2 to estimate the four parameters K, T 1 , T 2 and L of the model, by solving the four equations obtained from the magnitudes and arguments of: There are two techniques for generating oscillations at other frequencies different from ω osc : (i) by adding a phase lag to the original critical point G(jω osc ) and (ii) by using a relay with hysteresis. The phase lag can be introduced in G(s) in two ways: by including an additional transport delay [3][4][5][6] or by an integrator [7][8][9]. Regarding the hysteresis technique, the DF corresponding to a relay with hysteresis is: where δ = e A = e B = 0 and A osc > δ. As −1/N(A osc , δ) represents a horizontal line that can be moved down along the imaginary axis of the Nyquist plot changing δ. As the intersection of G(s) with −1/N(A osc , δ) means an oscillation, the critical point where the intersection occurs can be modified by manipulating the width of the hysteresis. Replacing Equation (9) in Equation (7) and Equation (8), it is possible to decrease the phase lag of the critical point as needed. From Equation (9), the theoretical phase lag obtained by applying hysteresis to the relay is defined by: This solution was originally described by [2]. However, the DF approach provides only an approximation of the true critical point at the frequencies where G(s) oscillates. These approximations reduce the accuracy of the results or even produce solutions in the complex plane. For instance, the DF analysis in a simple relay considers that the critical point has a phase lag of 180 • but in a real experiment the phase lag can differ some degrees from 0 • . Let see a practical example of this situation.
Example: The first order process: under relay control feedback with u A = −u B = 1 and e A = e B = 0 oscillates in a simulation at the frequency ω osc = 2.1056 and amplitude A osc = 0.6327. However, by the DF approach, it is known that: and argG(jω osc ) = −180 • . Supposing that the process gain K = 1 is known and using the information from the DF and Equations (5) and (6) to estimate the lag and the delay of G(s), we would obtain T = 0.8002 and L = 1.9842. Clearly, the estimation is not very accurate. Also, if the process to estimate had been a second order model with time delay, an additional point G(jω osc 2 ) would have been necessary by a second experiment with the relay. However, in this case also the estimation would contain similar errors. From this example, it is clear that it is essential to obtain the exact value of G(s) at different oscillation frequencies during one relay test.
Two strategies are the most relevant solutions found in the literature to measure accurately several critical points during just one relay experiment. The first one is based on the Laplace transform of y(t) and u(t); as both signals are periodic and piecewise, it can be written as: and so, it is possible to obtain the harmonics: It is worth noting that if the relay is not biased, that is, u A = u B , the signals and u(t) are symmetric in Equation (11) and thus harmonics of even order cancel out. Some works on system identification based on these approaches are [10][11][12][13][14]. If this strategy had been applied in the previous example to estimate the critical point, it would have been obtained |G(jω osc )| = 0.4289965 and argG(jω osc ) = 174.56 • . With these data, the estimation of the process parameters produces T = 0.9999 and L = 1.001. Notice the difference in the accuracy with respect to the previous situation.
A second solution to get a second harmonic is the shifting method described in [15,16]. By this approach, it is possible to determine two points G(jω osc ) and G(j2ω osc ) during one relay test by obtaining a second set of signals of frequency 2ω osc directly from y(t) and u(t). The generalization of this method to obtain the points G(jkω osc ) with k = 1 . . . n is the subject of this paper.
The original shifting and the extension presented here are characterized for their significantly short experiments (in fact, only a single experiment is needed). Once the convergence to the limit cycle has been achieved, only the information of one cycle of y(t) and u(t) is needed to make possible the identification of the model. It was demonstrated in [17] that if the experiment starts in a steady state and it is executed in the absence of load disturbances, then it is possible to use very short experiments, even under significant measurement noise. In general, the convergence is reached in one to three periods [17,18]. Similar works on identification based on short relay experiments are described in [19,20].
This paper presents the n-shifting method, a generalization of the basic shifting algorithm to obtain n points of G(s) in a short relay feedback experiment. Thus, the original method is extended to estimate transfer functions with more than 5 parameters, which is the maximum possible with the original approach.
The paper is organized as follows. Section 2 is dedicated to describing with more detail the original shifting technique. Section 3 demonstrates why it is possible to extend the shifting method to obtain n points from G(s) with data from just only one cycle of the oscillation generated by the relay in the system. The identification algorithm is explained in Section 4 for fitting models with n lags and one time delay and three illustrative examples are given in Section 5. Finally, conclusions and future extensions of the method are provided in Section 6.

The Shifting Technique
The original shifting procedure is based on generating, during a test with an asymmetric relay with hysteresis (see Figure 1), a new set of periodic signals y(t) and u(t) of frequency 2ω osc directly from the signals y(t) and u(t) of frequency ω osc where ω osc = 2π/T. These two new signals are obtained by the sum of the original signals and themselves but with a delay of T/2: The result of Equations (12) and (13) is that y(t) and u(t) are periodic and rectangular waveforms, respectively but of T/2 period. Figure 2 shows an example of the shifting of the signals corresponding to a process 2e −s /(s + 1) where e A = −e B = 1, u A = 1 and u B = −0.8.
Since y(t) and u(t) are periodic, it is possible to apply Equation (11) for obtaining the fundamental harmonics of G(s) when the system is oscillating at 2ω osc . So, two critical points can be generated from the information of only one cycle by the following expressions: At this point, it is important to remark that the relay must be asymmetric, that is, e A = −e B and/or u A = −u B to avoid the second critical point becoming infinite because u(t) = 0. If the relay is symmetrical the sum of the semi-periods of u(t) and u(t) in Equation (13) is zero because u(t) = −u(t − T/2). Also, if a third point of the Nyquist plot is required, the process gain K = G(0) can be obtained as follows: The integrals of the Equations (14), (15) and (16) can be computed numerically once the experiment is finished. During the experiment it is enough to save the information corresponding to ( ) and ( ) to be analysed and processed off-line to calculate the value of the three critical points.
Once the values of the three critical points are available, it is possible to estimate the parameters of a model composed of four or five parameters as a maximum. There are two alternatives to make the fitting. The first one is by solving the explicit formulas of the model parameters (see [12][13][14] for expressions of models of four parameters). A second alternative is by solving the following two complex nonlinear equations: The first version of the shifting method described in [15,16] did not take into consideration that the position of ( 2 ) could be located in the first or second quadrant of the Nyquist map (see Figure 3). Due to the high frequency noise, it could be preferable for system identification to have this point in the third quadrant or with a phase lag as close to 180° as possible. To get two points with a more convenient position in the map (with a phase lag from 180° to 90°) it is necessary to force a clockwise rolling of ( ) by adding an additional phase shift by an integrator or a delay. Such action produces a new system ′( ) whose intersection with the negative reciprocal of the relay DF occurs at a lower frequency than the original ( ). Once the point ′( ) is measured and ′( 2 ) Figure 2. Waveforms of four cycles of y(t) and u(t) (thick lines) and their corresponding shifting y(t) and u(t) (thin lines in red colour).
The integrals of the Equations (14), (15) and (16) can be computed numerically once the experiment is finished. During the experiment it is enough to save the information corresponding to y(t) and u(t) to be analysed and processed off-line to calculate the value of the three critical points.
Once the values of the three critical points are available, it is possible to estimate the parameters of a model composed of four or five parameters as a maximum. There are two alternatives to make the fitting. The first one is by solving the explicit formulas of the model parameters (see [12][13][14] for expressions of models of four parameters). A second alternative is by solving the following two complex nonlinear equations: The first version of the shifting method described in [15,16] did not take into consideration that the position of G(j2ω osc ) could be located in the first or second quadrant of the Nyquist map (see Figure 3). Due to the high frequency noise, it could be preferable for system identification to have this point in the third quadrant or with a phase lag as close to 180 • as possible. To get two points with a more convenient position in the map (with a phase lag from 180 • to 90 • ) it is necessary to force a clockwise rolling of G(s) by adding an additional phase shift by an integrator or a delay. Such action produces a new system G (s) whose intersection with the negative reciprocal of the relay DF occurs at a lower frequency than the original G(s). Once the point G (jω osc ) is measured and G (j2ω osc ) is calculated by shifting, the effect of the additional phase is discounted to obtain the original G(jω osc ) and G(j2ω osc ) and proceed with the fitting. So, for example, if the system G(s) is modified by adding a delay L add , Equation (17) must be replaced by: (18) to undo the phase shifting. The added phase is −ω osc L add for the point G(jω osc ) and −2ω osc L add for G(j2ω osc ). If an integrator is chosen, the lag is always −π/2. Figure 3 presents an example of phase shifting by inserting a delay of 0.4 that transforms the system G(s) = e −0.1s /(s + 1) into G (s) = e −0.5s /(s + 1). The oscillation occurs because G is intersecting with −1/N(A osc , δ) that corresponds to a relay with hysteresis. Notice that the two points of G(s) are located in the third and fourth quadrants of the map when the phase shifting is undone. If the experiment with G(s) had been done without the additional delay, the point G(j2ω osc ) would have been located in the second or first quadrant. is calculated by shifting, the effect of the additional phase is discounted to obtain the original ( ) and ( 2 ) and proceed with the fitting. So, for example, if the system ( ) is modified by adding a delay , Equation (17) must be replaced by: to undo the phase shifting. The added phase is − for the point ( ) and −2 for ( 2  ). If an integrator is chosen, the lag is always −π/2. Figure 3 presents an example of phase shifting by inserting a delay of 0.4 that transforms the system ( ) = . /( + 1) into ′( ) = . /( + 1). The oscillation occurs because ′ is intersecting with −1/ ( , ) that corresponds to a relay with hysteresis. Notice that the two points of ( ) are located in the third and fourth quadrants of the map when the phase shifting is undone. If the experiment with ( ) had been done without the additional delay, the point ( 2 ) would have been located in the second or first quadrant.
In this example,  This modification was included in the shifting method in [21] by adding an integrator to the block diagram and in [22,23] by a transport delay.

Extending the Basic Shifting Approach
The extension of the basic shifting method to generate n points means opening the possibility of estimating models of any order by one relay feedback experiment. This extension would consist in generating n sets of signals ( ) and ( ) by generalization of the Equations (12) and (13) to: in order to be able to apply the expression: This modification was included in the shifting method in [21] by adding an integrator to the block diagram and in [22,23] by a transport delay.

Extending the Basic Shifting Approach
The extension of the basic shifting method to generate n points means opening the possibility of estimating models of any order by one relay feedback experiment. This extension would consist in generating n sets of signals y n (t) and u n (t) by generalization of the Equations (12) and (13) to: This generalization is known as the n-shifting method. It must be noticed that with n = 2, the generalization corresponds to the original approach. Figure 4 shows the result of the n-shifting until n = 4 with the same process and relay setting as in Figure 2.
This generalization is known as the n-shifting method. It must be noticed that with n = 2, the generalization corresponds to the original approach. Figure 4 shows the result of the n-shifting until n = 4 with the same process and relay setting as in Figure 2. But in order to be able to apply the extension, it is necessary to prove that given a generic signal ( ) of period , the signal ̅ ( ) obtained by the sum of n signals ( ) delayed successively / with j = 0 … n−1 is a signal of period / . It leads to the formulation of the following theorem. Theorem 1. Given a periodic function ( ) having a period , that is, ( ) = ( + ), then: is a periodic function of period / . Proof of Theorem 1. If the ( ) is the Laplace transform of ( ) then the Laplace transform of ( − ) is ( ). So, the Laplace transform of (22) is: It is well known [24] that the Laplace transform of a periodic function ( ) = ( + ) is: Figure 4. Waveforms of four cycles of y(t) and u(t) (thick black) and their n-shifting signals y 2 (t), y 3 (t), y 4 (t) and u 2 (t), u 3 (t), u 4 (t) (thick red, thin black and thin red, respectively). But in order to be able to apply the extension, it is necessary to prove that given a generic signal f (t) of period T, the signal f (t) obtained by the sum of n signals f (t) delayed successively jT/n with j = 0 . . . n−1 is a signal of period T/n. It leads to the formulation of the following theorem.

Theorem 1.
Given a periodic function f (t) having a period T, that is, f (t) = f (t + T), then: is a periodic function of period T/n.
It is well known [24] that the Laplace transform of a periodic function f (t) = f (t + T) is: So, applying Equation (24) to the right term of Equation (23), we get: By applying the variable change e − T n s = x to Equation (25), we have: If: is applied to Equation (26), the following expression is obtained: Undoing the previous variable change in Equation (28), we get: that corresponds to the expression of the Laplace transform given in Equation (24) of a periodic function f (t) with period T/n.

The Identification Algorithm
Once proved that the generalization of the shifting method is feasible, it is necessary to define the identification method for fitting models of any order using the information provided by the n-shifting approach.
Obtaining the explicit expressions of the parameters of a model with more than four parameters (for instance, one zero, three lags and delay) can be complicated and tedious because of the length of the nonlinear expressions. The solution presented here is for fitting a model of n order plus time delay (n-OPTD): and it is based on solving a set of nonlinear equations. Without loss of generality, the procedure can be applied to models with other structures as non-minimum phase or with integrators. Continuing with Equation (30), the number of unknowns in this model is n + 1, L is not taken into account because it is obtained by solving the equation corresponding to the argument of Equation (30). Therefore, n + 1 experimental points G(0)), G(jω osc ), G(j2ω osc ), . . . , G(jnω osc ) must be obtained from a relay test by n-shifting. The magnitude expression of M n−OPTD (jnω osc ) is: Now, relating the magnitudes of the n points G(jω osc ), G(j2ω osc ), . . . , G(jnω osc ) with their corresponding magnitude expressions obtained from (31) and reordering the expressions, we obtain the following set of n nonlinear equations: with k = 1 . . . n. Once the T i s are obtained from Equation (32), we use these values and the argument of G(jω osc ) to get L. As it is possible to calculate: then: where it is considered that argG(jω osc ) < 0. Several valid models can be obtained on some occasions. The model selected must be the one that minimizes the following index based on the Euclidian distance: Summarizing, the generic procedure for fitting a model by the n-shifting approach is: 1. Make a short relay feedback experiment. The length of the experiment must be 2-3 cycles for stabilization and 1 cycle for recording y(t) and u(t).

2.
Calculate ω osc and save the last cycle of y(t) and u(t).
If multiple solutions are obtained, select the model that minimizes Equation (35).

Remark 1.
The phase lag of points G(j2ω osc ) to G(jnω osc ) in the Nyquist plot are conditioned by the relay because the hysteresis determines the position of G(jω osc ). The maximum decrease of phase lag of G(jω osc ) that the hysteresis can theoretically produce is 90 • . This could prevent the points of G at higher frequencies from being situated in the first and second quadrant, without thus being necessary to add any phase by a delay or an integrator.

Remark 2.
The solutions provided in the literature to force the phase shifting when needed do not allow to know the effect on the phase before starting the experiment. An on-the-fly way is starting the experiment with a small additional delay and increase it until: but that implies that G (jnω osc ) must be known on-line.

Remark 3.
The effect of the noise has not been considered in the procedure because results would be similar to those presented in [12,13] due to the same foundation of both methods. One of the advantages of the relay is that its hysteresis reduces the influence of the noise in e(t) by generating control actions u(t) free of noise. The recommendation is to adjust the hysteresis to 2-3 times the noise level [19,25]. However, if the noise level is too large, it could produce false control actions and the calculation of ω osc and G(jω osc ) would become difficult. This fact would demand a previous filtering of y(t).

Simulation Examples
The method is demonstrated on three examples of processes of high order. The experiments have been run in Matlab/Simulink by simulations and the integration step has been set to h = 0.001 s. The length of the experiments has been 3 cycles, leaving the first two cycles to allow the limit cycle to stabilize. So, just information from once cycle was used in the Equations (19)-(21). The nonlinear system given by Equations (32)  The method is applied to the estimation of a model of order n = 4, which means six parameters, as a study case for presenting the system of nonlinear equations to be solved. The process chosen is: so four points must be obtained by the n-shifting derived from the information of one cycle of y(t) and u(t). Applying Equation (32) with n = 4, it is obtained: |G(j2ω osc )| 2 + 1 + T 2 1 4ω 2 osc 1 + T 2 2 4ω 2 osc 1 + T 2 3 4ω 2 osc 1 + T 2 4 4ω 2 osc = 0 − K 2 |G(j3ω osc )| 2 + 1 + T 2 1 9ω 2 osc 1 + T 2 2 9ω 2 osc 1 + T 2 3 9ω 2 osc 1 + T 2 4 9ω 2 osc = 0 − K 2 |G(j4ω osc )| 2 + 1 + T 2 1 16ω 2 osc 1 + T 2 2 16ω 2 osc 1 + T 2 3 16ω 2 osc 1 + T 2 4 16ω 2 osc = 0.
Solving Equations (38) and (39), the model obtained is: Figure 5 presents the Nyquist plot of the process given by Equation (37) and the model given by Equation (39) with the points of G 1 (s) highlighted.

Example 2.
Fitting a model of lower order than the process using phase shifting.
In this second example, the 8th order process: is fitted using the model of a fourth order M 4-OPTD . An additional delay L add = 4 has been introduced to illustrate the effect of the phase shifting to move the G 2 points to a more convenient position in the Nyquist plot. Therefore, the simulation is run during three cycles with the system G 2 (s) = G 2 (s)e −sL add . The configuration of the relay during the experiment has been set to: e a = 0.5, e b = − 0.5, u a = 1, u b = − 0.8.
After running the experiment and applying the n-shifting approach to one cycle of data from ( ) and ( ), the following values were measured: ωosc = 0.892878 rad/s, G1(0)  Example 2: Fitting a model of lower order than the process using phase shifting In this second example, the 8 th order process: As the model to estimate is of order 4, the points to calculate are 4. Thus, the results of the n-shifting are: ω osc = 0.2238159 rad/s, G' 2 (0) = 1.000020.
It can be appreciated the phase shifting that the additional delay has introduced in the G 2 points when compared with the G 2 s. Solving Equations (32) and (34) with the G 2 points, many solutions are obtained but the only valid solution that corresponds to a stable system is the following: 1.00002 14.52s 4 + 21.96s 3 + 17.17s 2 + 6.5s + 1 e −1.504s .
The rest of the solutions correspond to unstable and/or non-causal systems (negative delays) and are thus discarded. Figures 6 and 7 present the Nyquist plot and the step response of the process G 2 and the model of Equation (41). The behaviour of the model in time and frequency is completely similar to the original process.

Example 3. Estimation of an unstable process of low order.
This last example is included to show that the proposed method is able to identify unstable plants from one relay feedback test. The process selected is a second order unstable plant: presented in [26], where T 1 = 2, T 2 = 0.5 and L = 0.5. The configuration of the relay is the one used by [26]: e a = 0.1, e b = − 0.1, u a = 1, u b = − 1. . (41) The rest of the solutions correspond to unstable and/or non-causal systems (negative delays) and are thus discarded. Figures 6 and 7 present the Nyquist plot and the step response of the process G2 and the model of Equation (41). The behaviour of the model in time and frequency is completely similar to the original process.  Step responses of of ( ) (in blue) and ( ) (in red).
Example 3: Estimation of an unstable process of low order This last example is included to show that the proposed method is able to identify unstable plants from one relay feedback test. The process selected is a second order unstable plant: presented in [26], where = 2, = 0.5 and = 0.5. The configuration of the relay is the one used by [26]: To select one, it is necessary to calculate the index given by Equation (35) for both models. The values are = 0.0362 and = 0.1607 , so the selected model is ( ). Figure 8 shows the Nyquist plot of ( ) and the two solutions where it is clearly appreciated that the model ( ) is the best option. After the simulation, the points obtained from applying the n-shifting are: ω osc = 0.592697416015431 rad/s, G 3 (0) = − 0.999979, G 3 (jω osc ) = − 0.593316 − 0.173648j, G 3 (j2ω osc ) = − 0.333983 − 0.014068j, G 3 (j3ω osc ) = − 0.191932 + 0.063947j, G 3 (j4ω osc ) = − 0.102077 + 0.085296j.
As in Example 2, the solutions obtained from Equations (32) and (34) with n = 4 are many but only two correspond to an unstable system. Those are: To select one, it is necessary to calculate the D index given by Equation (35) for both models. The values are D 1 = 0.0362 and D 2 = 0.1607, so the selected model is M 1 4−OPTD (s). Figure 8 shows the Nyquist plot of G 3 (s) and the two solutions where it is clearly appreciated that the model M 1 4−OPTD (s) is the best option.

Conclusions
This paper has presented the n-shifting approach, a generalization of the basic shifting approach presented in [15], that allows the fitting of transfer functions by a short relay feedback experiment. The basic shifting technique allows us to fit models of a maximum of five parameters by two points of G(s) obtained experimentally. However, by the n-shifting extension it is possible to obtain n points and make a fitting of a model of any order and structure.
In this work, the fitting of a generic n-order with time delay (n-OTD) model has been taken as a study case to present the approach but it can be extended to stable models with other structures (for example with integrators or zeros) or even unstable processes. The simulations presented show that the approach may be very appropriate to fit model templates in contexts where long experiments cannot be done. So, in this context, the models of the examples have been fitted with information from one cycle of simulation after two cycles for stabilization. To resume the approach, from n points derived from one cycle relay feedback experiment it is feasible to estimate a n-OPTD model composed of n + 2 parameters, that is, the static gain K, the transport delay L and the n time lags.
Apart from the fitting of models of high order, a practical use of the n-shifting approach is its application to the tuning of PID controllers by evolutionary algorithms. As the approach is able to provide an accurate estimation of n points of the output spectrum of a process, the tuning of PID controllers with a genetic algorithm using a n-order model is a direct application. The procedure for tuning would be (a) to estimate n points of an unknown process by the n-shifting approach, (b) to generate the n-order model with the points obtained in the previous step and (c) to tune the PID controller by a genetic algorithm using the n-order model with the goal of minimizing the integrated absolute error (IAE) after a set-point change or after a disturbance.
Regarding its limitations, the noise could be a problem in presence of relevant measurement noise. As said before, one of the features of an asymmetric relay is that it reduces the influence of the measurement noise by producing a control action free of noise in most of the cases. However, in presence of significant noise would be convenient a previous filtering of the process output because the signals ( ) are synthetized by the n-shifting of ( ) and the influence of noise could be transmitted to these signals. Another limitation of this estimation approach (and the approaches based on relays) is to obtain points of the process located in the low frequency range (the fourth quadrant in the Nyquist

Conclusions
This paper has presented the n-shifting approach, a generalization of the basic shifting approach presented in [15], that allows the fitting of transfer functions by a short relay feedback experiment. The basic shifting technique allows us to fit models of a maximum of five parameters by two points of G(s) obtained experimentally. However, by the n-shifting extension it is possible to obtain n points and make a fitting of a model of any order and structure.
In this work, the fitting of a generic n-order with time delay (n-OTD) model has been taken as a study case to present the approach but it can be extended to stable models with other structures (for example with integrators or zeros) or even unstable processes. The simulations presented show that the approach may be very appropriate to fit model templates in contexts where long experiments cannot be done. So, in this context, the models of the examples have been fitted with information from one cycle of simulation after two cycles for stabilization. To resume the approach, from n points derived from one cycle relay feedback experiment it is feasible to estimate a n-OPTD model composed of n + 2 parameters, that is, the static gain K, the transport delay L and the n time lags.
Apart from the fitting of models of high order, a practical use of the n-shifting approach is its application to the tuning of PID controllers by evolutionary algorithms. As the approach is able to provide an accurate estimation of n points of the output spectrum of a process, the tuning of PID controllers with a genetic algorithm using a n-order model is a direct application. The procedure for tuning would be (a) to estimate n points of an unknown process by the n-shifting approach, (b) to generate the n-order model with the points obtained in the previous step and (c) to tune the PID controller by a genetic algorithm using the n-order model with the goal of minimizing the integrated absolute error (IAE) after a set-point change or after a disturbance.
Regarding its limitations, the noise could be a problem in presence of relevant measurement noise. As said before, one of the features of an asymmetric relay is that it reduces the influence of the measurement noise by producing a control action free of noise in most of the cases. However, in presence of significant noise would be convenient a previous filtering of the process output because the signals y(t) are synthetized by the n-shifting of y(t) and the influence of noise could be transmitted to these signals. Another limitation of this estimation approach (and the approaches based on relays) is to obtain points of the process located in the low frequency range (the fourth quadrant in the Nyquist map). This is consequence that relays with hysteresis can produce oscillations with a minimum phase lag of around 90 • .
Further research lines are (i) the extension of the n-shifting approach to the identification in event-based control loops where a SSOD-PI (Symmetric-Send-On-Delta) controller would replace the asymmetric relay, (ii) the tuning of PID controllers by genetic algorithms using n-OPTD models obtained by the n-shifting approach and (iii) the full characterization of the waveforms u(t) that the n-shifting produce by changing the relay parameters and the value of n.